diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index a3ca703d15e..cf9d259d45a 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -1732,8 +1732,8 @@ are applied to the grid directly. In particular, these fields can be seen in the One can refer to input files in ``Examples/Tests/LoadExternalField`` for more information. Regarding how to prepare the openPMD data file, one can refer to the `openPMD-example-datasets `__. - Note that if both `B_ext_grid_init_style` and `E_ext_grid_init_style` are set to - `read_from_file`, the openPMD file specified by `warpx.read_fields_from_path` + Note that if both ``B_ext_grid_init_style`` and ``E_ext_grid_init_style`` are set to + ``read_from_file``, the openPMD file specified by ``warpx.read_fields_from_path`` should contain both B and E external fields data. * ``warpx.E_external_grid`` & ``warpx.B_external_grid`` (list of `3 floats`) @@ -1789,6 +1789,21 @@ are applied to the particles directly, at each timestep. As a results, these fie Note that the position is defined in Cartesian coordinates, as a function of (x,y,z), even for RZ. + * ``read_from_file``: load the external field from an openPMD file. + An additional parameter, indicating the path of an openPMD data file, ``warpx.read_fields_from_path`` + must be specified, from which the external E field data can be loaded into WarpX. + One can refer to input files in ``Examples/Tests/LoadExternalField`` for more information. + Regarding how to prepare the openPMD data file, one can refer to + the `openPMD-example-datasets `__. + Note that if both ``B_ext_particle_init_style`` and ``E_ext_particle_init_style`` are set to + ``read_from_file``, the openPMD file specified by ``warpx.read_fields_from_path`` + should contain both B and E external fields data. + + .. note:: + + When using ``read_from_file``, the fields loaded from the file will be interpolated + to the resolution of the grid used for the simulation. + * ``repeated_plasma_lens``: apply a series of plasma lenses. The properties of the lenses are defined in the lab frame by the input parameters: diff --git a/Examples/Tests/LoadExternalField/inputs_rz b/Examples/Tests/LoadExternalField/inputs_rz_grid_fields similarity index 100% rename from Examples/Tests/LoadExternalField/inputs_rz rename to Examples/Tests/LoadExternalField/inputs_rz_grid_fields diff --git a/Examples/Tests/LoadExternalField/inputs_rz_particle_fields b/Examples/Tests/LoadExternalField/inputs_rz_particle_fields new file mode 100644 index 00000000000..c5ddaecbddf --- /dev/null +++ b/Examples/Tests/LoadExternalField/inputs_rz_particle_fields @@ -0,0 +1,65 @@ +warpx.serialize_initial_conditions = 0 +warpx.do_dynamic_scheduling = 0 +particles.do_tiling = 0 + +particles.B_ext_particle_init_style = "read_from_file" +warpx.read_fields_from_path = "../../../../openPMD-example-datasets/example-femm-thetaMode.h5" + +warpx.grid_type = collocated +warpx.do_electrostatic = labframe + +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 300 +amr.n_cell = 40 40 +warpx.numprocs = 1 1 +amr.max_level = 0 +geometry.dims = RZ + +geometry.prob_lo = 0.0 0.0 +geometry.prob_hi = 1.0 5.0 + +################################# +###### Boundary Condition ####### +################################# +boundary.field_lo = none pec +boundary.field_hi = pec pec +boundary.potential_lo_x = 0 +boundary.potential_hi_x = 0 +boundary.potential_lo_y = 0 +boundary.potential_hi_y = 0 +boundary.potential_lo_z = 0 +boundary.potential_hi_z = 0 + +################################# +############ NUMERICS ########### +################################# +warpx.serialize_initial_conditions = 1 +warpx.verbose = 1 +warpx.const_dt = 4.40917904849092e-7 +warpx.use_filter = 0 + +# Order of particle shape factors +algo.particle_shape = 1 + +################################# +############ PLASMA ############# +################################# +particles.species_names = proton +proton.injection_style = "SingleParticle" +proton.single_particle_pos = 0.0 0.2 2.5 +proton.single_particle_u = 9.506735958279367e-05 0.0 0.00013435537232359165 +proton.single_particle_weight = 1.0 +proton.do_not_deposit = 1 +proton.mass = m_p +proton.charge = q_e + +# Diagnostics +diagnostics.diags_names = diag1 chk +diag1.intervals = 300 +diag1.diag_type = Full + +chk.intervals = 150 +chk.diag_type = Full +chk.format = checkpoint diff --git a/Regression/Checksum/benchmarks_json/LoadExternalFieldRZ.json b/Regression/Checksum/benchmarks_json/LoadExternalFieldRZGrid.json similarity index 100% rename from Regression/Checksum/benchmarks_json/LoadExternalFieldRZ.json rename to Regression/Checksum/benchmarks_json/LoadExternalFieldRZGrid.json diff --git a/Regression/Checksum/benchmarks_json/LoadExternalFieldRZParticles.json b/Regression/Checksum/benchmarks_json/LoadExternalFieldRZParticles.json new file mode 100644 index 00000000000..dc4ffa3beb2 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/LoadExternalFieldRZParticles.json @@ -0,0 +1,22 @@ +{ + "lev=0": { + "Br": 0.6380326770459374, + "Bt": 0.0, + "Bz": 4.850698699951099, + "Er": 0.0, + "Et": 0.0, + "Ez": 0.0, + "jr": 0.0, + "jt": 0.0, + "jz": 0.0 + }, + "proton": { + "particle_momentum_x": 1.2945379007125463e-23, + "particle_momentum_y": 7.901629565307941e-23, + "particle_momentum_z": 2.0004592432172922e-23, + "particle_position_x": 0.12402004585603355, + "particle_position_y": 4.363249204658946, + "particle_theta": 0.22200801011337173, + "particle_weight": 1.0 + } +} \ No newline at end of file diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 41db9a15bdc..259fa90bb3f 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -2266,9 +2266,28 @@ doVis = 0 compareParticles = 0 analysisRoutine = Examples/Tests/resampling/analysis_leveling_thinning.py -[LoadExternalFieldRZ] +[LoadExternalFieldRZGrid] buildDir = . -inputFile = Examples/Tests/LoadExternalField/inputs_rz +inputFile = Examples/Tests/LoadExternalField/inputs_rz_grid_fields +runtime_params = warpx.abort_on_warning_threshold=medium chk.file_prefix=LoadExternalFieldRZ_chk chk.file_min_digits=5 +dim = 2 +addToCompileString = USE_RZ=TRUE +cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_OPENPMD=ON +restartTest = 1 +restartFileNum = 150 +useMPI = 1 +numprocs = 1 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = proton +analysisRoutine = Examples/Tests/LoadExternalField/analysis_rz.py + +[LoadExternalFieldRZParticles] +buildDir = . +inputFile = Examples/Tests/LoadExternalField/inputs_rz_particle_fields runtime_params = warpx.abort_on_warning_threshold=medium chk.file_prefix=LoadExternalFieldRZ_chk chk.file_min_digits=5 dim = 2 addToCompileString = USE_RZ=TRUE diff --git a/Source/Initialization/ExternalField.cpp b/Source/Initialization/ExternalField.cpp index e7a272e2e95..c46a9b0c178 100644 --- a/Source/Initialization/ExternalField.cpp +++ b/Source/Initialization/ExternalField.cpp @@ -172,11 +172,6 @@ ExternalFieldParams::ExternalFieldParams(const amrex::ParmParse& pp_warpx) // // External fields from file // - - if (E_ext_grid_type == ExternalFieldType::read_from_file || - B_ext_grid_type == ExternalFieldType::read_from_file){ - const std::string read_fields_from_path="./"; - pp_warpx.query("read_fields_from_path", external_fields_path); - } + pp_warpx.query("read_fields_from_path", external_fields_path); //___________________________________________________________________________ } diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 68b38153f8d..c3ada061f65 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1359,7 +1359,10 @@ void WarpX::CheckKnownIssues() void WarpX::LoadExternalFieldsFromFile (int const lev) { - if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) { + // TODO: assert that we cannot have both external grid and particle field from file + + if ( (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) || + (mypc->m_B_ext_particle_s == "read_from_file") ){ #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)"); @@ -1372,7 +1375,8 @@ WarpX::LoadExternalFieldsFromFile (int const lev) ReadExternalFieldFromFile(m_p_ext_field_params->external_fields_path, Bfield_fp_external[lev][2].get(), "B", "z"); #endif } - if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) { + if ( (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) || + (mypc->m_E_ext_particle_s == "read_from_file") ){ #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)"); diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index dcf5fbe5621..1145256b8a1 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -18,6 +18,7 @@ #include "Utils/WarpXProfilerWrapper.H" #include "WarpXComm_K.H" #include "WarpXSumGuardCells.H" +#include "Particles/MultiParticleContainer.H" #include #include @@ -59,6 +60,24 @@ WarpX::UpdateAuxilaryData () } else { UpdateAuxilaryDataStagToNodal(); } + + // When loading particle fields from file: add the external fields: + // TODO reuse WarpX::AddExternalFields + // TODO replace fp_external with simply external + for (int lev = 0; lev <= finest_level; ++lev) { + // FIXME: RZ multimode has more than one component for all these + if (mypc->m_E_ext_particle_s == "read_from_file") { + amrex::MultiFab::Add(*Efield_aux[lev][0], *Efield_fp_external[lev][0], 0, 0, 1, guard_cells.ng_alloc_EB); + amrex::MultiFab::Add(*Efield_aux[lev][1], *Efield_fp_external[lev][1], 0, 0, 1, guard_cells.ng_alloc_EB); + amrex::MultiFab::Add(*Efield_aux[lev][2], *Efield_fp_external[lev][2], 0, 0, 1, guard_cells.ng_alloc_EB); + } + if (mypc->m_B_ext_particle_s == "read_from_file") { + amrex::MultiFab::Add(*Bfield_aux[lev][0], *Bfield_fp_external[lev][0], 0, 0, 1, guard_cells.ng_alloc_EB); + amrex::MultiFab::Add(*Bfield_aux[lev][1], *Bfield_fp_external[lev][1], 0, 0, 1, guard_cells.ng_alloc_EB); + amrex::MultiFab::Add(*Bfield_aux[lev][2], *Bfield_fp_external[lev][2], 0, 0, 1, guard_cells.ng_alloc_EB); + } + } + } void @@ -293,6 +312,18 @@ WarpX::UpdateAuxilaryDataStagToNodal () void WarpX::UpdateAuxilaryDataSameType () { + // Update aux field, including guard cells, up to ng_FieldGather + const amrex::IntVect& ng_src = guard_cells.ng_FieldGather; + + // Level 0: Copy from fine to aux + // TODO Check if aux and fp are aliases of each other. If so, we can skip the copy. + MultiFab::Copy(*Efield_aux[0][0], *Efield_fp[0][0], 0, 0, Efield_aux[0][0]->nComp(), ng_src); + MultiFab::Copy(*Efield_aux[0][1], *Efield_fp[0][1], 0, 0, Efield_aux[0][1]->nComp(), ng_src); + MultiFab::Copy(*Efield_aux[0][2], *Efield_fp[0][2], 0, 0, Efield_aux[0][2]->nComp(), ng_src); + MultiFab::Copy(*Bfield_aux[0][0], *Bfield_fp[0][0], 0, 0, Bfield_aux[0][0]->nComp(), ng_src); + MultiFab::Copy(*Bfield_aux[0][1], *Bfield_fp[0][1], 0, 0, Bfield_aux[0][1]->nComp(), ng_src); + MultiFab::Copy(*Bfield_aux[0][2], *Bfield_fp[0][2], 0, 0, Bfield_aux[0][2]->nComp(), ng_src); + for (int lev = 1; lev <= finest_level; ++lev) { const amrex::Periodicity& crse_period = Geom(lev-1).periodicity(); @@ -308,8 +339,6 @@ WarpX::UpdateAuxilaryDataSameType () dBy.setVal(0.0); dBz.setVal(0.0); - // Guard cells may not be up to date beyond ng_FieldGather - const amrex::IntVect& ng_src = guard_cells.ng_FieldGather; // Copy Bfield_aux to the dB MultiFabs, using up to ng_src (=ng_FieldGather) guard // cells from Bfield_aux and filling up to ng (=nGrow) guard cells in the dB MultiFabs @@ -379,8 +408,6 @@ WarpX::UpdateAuxilaryDataSameType () dEy.setVal(0.0); dEz.setVal(0.0); - // Guard cells may not be up to date beyond ng_FieldGather - const amrex::IntVect& ng_src = guard_cells.ng_FieldGather; // Copy Efield_aux to the dE MultiFabs, using up to ng_src (=ng_FieldGather) guard // cells from Efield_aux and filling up to ng (=nGrow) guard cells in the dE MultiFabs ablastr::utils::communication::ParallelCopy(dEx, *Efield_aux[lev - 1][0], 0, 0, diff --git a/Source/Particles/Gather/GetExternalFields.H b/Source/Particles/Gather/GetExternalFields.H index 7000d6d7c26..36c88989d0c 100644 --- a/Source/Particles/Gather/GetExternalFields.H +++ b/Source/Particles/Gather/GetExternalFields.H @@ -23,7 +23,7 @@ */ struct GetExternalEBField { - enum ExternalFieldInitType { None, Parser, RepeatedPlasmaLens, Unknown }; + enum ExternalFieldInitType { None, Parser, RepeatedPlasmaLens, FromFile, Unknown }; GetExternalEBField () = default; diff --git a/Source/Particles/Gather/GetExternalFields.cpp b/Source/Particles/Gather/GetExternalFields.cpp index 8e3c679aa3a..9ef7d407667 100644 --- a/Source/Particles/Gather/GetExternalFields.cpp +++ b/Source/Particles/Gather/GetExternalFields.cpp @@ -83,6 +83,13 @@ GetExternalEBField::GetExternalEBField (const WarpXParIter& a_pti, long a_offset m_repeated_plasma_lens_strengths_B = mypc.d_repeated_plasma_lens_strengths_B.data(); } + if (mypc.m_E_ext_particle_s == "read_from_file") { + m_Etype = ExternalFieldInitType::FromFile; + } + if (mypc.m_B_ext_particle_s == "read_from_file") { + m_Btype = ExternalFieldInitType::FromFile; + } + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_Etype != Unknown, "Unknown E_ext_particle_init_style"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_Btype != Unknown, "Unknown B_ext_particle_init_style"); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index ef81aef4482..dced5c77eba 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -2308,12 +2308,14 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm AllocInitMultiFab(current_fp[lev][2], amrex::convert(ba, jz_nodal_flag), dm, ncomps, ngJ, lev, "current_fp[z]", 0.0_rt); // Match external field MultiFabs to fine patch - if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) { + if ( (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) || + (mypc->m_B_ext_particle_s == "read_from_file") ){ AllocInitMultiFab(Bfield_fp_external[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, lev, "Bfield_fp_external[x]", 0.0_rt); AllocInitMultiFab(Bfield_fp_external[lev][1], amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, lev, "Bfield_fp_external[y]", 0.0_rt); AllocInitMultiFab(Bfield_fp_external[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, lev, "Bfield_fp_external[z]", 0.0_rt); } - if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) { + if ( (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) || + (mypc->m_E_ext_particle_s == "read_from_file") ){ AllocInitMultiFab(Efield_fp_external[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, lev, "Efield_fp_external[x]", 0.0_rt); AllocInitMultiFab(Efield_fp_external[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, lev, "Efield_fp_external[y]", 0.0_rt); AllocInitMultiFab(Efield_fp_external[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, lev, "Efield_fp_external[z]", 0.0_rt); @@ -2550,23 +2552,35 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm AllocInitMultiFab(Efield_aux[lev][1], nba, dm, ncomps, ngEB, lev, "Efield_aux[y]", 0.0_rt); AllocInitMultiFab(Efield_aux[lev][2], nba, dm, ncomps, ngEB, lev, "Efield_aux[z]", 0.0_rt); } else if (lev == 0) { - if (!WarpX::fft_do_time_averaging) { - // In this case, the aux grid is simply an alias of the fp grid - AliasInitMultiFab(Efield_aux[lev][0], *Efield_fp[lev][0], 0, ncomps, lev, "Efield_aux[x]", 0.0_rt); - AliasInitMultiFab(Efield_aux[lev][1], *Efield_fp[lev][1], 0, ncomps, lev, "Efield_aux[y]", 0.0_rt); - AliasInitMultiFab(Efield_aux[lev][2], *Efield_fp[lev][2], 0, ncomps, lev, "Efield_aux[z]", 0.0_rt); - - AliasInitMultiFab(Bfield_aux[lev][0], *Bfield_fp[lev][0], 0, ncomps, lev, "Bfield_aux[x]", 0.0_rt); - AliasInitMultiFab(Bfield_aux[lev][1], *Bfield_fp[lev][1], 0, ncomps, lev, "Bfield_aux[y]", 0.0_rt); - AliasInitMultiFab(Bfield_aux[lev][2], *Bfield_fp[lev][2], 0, ncomps, lev, "Bfield_aux[z]", 0.0_rt); - } else { - AliasInitMultiFab(Efield_aux[lev][0], *Efield_avg_fp[lev][0], 0, ncomps, lev, "Efield_aux[x]", 0.0_rt); - AliasInitMultiFab(Efield_aux[lev][1], *Efield_avg_fp[lev][1], 0, ncomps, lev, "Efield_aux[y]", 0.0_rt); - AliasInitMultiFab(Efield_aux[lev][2], *Efield_avg_fp[lev][2], 0, ncomps, lev, "Efield_aux[z]", 0.0_rt); - + if (WarpX::fft_do_time_averaging) { AliasInitMultiFab(Bfield_aux[lev][0], *Bfield_avg_fp[lev][0], 0, ncomps, lev, "Bfield_aux[x]", 0.0_rt); AliasInitMultiFab(Bfield_aux[lev][1], *Bfield_avg_fp[lev][1], 0, ncomps, lev, "Bfield_aux[y]", 0.0_rt); AliasInitMultiFab(Bfield_aux[lev][2], *Bfield_avg_fp[lev][2], 0, ncomps, lev, "Bfield_aux[z]", 0.0_rt); + + AliasInitMultiFab(Efield_aux[lev][0], *Efield_avg_fp[lev][0], 0, ncomps, lev, "Efield_aux[x]", 0.0_rt); + AliasInitMultiFab(Efield_aux[lev][1], *Efield_avg_fp[lev][1], 0, ncomps, lev, "Efield_aux[y]", 0.0_rt); + AliasInitMultiFab(Efield_aux[lev][2], *Efield_avg_fp[lev][2], 0, ncomps, lev, "Efield_aux[z]", 0.0_rt); + } else { + if (mypc->m_B_ext_particle_s == "read_from_file") { + AllocInitMultiFab(Bfield_aux[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, lev, "Bfield_aux[x]"); + AllocInitMultiFab(Bfield_aux[lev][1], amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, lev, "Bfield_aux[y]"); + AllocInitMultiFab(Bfield_aux[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, lev, "Bfield_aux[z]"); + } else { + // In this case, the aux grid is simply an alias of the fp grid (most common case in WarpX) + AliasInitMultiFab(Bfield_aux[lev][0], *Bfield_fp[lev][0], 0, ncomps, lev, "Bfield_aux[x]", 0.0_rt); + AliasInitMultiFab(Bfield_aux[lev][1], *Bfield_fp[lev][1], 0, ncomps, lev, "Bfield_aux[y]", 0.0_rt); + AliasInitMultiFab(Bfield_aux[lev][2], *Bfield_fp[lev][2], 0, ncomps, lev, "Bfield_aux[z]", 0.0_rt); + } + if (mypc->m_E_ext_particle_s == "read_from_file") { + AllocInitMultiFab(Efield_aux[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, lev, "Efield_aux[x]"); + AllocInitMultiFab(Efield_aux[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, lev, "Efield_aux[y]"); + AllocInitMultiFab(Efield_aux[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, lev, "Efield_aux[z]"); + } else { + // In this case, the aux grid is simply an alias of the fp grid (most common case in WarpX) + AliasInitMultiFab(Efield_aux[lev][0], *Efield_fp[lev][0], 0, ncomps, lev, "Efield_aux[x]", 0.0_rt); + AliasInitMultiFab(Efield_aux[lev][1], *Efield_fp[lev][1], 0, ncomps, lev, "Efield_aux[y]", 0.0_rt); + AliasInitMultiFab(Efield_aux[lev][2], *Efield_fp[lev][2], 0, ncomps, lev, "Efield_aux[z]", 0.0_rt); + } } } else { AllocInitMultiFab(Bfield_aux[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, lev, "Bfield_aux[x]");