aboutsummaryrefslogtreecommitdiff
path: root/Source/FieldSolver/SpectralSolver
diff options
context:
space:
mode:
authorGravatar Edoardo Zoni <59625522+EZoni@users.noreply.github.com> 2021-04-19 08:27:19 -0700
committerGravatar GitHub <noreply@github.com> 2021-04-19 08:27:19 -0700
commit0beb477484381a012d9e80f326c4cb51926b5e5c (patch)
treebe3b5e0d9374e2cb718582b4fdd64c13a1aa55ff /Source/FieldSolver/SpectralSolver
parente0d1e2494f1190f1b42088273409887a459987bd (diff)
downloadWarpX-0beb477484381a012d9e80f326c4cb51926b5e5c.tar.gz
WarpX-0beb477484381a012d9e80f326c4cb51926b5e5c.tar.zst
WarpX-0beb477484381a012d9e80f326c4cb51926b5e5c.zip
Divergence cleaning for PSATD in PML (#1600)
* Add div(E)/div(B) cleaning options for PSATD in PMLs * Pass missing flags to spectral solver in PML * Duplicate MPI exchange and communication functions for G * Use separate parameters for div cleaning in PMLs * Add asserts for features that are not implemented * Do not need to duplicate MPI exchange functions for G * Add short documentation for new input parameters * Set new parameters true by default with PSATD solver * Add CI test for PML div cleaning with PSATD * Use new syntax <diag_name>.intervals in new input file * Reset benchmark of new CI test * Always synchronize nodal points of G MultiFab * Fix few warnings in 2D build * Update Benchmark of pml_psatd_dive_divb_cleaning * Improve Documentation of warpx.do_pml_dive_cleaning Co-authored-by: Neïl Zaim <49716072+NeilZaim@users.noreply.github.com> * Improve Documentation of warpx.do_pml_divb_cleaning Co-authored-by: Neïl Zaim <49716072+NeilZaim@users.noreply.github.com> * Improve Abort Message * Clean Up Co-authored-by: Neïl Zaim <49716072+NeilZaim@users.noreply.github.com>
Diffstat (limited to 'Source/FieldSolver/SpectralSolver')
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H6
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp295
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralFieldData.H10
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralSolver.H6
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralSolver.cpp6
5 files changed, 246 insertions, 77 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
index 368ea27ad..1b9755826 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
@@ -21,7 +21,9 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm
const amrex::DistributionMapping& dm,
const int norder_x, const int norder_y,
const int norder_z, const bool nodal,
- const amrex::Real dt);
+ const amrex::Real dt,
+ const bool dive_cleaning,
+ const bool divb_cleaning);
void InitializeSpectralCoefficients(
const SpectralKSpace& spectral_kspace,
@@ -69,6 +71,8 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm
private:
SpectralRealCoefficients C_coef, S_ck_coef, inv_k2_coef;
amrex::Real m_dt;
+ bool m_dive_cleaning;
+ bool m_divb_cleaning;
};
#endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
index aafa3ba8d..3c1526357 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
@@ -18,10 +18,13 @@ using namespace amrex;
PMLPsatdAlgorithm::PMLPsatdAlgorithm(const SpectralKSpace& spectral_kspace,
const DistributionMapping& dm,
const int norder_x, const int norder_y,
- const int norder_z, const bool nodal, const Real dt)
+ const int norder_z, const bool nodal, const Real dt,
+ const bool dive_cleaning, const bool divb_cleaning)
// Initialize members of base class
: SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal),
- m_dt(dt)
+ m_dt(dt),
+ m_dive_cleaning(dive_cleaning),
+ m_divb_cleaning(divb_cleaning)
{
const BoxArray& ba = spectral_kspace.spectralspace_ba;
@@ -38,6 +41,9 @@ PMLPsatdAlgorithm::PMLPsatdAlgorithm(const SpectralKSpace& spectral_kspace,
void
PMLPsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const {
+ const bool dive_cleaning = m_dive_cleaning;
+ const bool divb_cleaning = m_divb_cleaning;
+
// Loop over boxes
for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){
@@ -61,8 +67,7 @@ PMLPsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const {
const amrex::Real dt = m_dt;
// Loop over indices within one box
- ParallelFor(bx,
- [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
// Record old values of the fields to be updated
using Idx = SpectralPMLIndex;
@@ -81,6 +86,56 @@ PMLPsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const {
const Complex Bzx = fields(i,j,k,Idx::Bzx);
const Complex Bzy = fields(i,j,k,Idx::Bzy);
+ Complex Ex, Ey, Ez;
+ Complex Bx, By, Bz;
+
+ // Used only if dive_cleaning = true and divb_cleaning = true
+ Complex Exx, Eyy, Ezz;
+ Complex Bxx, Byy, Bzz;
+ Complex Fx, Fy, Fz;
+ Complex Gx, Gy, Gz;
+ Complex F, G;
+
+ if (!dive_cleaning && !divb_cleaning)
+ {
+ Ex = Exy + Exz;
+ Ey = Eyx + Eyz;
+ Ez = Ezx + Ezy;
+
+ Bx = Bxy + Bxz;
+ By = Byx + Byz;
+ Bz = Bzx + Bzy;
+ }
+ else if (dive_cleaning && divb_cleaning)
+ {
+ Exx = fields(i,j,k,Idx::Exx);
+ Eyy = fields(i,j,k,Idx::Eyy);
+ Ezz = fields(i,j,k,Idx::Ezz);
+
+ Bxx = fields(i,j,k,Idx::Bxx);
+ Byy = fields(i,j,k,Idx::Byy);
+ Bzz = fields(i,j,k,Idx::Bzz);
+
+ Fx = fields(i,j,k,Idx::Fx);
+ Fy = fields(i,j,k,Idx::Fy);
+ Fz = fields(i,j,k,Idx::Fz);
+
+ Gx = fields(i,j,k,Idx::Gx);
+ Gy = fields(i,j,k,Idx::Gy);
+ Gz = fields(i,j,k,Idx::Gz);
+
+ Ex = Exx + Exy + Exz;
+ Ey = Eyx + Eyy + Eyz;
+ Ez = Ezx + Ezy + Ezz;
+
+ Bx = Bxx + Bxy + Bxz;
+ By = Byx + Byy + Byz;
+ Bz = Bzx + Bzy + Bzz;
+
+ F = Fx + Fy + Fz;
+ G = Gx + Gy + Gz;
+ }
+
// k vector values, and coefficients
const Real kx = modified_kx_arr[i];
#if (AMREX_SPACEDIM==3)
@@ -101,70 +156,174 @@ PMLPsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const {
if (k_norm != 0._rt) {
- const Real C = C_arr(i,j,k);
- const Real S_ck = S_ck_arr(i,j,k);
- const Real inv_k2 = inv_k2_arr(i,j,k);
-
- const Real C1 = (kx2 * C + ky2 + kz2) * inv_k2;
- const Real C2 = (kx2 + ky2 * C + kz2) * inv_k2;
- const Real C3 = (kx2 + ky2 + kz2 * C) * inv_k2;
- const Real C4 = kx2 * (C - 1._rt) * inv_k2;
- const Real C5 = ky2 * (C - 1._rt) * inv_k2;
- const Real C6 = kz2 * (C - 1._rt) * inv_k2;
- const Real C7 = ky * kz * (1._rt - C) * inv_k2;
- const Real C8 = kx * kz * (1._rt - C) * inv_k2;
- const Real C9 = kx * ky * (1._rt - C) * inv_k2;
- const Real C10 = c2 * kx * ky * kz * (dt - S_ck) * inv_k2;
- const Real C11 = c2 * ky2 * kz * (dt - S_ck) * inv_k2;
- const Real C12 = c2 * kz2 * ky * (dt - S_ck) * inv_k2;
- const Real C13 = c2 * kz2 * kx * (dt - S_ck) * inv_k2;
- const Real C14 = c2 * kx2 * kz * (dt - S_ck) * inv_k2;
- const Real C15 = c2 * kx2 * ky * (dt - S_ck) * inv_k2;
- const Real C16 = c2 * ky2 * kx * (dt - S_ck) * inv_k2;
- const Real C17 = c2 * kx * (ky2 * dt + (kz2 + kx2) * S_ck) * inv_k2;
- const Real C18 = c2 * kx * (kz2 * dt + (ky2 + kx2) * S_ck) * inv_k2;
- const Real C19 = c2 * ky * (kz2 * dt + (kx2 + ky2) * S_ck) * inv_k2;
- const Real C20 = c2 * ky * (kx2 * dt + (kz2 + ky2) * S_ck) * inv_k2;
- const Real C21 = c2 * kz * (kx2 * dt + (ky2 + kz2) * S_ck) * inv_k2;
- const Real C22 = c2 * kz * (ky2 * dt + (kx2 + kz2) * S_ck) * inv_k2;
-
- // Update E
- fields(i,j,k,Idx::Exy) = C2 * Exy + C5 * Exz + C9 * (Eyz + Eyx)
- + I*C10 * (Bxy + Bxz) + I*C11 * (Byz + Byx) + I*C19 * (Bzx + Bzy);
-
- fields(i,j,k,Idx::Exz) = C6 * Exy + C3 * Exz + C8 * (Ezx + Ezy)
- - I*C10 * (Bxy + Bxz) - I*C22 * (Byz + Byx) - I*C12 * (Bzx + Bzy);
-
- fields(i,j,k,Idx::Eyz) = C3 * Eyz + C6 * Eyx + C7 * (Ezx + Ezy)
- + I*C21 * (Bxy + Bxz) + I*C10 * (Byz + Byx) + I*C13 * (Bzx + Bzy);
-
- fields(i,j,k,Idx::Eyx) = C9 * (Exy + Exz) + C4 * Eyz + C1 * Eyx
- - I*C14 * (Bxy + Bxz) - I*C10 * (Byz + Byx) - I*C18 * (Bzx + Bzy);
-
- fields(i,j,k,Idx::Ezx) = C8 * (Exy + Exz) + C1 * Ezx + C4 * Ezy
- + I*C15 * (Bxy + Bxz) + I*C17 * (Byz + Byx) + I*C10 * (Bzx + Bzy);
-
- fields(i,j,k,Idx::Ezy) = C7 * (Eyz + Eyx) + C5 * Ezx + C2 * Ezy
- - I*C20 * (Bxy + Bxz) - I*C16 * (Byz + Byx) - I*C10 * (Bzx + Bzy);
-
- // Update B
- fields(i,j,k,Idx::Bxy) = C2 * Bxy + C5 * Bxz + C9 * (Byz + Byx)
- - I*C10/c2 * (Exy + Exz) - I*C11/c2 * (Eyz + Eyx) - I*C19/c2 * (Ezx + Ezy);
-
- fields(i,j,k,Idx::Bxz) = C6 * Bxy + C3 * Bxz + C8 * (Bzx + Bzy)
- + I*C10/c2 * (Exy + Exz) + I*C22/c2 * (Eyz + Eyx) + I*C12/c2 * (Ezx + Ezy);
-
- fields(i,j,k,Idx::Byz) = C3 * Byz + C6 * Byx + C7 * (Bzx + Bzy)
- - I*C21/c2 * (Exy + Exz) - I*C10/c2 * (Eyz + Eyx) - I*C13/c2 * (Ezx + Ezy);
-
- fields(i,j,k,Idx::Byx) = C9 * (Bxy + Bxz) + C4 * Byz + C1 * Byx
- + I*C14/c2 * (Exy + Exz) + I*C10/c2 * (Eyz + Eyx) + I*C18/c2 * (Ezx + Ezy);
-
- fields(i,j,k,Idx::Bzx) = C8 * (Bxy + Bxz) + C1 * Bzx + C4 * Bzy
- - I*C15/c2 * (Exy + Exz) - I*C17/c2 * (Eyz + Eyx) - I*C10/c2 * (Ezx + Ezy);
-
- fields(i,j,k,Idx::Bzy) = C7 * (Byz + Byx) + C5 * Bzx + C2 * Bzy
- + I*C20/c2 * (Exy + Exz) + I*C16/c2 * (Eyz + Eyx) + I*C10/c2 * (Ezx + Ezy);
+ const amrex::Real C = C_arr(i,j,k);
+ const amrex::Real S_ck = S_ck_arr(i,j,k);
+ const amrex::Real inv_k2 = inv_k2_arr(i,j,k);
+
+ const amrex::Real C1 = (kx2 * C + ky2 + kz2) * inv_k2;
+ const amrex::Real C2 = (kx2 + ky2 * C + kz2) * inv_k2;
+ const amrex::Real C3 = (kx2 + ky2 + kz2 * C) * inv_k2;
+ const amrex::Real C4 = kx2 * (C - 1._rt) * inv_k2;
+ const amrex::Real C5 = ky2 * (C - 1._rt) * inv_k2;
+ const amrex::Real C6 = kz2 * (C - 1._rt) * inv_k2;
+ const amrex::Real C7 = ky * kz * (1._rt - C) * inv_k2;
+ const amrex::Real C8 = kx * kz * (1._rt - C) * inv_k2;
+ const amrex::Real C9 = kx * ky * (1._rt - C) * inv_k2;
+
+ if (!dive_cleaning && !divb_cleaning)
+ {
+ const Complex C10 = I * c2 * kx * ky * kz * (dt - S_ck) * inv_k2;
+ const Complex C11 = I * c2 * ky2 * kz * (dt - S_ck) * inv_k2;
+ const Complex C12 = I * c2 * kz2 * ky * (dt - S_ck) * inv_k2;
+ const Complex C13 = I * c2 * kz2 * kx * (dt - S_ck) * inv_k2;
+ const Complex C14 = I * c2 * kx2 * kz * (dt - S_ck) * inv_k2;
+ const Complex C15 = I * c2 * kx2 * ky * (dt - S_ck) * inv_k2;
+ const Complex C16 = I * c2 * ky2 * kx * (dt - S_ck) * inv_k2;
+ const Complex C17 = I * c2 * kx * (ky2 * dt + (kz2 + kx2) * S_ck) * inv_k2;
+ const Complex C18 = I * c2 * kx * (kz2 * dt + (ky2 + kx2) * S_ck) * inv_k2;
+ const Complex C19 = I * c2 * ky * (kz2 * dt + (kx2 + ky2) * S_ck) * inv_k2;
+ const Complex C20 = I * c2 * ky * (kx2 * dt + (kz2 + ky2) * S_ck) * inv_k2;
+ const Complex C21 = I * c2 * kz * (kx2 * dt + (ky2 + kz2) * S_ck) * inv_k2;
+ const Complex C22 = I * c2 * kz * (ky2 * dt + (kx2 + kz2) * S_ck) * inv_k2;
+
+ const Complex C10_c2 = C10 / c2;
+ const Complex C11_c2 = C11 / c2;
+ const Complex C12_c2 = C12 / c2;
+ const Complex C13_c2 = C13 / c2;
+ const Complex C14_c2 = C14 / c2;
+ const Complex C15_c2 = C15 / c2;
+ const Complex C16_c2 = C16 / c2;
+ const Complex C17_c2 = C17 / c2;
+ const Complex C18_c2 = C18 / c2;
+ const Complex C19_c2 = C19 / c2;
+ const Complex C20_c2 = C20 / c2;
+ const Complex C21_c2 = C21 / c2;
+ const Complex C22_c2 = C22 / c2;
+
+ // Update E
+ fields(i,j,k,Idx::Exy) = C2 * Exy + C5 * Exz + C9 * Ey
+ + C10 * Bx + C11 * By + C19 * Bz;
+
+ fields(i,j,k,Idx::Exz) = C6 * Exy + C3 * Exz + C8 * Ez
+ - C10 * Bx - C22 * By - C12 * Bz;
+
+ fields(i,j,k,Idx::Eyz) = C3 * Eyz + C6 * Eyx + C7 * Ez
+ + C21 * Bx + C10 * By + C13 * Bz;
+
+ fields(i,j,k,Idx::Eyx) = C9 * Ex + C4 * Eyz + C1 * Eyx
+ - C14 * Bx - C10 * By - C18 * Bz;
+
+ fields(i,j,k,Idx::Ezx) = C8 * Ex + C1 * Ezx + C4 * Ezy
+ + C15 * Bx + C17 * By + C10 * Bz;
+
+ fields(i,j,k,Idx::Ezy) = C7 * Ey + C5 * Ezx + C2 * Ezy
+ - C20 * Bx - C16 * By - C10 * Bz;
+
+ // Update B
+ fields(i,j,k,Idx::Bxy) = C2 * Bxy + C5 * Bxz + C9 * By
+ - C10_c2 * Ex - C11_c2 * Ey - C19_c2 * Ez;
+
+ fields(i,j,k,Idx::Bxz) = C6 * Bxy + C3 * Bxz + C8 * Bz
+ + C10_c2 * Ex + C22_c2 * Ey + C12_c2 * Ez;
+
+ fields(i,j,k,Idx::Byz) = C3 * Byz + C6 * Byx + C7 * Bz
+ - C21_c2 * Ex - C10_c2 * Ey - C13_c2 * Ez;
+
+ fields(i,j,k,Idx::Byx) = C9 * Bx + C4 * Byz + C1 * Byx
+ + C14_c2 * Ex + C10_c2 * Ey + C18_c2 * Ez;
+
+ fields(i,j,k,Idx::Bzx) = C8 * Bx + C1 * Bzx + C4 * Bzy
+ - C15_c2 * Ex - C17_c2 * Ey - C10_c2 * Ez;
+
+ fields(i,j,k,Idx::Bzy) = C7 * By + C5 * Bzx + C2 * Bzy
+ + C20_c2 * Ex + C16_c2 * Ey + C10_c2 * Ez;
+ }
+ else if (dive_cleaning && divb_cleaning)
+ {
+ const Complex C23 = I * c2 * kx * S_ck;
+ const Complex C24 = I * c2 * ky * S_ck;
+ const Complex C25 = I * c2 * kz * S_ck;
+
+ const Complex C23_c2 = C23 / c2;
+ const Complex C24_c2 = C24 / c2;
+ const Complex C25_c2 = C25 / c2;
+
+ // Update E
+ fields(i,j,k,Idx::Exx) = C1 * Exx + C4 * Exy + C4 * Exz
+ - C9 * Ey - C8 * Ez + C23 * F;
+
+ fields(i,j,k,Idx::Exy) = C5 * Exx + C2 * Exy + C5 * Exz
+ + C9 * Ey + C24 * Bz - C7 * G;
+
+ fields(i,j,k,Idx::Exz) = C6 * Exx + C6 * Exy + C3 * Exz
+ + C8 * Ez - C25 * By + C7 * G;
+
+ fields(i,j,k,Idx::Eyx) = C9 * Ex + C1 * Eyx + C4 * Eyy
+ + C4 * Eyz - C23 * Bz + C8 * G;
+
+ fields(i,j,k,Idx::Eyy) = - C9 * Ex + C5 * Eyx + C2 * Eyy
+ + C5 * Eyz - C7 * Ez + C24 * F;
+
+ fields(i,j,k,Idx::Eyz) = C6 * Eyx + C6 * Eyy + C3 * Eyz
+ + C7 * Ez + C25 * Bx - C8 * G;
+
+ fields(i,j,k,Idx::Ezx) = C8 * Ex + C1 * Ezx + C4 * Ezy
+ + C4 * Ezz + C23 * By - C9 * G;
+
+ fields(i,j,k,Idx::Ezy) = C7 * Ey + C5 * Ezx + C2 * Ezy
+ + C5 * Ezz - C24 * Bx + C9 * G;
+
+ fields(i,j,k,Idx::Ezz) = - C8 * Ex - C7 * Ey + C6 * Ezx
+ + C6 * Ezy + C3 * Ezz + C25 * F;
+
+ // Update B
+ fields(i,j,k,Idx::Bxx) = C1 * Bxx + C4 * Bxy + C4 * Bxz
+ - C9 * By - C8 * Bz + C23_c2 * G;
+
+ fields(i,j,k,Idx::Bxy) = - C24_c2 * Ez + C5 * Bxx + C2 * Bxy
+ + C5 * Bxz + C9 * By + C7 * F;
+
+ fields(i,j,k,Idx::Bxz) = C25_c2 * Ey + C6 * Bxx + C6 * Bxy
+ + C3 * Bxz + C8 * Bz - C7 * F;
+
+ fields(i,j,k,Idx::Byx) = C23_c2 * Ez + C9 * Bx + C1 * Byx
+ + C4 * Byy + C4 * Byz - C8 * F;
+
+ fields(i,j,k,Idx::Byy) = - C9 * Bx + C5 * Byx + C2 * Byy
+ + C5 * Byz - C7 * Bz + C24_c2 * G;
+
+ fields(i,j,k,Idx::Byz) = - C25_c2 * Ex + C6 * Byx + C6 * Byy
+ + C3 * Byz + C7 * Bz + C8 * F;
+
+ fields(i,j,k,Idx::Bzx) = - C23_c2 * Ey + C8 * Bx + C1 * Bzx
+ + C4 * Bzy + C4 * Bzz + C9 * F;
+
+ fields(i,j,k,Idx::Bzy) = C24_c2 * Ex + C7 * By + C5 * Bzx
+ + C2 * Bzy + C5 * Bzz - C9 * F;
+
+ fields(i,j,k,Idx::Bzz) = - C8 * Bx - C7 * By + C6 * Bzx
+ + C6 * Bzy + C3 * Bzz + C25_c2 * G;
+
+ // Update F
+ fields(i,j,k,Idx::Fx) = C23_c2 * Ex + C8 * By - C9 * Bz
+ + C1 * Fx + C4 * Fy + C4 * Fz;
+
+ fields(i,j,k,Idx::Fy) = C24_c2 * Ey - C7 * Bx + C9 * Bz
+ + C5 * Fx + C2 * Fy + C5 * Fz;
+
+ fields(i,j,k,Idx::Fz) = C25_c2 * Ez + C7 * Bx - C8 * By
+ + C6 * Fx + C6 * Fy + C3 * Fz;
+
+ // Update G
+ fields(i,j,k,Idx::Gx) = - C8 * Ey + C9 * Ez + C23 * Bx
+ + C1 * Gx + C4 * Gy + C4 * Gz;
+
+ fields(i,j,k,Idx::Gy) = C7 * Ex - C9 * Ez + C24 * By
+ + C5 * Gx + C2 * Gy + C5 * Gz;
+
+ fields(i,j,k,Idx::Gz) = - C7 * Ex + C8 * Ey + C25 * Bz
+ + C6 * Gx + C6 * Gy + C3 * Gz;
+ }
}
});
}
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
index 764ecc8af..4fc734853 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
@@ -33,11 +33,13 @@ struct SpectralAvgFieldIndex {
// n_fields is automatically the total number of fields
};
-/* Index for the PML fields, when stored in spectral space */
+/** Index for the PML fields, when stored in spectral space,
+ * (n_fields is automatically set to the total number of fields)
+ * TODO How to include the diagonal components only when needed?
+ */
struct SpectralPMLIndex {
- enum { Exy=0, Exz, Eyx, Eyz, Ezx, Ezy,
- Bxy, Bxz, Byx, Byz, Bzx, Bzy, n_fields };
- // n_fields is automatically the total number of fields
+ enum {Exx=0, Exy, Exz, Eyx, Eyy, Eyz, Ezx, Ezy, Ezz,
+ Bxx , Bxy, Bxz, Byx, Byy, Byz, Bzx, Bzy, Bzz, Fx, Fy, Fz, Gx, Gy, Gz, n_fields};
};
/** \brief Class that stores the fields in spectral space, and performs the
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
index c9a7bbcc4..d5c7725fe 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
@@ -28,7 +28,7 @@ class SpectralSolver
// underlying classes `SpectralFieldData` and `PsatdAlgorithm`
// Constructor
- SpectralSolver( const int lev,
+ SpectralSolver (const int lev,
const amrex::BoxArray& realspace_ba,
const amrex::DistributionMapping& dm,
const int norder_x, const int norder_y,
@@ -40,7 +40,9 @@ class SpectralSolver
const bool pml=false,
const bool periodic_single_box=false,
const bool update_with_rho=false,
- const bool fft_do_time_averaging=false);
+ const bool fft_do_time_averaging=false,
+ const bool dive_cleaning = false,
+ const bool divb_cleaning = false);
/**
* \brief Transform the component `i_comp` of MultiFab `mf`
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
index ec1219f93..ae631ba6c 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -43,7 +43,9 @@ SpectralSolver::SpectralSolver(
const amrex::RealVect dx, const amrex::Real dt,
const bool pml, const bool periodic_single_box,
const bool update_with_rho,
- const bool fft_do_time_averaging) {
+ const bool fft_do_time_averaging,
+ const bool dive_cleaning,
+ const bool divb_cleaning) {
// Initialize all structures using the same distribution mapping dm
@@ -57,7 +59,7 @@ SpectralSolver::SpectralSolver(
if (pml) {
algorithm = std::make_unique<PMLPsatdAlgorithm>(
- k_space, dm, norder_x, norder_y, norder_z, nodal, dt);
+ k_space, dm, norder_x, norder_y, norder_z, nodal, dt, dive_cleaning, divb_cleaning);
}
else {
// Comoving PSATD algorithm