aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/BoundaryConditions/PML.cpp7
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt1
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H92
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp502
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H6
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp258
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package1
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralSolver.H1
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralSolver.cpp20
-rw-r--r--Source/Parallelization/GuardCellManager.H1
-rw-r--r--Source/Parallelization/GuardCellManager.cpp8
-rw-r--r--Source/Parallelization/WarpXComm.cpp60
-rw-r--r--Source/Parallelization/WarpXComm_K.H202
-rw-r--r--Source/WarpX.H8
-rw-r--r--Source/WarpX.cpp48
15 files changed, 1120 insertions, 95 deletions
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp
index 528dcb6cd..58ab7b9a4 100644
--- a/Source/BoundaryConditions/PML.cpp
+++ b/Source/BoundaryConditions/PML.cpp
@@ -546,10 +546,11 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/,
const RealVect dx{AMREX_D_DECL(geom->CellSize(0), geom->CellSize(1), geom->CellSize(2))};
// Get the cell-centered box, with guard cells
BoxArray realspace_ba = ba; // Copy box
- Array<Real,3> v_galilean_zero = {0,0,0};
+ Array<Real,3> v_galilean_zero = {0., 0., 0.};
+ Array<Real,3> v_comoving_zero = {0., 0., 0.};
realspace_ba.enclosedCells().grow(nge); // cell-centered + guard cells
spectral_solver_fp = std::make_unique<SpectralSolver>(realspace_ba, dm,
- nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, dx, dt, in_pml );
+ nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, dx, dt, in_pml );
#endif
if (cgeom)
@@ -619,7 +620,7 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/,
realspace_cba.enclosedCells().grow(nge); // cell-centered + guard cells
spectral_solver_cp = std::make_unique<SpectralSolver>(realspace_cba, cdm,
- nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, cdx, dt, in_pml );
+ nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, cdx, dt, in_pml );
#endif
}
}
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
index a050bd5c2..ab8e439ed 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
@@ -5,6 +5,7 @@ target_sources(WarpX
PsatdAlgorithm.cpp
SpectralBaseAlgorithm.cpp
AvgGalileanAlgorithm.cpp
+ ComovingPsatdAlgorithm.cpp
)
if(WarpX_DIMS STREQUAL RZ)
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H
new file mode 100644
index 000000000..2d9962e4d
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H
@@ -0,0 +1,92 @@
+#ifndef WARPX_COMOVING_PSATD_ALGORITHM_H_
+#define WARPX_COMOVING_PSATD_ALGORITHM_H_
+
+#include "SpectralBaseAlgorithm.H"
+
+#if WARPX_USE_PSATD
+
+/* \brief Class that updates the field in spectral space and stores the coefficients
+ * of the corresponding update equation, according to the comoving spectral scheme.
+ */
+class ComovingPsatdAlgorithm : public SpectralBaseAlgorithm
+{
+ public:
+
+ /**
+ * \brief Class constructor
+ */
+ ComovingPsatdAlgorithm (const SpectralKSpace& spectral_kspace,
+ const amrex::DistributionMapping& dm,
+ const int norder_x,
+ const int norder_y,
+ const int norder_z,
+ const bool nodal,
+ const amrex::Array<amrex::Real,3>& v_comoving,
+ const amrex::Real dt,
+ const bool update_with_rho);
+
+ /**
+ * \brief Override the update equations in Fourier space
+ */
+ virtual void pushSpectralFields (SpectralFieldData& f) const override final;
+
+ virtual int getRequiredNumberOfFields () const override final
+ {
+ return SpectralFieldIndex::n_fields;
+ };
+
+ /* \brief Initialize the coefficients needed in the update equations
+ */
+ void InitializeSpectralCoefficients (const SpectralKSpace& spectral_kspace,
+ const amrex::DistributionMapping& dm,
+ const amrex::Real dt);
+
+ /**
+ * \brief Virtual function for current correction in Fourier space.
+ * This function overrides the virtual function \c CurrentCorrection in the
+ * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+ * derived classes.
+ *
+ * \param[in,out] field_data All fields in Fourier space
+ * \param[in,out] current Array of unique pointers to \c MultiFab storing
+ * the three components of the current density
+ * \param[in] rho Unique pointer to \c MultiFab storing the charge density
+ */
+ virtual void CurrentCorrection (SpectralFieldData& field_data,
+ std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+ const std::unique_ptr<amrex::MultiFab>& rho) override final;
+
+ /**
+ * \brief Virtual function for Vay current deposition in Fourier space.
+ * This function overrides the virtual function \c VayDeposition in the
+ * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+ * derived classes.
+ *
+ * \param[in,out] field_data All fields in Fourier space
+ * \param[in,out] current Array of unique pointers to \c MultiFab storing
+ * the three components of the current density
+ */
+ virtual void VayDeposition (SpectralFieldData& field_data,
+ std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
+
+ private:
+
+ // Real and complex spectral coefficients
+ SpectralRealCoefficients C_coef, S_ck_coef;
+ SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef;
+
+ // k vectors
+ KVectorComponent kx_vec;
+#if (AMREX_SPACEDIM==3)
+ KVectorComponent ky_vec;
+#endif
+ KVectorComponent kz_vec;
+
+ // Additional member variables
+ amrex::Array<amrex::Real,3> m_v_comoving;
+ amrex::Real m_dt;
+ bool m_update_with_rho;
+};
+
+#endif // WARPX_USE_PSATD
+#endif // WARPX_COMOVING_PSATD_ALGORITHM_H_
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp
new file mode 100644
index 000000000..ca82d5a8c
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp
@@ -0,0 +1,502 @@
+#include "ComovingPsatdAlgorithm.H"
+#include "Utils/WarpXConst.H"
+
+#if WARPX_USE_PSATD
+
+using namespace amrex;
+
+ComovingPsatdAlgorithm::ComovingPsatdAlgorithm (const SpectralKSpace& spectral_kspace,
+ const DistributionMapping& dm,
+ const int norder_x, const int norder_y,
+ const int norder_z, const bool nodal,
+ const amrex::Array<amrex::Real, 3>& v_comoving,
+ const amrex::Real dt,
+ const bool update_with_rho)
+ // Members initialization
+ : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal),
+ // Initialize the infinite-order k vectors (the argument n_order = -1 selects
+ // the infinite order option, the argument nodal = false is then irrelevant)
+ kx_vec(spectral_kspace.getModifiedKComponent(dm, 0, -1, false)),
+#if (AMREX_SPACEDIM==3)
+ ky_vec(spectral_kspace.getModifiedKComponent(dm, 1, -1, false)),
+ kz_vec(spectral_kspace.getModifiedKComponent(dm, 2, -1, false)),
+#else
+ kz_vec(spectral_kspace.getModifiedKComponent(dm, 1, -1, false)),
+#endif
+ m_v_comoving(v_comoving),
+ m_dt(dt),
+ m_update_with_rho(update_with_rho)
+{
+ const BoxArray& ba = spectral_kspace.spectralspace_ba;
+
+ // Allocate arrays of real spectral coefficients
+ C_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+ S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+
+ // Allocate arrays of complex spectral coefficients
+ X1_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+ X2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+ X3_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+ X4_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+ Theta2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+
+ // Initialize real and complex spectral coefficients
+ InitializeSpectralCoefficients(spectral_kspace, dm, dt);
+}
+
+void
+ComovingPsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
+{
+ // Loop over boxes
+ for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi){
+
+ const amrex::Box& bx = f.fields[mfi].box();
+
+ // Extract arrays for the fields to be updated
+ amrex::Array4<Complex> fields = f.fields[mfi].array();
+
+ // Extract arrays for the coefficients
+ amrex::Array4<const amrex::Real> C_arr = C_coef [mfi].array();
+ amrex::Array4<const amrex::Real> S_ck_arr = S_ck_coef[mfi].array();
+ amrex::Array4<const Complex> X1_arr = X1_coef [mfi].array();
+ amrex::Array4<const Complex> X2_arr = X2_coef [mfi].array();
+ amrex::Array4<const Complex> X3_arr = X3_coef [mfi].array();
+ amrex::Array4<const Complex> X4_arr = X4_coef [mfi].array();
+
+ // Extract pointers for the k vectors
+ const amrex::Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr();
+#if (AMREX_SPACEDIM==3)
+ const amrex::Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr();
+#endif
+ const amrex::Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr();
+
+ // Loop over indices within one box
+ amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+ {
+ // Record old values of the fields to be updated
+ using Idx = SpectralFieldIndex;
+ const Complex Ex_old = fields(i,j,k,Idx::Ex);
+ const Complex Ey_old = fields(i,j,k,Idx::Ey);
+ const Complex Ez_old = fields(i,j,k,Idx::Ez);
+ const Complex Bx_old = fields(i,j,k,Idx::Bx);
+ const Complex By_old = fields(i,j,k,Idx::By);
+ const Complex Bz_old = fields(i,j,k,Idx::Bz);
+
+ // Shortcuts for the values of J and rho
+ const Complex Jx = fields(i,j,k,Idx::Jx);
+ const Complex Jy = fields(i,j,k,Idx::Jy);
+ const Complex Jz = fields(i,j,k,Idx::Jz);
+ const Complex rho_old = fields(i,j,k,Idx::rho_old);
+ const Complex rho_new = fields(i,j,k,Idx::rho_new);
+
+ // k vector values
+ const amrex::Real kx_mod = modified_kx_arr[i];
+#if (AMREX_SPACEDIM==3)
+ const amrex::Real ky_mod = modified_ky_arr[j];
+ const amrex::Real kz_mod = modified_kz_arr[k];
+#else
+ constexpr amrex::Real ky_mod = 0._rt;
+ const amrex::Real kz_mod = modified_kz_arr[j];
+#endif
+
+ // Physical constant c**2 and imaginary unit
+ constexpr amrex::Real c2 = PhysConst::c*PhysConst::c;
+ constexpr Complex I = Complex{0._rt,1._rt};
+
+ // The definition of these coefficients is explained in more detail
+ // in the function InitializeSpectralCoefficients below
+ const amrex::Real C = C_arr(i,j,k);
+ const amrex::Real S_ck = S_ck_arr(i,j,k);
+ const Complex X1 = X1_arr(i,j,k);
+ const Complex X2 = X2_arr(i,j,k);
+ const Complex X3 = X3_arr(i,j,k);
+ const Complex X4 = X4_arr(i,j,k);
+
+ // Update E
+ fields(i,j,k,Idx::Ex) = C*Ex_old + S_ck*c2*I*(ky_mod*Bz_old - kz_mod*By_old)
+ + X4*Jx - I*(X2*rho_new - X3*rho_old)*kx_mod;
+
+ fields(i,j,k,Idx::Ey) = C*Ey_old + S_ck*c2*I*(kz_mod*Bx_old - kx_mod*Bz_old)
+ + X4*Jy - I*(X2*rho_new - X3*rho_old)*ky_mod;
+
+ fields(i,j,k,Idx::Ez) = C*Ez_old + S_ck*c2*I*(kx_mod*By_old - ky_mod*Bx_old)
+ + X4*Jz - I*(X2*rho_new - X3*rho_old)*kz_mod;
+
+ // Update B
+ fields(i,j,k,Idx::Bx) = C*Bx_old - S_ck*I*(ky_mod*Ez_old - kz_mod*Ey_old)
+ + X1*I*(ky_mod*Jz - kz_mod*Jy);
+
+ fields(i,j,k,Idx::By) = C*By_old - S_ck*I*(kz_mod*Ex_old - kx_mod*Ez_old)
+ + X1*I*(kz_mod*Jx - kx_mod*Jz);
+
+ fields(i,j,k,Idx::Bz) = C*Bz_old - S_ck*I*(kx_mod*Ey_old - ky_mod*Ex_old)
+ + X1*I*(kx_mod*Jy - ky_mod*Jx);
+ });
+ }
+}
+
+void ComovingPsatdAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& spectral_kspace,
+ const amrex::DistributionMapping& dm,
+ const amrex::Real dt)
+{
+ const amrex::BoxArray& ba = spectral_kspace.spectralspace_ba;
+
+ // Loop over boxes and allocate the corresponding coefficients for each box
+ for (amrex::MFIter mfi(ba, dm); mfi.isValid(); ++mfi) {
+
+ const amrex::Box& bx = ba[mfi];
+
+ // Extract pointers for the k vectors
+ const amrex::Real* kx_mod = modified_kx_vec[mfi].dataPtr();
+ const amrex::Real* kx = kx_vec[mfi].dataPtr();
+#if (AMREX_SPACEDIM==3)
+ const amrex::Real* ky_mod = modified_ky_vec[mfi].dataPtr();
+ const amrex::Real* ky = ky_vec[mfi].dataPtr();
+#endif
+ const amrex::Real* kz_mod = modified_kz_vec[mfi].dataPtr();
+ const amrex::Real* kz = kz_vec[mfi].dataPtr();
+
+ // Extract arrays for the coefficients
+ amrex::Array4<amrex::Real> C = C_coef [mfi].array();
+ amrex::Array4<amrex::Real> S_ck = S_ck_coef [mfi].array();
+ amrex::Array4<Complex> X1 = X1_coef [mfi].array();
+ amrex::Array4<Complex> X2 = X2_coef [mfi].array();
+ amrex::Array4<Complex> X3 = X3_coef [mfi].array();
+ amrex::Array4<Complex> X4 = X4_coef [mfi].array();
+ amrex::Array4<Complex> T2 = Theta2_coef[mfi].array();
+
+ // Store comoving velocity
+ const amrex::Real vx = m_v_comoving[0];
+#if (AMREX_SPACEDIM==3)
+ const amrex::Real vy = m_v_comoving[1];
+#endif
+ const amrex::Real vz = m_v_comoving[2];
+
+ // Loop over indices within one box
+ amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+ {
+ // Calculate norm of finite-order k vector
+ const amrex::Real knorm_mod = std::sqrt(
+ std::pow(kx_mod[i], 2) +
+#if (AMREX_SPACEDIM==3)
+ std::pow(ky_mod[j], 2) +
+ std::pow(kz_mod[k], 2));
+#else
+ std::pow(kz_mod[j], 2));
+#endif
+ // Calculate norm of infinite-order k vector
+ const amrex::Real knorm = std::sqrt(
+ std::pow(kx[i], 2) +
+#if (AMREX_SPACEDIM==3)
+ std::pow(ky[j], 2) +
+ std::pow(kz[k], 2));
+#else
+ std::pow(kz[j], 2));
+#endif
+ // Physical constants c, c**2, and epsilon_0, and imaginary unit
+ constexpr amrex::Real c = PhysConst::c;
+ constexpr amrex::Real c2 = c*c;
+ constexpr amrex::Real ep0 = PhysConst::ep0;
+ constexpr Complex I = Complex{0._rt, 1._rt};
+
+ // Auxiliary coefficients used when update_with_rho=false
+ const amrex::Real dt2 = dt * dt;
+
+ // Calculate dot product of k vector with comoving velocity
+ const amrex::Real kv = kx[i]*vx +
+#if (AMREX_SPACEDIM==3)
+ ky[j]*vy + kz[k]*vz;
+#else
+ kz[j]*vz;
+#endif
+
+ if (knorm_mod != 0. && knorm != 0.) {
+
+ // Auxiliary coefficients
+ const amrex::Real om_mod = c * knorm_mod;
+ const amrex::Real om2_mod = om_mod * om_mod;
+ const amrex::Real om = c * knorm;
+ const amrex::Real om2 = om * om;
+ const Complex tmp1 = amrex::exp( I * om_mod * dt);
+ const Complex tmp2 = amrex::exp(- I * om_mod * dt);
+ const Complex tmp1_sqrt = amrex::exp( I * om_mod * dt * 0.5_rt);
+ const Complex tmp2_sqrt = amrex::exp(- I * om_mod * dt * 0.5_rt);
+
+ C (i,j,k) = std::cos(om_mod * dt);
+ S_ck(i,j,k) = std::sin(om_mod * dt) / om_mod;
+
+ const amrex::Real nu = - kv / om;
+ const Complex theta = amrex::exp( I * nu * om * dt * 0.5_rt);
+ const Complex theta_star = amrex::exp(- I * nu * om * dt * 0.5_rt);
+
+ T2(i,j,k) = theta * theta;
+
+ if ( (nu != om_mod/om) && (nu != -om_mod/om) && (nu != 0.) ) {
+
+ Complex x1 = om2 / (om2_mod - nu * nu * om2)
+ * (theta_star - theta * C(i,j,k) + I * nu * om * theta * S_ck(i,j,k));
+
+ // X1 multiplies i*(k \times J) in the update equation for B
+ X1(i,j,k) = x1 / (ep0 * om2);
+
+ // X2 multiplies rho_new in the update equation for E
+ // X3 multiplies rho_old in the update equation for E
+ X2(i,j,k) = c2 * (x1 * om2_mod - theta * (1._rt - C(i,j,k)) * om2)
+ / (theta_star - theta) / (ep0 * om2 * om2_mod);
+ X3(i,j,k) = c2 * (x1 * om2_mod - theta_star * (1._rt - C(i,j,k)) * om2)
+ / (theta_star - theta) / (ep0 * om2 * om2_mod);
+
+ // X4 multiplies J in the update equation for E
+ X4(i,j,k) = I * nu * om * X1(i,j,k) - theta * S_ck(i,j,k) / ep0;
+ }
+
+ // Limits for nu = 0
+ if (nu == 0.) {
+
+ // X1 multiplies i*(k \times J) in the update equation for B
+ X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2_mod);
+
+ // X2 multiplies rho_new in the update equation for E
+ // X3 multiplies rho_old in the update equation for E
+ X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2_mod);
+ X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2_mod);
+
+ // Coefficient multiplying J in update equation for E
+ X4(i,j,k) = - S_ck(i,j,k) / ep0;
+ }
+
+ // Limits for nu = omega_mod/omega
+ if (nu == om_mod/om) {
+
+ // X1 multiplies i*(k \times J) in the update equation for B
+ X1(i,j,k) = tmp1_sqrt * (1._rt - tmp2 * tmp2 - 2._rt * I * om_mod * dt) / (4._rt * ep0 * om2_mod);
+
+ // X2 multiplies rho_new in the update equation for E
+ // X3 multiplies rho_old in the update equation for E
+ X2(i,j,k) = c2 * (- 4._rt + 3._rt * tmp1 + tmp2 - 2._rt * I * om_mod * dt * tmp1)
+ / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt));
+ X3(i,j,k) = c2 * (2._rt - tmp2 - 3._rt * tmp1 + 2._rt * tmp1 * tmp1 - 2._rt * I * om_mod * dt * tmp1)
+ / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt));
+
+ // Coefficient multiplying J in update equation for E
+ X4(i,j,k) = tmp1_sqrt * (I - I * tmp2 * tmp2 - 2._rt * om_mod * dt) / (4._rt * ep0 * om_mod);
+ }
+
+ // Limits for nu = -omega_mod/omega
+ if (nu == -om_mod/om) {
+
+ // X1 multiplies i*(k \times J) in the update equation for B
+ X1(i,j,k) = tmp2_sqrt * (1._rt - tmp1 * tmp1 + 2._rt * I * om_mod * dt) / (4._rt * ep0 * om2_mod);
+
+ // X2 multiplies rho_new in the update equation for E
+ // X3 multiplies rho_old in the update equation for E
+ X2(i,j,k) = c2 * (- 3._rt + 4._rt * tmp1 - tmp1 * tmp1 - 2._rt * I * om_mod * dt)
+ / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt));
+ X3(i,j,k) = c2 * (3._rt - 2._rt * tmp2 - 2._rt * tmp1 + tmp1 * tmp1 - 2._rt * I * om_mod * dt)
+ / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt));
+
+ // Coefficient multiplying J in update equation for E
+ X4(i,j,k) = tmp2_sqrt * (- I + I * tmp1 * tmp1 - 2._rt * om_mod * dt) / (4._rt * ep0 * om_mod);
+ }
+ }
+
+ // Limits for omega = 0 only
+ else if (knorm_mod != 0. && knorm == 0.) {
+
+ const amrex::Real om_mod = c * knorm_mod;
+ const amrex::Real om2_mod = om_mod * om_mod;
+
+ C (i,j,k) = std::cos(om_mod * dt);
+ S_ck(i,j,k) = std::sin(om_mod * dt) / om_mod;
+ T2(i,j,k) = 1._rt;
+
+ // X1 multiplies i*(k \times J) in the update equation for B
+ X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2_mod);
+
+ // X2 multiplies rho_new in the update equation for E
+ // X3 multiplies rho_old in the update equation for E
+ X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2_mod);
+ X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2_mod);
+
+ // Coefficient multiplying J in update equation for E
+ X4(i,j,k) = - S_ck(i,j,k) / ep0;
+
+ }
+
+ // Limits for omega_mod = 0 only
+ else if (knorm_mod == 0. && knorm != 0.) {
+
+ const amrex::Real om = c * knorm;
+ const amrex::Real om2 = om * om;
+ const amrex::Real nu = - kv / om;
+ const Complex theta = amrex::exp(I * nu * om * dt * 0.5_rt);
+ const Complex theta_star = amrex::exp(- I * nu * om * dt * 0.5_rt);
+
+ C(i,j,k) = 1._rt;
+ S_ck(i,j,k) = dt;
+ T2(i,j,k) = theta * theta;
+
+ if (nu != 0.) {
+ // X1 multiplies i*(k \times J) in the update equation for B
+ X1(i,j,k) = (-theta_star + theta - I * nu * om * dt * theta)
+ / (ep0 * nu * nu * om2);
+
+ // X2 multiplies rho_new in the update equation for E
+ // X3 multiplies rho_old in the update equation for E
+ X2(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om * dt * T2(i,j,k)
+ + 0.5_rt * nu * nu * om2 * dt * dt * T2(i,j,k))
+ / (ep0 * nu * nu * om2 * (T2(i,j,k) - 1._rt));
+ X3(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om * dt * T2(i,j,k)
+ + 0.5_rt * nu * nu * om2 * dt * dt)
+ / (ep0 * nu * nu * om2 * (T2(i,j,k) - 1._rt));
+
+ // Coefficient multiplying J in update equation for E
+ X4(i,j,k) = I * (theta - theta_star) / (ep0 * nu * om);
+ }
+
+ else {
+ // X1 multiplies i*(k \times J) in the update equation for B
+ X1(i,j,k) = dt2 / (2._rt * ep0);
+
+ // X2 multiplies rho_new in the update equation for E
+ // X3 multiplies rho_old in the update equation for E
+ X2(i,j,k) = c2 * dt2 / (6._rt * ep0);
+ X3(i,j,k) = - c2 * dt2 / (3._rt * ep0);
+
+ // Coefficient multiplying J in update equation for E
+ X4(i,j,k) = -dt / ep0;
+ }
+ }
+
+ // Limits for omega_mod = 0 and omega = 0
+ else if (knorm_mod == 0. && knorm == 0.) {
+
+ C(i,j,k) = 1._rt;
+ S_ck(i,j,k) = dt;
+ T2(i,j,k) = 1._rt;
+
+ // X1 multiplies i*(k \times J) in the update equation for B
+ X1(i,j,k) = dt2 / (2._rt * ep0);
+
+ // X2 multiplies rho_new in the update equation for E
+ // X3 multiplies rho_old in the update equation for E
+ X2(i,j,k) = c2 * dt2 / (6._rt * ep0);
+ X3(i,j,k) = - c2 * dt2 / (3._rt * ep0);
+
+ // Coefficient multiplying J in update equation for E
+ X4(i,j,k) = -dt / ep0;
+ }
+ });
+ }
+}
+
+void
+ComovingPsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
+ std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+ const std::unique_ptr<amrex::MultiFab>& rho)
+{
+ // Profiling
+ BL_PROFILE("ComovingAlgorithm::CurrentCorrection");
+
+ using Idx = SpectralFieldIndex;
+
+ // Forward Fourier transform of J and rho
+ field_data.ForwardTransform(*current[0], Idx::Jx, 0);
+ field_data.ForwardTransform(*current[1], Idx::Jy, 0);
+ field_data.ForwardTransform(*current[2], Idx::Jz, 0);
+ field_data.ForwardTransform(*rho, Idx::rho_old, 0);
+ field_data.ForwardTransform(*rho, Idx::rho_new, 1);
+
+ // Loop over boxes
+ for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
+
+ const amrex::Box& bx = field_data.fields[mfi].box();
+
+ // Extract arrays for the fields to be updated
+ amrex::Array4<Complex> fields = field_data.fields[mfi].array();
+
+ // Extract pointers for the k vectors
+ const amrex::Real* const modified_kx_arr = modified_kx_vec[mfi].dataPtr();
+ const amrex::Real* const kx_arr = kx_vec[mfi].dataPtr();
+#if (AMREX_SPACEDIM==3)
+ const amrex::Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr();
+ const amrex::Real* const ky_arr = ky_vec[mfi].dataPtr();
+#endif
+ const amrex::Real* const modified_kz_arr = modified_kz_vec[mfi].dataPtr();
+ const amrex::Real* const kz_arr = kz_vec[mfi].dataPtr();
+
+ // Local copy of member variables before GPU loop
+ const amrex::Real dt = m_dt;
+
+ // Comoving velocity
+ const amrex::Real vx = m_v_comoving[0];
+ const amrex::Real vy = m_v_comoving[1];
+ const amrex::Real vz = m_v_comoving[2];
+
+ // Loop over indices within one box
+ amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+ {
+ // Shortcuts for the values of J and rho
+ const Complex Jx = fields(i,j,k,Idx::Jx);
+ const Complex Jy = fields(i,j,k,Idx::Jy);
+ const Complex Jz = fields(i,j,k,Idx::Jz);
+ const Complex rho_old = fields(i,j,k,Idx::rho_old);
+ const Complex rho_new = fields(i,j,k,Idx::rho_new);
+
+ // k vector values, and coefficients
+ const amrex::Real kx_mod = modified_kx_arr[i];
+ const amrex::Real kx = kx_arr[i];
+#if (AMREX_SPACEDIM==3)
+ const amrex::Real ky_mod = modified_ky_arr[j];
+ const amrex::Real kz_mod = modified_kz_arr[k];
+ const amrex::Real ky = ky_arr[j];
+ const amrex::Real kz = kz_arr[k];
+#else
+ constexpr amrex::Real ky_mod = 0._rt;
+ const amrex::Real kz_mod = modified_kz_arr[j];
+ constexpr amrex::Real ky = 0._rt;
+ const amrex::Real kz = kz_arr[j];
+#endif
+ constexpr Complex I = Complex{0._rt,1._rt};
+
+ const amrex::Real knorm_mod = std::sqrt(kx_mod * kx_mod + ky_mod * ky_mod + kz_mod * kz_mod);
+
+ // Correct J
+ if (knorm_mod != 0._rt)
+ {
+ const Complex kmod_dot_J = kx_mod * Jx + ky_mod * Jy + kz_mod * Jz;
+ const amrex::Real k_dot_v = kx * vx + ky * vy + kz * vz;
+
+ if ( k_dot_v != 0._rt ) {
+
+ const Complex theta = amrex::exp(- I * k_dot_v * dt * 0.5_rt);
+ const Complex den = 1._rt - theta * theta;
+
+ fields(i,j,k,Idx::Jx) = Jx - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kx_mod / (knorm_mod * knorm_mod);
+ fields(i,j,k,Idx::Jy) = Jy - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * ky_mod / (knorm_mod * knorm_mod);
+ fields(i,j,k,Idx::Jz) = Jz - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kz_mod / (knorm_mod * knorm_mod);
+
+ } else {
+
+ fields(i,j,k,Idx::Jx) = Jx - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kx_mod / (knorm_mod * knorm_mod);
+ fields(i,j,k,Idx::Jy) = Jy - (kmod_dot_J - I * (rho_new - rho_old) / dt) * ky_mod / (knorm_mod * knorm_mod);
+ fields(i,j,k,Idx::Jz) = Jz - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kz_mod / (knorm_mod * knorm_mod);
+ }
+ }
+ });
+ }
+
+ // Backward Fourier transform of J
+ field_data.BackwardTransform(*current[0], Idx::Jx, 0);
+ field_data.BackwardTransform(*current[1], Idx::Jy, 0);
+ field_data.BackwardTransform(*current[2], Idx::Jz, 0);
+}
+
+void
+ComovingPsatdAlgorithm::VayDeposition (SpectralFieldData& /*field_data*/,
+ std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/)
+{
+ amrex::Abort("Vay deposition not implemented for comoving PSATD");
+}
+
+#endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
index 83ae83d69..cc43f0489 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
@@ -59,6 +59,12 @@ class GalileanAlgorithm : public SpectralBaseAlgorithm
private:
SpectralRealCoefficients C_coef, S_ck_coef;
SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef;
+ // Centered modified finite-order k vectors
+ KVectorComponent modified_kx_vec_centered;
+#if (AMREX_SPACEDIM==3)
+ KVectorComponent modified_ky_vec_centered;
+#endif
+ KVectorComponent modified_kz_vec_centered;
amrex::Array<amrex::Real,3> m_v_galilean;
amrex::Real m_dt;
bool m_update_with_rho;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
index 0757d17cd..ea5aa7ec0 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
@@ -18,6 +18,16 @@ GalileanAlgorithm::GalileanAlgorithm(const SpectralKSpace& spectral_kspace,
const bool update_with_rho)
// Initialize members of base class
: SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal),
+ // Initialize the centered finite-order modified k vectors: these are computed
+ // always with the assumption of centered grids (argument nodal = true),
+ // for both nodal and staggered simulations
+ modified_kx_vec_centered(spectral_kspace.getModifiedKComponent(dm,0,norder_x,true)),
+#if (AMREX_SPACEDIM==3)
+ modified_ky_vec_centered(spectral_kspace.getModifiedKComponent(dm,1,norder_y,true)),
+ modified_kz_vec_centered(spectral_kspace.getModifiedKComponent(dm,2,norder_z,true)),
+#else
+ modified_kz_vec_centered(spectral_kspace.getModifiedKComponent(dm,1,norder_z,true)),
+#endif
m_v_galilean(v_galilean),
m_dt(dt),
m_update_with_rho(update_with_rho)
@@ -91,8 +101,8 @@ GalileanAlgorithm::pushSpectralFields (SpectralFieldData& f) const
const Real ky = modified_ky_arr[j];
const Real kz = modified_kz_arr[k];
#else
- constexpr Real ky = 0;
- const Real kz = modified_kz_arr[j];
+ constexpr Real ky = 0._rt;
+ const Real kz = modified_kz_arr[j];
#endif
// Physical constant c**2 and imaginary unit
constexpr Real c2 = PhysConst::c*PhysConst::c;
@@ -168,11 +178,14 @@ void GalileanAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& sp
const Box& bx = ba[mfi];
// Extract pointers for the k vectors
- const Real* modified_kx = modified_kx_vec[mfi].dataPtr();
+ const Real* kx = modified_kx_vec[mfi].dataPtr();
+ const Real* kx_c = modified_kx_vec_centered[mfi].dataPtr();
#if (AMREX_SPACEDIM==3)
- const Real* modified_ky = modified_ky_vec[mfi].dataPtr();
+ const Real* ky = modified_ky_vec[mfi].dataPtr();
+ const Real* ky_c = modified_ky_vec_centered[mfi].dataPtr();
#endif
- const Real* modified_kz = modified_kz_vec[mfi].dataPtr();
+ const Real* kz = modified_kz_vec[mfi].dataPtr();
+ const Real* kz_c = modified_kz_vec_centered[mfi].dataPtr();
// Extract arrays for the coefficients
Array4<Real> C = C_coef[mfi].array();
@@ -194,13 +207,22 @@ void GalileanAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& sp
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
// Calculate norm of vector
- const Real k_norm = std::sqrt(
- std::pow(modified_kx[i], 2) +
+ const Real knorm = std::sqrt(
+ std::pow(kx[i], 2) +
#if (AMREX_SPACEDIM==3)
- std::pow(modified_ky[j], 2) +
- std::pow(modified_kz[k], 2));
+ std::pow(ky[j], 2) +
+ std::pow(kz[k], 2));
#else
- std::pow(modified_kz[j], 2));
+ std::pow(kz[j], 2));
+#endif
+ // Calculate norm of vector
+ const Real knorm_c = std::sqrt(
+ std::pow(kx_c[i], 2) +
+#if (AMREX_SPACEDIM==3)
+ std::pow(ky_c[j], 2) +
+ std::pow(kz_c[k], 2));
+#else
+ std::pow(kz_c[j], 2));
#endif
// Physical constants c, c**2, and epsilon_0, and imaginary unit
constexpr Real c = PhysConst::c;
@@ -213,142 +235,222 @@ void GalileanAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& sp
const Real dt3 = dt * dt2;
Complex X2_old, X3_old;
- // Calculate dot product of k vector with Galilean velocity
- const Real kv = modified_kx[i]*vx +
+ // Calculate dot product of k vector with Galilean velocity:
+ // this has to be computed always with the centered finite-order modified k
+ // vectors, in order to work correctly for both nodal and staggered simulations
+ const Real kv = kx_c[i]*vx +
#if (AMREX_SPACEDIM==3)
- modified_ky[j]*vy + modified_kz[k]*vz;
+ ky_c[j]*vy + kz_c[k]*vz;
#else
- modified_kz[j]*vz;
+ kz_c[j]*vz;
#endif
// The coefficients in the following refer to the ones given in equations
// (12a)-(12d) of (Lehe et al, PRE 94, 2016), used to update B and E
// (equations (11a) and (11b) of the same reference, respectively)
- if (k_norm != 0.) {
+ if (knorm != 0. && knorm_c != 0.) {
// Auxiliary coefficients
- const Real k2 = k_norm * k_norm;
- const Real ck = c * k_norm;
- const Real ckdt = ck * dt;
- const Complex tmp1 = amrex::exp( I * ckdt); // limit of T2 for nu = 1
- const Complex tmp2 = amrex::exp(- I * ckdt); // limit of T2 for nu = -1
+ const Real om = c * knorm;
+ const Real om2 = om * om;
+ const Real om3 = om * om2;
+ const Real om_c = c * knorm_c;
+ const Real om2_c = om_c * om_c;
+ const Real om3_c = om_c * om2_c;
+ const Complex tmp1 = amrex::exp( I * om * dt);
+ const Complex tmp2 = amrex::exp(- I * om * dt);
// See equation (12a)
- C (i,j,k) = std::cos(ckdt);
- S_ck(i,j,k) = std::sin(ckdt) / ck;
+ C (i,j,k) = std::cos(om * dt);
+ S_ck(i,j,k) = std::sin(om * dt) / om;
// See equation (12b)
- const Real nu = kv / ck;
- const Complex theta = amrex::exp( I * 0.5_rt * kv * dt);
- const Complex theta_star = amrex::exp(- I * 0.5_rt * kv * dt);
+ const Real nu = kv / om_c;
+ const Complex theta = amrex::exp( I * nu * om_c * dt * 0.5_rt);
+ const Complex theta_star = amrex::exp(- I * nu * om_c * dt * 0.5_rt);
// This is exp(i*(k \dot v_gal)*dt)
T2(i,j,k) = theta * theta;
- if ( (nu != 1.) && (nu != 0.) ) {
+ if ( (nu != om/om_c) && (nu != -om/om_c) && (nu != 0.) ) {
// x1 is the coefficient chi_1 in equation (12c)
- Complex x1 = 1._rt / (1._rt - nu*nu)
- * (theta_star - C(i,j,k) * theta + I * kv * S_ck(i,j,k) * theta);
+ Complex x1 = om2_c / (om2 - nu * nu * om2_c)
+ * (theta_star - theta * C(i,j,k) + I * nu * om_c * theta * S_ck(i,j,k));
// X1 multiplies i*(k \times J) in the update equation for B
- X1(i,j,k) = theta * x1 / (ep0 * c2 * k2);
+ X1(i,j,k) = theta * x1 / (ep0 * om2_c);
if (update_with_rho) {
// X2 multiplies rho_new in the update equation for E
// X3 multiplies rho_old in the update equation for E
- X2(i,j,k) = (x1 - theta * (1._rt - C(i,j,k))) / (theta_star - theta) / (ep0 * k2);
- X3(i,j,k) = (x1 - theta_star * (1._rt - C(i,j,k))) / (theta_star - theta) / (ep0 * k2);
+ X2(i,j,k) = c2 * (x1 * om2 - theta * (1._rt - C(i,j,k)) * om2_c)
+ / (theta_star - theta) / (ep0 * om2_c * om2);
+ X3(i,j,k) = c2 * (x1 * om2 - theta_star * (1._rt - C(i,j,k)) * om2_c)
+ / (theta_star - theta) / (ep0 * om2_c * om2);
} else {
// X2_old is the coefficient chi_2 in equation (12d)
// X3_old is the coefficient chi_3 in equation (12d)
// X2 multiplies (k \dot E) in the update equation for E
// X3 multiplies (k \dot J) in the update equation for E
- X2_old = (x1 - theta * (1._rt - C(i,j,k))) / (theta_star - theta);
- X3_old = (x1 - theta_star * (1._rt - C(i,j,k))) / (theta_star - theta);
- X2(i,j,k) = T2(i,j,k) * (X2_old - X3_old) / k2;
- X3(i,j,k) = I * X2_old * (T2(i,j,k) - 1._rt) / (ep0 * k2 * kv);
+ X2_old = (x1 * om2 - theta * (1._rt - C(i,j,k)) * om2_c)
+ / (theta_star - theta);
+ X3_old = (x1 * om2 - theta_star * (1._rt - C(i,j,k)) * om2_c)
+ / (theta_star - theta);
+ X2(i,j,k) = c2 * T2(i,j,k) * (X2_old - X3_old)
+ / (om2_c * om2);
+ X3(i,j,k) = I * c2 * X2_old * (T2(i,j,k) - 1._rt)
+ / (ep0 * nu * om3_c * om2);
}
// X4 multiplies J in the update equation for E
- X4(i,j,k) = I * kv * X1(i,j,k) - T2(i,j,k) * S_ck(i,j,k) / ep0;
+ X4(i,j,k) = I * nu * om_c * X1(i,j,k) - T2(i,j,k) * S_ck(i,j,k) / ep0;
}
// Limits for nu = 0
if (nu == 0.) {
// X1 multiplies i*(k \times J) in the update equation for B
- X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * c2 * k2);
+ X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2);
if (update_with_rho) {
// X2 multiplies rho_new in the update equation for E
// X3 multiplies rho_old in the update equation for E
- X2(i,j,k) = (1._rt - S_ck(i,j,k) / dt) / (ep0 * k2);
- X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * k2);
+ X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2);
+ X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2);
} else {
// X2 multiplies (k \dot E) in the update equation for E
// X3 multiplies (k \dot J) in the update equation for E
- X2(i,j,k) = (1._rt - C(i,j,k)) / k2;
- X3(i,j,k) = (S_ck(i,j,k) / dt - 1._rt) * dt / (ep0 * k2);
+ X2(i,j,k) = c2 * (1._rt - C(i,j,k)) / om2;
+ X3(i,j,k) = c2 * (S_ck(i,j,k) / dt - 1._rt) * dt / (ep0 * om2);
}
// Coefficient multiplying J in update equation for E
X4(i,j,k) = - S_ck(i,j,k) / ep0;
}
- // Limits for nu = 1
- if (nu == 1.) {
+ // Limits for nu = omega/omega_c
+ if (nu == om/om_c) {
// X1 multiplies i*(k \times J) in the update equation for B
- X1(i,j,k) = (1._rt - tmp1 * tmp1 + 2._rt * I * ckdt) / (4._rt * ep0 * c2 * k2);
+ X1(i,j,k) = (1._rt - tmp1 * tmp1 + 2._rt * I * om * dt) / (4._rt * ep0 * om2);
if (update_with_rho) {
// X2 multiplies rho_new in the update equation for E
// X3 multiplies rho_old in the update equation for E
- X2(i,j,k) = (- 3._rt + 4._rt * tmp1 - tmp1 * tmp1 - 2._rt * I * ckdt)
- / (4._rt * ep0 * k2 * (tmp1 - 1._rt));
- X3(i,j,k) = (3._rt - 2._rt / tmp1 - 2._rt * tmp1 + tmp1 * tmp1 - 2._rt * I * ckdt)
- / (4._rt * ep0 * k2 * (tmp1 - 1._rt));
+ X2(i,j,k) = c2 * (- 3._rt + 4._rt * tmp1 - tmp1 * tmp1 - 2._rt * I * om * dt)
+ / (4._rt * ep0 * om2 * (tmp1 - 1._rt));
+ X3(i,j,k) = c2 * (3._rt - 2._rt * tmp2 - 2._rt * tmp1 + tmp1 * tmp1 - 2._rt * I * om * dt)
+ / (4._rt * ep0 * om2 * (tmp1 - 1._rt));
} else {
// X2 multiplies (k \dot E) in the update equation for E
// X3 multiplies (k \dot J) in the update equation for E
- X2(i,j,k) = (1._rt - C(i,j,k)) * tmp1 / k2;
- X3(i,j,k) = (2._rt * ckdt - I * tmp1 * tmp1 + 4._rt * I * tmp1 - 3._rt * I)
- / (4._rt * ep0 * ck * k2);
+ X2(i,j,k) = c2 * (1._rt - C(i,j,k)) * tmp1 / om2;
+ X3(i,j,k) = c2 * (2._rt * om * dt - I * tmp1 * tmp1 + 4._rt * I * tmp1 - 3._rt * I)
+ / (4._rt * ep0 * om3);
}
// Coefficient multiplying J in update equation for E
- X4(i,j,k) = (- I + I * tmp1 * tmp1 - 2._rt * ckdt) / (4._rt * ep0 * ck);
+ X4(i,j,k) = (- I + I * tmp1 * tmp1 - 2._rt * om * dt) / (4._rt * ep0 * om);
}
- // Limits for nu = -1
- if (nu == -1.) {
+ // Limits for nu = -omega/omega_c
+ if (nu == -om/om_c) {
// X1 multiplies i*(k \times J) in the update equation for B
- X1(i,j,k) = (1._rt - tmp2 * tmp2 - 2._rt * I * ckdt) / (4._rt * ep0 * c2 * k2);
+ X1(i,j,k) = (1._rt - tmp2 * tmp2 - 2._rt * I * om * dt) / (4._rt * ep0 * om2);
if (update_with_rho) {
// X2 multiplies rho_new in the update equation for E
// X3 multiplies rho_old in the update equation for E
- X2(i,j,k) = (- 4._rt + 3._rt * tmp1 + tmp2 - 2._rt * I * ckdt * tmp1)
- / (4._rt * ep0 * k2 * (tmp1 - 1._rt));
- X3(i,j,k) = (2._rt - tmp2 - 3._rt * tmp1 + 2._rt * tmp1 * tmp1 - 2._rt * I * ckdt * tmp1)
- / (4._rt * ep0 * k2 * (tmp1 - 1._rt));
+ X2(i,j,k) = c2 * (- 4._rt + 3._rt * tmp1 + tmp2 - 2._rt * I * om * dt * tmp1)
+ / (4._rt * ep0 * om2 * (tmp1 - 1._rt));
+ X3(i,j,k) = c2 * (2._rt - tmp2 - 3._rt * tmp1 + 2._rt * tmp1 * tmp1 - 2._rt * I * om * dt * tmp1)
+ / (4._rt * ep0 * om2 * (tmp1 - 1._rt));
} else {
// X2 multiplies (k \dot E) in the update equation for E
// X3 multiplies (k \dot J) in the update equation for E
- X2(i,j,k) = (1._rt - C(i,j,k)) * tmp2 / k2;
- X3(i,j,k) = (2._rt * ckdt + I * tmp2 * tmp2 - 4._rt * I * tmp2 + 3._rt * I)
- / (4._rt * ep0 * ck * k2);
+ X2(i,j,k) = c2 * (1._rt - C(i,j,k)) * tmp2 / om2;
+ X3(i,j,k) = c2 * (2._rt * om * dt + I * tmp2 * tmp2 - 4._rt * I * tmp2 + 3._rt * I)
+ / (4._rt * ep0 * om3);
}
// Coefficient multiplying J in update equation for E
- X4(i,j,k) = (I - I * tmp2 * tmp2 - 2._rt * ckdt) / (4._rt * ep0 * ck);
+ X4(i,j,k) = (I - I * tmp2 * tmp2 - 2._rt * om * dt) / (4._rt * ep0 * om);
+ }
+ }
+
+ // Limits for omega_c = 0 only
+ else if (knorm != 0. && knorm_c == 0.) {
+
+ const Real om = c * knorm;
+ const Real om2 = om * om;
+
+ C (i,j,k) = std::cos(om * dt);
+ S_ck(i,j,k) = std::sin(om * dt) / om;
+ T2(i,j,k) = 1._rt;
+
+ // X1 multiplies i*(k \times J) in the update equation for B
+ X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2);
+
+ if (update_with_rho) {
+ // X2 multiplies rho_new in the update equation for E
+ // X3 multiplies rho_old in the update equation for E
+ X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2);
+ X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2);
+ } else {
+ // X2 multiplies (k \dot E) in the update equation for E
+ // X3 multiplies (k \dot J) in the update equation for E
+ X2(i,j,k) = c2 * (1._rt - C(i,j,k)) / om2;
+ X3(i,j,k) = c2 * (S_ck(i,j,k) / dt - 1._rt) * dt / (ep0 * om2);
}
+
+ // Coefficient multiplying J in update equation for E
+ X4(i,j,k) = - S_ck(i,j,k) / ep0;
+
+ }
+
+ // Limits for omega = 0 only
+ else if (knorm == 0. && knorm_c != 0.) {
+
+ const Real om_c = c * knorm_c;
+ const Real om2_c = om_c * om_c;
+ const Real om3_c = om_c * om2_c;
+ const Real nu = kv / om_c;
+ const Complex theta = amrex::exp(I * nu * om_c * dt * 0.5_rt);
+
+ C(i,j,k) = 1._rt;
+ S_ck(i,j,k) = dt;
+ T2(i,j,k) = theta * theta;
+
+ // X1 multiplies i*(k \times J) in the update equation for B
+ X1(i,j,k) = (-1._rt + T2(i,j,k) - I * nu * om_c * dt * T2(i,j,k))
+ / (ep0 * nu * nu * om2_c);
+
+ if (update_with_rho) {
+ // X2 multiplies rho_new in the update equation for E
+ // X3 multiplies rho_old in the update equation for E
+ X2(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om_c * dt * T2(i,j,k)
+ + 0.5_rt * nu * nu * om2_c * dt * dt * T2(i,j,k))
+ / (ep0 * nu * nu * om2_c * (T2(i,j,k) - 1._rt));
+ X3(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om_c * dt * T2(i,j,k)
+ + 0.5_rt * nu * nu * om2_c * dt * dt)
+ / (ep0 * nu * nu * om2_c * (T2(i,j,k) - 1._rt));
+ } else {
+ // X2 multiplies (k \dot E) in the update equation for E
+ // X3 multiplies (k \dot J) in the update equation for E
+ X2(i,j,k) = c2 * dt * dt * T2(i,j,k) * 0.5_rt;
+ X3(i,j,k) = c2 * (2._rt * I - 2._rt * nu * om_c * dt * T2(i,j,k)
+ + I * nu * nu * om2_c * dt * dt * T2(i,j,k))
+ / (2._rt * ep0 * nu * nu * nu * om3_c);
+ }
+
+ // Coefficient multiplying J in update equation for E
+ X4(i,j,k) = I * (T2(i,j,k) - 1._rt) / (ep0 * nu * om_c);
}
- // Limits for k = 0
- else {
+ // Limits for omega = 0 and omega_c = 0
+ else if (knorm == 0. && knorm_c == 0.) {
// Limits of cos(c*k*dt) and sin(c*k*dt)/(c*k)
C(i,j,k) = 1._rt;
@@ -404,11 +506,14 @@ GalileanAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
amrex::Array4<Complex> fields = field_data.fields[mfi].array();
// Extract pointers for the k vectors
- const amrex::Real* const modified_kx_arr = modified_kx_vec[mfi].dataPtr();
+ const amrex::Real* const modified_kx_arr = modified_kx_vec[mfi].dataPtr();
+ const amrex::Real* const modified_kx_arr_c = modified_kx_vec_centered[mfi].dataPtr();
#if (AMREX_SPACEDIM==3)
- const amrex::Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr();
+ const amrex::Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr();
+ const amrex::Real* const modified_ky_arr_c = modified_ky_vec_centered[mfi].dataPtr();
#endif
- const amrex::Real* const modified_kz_arr = modified_kz_vec[mfi].dataPtr();
+ const amrex::Real* const modified_kz_arr = modified_kz_vec[mfi].dataPtr();
+ const amrex::Real* const modified_kz_arr_c = modified_kz_vec_centered[mfi].dataPtr();
// Local copy of member variables before GPU loop
const amrex::Real dt = m_dt;
@@ -429,13 +534,18 @@ GalileanAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
const Complex rho_new = fields(i,j,k,Idx::rho_new);
// k vector values, and coefficients
- const amrex::Real kx = modified_kx_arr[i];
+ const amrex::Real kx = modified_kx_arr[i];
+ const amrex::Real kx_c = modified_kx_arr_c[i];
#if (AMREX_SPACEDIM==3)
- const amrex::Real ky = modified_ky_arr[j];
- const amrex::Real kz = modified_kz_arr[k];
+ const amrex::Real ky = modified_ky_arr[j];
+ const amrex::Real kz = modified_kz_arr[k];
+ const amrex::Real ky_c = modified_ky_arr_c[j];
+ const amrex::Real kz_c = modified_kz_arr_c[k];
#else
- constexpr amrex::Real ky = 0._rt;
- const amrex::Real kz = modified_kz_arr[j];
+ constexpr amrex::Real ky = 0._rt;
+ const amrex::Real kz = modified_kz_arr[j];
+ constexpr amrex::Real ky_c = 0._rt;
+ const amrex::Real kz_c = modified_kz_arr_c[j];
#endif
constexpr Complex I = Complex{0._rt,1._rt};
@@ -445,7 +555,7 @@ GalileanAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
if (k_norm != 0._rt)
{
const Complex k_dot_J = kx * Jx + ky * Jy + kz * Jz;
- const amrex::Real k_dot_vg = kx * vgx + ky * vgy + kz * vgz;
+ const amrex::Real k_dot_vg = kx_c * vgx + ky_c * vgy + kz_c * vgz;
if ( k_dot_vg != 0._rt ) {
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
index 57f8c5e3f..71b97b8bc 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
@@ -3,6 +3,7 @@ CEXE_sources += PsatdAlgorithm.cpp
CEXE_sources += GalileanAlgorithm.cpp
CEXE_sources += AvgGalileanAlgorithm.cpp
CEXE_sources += PMLPsatdAlgorithm.cpp
+CEXE_sources += ComovingPsatdAlgorithm.cpp
ifeq ($(USE_RZ),TRUE)
CEXE_sources += SpectralBaseAlgorithmRZ.cpp
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
index 6e3338fa8..2626c6982 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
@@ -33,6 +33,7 @@ class SpectralSolver
const int norder_x, const int norder_y,
const int norder_z, const bool nodal,
const amrex::Array<amrex::Real,3>& v_galilean,
+ const amrex::Array<amrex::Real,3>& v_comoving,
const amrex::RealVect dx, const amrex::Real dt,
const bool pml=false,
const bool periodic_single_box=false,
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
index 0cfd899df..4c9cb6967 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -10,6 +10,7 @@
#include "SpectralAlgorithms/GalileanAlgorithm.H"
#include "SpectralAlgorithms/AvgGalileanAlgorithm.H"
#include "SpectralAlgorithms/PMLPsatdAlgorithm.H"
+#include "SpectralAlgorithms/ComovingPsatdAlgorithm.H"
#include "WarpX.H"
#include "Utils/WarpXProfilerWrapper.H"
#include "Utils/WarpXUtil.H"
@@ -39,6 +40,7 @@ SpectralSolver::SpectralSolver(
const int norder_x, const int norder_y,
const int norder_z, const bool nodal,
const amrex::Array<amrex::Real,3>& v_galilean,
+ const amrex::Array<amrex::Real,3>& v_comoving,
const amrex::RealVect dx, const amrex::Real dt,
const bool pml, const bool periodic_single_box,
const bool update_with_rho,
@@ -64,15 +66,21 @@ SpectralSolver::SpectralSolver(
k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt);
}
else {
- if ((v_galilean[0]==0) && (v_galilean[1]==0) && (v_galilean[2]==0)){
- // v_galilean is 0: use standard PSATD algorithm
- algorithm = std::make_unique<PsatdAlgorithm>(
- k_space, dm, norder_x, norder_y, norder_z, nodal, dt, update_with_rho);
- }
- else {
+ // Galilean PSATD algorithm
+ if (v_galilean[0] != 0. || v_galilean[1] != 0. || v_galilean[2] != 0.) {
algorithm = std::make_unique<GalileanAlgorithm>(
k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt, update_with_rho);
}
+ // Comoving PSATD algorithm
+ else if (v_comoving[0] != 0. || v_comoving[1] != 0. || v_comoving[2] != 0.) {
+ algorithm = std::make_unique<ComovingPsatdAlgorithm>(
+ k_space, dm, norder_x, norder_y, norder_z, nodal, v_comoving, dt, update_with_rho);
+ }
+ // Standard PSATD algorithm
+ else {
+ algorithm = std::make_unique<PsatdAlgorithm>(
+ k_space, dm, norder_x, norder_y, norder_z, nodal, dt, update_with_rho);
+ }
}
}
diff --git a/Source/Parallelization/GuardCellManager.H b/Source/Parallelization/GuardCellManager.H
index 5e9971544..d2f764fbb 100644
--- a/Source/Parallelization/GuardCellManager.H
+++ b/Source/Parallelization/GuardCellManager.H
@@ -47,6 +47,7 @@ public:
const int maxwell_solver_id,
const int max_level,
const amrex::Array<amrex::Real,3> v_galilean,
+ const amrex::Array<amrex::Real,3> v_comoving,
const bool safe_guard_cells);
// Guard cells allocated for MultiFabs E and B
diff --git a/Source/Parallelization/GuardCellManager.cpp b/Source/Parallelization/GuardCellManager.cpp
index 1b8d38f0c..85490f2e0 100644
--- a/Source/Parallelization/GuardCellManager.cpp
+++ b/Source/Parallelization/GuardCellManager.cpp
@@ -28,6 +28,7 @@ guardCellManager::Init(
const int maxwell_solver_id,
const int max_level,
const amrex::Array<amrex::Real,3> v_galilean,
+ const amrex::Array<amrex::Real,3> v_comoving,
const bool safe_guard_cells)
{
// When using subcycling, the particles on the fine level perform two pushes
@@ -37,10 +38,9 @@ guardCellManager::Init(
int ngy_tmp = (max_level > 0 && do_subcycling == 1) ? nox+1 : nox;
int ngz_tmp = (max_level > 0 && do_subcycling == 1) ? nox+1 : nox;
- if ((v_galilean[0]!=0) or
- (v_galilean[1]!=0) or
- (v_galilean[2]!=0)){
- // Add one guard cell in the case of the galilean algorithm
+ // Add one guard cell in the case of the Galilean or comoving algorithms
+ if (v_galilean[0] != 0. || v_galilean[1] != 0. || v_galilean[2] != 0. ||
+ v_comoving[0] != 0. || v_comoving[1] != 0. || v_comoving[2] != 0. ) {
ngx_tmp += 1;
ngy_tmp += 1;
ngz_tmp += 1;
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp
index e6d4c91d4..6f408613d 100644
--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -10,6 +10,9 @@
#include "WarpX.H"
#include "WarpXSumGuardCells.H"
#include "Utils/CoarsenMR.H"
+#ifdef WARPX_USE_PSATD
+#include "FieldSolver/SpectralSolver/SpectralKSpace.H"
+#endif
#include <algorithm>
#include <cstdlib>
@@ -70,6 +73,42 @@ WarpX::UpdateAuxilaryData ()
void
WarpX::UpdateAuxilaryDataStagToNodal ()
{
+#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;
+
+ // Compute real-space stencil coefficients along x
+ amrex::Vector<Real> h_stencil_coef_x = getFornbergStencilCoefficients(fg_nox, false);
+ amrex::Gpu::DeviceVector<Real> d_stencil_coef_x(h_stencil_coef_x.size());
+ amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+ h_stencil_coef_x.begin(),
+ h_stencil_coef_x.end(),
+ d_stencil_coef_x.begin());
+ amrex::Gpu::synchronize();
+ amrex::Real const* p_stencil_coef_x = d_stencil_coef_x.data();
+
+ // Compute real-space stencil coefficients along y
+ amrex::Vector<Real> h_stencil_coef_y = getFornbergStencilCoefficients(fg_noy, false);
+ amrex::Gpu::DeviceVector<Real> d_stencil_coef_y(h_stencil_coef_y.size());
+ amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+ h_stencil_coef_y.begin(),
+ h_stencil_coef_y.end(),
+ d_stencil_coef_y.begin());
+ amrex::Gpu::synchronize();
+ amrex::Real const* p_stencil_coef_y = d_stencil_coef_y.data();
+
+ // Compute real-space stencil coefficients along z
+ amrex::Vector<Real> h_stencil_coef_z = getFornbergStencilCoefficients(fg_noz, false);
+ amrex::Gpu::DeviceVector<Real> d_stencil_coef_z(h_stencil_coef_z.size());
+ amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+ h_stencil_coef_z.begin(),
+ h_stencil_coef_z.end(),
+ d_stencil_coef_z.begin());
+ amrex::Gpu::synchronize();
+ amrex::Real const* p_stencil_coef_z = d_stencil_coef_z.data();
+#endif
+
// For level 0, we only need to do the average.
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
@@ -90,19 +129,34 @@ WarpX::UpdateAuxilaryDataStagToNodal ()
Array4<Real const> const& ey_fp = Efield_fp[0][1]->const_array(mfi);
Array4<Real const> const& ez_fp = Efield_fp[0][2]->const_array(mfi);
- const Box& bx = mfi.fabbox();
- amrex::ParallelFor(bx,
- [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
+ Box bx = mfi.validbox();
+ // TODO It seems like it is necessary to loop over the valid box grown
+ // with 2 guard cells. Should this number of guard cells be expressed
+ // in terms of the parameters defined in the guardCellManager class?
+ bx.grow(2);
+ amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
{
+// FDTD variant
+#ifndef WARPX_USE_PSATD
warpx_interp_nd_bfield_x(j,k,l, bx_aux, bx_fp);
warpx_interp_nd_bfield_y(j,k,l, by_aux, by_fp);
warpx_interp_nd_bfield_z(j,k,l, bz_aux, bz_fp);
warpx_interp_nd_efield_x(j,k,l, ex_aux, ex_fp);
warpx_interp_nd_efield_y(j,k,l, ey_aux, ey_fp);
warpx_interp_nd_efield_z(j,k,l, ez_aux, ez_fp);
+// PSATD variant
+#else
+ warpx_interp_nd_bfield_x(j,k,l, bx_aux, bx_fp, fg_noy, fg_noz, p_stencil_coef_y, p_stencil_coef_z);
+ warpx_interp_nd_bfield_y(j,k,l, by_aux, by_fp, fg_nox, fg_noz, p_stencil_coef_x, p_stencil_coef_z);
+ warpx_interp_nd_bfield_z(j,k,l, bz_aux, bz_fp, fg_nox, fg_noy, p_stencil_coef_x, p_stencil_coef_y);
+ warpx_interp_nd_efield_x(j,k,l, ex_aux, ex_fp, fg_nox, p_stencil_coef_x);
+ warpx_interp_nd_efield_y(j,k,l, ey_aux, ey_fp, fg_noy, p_stencil_coef_y);
+ warpx_interp_nd_efield_z(j,k,l, ez_aux, ez_fp, fg_noz, p_stencil_coef_z);
+#endif
});
}
+ // NOTE: high-order interpolation is not implemented for mesh refinement
for (int lev = 1; lev <= finest_level; ++lev)
{
BoxArray const& nba = Bfield_aux[lev][0]->boxArray();
diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H
index b490636ee..7a8095f18 100644
--- a/Source/Parallelization/WarpXComm_K.H
+++ b/Source/Parallelization/WarpXComm_K.H
@@ -378,6 +378,10 @@ void warpx_interp_nd_bfield_z (int j, int k, int l,
Bza(j,k,l) = bg + (bf-bc);
}
+// With the FDTD Maxwell solver, this is the linear interpolation function used to
+// interpolate the Bx field from a staggered grid to a nodal grid, before gathering
+// the field from the nodes to the particle positions (momentum-conserving gathering)
+#ifndef WARPX_USE_PSATD
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void warpx_interp_nd_bfield_x (int j, int k, int l,
amrex::Array4<amrex::Real> const& Bxa,
@@ -391,6 +395,45 @@ void warpx_interp_nd_bfield_x (int j, int k, int l,
#endif
}
+// With the PSATD Maxwell solver, this is the arbitrary-order interpolation function used to
+// interpolate the Bx field from a staggered grid to a nodal grid, before gathering
+// the field from the nodes to the particle positions (momentum-conserving gathering)
+#else
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_bfield_x (int j, int k, int l,
+ amrex::Array4<amrex::Real> const& Bxa,
+ amrex::Array4<amrex::Real const> const& Bxf,
+ const int noy,
+ const int noz,
+ amrex::Real const* stencil_coef_y,
+ amrex::Real const* stencil_coef_z)
+{
+#if (AMREX_SPACEDIM == 2)
+ amrex::ignore_unused(noy);
+ amrex::ignore_unused(stencil_coef_y);
+ amrex::Real res = 0.;
+ for (int nz = 1; nz < noz/2+1; nz++) {
+ res += 0.5 * stencil_coef_z[nz-1] * (Bxf(j, k + nz - 1, l) + Bxf(j, k - nz, l));
+ }
+ Bxa(j,k,l) = res;
+#else
+ amrex::Real res = 0.;
+ for (int nz = 1; nz < noz/2+1; nz++) {
+ for (int ny = 1; ny < noy/2+1; ny++) {
+ res += 0.25 * stencil_coef_y[ny-1] * stencil_coef_z[nz-1] *
+ (Bxf(j, k + ny - 1, l + nz - 1) + Bxf(j, k - ny, l + nz - 1)
+ + Bxf(j, k + ny - 1, l - nz ) + Bxf(j, k - ny, l - nz ));
+ }
+ }
+ Bxa(j,k,l) = res;
+#endif
+}
+#endif
+
+// With the FDTD Maxwell solver, this is the linear interpolation function used to
+// interpolate the By field from a staggered grid to a nodal grid, before gathering
+// the field from the nodes to the particle positions (momentum-conserving gathering)
+#ifndef WARPX_USE_PSATD
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void warpx_interp_nd_bfield_y (int j, int k, int l,
amrex::Array4<amrex::Real> const& Bya,
@@ -404,6 +447,47 @@ void warpx_interp_nd_bfield_y (int j, int k, int l,
#endif
}
+// With the PSATD Maxwell solver, this is the arbitrary-order interpolation function used to
+// interpolate the By field from a staggered grid to a nodal grid, before gathering
+// the field from the nodes to the particle positions (momentum-conserving gathering)
+#else
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_bfield_y (int j, int k, int l,
+ amrex::Array4<amrex::Real> const& Bya,
+ amrex::Array4<amrex::Real const> const& Byf,
+ const int nox,
+ const int noz,
+ amrex::Real const* stencil_coef_x,
+ amrex::Real const* stencil_coef_z)
+{
+#if (AMREX_SPACEDIM == 2)
+ amrex::Real res = 0.;
+ for (int nz = 1; nz < noz/2+1; nz++) {
+ for (int nx = 1; nx < nox/2+1; nx++) {
+ res += 0.25 * stencil_coef_x[nx-1] * stencil_coef_z[nz-1] *
+ (Byf(j + nx - 1, k + nz - 1, l) + Byf(j - nx, k + nz - 1, l)
+ + Byf(j + nx - 1, k - nz , l) + Byf(j - nx, k - nz , l));
+ }
+ }
+ Bya(j,k,l) = res;
+#else
+ amrex::Real res = 0.;
+ for (int nz = 1; nz < noz/2+1; nz++) {
+ for (int nx = 1; nx < nox/2+1; nx++) {
+ res += 0.25 * stencil_coef_x[nx-1] * stencil_coef_z[nz-1] *
+ (Byf(j + nx - 1, k, l + nz - 1) + Byf(j - nx, k, l + nz - 1)
+ + Byf(j + nx - 1, k, l - nz ) + Byf(j - nx, k, l - nz ));
+ }
+ }
+ Bya(j,k,l) = res;
+#endif
+}
+#endif
+
+// With the FDTD Maxwell solver, this is the linear interpolation function used to
+// interpolate the Bz field from a staggered grid to a nodal grid, before gathering
+// the field from the nodes to the particle positions (momentum-conserving gathering)
+#ifndef WARPX_USE_PSATD
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void warpx_interp_nd_bfield_z (int j, int k, int l,
amrex::Array4<amrex::Real> const& Bza,
@@ -417,6 +501,41 @@ void warpx_interp_nd_bfield_z (int j, int k, int l,
#endif
}
+// With the PSATD Maxwell solver, this is the arbitrary-order interpolation function used to
+// interpolate the Bz field from a staggered grid to a nodal grid, before gathering
+// the field from the nodes to the particle positions (momentum-conserving gathering)
+#else
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_bfield_z (int j, int k, int l,
+ amrex::Array4<amrex::Real> const& Bza,
+ amrex::Array4<amrex::Real const> const& Bzf,
+ const int nox,
+ const int noy,
+ amrex::Real const* stencil_coef_x,
+ amrex::Real const* stencil_coef_y)
+{
+#if (AMREX_SPACEDIM == 2)
+ amrex::ignore_unused(noy);
+ amrex::ignore_unused(stencil_coef_y);
+ amrex::Real res = 0.;
+ for (int nx = 1; nx < nox/2+1; nx++) {
+ res += 0.5 * stencil_coef_x[nx-1] * (Bzf(j + nx - 1, k, l) + Bzf(j - nx, k, l));
+ }
+ Bza(j,k,l) = res;
+#else
+ amrex::Real res = 0.;
+ for (int ny = 1; ny < noy/2+1; ny++) {
+ for (int nx = 1; nx < nox/2+1; nx++) {
+ res += 0.25 * stencil_coef_x[nx-1] * stencil_coef_y[ny-1] *
+ (Bzf(j + nx - 1, k + ny - 1, l) + Bzf(j - nx, k + ny - 1, l)
+ + Bzf(j + nx - 1, k - ny , l) + Bzf(j - nx, k - ny , l));
+ }
+ }
+ Bza(j,k,l) = res;
+#endif
+}
+#endif
+
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void warpx_interp_nd_efield_x (int j, int k, int l,
amrex::Array4<amrex::Real> const& Exa,
@@ -621,6 +740,10 @@ void warpx_interp_nd_efield_z (int j, int k, int l,
Eza(j,k,l) = eg + (ef-ec);
}
+// With the FDTD Maxwell solver, this is the linear interpolation function used to
+// interpolate the Ex field from a staggered grid to a nodal grid, before gathering
+// the field from the nodes to the particle positions (momentum-conserving gathering)
+#ifndef WARPX_USE_PSATD
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void warpx_interp_nd_efield_x (int j, int k, int l,
amrex::Array4<amrex::Real> const& Exa,
@@ -629,6 +752,29 @@ void warpx_interp_nd_efield_x (int j, int k, int l,
Exa(j,k,l) = 0.5*(Exf(j-1,k,l) + Exf(j,k,l));
}
+// With the PSATD Maxwell solver, this is the arbitrary-order interpolation function used to
+// interpolate the Ex field from a staggered grid to a nodal grid, before gathering
+// the field from the nodes to the particle positions (momentum-conserving gathering)
+#else
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_efield_x (int j, int k, int l,
+ amrex::Array4<amrex::Real> const& Exa,
+ amrex::Array4<amrex::Real const> const& Exf,
+ const int nox,
+ amrex::Real const* stencil_coef_x)
+{
+ amrex::Real res = 0.;
+ for (int nx = 1; nx < nox/2+1; nx++) {
+ res += 0.5 * stencil_coef_x[nx-1] * (Exf(j + nx - 1, k, l) + Exf(j - nx, k, l));
+ }
+ Exa(j,k,l) = res;
+}
+#endif
+
+// With the FDTD Maxwell solver, this is the linear interpolation function used to
+// interpolate the Ey field from a staggered grid to a nodal grid, before gathering
+// the field from the nodes to the particle positions (momentum-conserving gathering)
+#ifndef WARPX_USE_PSATD
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void warpx_interp_nd_efield_y (int j, int k, int l,
amrex::Array4<amrex::Real> const& Eya,
@@ -642,6 +788,35 @@ void warpx_interp_nd_efield_y (int j, int k, int l,
#endif
}
+// With the PSATD Maxwell solver, this is the arbitrary-order interpolation function used to
+// interpolate the Ey field from a staggered grid to a nodal grid, before gathering
+// the field from the nodes to the particle positions (momentum-conserving gathering)
+#else
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_efield_y (int j, int k, int l,
+ amrex::Array4<amrex::Real> const& Eya,
+ amrex::Array4<amrex::Real const> const& Eyf,
+ const int noy,
+ amrex::Real const* stencil_coef_y)
+{
+#if (AMREX_SPACEDIM == 2)
+ amrex::ignore_unused(noy);
+ amrex::ignore_unused(stencil_coef_y);
+ Eya(j,k,l) = Eyf(j,k,l);
+#else
+ amrex::Real res = 0.;
+ for (int ny = 1; ny < noy/2+1; ny++) {
+ res += 0.5 * stencil_coef_y[ny-1] * (Eyf(j, k + ny - 1, l) + Eyf(j, k - ny, l));
+ }
+ Eya(j,k,l) = res;
+#endif
+}
+#endif
+
+// With the FDTD Maxwell solver, this is the linear interpolation function used to
+// interpolate the Ez field from a staggered grid to a nodal grid, before gathering
+// the field from the nodes to the particle positions (momentum-conserving gathering)
+#ifndef WARPX_USE_PSATD
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void warpx_interp_nd_efield_z (int j, int k, int l,
amrex::Array4<amrex::Real> const& Eza,
@@ -655,4 +830,31 @@ void warpx_interp_nd_efield_z (int j, int k, int l,
#endif
}
+// With the PSATD Maxwell solver, this is the arbitrary-order interpolation function used to
+// interpolate the Ez field from a staggered grid to a nodal grid, before gathering
+// the field from the nodes to the particle positions (momentum-conserving gathering)
+#else
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_efield_z (int j, int k, int l,
+ amrex::Array4<amrex::Real> const& Eza,
+ amrex::Array4<amrex::Real const> const& Ezf,
+ const int noz,
+ amrex::Real const* stencil_coef_z)
+{
+#if (AMREX_SPACEDIM == 2)
+ amrex::Real res = 0.;
+ for (int nz = 1; nz < noz/2+1; nz++) {
+ res += 0.5 * stencil_coef_z[nz-1] * (Ezf(j, k + nz - 1, l) + Ezf(j, k - nz, l));
+ }
+ Eza(j,k,l) = res;
+#else
+ amrex::Real res = 0.;
+ for (int nz = 1; nz < noz/2+1; nz++) {
+ res += 0.5 * stencil_coef_z[nz-1] * (Ezf(j, k, l + nz - 1) + Ezf(j, k, l - nz));
+ }
+ Eza(j,k,l) = res;
+#endif
+}
+#endif
+
#endif
diff --git a/Source/WarpX.H b/Source/WarpX.H
index b66f99b82..949247ac9 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -157,6 +157,12 @@ 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;
+
// Number of modes for the RZ multimode version
static long n_rz_azimuthal_modes;
static long ncomps;
@@ -274,6 +280,8 @@ public:
amrex::Array<amrex::Real,3> m_v_galilean = {{0}};
amrex::Array<amrex::Real,3> m_galilean_shift = {{0}};
+ amrex::Array<amrex::Real,3> m_v_comoving = {{0.}};
+
static int num_mirrors;
amrex::Vector<amrex::Real> mirror_z;
amrex::Vector<amrex::Real> mirror_z_width;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index a8afa17ec..3c057e044 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -81,6 +81,12 @@ 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;
+
bool WarpX::use_fdtd_nci_corr = false;
bool WarpX::galerkin_interpolation = true;
@@ -642,6 +648,22 @@ WarpX::ReadParameters ()
pp.query("noy", noy);
pp.query("noz", noz);
+#ifdef WARPX_USE_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.query("field_gathering_nox", field_gathering_nox);
+ pp.query("field_gathering_noy", field_gathering_noy);
+ pp.query("field_gathering_noz", field_gathering_noz);
+ }
+
+ 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");
+ }
+#endif
+
pp.query("galerkin_scheme",galerkin_interpolation);
AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nox == noy and nox == noz ,
@@ -688,25 +710,40 @@ WarpX::ReadParameters ()
pp.query("current_correction", current_correction);
pp.query("v_galilean", m_v_galilean);
+ pp.query("v_comoving", m_v_comoving);
pp.query("do_time_averaging", fft_do_time_averaging);
- // Scale the velocity by the speed of light
+ // Scale the velocity by the speed of light
for (int i=0; i<3; i++) m_v_galilean[i] *= PhysConst::c;
+ for (int i=0; i<3; i++) m_v_comoving[i] *= PhysConst::c;
+
+ // The comoving PSATD algorithm is not implemented nor tested with Esirkepov current deposition
+ if (current_deposition_algo == CurrentDepositionAlgo::Esirkepov) {
+ if (m_v_comoving[0] != 0. || m_v_comoving[1] != 0. || m_v_comoving[2] != 0.) {
+ amrex::Abort("Esirkepov current deposition cannot be used with the comoving PSATD algorithm");
+ }
+ }
# ifdef WARPX_DIM_RZ
update_with_rho = true; // Must be true for RZ PSATD
# else
- if (m_v_galilean[0] == 0. && m_v_galilean[1] == 0. && m_v_galilean[2] == 0.) {
+ if (m_v_galilean[0] == 0. && m_v_galilean[1] == 0. && m_v_galilean[2] == 0. &&
+ m_v_comoving[0] == 0. && m_v_comoving[1] == 0. && m_v_comoving[2] == 0.) {
update_with_rho = false; // standard PSATD
}
else {
- update_with_rho = true; // Galilean PSATD
+ update_with_rho = true; // Galilean PSATD or comoving PSATD
}
# endif
// Overwrite update_with_rho with value set in input file
pp.query("update_with_rho", update_with_rho);
+ if (m_v_comoving[0] != 0. || m_v_comoving[1] != 0. || m_v_comoving[2] != 0.) {
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(update_with_rho,
+ "psatd.update_with_rho must be equal to 1 for comoving PSATD");
+ }
+
# ifdef WARPX_DIM_RZ
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(update_with_rho,
"psatd.update_with_rho must be equal to 1 in RZ geometry");
@@ -883,6 +920,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d
maxwell_solver_id,
maxLevel(),
WarpX::m_v_galilean,
+ WarpX::m_v_comoving,
safe_guard_cells);
if (mypc->nSpeciesDepositOnMainGrid() && n_current_deposition_buffer == 0) {
@@ -1074,7 +1112,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
}
bool const pml_flag_false=false;
spectral_solver_fp[lev] = std::make_unique<SpectralSolver>( realspace_ba, dm,
- nox_fft, noy_fft, noz_fft, do_nodal, m_v_galilean, dx_vect, dt[lev],
+ nox_fft, noy_fft, noz_fft, do_nodal, m_v_galilean, m_v_comoving, dx_vect, dt[lev],
pml_flag_false, fft_periodic_single_box, update_with_rho, fft_do_time_averaging );
# endif
#endif
@@ -1192,7 +1230,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
# else
c_realspace_ba.grow(ngE); // add guard cells
spectral_solver_cp[lev] = std::make_unique<SpectralSolver>( c_realspace_ba, dm,
- nox_fft, noy_fft, noz_fft, do_nodal, m_v_galilean, cdx_vect, dt[lev],
+ nox_fft, noy_fft, noz_fft, do_nodal, m_v_galilean, m_v_comoving, cdx_vect, dt[lev],
pml_flag_false, fft_periodic_single_box, update_with_rho, fft_do_time_averaging );
# endif
#endif
generated by cgit v1.2.3 (git 2.39.1) at 2025-08-09 09:59:58 +0000