aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Docs/source/usage/parameters.rst10
-rw-r--r--Examples/Tests/spectral_staggered/inputs_hybrid_2d4
-rw-r--r--Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd_hybrid.json38
-rw-r--r--Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd_hybrid.json42
-rw-r--r--Regression/WarpX-tests.ini4
-rw-r--r--Source/Evolve/WarpXEvolve.cpp15
-rw-r--r--Source/Initialization/WarpXInitData.cpp5
-rw-r--r--Source/Parallelization/WarpXComm.cpp118
-rw-r--r--Source/Parallelization/WarpXComm_K.H148
-rw-r--r--Source/WarpX.H79
-rw-r--r--Source/WarpX.cpp200
11 files changed, 459 insertions, 204 deletions
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index f90ec5485..6903423bb 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -1228,8 +1228,14 @@ Numerics and algorithms
The default behavior should not normally be changed.
At present, this parameter is intended mainly for testing and development purposes.
-* ``interpolation.field_gathering_nox``, ``interpolation.field_gathering_noy``, ``interpolation.field_gathering_noz`` (default: ``2`` in all directions)
- The order of the interpolation used with staggered grids (``warpx.do_nodal = 0``) and momentum-conserving field gathering (``algo.field_gathering = momentum-conserving``) to interpolate the electric and magnetic fields from the cell centers to the cell nodes, before gathering the fields from the cell nodes to the particle positions. High-order interpolation (order 8 in each direction, at least) is necessary to ensure stability in typical LWFA boosted-frame simulations using the Galilean PSATD or comoving PSATD schemes. This arbitrary-order interpolation is used only when the PSATD solver is used for Maxwell's equations. With the FDTD solver, basic linear interpolation is used instead.
+* ``interpolation.field_centering_nox``, ``interpolation.field_centering_noy``, ``interpolation.field_centering_noz`` (default: ``2`` in all directions)
+ The order of interpolation used with staggered grids (``warpx.do_nodal = 0``) and momentum-conserving field gathering (``algo.field_gathering = momentum-conserving``) to interpolate the electric and magnetic fields from the cell centers to the cell nodes, before gathering the fields from the cell nodes to the particle positions. High-order interpolation (order 8 in each direction, at least) is necessary to ensure stability in typical LWFA boosted-frame simulations using the Galilean PSATD or comoving PSATD schemes. This finite-order interpolation is used only when the PSATD solver is used for Maxwell's equations. With the FDTD solver, basic linear interpolation is used instead.
+
+* ``interpolation.current_centering_nox``, ``interpolation.current_centering_noy``, ``interpolation.current_centering_noz`` (default: ``2`` in all directions)
+ The order of interpolation used to center the currents from nodal to staggered grids (if ``warpx.do_current_centering = 1``), before pushing the Maxwell fields on staggered grids. This finite-order interpolation is used only when the PSATD solver is used for Maxwell's equations. With the FDTD solver, basic linear interpolation is used instead.
+
+* ``warpx.do_current_centering`` (`0` or `1` ; default: 0)
+ If true, the current is deposited on a nodal grid and then centered to a staggered grid (Yee grid), using finite-order interpolation. If ``warpx.do_nodal = 1``, the Maxwell fields are pushed on a nodal grid, it is not necessary to center the currents to a staggered grid, and we set therefore ``warpx.do_current_centering = 0`` automatically, overwriting the user-defined input.
* ``warpx.do_dive_cleaning`` (`0` or `1` ; default: 0)
Whether to use modified Maxwell equations that progressively eliminate
diff --git a/Examples/Tests/spectral_staggered/inputs_hybrid_2d b/Examples/Tests/spectral_staggered/inputs_hybrid_2d
index 6e33e1171..c138a399f 100644
--- a/Examples/Tests/spectral_staggered/inputs_hybrid_2d
+++ b/Examples/Tests/spectral_staggered/inputs_hybrid_2d
@@ -36,8 +36,8 @@ warpx.use_filter = 1
warpx.serialize_ics = 1
warpx.verbose = 1
-interpolation.field_gathering_nox = 8
-interpolation.field_gathering_noz = 8
+interpolation.field_centering_nox = 8
+interpolation.field_centering_noz = 8
particles.species_names = electrons ions beam
particles.use_fdtd_nci_corr = 0
diff --git a/Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd_hybrid.json b/Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd_hybrid.json
index 5295441d7..c5d7e10ce 100644
--- a/Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd_hybrid.json
+++ b/Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd_hybrid.json
@@ -2,32 +2,32 @@
"electrons": {
"particle_cpu": 32768.0,
"particle_id": 1073774592.0,
- "particle_momentum_x": 2.695748345809732e-21,
- "particle_momentum_y": 2.3028925989947628e-21,
- "particle_momentum_z": 1.780767583794476e-16,
- "particle_position_x": 405588.6300079636,
- "particle_position_y": 20127109.122946262,
+ "particle_momentum_x": 2.527239041054006e-21,
+ "particle_momentum_y": 2.3029995759007875e-21,
+ "particle_momentum_z": 1.7807674965649623e-16,
+ "particle_position_x": 405588.62581602746,
+ "particle_position_y": 20127109.123043947,
"particle_weight": 6.917460794691972e+17
},
"ions": {
"particle_cpu": 32768.0,
"particle_id": 3221258240.0,
- "particle_momentum_x": 2.636863383941158e-18,
- "particle_momentum_y": 2.6255372112121483e-18,
- "particle_momentum_z": 3.2697610878128845e-13,
- "particle_position_x": 405588.6872870345,
- "particle_position_y": 20127109.122934863,
+ "particle_momentum_x": 2.6368724382965566e-18,
+ "particle_momentum_y": 2.6255372125997566e-18,
+ "particle_momentum_z": 3.2697610879128045e-13,
+ "particle_position_x": 405588.6872877858,
+ "particle_position_y": 20127109.122936234,
"particle_weight": 6.917460794691972e+17
},
"lev=0": {
- "Bx": 0.002181514694303784,
- "By": 0.015312399998106402,
- "Bz": 0.006258536700732526,
- "Ex": 4594891.40715035,
- "Ey": 2030809.91823095,
- "Ez": 362845.95726377127,
- "jx": 428.70648696668235,
- "jy": 819.0133385850143,
- "jz": 69647.2014472035
+ "Bx": 0.0021815053656360347,
+ "By": 0.004179957439645016,
+ "Bz": 0.006258472864965425,
+ "Ex": 1323685.1423194357,
+ "Ey": 2030803.1957688145,
+ "Ez": 306201.9811103587,
+ "jx": 407.69195463935387,
+ "jy": 818.9909953388023,
+ "jz": 4875.487262566359
}
} \ No newline at end of file
diff --git a/Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd_hybrid.json b/Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd_hybrid.json
index 7fda73008..99be8a0b0 100644
--- a/Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd_hybrid.json
+++ b/Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd_hybrid.json
@@ -2,34 +2,34 @@
"electrons": {
"particle_cpu": 0.0,
"particle_id": 536887296.0,
- "particle_momentum_x": 3.985020903846283e-20,
- "particle_momentum_y": 3.8527082179407666e-20,
- "particle_momentum_z": 8.904330117181464e-17,
- "particle_position_x": 158433.3336413721,
- "particle_position_y": 158432.19270800205,
- "particle_position_z": 15724303.898790203,
+ "particle_momentum_x": 3.094828314000131e-20,
+ "particle_momentum_y": 2.986276244935509e-20,
+ "particle_momentum_z": 8.904170376087041e-17,
+ "particle_position_x": 158433.3737028487,
+ "particle_position_y": 158432.1663227203,
+ "particle_position_z": 15724303.92810403,
"particle_weight": 4.082754265421834e+18
},
"ions": {
"particle_cpu": 0.0,
"particle_id": 1610629120.0,
- "particle_momentum_x": 1.3146876520749324e-18,
- "particle_momentum_y": 1.3118515839055124e-18,
- "particle_momentum_z": 1.6348803284082346e-13,
- "particle_position_x": 158433.36778651696,
- "particle_position_y": 158432.13062157098,
- "particle_position_z": 15724303.986645736,
+ "particle_momentum_x": 1.3143999865320922e-18,
+ "particle_momentum_y": 1.3115489851341607e-18,
+ "particle_momentum_z": 1.6348803506523476e-13,
+ "particle_position_x": 158433.36794520426,
+ "particle_position_y": 158432.13055551145,
+ "particle_position_z": 15724303.986785894,
"particle_weight": 4.082754265421834e+18
},
"lev=0": {
- "Bx": 0.10966781374215134,
- "By": 0.11084296672822039,
- "Bz": 0.5914307622478067,
- "Ex": 197527229.28904873,
- "Ey": 189459821.69744563,
- "Ez": 3216201.315611569,
- "jx": 85228.92279339595,
- "jy": 83520.36614829257,
- "jz": 350716.4444535297
+ "Bx": 0.056130529683723604,
+ "By": 0.05812528370481717,
+ "Bz": 0.5178262645776277,
+ "Ex": 171598317.0963192,
+ "Ey": 164544282.67461893,
+ "Ez": 1006956.810225577,
+ "jx": 76053.34653818703,
+ "jy": 74549.16508051571,
+ "jz": 20570.512532683482
}
} \ No newline at end of file
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index 8bcf8dbf3..1f3eb16ed 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -1782,7 +1782,7 @@ tolerance = 1e-6
[averaged_galilean_2d_psatd_hybrid]
buildDir = .
inputFile = Examples/Tests/averaged_galilean/inputs_avg_2d
-runtime_params = amr.max_grid_size_x = 128 amr.max_grid_size_y = 64 warpx.do_nodal = 0 algo.field_gathering = momentum-conserving interpolation.field_gathering_nox = 8 interpolation.field_gathering_noz = 8
+runtime_params = amr.max_grid_size_x = 128 amr.max_grid_size_y = 64 warpx.do_nodal = 0 algo.field_gathering = momentum-conserving interpolation.field_centering_nox = 8 interpolation.field_centering_noz = 8 warpx.do_current_centering = 1 interpolation.current_centering_nox = 8 interpolation.current_centering_noz = 8
dim = 2
addToCompileString = USE_PSATD=TRUE
restartTest = 0
@@ -1818,7 +1818,7 @@ tolerance = 1e-4
[averaged_galilean_3d_psatd_hybrid]
buildDir = .
inputFile = Examples/Tests/averaged_galilean/inputs_avg_3d
-runtime_params = warpx.do_nodal = 0 algo.field_gathering = momentum-conserving interpolation.field_gathering_nox = 8 interpolation.field_gathering_noy = 8 interpolation.field_gathering_noz = 8
+runtime_params = warpx.do_nodal = 0 algo.field_gathering = momentum-conserving interpolation.field_centering_nox = 8 interpolation.field_centering_noy = 8 interpolation.field_centering_noz = 8 warpx.do_current_centering = 1 interpolation.current_centering_nox = 8 interpolation.current_centering_noy = 8 interpolation.current_centering_noz = 8
dim = 3
addToCompileString = USE_PSATD=TRUE
restartTest = 0
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
index 39b51b261..e632409f6 100644
--- a/Source/Evolve/WarpXEvolve.cpp
+++ b/Source/Evolve/WarpXEvolve.cpp
@@ -21,7 +21,6 @@
#include <cmath>
#include <limits>
-
using namespace amrex;
void
@@ -303,6 +302,7 @@ WarpX::OneStep_nosub (Real cur_time)
if (warpx_py_particlescraper) warpx_py_particlescraper();
if (warpx_py_beforedeposition) warpx_py_beforedeposition();
PushParticlesandDepose(cur_time);
+
if (warpx_py_afterdeposition) warpx_py_afterdeposition();
// Synchronize J and rho
@@ -317,7 +317,6 @@ WarpX::OneStep_nosub (Real cur_time)
VayDeposition();
}
-
// At this point, J is up-to-date inside the domain, and E and B are
// up-to-date including enough guard cells for first step of the field
// solve.
@@ -601,17 +600,27 @@ WarpX::PushParticlesandDepose (amrex::Real cur_time)
void
WarpX::PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type)
{
+ // If warpx.do_current_centering = 1, the current is deposited on the nodal MultiFab current_fp_nodal
+ // and then centered onto the staggered MultiFab current_fp
+ amrex::MultiFab* current_x = (WarpX::do_current_centering) ? current_fp_nodal[lev][0].get()
+ : current_fp[lev][0].get();
+ amrex::MultiFab* current_y = (WarpX::do_current_centering) ? current_fp_nodal[lev][1].get()
+ : current_fp[lev][1].get();
+ amrex::MultiFab* current_z = (WarpX::do_current_centering) ? current_fp_nodal[lev][2].get()
+ : current_fp[lev][2].get();
+
mypc->Evolve(lev,
*Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2],
*Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2],
*Efield_avg_aux[lev][0],*Efield_avg_aux[lev][1],*Efield_avg_aux[lev][2],
*Bfield_avg_aux[lev][0],*Bfield_avg_aux[lev][1],*Bfield_avg_aux[lev][2],
- *current_fp[lev][0],*current_fp[lev][1],*current_fp[lev][2],
+ *current_x, *current_y, *current_z,
current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(),
rho_fp[lev].get(), charge_buf[lev].get(),
Efield_cax[lev][0].get(), Efield_cax[lev][1].get(), Efield_cax[lev][2].get(),
Bfield_cax[lev][0].get(), Bfield_cax[lev][1].get(), Bfield_cax[lev][2].get(),
cur_time, dt[lev], a_dt_type);
+
#ifdef WARPX_DIM_RZ
// This is called after all particles have deposited their current and charge.
ApplyInverseVolumeScalingToCurrentDensity(current_fp[lev][0].get(), current_fp[lev][1].get(), current_fp[lev][2].get(), lev);
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index 1852d871f..f07995c84 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -316,6 +316,11 @@ WarpX::InitLevelData (int lev, Real /*time*/)
}
}
+ if (WarpX::do_current_centering)
+ {
+ current_fp_nodal[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 (fft_do_time_averaging) {
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp
index e5a912c8b..4131b7307 100644
--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -56,6 +56,9 @@ WarpX::UpdateAuxilaryDataStagToNodal ()
const amrex::IntVect& Ey_stag = Emf[0][1]->ixType().toIntVect();
const amrex::IntVect& Ez_stag = Emf[0][2]->ixType().toIntVect();
+ // Destination MultiFab (aux) always has nodal index type when this function is called
+ const amrex::IntVect& dst_stag = amrex::IntVect::TheNodeVector();
+
// For level 0, we only need to do the average.
#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
@@ -82,44 +85,49 @@ WarpX::UpdateAuxilaryDataStagToNodal ()
Box bx = mfi.fabbox();
if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+
#ifdef WARPX_USE_PSATD
- const int fg_nox = WarpX::field_gathering_nox;
- const int fg_noy = WarpX::field_gathering_noy;
- const int fg_noz = WarpX::field_gathering_noz;
- // Device vectors for Fornberg stencil coefficients used for finite-order centering
- amrex::Real const * stencil_coeffs_x = WarpX::device_centering_stencil_coeffs_x.data();
- amrex::Real const * stencil_coeffs_y = WarpX::device_centering_stencil_coeffs_y.data();
- amrex::Real const * stencil_coeffs_z = WarpX::device_centering_stencil_coeffs_z.data();
+
+ // Order of finite-order centering of fields
+ const int fg_nox = WarpX::field_centering_nox;
+ const int fg_noy = WarpX::field_centering_noy;
+ const int fg_noz = WarpX::field_centering_noz;
+
+ // Device vectors of stencil coefficients used for finite-order centering of fields
+ amrex::Real const * stencil_coeffs_x = WarpX::device_field_centering_stencil_coeffs_x.data();
+ amrex::Real const * stencil_coeffs_y = WarpX::device_field_centering_stencil_coeffs_y.data();
+ amrex::Real const * stencil_coeffs_z = WarpX::device_field_centering_stencil_coeffs_z.data();
+
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
{
- warpx_interp<true>(j, k, l, bx_aux, bx_fp, Bx_stag, fg_nox, fg_noy, fg_noz,
+ warpx_interp<true>(j, k, l, bx_aux, bx_fp, dst_stag, Bx_stag, fg_nox, fg_noy, fg_noz,
stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z);
- warpx_interp<true>(j, k, l, by_aux, by_fp, By_stag, fg_nox, fg_noy, fg_noz,
+ warpx_interp<true>(j, k, l, by_aux, by_fp, dst_stag, By_stag, fg_nox, fg_noy, fg_noz,
stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z);
- warpx_interp<true>(j, k, l, bz_aux, bz_fp, Bz_stag, fg_nox, fg_noy, fg_noz,
+ warpx_interp<true>(j, k, l, bz_aux, bz_fp, dst_stag, Bz_stag, fg_nox, fg_noy, fg_noz,
stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z);
- warpx_interp<true>(j, k, l, ex_aux, ex_fp, Ex_stag, fg_nox, fg_noy, fg_noz,
+ warpx_interp<true>(j, k, l, ex_aux, ex_fp, dst_stag, Ex_stag, fg_nox, fg_noy, fg_noz,
stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z);
- warpx_interp<true>(j, k, l, ey_aux, ey_fp, Ey_stag, fg_nox, fg_noy, fg_noz,
+ warpx_interp<true>(j, k, l, ey_aux, ey_fp, dst_stag, Ey_stag, fg_nox, fg_noy, fg_noz,
stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z);
- warpx_interp<true>(j, k, l, ez_aux, ez_fp, Ez_stag, fg_nox, fg_noy, fg_noz,
+ warpx_interp<true>(j, k, l, ez_aux, ez_fp, dst_stag, Ez_stag, fg_nox, fg_noy, fg_noz,
stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z);
});
#endif
} else { // FDTD
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
{
- warpx_interp(j, k, l, bx_aux, bx_fp, Bx_stag);
- warpx_interp(j, k, l, by_aux, by_fp, By_stag);
- warpx_interp(j, k, l, bz_aux, bz_fp, Bz_stag);
- warpx_interp(j, k, l, ex_aux, ex_fp, Ex_stag);
- warpx_interp(j, k, l, ey_aux, ey_fp, Ey_stag);
- warpx_interp(j, k, l, ez_aux, ez_fp, Ez_stag);
+ warpx_interp(j, k, l, bx_aux, bx_fp, dst_stag, Bx_stag);
+ warpx_interp(j, k, l, by_aux, by_fp, dst_stag, By_stag);
+ warpx_interp(j, k, l, bz_aux, bz_fp, dst_stag, Bz_stag);
+ warpx_interp(j, k, l, ex_aux, ex_fp, dst_stag, Ex_stag);
+ warpx_interp(j, k, l, ey_aux, ey_fp, dst_stag, Ey_stag);
+ warpx_interp(j, k, l, ez_aux, ez_fp, dst_stag, Ez_stag);
});
}
}
@@ -360,6 +368,67 @@ WarpX::UpdateAuxilaryDataSameType ()
}
}
+void WarpX::UpdateCurrentNodalToStag (amrex::MultiFab& dst, amrex::MultiFab const& src)
+{
+ // If source and destination MultiFabs have the same index type, a simple copy is enough
+ // (for example, this happens with the current along y in 2D, which is always fully nodal)
+ if (dst.ixType() == src.ixType())
+ {
+ amrex::MultiFab::Copy(dst, src, 0, 0, dst.nComp(), dst.nGrowVect());
+ return;
+ }
+
+ amrex::IntVect const& dst_stag = dst.ixType().toIntVect();
+
+ // Source MultiFab always has nodal index type when this function is called
+ amrex::IntVect const& src_stag = amrex::IntVect::TheNodeVector();
+
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+
+ for (MFIter mfi(dst); mfi.isValid(); ++mfi)
+ {
+ // Loop over full box including ghost cells
+ // (input arrays will be padded with zeros beyond ghost cells
+ // for out-of-bound accesses due to large-stencil operations)
+ Box bx = mfi.fabbox();
+
+ amrex::Array4<amrex::Real const> const& src_arr = src.const_array(mfi);
+ amrex::Array4<amrex::Real> const& dst_arr = dst.array(mfi);
+
+ if (maxwell_solver_id == MaxwellSolverAlgo::PSATD)
+ {
+#ifdef WARPX_USE_PSATD
+
+ // Order of finite-order centering of currents
+ const int cc_nox = WarpX::current_centering_nox;
+ const int cc_noy = WarpX::current_centering_noy;
+ const int cc_noz = WarpX::current_centering_noz;
+
+ // Device vectors of stencil coefficients used for finite-order centering of currents
+ amrex::Real const * stencil_coeffs_x = WarpX::device_current_centering_stencil_coeffs_x.data();
+ amrex::Real const * stencil_coeffs_y = WarpX::device_current_centering_stencil_coeffs_y.data();
+ amrex::Real const * stencil_coeffs_z = WarpX::device_current_centering_stencil_coeffs_z.data();
+
+ amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
+ {
+ warpx_interp<true>(j, k, l, dst_arr, src_arr, dst_stag, src_stag, cc_nox, cc_noy, cc_noz,
+ stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z);
+ });
+#endif
+ }
+
+ else // FDTD
+ {
+ amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
+ {
+ warpx_interp<false>(j, k, l, dst_arr, src_arr, dst_stag, src_stag);
+ });
+ }
+ }
+}
+
void
WarpX::FillBoundaryB (IntVect ng)
{
@@ -710,6 +779,17 @@ WarpX::SyncCurrent ()
{
WARPX_PROFILE("WarpX::SyncCurrent()");
+ // If warpx.do_current_centering = 1, center currents from nodal grid to staggered grid
+ if (WarpX::do_current_centering)
+ {
+ for (int lev = 0; lev <= finest_level; lev++)
+ {
+ WarpX::UpdateCurrentNodalToStag(*current_fp[lev][0], *current_fp_nodal[lev][0]);
+ WarpX::UpdateCurrentNodalToStag(*current_fp[lev][1], *current_fp_nodal[lev][1]);
+ WarpX::UpdateCurrentNodalToStag(*current_fp[lev][2], *current_fp_nodal[lev][2]);
+ }
+ }
+
// Restrict fine patch current onto the coarse patch, before
// summing the guard cells of the fine patch
for (int lev = 1; lev <= finest_level; ++lev)
diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H
index 64c5d0f7c..d1600ec7a 100644
--- a/Source/Parallelization/WarpXComm_K.H
+++ b/Source/Parallelization/WarpXComm_K.H
@@ -487,24 +487,24 @@ void warpx_interp_nd_efield_z (int j, int k, int l,
}
/**
- * \brief Arbitrary-order interpolation function used to interpolate a given field from the Yee grid
- * to the nodal grid, before gathering the field from the nodes to the particle positions
- * (momentum-conserving field gathering). With the FDTD solver, this performs simple linear interpolation.
+ * \brief Arbitrary-order interpolation function used to center a given MultiFab between two grids
+ * with different staggerings. With the FDTD solver, this performs simple linear interpolation.
* With the PSATD solver, this performs arbitrary-order interpolation based on the Fornberg coefficients.
- * The result is stored in the output array \c arr_aux, allocated on the \c aux patch.
+ * The result is stored in the output array \c dst_arr.
*
* \param[in] j index along x of the output array
* \param[in] k index along y (in 3D) or z (in 2D) of the output array
* \param[in] l index along z (in 3D, \c l = 0 in 2D) of the output array
- * \param[in,out] dst_arr output array allocated on the \c aux patch where interpolated values are stored
- * \param[in] src_arr input array allocated on the fine patch
- * \param[in] src_stag \c IndexType of the input array (Yee staggering)
- * \param[in] nox order of finite-order interpolation along x
- * \param[in] noy order of finite-order interpolation along y
- * \param[in] noz order of finite-order interpolation along z
- * \param[in] stencil_coeffs_x array of ordered Fornberg coefficients for finite-order interpolation stencil along x
- * \param[in] stencil_coeffs_y array of ordered Fornberg coefficients for finite-order interpolation stencil along y
- * \param[in] stencil_coeffs_z array of ordered Fornberg coefficients for finite-order interpolation stencil along z
+ * \param[in,out] dst_arr output array where interpolated values are stored
+ * \param[in] src_arr input array storing the values used for interpolation
+ * \param[in] dst_stag \c IndexType of the output array
+ * \param[in] src_stag \c IndexType of the input array
+ * \param[in] nox order of finite-order centering along x
+ * \param[in] noy order of finite-order centering along y
+ * \param[in] noz order of finite-order centering along z
+ * \param[in] stencil_coeffs_x array of ordered Fornberg coefficients for finite-order centering stencil along x
+ * \param[in] stencil_coeffs_y array of ordered Fornberg coefficients for finite-order centering stencil along y
+ * \param[in] stencil_coeffs_z array of ordered Fornberg coefficients for finite-order centering stencil along z
*/
template< bool IS_PSATD = false >
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
@@ -513,6 +513,7 @@ void warpx_interp (const int j,
const int l,
amrex::Array4<amrex::Real > const& dst_arr,
amrex::Array4<amrex::Real const> const& src_arr,
+ const amrex::IntVect& dst_stag,
const amrex::IntVect& src_stag,
const int nox = 2,
const int noy = 2,
@@ -540,11 +541,18 @@ void warpx_interp (const int j,
amrex::ignore_unused(stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z);
#endif
+ // If dst_nodal = true , we are centering from a staggered grid to a nodal grid
+ // If dst_nodal = false, we are centering from a nodal grid to a staggered grid
+ const bool dst_nodal = (dst_stag == amrex::IntVect::TheNodeVector());
+
+ // See 1D examples below to understand the meaning of this integer shift
+ const int shift = (dst_nodal) ? 0 : 1;
+
// Staggering (s = 0 if cell-centered, s = 1 if nodal)
- const int sj = src_stag[0];
- const int sk = src_stag[1];
+ const int sj = (dst_nodal) ? src_stag[0] : dst_stag[0];
+ const int sk = (dst_nodal) ? src_stag[1] : dst_stag[1];
#if (AMREX_SPACEDIM == 3)
- const int sl = src_stag[2];
+ const int sl = (dst_nodal) ? src_stag[2] : dst_stag[2];
#endif
// Interpolate along j,k,l only if source MultiFab is staggered along j,k,l
@@ -572,12 +580,12 @@ void warpx_interp (const int j,
#endif
// Min and max for interpolation loop along j
- const int jmin = (interp_j) ? j - noj/2 : j;
- const int jmax = (interp_j) ? j + noj/2 - 1 : j;
+ const int jmin = (interp_j) ? j - noj/2 + shift : j;
+ const int jmax = (interp_j) ? j + noj/2 + shift - 1 : j;
// Min and max for interpolation loop along k
- const int kmin = (interp_k) ? k - nok/2 : k;
- const int kmax = (interp_k) ? k + nok/2 - 1 : k;
+ const int kmin = (interp_k) ? k - nok/2 + shift : k;
+ const int kmax = (interp_k) ? k + nok/2 + shift - 1 : k;
// Min and max for interpolation loop along l
#if (AMREX_SPACEDIM == 2)
@@ -585,44 +593,73 @@ void warpx_interp (const int j,
const int lmin = l;
const int lmax = l;
#elif (AMREX_SPACEDIM == 3)
- const int lmin = (interp_l) ? l - nol/2 : l;
- const int lmax = (interp_l) ? l + nol/2 - 1 : l;
+ const int lmin = (interp_l) ? l - nol/2 + shift : l;
+ const int lmax = (interp_l) ? l + nol/2 + shift - 1 : l;
#endif
+ // Number of interpolation points
+ const int nj = jmax - jmin;
+ const int nk = kmax - kmin;
+ const int nl = lmax - lmin;
+
+ // Example of 1D centering from nodal grid to nodal grid (simple copy):
+ //
+ // j
+ // --o-----o-----o-- result(j) = f(j)
+ // --o-----o-----o--
+ // j-1 j j+1
+ //
+ // Example of 1D linear centering from staggered grid to nodal grid:
+ //
+ // j
+ // --o-----o-----o-- result(j) = (f(j-1) + f(j)) / 2
+ // -----x-----x-----
+ // j-1 j
+ //
+ // Example of 1D linear centering from nodal grid to staggered grid:
+ // (note the shift of +1 in the indices with respect to the case above, see variable "shift")
+ //
+ // j
+ // --x-----x-----x-- result(j) = (f(j) + f(j+1)) / 2
+ // -----o-----o-----
+ // j j+1
+ //
+ // Example of 1D finite-order centering from staggered grid to nodal grid:
+ //
+ // j
+ // --o-----o-----o-----o-----o-----o-----o-- result(j) = c_0 * (f(j-1) + f(j) ) / 2
+ // -----x-----x-----x-----x-----x-----x----- + c_1 * (f(j-2) + f(j+1)) / 2
+ // j-3 j-2 j-1 j j+1 j+2 + c_2 * (f(j-3) + f(j+2)) / 2
+ // c_2 c_1 c_0 c_0 c_1 c_2 + ...
+ //
+ // Example of 1D finite-order centering from nodal grid to staggered grid:
+ // (note the shift of +1 in the indices with respect to the case above, see variable "shift")
+ //
+ // j
+ // --x-----x-----x-----x-----x-----x-----x-- result(j) = c_0 * (f(j) + f(j+1)) / 2
+ // -----o-----o-----o-----o-----o-----o----- + c_1 * (f(j-1) + f(j+2)) / 2
+ // j-2 j-1 j j+1 j+2 j+3 + c_2 * (f(j-2) + f(j+3)) / 2
+ // c_2 c_1 c_0 c_0 c_1 c_2 + ...
+
amrex::Real res = 0.0_rt;
- // FDTD (linear interpolation)
- if( !IS_PSATD ) {
- // Example of 1D linear interpolation from nodal grid to nodal grid:
- //
- // j
- // --o-----o-----o-- result(j) = f(j)
- // --o-----o-----o--
- // j-1 j j+1
- //
- // Example of 1D linear interpolation from staggered grid to nodal grid:
- //
- // j
- // --o-----o-----o-- result(j) = (f(j-1) + f(j)) / 2
- // -----x-----x-----
- // j-1 j
-
- for (int ll = lmin; ll <= lmax; ll++)
+ if( !IS_PSATD ) // FDTD (linear interpolation)
+ {
+ for (int ll = 0; ll <= nl; ll++)
{
- for (int kk = kmin; kk <= kmax; kk++)
+ for (int kk = 0; kk <= nk; kk++)
{
- for (int jj = jmin; jj <= jmax; jj++)
+ for (int jj = 0; jj <= nj; jj++)
{
- res += src_arr_zeropad(jj,kk,ll);
+ res += src_arr_zeropad(jmin+jj,kmin+kk,lmin+ll);
}
}
}
- // PSATD (finite-order interpolation)
- } else {
- const int nj = jmax - jmin;
- const int nk = kmax - kmin;
- const int nl = lmax - lmin;
+ }
+
+ else // PSATD (finite-order interpolation)
+ {
amrex::Real cj = 1.0_rt;
amrex::Real ck = 1.0_rt;
amrex::Real cl = 1.0_rt;
@@ -635,21 +672,6 @@ void warpx_interp (const int j,
amrex::Real const* scl = stencil_coeffs_z;
#endif
- // Example of 1D finite-order interpolation from nodal grid to nodal grid:
- //
- // j
- // --o-----o-----o-- result(j) = f(j)
- // --o-----o-----o--
- // j-1 j j+1
- //
- // Example of 1D finite-order interpolation from staggered grid to nodal grid:
- //
- // j
- // --o-----o-----o-----o-----o-----o-----o-- result(j) = c_0 * (f(j-1) + f(j) ) / 2
- // -----x-----x-----x-----x-----x-----x----- + c_1 * (f(j-2) + f(j+1)) / 2
- // j-3 j-2 j-1 j j+1 j+2 + c_2 * (f(j-3) + f(j+2)) / 2
- // c_2 c_1 c_0 c_0 c_1 c_2 + ...
-
for (int ll = 0; ll <= nl; ll++)
{
#if (AMREX_SPACEDIM == 3)
@@ -667,7 +689,7 @@ void warpx_interp (const int j,
}
}
}
- } // PSATD (finite-order interpolation)
+ }
dst_arr(j,k,l) = wj * wk * wl * res;
}
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 2c4b9c20e..4ba4c9d43 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -142,6 +142,9 @@ public:
static amrex::Vector<int> particle_boundary_hi;
+ // If true, the current is deposited on a nodal grid and then centered onto a staggered grid
+ static bool do_current_centering;
+
// PSATD: If true (overwritten by the user in the input file), the current correction
// defined in equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010 is applied
bool current_correction = false;
@@ -158,11 +161,15 @@ public:
static long noy;
static long noz;
- // For momentum-conserving field gathering, order of interpolation from the
- // staggered positions to the grid nodes
- static int field_gathering_nox;
- static int field_gathering_noy;
- static int field_gathering_noz;
+ // Order of finite-order centering of fields (staggered to nodal)
+ static int field_centering_nox;
+ static int field_centering_noy;
+ static int field_centering_noz;
+
+ // Order of finite-order centering of currents (nodal to staggered)
+ static int current_centering_nox;
+ static int current_centering_noy;
+ static int current_centering_noz;
// Number of modes for the RZ multimode version
static int n_rz_azimuthal_modes;
@@ -441,6 +448,16 @@ public:
void UpdateAuxilaryDataStagToNodal ();
void UpdateAuxilaryDataSameType ();
+ /**
+ * \brief This function is called if \c warpx.do_current_centering = 1 and
+ * it centers the currents from a nodal grid to a staggered grid (Yee) using
+ * finite-order interpolation based on the Fornberg coefficients.
+ *
+ * \param[in,out] dst destination \c MultiFab where the results of the finite-order centering are stored
+ * \param[in] src source \c MultiFab that contains the values of the nodal current to be centered
+ */
+ void UpdateCurrentNodalToStag (amrex::MultiFab& dst, amrex::MultiFab const& src);
+
// Fill boundary cells including coarse/fine boundaries
void FillBoundaryB (amrex::IntVect ng);
void FillBoundaryE (amrex::IntVect ng);
@@ -605,13 +622,14 @@ public:
void ApplyFilterandSumBoundaryRho (int lev, int glev, amrex::MultiFab& rho, int icomp, int ncomp);
#ifdef WARPX_USE_PSATD
- // Host and device vectors for ordered Fornberg stencil coefficients used for finite-order centering
- amrex::Vector<amrex::Real> host_centering_stencil_coeffs_x;
- amrex::Vector<amrex::Real> host_centering_stencil_coeffs_y;
- amrex::Vector<amrex::Real> host_centering_stencil_coeffs_z;
- amrex::Gpu::DeviceVector<amrex::Real> device_centering_stencil_coeffs_x;
- amrex::Gpu::DeviceVector<amrex::Real> device_centering_stencil_coeffs_y;
- amrex::Gpu::DeviceVector<amrex::Real> device_centering_stencil_coeffs_z;
+ // Device vectors of stencil coefficients used for finite-order centering of fields
+ amrex::Gpu::DeviceVector<amrex::Real> device_field_centering_stencil_coeffs_x;
+ amrex::Gpu::DeviceVector<amrex::Real> device_field_centering_stencil_coeffs_y;
+ amrex::Gpu::DeviceVector<amrex::Real> device_field_centering_stencil_coeffs_z;
+ // Device vectors of stencil coefficients used for finite-order centering of currents
+ amrex::Gpu::DeviceVector<amrex::Real> device_current_centering_stencil_coeffs_x;
+ amrex::Gpu::DeviceVector<amrex::Real> device_current_centering_stencil_coeffs_y;
+ amrex::Gpu::DeviceVector<amrex::Real> device_current_centering_stencil_coeffs_z;
#endif
protected:
@@ -760,9 +778,37 @@ private:
return gather_buffer_masks[lev].get();
}
- void ReorderFornbergCoefficients(amrex::Vector<amrex::Real>& ordered_coeffs,
- amrex::Vector<amrex::Real>& unordered_coeffs,
- const int order);
+#ifdef WARPX_USE_PSATD
+ /**
+ * \brief Re-orders the Fornberg coefficients so that they can be used more conveniently for
+ * finite-order centering operations. For example, for finite-order centering of order 6,
+ * the Fornberg coefficients \c (c_0,c_1,c_2) are re-ordered as \c (c_2,c_1,c_0,c_0,c_1,c_2).
+ *
+ * \param[in,out] ordered_coeffs host vector where the re-ordered Fornberg coefficients will be stored
+ * \param[in] unordered_coeffs host vector storing the original sequence of Fornberg coefficients
+ * \param[in] order order of the finite-order centering along a given direction
+ */
+ void ReorderFornbergCoefficients (amrex::Vector<amrex::Real>& ordered_coeffs,
+ amrex::Vector<amrex::Real>& unordered_coeffs,
+ const int order);
+ /**
+ * \brief Allocates and initializes the stencil coefficients used for the finite-order centering
+ * of fields and currents, and stores them in the given device vectors.
+ *
+ * \param[in,out] device_centering_stencil_coeffs_x device vector where the stencil coefficients along x will be stored
+ * \param[in,out] device_centering_stencil_coeffs_y device vector where the stencil coefficients along y will be stored
+ * \param[in,out] device_centering_stencil_coeffs_z device vector where the stencil coefficients along z will be stored
+ * \param[in] centering_nox order of the finite-order centering along x
+ * \param[in] centering_noy order of the finite-order centering along y
+ * \param[in] centering_noz order of the finite-order centering along z
+ */
+ void AllocateCenteringCoefficients (amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_x,
+ amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_y,
+ amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_z,
+ const int centering_nox,
+ const int centering_noy,
+ const int centering_noz);
+#endif
void AllocLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm,
const amrex::IntVect& ngE, const amrex::IntVect& ngJ,
@@ -805,6 +851,9 @@ private:
// store fine patch
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_store;
+ // Nodal MultiFab for nodal current deposition if warpx.do_current_centering = 1
+ amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>> current_fp_nodal;
+
// Coarse patch
amrex::Vector< std::unique_ptr<amrex::MultiFab> > F_cp;
amrex::Vector< std::unique_ptr<amrex::MultiFab> > rho_cp;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 92a2efb61..7417921ab 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -83,6 +83,8 @@ amrex::Vector<int> WarpX::field_boundary_hi(AMREX_SPACEDIM,0);
amrex::Vector<int> WarpX::particle_boundary_lo(AMREX_SPACEDIM,0);
amrex::Vector<int> WarpX::particle_boundary_hi(AMREX_SPACEDIM,0);
+bool WarpX::do_current_centering = false;
+
int WarpX::n_rz_azimuthal_modes = 1;
int WarpX::ncomps = 1;
@@ -90,11 +92,15 @@ long WarpX::nox = 1;
long WarpX::noy = 1;
long WarpX::noz = 1;
-// For momentum-conserving field gathering, order of interpolation from the
-// staggered positions to the grid nodes
-int WarpX::field_gathering_nox = 2;
-int WarpX::field_gathering_noy = 2;
-int WarpX::field_gathering_noz = 2;
+// Order of finite-order centering of fields (staggered to nodal)
+int WarpX::field_centering_nox = 2;
+int WarpX::field_centering_noy = 2;
+int WarpX::field_centering_noz = 2;
+
+// Order of finite-order centering of currents (nodal to staggered)
+int WarpX::current_centering_nox = 2;
+int WarpX::current_centering_noy = 2;
+int WarpX::current_centering_noz = 2;
bool WarpX::use_fdtd_nci_corr = false;
bool WarpX::galerkin_interpolation = true;
@@ -229,6 +235,11 @@ WarpX::WarpX ()
current_store.resize(nlevs_max);
+ if (do_current_centering)
+ {
+ current_fp_nodal.resize(nlevs_max);
+ }
+
F_cp.resize(nlevs_max);
rho_cp.resize(nlevs_max);
current_cp.resize(nlevs_max);
@@ -630,6 +641,21 @@ WarpX::ReadParameters ()
// Only needs to be set with WARPX_DIM_RZ, otherwise defaults to 1
pp_warpx.query("n_rz_azimuthal_modes", n_rz_azimuthal_modes);
+
+ // If true, the current is deposited on a nodal grid and then interpolated onto a Yee grid
+ pp_warpx.query("do_current_centering", do_current_centering);
+
+ // If do_nodal = 1, Maxwell's equations are solved on a nodal grid and
+ // the current should not be centered
+ if (do_nodal)
+ {
+ do_current_centering = false;
+ }
+
+ if ((maxLevel() > 0) && do_current_centering)
+ {
+ amrex::Abort("\nFinite-order centering of currents is not implemented with mesh refinement");
+ }
}
{
@@ -695,64 +721,50 @@ WarpX::ReadParameters ()
pp_interpolation.query("noz", noz);
#ifdef WARPX_USE_PSATD
+
if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+
// For momentum-conserving field gathering, read from input the order of
// interpolation from the staggered positions to the grid nodes
- if (field_gathering_algo == GatheringAlgo::MomentumConserving) {
- pp_interpolation.query("field_gathering_nox", field_gathering_nox);
- pp_interpolation.query("field_gathering_noy", field_gathering_noy);
- pp_interpolation.query("field_gathering_noz", field_gathering_noz);
+ if (WarpX::field_gathering_algo == GatheringAlgo::MomentumConserving) {
+ pp_interpolation.query("field_centering_nox", field_centering_nox);
+ pp_interpolation.query("field_centering_noy", field_centering_noy);
+ pp_interpolation.query("field_centering_noz", field_centering_noz);
+ }
+
+ // Read order of finite-order centering of currents (nodal to staggered)
+ if (WarpX::do_current_centering) {
+ pp_interpolation.query("current_centering_nox", current_centering_nox);
+ pp_interpolation.query("current_centering_noy", current_centering_noy);
+ pp_interpolation.query("current_centering_noz", current_centering_noz);
}
- if (maxLevel() > 0) {
+ if (maxLevel() > 0)
+ {
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
- field_gathering_nox == 2 && field_gathering_noy == 2 && field_gathering_noz == 2,
- "High-order interpolation (order > 2) is not implemented with mesh refinement");
+ field_centering_nox == 2 && field_centering_noy == 2 && field_centering_noz == 2,
+ "High-order centering of fields (order > 2) is not implemented with mesh refinement");
}
- // Host vectors for Fornberg stencil coefficients
- amrex::Vector<amrex::Real> host_Fornberg_stencil_coeffs_x;
- amrex::Vector<amrex::Real> host_Fornberg_stencil_coeffs_y;
- amrex::Vector<amrex::Real> host_Fornberg_stencil_coeffs_z;
-
- host_Fornberg_stencil_coeffs_x = getFornbergStencilCoefficients(field_gathering_nox, false);
- host_Fornberg_stencil_coeffs_y = getFornbergStencilCoefficients(field_gathering_noy, false);
- host_Fornberg_stencil_coeffs_z = getFornbergStencilCoefficients(field_gathering_noz, false);
-
- // Host vectors for ordered Fornberg stencil coefficients used for finite-order centering
- host_centering_stencil_coeffs_x.resize(field_gathering_nox);
- host_centering_stencil_coeffs_y.resize(field_gathering_noy);
- host_centering_stencil_coeffs_z.resize(field_gathering_noz);
-
- // Re-order Fornberg stencil coefficients
- // example for order 6: (c_0,c_1,c_2) becomes (c_2,c_1,c_0,c_0,c_1,c_2)
- ReorderFornbergCoefficients(host_centering_stencil_coeffs_x,
- host_Fornberg_stencil_coeffs_x, field_gathering_nox);
- ReorderFornbergCoefficients(host_centering_stencil_coeffs_y,
- host_Fornberg_stencil_coeffs_y, field_gathering_noy);
- ReorderFornbergCoefficients(host_centering_stencil_coeffs_z,
- host_Fornberg_stencil_coeffs_z, field_gathering_noz);
-
- // Device vectors for ordered Fornberg stencil coefficients used for finite-order centering
- device_centering_stencil_coeffs_x.resize(host_centering_stencil_coeffs_x.size());
- amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
- host_centering_stencil_coeffs_x.begin(),
- host_centering_stencil_coeffs_x.end(),
- device_centering_stencil_coeffs_x.begin());
-
- device_centering_stencil_coeffs_y.resize(host_centering_stencil_coeffs_y.size());
- amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
- host_centering_stencil_coeffs_y.begin(),
- host_centering_stencil_coeffs_y.end(),
- device_centering_stencil_coeffs_y.begin());
-
- device_centering_stencil_coeffs_z.resize(host_centering_stencil_coeffs_z.size());
- amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
- host_centering_stencil_coeffs_z.begin(),
- host_centering_stencil_coeffs_z.end(),
- device_centering_stencil_coeffs_z.begin());
-
- amrex::Gpu::synchronize();
+ if (WarpX::field_gathering_algo == GatheringAlgo::MomentumConserving)
+ {
+ AllocateCenteringCoefficients(device_field_centering_stencil_coeffs_x,
+ device_field_centering_stencil_coeffs_y,
+ device_field_centering_stencil_coeffs_z,
+ field_centering_nox,
+ field_centering_noy,
+ field_centering_noz);
+ }
+
+ if (WarpX::do_current_centering)
+ {
+ AllocateCenteringCoefficients(device_current_centering_stencil_coeffs_x,
+ device_current_centering_stencil_coeffs_y,
+ device_current_centering_stencil_coeffs_z,
+ current_centering_nox,
+ current_centering_noy,
+ current_centering_noz);
+ }
}
#endif
@@ -1031,6 +1043,11 @@ WarpX::ClearLevel (int lev)
current_store[lev][i].reset();
+ if (do_current_centering)
+ {
+ current_fp_nodal[lev][i].reset();
+ }
+
current_cp[lev][i].reset();
Efield_cp [lev][i].reset();
Bfield_cp [lev][i].reset();
@@ -1209,6 +1226,14 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
current_fp[lev][1] = std::make_unique<MultiFab>(amrex::convert(ba,jy_nodal_flag),dm,ncomps,ngJ,tag("current_fp[y]"));
current_fp[lev][2] = std::make_unique<MultiFab>(amrex::convert(ba,jz_nodal_flag),dm,ncomps,ngJ,tag("current_fp[z]"));
+ if (do_current_centering)
+ {
+ amrex::BoxArray const& nodal_ba = amrex::convert(ba, amrex::IntVect::TheNodeVector());
+ current_fp_nodal[lev][0] = std::make_unique<MultiFab>(nodal_ba, dm, ncomps, ngJ);
+ current_fp_nodal[lev][1] = std::make_unique<MultiFab>(nodal_ba, dm, ncomps, ngJ);
+ current_fp_nodal[lev][2] = std::make_unique<MultiFab>(nodal_ba, dm, ncomps, ngJ);
+ }
+
Bfield_avg_fp[lev][0] = std::make_unique<MultiFab>(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE,tag("Bfield_avg_fp[x]"));
Bfield_avg_fp[lev][1] = std::make_unique<MultiFab>(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE,tag("Bfield_avg_fp[y]"));
Bfield_avg_fp[lev][2] = std::make_unique<MultiFab>(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE,tag("Bfield_avg_fp[z]"));
@@ -1768,9 +1793,10 @@ WarpX::BuildBufferMasksInBox ( const amrex::Box tbx, amrex::IArrayBox &buffer_ma
#endif
}
-void WarpX::ReorderFornbergCoefficients(amrex::Vector<amrex::Real>& ordered_coeffs,
- amrex::Vector<amrex::Real>& unordered_coeffs,
- const int order)
+#ifdef WARPX_USE_PSATD
+void WarpX::ReorderFornbergCoefficients (amrex::Vector<amrex::Real>& ordered_coeffs,
+ amrex::Vector<amrex::Real>& unordered_coeffs,
+ const int order)
{
const int n = order / 2;
for (int i = 0; i < n; i++) {
@@ -1781,6 +1807,64 @@ void WarpX::ReorderFornbergCoefficients(amrex::Vector<amrex::Real>& ordered_coef
}
}
+void WarpX::AllocateCenteringCoefficients (amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_x,
+ amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_y,
+ amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_z,
+ const int centering_nox,
+ const int centering_noy,
+ const int centering_noz)
+{
+ // Vectors of Fornberg stencil coefficients
+ amrex::Vector<amrex::Real> Fornberg_stencil_coeffs_x;
+ amrex::Vector<amrex::Real> Fornberg_stencil_coeffs_y;
+ amrex::Vector<amrex::Real> Fornberg_stencil_coeffs_z;
+
+ // Host vectors of stencil coefficients used for finite-order centering
+ amrex::Vector<amrex::Real> host_centering_stencil_coeffs_x;
+ amrex::Vector<amrex::Real> host_centering_stencil_coeffs_y;
+ amrex::Vector<amrex::Real> host_centering_stencil_coeffs_z;
+
+ Fornberg_stencil_coeffs_x = getFornbergStencilCoefficients(centering_nox, false);
+ Fornberg_stencil_coeffs_y = getFornbergStencilCoefficients(centering_noy, false);
+ Fornberg_stencil_coeffs_z = getFornbergStencilCoefficients(centering_noz, false);
+
+ host_centering_stencil_coeffs_x.resize(centering_nox);
+ host_centering_stencil_coeffs_y.resize(centering_noy);
+ host_centering_stencil_coeffs_z.resize(centering_noz);
+
+ // Re-order Fornberg stencil coefficients:
+ // example for order 6: (c_0,c_1,c_2) becomes (c_2,c_1,c_0,c_0,c_1,c_2)
+ ReorderFornbergCoefficients(host_centering_stencil_coeffs_x,
+ Fornberg_stencil_coeffs_x, centering_nox);
+ ReorderFornbergCoefficients(host_centering_stencil_coeffs_y,
+ Fornberg_stencil_coeffs_y, centering_noy);
+ ReorderFornbergCoefficients(host_centering_stencil_coeffs_z,
+ Fornberg_stencil_coeffs_z, centering_noz);
+
+ // Device vectors of stencil coefficients used for finite-order centering
+
+ device_centering_stencil_coeffs_x.resize(host_centering_stencil_coeffs_x.size());
+ amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+ host_centering_stencil_coeffs_x.begin(),
+ host_centering_stencil_coeffs_x.end(),
+ device_centering_stencil_coeffs_x.begin());
+
+ device_centering_stencil_coeffs_y.resize(host_centering_stencil_coeffs_y.size());
+ amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+ host_centering_stencil_coeffs_y.begin(),
+ host_centering_stencil_coeffs_y.end(),
+ device_centering_stencil_coeffs_y.begin());
+
+ device_centering_stencil_coeffs_z.resize(host_centering_stencil_coeffs_z.size());
+ amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+ host_centering_stencil_coeffs_z.begin(),
+ host_centering_stencil_coeffs_z.end(),
+ device_centering_stencil_coeffs_z.begin());
+
+ amrex::Gpu::synchronize();
+}
+#endif
+
const iMultiFab*
WarpX::CurrentBufferMasks (int lev)
{
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-19 17:27:43 +0000