aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/Initialization/PlasmaInjector.cpp56
-rw-r--r--Source/Initialization/WarpXInitData.cpp258
-rw-r--r--Source/Laser/LaserParticleContainer.cpp3
-rw-r--r--Source/Laser/LaserProfiles.H195
-rw-r--r--Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp2
-rw-r--r--Source/Laser/LaserProfilesImpl/LaserProfileFromTXYEFile.cpp479
-rw-r--r--Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp2
-rw-r--r--Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp2
-rw-r--r--Source/Laser/LaserProfilesImpl/Make.package1
-rw-r--r--Source/Parser/Make.package1
-rw-r--r--Source/Parser/WarpXParserWrapper.H35
-rw-r--r--Source/Utils/WarpXMovingWindow.cpp77
-rw-r--r--Source/Utils/WarpXUtil.H110
-rw-r--r--Source/Utils/WarpXUtil.cpp41
-rw-r--r--Source/WarpX.H74
-rw-r--r--Source/WarpX.cpp15
16 files changed, 1257 insertions, 94 deletions
diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp
index f7c7e498f..5f75ed45a 100644
--- a/Source/Initialization/PlasmaInjector.cpp
+++ b/Source/Initialization/PlasmaInjector.cpp
@@ -3,6 +3,7 @@
#include <WarpXConst.H>
#include <WarpX_f.H>
#include <WarpX.H>
+#include <WarpXUtil.H>
#include <AMReX.H>
@@ -168,34 +169,6 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
}
}
-namespace {
-WarpXParser makeParser (std::string const& parse_function)
-{
- WarpXParser parser(parse_function);
- parser.registerVariables({"x","y","z"});
-
- ParmParse pp("my_constants");
- std::set<std::string> symbols = parser.symbols();
- symbols.erase("x");
- symbols.erase("y");
- symbols.erase("z"); // after removing variables, we are left with constants
- for (auto it = symbols.begin(); it != symbols.end(); ) {
- Real v;
- if (pp.query(it->c_str(), v)) {
- parser.setConstant(*it, v);
- it = symbols.erase(it);
- } else {
- ++it;
- }
- }
- for (auto const& s : symbols) { // make sure there no unknown symbols
- amrex::Abort("PlasmaInjector::makeParser: Unknown symbol "+s);
- }
-
- return parser;
-}
-}
-
// Depending on injection type at runtime, initialize inj_rho
// so that inj_rho->getDensity calls
// InjectorPosition[Constant or Custom or etc.].getDensity.
@@ -217,11 +190,7 @@ void PlasmaInjector::parseDensity (ParmParse& pp)
// Construct InjectorDensity with InjectorDensityPredefined.
inj_rho.reset(new InjectorDensity((InjectorDensityPredefined*)nullptr,species_name));
} else if (rho_prof_s == "parse_density_function") {
- std::vector<std::string> f;
- pp.getarr("density_function(x,y,z)", f);
- for (auto const& s : f) {
- str_density_function += s;
- }
+ Store_parserString(pp, "density_function(x,y,z)", str_density_function);
// Construct InjectorDensity with InjectorDensityParser.
inj_rho.reset(new InjectorDensity((InjectorDensityParser*)nullptr,
makeParser(str_density_function)));
@@ -339,21 +308,12 @@ void PlasmaInjector::parseMomentum (ParmParse& pp)
inj_mom.reset(new InjectorMomentum
((InjectorMomentumRadialExpansion*)nullptr, u_over_r));
} else if (mom_dist_s == "parse_momentum_function") {
- std::vector<std::string> f;
- pp.getarr("momentum_function_ux(x,y,z)", f);
- for (auto const& s : f) {
- str_momentum_function_ux += s;
- }
- f.clear();
- pp.getarr("momentum_function_uy(x,y,z)", f);
- for (auto const& s : f) {
- str_momentum_function_uy += s;
- }
- f.clear();
- pp.getarr("momentum_function_uz(x,y,z)", f);
- for (auto const& s : f) {
- str_momentum_function_uz += s;
- }
+ Store_parserString(pp, "momentum_function_ux(x,y,z)",
+ str_momentum_function_ux);
+ Store_parserString(pp, "momentum_function_uy(x,y,z)",
+ str_momentum_function_uy);
+ Store_parserString(pp, "momentum_function_uz(x,y,z)",
+ str_momentum_function_uz);
// Construct InjectorMomentum with InjectorMomentumParser.
inj_mom.reset(new InjectorMomentum((InjectorMomentumParser*)nullptr,
makeParser(str_momentum_function_ux),
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index 29c9a8955..be29a1cbc 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -1,4 +1,3 @@
-
#include <WarpX.H>
#include <WarpX_f.H>
#include <BilinearFilter.H>
@@ -10,6 +9,9 @@
#ifdef BL_USE_SENSEI_INSITU
#include <AMReX_AmrMeshInSituBridge.H>
#endif
+#include <GpuParser.H>
+#include <WarpXUtil.H>
+
using namespace amrex;
@@ -231,24 +233,161 @@ WarpX::PostRestart ()
mypc->PostRestart();
}
+
void
WarpX::InitLevelData (int lev, Real time)
{
+
+ ParmParse pp("warpx");
+
+ // default values of E_external_grid and B_external_grid
+ // are used to set the E and B field when "constant" or
+ // "parser" is not explicitly used in the input.
+ pp.query("B_ext_grid_init_style", B_ext_grid_s);
+ std::transform(B_ext_grid_s.begin(),
+ B_ext_grid_s.end(),
+ B_ext_grid_s.begin(),
+ ::tolower);
+
+ pp.query("E_ext_grid_init_style", E_ext_grid_s);
+ std::transform(E_ext_grid_s.begin(),
+ E_ext_grid_s.end(),
+ E_ext_grid_s.begin(),
+ ::tolower);
+
+ // if the input string is "constant", the values for the
+ // external grid must be provided in the input.
+ if (B_ext_grid_s == "constant")
+ pp.getarr("B_external_grid", B_external_grid);
+
+ // if the input string is "constant", the values for the
+ // external grid must be provided in the input.
+ if (E_ext_grid_s == "constant")
+ pp.getarr("E_external_grid", E_external_grid);
+
for (int i = 0; i < 3; ++i) {
current_fp[lev][i]->setVal(0.0);
- Efield_fp[lev][i]->setVal(E_external_grid[i]);
- Bfield_fp[lev][i]->setVal(B_external_grid[i]);
+ if (lev > 0)
+ current_cp[lev][i]->setVal(0.0);
+
+ if (B_ext_grid_s == "constant" || B_ext_grid_s == "default") {
+ Bfield_fp[lev][i]->setVal(B_external_grid[i]);
+ if (lev > 0) {
+ Bfield_aux[lev][i]->setVal(B_external_grid[i]);
+ Bfield_cp[lev][i]->setVal(B_external_grid[i]);
+ }
+ }
+ if (E_ext_grid_s == "constant" || E_ext_grid_s == "default") {
+ Efield_fp[lev][i]->setVal(E_external_grid[i]);
+ if (lev > 0) {
+ Efield_aux[lev][i]->setVal(E_external_grid[i]);
+ Efield_cp[lev][i]->setVal(E_external_grid[i]);
+ }
+ }
}
- if (lev > 0) {
- for (int i = 0; i < 3; ++i) {
- Efield_aux[lev][i]->setVal(E_external_grid[i]);
- Bfield_aux[lev][i]->setVal(B_external_grid[i]);
+ // if the input string for the B-field is "parse_b_ext_grid_function",
+ // then the analytical expression or function must be
+ // provided in the input file.
+ if (B_ext_grid_s == "parse_b_ext_grid_function") {
- current_cp[lev][i]->setVal(0.0);
- Efield_cp[lev][i]->setVal(E_external_grid[i]);
- Bfield_cp[lev][i]->setVal(B_external_grid[i]);
- }
+#ifdef WARPX_DIM_RZ
+ amrex::Abort("E and B parser for external fields does not work with RZ -- TO DO");
+#endif
+ Store_parserString(pp, "Bx_external_grid_function(x,y,z)",
+ str_Bx_ext_grid_function);
+ Store_parserString(pp, "By_external_grid_function(x,y,z)",
+ str_By_ext_grid_function);
+ Store_parserString(pp, "Bz_external_grid_function(x,y,z)",
+ str_Bz_ext_grid_function);
+
+ Bxfield_parser.reset(new ParserWrapper(
+ makeParser(str_Bx_ext_grid_function)));
+ Byfield_parser.reset(new ParserWrapper(
+ makeParser(str_By_ext_grid_function)));
+ Bzfield_parser.reset(new ParserWrapper(
+ makeParser(str_Bz_ext_grid_function)));
+
+ // Initialize Bfield_fp with external function
+ InitializeExternalFieldsOnGridUsingParser(Bfield_fp[lev][0].get(),
+ Bfield_fp[lev][1].get(),
+ Bfield_fp[lev][2].get(),
+ Bxfield_parser.get(),
+ Byfield_parser.get(),
+ Bzfield_parser.get(),
+ Bx_nodal_flag, By_nodal_flag,
+ Bz_nodal_flag, lev);
+ if (lev > 0) {
+ InitializeExternalFieldsOnGridUsingParser(Bfield_aux[lev][0].get(),
+ Bfield_aux[lev][1].get(),
+ Bfield_aux[lev][2].get(),
+ Bxfield_parser.get(),
+ Byfield_parser.get(),
+ Bzfield_parser.get(),
+ Bx_nodal_flag, By_nodal_flag,
+ Bz_nodal_flag, lev);
+
+ InitializeExternalFieldsOnGridUsingParser(Bfield_cp[lev][0].get(),
+ Bfield_cp[lev][1].get(),
+ Bfield_cp[lev][2].get(),
+ Bxfield_parser.get(),
+ Byfield_parser.get(),
+ Bzfield_parser.get(),
+ Bx_nodal_flag, By_nodal_flag,
+ Bz_nodal_flag, lev);
+ }
+ }
+
+ // if the input string for the E-field is "parse_e_ext_grid_function",
+ // then the analytical expression or function must be
+ // provided in the input file.
+ if (E_ext_grid_s == "parse_e_ext_grid_function") {
+
+#ifdef WARPX_DIM_RZ
+ amrex::Abort("E and B parser for external fields does not work with RZ -- TO DO");
+#endif
+ Store_parserString(pp, "Ex_external_grid_function(x,y,z)",
+ str_Ex_ext_grid_function);
+ Store_parserString(pp, "Ey_external_grid_function(x,y,z)",
+ str_Ey_ext_grid_function);
+ Store_parserString(pp, "Ez_external_grid_function(x,y,z)",
+ str_Ez_ext_grid_function);
+
+ Exfield_parser.reset(new ParserWrapper(
+ makeParser(str_Ex_ext_grid_function)));
+ Eyfield_parser.reset(new ParserWrapper(
+ makeParser(str_Ey_ext_grid_function)));
+ Ezfield_parser.reset(new ParserWrapper(
+ makeParser(str_Ez_ext_grid_function)));
+
+ // Initialize Efield_fp with external function
+ InitializeExternalFieldsOnGridUsingParser(Efield_fp[lev][0].get(),
+ Efield_fp[lev][1].get(),
+ Efield_fp[lev][2].get(),
+ Exfield_parser.get(),
+ Eyfield_parser.get(),
+ Ezfield_parser.get(),
+ Ex_nodal_flag, Ey_nodal_flag,
+ Ez_nodal_flag, lev);
+ if (lev > 0) {
+ InitializeExternalFieldsOnGridUsingParser(Efield_aux[lev][0].get(),
+ Efield_aux[lev][1].get(),
+ Efield_aux[lev][2].get(),
+ Exfield_parser.get(),
+ Eyfield_parser.get(),
+ Ezfield_parser.get(),
+ Ex_nodal_flag, Ey_nodal_flag,
+ Ez_nodal_flag, lev);
+
+ InitializeExternalFieldsOnGridUsingParser(Efield_cp[lev][0].get(),
+ Efield_cp[lev][1].get(),
+ Efield_cp[lev][2].get(),
+ Exfield_parser.get(),
+ Eyfield_parser.get(),
+ Ezfield_parser.get(),
+ Ex_nodal_flag, Ey_nodal_flag,
+ Ez_nodal_flag, lev);
+ }
}
if (F_fp[lev]) {
@@ -306,3 +445,100 @@ WarpX::InitLevelDataFFT (int lev, Real time)
}
#endif
+
+void
+WarpX::InitializeExternalFieldsOnGridUsingParser (
+ MultiFab *mfx, MultiFab *mfy, MultiFab *mfz,
+ ParserWrapper *xfield_parser, ParserWrapper *yfield_parser,
+ ParserWrapper *zfield_parser, IntVect x_nodal_flag,
+ IntVect y_nodal_flag, IntVect z_nodal_flag,
+ const int lev)
+{
+
+ const auto dx_lev = geom[lev].CellSizeArray();
+ const RealBox& real_box = geom[lev].ProbDomain();
+ for ( MFIter mfi(*mfx, TilingIfNotGPU()); mfi.isValid(); ++mfi)
+ {
+ const Box& tbx = mfi.tilebox(x_nodal_flag);
+ const Box& tby = mfi.tilebox(y_nodal_flag);
+ const Box& tbz = mfi.tilebox(z_nodal_flag);
+
+ auto const& mfxfab = mfx->array(mfi);
+ auto const& mfyfab = mfy->array(mfi);
+ auto const& mfzfab = mfz->array(mfi);
+
+ auto const& mfx_IndexType = (*mfx).ixType();
+ auto const& mfy_IndexType = (*mfy).ixType();
+ auto const& mfz_IndexType = (*mfz).ixType();
+
+ // Initialize IntVect based on the index type of multiFab
+ // 0 if cell-centered, 1 if node-centered.
+ IntVect mfx_type(AMREX_D_DECL(0,0,0));
+ IntVect mfy_type(AMREX_D_DECL(0,0,0));
+ IntVect mfz_type(AMREX_D_DECL(0,0,0));
+
+ for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+ mfx_type[idim] = mfx_IndexType.nodeCentered(idim);
+ mfy_type[idim] = mfy_IndexType.nodeCentered(idim);
+ mfz_type[idim] = mfz_IndexType.nodeCentered(idim);
+ }
+
+ amrex::ParallelFor (tbx, tby, tbz,
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+ // Shift required in the x-, y-, or z- position
+ // depending on the index type of the multifab
+ Real fac_x = (1.0 - mfx_type[0]) * dx_lev[0]*0.5;
+ Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
+#if (AMREX_SPACEDIM==2)
+ Real y = 0.0;
+ Real fac_z = (1.0 - mfx_type[1]) * dx_lev[1]*0.5;
+ Real z = j*dx_lev[1] + real_box.lo(1) + fac_z;
+#else
+ Real fac_y = (1.0 - mfx_type[1]) * dx_lev[1]*0.5;
+ Real y = j*dx_lev[1] + real_box.lo(1) + fac_y;
+ Real fac_z = (1.0 - mfx_type[2]) * dx_lev[2]*0.5;
+ Real z = k*dx_lev[2] + real_box.lo(2) + fac_z;
+#endif
+ // Initialize the x-component of the field.
+ mfxfab(i,j,k) = xfield_parser->getField(x,y,z);
+ },
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+ Real fac_x = (1.0 - mfy_type[0]) * dx_lev[0]*0.5;
+ Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
+#if (AMREX_SPACEDIM==2)
+ Real y = 0.0;
+ Real fac_z = (1.0 - mfx_type[1]) * dx_lev[1]*0.5;
+ Real z = j*dx_lev[1] + real_box.lo(1) + fac_z;
+#elif (AMREX_SPACEDIM==3)
+ Real fac_y = (1.0 - mfx_type[1]) * dx_lev[1]*0.5;
+ Real y = j*dx_lev[1] + real_box.lo(1) + fac_y;
+ Real fac_z = (1.0 - mfx_type[2]) * dx_lev[2]*0.5;
+ Real z = k*dx_lev[2] + real_box.lo(2) + fac_z;
+#endif
+ // Initialize the y-component of the field.
+ mfyfab(i,j,k) = yfield_parser->getField(x,y,z);
+ },
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+ Real fac_x = (1.0 - mfz_type[0]) * dx_lev[0]*0.5;
+ Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
+#if (AMREX_SPACEDIM==2)
+ Real y = 0.0;
+ Real fac_z = (1.0 - mfx_type[1]) * dx_lev[1]*0.5;
+ Real z = j*dx_lev[1] + real_box.lo(1) + fac_z;
+#elif (AMREX_SPACEDIM==3)
+ Real fac_y = (1.0 - mfx_type[1]) * dx_lev[1]*0.5;
+ Real y = j*dx_lev[1] + real_box.lo(1) + fac_y;
+ Real fac_z = (1.0 - mfz_type[2]) * dx_lev[2]*0.5;
+ Real z = k*dx_lev[2] + real_box.lo(2) + fac_z;
+#endif
+ // Initialize the z-component of the field.
+ mfzfab(i,j,k) = zfield_parser->getField(x,y,z);
+ },
+ /* To allocate shared memory for the GPU threads. */
+ /* But, for now only 3 doubles (x,y,z) are allocated. */
+ amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 3
+ );
+
+ }
+
+}
diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp
index 1d9689979..ed9f5eda0 100644
--- a/Source/Laser/LaserParticleContainer.cpp
+++ b/Source/Laser/LaserParticleContainer.cpp
@@ -386,6 +386,9 @@ LaserParticleContainer::Evolve (int lev,
t_lab = 1./WarpX::gamma_boost*t + WarpX::beta_boost*Z0_lab/PhysConst::c;
}
+ // Update laser profile
+ m_up_laser_profile->update(t);
+
BL_ASSERT(OnSameGrids(lev,jx));
MultiFab* cost = WarpX::getCosts(lev);
diff --git a/Source/Laser/LaserProfiles.H b/Source/Laser/LaserProfiles.H
index e0ec0dc28..f97bf915e 100644
--- a/Source/Laser/LaserProfiles.H
+++ b/Source/Laser/LaserProfiles.H
@@ -5,6 +5,7 @@
#include <WarpXParser.H>
#include <AMReX_ParmParse.H>
#include <AMReX_Vector.H>
+#include <AMReX_Gpu.H>
#include <WarpXParser.H>
@@ -13,6 +14,7 @@
#include <memory>
#include <functional>
#include <limits>
+#include <utility>
namespace WarpXLaserProfiles {
@@ -31,8 +33,8 @@ struct CommonLaserParameters
/** Abstract interface for laser profile classes
*
- * Each new laser profile should inherit this interface and implement two
- * methods: init and fill_amplitude (described below).
+ * Each new laser profile should inherit this interface and implement three
+ * methods: init, update and fill_amplitude (described below).
*
* The declaration of a LaserProfile class should be placed in this file,
* while the implementation of the methods should be in a dedicated file in
@@ -60,6 +62,17 @@ public:
const amrex::ParmParse& ppc,
CommonLaserParameters params) = 0;
+ /** Update Laser Profile
+ *
+ * Some laser profiles might need to perform an "update" operation per
+ * time step.
+ *
+ * @param[in] t Current physical time in the simulation (seconds)
+ */
+ virtual void
+ update (
+ amrex::Real t) = 0;
+
/** Fill Electric Field Amplitude for each particle of the antenna.
*
* Xp, Yp and amplitude must be arrays with the same length
@@ -76,7 +89,7 @@ public:
amrex::Real const * AMREX_RESTRICT const Xp,
amrex::Real const * AMREX_RESTRICT const Yp,
amrex::Real t,
- amrex::Real * AMREX_RESTRICT const amplitude) = 0;
+ amrex::Real * AMREX_RESTRICT const amplitude) const = 0;
virtual ~ILaserProfile(){};
};
@@ -94,13 +107,17 @@ public:
const amrex::ParmParse& ppc,
CommonLaserParameters params) override final;
+ //No update needed
+ void
+ update (amrex::Real /*t */) override final {}
+
void
fill_amplitude (
const int np,
amrex::Real const * AMREX_RESTRICT const Xp,
amrex::Real const * AMREX_RESTRICT const Yp,
amrex::Real t,
- amrex::Real * AMREX_RESTRICT const amplitude) override final;
+ amrex::Real * AMREX_RESTRICT const amplitude) const override final;
private:
struct {
@@ -132,13 +149,17 @@ public:
const amrex::ParmParse& ppc,
CommonLaserParameters params) override final;
+ //No update needed
+ void
+ update (amrex::Real /*t */) override final {}
+
void
fill_amplitude (
const int np,
amrex::Real const * AMREX_RESTRICT const Xp,
amrex::Real const * AMREX_RESTRICT const Yp,
amrex::Real t,
- amrex::Real * AMREX_RESTRICT const amplitude) override final;
+ amrex::Real * AMREX_RESTRICT const amplitude) const override final;
private:
struct {
@@ -163,13 +184,17 @@ public:
const amrex::ParmParse& ppc,
CommonLaserParameters params) override final;
+ //No update needed
+ void
+ update (amrex::Real /*t */) override final {}
+
void
fill_amplitude (
const int np,
amrex::Real const * AMREX_RESTRICT const Xp,
amrex::Real const * AMREX_RESTRICT const Yp,
amrex::Real t,
- amrex::Real * AMREX_RESTRICT const amplitude) override final;
+ amrex::Real * AMREX_RESTRICT const amplitude) const override final;
private:
struct{
@@ -180,6 +205,160 @@ private:
};
/**
+ * Laser profile read from an 'TXYE' file
+ * The binary file must contain:
+ * - 3 unsigned integers (4 bytes): nt (points along t), nx (points along x) and ny (points along y)
+ * - nt*nx*ny doubles (8 bytes) in row major order : field amplitude
+ */
+class FromTXYEFileLaserProfile : public ILaserProfile
+{
+
+public:
+ void
+ init (
+ const amrex::ParmParse& ppl,
+ const amrex::ParmParse& ppc,
+ CommonLaserParameters params) override final;
+
+ /** \brief Reads new field data chunk from file if needed
+ *
+ * @param[in] t simulation time (seconds)
+ */
+ void
+ update (amrex::Real t) override final;
+
+ /** \brief compute field amplitude at particles' position for a laser beam
+ * loaded from an E(x,y,t) file.
+ *
+ * Both Xp and Yp are given in laser plane coordinate.
+ * For each particle with position Xp and Yp, this routine computes the
+ * amplitude of the laser electric field, stored in array amplitude.
+ *
+ * \param np: number of laser particles
+ * \param Xp: pointer to first component of positions of laser particles
+ * \param Yp: pointer to second component of positions of laser particles
+ * \param t: Current physical time
+ * \param amplitude: pointer to array of field amplitude.
+ */
+ void
+ fill_amplitude (
+ const int np,
+ amrex::Real const * AMREX_RESTRICT const Xp,
+ amrex::Real const * AMREX_RESTRICT const Yp,
+ amrex::Real t,
+ amrex::Real * AMREX_RESTRICT const amplitude) const override final;
+
+ /** \brief Function to fill the amplitude in case of a uniform grid.
+ * This function cannot be private due to restrictions related to
+ * the use of extended __device__ lambda
+ *
+ * \param idx_t_left index of the last time coordinate < t
+ * \param np: number of laser particles
+ * \param Xp: pointer to first component of positions of laser particles
+ * \param Yp: pointer to second component of positions of laser particles
+ * \param t: Current physical time
+ * \param amplitude: pointer to array of field amplitude.
+ */
+ void internal_fill_amplitude_uniform(
+ const int idx_t_left,
+ const int np,
+ amrex::Real const * AMREX_RESTRICT const Xp,
+ amrex::Real const * AMREX_RESTRICT const Yp,
+ amrex::Real t,
+ amrex::Real * AMREX_RESTRICT const amplitude) const;
+
+ /** \brief Function to fill the amplitude in case of a non-uniform grid.
+ * This function cannot be private due to restrictions related to
+ * the use of extended __device__ lambda
+ *
+ * \param idx_t_left index of the last time coordinate < t
+ * \param np: number of laser particles
+ * \param Xp: pointer to first component of positions of laser particles
+ * \param Yp: pointer to second component of positions of laser particles
+ * \param t: Current physical time
+ * \param amplitude: pointer to array of field amplitude.
+ */
+ void internal_fill_amplitude_nonuniform(
+ const int idx_t_left,
+ const int np,
+ amrex::Real const * AMREX_RESTRICT const Xp,
+ amrex::Real const * AMREX_RESTRICT const Yp,
+ amrex::Real t,
+ amrex::Real * AMREX_RESTRICT const amplitude) const;
+
+private:
+ /** \brief parse a field file in the binary 'txye' format (whose details are given below).
+ *
+ * A 'txye' file should be a binary file with the following format:
+ * -flag to indicate if the grid is uniform or not (1 byte, 0 means non-uniform, !=0 means uniform)
+ * -np, number of timesteps (uint32_t, must be >=2)
+ * -nx, number of points along x (uint32_t, must be >=2)
+ * -ny, number of points along y (uint32_t, must be 1 for 2D simulations and >=2 for 3D simulations)
+ * -timesteps (double[2] if grid is uniform, double[np] otherwise)
+ * -x_coords (double[2] if grid is uniform, double[nx] otherwise)
+ * -y_coords (double[1] if 2D, double[2] if 3D & uniform grid, double[ny] if 3D & non uniform grid)
+ * -field_data (double[nt * nx * ny], with nt being the slowest coordinate).
+ * The spatiotemporal grid must be rectangular.
+ *
+ * \param txye_file_name: name of the file to parse
+ */
+ void parse_txye_file(std::string txye_file_name);
+
+ /** \brief Finds left and right time indices corresponding to time t
+ *
+ *
+ * \param t: simulation time
+ */
+ std::pair<int,int> find_left_right_time_indices(amrex::Real t) const;
+
+ /** \brief Load field data within the temporal range [t_begin, t_end)
+ *
+ * Must be called after having parsed a data file with parse_txye_file.
+ *
+ * \param t_begin: left limit of the timestep range to read
+ * \param t_end: right limit of the timestep range to read (t_end is not read)
+ */
+ void read_data_t_chuck(int t_begin, int t_end);
+
+ /**
+ * \brief m_params contains all the internal parameters
+ * used by this laser profile
+ */
+ struct{
+ /** Name of the file containing the data */
+ std::string txye_file_name;
+ /** Flag to store if the grid is uniform */
+ bool is_grid_uniform = false;
+ /** Dimensions of E_data. nt, nx must be >=2.
+ * If DIM=3, ny must be >=2 as well.
+ * If DIM=2, ny must be 1 */
+ int nt, nx, ny;
+ /** Vector of temporal coordinates. For a non-uniform grid, it contains
+ * all values of time. For a uniform grid, it contains only the start and stop
+ * times and intermediate times are obtained with nt */
+ amrex::Gpu::ManagedVector<amrex::Real> t_coords;
+ /** Vector or x coordinates. For a non-uniform grid, it contains all values
+ * of space dimension x. For a uniform grid, it contains only the min and max
+ * x coordinates, and intermediate positions are obtained with nx */
+ amrex::Gpu::ManagedVector<amrex::Real> x_coords;
+ /** Vector or y coordinates. For a non-uniform grid, it contains all values
+ * of space dimension y. For a uniform grid, it contains only the min and max
+ * y coordinates, and intermediate positions are obtained with ny */
+ amrex::Gpu::ManagedVector<amrex::Real> y_coords;
+ /** Size of the timestep range to load */
+ int time_chunk_size;
+ /** Index of the first timestep in memory */
+ int first_time_index;
+ /** Index of the last timestep in memory */
+ int last_time_index;
+ /** Field data */
+ amrex::Gpu::ManagedVector<amrex::Real> E_data;
+ } m_params;
+
+ CommonLaserParameters m_common_params;
+};
+
+/**
* Maps laser profile names to lambdas returing unique pointers
* to the corresponding laser profile objects.
*/
@@ -195,7 +374,9 @@ laser_profiles_dictionary =
{"harris",
[] () {return std::make_unique<HarrisLaserProfile>();} },
{"parse_field_function",
- [] () {return std::make_unique<FieldFunctionLaserProfile>();} }
+ [] () {return std::make_unique<FieldFunctionLaserProfile>();} },
+ {"from_txye_file",
+ [] () {return std::make_unique<FromTXYEFileLaserProfile>();} }
};
} //WarpXLaserProfiles
diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp
index 3c9d7379a..d34bc6aba 100644
--- a/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp
+++ b/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp
@@ -37,7 +37,7 @@ FieldFunctionLaserProfile::init (
void
FieldFunctionLaserProfile::fill_amplitude (
const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp,
- Real t, Real * AMREX_RESTRICT const amplitude)
+ Real t, Real * AMREX_RESTRICT const amplitude) const
{
for (int i = 0; i < np; ++i) {
amplitude[i] = m_parser.eval(Xp[i], Yp[i], t);
diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileFromTXYEFile.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileFromTXYEFile.cpp
new file mode 100644
index 000000000..8f44c2d5f
--- /dev/null
+++ b/Source/Laser/LaserProfilesImpl/LaserProfileFromTXYEFile.cpp
@@ -0,0 +1,479 @@
+#include <LaserProfiles.H>
+
+#include <WarpX_Complex.H>
+#include <WarpXConst.H>
+#include <WarpXUtil.H>
+
+#include <AMReX_Print.H>
+#include <AMReX_ParallelDescriptor.H>
+
+#include <limits>
+#include <iostream>
+#include <fstream>
+#include <cstdint>
+#include <algorithm>
+
+using namespace amrex;
+using namespace WarpXLaserProfiles;
+
+void
+FromTXYEFileLaserProfile::init (
+ const amrex::ParmParse& ppl,
+ const amrex::ParmParse& /* ppc */,
+ CommonLaserParameters params)
+{
+ if (!std::numeric_limits< double >::is_iec559)
+ {
+ Print() << R"(Warning: double does not comply with IEEE 754: bad
+ things will happen parsing the X, Y and T profiles for the laser!)";
+ }
+
+ // Parse the TXYE file
+ ppl.get("txye_file_name", m_params.txye_file_name);
+ if(m_params.txye_file_name.empty())
+ {
+ Abort("txye_file_name must be provided for txye_file laser profile!");
+ }
+ parse_txye_file(m_params.txye_file_name);
+
+ //Set time_chunk_size
+ m_params.time_chunk_size = m_params.nt;
+ int temp = 1;
+ if(ppl.query("time_chunk_size", temp)){
+ m_params.time_chunk_size = min(
+ temp, m_params.time_chunk_size);
+ }
+ if(m_params.time_chunk_size < 2){
+ Abort("Error! time_chunk_size must be >= 2!");
+ }
+
+ //Allocate memory for E_data Vector
+ const int data_size = m_params.time_chunk_size*
+ m_params.nx*m_params.ny;
+ m_params.E_data = Gpu::ManagedVector<amrex::Real>(data_size);
+
+ //Read first time chunck
+ read_data_t_chuck(0, m_params.time_chunk_size);
+
+ //Copy common params
+ m_common_params = params;
+}
+
+void
+FromTXYEFileLaserProfile::update (amrex::Real t)
+{
+ if(t >= m_params.t_coords.back())
+ return;
+
+ const auto idx_times = find_left_right_time_indices(t);
+ const auto idx_t_left = idx_times.first;
+ const auto idx_t_right = idx_times.second;
+
+ //Load data chunck if needed
+ if(idx_t_right > m_params.last_time_index){
+ read_data_t_chuck(idx_t_left, idx_t_left+m_params.time_chunk_size);
+ }
+}
+
+void
+FromTXYEFileLaserProfile::fill_amplitude (
+ const int np,
+ Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp,
+ Real t, Real * AMREX_RESTRICT const amplitude) const
+{
+ //Amplitude is 0 if time is out of range
+ if(t < m_params.t_coords.front() || t > m_params.t_coords.back()){
+ amrex::ParallelFor(np,
+ [=] AMREX_GPU_DEVICE (int i) {
+ amplitude[i] = 0.0_rt;});
+ return;
+ }
+
+ //Find left and right time indices
+ int idx_t_left, idx_t_right;
+ std::tie(idx_t_left, idx_t_right) = find_left_right_time_indices(t);
+
+ if(idx_t_left < m_params.first_time_index){
+ Abort("Something bad has happened with the simulation time");
+ }
+
+ if(m_params.is_grid_uniform){
+ internal_fill_amplitude_uniform(
+ idx_t_left, np, Xp, Yp, t, amplitude);
+ }
+ else{
+ internal_fill_amplitude_nonuniform(
+ idx_t_left, np, Xp, Yp, t, amplitude);
+ }
+}
+
+void
+FromTXYEFileLaserProfile::parse_txye_file(std::string txye_file_name)
+{
+ if(ParallelDescriptor::IOProcessor()){
+ std::ifstream inp(txye_file_name, std::ios::binary);
+ if(!inp) Abort("Failed to open txye file");
+
+ //Uniform grid flag
+ char flag;
+ inp.read(&flag, 1);
+ if(!inp) Abort("Failed to read sizes from txye file");
+ m_params.is_grid_uniform=flag;
+
+ //Grid points along t, x and y
+ auto const three_uint32_size = sizeof(uint32_t)*3;
+ char buf[three_uint32_size];
+ inp.read(buf, three_uint32_size);
+ if(!inp) Abort("Failed to read sizes from txye file");
+ m_params.nt = reinterpret_cast<uint32_t*>(buf)[0];
+ m_params.nx = reinterpret_cast<uint32_t*>(buf)[1];
+ m_params.ny = reinterpret_cast<uint32_t*>(buf)[2];
+ if(m_params.nt <= 1) Abort("nt in txye file must be >=2");
+ if(m_params.nx <= 1) Abort("nx in txye file must be >=2");
+#if (AMREX_SPACEDIM == 3)
+ if(m_params.ny <= 1) Abort("ny in txye file must be >=2 in 3D");
+#elif(AMREX_SPACEDIM == 2)
+ if(m_params.ny != 1) Abort("ny in txye file must be 1 in 2D");
+#endif
+
+ //Coordinates
+ Vector<double> buf_t, buf_x, buf_y;
+ if(m_params.is_grid_uniform){
+ buf_t.resize(2);
+ buf_x.resize(2);
+#if (AMREX_SPACEDIM == 3)
+ buf_y.resize(2);
+#elif(AMREX_SPACEDIM == 2)
+ buf_y.resize(1);
+#endif
+ }
+ else{
+ buf_t.resize(m_params.nt);
+ buf_x.resize(m_params.nx);
+ buf_y.resize(m_params.ny);
+ }
+ inp.read(reinterpret_cast<char*>(buf_t.dataPtr()),
+ buf_t.size()*sizeof(double));
+ if(!inp)
+ Abort("Failed to read coords from txye file");
+ if (!std::is_sorted(buf_t.begin(), buf_t.end()))
+ Abort("Coordinates are not sorted in txye file");
+ inp.read(reinterpret_cast<char*>(buf_x.dataPtr()),
+ buf_x.size()*sizeof(double));
+ if(!inp)
+ Abort("Failed to read coords from txye file");
+ if (!std::is_sorted(buf_x.begin(), buf_x.end()))
+ Abort("Coordinates are not sorted in txye file");
+ inp.read(reinterpret_cast<char*>(buf_y.dataPtr()),
+ buf_y.size()*sizeof(double));
+ if(!inp)
+ Abort("Failed to read coords from txye file");
+ if (!std::is_sorted(buf_y.begin(), buf_y.end()))
+ Abort("Coordinates are not sorted in txye file");
+ m_params.t_coords = Gpu::ManagedVector<amrex::Real>(buf_t.size());
+ m_params.x_coords = Gpu::ManagedVector<amrex::Real>(buf_x.size());
+ m_params.y_coords = Gpu::ManagedVector<amrex::Real>(buf_y.size());
+ // Convert from double to amrex::Real
+ std::transform(buf_t.begin(), buf_t.end(), m_params.t_coords.begin(),
+ [](auto x) {return static_cast<amrex::Real>(x);} );
+ std::transform(buf_x.begin(), buf_x.end(), m_params.x_coords.begin(),
+ [](auto x) {return static_cast<amrex::Real>(x);} );
+ std::transform(buf_y.begin(), buf_y.end(), m_params.y_coords.begin(),
+ [](auto x) {return static_cast<amrex::Real>(x);} );
+ }
+
+ //Broadcast grid uniformity
+ char is_grid_uniform = m_params.is_grid_uniform;
+ ParallelDescriptor::Bcast(&is_grid_uniform, 1,
+ ParallelDescriptor::IOProcessorNumber());
+ ParallelDescriptor::Barrier();
+ m_params.is_grid_uniform = is_grid_uniform;
+
+ //Broadcast grid size and coordinate sizes
+ //When a non-uniform grid is used, nt, nx and ny are identical
+ //to t_coords.size(), x_coords.size() and y_coords.size().
+ //When a uniform grid is used, nt,nx and ny store the number of points
+ //used for the mesh, while t_coords, x_coords and y_coords store the
+ //extrems in each direaction. Thus t_coords and x_coords in this case
+ //have size 2 and y_coords has size 1 in 2D and size 2 in 3D.
+ int t_sizes[6] = {m_params.nt, m_params.nx, m_params.ny,
+ static_cast<int>(m_params.t_coords.size()),
+ static_cast<int>(m_params.x_coords.size()),
+ static_cast<int>(m_params.y_coords.size())};
+ ParallelDescriptor::Bcast(t_sizes, 6,
+ ParallelDescriptor::IOProcessorNumber());
+ ParallelDescriptor::Barrier();
+ m_params.nt = t_sizes[0]; m_params.nx = t_sizes[1]; m_params.ny = t_sizes[2];
+
+ //Broadcast coordinates
+ if(!ParallelDescriptor::IOProcessor()){
+ m_params.t_coords = Gpu::ManagedVector<amrex::Real>(t_sizes[3]);
+ m_params.x_coords = Gpu::ManagedVector<amrex::Real>(t_sizes[4]);
+ m_params.y_coords = Gpu::ManagedVector<amrex::Real>(t_sizes[5]);
+ }
+ ParallelDescriptor::Bcast(m_params.t_coords.dataPtr(),
+ m_params.t_coords.size(), ParallelDescriptor::IOProcessorNumber());
+ ParallelDescriptor::Bcast(m_params.x_coords.dataPtr(),
+ m_params.x_coords.size(), ParallelDescriptor::IOProcessorNumber());
+ ParallelDescriptor::Bcast(m_params.y_coords.dataPtr(),
+ m_params.y_coords.size(), ParallelDescriptor::IOProcessorNumber());
+ ParallelDescriptor::Barrier();
+}
+
+std::pair<int,int>
+FromTXYEFileLaserProfile::find_left_right_time_indices(amrex::Real t) const
+{
+ int idx_t_right;
+ if(m_params.is_grid_uniform){
+ const auto t_min = m_params.t_coords.front();
+ const auto t_max = m_params.t_coords.back();
+ const auto temp_idx_t_right = static_cast<int>(
+ ceil( (m_params.nt-1)*(t-t_min)/(t_max-t_min)));
+ idx_t_right = max(min(temp_idx_t_right, m_params.nt-1),1);
+ }
+ else{
+ idx_t_right = std::distance(m_params.t_coords.begin(),
+ std::upper_bound(m_params.t_coords.begin(),
+ m_params.t_coords.end(), t));
+ }
+ return std::make_pair(idx_t_right-1, idx_t_right);
+}
+
+void
+FromTXYEFileLaserProfile::read_data_t_chuck(int t_begin, int t_end)
+{
+ amrex::Print() <<
+ "Reading [" << t_begin << ", " << t_end <<
+ ") data chunk from " << m_params.txye_file_name << "\n";
+
+ //Indices of the first and last timestep to read
+ auto i_first = max(0, t_begin);
+ auto i_last = min(t_end-1, m_params.nt-1);
+ if(i_last-i_first+1 > m_params.E_data.size())
+ Abort("Data chunk to read from file is too large");
+
+ if(ParallelDescriptor::IOProcessor()){
+ //Read data chunk
+ std::ifstream inp(m_params.txye_file_name, std::ios::binary);
+ if(!inp) Abort("Failed to open txye file");
+ auto skip_amount = 1 +
+ 3*sizeof(uint32_t) +
+ m_params.t_coords.size()*sizeof(double) +
+ m_params.x_coords.size()*sizeof(double) +
+ m_params.y_coords.size()*sizeof(double) +
+ sizeof(double)*t_begin*m_params.nx*m_params.ny;
+ inp.ignore(skip_amount);
+ if(!inp) Abort("Failed to read field data from txye file");
+ const int read_size = (i_last - i_first + 1)*
+ m_params.nx*m_params.ny;
+ Vector<double> buf_e(read_size);
+ inp.read(reinterpret_cast<char*>(buf_e.dataPtr()), read_size*sizeof(double));
+ if(!inp) Abort("Failed to read field data from txye file");
+ std::transform(buf_e.begin(), buf_e.end(), m_params.E_data.begin(),
+ [](auto x) {return static_cast<amrex::Real>(x);} );
+ }
+
+ //Broadcast E_data
+ ParallelDescriptor::Bcast(m_params.E_data.dataPtr(),
+ m_params.E_data.size(), ParallelDescriptor::IOProcessorNumber());
+ ParallelDescriptor::Barrier();
+
+ //Update first and last indices
+ m_params.first_time_index = i_first;
+ m_params.last_time_index = i_last;
+}
+
+void
+FromTXYEFileLaserProfile::internal_fill_amplitude_uniform(
+ const int idx_t_left,
+ const int np,
+ Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp,
+ Real t, Real * AMREX_RESTRICT const amplitude) const
+{
+ // Copy member variables to tmp copies
+ // and get pointers to underlying data for GPU.
+ const auto tmp_e_max = m_common_params.e_max;
+ const auto tmp_x_min = m_params.x_coords.front();
+ const auto tmp_x_max = m_params.x_coords.back();
+ const auto tmp_y_min = m_params.y_coords.front();
+ const auto tmp_y_max = m_params.y_coords.back();
+ const auto tmp_nx = m_params.nx;
+#if (AMREX_SPACEDIM == 3)
+ const auto tmp_ny = m_params.ny;
+#endif
+ const auto p_E_data = m_params.E_data.dataPtr();
+ const auto tmp_idx_first_time = m_params.first_time_index;
+ const int idx_t_right = idx_t_left+1;
+ const auto t_left = idx_t_left*
+ (m_params.t_coords.back()-m_params.t_coords.front())/(m_params.nt-1) +
+ m_params.t_coords.front();
+ const auto t_right = idx_t_right*
+ (m_params.t_coords.back()-m_params.t_coords.front())/(m_params.nt-1) +
+ m_params.t_coords.front();
+
+ // Loop through the macroparticle to calculate the proper amplitude
+ amrex::ParallelFor(
+ np,
+ [=] AMREX_GPU_DEVICE (int i) {
+ //Amplitude is zero if we are out of bounds
+ if (Xp[i] <= tmp_x_min || Xp[i] >= tmp_x_max){
+ amplitude[i] = 0.0_rt;
+ return;
+ }
+#if (AMREX_SPACEDIM == 3)
+ if (Yp[i] <= tmp_y_min || Yp[i] >= tmp_y_max){
+ amplitude[i] = 0.0_rt;
+ return;
+ }
+#endif
+ //Find indices and coordinates along x
+ const int temp_idx_x_right = static_cast<int>(
+ ceil((tmp_nx-1)*(Xp[i]- tmp_x_min)/(tmp_x_max-tmp_x_min)));
+ const int idx_x_right =
+ max(min(temp_idx_x_right,tmp_nx-1),static_cast<int>(1));
+ const int idx_x_left = idx_x_right - 1;
+ const auto x_0 =
+ idx_x_left*(tmp_x_max-tmp_x_min)/(tmp_nx-1) + tmp_x_min;
+ const auto x_1 =
+ idx_x_right*(tmp_x_max-tmp_x_min)/(tmp_nx-1) + tmp_x_min;
+
+#if (AMREX_SPACEDIM == 2)
+ //Interpolate amplitude
+ const auto idx = [=](int i, int j){
+ return (i-tmp_idx_first_time) * tmp_nx + j;
+ };
+ amplitude[i] = WarpXUtilAlgo::bilinear_interp(
+ t_left, t_right,
+ x_0, x_1,
+ p_E_data[idx(idx_t_left, idx_x_left)],
+ p_E_data[idx(idx_t_left, idx_x_right)],
+ p_E_data[idx(idx_t_right, idx_x_left)],
+ p_E_data[idx(idx_t_right, idx_x_right)],
+ t, Xp[i])*tmp_e_max;
+
+#elif (AMREX_SPACEDIM == 3)
+ //Find indices and coordinates along y
+ const int temp_idx_y_right = static_cast<int>(
+ ceil((tmp_ny-1)*(Yp[i]- tmp_y_min)/(tmp_y_max-tmp_y_min)));
+ const int idx_y_right =
+ max(min(temp_idx_y_right,tmp_ny-1),static_cast<int>(1));
+ const int idx_y_left = idx_y_right - 1;
+ const auto y_0 =
+ idx_y_left*(tmp_y_max-tmp_y_min)/(tmp_ny-1) + tmp_y_min;
+ const auto y_1 =
+ idx_y_right*(tmp_y_max-tmp_y_min)/(tmp_ny-1) + tmp_y_min;
+
+ //Interpolate amplitude
+ const auto idx = [=](int i, int j, int k){
+ return
+ (i-tmp_idx_first_time)*tmp_nx*tmp_ny+
+ j*tmp_ny + k;
+ };
+ amplitude[i] = WarpXUtilAlgo::trilinear_interp(
+ t_left, t_right,
+ x_0, x_1,
+ y_0, y_1,
+ p_E_data[idx(idx_t_left, idx_x_left, idx_y_left)],
+ p_E_data[idx(idx_t_left, idx_x_left, idx_y_right)],
+ p_E_data[idx(idx_t_left, idx_x_right, idx_y_left)],
+ p_E_data[idx(idx_t_left, idx_x_right, idx_y_right)],
+ p_E_data[idx(idx_t_right, idx_x_left, idx_y_left)],
+ p_E_data[idx(idx_t_right, idx_x_left, idx_y_right)],
+ p_E_data[idx(idx_t_right, idx_x_right, idx_y_left)],
+ p_E_data[idx(idx_t_right, idx_x_right, idx_y_right)],
+ t, Xp[i], Yp[i])*tmp_e_max;
+#endif
+ }
+ );
+}
+
+void
+FromTXYEFileLaserProfile::internal_fill_amplitude_nonuniform(
+ const int idx_t_left,
+ const int np,
+ Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp,
+ Real t, Real * AMREX_RESTRICT const amplitude) const
+{
+ // Copy member variables to tmp copies
+ // and get pointers to underlying data for GPU.
+ const auto tmp_e_max = m_common_params.e_max;
+ const auto tmp_x_min = m_params.x_coords.front();
+ const auto tmp_x_max = m_params.x_coords.back();
+ const auto tmp_y_min = m_params.y_coords.front();
+ const auto tmp_y_max = m_params.y_coords.back();
+ const auto p_x_coords = m_params.x_coords.dataPtr();
+ const int tmp_x_coords_size = static_cast<int>(m_params.x_coords.size());
+ const auto p_y_coords = m_params.y_coords.dataPtr();
+ const int tmp_y_coords_size = static_cast<int>(m_params.y_coords.size());
+ const auto p_E_data = m_params.E_data.dataPtr();
+ const auto tmp_idx_first_time = m_params.first_time_index;
+ const int idx_t_right = idx_t_left+1;
+ const auto t_left = m_params.t_coords[idx_t_left];
+ const auto t_right = m_params.t_coords[idx_t_right];
+
+ // Loop through the macroparticle to calculate the proper amplitude
+ amrex::ParallelFor(
+ np,
+ [=] AMREX_GPU_DEVICE (int i) {
+ //Amplitude is zero if we are out of bounds
+ if (Xp[i] <= tmp_x_min || Xp[i] >= tmp_x_max){
+ amplitude[i] = 0.0_rt;
+ return;
+ }
+#if (AMREX_SPACEDIM == 3)
+ if (Yp[i] <= tmp_y_min || Yp[i] >= tmp_y_max){
+ amplitude[i] = 0.0_rt;
+ return;
+ }
+#endif
+
+ //Find indices along x
+ auto const p_x_right = WarpXUtilAlgo::upper_bound(
+ p_x_coords, p_x_coords+tmp_x_coords_size, Xp[i]);
+ const int idx_x_right = p_x_right - p_x_coords;
+ const int idx_x_left = idx_x_right - 1;
+
+#if (AMREX_SPACEDIM == 2)
+ //Interpolate amplitude
+ const auto idx = [=](int i, int j){
+ return (i-tmp_idx_first_time) * tmp_x_coords_size + j;
+ };
+ amplitude[i] = WarpXUtilAlgo::bilinear_interp(
+ t_left, t_right,
+ p_x_coords[idx_x_left], p_x_coords[idx_x_right],
+ p_E_data[idx(idx_t_left, idx_x_left)],
+ p_E_data[idx(idx_t_left, idx_x_right)],
+ p_E_data[idx(idx_t_right, idx_x_left)],
+ p_E_data[idx(idx_t_right, idx_x_right)],
+ t, Xp[i])*tmp_e_max;
+
+#elif (AMREX_SPACEDIM == 3)
+ //Find indices along y
+ auto const p_y_right = WarpXUtilAlgo::upper_bound(
+ p_y_coords, p_y_coords+tmp_y_coords_size, Yp[i]);
+ const int idx_y_right = p_y_right - p_y_coords;
+ const int idx_y_left = idx_y_right - 1;
+
+ //Interpolate amplitude
+ const auto idx = [=](int i, int j, int k){
+ return
+ (i-tmp_idx_first_time)*tmp_x_coords_size*tmp_y_coords_size+
+ j*tmp_y_coords_size + k;
+ };
+ amplitude[i] = WarpXUtilAlgo::trilinear_interp(
+ t_left, t_right,
+ p_x_coords[idx_x_left], p_x_coords[idx_x_right],
+ p_y_coords[idx_y_left], p_y_coords[idx_y_right],
+ p_E_data[idx(idx_t_left, idx_x_left, idx_y_left)],
+ p_E_data[idx(idx_t_left, idx_x_left, idx_y_right)],
+ p_E_data[idx(idx_t_left, idx_x_right, idx_y_left)],
+ p_E_data[idx(idx_t_left, idx_x_right, idx_y_right)],
+ p_E_data[idx(idx_t_right, idx_x_left, idx_y_left)],
+ p_E_data[idx(idx_t_right, idx_x_left, idx_y_right)],
+ p_E_data[idx(idx_t_right, idx_x_right, idx_y_left)],
+ p_E_data[idx(idx_t_right, idx_x_right, idx_y_right)],
+ t, Xp[i], Yp[i])*tmp_e_max;
+#endif
+ }
+ );
+}
diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp
index a0b5dd855..18bdec4a7 100644
--- a/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp
+++ b/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp
@@ -72,7 +72,7 @@ GaussianLaserProfile::init (
void
GaussianLaserProfile::fill_amplitude (
const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp,
- Real t, Real * AMREX_RESTRICT const amplitude)
+ Real t, Real * AMREX_RESTRICT const amplitude) const
{
Complex I(0,1);
// Calculate a few factors which are independent of the macroparticle
diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp
index 55374c5ea..7fe75cf56 100644
--- a/Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp
+++ b/Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp
@@ -35,7 +35,7 @@ HarrisLaserProfile::init (
void
HarrisLaserProfile::fill_amplitude (
const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp,
- Real t, Real * AMREX_RESTRICT const amplitude)
+ Real t, Real * AMREX_RESTRICT const amplitude) const
{
// This function uses the Harris function as the temporal profile of the pulse
const Real omega0 =
diff --git a/Source/Laser/LaserProfilesImpl/Make.package b/Source/Laser/LaserProfilesImpl/Make.package
index 32284c4e4..2fef27b9f 100644
--- a/Source/Laser/LaserProfilesImpl/Make.package
+++ b/Source/Laser/LaserProfilesImpl/Make.package
@@ -1,6 +1,7 @@
CEXE_sources += LaserProfileGaussian.cpp
CEXE_sources += LaserProfileHarris.cpp
CEXE_sources += LaserProfileFieldFunction.cpp
+CEXE_sources += LaserProfileFromTXYEFile.cpp
INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Laser/LaserProfilesImpl
VPATH_LOCATIONS += $(WARPX_HOME)/Source/Laser/LaserProfilesImpl
diff --git a/Source/Parser/Make.package b/Source/Parser/Make.package
index 5ce02cbda..15115c138 100644
--- a/Source/Parser/Make.package
+++ b/Source/Parser/Make.package
@@ -5,6 +5,7 @@ CEXE_sources += WarpXParser.cpp
CEXE_headers += WarpXParser.H
CEXE_headers += GpuParser.H
CEXE_sources += GpuParser.cpp
+CEXE_headers += WarpXParserWrapper.H
INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Parser
VPATH_LOCATIONS += $(WARPX_HOME)/Source/Parser
diff --git a/Source/Parser/WarpXParserWrapper.H b/Source/Parser/WarpXParserWrapper.H
new file mode 100644
index 000000000..2dd7f72c7
--- /dev/null
+++ b/Source/Parser/WarpXParserWrapper.H
@@ -0,0 +1,35 @@
+#ifndef WARPX_PARSER_WRAPPER_H_
+#define WARPX_PARSER_WRAPPER_H_
+
+#include <WarpXParser.H>
+#include <AMReX_Gpu.H>
+#include <GpuParser.H>
+/**
+ * \brief
+ * The ParserWrapper struct is constructed to safely use the GpuParser class
+ * that can essentially be though of as a raw pointer. The GpuParser does
+ * not have an explicit destructor and the AddPlasma subroutines handle the GpuParser
+ * in a safe way. The ParserWrapper struct is used to avoid memory leak
+ * in the EB parser functions.
+ */
+struct ParserWrapper
+ : public amrex::Gpu::Managed
+{
+ ParserWrapper (WarpXParser const& a_parser) noexcept
+ : m_parser(a_parser) {}
+
+ ~ParserWrapper() {
+ m_parser.clear();
+ }
+
+ AMREX_GPU_HOST_DEVICE
+ amrex::Real
+ getField (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
+ {
+ return m_parser(x,y,z);
+ }
+
+ GpuParser m_parser;
+};
+
+#endif
diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp
index c577da7f3..e05a64bfe 100644
--- a/Source/Utils/WarpXMovingWindow.cpp
+++ b/Source/Utils/WarpXMovingWindow.cpp
@@ -1,4 +1,5 @@
#include <WarpX.H>
+#include <WarpXUtil.H>
#include <WarpXConst.H>
using namespace amrex;
@@ -97,10 +98,25 @@ WarpX::MoveWindow (bool move_j)
// Shift each component of vector fields (E, B, j)
for (int dim = 0; dim < 3; ++dim) {
-
// Fine grid
- shiftMF(*Bfield_fp[lev][dim], geom[lev], num_shift, dir, B_external_grid[dim]);
- shiftMF(*Efield_fp[lev][dim], geom[lev], num_shift, dir, E_external_grid[dim]);
+ ParserWrapper *Bfield_parser;
+ ParserWrapper *Efield_parser;
+ bool use_Bparser = false;
+ bool use_Eparser = false;
+ if (B_ext_grid_s == "parse_b_ext_grid_function") {
+ use_Bparser = true;
+ if (dim == 0) Bfield_parser = Bxfield_parser.get();
+ if (dim == 1) Bfield_parser = Byfield_parser.get();
+ if (dim == 2) Bfield_parser = Bzfield_parser.get();
+ }
+ if (E_ext_grid_s == "parse_e_ext_grid_function") {
+ use_Eparser = true;
+ if (dim == 0) Efield_parser = Exfield_parser.get();
+ if (dim == 1) Efield_parser = Eyfield_parser.get();
+ if (dim == 2) Efield_parser = Ezfield_parser.get();
+ }
+ shiftMF(*Bfield_fp[lev][dim], geom[lev], num_shift, dir, B_external_grid[dim], use_Bparser, Bfield_parser);
+ shiftMF(*Efield_fp[lev][dim], geom[lev], num_shift, dir, E_external_grid[dim], use_Eparser, Efield_parser);
if (move_j) {
shiftMF(*current_fp[lev][dim], geom[lev], num_shift, dir);
}
@@ -112,9 +128,9 @@ WarpX::MoveWindow (bool move_j)
}
if (lev > 0) {
- // Coarse grid
- shiftMF(*Bfield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, B_external_grid[dim]);
- shiftMF(*Efield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, E_external_grid[dim]);
+ // coarse grid
+ shiftMF(*Bfield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, B_external_grid[dim], use_Bparser, Bfield_parser);
+ shiftMF(*Efield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, E_external_grid[dim], use_Eparser, Efield_parser);
shiftMF(*Bfield_aux[lev][dim], geom[lev], num_shift, dir);
shiftMF(*Efield_aux[lev][dim], geom[lev], num_shift, dir);
if (move_j) {
@@ -132,7 +148,7 @@ WarpX::MoveWindow (bool move_j)
// Shift scalar component F for dive cleaning
if (do_dive_cleaning) {
// Fine grid
- shiftMF(*F_fp[lev], geom[lev], num_shift, dir);
+ shiftMF(*F_fp[lev], geom[lev], num_shift, dir);
if (do_pml && pml[lev]->ok()) {
MultiFab* pml_F = pml[lev]->GetF_fp();
shiftMF(*pml_F, geom[lev], num_shift, dir);
@@ -204,7 +220,7 @@ WarpX::MoveWindow (bool move_j)
void
WarpX::shiftMF (MultiFab& mf, const Geometry& geom, int num_shift, int dir,
- amrex::Real external_field)
+ amrex::Real external_field, bool useparser, ParserWrapper *field_parser)
{
BL_PROFILE("WarpX::shiftMF()");
const BoxArray& ba = mf.boxArray();
@@ -246,20 +262,57 @@ WarpX::shiftMF (MultiFab& mf, const Geometry& geom, int num_shift, int dir,
shiftiv[dir] = num_shift;
Dim3 shift = shiftiv.dim3();
+ const RealBox& real_box = geom.ProbDomain();
+ const auto dx = geom.CellSizeArray();
+
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
+
+
for (MFIter mfi(tmpmf); mfi.isValid(); ++mfi )
{
auto const& dstfab = mf.array(mfi);
auto const& srcfab = tmpmf.array(mfi);
const Box& outbox = mfi.fabbox() & adjBox;
+
if (outbox.ok()) {
- AMREX_PARALLEL_FOR_4D ( outbox, nc, i, j, k, n,
- {
- srcfab(i,j,k,n) = external_field;
- });
+ if (useparser == false) {
+ AMREX_PARALLEL_FOR_4D ( outbox, nc, i, j, k, n,
+ {
+ srcfab(i,j,k,n) = external_field;
+ });
+ } else if (useparser == true) {
+ // index type of the src mf
+ auto const& mf_IndexType = (tmpmf).ixType();
+ IntVect mf_type(AMREX_D_DECL(0,0,0));
+ for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+ mf_type[idim] = mf_IndexType.nodeCentered(idim);
+ }
+
+ amrex::ParallelFor (outbox, nc,
+ [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+ {
+ // Compute x,y,z co-ordinates based on index type of mf
+ Real fac_x = (1.0 - mf_type[0]) * dx[0]*0.5;
+ Real x = i*dx[0] + real_box.lo(0) + fac_x;
+#if (AMREX_SPACEDIM==2)
+ Real y = 0.0;
+ Real fac_z = (1.0 - mf_type[1]) * dx[1]*0.5;
+ Real z = j*dx[1] + real_box.lo(1) + fac_z;
+#else
+ Real fac_y = (1.0 - mf_type[1]) * dx[1]*0.5;
+ Real y = j*dx[1] + real_box.lo(1) + fac_y;
+ Real fac_z = (1.0 - mf_type[2]) * dx[2]*0.5;
+ Real z = k*dx[2] + real_box.lo(2) + fac_z;
+#endif
+ srcfab(i,j,k,n) = field_parser->getField(x,y,z);
+ }
+ , amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double)*3
+ );
+ }
+
}
Box dstBox = mf[mfi].box();
diff --git a/Source/Utils/WarpXUtil.H b/Source/Utils/WarpXUtil.H
index 195e309cf..e7b2ef196 100644
--- a/Source/Utils/WarpXUtil.H
+++ b/Source/Utils/WarpXUtil.H
@@ -4,6 +4,8 @@
#include <AMReX_REAL.H>
#include <AMReX_Vector.H>
#include <AMReX_MultiFab.H>
+#include <AMReX_ParmParse.H>
+#include <WarpXParser.H>
#include <string>
@@ -15,14 +17,108 @@ void ConvertLabParamsToBoost();
void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin,
amrex::Real zmax);
+/**
+ * \brief Parse a string (typically a mathematical expression) from the
+ * input file and store it into a variable.
+ *
+ * \param ParmParse pp used to read the query_string pp.<function>=string
+ * \param parmparse_string String used to initialize ParmParse
+ * \param query_string ParmParse.query will look for this string
+ * \param stored_string variable in which the string to parse is stored
+ */
+void Store_parserString(amrex::ParmParse &pp, std::string query_string,
+ std::string& stored_string);
+
namespace WarpXUtilIO{
- /**
- * A helper function to write binary data on disk.
- * @param[in] filename where to write
- * @param[in] data Vector containing binary data to write on disk
- * return true if it succeeds, false otherwise
- */
- bool WriteBinaryDataOnFile(std::string filename, const amrex::Vector<char>& data);
+/**
+ * A helper function to write binary data on disk.
+ * @param[in] filename where to write
+ * @param[in] data Vector containing binary data to write on disk
+ * return true if it succeeds, false otherwise
+ */
+bool WriteBinaryDataOnFile(std::string filename, const amrex::Vector<char>& data);
+}
+
+namespace WarpXUtilAlgo{
+
+/** \brief Returns a pointer to the first element in the range [first, last) that is greater than val
+ *
+ * A re-implementation of the upper_bound algorithm suitable for GPU kernels.
+ *
+ * @param first: pointer to left limit of the range to consider
+ * @param last: pointer to right limit of the range to consider
+ * @param val: value to compare the elements of [first, last) to
+ */
+template<typename T> AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+const T* upper_bound(const T* first, const T* last, const T& val)
+{
+ const T* it;
+ size_t count, step;
+ count = last-first;
+ while(count>0){
+ it = first;
+ step = count/2;
+ it += step;
+ if (!(val<*it)){
+ first = ++it;
+ count -= step + 1;
+ }
+ else{
+ count = step;
+ }
+ }
+ return first;
}
+/** \brief Performs a linear interpolation
+ *
+ * Performs a linear interpolation at x given the 2 points
+ * (x0, f0) and (x1, f1)
+ */
+template<typename T> AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+T linear_interp(T x0, T x1, T f0, T f1, T x)
+{
+ return ((x1-x)*f0 + (x-x0)*f1)/(x1-x0);
+}
+
+/** \brief Performs a bilinear interpolation
+ *
+ * Performs a bilinear interpolation at (x,y) given the 4 points
+ * (x0, y0, f00), (x0, y1, f01), (x1, y0, f10), (x1, y1, f11).
+ */
+template<typename T> AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+T bilinear_interp(T x0, T x1, T y0, T y1, T f00, T f01, T f10, T f11, T x, T y)
+{
+ const T fx0 = linear_interp(x0, x1, f00, f10, x);
+ const T fx1 = linear_interp(x0, x1, f01, f11, x);
+ return linear_interp(y0, y1, fx0, fx1, y);
+}
+
+/** \brief Performs a trilinear interpolation
+ *
+ * Performs a trilinear interpolation at (x,y,z) given the 8 points
+ * (x0, y0, z0, f000), (x0, y0, z1, f001), (x0, y1, z0, f010), (x0, y1, z1, f011),
+ * (x1, y0, z0, f100), (x1, y0, z1, f101), (x1, y1, z0, f110), (x1, y1, z1, f111)
+ */
+template<typename T> AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+T trilinear_interp(T x0, T x1,T y0, T y1, T z0, T z1,
+ T f000, T f001, T f010, T f011, T f100, T f101, T f110, T f111,
+ T x, T y, T z)
+{
+ const T fxy0 = bilinear_interp(
+ x0, x1, y0, y1, f000, f010, f100, f110, x, y);
+ const T fxy1 = bilinear_interp(
+ x0, x1, y0, y1, f001, f011, f101, f111, x, y);
+ return linear_interp(z0, z1, fxy0, fxy1, z);
+}
+
+}
+
+/**
+* \brief Initialize a WarpXParser object from a string containing a math expression
+*
+* \param parse_function String to read to initialize the parser.
+*/
+WarpXParser makeParser (std::string const& parse_function);
+
#endif //WARPX_UTILS_H_
diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp
index 8764a09c6..e9fb958fd 100644
--- a/Source/Utils/WarpXUtil.cpp
+++ b/Source/Utils/WarpXUtil.cpp
@@ -164,3 +164,44 @@ namespace WarpXUtilIO{
}
}
+void Store_parserString(amrex::ParmParse& pp, std::string query_string,
+ std::string& stored_string)
+{
+
+ char cstr[query_string.size()+1];
+ strcpy(cstr, query_string.c_str());
+
+ std::vector<std::string> f;
+ pp.getarr(cstr, f);
+ stored_string.clear();
+ for (auto const& s : f) {
+ stored_string += s;
+ }
+ f.clear();
+
+}
+
+
+WarpXParser makeParser (std::string const& parse_function)
+{
+ WarpXParser parser(parse_function);
+ parser.registerVariables({"x","y","z"});
+ ParmParse pp("my_constants");
+ std::set<std::string> symbols = parser.symbols();
+ symbols.erase("x");
+ symbols.erase("y");
+ symbols.erase("z");
+ for (auto it = symbols.begin(); it != symbols.end(); ) {
+ Real v;
+ if (pp.query(it->c_str(), v)) {
+ parser.setConstant(*it, v);
+ it = symbols.erase(it);
+ } else {
+ ++it;
+ }
+ }
+ for (auto const& s : symbols) {
+ amrex::Abort("makeParser::Unknown symbol "+s);
+ }
+ return parser;
+}
diff --git a/Source/WarpX.H b/Source/WarpX.H
index decde3dc2..3bab73833 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -24,6 +24,7 @@
#include <AMReX_LayoutData.H>
#include <AMReX_Interpolater.H>
#include <AMReX_FillPatchUtil.H>
+#include <WarpXParserWrapper.H>
#ifdef _OPENMP
# include <omp.h>
@@ -73,8 +74,10 @@ public:
MultiParticleContainer& GetPartContainer () { return *mypc; }
- static void shiftMF(amrex::MultiFab& mf, const amrex::Geometry& geom, int num_shift, int dir,
- amrex::Real external_field = 0.);
+ static void shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom,
+ int num_shift, int dir, amrex::Real external_field=0.0,
+ bool useparser = false,
+ ParserWrapper *field_parser=nullptr);
static void GotoNextLine (std::istream& is);
@@ -86,6 +89,28 @@ public:
static amrex::Vector<amrex::Real> E_external_grid;
static amrex::Vector<amrex::Real> B_external_grid;
+ // Initialization Type for External E and B
+ static std::string B_ext_grid_s;
+ static std::string E_ext_grid_s;
+
+ // Parser for B_external on the grid
+ static std::string str_Bx_ext_grid_function;
+ static std::string str_By_ext_grid_function;
+ static std::string str_Bz_ext_grid_function;
+ // Parser for E_external on the grid
+ static std::string str_Ex_ext_grid_function;
+ static std::string str_Ey_ext_grid_function;
+ static std::string str_Ez_ext_grid_function;
+
+ // ParserWrapper for B_external on the grid
+ std::unique_ptr<ParserWrapper> Bxfield_parser;
+ std::unique_ptr<ParserWrapper> Byfield_parser;
+ std::unique_ptr<ParserWrapper> Bzfield_parser;
+ // ParserWrapper for E_external on the grid
+ std::unique_ptr<ParserWrapper> Exfield_parser;
+ std::unique_ptr<ParserWrapper> Eyfield_parser;
+ std::unique_ptr<ParserWrapper> Ezfield_parser;
+
// Algorithms
static long current_deposition_algo;
static long charge_deposition_algo;
@@ -323,6 +348,49 @@ public:
protected:
+ /**
+ * \brief
+ * This function initializes E, B, rho, and F, at all the levels
+ * of the multifab. rho and F are initialized with 0.
+ * The E and B fields are initialized using user-defined inputs.
+ * The initialization type is set using "B_ext_grid_init_style"
+ * and "E_ext_grid_init_style". The initialization style is set to "default"
+ * if not explicitly defined by the user, and the E and B fields are
+ * initialized with E_external_grid and B_external_grid, respectively, each with
+ * a default value of 0.
+ * If the initialization type for the E and B field is "constant",
+ * then, the E and B fields at all the levels are initialized with
+ * user-defined values for E_external_grid and B_external_grid.
+ * If the initialization type for B-field is set to
+ * "parse_B_ext_grid_function", then, the parser is used to read
+ * Bx_external_grid_function(x,y,z), By_external_grid_function(x,y,z),
+ * and Bz_external_grid_function(x,y,z).
+ * Similarly, if the E-field initialization type is set to
+ * "parse_E_ext_grid_function", then, the parser is used to read
+ * Ex_external_grid_function(x,y,z), Ey_external_grid_function(x,y,z),
+ * and Ex_external_grid_function(x,y,z). The parser for the E and B
+ * initialization assumes that the function has three independent
+ * variables, at max, namely, x, y, z. However, any number of constants
+ * can be used in the function used to define the E and B fields on the grid.
+ */
+ void InitLevelData (int lev, amrex::Real time);
+
+ /**
+ * \brief
+ * This function initializes the E and B fields on each level
+ * using the parser and the user-defined function for the external fields.
+ * The subroutine will parse the x_/y_z_external_grid_function and
+ * then, the B or E multifab is initialized based on the (x,y,z) position
+ * on the staggered yee-grid or cell-centered grid.
+ */
+ void InitializeExternalFieldsOnGridUsingParser (
+ amrex::MultiFab *mfx, amrex::MultiFab *mfy, amrex::MultiFab *mfz,
+ ParserWrapper *xfield_parser, ParserWrapper *yfield_parser,
+ ParserWrapper *zfield_parser, amrex::IntVect x_nodal_flag,
+ amrex::IntVect y_nodal_flag, amrex::IntVect z_nodal_flag,
+ const int lev);
+
+
//! Tagging cells for refinement
virtual void ErrorEst (int lev, amrex::TagBoxArray& tags, amrex::Real time, int /*ngrow*/) final;
@@ -410,7 +478,6 @@ private:
void AllocLevelData (int lev, const amrex::BoxArray& new_grids,
const amrex::DistributionMapping& new_dmap);
- void InitLevelData (int lev, amrex::Real time);
void InitFromCheckpoint ();
void PostRestart ();
@@ -665,4 +732,5 @@ private:
int insitu_pin_mesh;
};
+
#endif
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 377d103d1..48b4bbd55 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -30,6 +30,18 @@ Vector<Real> WarpX::E_external_particle(3, 0.0);
Vector<Real> WarpX::E_external_grid(3, 0.0);
Vector<Real> WarpX::B_external_grid(3, 0.0);
+std::string WarpX::B_ext_grid_s = "default";
+std::string WarpX::E_ext_grid_s = "default";
+
+// Parser for B_external on the grid
+std::string WarpX::str_Bx_ext_grid_function;
+std::string WarpX::str_By_ext_grid_function;
+std::string WarpX::str_Bz_ext_grid_function;
+// Parser for E_external on the grid
+std::string WarpX::str_Ex_ext_grid_function;
+std::string WarpX::str_Ey_ext_grid_function;
+std::string WarpX::str_Ez_ext_grid_function;
+
int WarpX::do_moving_window = 0;
int WarpX::moving_window_dir = -1;
Real WarpX::moving_window_v = std::numeric_limits<amrex::Real>::max();
@@ -309,9 +321,6 @@ WarpX::ReadParameters ()
pp.queryarr("B_external_particle", B_external_particle);
pp.queryarr("E_external_particle", E_external_particle);
- pp.queryarr("E_external_grid", E_external_grid);
- pp.queryarr("B_external_grid", B_external_grid);
-
pp.query("do_moving_window", do_moving_window);
if (do_moving_window)
{
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-14 10:12:05 +0000