Add ability to load external particle fields from file

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 <https://github.com/openPMD/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 <https://github.com/openPMD/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:
 

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

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
+  }
+}

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

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);
     //___________________________________________________________________________
 }

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)");

Source/Parallelization/WarpXComm.cpp

@@ -18,6 +18,7 @@
 #include "Utils/WarpXProfilerWrapper.H"
 #include "WarpXComm_K.H"
 #include "WarpXSumGuardCells.H"
+#include "Particles/MultiParticleContainer.H"
 
 #include <ablastr/coarsen/average.H>
 #include <ablastr/utils/Communication.H>
@@ -59,6 +60,24 @@
     } 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,29 @@
 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.
+    if (WarpX::fft_do_time_averaging)
+    {
+        MultiFab::Copy(*Efield_aux[0][0], *Efield_avg_fp[0][0], 0, 0, Efield_aux[0][0]->nComp(), ng_src);
+        MultiFab::Copy(*Efield_aux[0][1], *Efield_avg_fp[0][1], 0, 0, Efield_aux[0][1]->nComp(), ng_src);
+        MultiFab::Copy(*Efield_aux[0][2], *Efield_avg_fp[0][2], 0, 0, Efield_aux[0][2]->nComp(), ng_src);
+        MultiFab::Copy(*Bfield_aux[0][0], *Bfield_avg_fp[0][0], 0, 0, Bfield_aux[0][0]->nComp(), ng_src);
+        MultiFab::Copy(*Bfield_aux[0][1], *Bfield_avg_fp[0][1], 0, 0, Bfield_aux[0][1]->nComp(), ng_src);
+        MultiFab::Copy(*Bfield_aux[0][2], *Bfield_avg_fp[0][2], 0, 0, Bfield_aux[0][2]->nComp(), ng_src);
+    }
+    else
+    {
+        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 +350,6 @@
             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 +419,6 @@
             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,

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;
 

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");
 

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]");