diff --git a/GNUmakefile b/GNUmakefile index d7fd5b6a9..b5dcabdb1 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -1,4 +1,4 @@ -AMREX_HOME ?= ../amrex +AMREX_HOME = ../../amrex DEBUG = FALSE #DEBUG = TRUE diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp index 249f204a9..97a681a96 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp @@ -115,8 +115,10 @@ FlushFormatCheckpoint::WriteToFile ( if (warpx.getis_synchronized() || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon - || warpx.use_lumped_inductor == 1) { + || warpx.use_lumped_inductor == 1 + || warpx.use_josephson_junction == 1) { // Need to save j if synchronized because after restart we need j to evolve E by dt/2. + // London / lumped inductor / Josephson junction all treat j as persistent state. VisMF::Write(warpx.getcurrent_fp(lev, 0), amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jx_fp")); VisMF::Write(warpx.getcurrent_fp(lev, 1), @@ -180,7 +182,8 @@ FlushFormatCheckpoint::WriteToFile ( if (warpx.getis_synchronized() || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon - || warpx.use_lumped_inductor == 1) { + || warpx.use_lumped_inductor == 1 + || warpx.use_josephson_junction == 1) { // Need to save j if synchronized because after restart we need j to evolve E by dt/2. VisMF::Write(warpx.getcurrent_cp(lev, 0), amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jx_cp")); @@ -205,6 +208,12 @@ FlushFormatCheckpoint::WriteToFile ( } } + // Josephson junction state (phi, Js). phi is the junction's only real state and + // is NOT reconstructible from the fields, so without this a restart resets it to 0. + if (warpx.use_josephson_junction == 1) { + warpx.getJosephsonJunction().WriteCheckpoint(checkpointname); + } + CheckpointParticles(checkpointname, particle_diags); WriteDMaps(checkpointname, nlev); diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp.orig b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp.orig new file mode 100644 index 000000000..249f204a9 --- /dev/null +++ b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp.orig @@ -0,0 +1,281 @@ +#include "FlushFormatCheckpoint.H" + +#include "BoundaryConditions/PML.H" +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) +# include "BoundaryConditions/PML_RZ.H" +#endif +#include "Diagnostics/ParticleDiag/ParticleDiag.H" +#include "Particles/WarpXParticleContainer.H" +#include "Utils/TextMsg.H" +#include "Utils/WarpXProfilerWrapper.H" +#include "WarpX.H" + +#include +#include +#include +#include +#include +#include +#include + +using namespace amrex; + +namespace +{ + const std::string default_level_prefix {"Level_"}; +} + +void +FlushFormatCheckpoint::WriteToFile ( + const amrex::Vector /*varnames*/, + const amrex::Vector& /*mf*/, + amrex::Vector& geom, + const amrex::Vector iteration, const double /*time*/, + const amrex::Vector& particle_diags, int nlev, + const std::string prefix, int file_min_digits, + bool /*plot_raw_fields*/, + bool /*plot_raw_fields_guards*/, + const bool /*use_pinned_pc*/, + bool /*isBTD*/, int /*snapshotID*/, + int /*bufferID*/, int /*numBuffers*/, + const amrex::Geometry& /*full_BTD_snapshot*/, + bool /*isLastBTDFlush*/, const amrex::Vector& /* totalParticlesFlushedAlready*/) const +{ + WARPX_PROFILE("FlushFormatCheckpoint::WriteToFile()"); + + auto & warpx = WarpX::GetInstance(); + + VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); + VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1); + + const std::string& checkpointname = amrex::Concatenate(prefix, iteration[0], file_min_digits); + + amrex::Print() << Utils::TextMsg::Info( + "Writing checkpoint " + checkpointname); + + // const int nlevels = finestLevel()+1; + amrex::PreBuildDirectorHierarchy(checkpointname, default_level_prefix, nlev, true); + + WriteWarpXHeader(checkpointname, geom); + + WriteJobInfo(checkpointname); + + for (int lev = 0; lev < nlev; ++lev) + { + VisMF::Write(warpx.getEfield_fp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Ex_fp")); + VisMF::Write(warpx.getEfield_fp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Ey_fp")); + VisMF::Write(warpx.getEfield_fp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Ez_fp")); + VisMF::Write(warpx.getBfield_fp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Bx_fp")); + VisMF::Write(warpx.getBfield_fp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "By_fp")); + VisMF::Write(warpx.getBfield_fp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Bz_fp")); + +#ifdef WARPX_MAG_LLG + VisMF::Write(warpx.getHfield_fp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hx_fp")); + VisMF::Write(warpx.getHfield_fp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hy_fp")); + VisMF::Write(warpx.getHfield_fp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hz_fp")); + VisMF::Write(warpx.getMfield_fp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Mx_fp")); + VisMF::Write(warpx.getMfield_fp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "My_fp")); + VisMF::Write(warpx.getMfield_fp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Mz_fp")); + VisMF::Write(warpx.getH_biasfield_fp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hxbias_fp")); + VisMF::Write(warpx.getH_biasfield_fp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hybias_fp")); + VisMF::Write(warpx.getH_biasfield_fp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hzbias_fp")); +#endif + + if (WarpX::fft_do_time_averaging) + { + VisMF::Write(warpx.getEfield_avg_fp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Ex_avg_fp")); + VisMF::Write(warpx.getEfield_avg_fp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Ey_avg_fp")); + VisMF::Write(warpx.getEfield_avg_fp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Ez_avg_fp")); + + VisMF::Write(warpx.getBfield_avg_fp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Bx_avg_fp")); + VisMF::Write(warpx.getBfield_avg_fp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "By_avg_fp")); + VisMF::Write(warpx.getBfield_avg_fp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Bz_avg_fp")); + } + + if (warpx.getis_synchronized() + || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon + || warpx.use_lumped_inductor == 1) { + // Need to save j if synchronized because after restart we need j to evolve E by dt/2. + VisMF::Write(warpx.getcurrent_fp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jx_fp")); + VisMF::Write(warpx.getcurrent_fp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jy_fp")); + VisMF::Write(warpx.getcurrent_fp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jz_fp")); + } + + if (lev > 0) + { + VisMF::Write(warpx.getEfield_cp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Ex_cp")); + VisMF::Write(warpx.getEfield_cp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Ey_cp")); + VisMF::Write(warpx.getEfield_cp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Ez_cp")); + VisMF::Write(warpx.getBfield_cp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Bx_cp")); + VisMF::Write(warpx.getBfield_cp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "By_cp")); + VisMF::Write(warpx.getBfield_cp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Bz_cp")); + +#ifdef WARPX_MAG_LLG + VisMF::Write(warpx.getHfield_cp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hx_cp")); + VisMF::Write(warpx.getHfield_cp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hy_cp")); + VisMF::Write(warpx.getHfield_cp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hz_cp")); + VisMF::Write(warpx.getMfield_cp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Mx_cp")); + VisMF::Write(warpx.getMfield_cp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "My_cp")); + VisMF::Write(warpx.getMfield_cp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Mz_cp")); + VisMF::Write(warpx.getH_biasfield_cp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hxbias_fp")); + VisMF::Write(warpx.getH_biasfield_cp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hybias_fp")); + VisMF::Write(warpx.getH_biasfield_cp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hzbias_fp")); +#endif + + if (WarpX::fft_do_time_averaging) + { + VisMF::Write(warpx.getEfield_avg_cp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Ex_avg_cp")); + VisMF::Write(warpx.getEfield_avg_cp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Ey_avg_cp")); + VisMF::Write(warpx.getEfield_avg_cp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Ez_avg_cp")); + + VisMF::Write(warpx.getBfield_avg_cp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Bx_avg_cp")); + VisMF::Write(warpx.getBfield_avg_cp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "By_avg_cp")); + VisMF::Write(warpx.getBfield_avg_cp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Bz_avg_cp")); + } + + if (warpx.getis_synchronized() + || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon + || warpx.use_lumped_inductor == 1) { + // Need to save j if synchronized because after restart we need j to evolve E by dt/2. + VisMF::Write(warpx.getcurrent_cp(lev, 0), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jx_cp")); + VisMF::Write(warpx.getcurrent_cp(lev, 1), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jy_cp")); + VisMF::Write(warpx.getcurrent_cp(lev, 2), + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jz_cp")); + } + } + + if (warpx.DoPML()) { + if (warpx.GetPML(lev)) { + warpx.GetPML(lev)->CheckPoint( + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "pml")); + } +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + if (warpx.GetPML_RZ(lev)) { + warpx.GetPML_RZ(lev)->CheckPoint( + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "pml_rz")); + } +#endif + } + } + + CheckpointParticles(checkpointname, particle_diags); + + WriteDMaps(checkpointname, nlev); + + VisMF::SetHeaderVersion(current_version); + +} + +void +FlushFormatCheckpoint::CheckpointParticles ( + const std::string& dir, + const amrex::Vector& particle_diags) const +{ + for (auto& part_diag: particle_diags) { + WarpXParticleContainer* pc = part_diag.getParticleContainer(); + + Vector real_names; + Vector int_names; + Vector int_flags; + Vector real_flags; + + real_names.push_back("weight"); + + real_names.push_back("momentum_x"); + real_names.push_back("momentum_y"); + real_names.push_back("momentum_z"); + +#ifdef WARPX_DIM_RZ + real_names.push_back("theta"); +#endif + + // get the names of the real comps + real_names.resize(pc->NumRealComps()); + auto runtime_rnames = pc->getParticleRuntimeComps(); + for (auto const& x : runtime_rnames) { real_names[x.second+PIdx::nattribs] = x.first; } + + // and the int comps + int_names.resize(pc->NumIntComps()); + auto runtime_inames = pc->getParticleRuntimeiComps(); + for (auto const& x : runtime_inames) { int_names[x.second+0] = x.first; } + + pc->Checkpoint(dir, part_diag.getSpeciesName(), true, + real_names, int_names); + } +} + +void +FlushFormatCheckpoint::WriteDMaps (const std::string& dir, int nlev) const +{ + if (ParallelDescriptor::IOProcessor()) { + auto & warpx = WarpX::GetInstance(); + for (int lev = 0; lev < nlev; ++lev) { + std::string DMFileName = dir; + if (!DMFileName.empty() && DMFileName[DMFileName.size()-1] != '/') {DMFileName += '/';} + DMFileName = amrex::Concatenate(DMFileName + "Level_", lev, 1); + DMFileName += "/DM"; + + std::ofstream DMFile; + DMFile.open(DMFileName.c_str(), std::ios::out|std::ios::trunc); + + if (!DMFile.good()) { amrex::FileOpenFailed(DMFileName); } + + DMFile << ParallelDescriptor::NProcs() << "\n"; + warpx.DistributionMap(lev).writeOn(DMFile); + + DMFile.flush(); + DMFile.close(); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + DMFile.good(), + "FlushFormatCheckpoint::WriteDMaps: problem writing DMFile" + ); + } + } +} diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp index f9814c84a..489c44316 100644 --- a/Source/Diagnostics/FullDiagnostics.cpp +++ b/Source/Diagnostics/FullDiagnostics.cpp @@ -716,6 +716,18 @@ FullDiagnostics::InitializeFieldFunctors (int lev) m_all_field_functors[lev][comp] = std::make_unique(warpx.getInductor().m_inductor_y_mf.get(), lev, m_crse_ratio); } else if ( m_varnames[comp] == "inductorz" ){ m_all_field_functors[lev][comp] = std::make_unique(warpx.getInductor().m_inductor_z_mf.get(), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "josephson_phi_x" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getJosephsonJunction().m_phi_x_mf.get(), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "josephson_phi_y" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getJosephsonJunction().m_phi_y_mf.get(), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "josephson_phi_z" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getJosephsonJunction().m_phi_z_mf.get(), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "josephson_Ic_x" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getJosephsonJunction().m_Ic_x_mf.get(), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "josephson_Ic_y" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getJosephsonJunction().m_Ic_y_mf.get(), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "josephson_Ic_z" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getJosephsonJunction().m_Ic_z_mf.get(), lev, m_crse_ratio); } else { diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index ed617acaf..025ac5393 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -359,7 +359,8 @@ WarpX::InitFromCheckpoint () if (is_synchronized || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon - || use_lumped_inductor == 1) + || use_lumped_inductor == 1 + || use_josephson_junction == 1) { VisMF::Read(*current_fp[lev][0], amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jx_fp")); @@ -426,7 +427,8 @@ WarpX::InitFromCheckpoint () if (is_synchronized || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon - || use_lumped_inductor == 1) { + || use_lumped_inductor == 1 + || use_josephson_junction == 1) { VisMF::Read(*current_cp[lev][0], amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jx_cp")); VisMF::Read(*current_cp[lev][1], diff --git a/Source/Diagnostics/WarpXIO.cpp.orig b/Source/Diagnostics/WarpXIO.cpp.orig new file mode 100644 index 000000000..ed617acaf --- /dev/null +++ b/Source/Diagnostics/WarpXIO.cpp.orig @@ -0,0 +1,459 @@ +/* Copyright 2019-2020 Andrew Myers, Ann Almgren, Axel Huebl + * Burlen Loring, David Grote, Gunther H. Weber + * Junmin Gu, Maxence Thevenet, Remi Lehe + * Revathi Jambunathan, Weiqun Zhang + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "BoundaryConditions/PML.H" +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) +# include "BoundaryConditions/PML_RZ.H" +#endif +#include "FieldIO.H" +#include "Particles/MultiParticleContainer.H" +#include "Utils/TextMsg.H" +#include "Utils/WarpXProfilerWrapper.H" +#include "WarpX.H" +#include "Diagnostics/MultiDiagnostics.H" + +#include + +#ifdef AMREX_USE_SENSEI_INSITU +# include +#endif +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +using namespace amrex; + +namespace +{ + const std::string level_prefix {"Level_"}; +} + +void +WarpX::GotoNextLine (std::istream& is) +{ + constexpr std::streamsize bl_ignore_max { 100000 }; + is.ignore(bl_ignore_max, '\n'); +} + +amrex::DistributionMapping +WarpX::GetRestartDMap (const std::string& chkfile, const amrex::BoxArray& ba, int lev) const { + std::string DMFileName = chkfile; + if (!DMFileName.empty() && DMFileName[DMFileName.size()-1] != '/') {DMFileName += '/';} + DMFileName = amrex::Concatenate(DMFileName + "Level_", lev, 1); + DMFileName += "/DM"; + + if (!amrex::FileExists(DMFileName)) { + return amrex::DistributionMapping{ba, ParallelDescriptor::NProcs()}; + } + + Vector fileCharPtr; + ParallelDescriptor::ReadAndBcastFile(DMFileName, fileCharPtr); + std::string fileCharPtrString(fileCharPtr.dataPtr()); + std::istringstream DMFile(fileCharPtrString, std::istringstream::in); + if ( ! DMFile.good()) amrex::FileOpenFailed(DMFileName); + DMFile.exceptions(std::ios_base::failbit | std::ios_base::badbit); + + int nprocs_in_checkpoint; + DMFile >> nprocs_in_checkpoint; + if (nprocs_in_checkpoint != ParallelDescriptor::NProcs()) { + return amrex::DistributionMapping{ba, ParallelDescriptor::NProcs()}; + } + + amrex::DistributionMapping dm; + dm.readFrom(DMFile); + if (dm.size() != ba.size()) { + return amrex::DistributionMapping{ba, ParallelDescriptor::NProcs()}; + } + + return dm; +} + +void +WarpX::InitFromCheckpoint () +{ + WARPX_PROFILE("WarpX::InitFromCheckpoint()"); + + amrex::Print()<< Utils::TextMsg::Info( + "restart from checkpoint " + restart_chkfile); + + // Header + { + std::string File(restart_chkfile + "/WarpXHeader"); + + VisMF::IO_Buffer io_buffer(VisMF::GetIOBufferSize()); + + Vector fileCharPtr; + ParallelDescriptor::ReadAndBcastFile(File, fileCharPtr); + std::string fileCharPtrString(fileCharPtr.dataPtr()); + std::istringstream is(fileCharPtrString, std::istringstream::in); + is.exceptions(std::ios_base::failbit | std::ios_base::badbit); + + std::string line, word; + + std::getline(is, line); + + int nlevs; + is >> nlevs; + GotoNextLine(is); + finest_level = nlevs-1; + + std::getline(is, line); + { + std::istringstream lis(line); + lis.exceptions(std::ios_base::failbit | std::ios_base::badbit); + for (auto& istep_lev : istep) { + lis >> word; + istep_lev = std::stoi(word); + } + } + + std::getline(is, line); + { + std::istringstream lis(line); + lis.exceptions(std::ios_base::failbit | std::ios_base::badbit); + for (auto& nsub : nsubsteps) { + lis >> word; + nsub = std::stoi(word); + } + } + + std::getline(is, line); + { + std::istringstream lis(line); + lis.exceptions(std::ios_base::failbit | std::ios_base::badbit); + for (auto& t_new_lev : t_new) { + lis >> word; + t_new_lev = static_cast(std::stod(word)); + } + } + + std::getline(is, line); + { + std::istringstream lis(line); + lis.exceptions(std::ios_base::failbit | std::ios_base::badbit); + for (auto& t_old_lev : t_old) { + lis >> word; + t_old_lev = static_cast(std::stod(word)); + } + } + + std::getline(is, line); + { + std::istringstream lis(line); + lis.exceptions(std::ios_base::failbit | std::ios_base::badbit); + for (auto& dt_lev : dt) { + lis >> word; + dt_lev = static_cast(std::stod(word)); + } + } + + amrex::Real moving_window_x_checkpoint; + is >> moving_window_x_checkpoint; + GotoNextLine(is); + + is >> is_synchronized; + GotoNextLine(is); + + amrex::Vector prob_lo( AMREX_SPACEDIM ); + std::getline(is, line); + { + std::istringstream lis(line); + lis.exceptions(std::ios_base::failbit | std::ios_base::badbit); + for (auto& prob_lo_comp : prob_lo) { + lis >> word; + prob_lo_comp = static_cast(std::stod(word)); + } + } + + amrex::Vector prob_hi( AMREX_SPACEDIM ); + std::getline(is, line); + { + std::istringstream lis(line); + lis.exceptions(std::ios_base::failbit | std::ios_base::badbit); + for (auto& prob_hi_comp : prob_hi) { + lis >> word; + prob_hi_comp = static_cast(std::stod(word)); + } + } + + ResetProbDomain(RealBox(prob_lo.data(),prob_hi.data())); + + for (int lev = 0; lev < nlevs; ++lev) { + BoxArray ba; + ba.readFrom(is); + GotoNextLine(is); + DistributionMapping dm = GetRestartDMap(restart_chkfile, ba, lev); + SetBoxArray(lev, ba); + SetDistributionMap(lev, dm); + AllocLevelData(lev, ba, dm); + } + + mypc->ReadHeader(is); + is >> current_injection_position; + GotoNextLine(is); + + int do_moving_window_before_restart; + is >> do_moving_window_before_restart; + GotoNextLine(is); + + if (do_moving_window_before_restart) { + moving_window_x = moving_window_x_checkpoint; + } + + is >> time_of_last_gal_shift; + GotoNextLine(is); + + + auto & warpx = WarpX::GetInstance(); + for (int idiag = 0; idiag < warpx.GetMultiDiags().GetTotalDiags(); ++idiag) + { + if( warpx.GetMultiDiags().diagstypes(idiag) == DiagTypes::BackTransformed ) + { + auto& diag = warpx.GetMultiDiags().GetDiag(idiag); + if (diag.getnumbuffers() > 0) { + diag.InitDataBeforeRestart(); + for (int i_buffer=0; i_buffer> tlab; + diag.settlab(i_buffer, tlab); + int kindex_hi; + is >> kindex_hi; + diag.set_buffer_k_index_hi(i_buffer, kindex_hi); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + amrex::Real snapshot_lo; + is >> snapshot_lo; + diag.setSnapshotDomainLo(i_buffer, idim, snapshot_lo); + } + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + amrex::Real snapshot_hi; + is >> snapshot_hi; + diag.setSnapshotDomainHi(i_buffer, idim, snapshot_hi); + } + + int flush_counter; + is >> flush_counter; + diag.set_flush_counter(i_buffer, flush_counter); + + int last_valid_Zslice; + is >> last_valid_Zslice; + diag.set_last_valid_Zslice(i_buffer, last_valid_Zslice); + + int snapshot_full_flag; + is >> snapshot_full_flag; + diag.set_snapshot_full(i_buffer, snapshot_full_flag); + + } + diag.InitDataAfterRestart(); + } else { + diag.InitData(); + } + } else { + warpx.GetMultiDiags().GetDiag(idiag).InitData(); + } + } + } + + const int nlevs = finestLevel()+1; + + // Initialize the field data + for (int lev = 0; lev < nlevs; ++lev) + { + for (int i = 0; i < 3; ++i) { + current_fp[lev][i]->setVal(0.0); + Efield_fp[lev][i]->setVal(0.0); + Bfield_fp[lev][i]->setVal(0.0); +#ifdef WARPX_MAG_LLG + Mfield_fp[lev][i]->setVal(0.0); +#endif + } + + if (lev > 0) { + for (int i = 0; i < 3; ++i) { + Efield_aux[lev][i]->setVal(0.0); + Bfield_aux[lev][i]->setVal(0.0); +#ifdef WARPX_MAG_LLG + Mfield_aux[lev][i]->setVal(0.0); +#endif + current_cp[lev][i]->setVal(0.0); + Efield_cp[lev][i]->setVal(0.0); + Bfield_cp[lev][i]->setVal(0.0); +#ifdef WARPX_MAG_LLG + Mfield_cp[lev][i]->setVal(0.0); +#endif + } + } + + VisMF::Read(*Efield_fp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Ex_fp")); + VisMF::Read(*Efield_fp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Ey_fp")); + VisMF::Read(*Efield_fp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Ez_fp")); + + VisMF::Read(*Bfield_fp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Bx_fp")); + VisMF::Read(*Bfield_fp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "By_fp")); + VisMF::Read(*Bfield_fp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Bz_fp")); + +#ifdef WARPX_MAG_LLG + VisMF::Read(*Hfield_fp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hx_fp")); + VisMF::Read(*Hfield_fp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hy_fp")); + VisMF::Read(*Hfield_fp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hz_fp")); + VisMF::Read(*Mfield_fp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Mx_fp")); + VisMF::Read(*Mfield_fp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "My_fp")); + VisMF::Read(*Mfield_fp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Mz_fp")); + VisMF::Read(*H_biasfield_fp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hxbias_fp")); + VisMF::Read(*H_biasfield_fp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hybias_fp")); + VisMF::Read(*H_biasfield_fp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hzbias_fp")); +#endif + if (WarpX::fft_do_time_averaging) + { + VisMF::Read(*Efield_avg_fp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Ex_avg_fp")); + VisMF::Read(*Efield_avg_fp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Ey_avg_fp")); + VisMF::Read(*Efield_avg_fp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Ez_avg_fp")); + + VisMF::Read(*Bfield_avg_fp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Bx_avg_fp")); + VisMF::Read(*Bfield_avg_fp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "By_avg_fp")); + VisMF::Read(*Bfield_avg_fp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Bz_avg_fp")); + } + + if (is_synchronized + || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon + || use_lumped_inductor == 1) + { + VisMF::Read(*current_fp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jx_fp")); + VisMF::Read(*current_fp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jy_fp")); + VisMF::Read(*current_fp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jz_fp")); + } + + if (lev > 0) + { + VisMF::Read(*Efield_cp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Ex_cp")); + VisMF::Read(*Efield_cp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Ey_cp")); + VisMF::Read(*Efield_cp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Ez_cp")); + + VisMF::Read(*Bfield_cp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Bx_cp")); + VisMF::Read(*Bfield_cp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "By_cp")); + VisMF::Read(*Bfield_cp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Bz_cp")); + +#ifdef WARPX_MAG_LLG + VisMF::Read(*Hfield_cp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hx_cp")); + VisMF::Read(*Hfield_cp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hy_cp")); + VisMF::Read(*Hfield_cp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hz_cp")); + + VisMF::Read(*Mfield_cp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Mx_cp")); + VisMF::Read(*Mfield_cp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "My_cp")); + VisMF::Read(*Mfield_cp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Mz_cp")); + + VisMF::Read(*H_biasfield_cp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hxbias_cp")); + VisMF::Read(*H_biasfield_cp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hybias_cp")); + VisMF::Read(*H_biasfield_cp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hzbias_cp")); +#endif + if (WarpX::fft_do_time_averaging) + { + VisMF::Read(*Efield_avg_cp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Ex_avg_cp")); + VisMF::Read(*Efield_avg_cp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Ey_avg_cp")); + VisMF::Read(*Efield_avg_cp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Ez_avg_cp")); + + VisMF::Read(*Bfield_avg_cp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Bx_avg_cp")); + VisMF::Read(*Bfield_avg_cp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "By_avg_cp")); + VisMF::Read(*Bfield_avg_cp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Bz_avg_cp")); + } + + if (is_synchronized + || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon + || use_lumped_inductor == 1) { + VisMF::Read(*current_cp[lev][0], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jx_cp")); + VisMF::Read(*current_cp[lev][1], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jy_cp")); + VisMF::Read(*current_cp[lev][2], + amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jz_cp")); + } + } + } + + InitPML(); + if (do_pml) + { + for (int lev = 0; lev < nlevs; ++lev) { + if (pml[lev]) + pml[lev]->Restart(amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "pml")); +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + if (pml_rz[lev]) + pml_rz[lev]->Restart(amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "pml_rz")); +#endif + } + } + + InitializeEBGridData(maxLevel()); + + // Initialize particles + mypc->AllocData(); + mypc->Restart(restart_chkfile); + +} diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 65a5564f5..a160a356b 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -154,6 +154,10 @@ WarpX::Evolve (int numsteps) m_inductor->EvolveInductorJ(-0.5_rt*dt[0]); // J^(n) to J^(n-1/2) using E^(n) FillBoundaryJ(guard_cells.ng_alloc_EB); } + if (use_josephson_junction) { + m_jj->EvolveJunctionJ(-0.5_rt*dt[0]); // phi^(0) to phi^(-1/2), add JJ J^(-1/2) + FillBoundaryJ(guard_cells.ng_alloc_EB); + } is_synchronized = false; } else { if (electrostatic_solver_id == ElectrostaticSolverAlgo::None) { @@ -431,6 +435,10 @@ WarpX::OneStep_nosub (Real cur_time) m_inductor->EvolveInductorJ(dt[0]); // J^(n-1/2) to J^(n+1/2) using E^(n) FillBoundaryJ(guard_cells.ng_alloc_EB); } + if (use_josephson_junction == 1) { + m_jj->EvolveJunctionJ(dt[0]); // phi^(n-1/2) to phi^(n+1/2), add JJ J^(n+1/2) + FillBoundaryJ(guard_cells.ng_alloc_EB); + } ExecutePythonCallback("afterdeposition"); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp index 1ff8f7f60..cfa6ec973 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp @@ -133,12 +133,37 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( amrex::ignore_unused(edge_lengths); #endif + auto & warpx = WarpX::GetInstance(); + const int lev = 0; + const amrex::GpuArray dx = warpx.Geom(lev).CellSizeArray(); + + const int use_lumped_resistor = warpx.use_lumped_resistor; + const int use_lumped_capacitor = warpx.use_lumped_capacitor; + amrex::MultiFab& sigma_mf = macroscopic_properties->getsigma_mf(); amrex::MultiFab& epsilon_mf = macroscopic_properties->getepsilon_mf(); #ifndef WARPX_MAG_LLG amrex::MultiFab& mu_mf = macroscopic_properties->getmu_mf(); #endif + amrex::MultiFab* lumped_resistor_x_mf = nullptr; + amrex::MultiFab* lumped_resistor_y_mf = nullptr; + amrex::MultiFab* lumped_resistor_z_mf = nullptr; + amrex::MultiFab* lumped_capacitor_x_mf = nullptr; + amrex::MultiFab* lumped_capacitor_y_mf = nullptr; + amrex::MultiFab* lumped_capacitor_z_mf = nullptr; + + if (use_lumped_resistor) { + lumped_resistor_x_mf = macroscopic_properties->get_pointer_lumped_resistor_x(); + lumped_resistor_y_mf = macroscopic_properties->get_pointer_lumped_resistor_y(); + lumped_resistor_z_mf = macroscopic_properties->get_pointer_lumped_resistor_z(); + } + if (use_lumped_capacitor) { + lumped_capacitor_x_mf = macroscopic_properties->get_pointer_lumped_capacitor_x(); + lumped_capacitor_y_mf = macroscopic_properties->get_pointer_lumped_capacitor_y(); + lumped_capacitor_z_mf = macroscopic_properties->get_pointer_lumped_capacitor_z(); + } + // Index type required for calling ablastr::coarsen::sample::Interp to interpolate macroscopic // properties from their respective staggering to the Ex, Ey, Ez locations amrex::GpuArray const& sigma_stag = macroscopic_properties->sigma_IndexType; @@ -180,6 +205,20 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( amrex::Array4 const& mu_arr = mu_mf.array(mfi); #endif + // Lumped element arrays (default-constructed; populated only when feature enabled). + amrex::Array4 resistor_x_arr, resistor_y_arr, resistor_z_arr; + amrex::Array4 capacitor_x_arr, capacitor_y_arr, capacitor_z_arr; + if (use_lumped_resistor) { + resistor_x_arr = lumped_resistor_x_mf->array(mfi); + resistor_y_arr = lumped_resistor_y_mf->array(mfi); + resistor_z_arr = lumped_resistor_z_mf->array(mfi); + } + if (use_lumped_capacitor) { + capacitor_x_arr = lumped_capacitor_x_mf->array(mfi); + capacitor_y_arr = lumped_capacitor_y_mf->array(mfi); + capacitor_z_arr = lumped_capacitor_z_mf->array(mfi); + } + // Extract stencil coefficients Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); int const n_coefs_x = m_stencil_coefs_x.size(); @@ -214,18 +253,25 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( // Skip field push if this cell is fully covered by embedded boundaries if (lx(i, j, k) <= 0) return; #endif - // Interpolate conductivity, sigma, to Ex position on the grid amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, Ex_stag, macro_cr, i, j, k, scomp); - // Interpolated permittivity, epsilon, to Ex position on the grid amrex::Real const epsilon_interp = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag, Ex_stag, macro_cr, i, j, k, scomp); - amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt); - amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); - Ex(i, j, k) = alpha * Ex(i, j, k) - + beta * ( - T_Algo::DownwardDz(Hy, coefs_z, n_coefs_z, i, j, k,0) - + T_Algo::DownwardDy(Hz, coefs_y, n_coefs_y, i, j, k,0) - ) - beta * jx(i, j, k); + // fac1 = (dt/eps)*(sigma + lumped-resistor equivalent conductance) + amrex::Real const fac1 = (dt/epsilon_interp) * + ( sigma_interp + ((use_lumped_resistor && resistor_x_arr(i,j,k)!=0._rt) + ? dx[0] / (dx[1]*dx[2]*resistor_x_arr(i,j,k)) + : 0._rt) ); + // fac2 = lumped-capacitor / natural cell capacitance + amrex::Real const fac2 = (use_lumped_capacitor && capacitor_x_arr(i,j,k)!=0._rt) + ? capacitor_x_arr(i,j,k) * dx[0] / (dx[1]*dx[2]*epsilon_interp) + : 0._rt; + amrex::Real const alpha_c = T_MacroAlgo::alpha_compact(fac1, fac2); + amrex::Real const beta_c = T_MacroAlgo::beta_compact(fac1, fac2); + Ex(i, j, k) = alpha_c * Ex(i, j, k) + + (dt/epsilon_interp) * beta_c * ( - T_Algo::DownwardDz(Hy, coefs_z, n_coefs_z, i, j, k,0) + + T_Algo::DownwardDy(Hz, coefs_y, n_coefs_y, i, j, k,0)) + - (dt/epsilon_interp) * beta_c * jx(i, j, k); }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ @@ -238,18 +284,23 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( if (lx(i, j, k)<=0 || lx(i-1, j, k)<=0 || lz(i, j, k)<=0 || lz(i, j-1, k)<=0) return; #endif #endif - // Interpolate conductivity, sigma, to Ey position on the grid amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, Ey_stag, macro_cr, i, j, k, scomp); - // Interpolated permittivity, epsilon, to Ey position on the grid amrex::Real const epsilon_interp = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag, Ey_stag, macro_cr, i, j, k, scomp); - amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt); - amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); - Ey(i, j, k) = alpha * Ey(i, j, k) - + beta * ( - T_Algo::DownwardDx(Hz, coefs_x, n_coefs_x, i, j, k,0) - + T_Algo::DownwardDz(Hx, coefs_z, n_coefs_z, i, j, k,0) - ) - beta * jy(i, j, k); + amrex::Real const fac1 = (dt/epsilon_interp) * + ( sigma_interp + ((use_lumped_resistor && resistor_y_arr(i,j,k)!=0._rt) + ? dx[1] / (dx[0]*dx[2]*resistor_y_arr(i,j,k)) + : 0._rt) ); + amrex::Real const fac2 = (use_lumped_capacitor && capacitor_y_arr(i,j,k)!=0._rt) + ? capacitor_y_arr(i,j,k) * dx[1] / (dx[0]*dx[2]*epsilon_interp) + : 0._rt; + amrex::Real const alpha_c = T_MacroAlgo::alpha_compact(fac1, fac2); + amrex::Real const beta_c = T_MacroAlgo::beta_compact(fac1, fac2); + Ey(i, j, k) = alpha_c * Ey(i, j, k) + + (dt/epsilon_interp) * beta_c * ( - T_Algo::DownwardDx(Hz, coefs_x, n_coefs_x, i, j, k,0) + + T_Algo::DownwardDz(Hx, coefs_z, n_coefs_z, i, j, k,0)) + - (dt/epsilon_interp) * beta_c * jy(i, j, k); }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ @@ -257,18 +308,23 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( // Skip field push if this cell is fully covered by embedded boundaries if (lz(i,j,k) <= 0) return; #endif - // Interpolate conductivity, sigma, to Ez position on the grid amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, Ez_stag, macro_cr, i, j, k, scomp); - // Interpolated permittivity, epsilon, to Ez position on the grid amrex::Real const epsilon_interp = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag, Ez_stag, macro_cr, i, j, k, scomp); - amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt); - amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); - Ez(i, j, k) = alpha * Ez(i, j, k) - + beta * ( - T_Algo::DownwardDy(Hx, coefs_y, n_coefs_y, i, j, k,0) - + T_Algo::DownwardDx(Hy, coefs_x, n_coefs_x, i, j, k,0) - ) - beta * jz(i, j, k); + amrex::Real const fac1 = (dt/epsilon_interp) * + ( sigma_interp + ((use_lumped_resistor && resistor_z_arr(i,j,k)!=0._rt) + ? dx[2] / (dx[0]*dx[1]*resistor_z_arr(i,j,k)) + : 0._rt) ); + amrex::Real const fac2 = (use_lumped_capacitor && capacitor_z_arr(i,j,k)!=0._rt) + ? capacitor_z_arr(i,j,k) * dx[2] / (dx[0]*dx[1]*epsilon_interp) + : 0._rt; + amrex::Real const alpha_c = T_MacroAlgo::alpha_compact(fac1, fac2); + amrex::Real const beta_c = T_MacroAlgo::beta_compact(fac1, fac2); + Ez(i, j, k) = alpha_c * Ez(i, j, k) + + (dt/epsilon_interp) * beta_c * ( - T_Algo::DownwardDy(Hx, coefs_y, n_coefs_y, i, j, k,0) + + T_Algo::DownwardDx(Hy, coefs_x, n_coefs_x, i, j, k,0)) + - (dt/epsilon_interp) * beta_c * jz(i, j, k); } ); } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H index 325c995e7..6fad655c5 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H @@ -47,6 +47,22 @@ public: amrex::MultiFab& getmu_mf () {return (*m_mu_mf);} amrex::MultiFab * get_pointer_mu () {return m_mu_mf.get();} + /** return MultiFabs holding the per-edge lumped resistor values [Ohm] */ + amrex::MultiFab& getlumped_resistor_x_mf () {return (*m_lumped_resistor_x_mf);} + amrex::MultiFab * get_pointer_lumped_resistor_x () {return m_lumped_resistor_x_mf.get();} + amrex::MultiFab& getlumped_resistor_y_mf () {return (*m_lumped_resistor_y_mf);} + amrex::MultiFab * get_pointer_lumped_resistor_y () {return m_lumped_resistor_y_mf.get();} + amrex::MultiFab& getlumped_resistor_z_mf () {return (*m_lumped_resistor_z_mf);} + amrex::MultiFab * get_pointer_lumped_resistor_z () {return m_lumped_resistor_z_mf.get();} + + /** return MultiFabs holding the per-edge lumped capacitor values [Farad] */ + amrex::MultiFab& getlumped_capacitor_x_mf () {return (*m_lumped_capacitor_x_mf);} + amrex::MultiFab * get_pointer_lumped_capacitor_x () {return m_lumped_capacitor_x_mf.get();} + amrex::MultiFab& getlumped_capacitor_y_mf () {return (*m_lumped_capacitor_y_mf);} + amrex::MultiFab * get_pointer_lumped_capacitor_y () {return m_lumped_capacitor_y_mf.get();} + amrex::MultiFab& getlumped_capacitor_z_mf () {return (*m_lumped_capacitor_z_mf);} + amrex::MultiFab * get_pointer_lumped_capacitor_z () {return m_lumped_capacitor_z_mf.get();} + /** Gpu Vector with index type of coarsening ratio with default value (1,1,1) */ amrex::GpuArray macro_cr_ratio; /** Initializes the Multifabs storing macroscopic properties @@ -86,6 +102,21 @@ public: /** Gpu Vector with index type of the Bz multifab */ amrex::GpuArray Bz_IndexType; + /** Gpu Vector with index type of the lumped resistor multifab */ + amrex::GpuArray lumped_resistor_x_IndexType; + amrex::GpuArray lumped_resistor_y_IndexType; + amrex::GpuArray lumped_resistor_z_IndexType; + + /** Gpu Vector with index type of the lumped capacitor multifab */ + amrex::GpuArray lumped_capacitor_x_IndexType; + amrex::GpuArray lumped_capacitor_y_IndexType; + amrex::GpuArray lumped_capacitor_z_IndexType; + + /** Gpu Vector with index type of jx/jy/jz multifab (edge-centered) */ + amrex::GpuArray jx_IndexType; + amrex::GpuArray jy_IndexType; + amrex::GpuArray jz_IndexType; + /** Stores initialization type for conductivity : constant or parser */ std::string m_sigma_s = "constant"; /** Stores initialization type for permittivity : constant or parser */ @@ -117,6 +148,14 @@ public: std::unique_ptr m_epsilon_parser; std::unique_ptr m_mu_parser; + std::unique_ptr m_lumped_resistor_x_parser; + std::unique_ptr m_lumped_resistor_y_parser; + std::unique_ptr m_lumped_resistor_z_parser; + + std::unique_ptr m_lumped_capacitor_x_parser; + std::unique_ptr m_lumped_capacitor_y_parser; + std::unique_ptr m_lumped_capacitor_z_parser; + #ifdef WARPX_MAG_LLG /** Gpu Vector with index type of the Hx multifab */ amrex::GpuArray Hx_IndexType; @@ -330,6 +369,22 @@ private: std::string m_str_epsilon_function; std::string m_str_mu_function; + /** Multifabs for lumped resistor (per-edge Ohm values) */ + std::unique_ptr m_lumped_resistor_x_mf; + std::unique_ptr m_lumped_resistor_y_mf; + std::unique_ptr m_lumped_resistor_z_mf; + std::string m_str_lumped_resistor_x_function; + std::string m_str_lumped_resistor_y_function; + std::string m_str_lumped_resistor_z_function; + + /** Multifabs for lumped capacitor (per-edge Farad values) */ + std::unique_ptr m_lumped_capacitor_x_mf; + std::unique_ptr m_lumped_capacitor_y_mf; + std::unique_ptr m_lumped_capacitor_z_mf; + std::string m_str_lumped_capacitor_x_function; + std::string m_str_lumped_capacitor_y_function; + std::string m_str_lumped_capacitor_z_function; + }; /** @@ -362,6 +417,24 @@ struct LaxWendroffAlgo { return beta; } + /** Compact form for combined dissipation (fac1) and added capacitance (fac2). + * Used when lumped resistor/capacitor terms are present. */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real alpha_compact (amrex::Real const fac1, + amrex::Real const fac2) { + using namespace amrex; + amrex::Real alpha_compact = (1._rt - 0.5_rt * fac1 + fac2)/(1._rt + 0.5_rt * fac1 + fac2); + return alpha_compact; + } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real beta_compact (amrex::Real const fac1, + amrex::Real const fac2) { + using namespace amrex; + amrex::Real beta_compact = 1._rt / (1._rt + 0.5_rt * fac1 + fac2); + return beta_compact; + } + }; /** @@ -394,6 +467,24 @@ struct BackwardEulerAlgo { return beta; } + /** Compact form for combined dissipation (fac1) and added capacitance (fac2). + * Used when lumped resistor/capacitor terms are present. */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real alpha_compact (amrex::Real const fac1, + amrex::Real const fac2) { + using namespace amrex; + amrex::Real alpha_compact = (1._rt + fac2)/(1._rt + fac1 + fac2); + return alpha_compact; + } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real beta_compact (amrex::Real const fac1, + amrex::Real const fac2) { + using namespace amrex; + amrex::Real beta_compact = 1._rt / (1._rt + fac1 + fac2); + return beta_compact; + } + }; #endif // WARPX_MACROSCOPIC_PROPERTIES_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp index 1a1774cbb..0a22fddfb 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp @@ -164,6 +164,36 @@ MacroscopicProperties::ReadParameters () utils::parser::makeParser(m_str_mu_function,{"x","y","z"})); } + // Lumped resistor R(x,y,z) [Ohm], one parser per edge direction + if (WarpX::use_lumped_resistor == 1) { + utils::parser::Store_parserString(pp_macroscopic, "lumped_resistor_x_function(x,y,z)", m_str_lumped_resistor_x_function); + m_lumped_resistor_x_parser = std::make_unique( + utils::parser::makeParser(m_str_lumped_resistor_x_function,{"x","y","z"})); + + utils::parser::Store_parserString(pp_macroscopic, "lumped_resistor_y_function(x,y,z)", m_str_lumped_resistor_y_function); + m_lumped_resistor_y_parser = std::make_unique( + utils::parser::makeParser(m_str_lumped_resistor_y_function,{"x","y","z"})); + + utils::parser::Store_parserString(pp_macroscopic, "lumped_resistor_z_function(x,y,z)", m_str_lumped_resistor_z_function); + m_lumped_resistor_z_parser = std::make_unique( + utils::parser::makeParser(m_str_lumped_resistor_z_function,{"x","y","z"})); + } + + // Lumped capacitor C(x,y,z) [F], one parser per edge direction + if (WarpX::use_lumped_capacitor == 1) { + utils::parser::Store_parserString(pp_macroscopic, "lumped_capacitor_x_function(x,y,z)", m_str_lumped_capacitor_x_function); + m_lumped_capacitor_x_parser = std::make_unique( + utils::parser::makeParser(m_str_lumped_capacitor_x_function,{"x","y","z"})); + + utils::parser::Store_parserString(pp_macroscopic, "lumped_capacitor_y_function(x,y,z)", m_str_lumped_capacitor_y_function); + m_lumped_capacitor_y_parser = std::make_unique( + utils::parser::makeParser(m_str_lumped_capacitor_y_function,{"x","y","z"})); + + utils::parser::Store_parserString(pp_macroscopic, "lumped_capacitor_z_function(x,y,z)", m_str_lumped_capacitor_z_function); + m_lumped_capacitor_z_parser = std::make_unique( + utils::parser::makeParser(m_str_lumped_capacitor_z_function,{"x","y","z"})); + } + #ifdef WARPX_MAG_LLG auto &warpx = WarpX::GetInstance(); pp_macroscopic.get("mag_Ms_init_style", m_mag_Ms_s); @@ -353,6 +383,42 @@ MacroscopicProperties::InitData () } } + // Lumped elements live on the same staggering as J (edge-centered). + amrex::IntVect jx_stag = warpx.get_pointer_current_fp(lev,0)->ixType().toIntVect(); + amrex::IntVect jy_stag = warpx.get_pointer_current_fp(lev,1)->ixType().toIntVect(); + amrex::IntVect jz_stag = warpx.get_pointer_current_fp(lev,2)->ixType().toIntVect(); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + jx_IndexType[idim] = jx_stag[idim]; + jy_IndexType[idim] = jy_stag[idim]; + jz_IndexType[idim] = jz_stag[idim]; + lumped_resistor_x_IndexType[idim] = jx_stag[idim]; + lumped_resistor_y_IndexType[idim] = jy_stag[idim]; + lumped_resistor_z_IndexType[idim] = jz_stag[idim]; + lumped_capacitor_x_IndexType[idim] = jx_stag[idim]; + lumped_capacitor_y_IndexType[idim] = jy_stag[idim]; + lumped_capacitor_z_IndexType[idim] = jz_stag[idim]; + } + + if (warpx.use_lumped_resistor == 1) { + m_lumped_resistor_x_mf = std::make_unique(amrex::convert(ba,jx_stag), dmap, 1, ng_EB_alloc); + m_lumped_resistor_y_mf = std::make_unique(amrex::convert(ba,jy_stag), dmap, 1, ng_EB_alloc); + m_lumped_resistor_z_mf = std::make_unique(amrex::convert(ba,jz_stag), dmap, 1, ng_EB_alloc); + + InitializeMacroMultiFabUsingParser(m_lumped_resistor_x_mf.get(), m_lumped_resistor_x_parser->compile<3>(), lev); + InitializeMacroMultiFabUsingParser(m_lumped_resistor_y_mf.get(), m_lumped_resistor_y_parser->compile<3>(), lev); + InitializeMacroMultiFabUsingParser(m_lumped_resistor_z_mf.get(), m_lumped_resistor_z_parser->compile<3>(), lev); + } + + if (warpx.use_lumped_capacitor == 1) { + m_lumped_capacitor_x_mf = std::make_unique(amrex::convert(ba,jx_stag), dmap, 1, ng_EB_alloc); + m_lumped_capacitor_y_mf = std::make_unique(amrex::convert(ba,jy_stag), dmap, 1, ng_EB_alloc); + m_lumped_capacitor_z_mf = std::make_unique(amrex::convert(ba,jz_stag), dmap, 1, ng_EB_alloc); + + InitializeMacroMultiFabUsingParser(m_lumped_capacitor_x_mf.get(), m_lumped_capacitor_x_parser->compile<3>(), lev); + InitializeMacroMultiFabUsingParser(m_lumped_capacitor_y_mf.get(), m_lumped_capacitor_y_parser->compile<3>(), lev); + InitializeMacroMultiFabUsingParser(m_lumped_capacitor_z_mf.get(), m_lumped_capacitor_z_parser->compile<3>(), lev); + } + #ifdef WARPX_MAG_LLG // all magnetic macroparameters are stored on faces diff --git a/Source/FieldSolver/LumpedElement/CMakeLists.txt b/Source/FieldSolver/LumpedElement/CMakeLists.txt index a1f411f4a..5d4264a2e 100644 --- a/Source/FieldSolver/LumpedElement/CMakeLists.txt +++ b/Source/FieldSolver/LumpedElement/CMakeLists.txt @@ -1,4 +1,5 @@ target_sources(WarpX PRIVATE Inductor.cpp + JosephsonJunction.cpp ) diff --git a/Source/FieldSolver/LumpedElement/JosephsonJunction.H b/Source/FieldSolver/LumpedElement/JosephsonJunction.H new file mode 100644 index 000000000..61fa6a8ea --- /dev/null +++ b/Source/FieldSolver/LumpedElement/JosephsonJunction.H @@ -0,0 +1,97 @@ +/* + * This file is part of MicroEleX + * + * License: BSD-3-Clause-LBNL + */ + + +#ifndef JOSEPHSON_JUNCTION_H +#define JOSEPHSON_JUNCTION_H + +#include +#include +#include +#include +#include +#include +#include +#include + + +/** \brief Josephson junction (nonlinear inductor). + * + * Sets the supercurrent term J_J = Ic * sin(phi) / A on every edge where + * Ic(x,y,z) != 0. The phase phi obeys the second Josephson relation + * d phi / dt = (2e / hbar) * V, V = E_edge * edge_length + * advanced with leapfrog + * phi^{n+1/2} = phi^{n-1/2} + dt * (2e/hbar) * V^n. + * Because the current MultiFab is a persistent state in this solver (it is not + * re-zeroed each step), the junction updates J by the *difference* between the + * new and previously-added supercurrent so that J holds the instantaneous + * Ic*sin(phi)/A, not a running sum. This also composes with any London / + * inductor current sharing the same edge. + * + * Linear resistive shunt (R) and junction capacitance (C) of an RCSJ model are + * NOT handled in this class — set them via macroscopic.lumped_resistor_* + * and macroscopic.lumped_capacitor_* on the same cell (see MacroscopicProperties). + */ +class JosephsonJunction { + +public: + JosephsonJunction (); + + void ReadParameters (); + void InitData (); + + /** Advance phi by dt using E^n, then add Ic*sin(phi)/A to the J edges. + * Pass -0.5*dt for the leapfrog initialization push, dt for the main loop. */ + void EvolveJunctionJ (amrex::Real dt); + + void InitializeJunctionMultiFabUsingParser (amrex::MultiFab *mf, + amrex::ParserExecutor<3> const& parser, + int lev); + + /** Persist the junction state (phi, Js) across a restart. + * Without these, a restart silently resets phi to 0 while the EM fields are + * restored from the checkpoint, leaving the junction phase inconsistent with + * the fields (the JJ loses its history and the post-restart dynamics are wrong). + * Call WriteCheckpoint when the checkpoint is flushed, and ReadCheckpoint AFTER + * InitData() on restart (InitData allocates and zeroes these MultiFabs). */ + void WriteCheckpoint (const std::string& dir) const; + void ReadCheckpoint (const std::string& dir); + + // Parser string (from input file) + std::string m_str_Ic_x_function; + std::string m_str_Ic_y_function; + std::string m_str_Ic_z_function; + + std::unique_ptr m_Ic_x_parser; + std::unique_ptr m_Ic_y_parser; + std::unique_ptr m_Ic_z_parser; + + /** Critical current Ic [A] per edge. Non-zero entries mark junction edges. */ + std::unique_ptr m_Ic_x_mf; + std::unique_ptr m_Ic_y_mf; + std::unique_ptr m_Ic_z_mf; + + /** Persistent phase state phi [rad], same staggering as J. */ + std::unique_ptr m_phi_x_mf; + std::unique_ptr m_phi_y_mf; + std::unique_ptr m_phi_z_mf; + + /** Supercurrent density Ic*sin(phi)/A [A/m^2] added to J on the previous + * push. J is a persistent state in this solver (never re-zeroed), so we + * subtract this on the next push and add the new value — J then carries the + * instantaneous supercurrent rather than a running sum over time steps. */ + std::unique_ptr m_Js_x_mf; + std::unique_ptr m_Js_y_mf; + std::unique_ptr m_Js_z_mf; + + amrex::GpuArray jx_IndexType; + amrex::GpuArray jy_IndexType; + amrex::GpuArray jz_IndexType; + +}; + + +#endif diff --git a/Source/FieldSolver/LumpedElement/JosephsonJunction.H.orig b/Source/FieldSolver/LumpedElement/JosephsonJunction.H.orig new file mode 100644 index 000000000..55f30fe64 --- /dev/null +++ b/Source/FieldSolver/LumpedElement/JosephsonJunction.H.orig @@ -0,0 +1,75 @@ +/* + * This file is part of MicroEleX + * + * License: BSD-3-Clause-LBNL + */ + + +#ifndef JOSEPHSON_JUNCTION_H +#define JOSEPHSON_JUNCTION_H + +#include +#include +#include +#include +#include +#include +#include +#include + + +/** \brief Josephson junction (nonlinear inductor). + * + * Adds the supercurrent term J_J = Ic * sin(phi) / A on every edge where + * Ic(x,y,z) != 0. The phase phi obeys the second Josephson relation + * d phi / dt = (2e / hbar) * V, V = E_edge * edge_length + * advanced with leapfrog + * phi^{n+1/2} = phi^{n-1/2} + dt * (2e/hbar) * V^n. + * + * Linear resistive shunt (R) and junction capacitance (C) of an RCSJ model are + * NOT handled in this class — set them via macroscopic.lumped_resistor_* + * and macroscopic.lumped_capacitor_* on the same cell (see MacroscopicProperties). + */ +class JosephsonJunction { + +public: + JosephsonJunction (); + + void ReadParameters (); + void InitData (); + + /** Advance phi by dt using E^n, then add Ic*sin(phi)/A to the J edges. + * Pass -0.5*dt for the leapfrog initialization push, dt for the main loop. */ + void EvolveJunctionJ (amrex::Real dt); + + void InitializeJunctionMultiFabUsingParser (amrex::MultiFab *mf, + amrex::ParserExecutor<3> const& parser, + int lev); + + // Parser string (from input file) + std::string m_str_Ic_x_function; + std::string m_str_Ic_y_function; + std::string m_str_Ic_z_function; + + std::unique_ptr m_Ic_x_parser; + std::unique_ptr m_Ic_y_parser; + std::unique_ptr m_Ic_z_parser; + + /** Critical current Ic [A] per edge. Non-zero entries mark junction edges. */ + std::unique_ptr m_Ic_x_mf; + std::unique_ptr m_Ic_y_mf; + std::unique_ptr m_Ic_z_mf; + + /** Persistent phase state phi [rad], same staggering as J. */ + std::unique_ptr m_phi_x_mf; + std::unique_ptr m_phi_y_mf; + std::unique_ptr m_phi_z_mf; + + amrex::GpuArray jx_IndexType; + amrex::GpuArray jy_IndexType; + amrex::GpuArray jz_IndexType; + +}; + + +#endif diff --git a/Source/FieldSolver/LumpedElement/JosephsonJunction.cpp b/Source/FieldSolver/LumpedElement/JosephsonJunction.cpp new file mode 100644 index 000000000..7cc395f6f --- /dev/null +++ b/Source/FieldSolver/LumpedElement/JosephsonJunction.cpp @@ -0,0 +1,267 @@ +#include "JosephsonJunction.H" +#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" +#include "Utils/WarpXUtil.H" +#include "Utils/WarpXConst.H" +#include "WarpX.H" +#include +#include "Utils/Parser/IntervalsParser.H" +#include "Utils/Parser/ParserUtils.H" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +JosephsonJunction::JosephsonJunction () +{ + amrex::Print() << " JosephsonJunction (nonlinear inductor) class is constructed \n"; + ReadParameters(); +} + +void +JosephsonJunction::ReadParameters () +{ + amrex::ParmParse pp_jj("josephson"); + + utils::parser::Store_parserString(pp_jj, "Ic_x_function(x,y,z)", m_str_Ic_x_function); + m_Ic_x_parser = std::make_unique( + utils::parser::makeParser(m_str_Ic_x_function, {"x", "y", "z"})); + + utils::parser::Store_parserString(pp_jj, "Ic_y_function(x,y,z)", m_str_Ic_y_function); + m_Ic_y_parser = std::make_unique( + utils::parser::makeParser(m_str_Ic_y_function, {"x", "y", "z"})); + + utils::parser::Store_parserString(pp_jj, "Ic_z_function(x,y,z)", m_str_Ic_z_function); + m_Ic_z_parser = std::make_unique( + utils::parser::makeParser(m_str_Ic_z_function, {"x", "y", "z"})); +} + +#if( AMREX_SPACEDIM == 3) +void +JosephsonJunction::InitData() +{ + auto& warpx = WarpX::GetInstance(); + + const int lev = 0; + amrex::BoxArray ba = warpx.boxArray(lev); + amrex::DistributionMapping dmap = warpx.DistributionMap(lev); + const amrex::IntVect ng_EB_alloc = warpx.getngEB(); + + amrex::IntVect jx_stag = warpx.get_pointer_current_fp(lev,0)->ixType().toIntVect(); + amrex::IntVect jy_stag = warpx.get_pointer_current_fp(lev,1)->ixType().toIntVect(); + amrex::IntVect jz_stag = warpx.get_pointer_current_fp(lev,2)->ixType().toIntVect(); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + jx_IndexType[idim] = jx_stag[idim]; + jy_IndexType[idim] = jy_stag[idim]; + jz_IndexType[idim] = jz_stag[idim]; + } + + const int ncomps = 1; + m_Ic_x_mf = std::make_unique(amrex::convert(ba,jx_stag), dmap, ncomps, ng_EB_alloc); + m_Ic_y_mf = std::make_unique(amrex::convert(ba,jy_stag), dmap, ncomps, ng_EB_alloc); + m_Ic_z_mf = std::make_unique(amrex::convert(ba,jz_stag), dmap, ncomps, ng_EB_alloc); + m_phi_x_mf = std::make_unique(amrex::convert(ba,jx_stag), dmap, ncomps, ng_EB_alloc); + m_phi_y_mf = std::make_unique(amrex::convert(ba,jy_stag), dmap, ncomps, ng_EB_alloc); + m_phi_z_mf = std::make_unique(amrex::convert(ba,jz_stag), dmap, ncomps, ng_EB_alloc); + m_Js_x_mf = std::make_unique(amrex::convert(ba,jx_stag), dmap, ncomps, ng_EB_alloc); + m_Js_y_mf = std::make_unique(amrex::convert(ba,jy_stag), dmap, ncomps, ng_EB_alloc); + m_Js_z_mf = std::make_unique(amrex::convert(ba,jz_stag), dmap, ncomps, ng_EB_alloc); + + InitializeJunctionMultiFabUsingParser(m_Ic_x_mf.get(), m_Ic_x_parser->compile<3>(), lev); + InitializeJunctionMultiFabUsingParser(m_Ic_y_mf.get(), m_Ic_y_parser->compile<3>(), lev); + InitializeJunctionMultiFabUsingParser(m_Ic_z_mf.get(), m_Ic_z_parser->compile<3>(), lev); + + m_phi_x_mf->setVal(0.0); + m_phi_y_mf->setVal(0.0); + m_phi_z_mf->setVal(0.0); + m_Js_x_mf->setVal(0.0); + m_Js_y_mf->setVal(0.0); + m_Js_z_mf->setVal(0.0); +} + +void +JosephsonJunction::EvolveJunctionJ (amrex::Real dt) +{ + using namespace amrex::literals; + + amrex::Print() << " evolve Josephson junction: advance phi, add Ic*sin(phi) to J\n"; + auto & warpx = WarpX::GetInstance(); + const int lev = 0; + + const amrex::GpuArray dx = warpx.Geom(lev).CellSizeArray(); + + // 2 e / hbar = Josephson constant ~ 3.0394e15 rad/(V*s) + constexpr amrex::Real two_e_over_hbar = + 2.0_rt * PhysConst::q_e / PhysConst::hbar; + + amrex::MultiFab * jx = warpx.get_pointer_current_fp(lev, 0); + amrex::MultiFab * jy = warpx.get_pointer_current_fp(lev, 1); + amrex::MultiFab * jz = warpx.get_pointer_current_fp(lev, 2); + + amrex::MultiFab * Ex = warpx.get_pointer_Efield_fp(lev, 0); + amrex::MultiFab * Ey = warpx.get_pointer_Efield_fp(lev, 1); + amrex::MultiFab * Ez = warpx.get_pointer_Efield_fp(lev, 2); + + for (amrex::MFIter mfi(*jx, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + amrex::Array4 const& jx_arr = jx->array(mfi); + amrex::Array4 const& jy_arr = jy->array(mfi); + amrex::Array4 const& jz_arr = jz->array(mfi); + amrex::Array4 const& Ex_arr = Ex->array(mfi); + amrex::Array4 const& Ey_arr = Ey->array(mfi); + amrex::Array4 const& Ez_arr = Ez->array(mfi); + amrex::Array4 const& Ic_x_arr = m_Ic_x_mf->array(mfi); + amrex::Array4 const& Ic_y_arr = m_Ic_y_mf->array(mfi); + amrex::Array4 const& Ic_z_arr = m_Ic_z_mf->array(mfi); + amrex::Array4 const& phi_x_arr = m_phi_x_mf->array(mfi); + amrex::Array4 const& phi_y_arr = m_phi_y_mf->array(mfi); + amrex::Array4 const& phi_z_arr = m_phi_z_mf->array(mfi); + amrex::Array4 const& Js_x_arr = m_Js_x_mf->array(mfi); + amrex::Array4 const& Js_y_arr = m_Js_y_mf->array(mfi); + amrex::Array4 const& Js_z_arr = m_Js_z_mf->array(mfi); + amrex::Box const& tjx = mfi.tilebox(jx->ixType().toIntVect()); + amrex::Box const& tjy = mfi.tilebox(jy->ixType().toIntVect()); + amrex::Box const& tjz = mfi.tilebox(jz->ixType().toIntVect()); + + amrex::ParallelFor(tjx, tjy, tjz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (Ic_x_arr(i,j,k) != 0.) { + const amrex::Real V = Ex_arr(i,j,k) * dx[0]; + phi_x_arr(i,j,k) += dt * two_e_over_hbar * V; + const amrex::Real A = dx[1] * dx[2]; + const amrex::Real Js_new = Ic_x_arr(i,j,k) * std::sin(phi_x_arr(i,j,k)) / A; + jx_arr(i,j,k) += Js_new - Js_x_arr(i,j,k); + Js_x_arr(i,j,k) = Js_new; + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (Ic_y_arr(i,j,k) != 0.) { + const amrex::Real V = Ey_arr(i,j,k) * dx[1]; + phi_y_arr(i,j,k) += dt * two_e_over_hbar * V; + const amrex::Real A = dx[0] * dx[2]; + const amrex::Real Js_new = Ic_y_arr(i,j,k) * std::sin(phi_y_arr(i,j,k)) / A; + jy_arr(i,j,k) += Js_new - Js_y_arr(i,j,k); + Js_y_arr(i,j,k) = Js_new; + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (Ic_z_arr(i,j,k) != 0.) { + const amrex::Real V = Ez_arr(i,j,k) * dx[2]; + phi_z_arr(i,j,k) += dt * two_e_over_hbar * V; + const amrex::Real A = dx[0] * dx[1]; + const amrex::Real Js_new = Ic_z_arr(i,j,k) * std::sin(phi_z_arr(i,j,k)) / A; + jz_arr(i,j,k) += Js_new - Js_z_arr(i,j,k); + Js_z_arr(i,j,k) = Js_new; + } + }); + } +} + +void +JosephsonJunction::InitializeJunctionMultiFabUsingParser (amrex::MultiFab *mf, + amrex::ParserExecutor<3> const& parser, + const int lev) +{ + using namespace amrex::literals; + + WarpX& warpx = WarpX::GetInstance(); + const amrex::GpuArray dx = warpx.Geom(lev).CellSizeArray(); + const amrex::RealBox& real_box = warpx.Geom(lev).ProbDomain(); + amrex::IntVect iv = mf->ixType().toIntVect(); + for (amrex::MFIter mfi(*mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const amrex::Box& tb = mfi.tilebox(iv, mf->nGrowVect()); + amrex::Array4 const& fab = mf->array(mfi); + amrex::ParallelFor(tb, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real fac_x = (1._rt - iv[0]) * dx[0] * 0.5_rt; + amrex::Real x = i * dx[0] + real_box.lo(0) + fac_x; + amrex::Real fac_y = (1._rt - iv[1]) * dx[1] * 0.5_rt; + amrex::Real y = j * dx[1] + real_box.lo(1) + fac_y; + amrex::Real fac_z = (1._rt - iv[2]) * dx[2] * 0.5_rt; + amrex::Real z = k * dx[2] + real_box.lo(2) + fac_z; + fab(i,j,k) = parser(x,y,z); + }); + } +} + + +#else +void +JosephsonJunction::InitData() +{ + amrex::Abort("JosephsonJunction only works with 3D"); +} + +void +JosephsonJunction::EvolveJunctionJ (amrex::Real) +{ + amrex::Abort("JosephsonJunction only works with 3D"); +} + +void +JosephsonJunction::InitializeJunctionMultiFabUsingParser (amrex::MultiFab *, + amrex::ParserExecutor<3> const&, + const int) +{ + amrex::Abort("JosephsonJunction only works with 3D"); +} +#endif + + +// --------------------------------------------------------------------------- +// Checkpoint the persistent junction state (phi, Js). +// phi is the junction's only real state; Js is the supercurrent this class +// last added to J (needed by the delta update in EvolveJunctionJ). Both are +// dimension-agnostic MultiFab I/O, so they live outside the 3D-only block. +// --------------------------------------------------------------------------- +void +JosephsonJunction::WriteCheckpoint (const std::string& dir) const +{ + if (!m_phi_x_mf) { return; } // junction never initialized + const int lev = 0; + const std::string lp{"Level_"}; + amrex::VisMF::Write(*m_phi_x_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_x")); + amrex::VisMF::Write(*m_phi_y_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_y")); + amrex::VisMF::Write(*m_phi_z_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_z")); + amrex::VisMF::Write(*m_Js_x_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_Js_x")); + amrex::VisMF::Write(*m_Js_y_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_Js_y")); + amrex::VisMF::Write(*m_Js_z_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_Js_z")); +} + +void +JosephsonJunction::ReadCheckpoint (const std::string& dir) +{ + if (!m_phi_x_mf) { return; } + const int lev = 0; + const std::string lp{"Level_"}; + // Checkpoints written before the junction state was saved have no jj_* files. + // Warn and keep the freshly-zeroed phi rather than aborting -- but the restart + // is then NOT physically continuous (phi loses its history). + if (!amrex::FileExists(amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_x") + "_H")) { + amrex::Print() << "WARNING: checkpoint '" << dir << "' contains no Josephson junction " + << "state (jj_phi_*). phi restarts from 0 while the fields are restored, " + << "so the junction is inconsistent with the fields -- rerun from scratch " + << "or use a checkpoint written with JJ state.\n"; + return; + } + amrex::VisMF::Read(*m_phi_x_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_x")); + amrex::VisMF::Read(*m_phi_y_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_y")); + amrex::VisMF::Read(*m_phi_z_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_z")); + amrex::VisMF::Read(*m_Js_x_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_Js_x")); + amrex::VisMF::Read(*m_Js_y_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_Js_y")); + amrex::VisMF::Read(*m_Js_z_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_Js_z")); + amrex::Print() << " Josephson junction state (phi, Js) restored from checkpoint\n"; +} diff --git a/Source/FieldSolver/LumpedElement/JosephsonJunction.cpp.orig b/Source/FieldSolver/LumpedElement/JosephsonJunction.cpp.orig new file mode 100644 index 000000000..a275bec46 --- /dev/null +++ b/Source/FieldSolver/LumpedElement/JosephsonJunction.cpp.orig @@ -0,0 +1,203 @@ +#include "JosephsonJunction.H" +#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" +#include "Utils/WarpXUtil.H" +#include "Utils/WarpXConst.H" +#include "WarpX.H" +#include +#include "Utils/Parser/IntervalsParser.H" +#include "Utils/Parser/ParserUtils.H" +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +JosephsonJunction::JosephsonJunction () +{ + amrex::Print() << " JosephsonJunction (nonlinear inductor) class is constructed \n"; + ReadParameters(); +} + +void +JosephsonJunction::ReadParameters () +{ + amrex::ParmParse pp_jj("josephson"); + + utils::parser::Store_parserString(pp_jj, "Ic_x_function(x,y,z)", m_str_Ic_x_function); + m_Ic_x_parser = std::make_unique( + utils::parser::makeParser(m_str_Ic_x_function, {"x", "y", "z"})); + + utils::parser::Store_parserString(pp_jj, "Ic_y_function(x,y,z)", m_str_Ic_y_function); + m_Ic_y_parser = std::make_unique( + utils::parser::makeParser(m_str_Ic_y_function, {"x", "y", "z"})); + + utils::parser::Store_parserString(pp_jj, "Ic_z_function(x,y,z)", m_str_Ic_z_function); + m_Ic_z_parser = std::make_unique( + utils::parser::makeParser(m_str_Ic_z_function, {"x", "y", "z"})); +} + +#if( AMREX_SPACEDIM == 3) +void +JosephsonJunction::InitData() +{ + auto& warpx = WarpX::GetInstance(); + + const int lev = 0; + amrex::BoxArray ba = warpx.boxArray(lev); + amrex::DistributionMapping dmap = warpx.DistributionMap(lev); + const amrex::IntVect ng_EB_alloc = warpx.getngEB(); + + amrex::IntVect jx_stag = warpx.get_pointer_current_fp(lev,0)->ixType().toIntVect(); + amrex::IntVect jy_stag = warpx.get_pointer_current_fp(lev,1)->ixType().toIntVect(); + amrex::IntVect jz_stag = warpx.get_pointer_current_fp(lev,2)->ixType().toIntVect(); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + jx_IndexType[idim] = jx_stag[idim]; + jy_IndexType[idim] = jy_stag[idim]; + jz_IndexType[idim] = jz_stag[idim]; + } + + const int ncomps = 1; + m_Ic_x_mf = std::make_unique(amrex::convert(ba,jx_stag), dmap, ncomps, ng_EB_alloc); + m_Ic_y_mf = std::make_unique(amrex::convert(ba,jy_stag), dmap, ncomps, ng_EB_alloc); + m_Ic_z_mf = std::make_unique(amrex::convert(ba,jz_stag), dmap, ncomps, ng_EB_alloc); + m_phi_x_mf = std::make_unique(amrex::convert(ba,jx_stag), dmap, ncomps, ng_EB_alloc); + m_phi_y_mf = std::make_unique(amrex::convert(ba,jy_stag), dmap, ncomps, ng_EB_alloc); + m_phi_z_mf = std::make_unique(amrex::convert(ba,jz_stag), dmap, ncomps, ng_EB_alloc); + + InitializeJunctionMultiFabUsingParser(m_Ic_x_mf.get(), m_Ic_x_parser->compile<3>(), lev); + InitializeJunctionMultiFabUsingParser(m_Ic_y_mf.get(), m_Ic_y_parser->compile<3>(), lev); + InitializeJunctionMultiFabUsingParser(m_Ic_z_mf.get(), m_Ic_z_parser->compile<3>(), lev); + + m_phi_x_mf->setVal(0.0); + m_phi_y_mf->setVal(0.0); + m_phi_z_mf->setVal(0.0); +} + +void +JosephsonJunction::EvolveJunctionJ (amrex::Real dt) +{ + using namespace amrex::literals; + + amrex::Print() << " evolve Josephson junction: advance phi, add Ic*sin(phi) to J\n"; + auto & warpx = WarpX::GetInstance(); + const int lev = 0; + + const amrex::GpuArray dx = warpx.Geom(lev).CellSizeArray(); + + // 2 e / hbar = Josephson constant ~ 3.0394e15 rad/(V*s) + constexpr amrex::Real two_e_over_hbar = + 2.0_rt * PhysConst::q_e / PhysConst::hbar; + + amrex::MultiFab * jx = warpx.get_pointer_current_fp(lev, 0); + amrex::MultiFab * jy = warpx.get_pointer_current_fp(lev, 1); + amrex::MultiFab * jz = warpx.get_pointer_current_fp(lev, 2); + + amrex::MultiFab * Ex = warpx.get_pointer_Efield_fp(lev, 0); + amrex::MultiFab * Ey = warpx.get_pointer_Efield_fp(lev, 1); + amrex::MultiFab * Ez = warpx.get_pointer_Efield_fp(lev, 2); + + for (amrex::MFIter mfi(*jx, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + amrex::Array4 const& jx_arr = jx->array(mfi); + amrex::Array4 const& jy_arr = jy->array(mfi); + amrex::Array4 const& jz_arr = jz->array(mfi); + amrex::Array4 const& Ex_arr = Ex->array(mfi); + amrex::Array4 const& Ey_arr = Ey->array(mfi); + amrex::Array4 const& Ez_arr = Ez->array(mfi); + amrex::Array4 const& Ic_x_arr = m_Ic_x_mf->array(mfi); + amrex::Array4 const& Ic_y_arr = m_Ic_y_mf->array(mfi); + amrex::Array4 const& Ic_z_arr = m_Ic_z_mf->array(mfi); + amrex::Array4 const& phi_x_arr = m_phi_x_mf->array(mfi); + amrex::Array4 const& phi_y_arr = m_phi_y_mf->array(mfi); + amrex::Array4 const& phi_z_arr = m_phi_z_mf->array(mfi); + amrex::Box const& tjx = mfi.tilebox(jx->ixType().toIntVect()); + amrex::Box const& tjy = mfi.tilebox(jy->ixType().toIntVect()); + amrex::Box const& tjz = mfi.tilebox(jz->ixType().toIntVect()); + + amrex::ParallelFor(tjx, tjy, tjz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (Ic_x_arr(i,j,k) != 0.) { + const amrex::Real V = Ex_arr(i,j,k) * dx[0]; + phi_x_arr(i,j,k) += dt * two_e_over_hbar * V; + const amrex::Real A = dx[1] * dx[2]; + jx_arr(i,j,k) += Ic_x_arr(i,j,k) * std::sin(phi_x_arr(i,j,k)) / A; + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (Ic_y_arr(i,j,k) != 0.) { + const amrex::Real V = Ey_arr(i,j,k) * dx[1]; + phi_y_arr(i,j,k) += dt * two_e_over_hbar * V; + const amrex::Real A = dx[0] * dx[2]; + jy_arr(i,j,k) += Ic_y_arr(i,j,k) * std::sin(phi_y_arr(i,j,k)) / A; + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (Ic_z_arr(i,j,k) != 0.) { + const amrex::Real V = Ez_arr(i,j,k) * dx[2]; + phi_z_arr(i,j,k) += dt * two_e_over_hbar * V; + const amrex::Real A = dx[0] * dx[1]; + jz_arr(i,j,k) += Ic_z_arr(i,j,k) * std::sin(phi_z_arr(i,j,k)) / A; + } + }); + } +} + +void +JosephsonJunction::InitializeJunctionMultiFabUsingParser (amrex::MultiFab *mf, + amrex::ParserExecutor<3> const& parser, + const int lev) +{ + using namespace amrex::literals; + + WarpX& warpx = WarpX::GetInstance(); + const amrex::GpuArray dx = warpx.Geom(lev).CellSizeArray(); + const amrex::RealBox& real_box = warpx.Geom(lev).ProbDomain(); + amrex::IntVect iv = mf->ixType().toIntVect(); + for (amrex::MFIter mfi(*mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const amrex::Box& tb = mfi.tilebox(iv, mf->nGrowVect()); + amrex::Array4 const& fab = mf->array(mfi); + amrex::ParallelFor(tb, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real fac_x = (1._rt - iv[0]) * dx[0] * 0.5_rt; + amrex::Real x = i * dx[0] + real_box.lo(0) + fac_x; + amrex::Real fac_y = (1._rt - iv[1]) * dx[1] * 0.5_rt; + amrex::Real y = j * dx[1] + real_box.lo(1) + fac_y; + amrex::Real fac_z = (1._rt - iv[2]) * dx[2] * 0.5_rt; + amrex::Real z = k * dx[2] + real_box.lo(2) + fac_z; + fab(i,j,k) = parser(x,y,z); + }); + } +} + + +#else +void +JosephsonJunction::InitData() +{ + amrex::Abort("JosephsonJunction only works with 3D"); +} + +void +JosephsonJunction::EvolveJunctionJ (amrex::Real) +{ + amrex::Abort("JosephsonJunction only works with 3D"); +} + +void +JosephsonJunction::InitializeJunctionMultiFabUsingParser (amrex::MultiFab *, + amrex::ParserExecutor<3> const&, + const int) +{ + amrex::Abort("JosephsonJunction only works with 3D"); +} +#endif diff --git a/Source/FieldSolver/LumpedElement/Make.package b/Source/FieldSolver/LumpedElement/Make.package index a45b200ca..087c0681a 100644 --- a/Source/FieldSolver/LumpedElement/Make.package +++ b/Source/FieldSolver/LumpedElement/Make.package @@ -1,3 +1,3 @@ CEXE_sources += Inductor.cpp - +CEXE_sources += JosephsonJunction.cpp VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/LumpedElement diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 6a32c40b8..f7d2ab1f5 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -397,6 +397,16 @@ WarpX::InitData () if (WarpX::em_solver_medium==1) { m_macroscopic_properties->InitData(); } + // Lumped-element InitData must run BEFORE diagnostics InitData so that + // their multifabs exist when CellCenterFunctor captures the pointer. + if (use_lumped_inductor) { + amrex::Print() << " calling inductor (early) \n"; + m_inductor->InitData(); + } + if (use_josephson_junction) { + amrex::Print() << " calling Josephson junction (early) \n"; + m_jj->InitData(); + } InitDiagnostics(); } else @@ -407,6 +417,19 @@ WarpX::InitData () m_macroscopic_properties->InitData(); } PostRestart(); + // Same ordering constraint as the fresh-start path. + if (use_lumped_inductor) { + amrex::Print() << " calling inductor (early, restart) \n"; + m_inductor->InitData(); + } + if (use_josephson_junction) { + amrex::Print() << " calling Josephson junction (early, restart) \n"; + m_jj->InitData(); + // InitData() allocates and ZEROES phi/Js. Restore them from the checkpoint + // so the junction continues its phase; otherwise phi restarts at 0 while the + // fields are restored, which silently corrupts the post-restart dynamics. + m_jj->ReadCheckpoint(restart_chkfile); + } reduced_diags->InitData(); multi_diags->InitData(); } @@ -427,10 +450,9 @@ WarpX::InitData () m_london->InitData(); } - if (use_lumped_inductor) { - amrex::Print() << " calling inductor \n"; - m_inductor->InitData(); - } + // NOTE: Inductor and JosephsonJunction InitData were moved earlier + // (into the if/else branches above) so that diagnostics InitData sees + // their multifabs already allocated. Do not re-call here. if (ParallelDescriptor::IOProcessor()) { std::cout << "\nGrids Summary:\n"; diff --git a/Source/Initialization/WarpXInitData.cpp.orig b/Source/Initialization/WarpXInitData.cpp.orig new file mode 100644 index 000000000..4a792717a --- /dev/null +++ b/Source/Initialization/WarpXInitData.cpp.orig @@ -0,0 +1,1866 @@ +/* Copyright 2019-2020 Andrew Myers, Ann Almgren, Aurore Blelly + * Axel Huebl, Burlen Loring, Maxence Thevenet + * Michael Rowan, Remi Lehe, Revathi Jambunathan + * Weiqun Zhang + * + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "WarpX.H" + +#include "BoundaryConditions/PML.H" +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) +# include "BoundaryConditions/PML_RZ.H" +#endif +#include "Diagnostics/MultiDiagnostics.H" +#include "Diagnostics/ReducedDiags/MultiReducedDiags.H" +#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" +#include "Filter/BilinearFilter.H" +#include "Filter/NCIGodfreyFilter.H" +#include "Particles/MultiParticleContainer.H" +#include "Utils/Algorithms/LinearInterpolation.H" +#include "Utils/Logo/GetLogo.H" +#include "Utils/MPIInitHelpers.H" +#include "Utils/Parser/ParserUtils.H" +#include "Utils/TextMsg.H" +#include "Utils/WarpXAlgorithmSelection.H" +#include "Utils/WarpXConst.H" +#include "Utils/WarpXProfilerWrapper.H" +#include "Utils/WarpXUtil.H" +#include "Python/WarpX_py.H" + +#include +#include +#include + +#include +#include +#ifdef AMREX_USE_SENSEI_INSITU +# include +#endif +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#ifdef WARPX_USE_OPENPMD +# include +#endif + +#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" + +using namespace amrex; + +void +WarpX::PostProcessBaseGrids (BoxArray& ba0) const +{ + if (numprocs != 0) { + const Box& dom = Geom(0).Domain(); + const IntVect& domlo = dom.smallEnd(); + const IntVect& domlen = dom.size(); + const IntVect sz = domlen / numprocs; + const IntVect extra = domlen - sz*numprocs; + BoxList bl; +#if defined(WARPX_DIM_3D) + for (int k = 0; k < numprocs[2]; ++k) { + // The first extra[2] blocks get one extra cell with a total of + // sz[2]+1. The rest get sz[2] cells. The docomposition in y + // and x directions are similar. + int klo = (k < extra[2]) ? k*(sz[2]+1) : (k*sz[2]+extra[2]); + int khi = (k < extra[2]) ? klo+(sz[2]+1)-1 : klo+sz[2]-1; + klo += domlo[2]; + khi += domlo[2]; +#endif +#if (AMREX_SPACEDIM >= 2) + for (int j = 0; j < numprocs[1]; ++j) { + int jlo = (j < extra[1]) ? j*(sz[1]+1) : (j*sz[1]+extra[1]); + int jhi = (j < extra[1]) ? jlo+(sz[1]+1)-1 : jlo+sz[1]-1; + jlo += domlo[1]; + jhi += domlo[1]; +#endif + for (int i = 0; i < numprocs[0]; ++i) { + int ilo = (i < extra[0]) ? i*(sz[0]+1) : (i*sz[0]+extra[0]); + int ihi = (i < extra[0]) ? ilo+(sz[0]+1)-1 : ilo+sz[0]-1; + ilo += domlo[0]; + ihi += domlo[0]; + bl.push_back(Box(IntVect(AMREX_D_DECL(ilo,jlo,klo)), + IntVect(AMREX_D_DECL(ihi,jhi,khi)))); + AMREX_D_TERM(},},}) + ba0 = BoxArray(std::move(bl)); + } +} + +void +WarpX::PrintMainPICparameters () +{ + amrex::Print() << "-------------------------------------------------------------------------------\n"; + amrex::Print() << "--------------------------- MAIN EM PIC PARAMETERS ----------------------------\n"; + amrex::Print() << "-------------------------------------------------------------------------------\n"; + + // print warpx build information + if constexpr (std::is_same::value) { + amrex::Print() << "Precision: | SINGLE" << "\n"; + } + else { + amrex::Print() << "Precision: | DOUBLE" << "\n"; + } + if constexpr (std::is_same::value) { + amrex::Print() << "Particle precision: | SINGLE" << "\n"; + } + else { + amrex::Print() << "Particle precision: | DOUBLE" << "\n"; + } + + // Print geometry dimensionality + amrex::ParmParse pp_geometry("geometry"); + std::string dims; + pp_geometry.query( "dims", dims ); + if (dims=="1") { + amrex::Print() << "Geometry: | 1D (Z)" << "\n"; + } + else if (dims=="2") { + amrex::Print() << "Geometry: | 2D (XZ)" << "\n"; + } + else if (dims=="3") { + amrex::Print() << "Geometry: | 3D (XYZ)" << "\n"; + } + else if (dims=="RZ") { + amrex::Print() << "Geometry: | 2D (RZ)" << "\n"; + } + + #ifdef WARPX_DIM_RZ + amrex::Print() << " | - n_rz_azimuthal_modes = " << + WarpX::n_rz_azimuthal_modes << "\n"; + #endif // WARPX_USE_RZ + //Print solver's operation mode (e.g., EM or electrostatic) + if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame) { + amrex::Print() << "Operation mode: | Electrostatic" << "\n"; + amrex::Print() << " | - laboratory frame" << "\n"; + } + else if (electrostatic_solver_id == ElectrostaticSolverAlgo::Relativistic){ + amrex::Print() << "Operation mode: | Electrostatic" << "\n"; + amrex::Print() << " | - relativistic" << "\n"; + } + else if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic){ + amrex::Print() << "Operation mode: | Electrostatic" << "\n"; + amrex::Print() << " | - laboratory frame, electrostatic + magnetostatic" << "\n"; + } + else{ + amrex::Print() << "Operation mode: | Electromagnetic" << "\n"; + } + if (em_solver_medium == MediumForEM::Vacuum ){ + amrex::Print() << " | - vacuum" << "\n"; + } + else if (em_solver_medium == MediumForEM::Macroscopic ){ + amrex::Print() << " | - macroscopic" << "\n"; + } + if ( (em_solver_medium == MediumForEM::Macroscopic) && + (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff)){ + amrex::Print() << " | - Lax-Wendroff algorithm\n"; + } + else if ((em_solver_medium == MediumForEM::Macroscopic) && + (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler)){ + amrex::Print() << " | - Backward Euler algorithm\n"; + } + amrex::Print() << "-------------------------------------------------------------------------------\n"; + // Print type of current deposition + if (current_deposition_algo == CurrentDepositionAlgo::Direct){ + amrex::Print() << "Current Deposition: | direct \n"; + } + else if (current_deposition_algo == CurrentDepositionAlgo::Vay){ + amrex::Print() << "Current Deposition: | Vay \n"; + } + else if (current_deposition_algo == CurrentDepositionAlgo::Esirkepov){ + amrex::Print() << "Current Deposition: | Esirkepov \n"; + } + // Print type of particle pusher + if (particle_pusher_algo == ParticlePusherAlgo::Vay){ + amrex::Print() << "Particle Pusher: | Vay \n"; + } + else if (particle_pusher_algo == ParticlePusherAlgo::HigueraCary){ + amrex::Print() << "Particle Pusher: | Higuera-Cary \n"; + } + else if (particle_pusher_algo == ParticlePusherAlgo::Boris){ + amrex::Print() << "Particle Pusher: | Boris \n"; + } + // Print type of charge deposition + if (charge_deposition_algo == ChargeDepositionAlgo::Standard){ + amrex::Print() << "Charge Deposition: | standard \n"; + } + // Print field gathering algorithm + if (field_gathering_algo == GatheringAlgo::MomentumConserving){ + amrex::Print() << "Field Gathering: | momentum-conserving \n"; + } + else{ + amrex::Print() << "Field Gathering: | energy-conserving \n"; + } + // Print particle's shape factors + amrex::Print() << "Particle Shape Factor:| " << WarpX::nox << "\n"; + amrex::Print() << "-------------------------------------------------------------------------------\n"; + // Print solver's type: Yee, CKC, ECT + if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::Yee){ + amrex::Print() << "Maxwell Solver: | Yee \n"; + } + else if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::CKC){ + amrex::Print() << "Maxwell Solver: | CKC \n"; + } + else if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT){ + amrex::Print() << "Maxwell Solver: | ECT \n"; + } + #ifdef WARPX_USE_PSATD + // Print PSATD solver's configuration + if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD){ + amrex::Print() << "Maxwell Solver: | PSATD \n"; + } + if ((m_v_galilean[0]!=0) or (m_v_galilean[1]!=0) or (m_v_galilean[2]!=0)) { + amrex::Print() << " | - Galilean \n" << + " | - v_galilean = (" << m_v_galilean[0] << "," << + m_v_galilean[1] << "," << m_v_galilean[2] << ")\n"; + } + if ((m_v_comoving[0]!=0) or (m_v_comoving[1]!=0) or (m_v_comoving[2]!=0)) { + amrex::Print() << " | - comoving \n" << + " | - v_comoving = (" << m_v_comoving[0] << "," << + m_v_comoving[1] << "," << m_v_comoving[2] << ")\n"; + } + if (WarpX::update_with_rho==1) { + amrex::Print() << " | - update with rho is ON \n"; + } + if (current_correction==1) { + amrex::Print() << " | - current correction is ON \n"; + } + if (WarpX::do_dive_cleaning==1) { + amrex::Print() << " | - div(E) cleaning is ON \n"; + } + if (WarpX::do_divb_cleaning==1) { + amrex::Print() << " | - div(B) cleaning is ON \n"; + } + if (do_multi_J == 1){ + amrex::Print() << " | - multi-J deposition is ON \n"; + amrex::Print() << " | - do_multi_J_n_depositions = " + << WarpX::do_multi_J_n_depositions << "\n"; + } + if (fft_do_time_averaging == 1){ + amrex::Print()<<" | - time-averaged is ON \n"; + } + #endif // WARPX_USE_PSATD + + if (grid_type == GridType::Collocated){ + amrex::Print() << " | - collocated grid \n"; + } + #ifdef WARPX_USE_PSATD + if ( (grid_type == GridType::Staggered) && (field_gathering_algo == GatheringAlgo::EnergyConserving) ){ + amrex::Print()<<" | - staggered grid " << "\n"; + } + else if ( (grid_type == GridType::Hybrid) && (field_gathering_algo == GatheringAlgo::MomentumConserving) ){ + amrex::Print()<<" | - hybrid grid " << "\n"; + if (dims=="3"){ + amrex::Print() << " | - field_centering_nox = " << WarpX::field_centering_nox << "\n"; + amrex::Print() << " | - field_centering_noy = " << WarpX::field_centering_noy << "\n"; + amrex::Print() << " | - field_centering_noz = " << WarpX::field_centering_noz << "\n"; + amrex::Print() << " | - current_centering_nox = " << WarpX::current_centering_nox << "\n"; + amrex::Print() << " | - current_centering_noy = " << WarpX::current_centering_noy << "\n"; + amrex::Print() << " | - current_centering_noz = " << WarpX::current_centering_noz << "\n"; + } + else if (dims=="2"){ + amrex::Print() << " | - field_centering_nox = " << WarpX::field_centering_nox << "\n"; + amrex::Print() << " | - field_centering_noz = " << WarpX::field_centering_noz << "\n"; + amrex::Print() << " | - current_centering_nox = " << WarpX::current_centering_nox << "\n"; + amrex::Print() << " | - current_centering_noz = " << WarpX::current_centering_noz << "\n"; + } + else if (dims=="1"){ + amrex::Print() << " | - field_centering_noz = " << WarpX::field_centering_noz << "\n"; + amrex::Print() << " | - current_centering_noz = " << WarpX::current_centering_noz << "\n"; + } + } + if (WarpX::use_hybrid_QED == true){ + amrex::Print() << " | - use_hybrid_QED = true \n"; + } + + if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD){ + // Print solver's order + std::string psatd_nox_fft, psatd_noy_fft, psatd_noz_fft; + psatd_nox_fft = (nox_fft == -1) ? "inf" : std::to_string(nox_fft); + psatd_noy_fft = (noy_fft == -1) ? "inf" : std::to_string(noy_fft); + psatd_noz_fft = (noz_fft == -1) ? "inf" : std::to_string(noz_fft); + + if (dims=="3" ){ + amrex::Print() << "Spectral order: | - psatd.nox = " << psatd_nox_fft << "\n"; + amrex::Print() << " | - psatd.noy = " << psatd_noy_fft << "\n"; + amrex::Print() << " | - psatd.noz = " << psatd_noz_fft << "\n"; + } + else if (dims=="2" and WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD){ + amrex::Print() << "Spectral order: | - psatd.nox = " << psatd_nox_fft << "\n"; + amrex::Print() << " | - psatd.noz = " << psatd_noz_fft << "\n"; + } + else if (dims=="1" and WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD){ + amrex::Print() << "Spectral order: | - psatd.noz = " << psatd_noz_fft << "\n"; + } + } + // Print guard cells number + amrex::Print() << "Guard cells | - ng_alloc_EB = " << guard_cells.ng_alloc_EB << "\n"; + amrex::Print() << " (allocated for E/B) | \n"; + + #endif // WARPX_USE_PSATD + amrex::Print() << "-------------------------------------------------------------------------------" << "\n"; + //Print main boosted frame algorithm's parameters + if (WarpX::gamma_boost!=1){ + amrex::Print() << "Boosted Frame: | ON \n"; + amrex::Print() << " | - gamma_boost = " << WarpX::gamma_boost << "\n"; + amrex::Print() << " | - boost_direction = (" << WarpX::boost_direction[0] << + "," << WarpX::boost_direction[1] << "," << WarpX::boost_direction[2] << ")\n"; + amrex::Print() << "------------------------------------------------------------------------------- \n"; + } + //Print moving window details + if (WarpX::do_moving_window == 1){ + amrex::Print() << "Moving window: | ON \n"; + if (WarpX::moving_window_dir == 0){ + amrex::Print() << " | - moving_window_dir = x \n"; + } + #if defined(WARPX_DIM_3D) + else if (WarpX::moving_window_dir == 1){ + amrex::Print() << " | - moving_window_dir = y \n"; + } + #endif + else if (WarpX::moving_window_dir == WARPX_ZINDEX) { + amrex::Print() << " | - moving_window_dir = z \n"; + } + amrex::Print() << " | - moving_window_v = " << WarpX::moving_window_v << "\n"; + amrex::Print() << "------------------------------------------------------------------------------- \n"; + } +} + +void +WarpX::WriteUsedInputsFile () const +{ + std::string filename = "warpx_used_inputs"; + ParmParse pp_warpx("warpx"); + pp_warpx.queryAdd("used_inputs_file", filename); + + ablastr::utils::write_used_inputs_file(filename); +} + +void +WarpX::InitData () +{ + WARPX_PROFILE("WarpX::InitData()"); + utils::warpx_check_mpi_thread_level(); + +#ifdef WARPX_QED + Print() << "PICSAR (" << WarpX::PicsarVersion() << ")\n"; +#endif + + Print() << "WarpX (" << WarpX::Version() << ")\n"; + + Print() << utils::logo::get_logo(); + + // WarpX::computeMaxStepBoostAccelerator + // needs to start from the initial zmin_domain_boost, + // even if restarting from a checkpoint file + if (do_compute_max_step_from_zmax) { + zmin_domain_boost_step_0 = geom[0].ProbLo(WARPX_ZINDEX); + } + if (restart_chkfile.empty()) + { + ComputeDt(); + WarpX::PrintDtDxDyDz(); + InitFromScratch(); + if (WarpX::em_solver_medium==1) { + m_macroscopic_properties->InitData(); + } + // Lumped-element InitData must run BEFORE diagnostics InitData so that + // their multifabs exist when CellCenterFunctor captures the pointer. + if (use_lumped_inductor) { + amrex::Print() << " calling inductor (early) \n"; + m_inductor->InitData(); + } + if (use_josephson_junction) { + amrex::Print() << " calling Josephson junction (early) \n"; + m_jj->InitData(); + } + InitDiagnostics(); + } + else + { + InitFromCheckpoint(); + WarpX::PrintDtDxDyDz(); + if (WarpX::em_solver_medium==1) { + m_macroscopic_properties->InitData(); + } + PostRestart(); + // Same ordering constraint as the fresh-start path. + if (use_lumped_inductor) { + amrex::Print() << " calling inductor (early, restart) \n"; + m_inductor->InitData(); + } + if (use_josephson_junction) { + amrex::Print() << " calling Josephson junction (early, restart) \n"; + m_jj->InitData(); + } + reduced_diags->InitData(); + multi_diags->InitData(); + } + + ComputeMaxStep(); + + ComputePMLFactors(); + + if (WarpX::use_fdtd_nci_corr) { + WarpX::InitNCICorrector(); + } + + BuildBufferMasks(); + + + if (WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon) { + amrex::Print() << " calling london \n"; + m_london->InitData(); + } + + // NOTE: Inductor and JosephsonJunction InitData were moved earlier + // (into the if/else branches above) so that diagnostics InitData sees + // their multifabs already allocated. Do not re-call here. + + if (ParallelDescriptor::IOProcessor()) { + std::cout << "\nGrids Summary:\n"; + printGridSummary(std::cout, 0, finestLevel()); + } + + // Check that the number of guard cells is smaller than the number of valid cells for all MultiFabs + // (example: a box with 16 valid cells and 32 guard cells in z will not be considered valid) + CheckGuardCells(); + + PrintMainPICparameters(); + WriteUsedInputsFile(); + + if (restart_chkfile.empty()) + { + // Loop through species and calculate their space-charge field + bool const reset_fields = false; // Do not erase previous user-specified values on the grid + ExecutePythonCallback("beforeInitEsolve"); + ComputeSpaceChargeField(reset_fields); + ExecutePythonCallback("afterInitEsolve"); + if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) + ComputeMagnetostaticField(); + + // Set up an invariant condition through the rest of + // execution, that any code besides the field solver that + // looks at field values will see the composite of the field + // solution and any external field + AddExternalFields(); + + // Write full diagnostics before the first iteration. + multi_diags->FilterComputePackFlush( -1 ); + + // Write reduced diagnostics before the first iteration. + if (reduced_diags->m_plot_rd != 0) + { + reduced_diags->ComputeDiags(-1); + reduced_diags->WriteToFile(-1); + } + } + + PerformanceHints(); + + CheckKnownIssues(); +} + +void +WarpX::AddExternalFields () { + for (int lev = 0; lev <= finest_level; ++lev) { + // FIXME: RZ multimode has more than one component for all these + if (add_external_E_field) { + amrex::MultiFab::Add(*Efield_fp[lev][0], *Efield_fp_external[lev][0], 0, 0, 1, guard_cells.ng_alloc_EB); + amrex::MultiFab::Add(*Efield_fp[lev][1], *Efield_fp_external[lev][1], 0, 0, 1, guard_cells.ng_alloc_EB); + amrex::MultiFab::Add(*Efield_fp[lev][2], *Efield_fp_external[lev][2], 0, 0, 1, guard_cells.ng_alloc_EB); + } + if (add_external_B_field) { + amrex::MultiFab::Add(*Bfield_fp[lev][0], *Bfield_fp_external[lev][0], 0, 0, 1, guard_cells.ng_alloc_EB); + amrex::MultiFab::Add(*Bfield_fp[lev][1], *Bfield_fp_external[lev][1], 0, 0, 1, guard_cells.ng_alloc_EB); + amrex::MultiFab::Add(*Bfield_fp[lev][2], *Bfield_fp_external[lev][2], 0, 0, 1, guard_cells.ng_alloc_EB); + } + } +} + +void +WarpX::InitDiagnostics () { + multi_diags->InitData(); + reduced_diags->InitData(); +} + +void +WarpX::InitFromScratch () +{ + const Real time = 0.0; + + AmrCore::InitFromScratch(time); // This will call MakeNewLevelFromScratch + + mypc->AllocData(); + mypc->InitData(); + + InitPML(); +} + +void +WarpX::InitPML () +{ + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (WarpX::field_boundary_lo[idim] == FieldBoundaryType::PML) { + do_pml = 1; + do_pml_Lo[0][idim] = 1; // on level 0 + } + if (WarpX::field_boundary_hi[idim] == FieldBoundaryType::PML) { + do_pml = 1; + do_pml_Hi[0][idim] = 1; // on level 0 + } + } + if (finest_level > 0) do_pml = 1; + if (do_pml) + { +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + do_pml_Lo[0][0] = 0; // no PML at r=0, in cylindrical geometry + pml_rz[0] = std::make_unique(0, boxArray(0), DistributionMap(0), &Geom(0), pml_ncell, do_pml_in_domain); +#else + // Note: fill_guards_fields and fill_guards_current are both set to + // zero (amrex::IntVect(0)) (what we do with damping BCs does not apply + // to the PML, for example in the presence of mesh refinement patches) + pml[0] = std::make_unique(0, boxArray(0), DistributionMap(0), &Geom(0), nullptr, + pml_ncell, pml_delta, amrex::IntVect::TheZeroVector(), + dt[0], nox_fft, noy_fft, noz_fft, grid_type, + do_moving_window, pml_has_particles, do_pml_in_domain, + psatd_solution_type, J_in_time, rho_in_time, + do_pml_dive_cleaning, do_pml_divb_cleaning, + amrex::IntVect(0), amrex::IntVect(0), + guard_cells.ng_FieldSolver.max(), + v_particle_pml, + do_pml_Lo[0], do_pml_Hi[0]); +#endif + + for (int lev = 1; lev <= finest_level; ++lev) + { + do_pml_Lo[lev] = amrex::IntVect::TheUnitVector(); + do_pml_Hi[lev] = amrex::IntVect::TheUnitVector(); + // check if fine patch edges co-incide with domain boundary + amrex::Box levelBox = boxArray(lev).minimalBox(); + // Domain box at level, lev + amrex::Box DomainBox = Geom(lev).Domain(); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (levelBox.smallEnd(idim) == DomainBox.smallEnd(idim)) + do_pml_Lo[lev][idim] = do_pml_Lo[0][idim]; + if (levelBox.bigEnd(idim) == DomainBox.bigEnd(idim)) + do_pml_Hi[lev][idim] = do_pml_Hi[0][idim]; + } + +#ifdef WARPX_DIM_RZ + //In cylindrical geometry, if the edge of the patch is at r=0, do not add PML + if ((max_level > 0) && (fine_tag_lo[0]==0.)) { + do_pml_Lo[lev][0] = 0; + } +#endif + // Note: fill_guards_fields and fill_guards_current are both set to + // zero (amrex::IntVect(0)) (what we do with damping BCs does not apply + // to the PML, for example in the presence of mesh refinement patches) + pml[lev] = std::make_unique(lev, boxArray(lev), DistributionMap(lev), + &Geom(lev), &Geom(lev-1), + pml_ncell, pml_delta, refRatio(lev-1), + dt[lev], nox_fft, noy_fft, noz_fft, grid_type, + do_moving_window, pml_has_particles, do_pml_in_domain, + psatd_solution_type, J_in_time, rho_in_time, do_pml_dive_cleaning, do_pml_divb_cleaning, + amrex::IntVect(0), amrex::IntVect(0), + guard_cells.ng_FieldSolver.max(), + v_particle_pml, + do_pml_Lo[lev], do_pml_Hi[lev]); + } + } +} + +void +WarpX::ComputePMLFactors () +{ + if (do_pml) + { + for (int lev = 0; lev <= finest_level; ++lev) + { + if (pml[lev]) + pml[lev]->ComputePMLFactors(dt[lev]); + } + } +} + +void +WarpX::ComputeMaxStep () +{ + if (do_compute_max_step_from_zmax) { + computeMaxStepBoostAccelerator(); + } +} + +/* \brief computes max_step for wakefield simulation in boosted frame. + * \param geom: Geometry object that contains simulation domain. + * + * max_step is set so that the simulation stop when the lower corner of the + * simulation box passes input parameter zmax_plasma_to_compute_max_step. + */ +void +WarpX::computeMaxStepBoostAccelerator() { + // Sanity checks: can use zmax_plasma_to_compute_max_step only if + // the moving window and the boost are all in z direction. + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + WarpX::moving_window_dir == WARPX_ZINDEX, + "Can use zmax_plasma_to_compute_max_step only if " + "moving window along z. TODO: all directions."); + if (gamma_boost > 1){ + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + (WarpX::boost_direction[0]-0)*(WarpX::boost_direction[0]-0) + + (WarpX::boost_direction[1]-0)*(WarpX::boost_direction[1]-0) + + (WarpX::boost_direction[2]-1)*(WarpX::boost_direction[2]-1) < 1.e-12, + "Can use zmax_plasma_to_compute_max_step in boosted frame only if " + "warpx.boost_direction = z. TODO: all directions."); + } + + // Lower end of the simulation domain. All quantities are given in boosted + // frame except zmax_plasma_to_compute_max_step. + + // End of the plasma: Transform input argument + // zmax_plasma_to_compute_max_step to boosted frame. + const Real len_plasma_boost = zmax_plasma_to_compute_max_step/gamma_boost; + // Plasma velocity + const Real v_plasma_boost = -beta_boost * PhysConst::c; + // Get time at which the lower end of the simulation domain passes the + // upper end of the plasma (in the z direction). + const Real interaction_time_boost = (len_plasma_boost-zmin_domain_boost_step_0)/ + (moving_window_v-v_plasma_boost); + // Divide by dt, and update value of max_step. + int computed_max_step; + if (do_subcycling){ + computed_max_step = static_cast(interaction_time_boost/dt[0]); + } else { + computed_max_step = + static_cast(interaction_time_boost/dt[maxLevel()]); + } + max_step = computed_max_step; + Print()<<"max_step computed in computeMaxStepBoostAccelerator: " + <( + godfrey_coeff_set::Ex_Ey_Bz, cdtodz, nodal_gather); + // Same filter for fields Bx, By and Ez + nci_godfrey_filter_bxbyez[lev] = std::make_unique( + godfrey_coeff_set::Bx_By_Ez, cdtodz, nodal_gather); + // Compute Godfrey filters stencils + nci_godfrey_filter_exeybz[lev]->ComputeStencils(); + nci_godfrey_filter_bxbyez[lev]->ComputeStencils(); + } + } +#endif +} + +void +WarpX::InitFilter (){ + if (WarpX::use_filter){ + WarpX::bilinear_filter.npass_each_dir = WarpX::filter_npass_each_dir.toArray(); + WarpX::bilinear_filter.ComputeStencils(); + } +} + +void +WarpX::PostRestart () +{ + mypc->PostRestart(); +} + + +void +WarpX::InitLevelData (int lev, Real /*time*/) +{ + + ParmParse pp_warpx("warpx"); + + // default values of E_external_grid and B_external_grid + // are used to set the E and B field when "constant" or + // "parser" is not explicitly used in the input. + pp_warpx.query("B_ext_grid_init_style", B_ext_grid_s); + std::transform(B_ext_grid_s.begin(), + B_ext_grid_s.end(), + B_ext_grid_s.begin(), + ::tolower); + +#ifdef WARPX_MAG_LLG + if (pp_warpx.query("B_ext_grid_init_style", B_ext_grid_s) ) { + amrex::Abort("ERROR: Initialization of B field is not allowed in the LLG simulation! \nThe initial magnetic field must be H and M! \n"); + } +#endif + + pp_warpx.query("E_ext_grid_init_style", E_ext_grid_s); + std::transform(E_ext_grid_s.begin(), + E_ext_grid_s.end(), + E_ext_grid_s.begin(), + ::tolower); +#ifdef WARPX_MAG_LLG + pp_warpx.query("M_ext_grid_init_style", M_ext_grid_s); // user-defined initial M + std::transform(M_ext_grid_s.begin(), + M_ext_grid_s.end(), + M_ext_grid_s.begin(), + ::tolower); + + pp_warpx.query("H_ext_grid_init_style", H_ext_grid_s); // user-defined initial H + std::transform(H_ext_grid_s.begin(), + H_ext_grid_s.end(), + H_ext_grid_s.begin(), + ::tolower); + + pp_warpx.query("H_bias_ext_grid_init_style", H_bias_ext_grid_s); // user-defined initial M + std::transform(H_bias_ext_grid_s.begin(), + H_bias_ext_grid_s.end(), + H_bias_ext_grid_s.begin(), + ::tolower); +#endif + + // * Functions with the string "arr" in their names get an Array of + // values from the given entry in the table. The array argument is + // resized (if necessary) to hold all the values requested. + // + // * Functions without the string "arr" in their names get single + // values from the given entry in the table. + + // if the input string is "constant", the values for the + // external grid must be provided in the input. + if (B_ext_grid_s == "constant") + utils::parser::getArrWithParser(pp_warpx, "B_external_grid", B_external_grid); + + // if the input string is "constant", the values for the + // external grid must be provided in the input. + if (E_ext_grid_s == "constant") + utils::parser::getArrWithParser(pp_warpx, "E_external_grid", E_external_grid); + + +#ifdef WARPX_MAG_LLG + if (M_ext_grid_s == "constant") + utils::parser::getArrWithParser(pp_warpx, "M_external_grid", M_external_grid); + + if (H_ext_grid_s == "constant") + utils::parser::getArrWithParser(pp_warpx, "H_external_grid", H_external_grid); + + if (H_bias_ext_grid_s == "constant") + utils::parser::getArrWithParser(pp_warpx,"H_bias_external_grid", H_bias_external_grid); +#endif + // initialize the averaged fields only if the averaged algorithm + // is activated ('psatd.do_time_averaging=1') + ParmParse pp_psatd("psatd"); + pp_psatd.query("do_time_averaging", fft_do_time_averaging ); + + for (int i = 0; i < 3; ++i) { + current_fp[lev][i]->setVal(0.0); + if (lev > 0) + current_cp[lev][i]->setVal(0.0); + Bfield_sc_fp[lev][i]->setVal(0.0); + + if (B_ext_grid_s == "constant" || B_ext_grid_s == "default") { + Bfield_fp[lev][i]->setVal(B_external_grid[i]); + if (fft_do_time_averaging) { + Bfield_avg_fp[lev][i]->setVal(B_external_grid[i]); + } + + if (lev > 0) { + Bfield_aux[lev][i]->setVal(B_external_grid[i]); + Bfield_cp[lev][i]->setVal(B_external_grid[i]); + if (fft_do_time_averaging) { + Bfield_avg_cp[lev][i]->setVal(B_external_grid[i]); + } + } + } + if (E_ext_grid_s == "constant" || E_ext_grid_s == "default") { + Efield_fp[lev][i]->setVal(E_external_grid[i]); + if (fft_do_time_averaging) { + Efield_avg_fp[lev][i]->setVal(E_external_grid[i]); + } + + if (lev > 0) { + Efield_aux[lev][i]->setVal(E_external_grid[i]); + Efield_cp[lev][i]->setVal(E_external_grid[i]); + if (fft_do_time_averaging) { + Efield_avg_cp[lev][i]->setVal(E_external_grid[i]); + } + } + } + +#ifdef WARPX_MAG_LLG + if (M_ext_grid_s == "constant" || M_ext_grid_s == "default"){ + // this if condition finds out if the user-input is constant + // if not, set initial value to default, default = 0.0 + + // Set the value of num_comp components in the valid region of + // each FAB in the FabArray, starting at component comp to val. + // Also set the value of nghost boundary cells. + // template ::value>::type > + // void setVal (value_type val, + // int comp, + // int num_comp, + // int nghost = 0); + + int nghost = 1; + for (int icomp = 0; icomp < 3; ++icomp){ // icomp is the index of components at each i face + Mfield_fp[lev][i]->setVal(M_external_grid[icomp], icomp, 1, nghost); + } + } + + if (H_ext_grid_s == "constant" || H_ext_grid_s == "default") { + Hfield_fp[lev][i]->setVal(H_external_grid[i]); + if (lev > 0) { + Hfield_aux[lev][i]->setVal(H_external_grid[i]); + Hfield_cp[lev][i]->setVal(H_external_grid[i]); + } + } + + if (H_bias_ext_grid_s == "constant" || H_bias_ext_grid_s == "default") { + H_biasfield_fp[lev][i]->setVal(H_bias_external_grid[i]); + if (lev > 0) { + H_biasfield_aux[lev][i]->setVal(H_bias_external_grid[i]); + H_biasfield_cp[lev][i]->setVal(H_bias_external_grid[i]); + } + } + +#endif + } + +#ifdef AMREX_USE_EB + InitializeEBGridData(lev); +#endif + + // if the input string for the B-field is "parse_b_ext_grid_function", + // then the analytical expression or function must be + // provided in the input file. + if (B_ext_grid_s == "parse_b_ext_grid_function") { + +#ifdef WARPX_DIM_RZ + WARPX_ABORT_WITH_MESSAGE( + "E and B parser for external fields does not work with RZ -- TO DO"); +#endif + utils::parser::Store_parserString(pp_warpx, "Bx_external_grid_function(x,y,z)", + str_Bx_ext_grid_function); + utils::parser::Store_parserString(pp_warpx, "By_external_grid_function(x,y,z)", + str_By_ext_grid_function); + utils::parser::Store_parserString(pp_warpx, "Bz_external_grid_function(x,y,z)", + str_Bz_ext_grid_function); + Bxfield_parser = std::make_unique( + utils::parser::makeParser(str_Bx_ext_grid_function,{"x","y","z"})); + Byfield_parser = std::make_unique( + utils::parser::makeParser(str_By_ext_grid_function,{"x","y","z"})); + Bzfield_parser = std::make_unique( + utils::parser::makeParser(str_Bz_ext_grid_function,{"x","y","z"})); + + // Initialize Bfield_fp with external function + InitializeExternalFieldsOnGridUsingParser(Bfield_fp[lev][0].get(), + Bfield_fp[lev][1].get(), + Bfield_fp[lev][2].get(), + Bxfield_parser->compile<3>(), + Byfield_parser->compile<3>(), + Bzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'B', + lev, PatchType::fine); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(Bfield_aux[lev][0].get(), + Bfield_aux[lev][1].get(), + Bfield_aux[lev][2].get(), + Bxfield_parser->compile<3>(), + Byfield_parser->compile<3>(), + Bzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'B', + lev, PatchType::fine); + + InitializeExternalFieldsOnGridUsingParser(Bfield_cp[lev][0].get(), + Bfield_cp[lev][1].get(), + Bfield_cp[lev][2].get(), + Bxfield_parser->compile<3>(), + Byfield_parser->compile<3>(), + Bzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'B', + lev, PatchType::coarse); + } + } + + // if the input string for the E-field is "parse_e_ext_grid_function", + // then the analytical expression or function must be + // provided in the input file. + if (E_ext_grid_s == "parse_e_ext_grid_function") { + +#ifdef WARPX_DIM_RZ + WARPX_ABORT_WITH_MESSAGE( + "E and B parser for external fields does not work with RZ -- TO DO"); +#endif + utils::parser::Store_parserString(pp_warpx, "Ex_external_grid_function(x,y,z)", + str_Ex_ext_grid_function); + utils::parser::Store_parserString(pp_warpx, "Ey_external_grid_function(x,y,z)", + str_Ey_ext_grid_function); + utils::parser::Store_parserString(pp_warpx, "Ez_external_grid_function(x,y,z)", + str_Ez_ext_grid_function); + + Exfield_parser = std::make_unique( + utils::parser::makeParser(str_Ex_ext_grid_function,{"x","y","z"})); + Eyfield_parser = std::make_unique( + utils::parser::makeParser(str_Ey_ext_grid_function,{"x","y","z"})); + Ezfield_parser = std::make_unique( + utils::parser::makeParser(str_Ez_ext_grid_function,{"x","y","z"})); + + // Initialize Efield_fp with external function + InitializeExternalFieldsOnGridUsingParser(Efield_fp[lev][0].get(), + Efield_fp[lev][1].get(), + Efield_fp[lev][2].get(), + Exfield_parser->compile<3>(), + Eyfield_parser->compile<3>(), + Ezfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'E', + lev, PatchType::fine); + +#ifdef AMREX_USE_EB + // We initialize ECTRhofield consistently with the Efield + if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { + m_fdtd_solver_fp[lev]->EvolveECTRho(Efield_fp[lev], m_edge_lengths[lev], + m_face_areas[lev], ECTRhofield[lev], lev); + + } +#endif + + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(Efield_aux[lev][0].get(), + Efield_aux[lev][1].get(), + Efield_aux[lev][2].get(), + Exfield_parser->compile<3>(), + Eyfield_parser->compile<3>(), + Ezfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'E', + lev, PatchType::fine); + + InitializeExternalFieldsOnGridUsingParser(Efield_cp[lev][0].get(), + Efield_cp[lev][1].get(), + Efield_cp[lev][2].get(), + Exfield_parser->compile<3>(), + Eyfield_parser->compile<3>(), + Ezfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'E', + lev, PatchType::coarse); +#ifdef AMREX_USE_EB + if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { + // We initialize ECTRhofield consistently with the Efield + m_fdtd_solver_cp[lev]->EvolveECTRho(Efield_cp[lev], m_edge_lengths[lev], + m_face_areas[lev], ECTRhofield[lev], lev); + + } +#endif + } + } + +#ifdef WARPX_MAG_LLG + // if the input string for the Hbias-field is "parse_h_bias_ext_grid_function", + // then the analytical expression or function must be + // provided in the input file. + if (H_bias_ext_grid_s == "parse_h_bias_ext_grid_function") { + +#ifdef WARPX_DIM_RZ + amrex::Abort("H bias parser for external fields does not work with RZ -- TO DO"); +#endif + utils::parser::Store_parserString(pp_warpx, "Hx_bias_external_grid_function(x,y,z)", + str_Hx_bias_ext_grid_function); + utils::parser::Store_parserString(pp_warpx, "Hy_bias_external_grid_function(x,y,z)", + str_Hy_bias_ext_grid_function); + utils::parser::Store_parserString(pp_warpx, "Hz_bias_external_grid_function(x,y,z)", + str_Hz_bias_ext_grid_function); + + Hx_biasfield_parser = std::make_unique( + utils::parser::makeParser(str_Hx_bias_ext_grid_function,{"x","y","z"})); + Hy_biasfield_parser = std::make_unique( + utils::parser::makeParser(str_Hy_bias_ext_grid_function,{"x","y","z"})); + Hz_biasfield_parser = std::make_unique( + utils::parser::makeParser(str_Hz_bias_ext_grid_function,{"x","y","z"})); + + // Initialize Efield_fp with external function + InitializeExternalFieldsOnGridUsingParser(H_biasfield_fp[lev][0].get(), + H_biasfield_fp[lev][1].get(), + H_biasfield_fp[lev][2].get(), + Hx_biasfield_parser->compile<3>(), + Hy_biasfield_parser->compile<3>(), + Hz_biasfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'H', + lev, PatchType::fine); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(H_biasfield_aux[lev][0].get(), + H_biasfield_aux[lev][1].get(), + H_biasfield_aux[lev][2].get(), + Hx_biasfield_parser->compile<3>(), + Hy_biasfield_parser->compile<3>(), + Hz_biasfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'H', + lev, PatchType::coarse); + + InitializeExternalFieldsOnGridUsingParser(H_biasfield_cp[lev][0].get(), + H_biasfield_cp[lev][1].get(), + H_biasfield_cp[lev][2].get(), + Hx_biasfield_parser->compile<3>(), + Hy_biasfield_parser->compile<3>(), + Hz_biasfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'H', + lev, PatchType::coarse); + } + } + + if (H_ext_grid_s == "parse_h_ext_grid_function") { + +#ifdef WARPX_DIM_RZ + amrex::Abort("H parser for external fields does not work with RZ -- TO DO"); +#endif + utils::parser::Store_parserString(pp_warpx, "Hx_external_grid_function(x,y,z)", + str_Hx_ext_grid_function); + utils::parser::Store_parserString(pp_warpx, "Hy_external_grid_function(x,y,z)", + str_Hy_ext_grid_function); + utils::parser::Store_parserString(pp_warpx, "Hz_external_grid_function(x,y,z)", + str_Hz_ext_grid_function); + + Hxfield_parser = std::make_unique( + utils::parser::makeParser(str_Hx_ext_grid_function,{"x","y","z"})); + Hyfield_parser = std::make_unique( + utils::parser::makeParser(str_Hy_ext_grid_function,{"x","y","z"})); + Hzfield_parser = std::make_unique( + utils::parser::makeParser(str_Hz_ext_grid_function,{"x","y","z"})); + + // Initialize Hfield_fp with external function + InitializeExternalFieldsOnGridUsingParser(Hfield_fp[lev][0].get(), + Hfield_fp[lev][1].get(), + Hfield_fp[lev][2].get(), + Hxfield_parser->compile<3>(), + Hyfield_parser->compile<3>(), + Hzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'H', + lev, PatchType::fine); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(Hfield_aux[lev][0].get(), + Hfield_aux[lev][1].get(), + Hfield_aux[lev][2].get(), + Hxfield_parser->compile<3>(), + Hyfield_parser->compile<3>(), + Hzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'H', + lev, PatchType::coarse); + + InitializeExternalFieldsOnGridUsingParser(Hfield_cp[lev][0].get(), + Hfield_cp[lev][1].get(), + Hfield_cp[lev][2].get(), + Hxfield_parser->compile<3>(), + Hyfield_parser->compile<3>(), + Hzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'H', + lev, PatchType::coarse); + } + } + + if (M_ext_grid_s == "parse_m_ext_grid_function") { +#ifdef WARPX_DIM_RZ + amrex::Abort("M-field parser for external fields does not work with RZ"); +#endif + utils::parser::Store_parserString(pp_warpx, "Mx_external_grid_function(x,y,z)", + str_Mx_ext_grid_function); + utils::parser::Store_parserString(pp_warpx, "My_external_grid_function(x,y,z)", + str_My_ext_grid_function); + utils::parser::Store_parserString(pp_warpx, "Mz_external_grid_function(x,y,z)", + str_Mz_ext_grid_function); + + Mxfield_parser = std::make_unique( + utils::parser::makeParser(str_Mx_ext_grid_function,{"x","y","z"})); + Myfield_parser = std::make_unique( + utils::parser::makeParser(str_My_ext_grid_function,{"x","y","z"})); + Mzfield_parser = std::make_unique( + utils::parser::makeParser(str_Mz_ext_grid_function,{"x","y","z"})); + + // Initialize Mfield_fp with external function directly on the faces + InitializeExternalFieldsOnGridUsingParser(Mfield_fp[lev][0].get(), + Mfield_fp[lev][1].get(), + Mfield_fp[lev][2].get(), + Mxfield_parser->compile<3>(), + Myfield_parser->compile<3>(), + Mzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'M', + lev, PatchType::fine); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(Mfield_aux[lev][0].get(), + Mfield_aux[lev][1].get(), + Mfield_aux[lev][2].get(), + Mxfield_parser->compile<3>(), + Myfield_parser->compile<3>(), + Mzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'M', + lev, PatchType::coarse); + + InitializeExternalFieldsOnGridUsingParser(Mfield_cp[lev][0].get(), + Mfield_cp[lev][1].get(), + Mfield_cp[lev][2].get(), + Mxfield_parser->compile<3>(), + Myfield_parser->compile<3>(), + Mzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'M', + lev, PatchType::coarse); + } + } + +#endif //closes #ifdef WARPX_MAG_LLG + + // Reading external fields from data file + if (add_external_B_field) { + std::string read_fields_from_path="./"; + pp_warpx.query("read_fields_from_path", read_fields_from_path); +#if defined(WARPX_DIM_RZ) + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(n_rz_azimuthal_modes == 1, + "External field reading is not implemented for more than one RZ mode (see #3829)"); + ReadExternalFieldFromFile(read_fields_from_path, Bfield_fp_external[lev][0].get(), "B", "r"); + ReadExternalFieldFromFile(read_fields_from_path, Bfield_fp_external[lev][1].get(), "B", "t"); + ReadExternalFieldFromFile(read_fields_from_path, Bfield_fp_external[lev][2].get(), "B", "z"); +#else + ReadExternalFieldFromFile(read_fields_from_path, Bfield_fp_external[lev][0].get(), "B", "x"); + ReadExternalFieldFromFile(read_fields_from_path, Bfield_fp_external[lev][1].get(), "B", "y"); + ReadExternalFieldFromFile(read_fields_from_path, Bfield_fp_external[lev][2].get(), "B", "z"); +#endif + } + if (add_external_E_field) { + std::string read_fields_from_path="./"; + pp_warpx.query("read_fields_from_path", read_fields_from_path); +#if defined(WARPX_DIM_RZ) + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(n_rz_azimuthal_modes == 1, + "External field reading is not implemented for more than one RZ mode (see #3829)"); + ReadExternalFieldFromFile(read_fields_from_path, Efield_fp_external[lev][0].get(), "E", "r"); + ReadExternalFieldFromFile(read_fields_from_path, Efield_fp_external[lev][1].get(), "E", "t"); + ReadExternalFieldFromFile(read_fields_from_path, Efield_fp_external[lev][2].get(), "E", "z"); +#else + ReadExternalFieldFromFile(read_fields_from_path, Efield_fp_external[lev][0].get(), "E", "x"); + ReadExternalFieldFromFile(read_fields_from_path, Efield_fp_external[lev][1].get(), "E", "y"); + ReadExternalFieldFromFile(read_fields_from_path, Efield_fp_external[lev][2].get(), "E", "z"); +#endif + } + + if (costs[lev]) { + const auto iarr = costs[lev]->IndexArray(); + for (int i : iarr) { + (*costs[lev])[i] = 0.0; + WarpX::setLoadBalanceEfficiency(lev, -1); + } + } +} + +#ifdef WARPX_MAG_LLG +void WarpX::AverageParsedMtoFaces(MultiFab& Mx_cc, + MultiFab& My_cc, + MultiFab& Mz_cc, + MultiFab& Mx_face, + MultiFab& My_face, + MultiFab& Mz_face) +{ + // average Mx, My, Mz to faces + for (MFIter mfi(Mx_face, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + amrex::IntVect x_nodal_flag = Mx_face.ixType().toIntVect(); + amrex::IntVect y_nodal_flag = My_face.ixType().toIntVect(); + amrex::IntVect z_nodal_flag = Mz_face.ixType().toIntVect(); + const amrex::Box& tbx = mfi.tilebox( x_nodal_flag, Mx_face.nGrowVect() ); + const amrex::Box& tby = mfi.tilebox( y_nodal_flag, My_face.nGrowVect() ); + const amrex::Box& tbz = mfi.tilebox( z_nodal_flag, Mz_face.nGrowVect() ); + + auto const& mx_cc = Mx_cc.array(mfi); + auto const& my_cc = My_cc.array(mfi); + auto const& mz_cc = Mz_cc.array(mfi); + + auto const& mx_face = Mx_face.array(mfi); + auto const& my_face = My_face.array(mfi); + auto const& mz_face = Mz_face.array(mfi); + + amrex::ParallelFor (tbx, tby, tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + mx_face(i,j,k,0) = 0.5*(mx_cc(i-1,j,k) + mx_cc(i,j,k)); + mx_face(i,j,k,1) = 0.5*(my_cc(i-1,j,k) + my_cc(i,j,k)); + mx_face(i,j,k,2) = 0.5*(mz_cc(i-1,j,k) + mz_cc(i,j,k)); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + my_face(i,j,k,0) = 0.5*(mx_cc(i,j-1,k) + mx_cc(i,j,k)); + my_face(i,j,k,1) = 0.5*(my_cc(i,j-1,k) + my_cc(i,j,k)); + my_face(i,j,k,2) = 0.5*(mz_cc(i,j-1,k) + mz_cc(i,j,k)); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + mz_face(i,j,k,0) = 0.5*(mx_cc(i,j,k-1) + mx_cc(i,j,k)); + mz_face(i,j,k,1) = 0.5*(my_cc(i,j,k-1) + my_cc(i,j,k)); + mz_face(i,j,k,2) = 0.5*(mz_cc(i,j,k-1) + mz_cc(i,j,k)); + }); + } +} +#endif + +void +WarpX::InitializeExternalFieldsOnGridUsingParser ( + MultiFab *mfx, MultiFab *mfy, MultiFab *mfz, + ParserExecutor<3> const& xfield_parser, ParserExecutor<3> const& yfield_parser, + ParserExecutor<3> const& zfield_parser, + std::array< std::unique_ptr, 3 > const& edge_lengths, + std::array< std::unique_ptr, 3 > const& face_areas, + const char field, + const int lev, PatchType patch_type) +{ + + auto dx_lev = geom[lev].CellSizeArray(); + amrex::IntVect refratio = (lev > 0 ) ? WarpX::RefRatio(lev-1) : amrex::IntVect(1); + if (patch_type == PatchType::coarse) { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + dx_lev[idim] = dx_lev[idim] * refratio[idim]; + } + } + + const RealBox& real_box = geom[lev].ProbDomain(); + amrex::IntVect x_nodal_flag = mfx->ixType().toIntVect(); + amrex::IntVect y_nodal_flag = mfy->ixType().toIntVect(); + amrex::IntVect z_nodal_flag = mfz->ixType().toIntVect(); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + mfx->nComp() == mfy->nComp() and mfx->nComp() == mfz->nComp(), + "The number of components for the three Multifabs must be equal"); + // Number of multifab components +#ifdef WARPX_MAG_LLG + int ncomp = mfx->nComp(); +#endif + + for ( MFIter mfi(*mfx, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + + const amrex::Box& tbx = mfi.tilebox( x_nodal_flag, mfx->nGrowVect() ); + const amrex::Box& tby = mfi.tilebox( y_nodal_flag, mfy->nGrowVect() ); + const amrex::Box& tbz = mfi.tilebox( z_nodal_flag, mfz->nGrowVect() ); + + auto const& mfxfab = mfx->array(mfi); + auto const& mfyfab = mfy->array(mfi); + auto const& mfzfab = mfz->array(mfi); + +#ifdef AMREX_USE_EB + amrex::Array4 const& lx = edge_lengths[0]->array(mfi); + amrex::Array4 const& ly = edge_lengths[1]->array(mfi); + amrex::Array4 const& lz = edge_lengths[2]->array(mfi); + amrex::Array4 const& Sx = face_areas[0]->array(mfi); + amrex::Array4 const& Sy = face_areas[1]->array(mfi); + amrex::Array4 const& Sz = face_areas[2]->array(mfi); + +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + const amrex::Dim3 lx_lo = amrex::lbound(lx); + const amrex::Dim3 lx_hi = amrex::ubound(lx); + const amrex::Dim3 lz_lo = amrex::lbound(lz); + const amrex::Dim3 lz_hi = amrex::ubound(lz); +#endif + +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + amrex::ignore_unused(ly, Sx, Sz); +#elif defined(WARPX_DIM_1D_Z) + amrex::ignore_unused(lx, ly, lz, Sx, Sy, Sz); +#endif + +#else + amrex::ignore_unused(edge_lengths, face_areas, field); +#endif + + amrex::ParallelFor (tbx, tby, tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { +#ifdef AMREX_USE_EB +#ifdef WARPX_DIM_3D + if((field=='E' and lx(i, j, k)<=0) or (field=='B' and Sx(i, j, k)<=0)) return; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + //In XZ and RZ Ex is associated with a x-edge, while Bx is associated with a z-edge + if((field=='E' and lx(i, j, k)<=0) or (field=='B' and lz(i, j, k)<=0)) return; +#endif +#endif + // Shift required in the x-, y-, or z- position + // depending on the index type of the multifab +#if defined(WARPX_DIM_1D_Z) + amrex::Real x = 0._rt; + amrex::Real y = 0._rt; + amrex::Real fac_z = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real z = j*dx_lev[0] + real_box.lo(0) + fac_z; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + amrex::Real fac_x = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; + amrex::Real y = 0._rt; + amrex::Real fac_z = (1._rt - x_nodal_flag[1]) * dx_lev[1] * 0.5_rt; + amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#else + amrex::Real fac_x = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; + amrex::Real fac_y = (1._rt - x_nodal_flag[1]) * dx_lev[1] * 0.5_rt; + amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; + amrex::Real fac_z = (1._rt - x_nodal_flag[2]) * dx_lev[2] * 0.5_rt; + amrex::Real z = k*dx_lev[2] + real_box.lo(2) + fac_z; +#endif +#ifdef WARPX_MAG_LLG + if (ncomp > 1) { + // This condition is specific to Mfield, where, + // x-, y-, and z-components are stored on the x-face + mfxfab(i,j,k,0) = xfield_parser(x,y,z); + mfxfab(i,j,k,1) = yfield_parser(x,y,z); + mfxfab(i,j,k,2) = zfield_parser(x,y,z); + } else +#endif + { + mfxfab(i,j,k) = xfield_parser(x,y,z); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { +#ifdef AMREX_USE_EB +#ifdef WARPX_DIM_3D + if((field=='E' and ly(i, j, k)<=0) or (field=='B' and Sy(i, j, k)<=0)) return; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + //In XZ and RZ Ey is associated with a mesh node, so we need to check if the mesh node is covered + if((field=='E' and (lx(std::min(i , lx_hi.x), std::min(j , lx_hi.y), k)<=0 + || lx(std::max(i-1, lx_lo.x), std::min(j , lx_hi.y), k)<=0 + || lz(std::min(i , lz_hi.x), std::min(j , lz_hi.y), k)<=0 + || lz(std::min(i , lz_hi.x), std::max(j-1, lz_lo.y), k)<=0)) or + (field=='B' and Sy(i,j,k)<=0)) return; +#endif +#endif +#if defined(WARPX_DIM_1D_Z) + amrex::Real x = 0._rt; + amrex::Real y = 0._rt; + amrex::Real fac_z = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real z = j*dx_lev[0] + real_box.lo(0) + fac_z; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + amrex::Real fac_x = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; + amrex::Real y = 0._rt; + amrex::Real fac_z = (1._rt - y_nodal_flag[1]) * dx_lev[1] * 0.5_rt; + amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#elif defined(WARPX_DIM_3D) + amrex::Real fac_x = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; + amrex::Real fac_y = (1._rt - y_nodal_flag[1]) * dx_lev[1] * 0.5_rt; + amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; + amrex::Real fac_z = (1._rt - y_nodal_flag[2]) * dx_lev[2] * 0.5_rt; + amrex::Real z = k*dx_lev[2] + real_box.lo(2) + fac_z; +#endif +#ifdef WARPX_MAG_LLG + if (ncomp > 1) { + // This condition is specific to Mfield, where, + // x-, y-, and z-components are stored on the y-face + mfyfab(i,j,k,0) = xfield_parser(x,y,z); + mfyfab(i,j,k,1) = yfield_parser(x,y,z); + mfyfab(i,j,k,2) = zfield_parser(x,y,z); + } else +#endif + { + mfyfab(i,j,k) = yfield_parser(x,y,z); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { +#ifdef AMREX_USE_EB +#ifdef WARPX_DIM_3D + if((field=='E' and lz(i, j, k)<=0) or (field=='B' and Sz(i, j, k)<=0)) return; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + //In XZ and RZ Ez is associated with a z-edge, while Bz is associated with a x-edge + if((field=='E' and lz(i, j, k)<=0) or (field=='B' and lx(i, j, k)<=0)) return; +#endif +#endif +#if defined(WARPX_DIM_1D_Z) + amrex::Real x = 0._rt; + amrex::Real y = 0._rt; + amrex::Real fac_z = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real z = j*dx_lev[0] + real_box.lo(0) + fac_z; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + amrex::Real fac_x = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; + amrex::Real y = 0._rt; + amrex::Real fac_z = (1._rt - z_nodal_flag[1]) * dx_lev[1] * 0.5_rt; + amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#elif defined(WARPX_DIM_3D) + amrex::Real fac_x = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; + amrex::Real fac_y = (1._rt - z_nodal_flag[1]) * dx_lev[1] * 0.5_rt; + amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; + amrex::Real fac_z = (1._rt - z_nodal_flag[2]) * dx_lev[2] * 0.5_rt; + amrex::Real z = k*dx_lev[2] + real_box.lo(2) + fac_z; +#endif +#ifdef WARPX_MAG_LLG + if (ncomp > 1) { + // This condition is specific to Mfield, where, + // x-, y-, and z-components are stored on the z-face + mfzfab(i,j,k,0) = xfield_parser(x,y,z); + mfzfab(i,j,k,1) = yfield_parser(x,y,z); + mfzfab(i,j,k,2) = zfield_parser(x,y,z); + } else +#endif + { + mfzfab(i,j,k) = zfield_parser(x,y,z); + } + } + ); + } +} + +void +WarpX::PerformanceHints () +{ + // Check requested MPI ranks and available boxes + amrex::Long total_nboxes = 0; // on all MPI ranks + for (int ilev = 0; ilev <= finestLevel(); ++ilev) { + total_nboxes += boxArray(ilev).size(); + } + auto const nprocs = ParallelDescriptor::NProcs(); + + // Check: are there more MPI ranks than Boxes? + if (nprocs > total_nboxes) { + std::stringstream warnMsg; + warnMsg << "Too many resources / too little work!\n" + << " It looks like you requested more compute resources than " + << "there are total number of boxes of cells available (" + << total_nboxes << "). " + << "You started with (" << nprocs + << ") MPI ranks, so (" << nprocs - total_nboxes + << ") rank(s) will have no work.\n" +#ifdef AMREX_USE_GPU + << " On GPUs, consider using 1-8 boxes per GPU that together fill " + << "each GPU's memory sufficiently. If you do not rely on dynamic " + << "load-balancing, then one large box per GPU is ideal.\n" +#endif + << "Consider decreasing the amr.blocking_factor and " + << "amr.max_grid_size parameters and/or using fewer MPI ranks.\n" + << " More information:\n" + << " https://warpx.readthedocs.io/en/latest/usage/workflows/parallelization.html\n"; + + ablastr::warn_manager::WMRecordWarning( + "Performance", warnMsg.str(), ablastr::warn_manager::WarnPriority::high); + } + +#ifdef AMREX_USE_GPU + // Check: Are there more than 12 boxes per GPU? + if (total_nboxes > nprocs * 12) { + std::stringstream warnMsg; + warnMsg << "Too many boxes per GPU!\n" + << " It looks like you split your simulation domain " + << "in too many boxes (" << total_nboxes << "), which " + << "results in an average number of (" + << amrex::Long(total_nboxes/nprocs) << ") per GPU. " + << "This causes severe overhead in the communication of " + << "border/guard regions.\n" + << " On GPUs, consider using 1-8 boxes per GPU that together fill " + << "each GPU's memory sufficiently. If you do not rely on dynamic " + << "load-balancing, then one large box per GPU is ideal.\n" + << "Consider increasing the amr.blocking_factor and " + << "amr.max_grid_size parameters and/or using more MPI ranks.\n" + << " More information:\n" + << " https://warpx.readthedocs.io/en/latest/usage/workflows/parallelization.html\n"; + + ablastr::warn_manager::WMRecordWarning( + "Performance", warnMsg.str(), ablastr::warn_manager::WarnPriority::high); + } +#endif + + // TODO: warn if some ranks have disproportionally more work than all others + // tricky: it can be ok to assign "vacuum" boxes to some ranks w/o slowing down + // all other ranks; we need to measure this with our load-balancing + // routines and issue a warning only of some ranks stall all other ranks + // TODO: check MPI-rank to GPU ratio (should be 1:1) + // TODO: check memory per MPI rank, especially if GPUs are underutilized + // TODO: CPU tiling hints with OpenMP +} + +void WarpX::CheckGuardCells() +{ + for (int lev = 0; lev <= finest_level; ++lev) + { + for (int dim = 0; dim < 3; ++dim) + { + CheckGuardCells(*Efield_fp[lev][dim]); + CheckGuardCells(*Bfield_fp[lev][dim]); + CheckGuardCells(*current_fp[lev][dim]); + + if (WarpX::fft_do_time_averaging) + { + CheckGuardCells(*Efield_avg_fp[lev][dim]); + CheckGuardCells(*Bfield_avg_fp[lev][dim]); + } + } + + if (rho_fp[lev]) + { + CheckGuardCells(*rho_fp[lev]); + } + + if (F_fp[lev]) + { + CheckGuardCells(*F_fp[lev]); + } + + if (G_fp[lev]) + { + CheckGuardCells(*G_fp[lev]); + } + + // MultiFabs on coarse patch + if (lev > 0) + { + for (int dim = 0; dim < 3; ++dim) + { + CheckGuardCells(*Efield_cp[lev][dim]); + CheckGuardCells(*Bfield_cp[lev][dim]); + CheckGuardCells(*current_cp[lev][dim]); + + if (WarpX::fft_do_time_averaging) + { + CheckGuardCells(*Efield_avg_cp[lev][dim]); + CheckGuardCells(*Bfield_avg_cp[lev][dim]); + } + } + + if (rho_cp[lev]) + { + CheckGuardCells(*rho_cp[lev]); + } + + if (F_cp[lev]) + { + CheckGuardCells(*F_cp[lev]); + } + + if (G_cp[lev]) + { + CheckGuardCells(*G_cp[lev]); + } + } + } +} + +void WarpX::CheckGuardCells(amrex::MultiFab const& mf) +{ + for (amrex::MFIter mfi(mf); mfi.isValid(); ++mfi) + { + const amrex::IntVect vc = mfi.validbox().enclosedCells().size(); + const amrex::IntVect gc = mf.nGrowVect(); + + std::stringstream ss_msg; + ss_msg << "MultiFab " << mf.tags()[1].c_str() << ":" << + " the number of guard cells " << gc << + " is larger than or equal to the number of valid cells " + << vc << ", please reduce the number of guard cells" << + " or increase the grid size by changing domain decomposition."; + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(vc.allGT(gc), ss_msg.str()); + } +} + +void WarpX::InitializeEBGridData (int lev) +{ +#ifdef AMREX_USE_EB + if (lev == maxLevel()) { + + // Throw a warning if EB is on and particle_shape > 1 + bool flag_eb_on = not fieldEBFactory(lev).isAllRegular(); + + if ((nox > 1 or noy > 1 or noz > 1) and flag_eb_on) + { + ablastr::warn_manager::WMRecordWarning("Particles", + "when algo.particle_shape > 1, numerical artifacts will be present when\n" + "particles are close to embedded boundaries"); + } + + if (WarpX::electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD ) { + + auto const eb_fact = fieldEBFactory(lev); + + ComputeEdgeLengths(m_edge_lengths[lev], eb_fact); + ScaleEdges(m_edge_lengths[lev], CellSize(lev)); + ComputeFaceAreas(m_face_areas[lev], eb_fact); + ScaleAreas(m_face_areas[lev], CellSize(lev)); + + if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { + MarkCells(); + ComputeFaceExtensions(); + } + } + + ComputeDistanceToEB(); + + } +#else + amrex::ignore_unused(lev); +#endif +} + +void WarpX::CheckKnownIssues() +{ + if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD && + (std::any_of(do_pml_Lo[0].begin(),do_pml_Lo[0].end(),[](const auto& ee){return ee;}) || + std::any_of(do_pml_Hi[0].begin(),do_pml_Hi[0].end(),[](const auto& ee){return ee;})) ) + { + ablastr::warn_manager::WMRecordWarning( + "PML", + "Using PSATD together with PML may lead to instabilities if the plasma touches the PML region. " + "It is recommended to leave enough empty space between the plasma boundary and the PML region.", + ablastr::warn_manager::WarnPriority::low); + } +} + +#if defined(WARPX_USE_OPENPMD) && !defined(WARPX_DIM_1D_Z) && !defined(WARPX_DIM_XZ) +void +WarpX::ReadExternalFieldFromFile ( + std::string read_fields_from_path, amrex::MultiFab* mf, + std::string F_name, std::string F_component) +{ + // Get WarpX domain info + auto& warpx = WarpX::GetInstance(); + amrex::Geometry const& geom0 = warpx.Geom(0); + const amrex::RealBox& real_box = geom0.ProbDomain(); + const auto dx = geom0.CellSizeArray(); + amrex::IntVect nodal_flag = mf->ixType().toIntVect(); + + // Read external field openPMD data + auto series = openPMD::Series(read_fields_from_path, openPMD::Access::READ_ONLY); + auto iseries = series.iterations.begin()->second; + auto F = iseries.meshes[F_name]; + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(F.getAttribute("dataOrder").get() == "C", + "Reading from files with non-C dataOrder is not implemented"); + + auto axisLabels = F.getAttribute("axisLabels").get>(); + auto fileGeom = F.getAttribute("geometry").get(); + +#if defined(WARPX_DIM_3D) + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(fileGeom == "cartesian", "3D can only read from files with cartesian geometry"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(axisLabels[0] == "x" && axisLabels[1] == "y" && axisLabels[2] == "z", + "3D expects axisLabels {x, y, z}"); +#elif defined(WARPX_DIM_XZ) + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(fileGeom == "cartesian", "XZ can only read from files with cartesian geometry"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(axisLabels[0] == "x" && axisLabels[1] == "z", + "XZ expects axisLabels {x, z}"); +#elif defined(WARPX_DIM_1D_Z) + WARPX_ABORT_WITH_MESSAGE( + "Reading from openPMD for external fields is not known to work with 1D3V (see #3830)"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(fileGeom == "cartesian", "1D3V can only read from files with cartesian geometry"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(axisLabels[0] == "z"); +#elif defined(WARPX_DIM_RZ) + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(fileGeom == "thetaMode", "RZ can only read from files with 'thetaMode' geometry"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(axisLabels[0] == "r" && axisLabels[1] == "z", + "RZ expects axisLabels {r, z}"); +#endif + + auto offset = F.gridGlobalOffset(); + amrex::Real offset0 = offset[0]; + amrex::Real offset1 = offset[1]; +#if defined(WARPX_DIM_3D) + amrex::Real offset2 = offset[2]; +#endif + auto d = F.gridSpacing(); + +#if defined(WARPX_DIM_RZ) + amrex::Real file_dr = d[0]; + amrex::Real file_dz = d[1]; +#elif defined(WARPX_DIM_3D) + amrex::Real file_dx = d[0]; + amrex::Real file_dy = d[1]; + amrex::Real file_dz = d[2]; +#endif + + auto FC = F[F_component]; + auto extent = FC.getExtent(); + int extent0 = extent[0]; + int extent1 = extent[1]; + int extent2 = extent[2]; + + // Determine the chunk data that will be loaded. + // Now, the full range of data is loaded. + // Loading chunk data can speed up the process. + // Thus, `chunk_offset` and `chunk_extent` should be modified accordingly in another PR. + openPMD::Offset chunk_offset = {0,0,0}; + openPMD::Extent chunk_extent = {extent[0], extent[1], extent[2]}; + + auto FC_chunk_data = FC.loadChunk(chunk_offset,chunk_extent); + series.flush(); + auto FC_data_host = FC_chunk_data.get(); + + // Load data to GPU + size_t total_extent = size_t(extent[0]) * extent[1] * extent[2]; + amrex::Gpu::DeviceVector FC_data_gpu(total_extent); + auto FC_data = FC_data_gpu.data(); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, FC_data_host, FC_data_host + total_extent, FC_data); + + // Loop over boxes + for (MFIter mfi(*mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + amrex::Box box = mfi.growntilebox(); + amrex::Box tb = mfi.tilebox(nodal_flag, mf->nGrowVect()); + auto const& mffab = mf->array(mfi); + + // Start ParallelFor + amrex::ParallelFor (tb, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + // i,j,k denote x,y,z indices in 3D xyz. + // i,j denote r,z indices in 2D rz; k is just 0 + + // ii is used for 2D RZ mode + int ii = i; +#if defined(WARPX_DIM_RZ) + // In 2D RZ, i denoting r can be < 0 + // but mirrored values should be assigned. + // Namely, mffab(i) = FC_data[-i] when i<0. + if (i<0) {ii = -i;} +#endif + + // Physical coordinates of the grid point + // 0,1,2 denote x,y,z in 3D xyz. + // 0,1 denote r,z in 2D rz. + amrex::Real x0, x1; + if ( box.type(0)==amrex::IndexType::CellIndex::NODE ) + { x0 = real_box.lo(0) + ii*dx[0]; } + else { x0 = real_box.lo(0) + ii*dx[0] + 0.5*dx[0]; } + if ( box.type(1)==amrex::IndexType::CellIndex::NODE ) + { x1 = real_box.lo(1) + j*dx[1]; } + else { x1 = real_box.lo(1) + j*dx[1] + 0.5*dx[1]; } + +#if defined(WARPX_DIM_RZ) + // Get index of the external field array + int const ir = floor( (x0-offset0)/file_dr ); + int const iz = floor( (x1-offset1)/file_dz ); + + // Get coordinates of external grid point + amrex::Real const xx0 = offset0 + ir * file_dr; + amrex::Real const xx1 = offset1 + iz * file_dz; + +#elif defined(WARPX_DIM_3D) + amrex::Real x2; + if ( box.type(2)==amrex::IndexType::CellIndex::NODE ) + { x2 = real_box.lo(2) + k*dx[2]; } + else { x2 = real_box.lo(2) + k*dx[2] + 0.5*dx[2]; } + + // Get index of the external field array + int const ix = floor( (x0-offset0)/file_dx ); + int const iy = floor( (x1-offset1)/file_dy ); + int const iz = floor( (x2-offset2)/file_dz ); + + // Get coordinates of external grid point + amrex::Real const xx0 = offset0 + ix * file_dx; + amrex::Real const xx1 = offset1 + iy * file_dy; + amrex::Real const xx2 = offset2 + iz * file_dz; +#endif + +#if defined(WARPX_DIM_RZ) + amrex::Array4 fc_array(FC_data, {0,0,0}, {extent0, extent2, extent1}, 1); + double + f00 = fc_array(0, iz , ir ), + f01 = fc_array(0, iz , ir+1), + f10 = fc_array(0, iz+1, ir ), + f11 = fc_array(0, iz+1, ir+1); + mffab(i,j,k) = utils::algorithms::bilinear_interp + (xx0, xx0+file_dr, xx1, xx1+file_dz, + f00, f01, f10, f11, + x0, x1); +#elif defined(WARPX_DIM_3D) + amrex::Array4 fc_array(FC_data, {0,0,0}, {extent2, extent1, extent0}, 1); + double + f000 = fc_array(iz , iy , ix ), + f001 = fc_array(iz+1, iy , ix ), + f010 = fc_array(iz , iy+1, ix ), + f011 = fc_array(iz+1, iy+1, ix ), + f100 = fc_array(iz , iy , ix+1), + f101 = fc_array(iz+1, iy , ix+1), + f110 = fc_array(iz , iy+1, ix+1), + f111 = fc_array(iz+1, iy+1, ix+1); + mffab(i,j,k) = utils::algorithms::trilinear_interp + (xx0, xx0+file_dx, xx1, xx1+file_dy, xx2, xx2+file_dz, + f000, f001, f010, f011, f100, f101, f110, f111, + x0, x1, x2); +#endif + + } + + ); // End ParallelFor + + } // End loop over boxes. + +} // End function WarpX::ReadExternalFieldFromFile +#else // WARPX_USE_OPENPMD && !WARPX_DIM_1D_Z && !defined(WARPX_DIM_XZ) +void +WarpX::ReadExternalFieldFromFile (std::string , amrex::MultiFab* ,std::string, std::string) +{ +#if defined(WARPX_DIM_1D) + WARPX_ABORT_WITH_MESSAGE("Reading fields from openPMD files is not supported in 1D"); +#elif defined(WARPX_DIM_XZ) + WARPX_ABORT_WITH_MESSAGE("Reading from openPMD for external fields is not known to work with XZ (see #3828)"); +#elif !defined(WARPX_USE_OPENPMD) + WARPX_ABORT_WITH_MESSAGE("OpenPMD field reading requires OpenPMD support to be enabled"); +#endif +} +#endif // WARPX_USE_OPENPMD diff --git a/Source/WarpX.H b/Source/WarpX.H index 870c7eeb5..afb335463 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -40,6 +40,7 @@ #include "Utils/WarpXAlgorithmSelection.H" #include "FieldSolver/London/London.H" #include "FieldSolver/LumpedElement/Inductor.H" +#include "FieldSolver/LumpedElement/JosephsonJunction.H" #include #include @@ -102,6 +103,7 @@ public: MacroscopicProperties& GetMacroscopicProperties () { return *m_macroscopic_properties; } London& getLondon () { return *m_london; } Inductor& getInductor () { return *m_inductor; } + JosephsonJunction& getJosephsonJunction () { return *m_jj; } MultiDiagnostics& GetMultiDiags () {return *multi_diags;} ParticleBoundaryBuffer& GetParticleBoundaryBuffer () { return *m_particle_boundary_buffer; } @@ -320,6 +322,9 @@ public: static int use_PEC_mask; static int use_lumped_inductor; + static int use_lumped_resistor; + static int use_lumped_capacitor; + static int use_josephson_junction; //! Integer that corresponds to the order of the PSATD solution //! (whether the PSATD equations are derived from first-order or @@ -1742,6 +1747,8 @@ private: std::unique_ptr m_london; // Lumped inductor std::unique_ptr m_inductor; + // Josephson junction (nonlinear inductor) + std::unique_ptr m_jj; #ifdef WARPX_MAG_LLG // time advancement scheme of M field diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 46d425e94..82f440adb 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -225,6 +225,9 @@ amrex::Vector WarpX::particle_boundary_hi(AMREX_SPACEDIM,P int WarpX::yee_coupled_solver_algo; int WarpX::use_PEC_mask = 0; int WarpX::use_lumped_inductor = 0; +int WarpX::use_lumped_resistor = 0; +int WarpX::use_lumped_capacitor = 0; +int WarpX::use_josephson_junction = 0; bool WarpX::do_current_centering = false; @@ -495,6 +498,11 @@ WarpX::WarpX () m_inductor = std::make_unique(); } + // Josephson junction (nonlinear inductor) + if (use_josephson_junction) { + m_jj = std::make_unique(); + } + // Set default values for particle and cell weights for costs update; // Default values listed here for the case AMREX_USE_GPU are determined // from single-GPU tests on Summit. @@ -1033,6 +1041,9 @@ WarpX::ReadParameters () ); } + pp_warpx.query("use_lumped_resistor", use_lumped_resistor); + pp_warpx.query("use_lumped_capacitor", use_lumped_capacitor); + #ifdef WARPX_MAG_LLG // Read the value of the time advancement scheme of M field pp_warpx.query("mag_time_scheme_order", mag_time_scheme_order); @@ -1267,6 +1278,7 @@ WarpX::ReadParameters () pp_algo.query("use_PEC_mask",use_PEC_mask); pp_algo.query("use_lumped_inductor",use_lumped_inductor); + pp_algo.query("use_josephson_junction",use_josephson_junction); // Load balancing parameters std::vector load_balance_intervals_string_vec = {"0"};