Skip to content

Commit

Permalink
Treat external field types consistently (#5104)
Browse files Browse the repository at this point in the history 
* use external field multifabs consistently for different external field types

* make external field MF vectors length equal to `nlevs_max`

* avoid adding external grid fields multiple times

* load parse ext grid fields after restart

* rename `LoadExternalFieldsFromFile` to `LoadExternalFields`

* Update inline comment on why external fields are added to ES solutions

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

---------

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>
  • Loading branch information
@roelof-groenewald @RemiLehe
roelof-groenewald and RemiLehe authored Aug 2, 2024
1 parent 6d54864 commit f1f3c38
Show file tree
Hide file tree
Showing 4 changed files with 82 additions and 70 deletions.
7 changes: 6 additions & 1 deletion Source/Evolve/WarpXEvolve.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -246,7 +246,12 @@ WarpX::Evolve (int numsteps)
// This is currently a lab frame calculation.
ComputeMagnetostaticField();
}
AddExternalFields();
// Since the fields were reset above, the external fields are added
// back on to the fine patch fields. This make it so that the net fields
// are the sum of the field solution and any external field.
for (int lev = 0; lev <= max_level; ++lev) {
AddExternalFields(lev);
}
} else if (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) {
// Hybrid-PIC case:
// The particles are now at p^{n+1/2} and x^{n+1}. The fields
Expand Down
131 changes: 69 additions & 62 deletions Source/Initialization/WarpXInitData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -513,12 +513,6 @@ WarpX::InitData ()
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();
}

if (restart_chkfile.empty() || write_diagnostics_on_restart) {
Expand All @@ -539,15 +533,27 @@ WarpX::InitData ()
}

