aboutsummaryrefslogtreecommitdiff
path: root/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp')
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp176
1 files changed, 8 insertions, 168 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
index 2b1e09c2d..0515aa598 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
@@ -39,7 +39,6 @@ PsatdAlgorithm::PsatdAlgorithm(
const amrex::Real dt,
const bool update_with_rho,
const bool time_averaging,
- const bool do_multi_J,
const bool dive_cleaning,
const bool divb_cleaning)
// Initializer list
@@ -59,7 +58,6 @@ PsatdAlgorithm::PsatdAlgorithm(
m_dt(dt),
m_update_with_rho(update_with_rho),
m_time_averaging(time_averaging),
- m_do_multi_J(do_multi_J),
m_dive_cleaning(dive_cleaning),
m_divb_cleaning(divb_cleaning)
{
@@ -84,7 +82,7 @@ PsatdAlgorithm::PsatdAlgorithm(
InitializeSpectralCoefficients(spectral_kspace, dm, dt);
// Allocate these coefficients only with time averaging
- if (time_averaging && !do_multi_J)
+ if (time_averaging)
{
Psi1_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
Psi2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
@@ -94,14 +92,6 @@ PsatdAlgorithm::PsatdAlgorithm(
Y4_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
InitializeSpectralCoefficientsAveraging(spectral_kspace, dm, dt);
}
- // Allocate these coefficients only with time averaging
- // and with the assumption that J is linear in time (always with multi-J algorithm)
- else if (time_averaging && do_multi_J)
- {
- X5_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
- X6_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
- InitializeSpectralCoefficientsAvgLin(spectral_kspace, dm, dt);
- }
if (dive_cleaning && m_is_galilean)
{
@@ -121,12 +111,11 @@ PsatdAlgorithm::PsatdAlgorithm(
void
PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
{
- const bool update_with_rho = m_update_with_rho;
- const bool time_averaging = m_time_averaging;
- const bool do_multi_J = m_do_multi_J;
- const bool dive_cleaning = m_dive_cleaning;
- const bool divb_cleaning = m_divb_cleaning;
- const bool is_galilean = m_is_galilean;
+ const bool update_with_rho = m_update_with_rho;
+ const bool time_averaging = m_time_averaging;
+ const bool dive_cleaning = m_dive_cleaning;
+ const bool divb_cleaning = m_divb_cleaning;
+ const bool is_galilean = m_is_galilean;
const amrex::Real dt = m_dt;
@@ -163,7 +152,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
amrex::Array4<const Complex> Y3_arr;
amrex::Array4<const Complex> Y4_arr;
- if (time_averaging && !do_multi_J)
+ if (time_averaging)
{
Psi1_arr = Psi1_coef[mfi].array();
Psi2_arr = Psi2_coef[mfi].array();
@@ -173,14 +162,6 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
Y4_arr = Y4_coef[mfi].array();
}
- Array4<const Complex> X5_arr;
- Array4<const Complex> X6_arr;
- if (time_averaging && do_multi_J)
- {
- X5_arr = X5_coef[mfi].array();
- X6_arr = X6_coef[mfi].array();
- }
-
// Extract pointers for the k vectors
const amrex::Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr();
#if defined(WARPX_DIM_3D)
@@ -229,7 +210,6 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
#endif
// Physical constants and imaginary unit
constexpr Real c2 = PhysConst::c * PhysConst::c;
- constexpr Real ep0 = PhysConst::ep0;
constexpr Real inv_ep0 = 1._rt / PhysConst::ep0;
constexpr Complex I = Complex{0._rt, 1._rt};
@@ -320,78 +300,8 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
fields(i,j,k,Idx.G) = C * G_old + I * c2 * S_ck * k_dot_B;
}
- if (do_multi_J)
- {
- const Complex Jx_new = fields(i,j,k,Idx.Jx_new);
- const Complex Jy_new = fields(i,j,k,Idx.Jy_new);
- const Complex Jz_new = fields(i,j,k,Idx.Jz_new);
-
- fields(i,j,k,Idx.Ex) += -X1 * (Jx_new - Jx) / dt;
- fields(i,j,k,Idx.Ey) += -X1 * (Jy_new - Jy) / dt;
- fields(i,j,k,Idx.Ez) += -X1 * (Jz_new - Jz) / dt;
-
- fields(i,j,k,Idx.Bx) += I * X2/c2 * (ky * (Jz_new - Jz) - kz * (Jy_new - Jy));
- fields(i,j,k,Idx.By) += I * X2/c2 * (kz * (Jx_new - Jx) - kx * (Jz_new - Jz));
- fields(i,j,k,Idx.Bz) += I * X2/c2 * (kx * (Jy_new - Jy) - ky * (Jx_new - Jx));
-
- if (dive_cleaning)
- {
- const Complex k_dot_dJ = kx * (Jx_new - Jx) + ky * (Jy_new - Jy) + kz * (Jz_new - Jz);
-
- fields(i,j,k,Idx.F) += -I * X2/c2 * k_dot_dJ;
- }
-
- if (time_averaging)
- {
- const Complex X5 = X5_arr(i,j,k);
- const Complex X6 = X6_arr(i,j,k);
-
- // TODO: Here the code is *accumulating* the average,
- // because it is meant to be used with sub-cycling
- // maybe this should be made more generic
-
- fields(i,j,k,Idx.Ex_avg) += S_ck * Ex_old
- + I * c2 * ep0 * X1 * (ky * Bz_old - kz * By_old)
- + I * X5 * rho_old * kx + I * X6 * rho_new * kx + X3/c2 * Jx - X2/c2 * Jx_new;
-
- fields(i,j,k,Idx.Ey_avg) += S_ck * Ey_old
- + I * c2 * ep0 * X1 * (kz * Bx_old - kx * Bz_old)
- + I * X5 * rho_old * ky + I * X6 * rho_new * ky + X3/c2 * Jy - X2/c2 * Jy_new;
-
- fields(i,j,k,Idx.Ez_avg) += S_ck * Ez_old
- + I * c2 * ep0 * X1 * (kx * By_old - ky * Bx_old)
- + I * X5 * rho_old * kz + I * X6 * rho_new * kz + X3/c2 * Jz - X2/c2 * Jz_new;
-
- fields(i,j,k,Idx.Bx_avg) += S_ck * Bx_old
- - I * ep0 * X1 * (ky * Ez_old - kz * Ey_old)
- - I * X5/c2 * (ky * Jz - kz * Jy) - I * X6/c2 * (ky * Jz_new - kz * Jy_new);
-
- fields(i,j,k,Idx.By_avg) += S_ck * By_old
- - I * ep0 * X1 * (kz * Ex_old - kx * Ez_old)
- - I * X5/c2 * (kz * Jx - kx * Jz) - I * X6/c2 * (kz * Jx_new - kx * Jz_new);
-
- fields(i,j,k,Idx.Bz_avg) += S_ck * Bz_old
- - I * ep0 * X1 * (kx * Ey_old - ky * Ex_old)
- - I * X5/c2 * (kx * Jy - ky * Jx) - I * X6/c2 * (kx * Jy_new - ky * Jx_new);
-
- if (dive_cleaning)
- {
- fields(i,j,k,Idx.Ex_avg) += I * c2 * ep0 * X1 * F_old * kx;
- fields(i,j,k,Idx.Ey_avg) += I * c2 * ep0 * X1 * F_old * ky;
- fields(i,j,k,Idx.Ez_avg) += I * c2 * ep0 * X1 * F_old * kz;
- }
-
- if (divb_cleaning)
- {
- fields(i,j,k,Idx.Bx_avg) += I * ep0 * X1 * G_old * kx;
- fields(i,j,k,Idx.By_avg) += I * ep0 * X1 * G_old * ky;
- fields(i,j,k,Idx.Bz_avg) += I * ep0 * X1 * G_old * kz;
- }
- }
- }
-
// Additional update equations for averaged Galilean algorithm
- if (time_averaging && !do_multi_J)
+ if (time_averaging)
{
// These coefficients are initialized in the function InitializeSpectralCoefficients below
const Complex Psi1 = Psi1_arr(i,j,k);
@@ -822,76 +732,6 @@ void PsatdAlgorithm::InitializeSpectralCoefficientsAveraging (
}
}
-void PsatdAlgorithm::InitializeSpectralCoefficientsAvgLin (
- const SpectralKSpace& spectral_kspace,
- const amrex::DistributionMapping& dm,
- const amrex::Real dt)
-{
- const BoxArray& ba = spectral_kspace.spectralspace_ba;
-
- // Loop over boxes and allocate the corresponding coefficients for each box
- for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi)
- {
- const Box& bx = ba[mfi];
-
- // Extract pointers for the k vectors
- const Real* kx_s = modified_kx_vec[mfi].dataPtr();
-#if defined(WARPX_DIM_3D)
- const Real* ky_s = modified_ky_vec[mfi].dataPtr();
-#endif
- const Real* kz_s = modified_kz_vec[mfi].dataPtr();
-
- Array4<Real const> C = C_coef[mfi].array();
- Array4<Real const> S_ck = S_ck_coef[mfi].array();
-
- Array4<Complex> X5 = X5_coef[mfi].array();
- Array4<Complex> X6 = X6_coef[mfi].array();
-
- // Loop over indices within one box
- ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
- {
- // Calculate norm of k vector
- const Real knorm_s = std::sqrt(
- std::pow(kx_s[i], 2) +
-#if defined(WARPX_DIM_3D)
- std::pow(ky_s[j], 2) + std::pow(kz_s[k], 2));
-#else
- std::pow(kz_s[j], 2));
-#endif
- // Physical constants and imaginary unit
- constexpr Real c = PhysConst::c;
- constexpr Real c2 = c*c;
- constexpr Real ep0 = PhysConst::ep0;
-
- // Auxiliary coefficients
- const Real dt3 = dt * dt * dt;
-
- const Real om_s = c * knorm_s;
- const Real om2_s = om_s * om_s;
- const Real om4_s = om2_s * om2_s;
-
- if (om_s != 0.)
- {
- X5(i,j,k) = c2 / ep0 * (S_ck(i,j,k) / om2_s - (1._rt - C(i,j,k)) / (om4_s * dt)
- - 0.5_rt * dt / om2_s);
- }
- else
- {
- X5(i,j,k) = - c2 * dt3 / (8._rt * ep0);
- }
-
- if (om_s != 0.)
- {
- X6(i,j,k) = c2 / ep0 * ((1._rt - C(i,j,k)) / (om4_s * dt) - 0.5_rt * dt / om2_s);
- }
- else
- {
- X6(i,j,k) = - c2 * dt3 / (24._rt * ep0);
- }
- });
- }
-}
-
void
PsatdAlgorithm::CurrentCorrection (
const int lev,
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-11 22:49:02 +0000