aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/Evolve/WarpXEvolve.cpp1
-rw-r--r--Source/Initialization/WarpXInitData.cpp252
-rw-r--r--Source/Parallelization/WarpXRegrid.cpp6
-rw-r--r--Source/WarpX.H15
-rw-r--r--Source/WarpX.cpp34
5 files changed, 308 insertions, 0 deletions
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
index 21deb66d1..269270bd9 100644
--- a/Source/Evolve/WarpXEvolve.cpp
+++ b/Source/Evolve/WarpXEvolve.cpp
@@ -324,6 +324,7 @@ WarpX::Evolve (int numsteps)
// This is currently a lab frame calculation.
ComputeMagnetostaticField();
}
+ AddExternalFields();
ExecutePythonCallback("afterEsolve");
}
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index cc2225375..062f8a83b 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -20,6 +20,7 @@
#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"
@@ -72,6 +73,10 @@
#include <string>
#include <utility>
+#ifdef WARPX_USE_OPENPMD
+# include <openPMD/openPMD.hpp>
+#endif
+
#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H"
using namespace amrex;
@@ -435,6 +440,12 @@ 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();
+
// Write full diagnostics before the first iteration.
multi_diags->FilterComputePackFlush( -1 );
@@ -452,6 +463,23 @@ 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 (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();
@@ -870,6 +898,38 @@ WarpX::InitLevelData (int lev, Real /*time*/)
}
}
+ // 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) {
@@ -1242,3 +1302,195 @@ void WarpX::CheckKnownIssues()
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<std::string>() == "C",
+ "Reading from files with non-C dataOrder is not implemented");
+
+ auto axisLabels = F.getAttribute("axisLabels").get<std::vector<std::string>>();
+ auto fileGeom = F.getAttribute("geometry").get<std::string>();
+
+#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)
+ amrex::Abort(Utils::TextMsg::Err(
+ "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<long double>();
+
+#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<double>(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<double> 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<double> 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<double>
+ (xx0, xx0+file_dr, xx1, xx1+file_dz,
+ f00, f01, f10, f11,
+ x0, x1);
+#elif defined(WARPX_DIM_3D)
+ amrex::Array4<double> 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<double>
+ (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)
+ Abort(Utils::TextMsg::Err("Reading fields from openPMD files is not supported in 1D");
+#elif defined(WARPX_DIM_XZ)
+ Abort(Utils::TextMsg::Err(
+ "Reading from openPMD for external fields is not known to work with XZ (see #3828)"));
+#elif !defined(WARPX_USE_OPENPMD)
+ Abort(Utils::TextMsg::Err("OpenPMD field reading requires OpenPMD support to be enabled"));
+#endif
+}
+#endif // WARPX_USE_OPENPMD
diff --git a/Source/Parallelization/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp
index 3c90a07cb..16ecd49ce 100644
--- a/Source/Parallelization/WarpXRegrid.cpp
+++ b/Source/Parallelization/WarpXRegrid.cpp
@@ -169,6 +169,12 @@ WarpX::RemakeLevel (int lev, Real /*time*/, const BoxArray& ba, const Distributi
{
RemakeMultiFab(Bfield_fp[lev][idim], dm, true);
RemakeMultiFab(Efield_fp[lev][idim], dm, true);
+ if (add_external_B_field) {
+ RemakeMultiFab(Bfield_fp_external[lev][idim], dm, true);
+ }
+ if (add_external_E_field) {
+ RemakeMultiFab(Efield_fp_external[lev][idim], dm, true);
+ }
RemakeMultiFab(current_fp[lev][idim], dm, false);
RemakeMultiFab(current_store[lev][idim], dm, false);
if (current_deposition_algo == CurrentDepositionAlgo::Vay) {
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 7d06b8b19..7b013f1c6 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -122,6 +122,11 @@ public:
//! Initialization type for external electric field on the grid
static std::string E_ext_grid_s;
+ //! Whether to apply the effect of an externally-defined electric field
+ static bool add_external_E_field;
+ //! Whether to apply the effect of an externally-defined magnetic field
+ static bool add_external_B_field;
+
//! String storing parser function to initialize x-component of the magnetic field on the grid
static std::string str_Bx_ext_grid_function;
//! String storing parser function to initialize y-component of the magnetic field on the grid
@@ -976,6 +981,10 @@ public:
const char field,
const int lev, PatchType patch_type);
+ void ReadExternalFieldFromFile (
+ std::string read_fields_from_path, amrex::MultiFab* mf,
+ std::string F_name, std::string F_component);
+
/**
* \brief
* This function initializes and calculates grid quantities used along with
@@ -1099,6 +1108,8 @@ 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 OneStep_nosub (amrex::Real t);
void OneStep_sub1 (amrex::Real t);
@@ -1311,6 +1322,10 @@ private:
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > vector_potential_grad_buf_e_stag;
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > vector_potential_grad_buf_b_stag;
+ // Same as Bfield_fp/Efield_fp for reading external field data
+ amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_fp_external;
+ amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_fp_external;
+
//! EB: Lengths of the mesh edges
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > m_edge_lengths;
//! EB: Areas of the mesh faces
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 518d86def..dfc340861 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -88,6 +88,8 @@ Vector<Real> WarpX::B_external_grid(3, 0.0);
std::string WarpX::authors = "";
std::string WarpX::B_ext_grid_s = "default";
std::string WarpX::E_ext_grid_s = "default";
+bool WarpX::add_external_E_field = false;
+bool WarpX::add_external_B_field = false;
// Parser for B_external on the grid
std::string WarpX::str_Bx_ext_grid_function;
@@ -310,6 +312,10 @@ WarpX::WarpX ()
Bfield_avg_fp.resize(nlevs_max);
}
+ // Same as Bfield_fp/Efield_fp for reading external field data
+ Bfield_fp_external.resize(1);
+ Efield_fp_external.resize(1);
+
m_edge_lengths.resize(nlevs_max);
m_face_areas.resize(nlevs_max);
m_distance_to_eb.resize(nlevs_max);
@@ -671,6 +677,22 @@ WarpX::ReadParameters ()
moving_window_v *= PhysConst::c;
}
+ pp_warpx.query("B_ext_grid_init_style", WarpX::B_ext_grid_s);
+ pp_warpx.query("E_ext_grid_init_style", WarpX::E_ext_grid_s);
+
+ if (WarpX::B_ext_grid_s == "read_from_file")
+ {
+ WARPX_ALWAYS_ASSERT_WITH_MESSAGE(max_level == 0,
+ "External field reading is not implemented for more than one level");
+ add_external_B_field = true;
+ }
+ if (WarpX::E_ext_grid_s == "read_from_file")
+ {
+ WARPX_ALWAYS_ASSERT_WITH_MESSAGE(max_level == 0,
+ "External field reading is not implemented for more than one level");
+ add_external_E_field = true;
+ }
+
electrostatic_solver_id = GetAlgorithmInteger(pp_warpx, "do_electrostatic");
// if an electrostatic solver is used, set the Maxwell solver to None
if (electrostatic_solver_id != ElectrostaticSolverAlgo::None) {
@@ -2074,6 +2096,18 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
AllocInitMultiFab(current_fp[lev][1], amrex::convert(ba, jy_nodal_flag), dm, ncomps, ngJ, tag("current_fp[y]"), 0.0_rt);
AllocInitMultiFab(current_fp[lev][2], amrex::convert(ba, jz_nodal_flag), dm, ncomps, ngJ, tag("current_fp[z]"), 0.0_rt);
+ // Match external field MultiFabs to fine patch
+ if (add_external_B_field) {
+ AllocInitMultiFab(Bfield_fp_external[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, tag("Bfield_fp_external[x]"));
+ AllocInitMultiFab(Bfield_fp_external[lev][1], amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, tag("Bfield_fp_external[y]"));
+ AllocInitMultiFab(Bfield_fp_external[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, tag("Bfield_fp_external[z]"));
+ }
+ if (add_external_E_field) {
+ AllocInitMultiFab(Efield_fp_external[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, tag("Efield_fp_external[x]"));
+ AllocInitMultiFab(Efield_fp_external[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, tag("Efield_fp_external[y]"));
+ AllocInitMultiFab(Efield_fp_external[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, tag("Efield_fp_external[z]"));
+ }
+
if (do_current_centering)
{
amrex::BoxArray const& nodal_ba = amrex::convert(ba, amrex::IntVect::TheNodeVector());
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-17 23:01:21 +0000