diff options
author | 2021-04-19 08:27:19 -0700 | |
---|---|---|
committer | 2021-04-19 08:27:19 -0700 | |
commit | 0beb477484381a012d9e80f326c4cb51926b5e5c (patch) | |
tree | be3b5e0d9374e2cb718582b4fdd64c13a1aa55ff /Source/FieldSolver/SpectralSolver | |
parent | e0d1e2494f1190f1b42088273409887a459987bd (diff) | |
download | WarpX-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')
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 |