Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Add ability to load external particle fields from file #4981

Closed
wants to merge 11 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 17 additions & 2 deletions Docs/source/usage/parameters.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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`)
Expand Down Expand Up @@ -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:

Expand Down
65 changes: 65 additions & 0 deletions Examples/Tests/LoadExternalField/inputs_rz_particle_fields
Original file line number Diff line number Diff line change
@@ -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
Original file line number Diff line number Diff line change
@@ -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
}
}
23 changes: 21 additions & 2 deletions Regression/WarpX-tests.ini
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
7 changes: 1 addition & 6 deletions Source/Initialization/ExternalField.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
//___________________________________________________________________________
}
8 changes: 6 additions & 2 deletions Source/Initialization/WarpXInitData.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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)");
Expand All @@ -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)");
Expand Down
35 changes: 31 additions & 4 deletions Source/Parallelization/WarpXComm.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down Expand Up @@ -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

Check notice

Code scanning / CodeQL

FIXME comment Note

FIXME comment: 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
Expand Down Expand Up @@ -293,6 +312,18 @@
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.
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not needed anymore after this PR is merged:
AMReX-Codes/amrex#3986

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);
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here, we need to check whether WarpX::fft_do_time_averaging is true. If yes,

MultiFab::Copy(*Bfield_aux[0][0], *Bfield_avg_fp[0][0], 0, 0, Bfield_aux[0][0]->nComp(), ng_src);

and otherwise, keep this code


for (int lev = 1; lev <= finest_level; ++lev)
{
const amrex::Periodicity& crse_period = Geom(lev-1).periodicity();
Expand All @@ -308,8 +339,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

Expand Down Expand Up @@ -379,8 +408,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,
Expand Down
2 changes: 1 addition & 1 deletion Source/Particles/Gather/GetExternalFields.H
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
*/
struct GetExternalEBField
{
enum ExternalFieldInitType { None, Parser, RepeatedPlasmaLens, Unknown };
enum ExternalFieldInitType { None, Parser, RepeatedPlasmaLens, FromFile, Unknown };

GetExternalEBField () = default;

Expand Down
7 changes: 7 additions & 0 deletions Source/Particles/Gather/GetExternalFields.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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");

Expand Down
46 changes: 30 additions & 16 deletions Source/WarpX.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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]");
Expand Down
Loading