void
WarpX::AddExternalFields () {
for (int lev = 0; lev <= finest_level; ++lev) {
// FIXME: RZ multimode has more than one component for all these
if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) {
WarpX::AddExternalFields (int const lev) {
// FIXME: RZ multimode has more than one component for all these
if (m_p_ext_field_params->E_ext_grid_type != ExternalFieldType::default_zero) {
if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::constant) {
Efield_fp[lev][0]->plus(m_p_ext_field_params->E_external_grid[0], guard_cells.ng_alloc_EB.min());
Efield_fp[lev][1]->plus(m_p_ext_field_params->E_external_grid[1], guard_cells.ng_alloc_EB.min());
Efield_fp[lev][2]->plus(m_p_ext_field_params->E_external_grid[2], guard_cells.ng_alloc_EB.min());
}
else {
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 (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) {
}
if (m_p_ext_field_params->B_ext_grid_type != ExternalFieldType::default_zero) {
if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::constant) {
Bfield_fp[lev][0]->plus(m_p_ext_field_params->B_external_grid[0], guard_cells.ng_alloc_EB.min());
Bfield_fp[lev][1]->plus(m_p_ext_field_params->B_external_grid[1], guard_cells.ng_alloc_EB.min());
Bfield_fp[lev][2]->plus(m_p_ext_field_params->B_external_grid[2], guard_cells.ng_alloc_EB.min());
}
else {
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);
Expand Down Expand Up @@ -788,7 +794,7 @@ WarpX::PostRestart ()
{
mypc->PostRestart();
for (int lev = 0; lev <= maxLevel(); ++lev) {
LoadExternalFieldsFromFile(lev);
LoadExternalFields(lev);
}
}

Expand All @@ -810,7 +816,6 @@ WarpX::InitLevelData (int lev, Real /*time*/)
m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::default_zero;
if ( is_B_ext_const && (lev <= maxlevel_extEMfield_init) )
{
Bfield_fp[lev][i]->setVal(m_p_ext_field_params->B_external_grid[i]);
if (fft_do_time_averaging) {
Bfield_avg_fp[lev][i]->setVal(m_p_ext_field_params->B_external_grid[i]);
}
Expand All @@ -831,7 +836,6 @@ WarpX::InitLevelData (int lev, Real /*time*/)
m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::default_zero;
if ( is_E_ext_const && (lev <= maxlevel_extEMfield_init) )
{
Efield_fp[lev][i]->setVal(m_p_ext_field_params->E_external_grid[i]);
if (fft_do_time_averaging) {
Efield_avg_fp[lev][i]->setVal(m_p_ext_field_params->E_external_grid[i]);
}
Expand All @@ -856,13 +860,12 @@ WarpX::InitLevelData (int lev, Real /*time*/)
// Externally imposed fields are only initialized until the user-defined maxlevel_extEMfield_init.
// The default maxlevel_extEMfield_init value is the total number of levels in the simulation
if ((m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::parse_ext_grid_function)
&& (lev <= maxlevel_extEMfield_init)) {
&& (lev > 0) && (lev <= maxlevel_extEMfield_init)) {

// Initialize Bfield_fp with external function
InitializeExternalFieldsOnGridUsingParser(
Bfield_fp[lev][0].get(),
Bfield_fp[lev][1].get(),
Bfield_fp[lev][2].get(),
Bfield_aux[lev][0].get(),
Bfield_aux[lev][1].get(),
Bfield_aux[lev][2].get(),
m_p_ext_field_params->Bxfield_parser->compile<3>(),
m_p_ext_field_params->Byfield_parser->compile<3>(),
m_p_ext_field_params->Bzfield_parser->compile<3>(),
Expand All @@ -871,31 +874,17 @@ WarpX::InitLevelData (int lev, Real /*time*/)
'B',
lev, PatchType::fine);

if (lev > 0) {
InitializeExternalFieldsOnGridUsingParser(
Bfield_aux[lev][0].get(),
Bfield_aux[lev][1].get(),
Bfield_aux[lev][2].get(),
m_p_ext_field_params->Bxfield_parser->compile<3>(),
m_p_ext_field_params->Byfield_parser->compile<3>(),
m_p_ext_field_params->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(),
m_p_ext_field_params->Bxfield_parser->compile<3>(),
m_p_ext_field_params->Byfield_parser->compile<3>(),
m_p_ext_field_params->Bzfield_parser->compile<3>(),
m_edge_lengths[lev],
m_face_areas[lev],
'B',
lev, PatchType::coarse);
}
InitializeExternalFieldsOnGridUsingParser(
Bfield_cp[lev][0].get(),
Bfield_cp[lev][1].get(),
Bfield_cp[lev][2].get(),
m_p_ext_field_params->Bxfield_parser->compile<3>(),
m_p_ext_field_params->Byfield_parser->compile<3>(),
m_p_ext_field_params->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",
Expand All @@ -906,19 +895,6 @@ WarpX::InitLevelData (int lev, Real /*time*/)
if ((m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::parse_ext_grid_function)
&& (lev <= maxlevel_extEMfield_init)) {

// Initialize Efield_fp with external function
InitializeExternalFieldsOnGridUsingParser(
Efield_fp[lev][0].get(),
Efield_fp[lev][1].get(),
Efield_fp[lev][2].get(),
m_p_ext_field_params->Exfield_parser->compile<3>(),
m_p_ext_field_params->Eyfield_parser->compile<3>(),
m_p_ext_field_params->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) {
Expand Down Expand Up @@ -963,7 +939,10 @@ WarpX::InitLevelData (int lev, Real /*time*/)
}
}

LoadExternalFieldsFromFile(lev);
// load external grid fields into E/Bfield_fp_external multifabs
LoadExternalFields(lev);
// add the external fields to the fine patch fields as initial conditions for the fields
AddExternalFields(lev);

if (costs[lev]) {
const auto iarr = costs[lev]->IndexArray();
Expand Down Expand Up @@ -1354,7 +1333,7 @@ void WarpX::CheckKnownIssues()
}

void
WarpX::LoadExternalFieldsFromFile (int const lev)
WarpX::LoadExternalFields (int const lev)
{
// External fields from file are currently not compatible with the moving window
// In order to support the moving window, the MultiFab containing the external
Expand All @@ -1370,8 +1349,21 @@ WarpX::LoadExternalFieldsFromFile (int const lev)
}

// External grid fields

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::parse_ext_grid_function) {
// Initialize Bfield_fp_external with external function
InitializeExternalFieldsOnGridUsingParser(
Bfield_fp_external[lev][0].get(),
Bfield_fp_external[lev][1].get(),
Bfield_fp_external[lev][2].get(),
m_p_ext_field_params->Bxfield_parser->compile<3>(),
m_p_ext_field_params->Byfield_parser->compile<3>(),
m_p_ext_field_params->Bzfield_parser->compile<3>(),
m_edge_lengths[lev],
m_face_areas[lev],
'B',
lev, PatchType::fine);
}
else if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::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 @@ -1384,7 +1376,22 @@ 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::parse_ext_grid_function) {
// Initialize Efield_fp_external with external function
InitializeExternalFieldsOnGridUsingParser(
Efield_fp_external[lev][0].get(),
Efield_fp_external[lev][1].get(),
Efield_fp_external[lev][2].get(),
m_p_ext_field_params->Exfield_parser->compile<3>(),
m_p_ext_field_params->Eyfield_parser->compile<3>(),
m_p_ext_field_params->Ezfield_parser->compile<3>(),
m_edge_lengths[lev],
m_face_areas[lev],
'E',
lev, PatchType::fine);
}
else if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::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
4 changes: 2 additions & 2 deletions Source/WarpX.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1050,7 +1050,7 @@ public:
* \brief Load field values from a user-specified openPMD file,
* for the fields Ex, Ey, Ez, Bx, By, Bz
*/
void LoadExternalFieldsFromFile (int lev);
void LoadExternalFields (int lev);

/**
* \brief Load field values from a user-specified openPMD file
Expand Down Expand Up @@ -1293,7 +1293,7 @@ private:
void FillBoundaryB_avg (int lev, PatchType patch_type, amrex::IntVect ng);
void FillBoundaryE_avg (int lev, PatchType patch_type, amrex::IntVect ng);

void AddExternalFields ();
void AddExternalFields (int lev);

void OneStep_nosub (amrex::Real cur_time);
void OneStep_sub1 (amrex::Real cur_time);
Expand Down
10 changes: 5 additions & 5 deletions Source/WarpX.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -285,7 +285,7 @@ WarpX::WarpX ()

// Loop over species (particles and lasers)
// and set current injection position per species
if (do_moving_window){
if (do_moving_window){
const int n_containers = mypc->nContainers();
for (int i=0; i<n_containers; i++)
{
Expand Down Expand Up @@ -341,8 +341,8 @@ WarpX::WarpX ()
}

// Same as Bfield_fp/Efield_fp for reading external field data
Bfield_fp_external.resize(1);
Efield_fp_external.resize(1);
Bfield_fp_external.resize(nlevs_max);
Efield_fp_external.resize(nlevs_max);
B_external_particle_field.resize(1);
E_external_particle_field.resize(1);

Expand Down Expand Up @@ -2596,7 +2596,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
}

// The external fields that are read from file
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::default_zero && m_p_ext_field_params->B_ext_grid_type != ExternalFieldType::constant) {
// These fields will be added directly to the grid, i.e. to fp, and need to match the index type
AllocInitMultiFab(Bfield_fp_external[lev][0], amrex::convert(ba, Bfield_fp[lev][0]->ixType()),
dm, ncomps, ngEB, lev, "Bfield_fp_external[x]", 0.0_rt);
Expand All @@ -2614,7 +2614,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
AllocInitMultiFab(B_external_particle_field[lev][2], amrex::convert(ba, Bfield_aux[lev][2]->ixType()),
dm, ncomps, ngEB, lev, "B_external_particle_field[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::default_zero && m_p_ext_field_params->E_ext_grid_type != ExternalFieldType::constant) {
// These fields will be added directly to the grid, i.e. to fp, and need to match the index type
AllocInitMultiFab(Efield_fp_external[lev][0], amrex::convert(ba, Efield_fp[lev][0]->ixType()),
dm, ncomps, ngEB, lev, "Efield_fp_external[x]", 0.0_rt);
Expand Down

0 comments on commit f1f3c38

Please sign in to comment.