diff options
Diffstat (limited to 'Source/FieldSolver/SpectralSolver')
5 files changed, 307 insertions, 51 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index 7b917284b..b7f349bb0 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -41,6 +41,8 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm * \param[in] dt time step of the simulation * \param[in] update_with_rho whether the update equation for E uses rho or not * \param[in] time_averaging whether to use time averaging for large time steps + * \param[in] J_linear_in_time whether to use two currents computed at the beginning and the end + * of the time interval (instead of using one current computed at half time) */ PsatdAlgorithm ( const SpectralKSpace& spectral_kspace, @@ -52,7 +54,8 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm const amrex::Array<amrex::Real,3>& v_galilean, const amrex::Real dt, const bool update_with_rho, - const bool time_averaging); + const bool time_averaging, + const bool J_linear_in_time); /** * \brief Updates the E and B fields in spectral space, according to the relevant PSATD equations @@ -66,12 +69,22 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm */ virtual int getRequiredNumberOfFields () const override final { - if (m_time_averaging) { - return SpectralAvgFieldIndex::n_fields; - } else { - return SpectralFieldIndex::n_fields; + if (m_J_linear_in_time) + { + return SpectralFieldIndexJLinearInTime::n_fields; } - } + else + { + if (m_time_averaging) + { + return SpectralFieldIndexTimeAveraging::n_fields; + } + else + { + return SpectralFieldIndex::n_fields; + } + } + }; /** * \brief Initializes the coefficients used in \c pushSpectralFields to update the E and B fields @@ -86,6 +99,19 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm const amrex::Real dt); /** + * \brief Initialize additional coefficients used in \c pushSpectralFields to update E,B, + * required only when using time averaging with the assumption that J is linear in time + * + * \param[in] spectral_kspace spectral space + * \param[in] dm distribution mapping + * \param[in] dt time step of the simulation + */ + void InitializeSpectralCoefficientsAvgLin ( + const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const amrex::Real dt); + + /** * \brief Initializes additional coefficients used in \c pushSpectralFields to update the E and B fields, * required only when using time averaging with large time steps * @@ -138,6 +164,8 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm SpectralRealCoefficients C_coef, S_ck_coef; SpectralComplexCoefficients T2_coef, X1_coef, X2_coef, X3_coef, X4_coef; + SpectralComplexCoefficients X5_coef, X6_coef; + // These real and complex coefficients are allocated only with averaged Galilean PSATD SpectralComplexCoefficients Psi1_coef, Psi2_coef, Y1_coef, Y2_coef, Y3_coef, Y4_coef; @@ -153,6 +181,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm amrex::Real m_dt; bool m_update_with_rho; bool m_time_averaging; + bool m_J_linear_in_time; bool m_is_galilean; }; #endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index 8a5b791ab..b454d79ba 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -36,7 +36,8 @@ PsatdAlgorithm::PsatdAlgorithm( const amrex::Array<amrex::Real,3>& v_galilean, const amrex::Real dt, const bool update_with_rho, - const bool time_averaging) + const bool time_averaging, + const bool J_linear_in_time) // Initializer list : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal), // Initialize the centered finite-order modified k vectors: @@ -52,7 +53,8 @@ PsatdAlgorithm::PsatdAlgorithm( m_v_galilean(v_galilean), m_dt(dt), m_update_with_rho(update_with_rho), - m_time_averaging(time_averaging) + m_time_averaging(time_averaging), + m_J_linear_in_time(J_linear_in_time) { const amrex::BoxArray& ba = spectral_kspace.spectralspace_ba; @@ -65,6 +67,7 @@ PsatdAlgorithm::PsatdAlgorithm( X2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); X3_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + // Allocate these coefficients only with Galilean PSATD if (m_is_galilean) { X4_coef = SpectralComplexCoefficients(ba, dm, 1, 0); @@ -73,8 +76,8 @@ PsatdAlgorithm::PsatdAlgorithm( InitializeSpectralCoefficients(spectral_kspace, dm, dt); - // Allocate these coefficients only with averaged Galilean PSATD - if (time_averaging) + // Allocate these coefficients only with time averaging + if (time_averaging && !J_linear_in_time) { Psi1_coef = SpectralComplexCoefficients(ba, dm, 1, 0); Psi2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); @@ -82,9 +85,16 @@ PsatdAlgorithm::PsatdAlgorithm( Y3_coef = SpectralComplexCoefficients(ba, dm, 1, 0); Y2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); 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 + else if (time_averaging && J_linear_in_time) + { + X5_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + X6_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + InitializeSpectralCoefficientsAvgLin(spectral_kspace, dm, dt); + } } void @@ -92,8 +102,11 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const { const bool update_with_rho = m_update_with_rho; const bool time_averaging = m_time_averaging; + const bool J_linear_in_time = m_J_linear_in_time; const bool is_galilean = m_is_galilean; + const amrex::Real dt = m_dt; + // Loop over boxes for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi) { @@ -125,7 +138,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const amrex::Array4<const Complex> Y3_arr; amrex::Array4<const Complex> Y4_arr; - if (time_averaging) + if (time_averaging && !J_linear_in_time) { Psi1_arr = Psi1_coef[mfi].array(); Psi2_arr = Psi2_coef[mfi].array(); @@ -135,6 +148,14 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const Y4_arr = Y4_coef[mfi].array(); } + Array4<const Complex> X5_arr; + Array4<const Complex> X6_arr; + if (time_averaging && J_linear_in_time) + { + 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 (AMREX_SPACEDIM == 3) @@ -146,7 +167,8 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { using Idx = SpectralFieldIndex; - using AvgIdx = SpectralAvgFieldIndex; + using IdxAvg = SpectralFieldIndexTimeAveraging; + using IdxLin = SpectralFieldIndexJLinearInTime; // Record old values of the fields to be updated const Complex Ex_old = fields(i,j,k,Idx::Ex); @@ -173,7 +195,9 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const const amrex::Real kz = modified_kz_arr[j]; #endif // Physical constants and imaginary unit - const amrex::Real c2 = std::pow(PhysConst::c, 2); + 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}; // These coefficients are initialized in the function InitializeSpectralCoefficients @@ -239,9 +263,83 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const - I * T2 * S_ck * (kx * Ey_old - ky * Ex_old) + I * X1 * (kx * Jy - ky * Jx); - // Additional update equations for averaged Galilean algorithm + if (J_linear_in_time) + { + const Complex Jx_new = fields(i,j,k,IdxLin::Jx_new); + const Complex Jy_new = fields(i,j,k,IdxLin::Jy_new); + const Complex Jz_new = fields(i,j,k,IdxLin::Jz_new); + + const Complex F_old = fields(i,j,k,IdxLin::F); + const Complex G_old = fields(i,j,k,IdxLin::G); + + fields(i,j,k,Idx::Ex) += -X1 * (Jx_new - Jx) / dt + I * c2 * S_ck * F_old * kx; + + fields(i,j,k,Idx::Ey) += -X1 * (Jy_new - Jy) / dt + I * c2 * S_ck * F_old * ky; + + fields(i,j,k,Idx::Ez) += -X1 * (Jz_new - Jz) / dt + I * c2 * S_ck * F_old * kz; + + fields(i,j,k,Idx::Bx) += I * X2/c2 * (ky * (Jz_new - Jz) - kz * (Jy_new - Jy)); + + I * c2 * S_ck * G_old * kx; + + fields(i,j,k,Idx::By) += I * X2/c2 * (kz * (Jx_new - Jx) - kx * (Jz_new - Jz)); + + I * c2 * S_ck * G_old * ky; + + fields(i,j,k,Idx::Bz) += I * X2/c2 * (kx * (Jy_new - Jy) - ky * (Jx_new - Jx)); + + I * c2 * S_ck * G_old * kz; + + const Complex k_dot_J = kx * Jx + ky * Jy + kz * Jz; + const Complex k_dot_dJ = kx * (Jx_new - Jx) + ky * (Jy_new - Jy) + kz * (Jz_new - Jz); + const Complex k_dot_E = kx * Ex_old + ky * Ey_old + kz * Ez_old; + const Complex k_dot_B = kx * Bx_old + ky * By_old + kz * Bz_old; - if (time_averaging) + fields(i,j,k,IdxLin::F) = C * F_old + S_ck * (I * k_dot_E - rho_old * inv_ep0) + - X1 * ((rho_new - rho_old) / dt + I * k_dot_J) - I * X2/c2 * k_dot_dJ; + + fields(i,j,k,IdxLin::G) = C * G_old + I * S_ck * k_dot_B; + + 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,IdxLin::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 + I * c2 * ep0 * X1 * F_old * kx; + + fields(i,j,k,IdxLin::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 + I * c2 * ep0 * X1 * F_old * ky; + + fields(i,j,k,IdxLin::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 + I * c2 * ep0 * X1 * F_old * kz; + + fields(i,j,k,IdxLin::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); + + I * c2 * ep0 * X1 * G_old * kx; + + fields(i,j,k,IdxLin::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); + + I * c2 * ep0 * X1 * G_old * ky; + + fields(i,j,k,IdxLin::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); + + I * c2 * ep0 * X1 * G_old * kz; + } + } + + // Additional update equations for averaged Galilean algorithm + if (time_averaging && !J_linear_in_time) { // These coefficients are initialized in the function InitializeSpectralCoefficients below const Complex Psi1 = Psi1_arr(i,j,k); @@ -251,27 +349,27 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const const Complex Y2 = Y2_arr(i,j,k); const Complex Y4 = Y4_arr(i,j,k); - fields(i,j,k,AvgIdx::Ex_avg) = Psi1 * Ex_old + fields(i,j,k,IdxAvg::Ex_avg) = Psi1 * Ex_old - I * c2 * Psi2 * (ky * Bz_old - kz * By_old) + Y4 * Jx + (Y2 * rho_new + Y3 * rho_old) * kx; - fields(i,j,k,AvgIdx::Ey_avg) = Psi1 * Ey_old + fields(i,j,k,IdxAvg::Ey_avg) = Psi1 * Ey_old - I * c2 * Psi2 * (kz * Bx_old - kx * Bz_old) + Y4 * Jy + (Y2 * rho_new + Y3 * rho_old) * ky; - fields(i,j,k,AvgIdx::Ez_avg) = Psi1 * Ez_old + fields(i,j,k,IdxAvg::Ez_avg) = Psi1 * Ez_old - I * c2 * Psi2 * (kx * By_old - ky * Bx_old) + Y4 * Jz + (Y2 * rho_new + Y3 * rho_old) * kz; - fields(i,j,k,AvgIdx::Bx_avg) = Psi1 * Bx_old + fields(i,j,k,IdxAvg::Bx_avg) = Psi1 * Bx_old + I * Psi2 * (ky * Ez_old - kz * Ey_old) + I * Y1 * (ky * Jz - kz * Jy); - fields(i,j,k,AvgIdx::By_avg) = Psi1 * By_old + fields(i,j,k,IdxAvg::By_avg) = Psi1 * By_old + I * Psi2 * (kz * Ex_old - kx * Ez_old) + I * Y1 * (kz * Jx - kx * Jz); - fields(i,j,k,AvgIdx::Bz_avg) = Psi1 * Bz_old + fields(i,j,k,IdxAvg::Bz_avg) = Psi1 * Bz_old + I * Psi2 * (kx * Ey_old - ky * Ex_old) + I * Y1 * (kx * Jy - ky * Jx); } @@ -672,6 +770,76 @@ 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 (AMREX_SPACEDIM==3) + 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 (AMREX_SPACEDIM==3) + 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, diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index a1fed5cce..4999a268d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -37,8 +37,14 @@ struct SpectralFieldIndex { enum { Ex=0, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new, n_fields, divE=3 }; }; +struct SpectralFieldIndexJLinearInTime { + enum { Ex=0, Ey, Ez, Bx, By, Bz, Jx_old, Jy_old, Jz_old, rho_old, rho_new, + Ex_avg, Ey_avg, Ez_avg, Bx_avg, By_avg, Bz_avg, + Jx_new, Jy_new, Jz_new, F, G, n_fields, divE=3 }; +}; + /* Index for the regular fields + averaged fields, when stored in spectral space */ -struct SpectralAvgFieldIndex { +struct SpectralFieldIndexTimeAveraging { enum { Ex=0, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new, Ex_avg, Ey_avg, Ez_avg, Bx_avg, By_avg, Bz_avg,n_fields }; // n_fields is automatically the total number of fields }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 996053a88..e5d84a2af 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -32,12 +32,40 @@ class SpectralSolver { public: - // Inline definition of the member functions of `SpectralSolver`, - // except the constructor (see `SpectralSolver.cpp`) - // The body of these functions is short, since the work is done in the - // underlying classes `SpectralFieldData` and `PsatdAlgorithm` - // Constructor + /** + * \brief Constructor of the class SpectralSolver + * + * Select the spectral algorithm to be used, allocate the corresponding coefficients + * for the discrete field update equations, and prepare the structures that store + * the fields in spectral space. + * + * \param[in] lev mesh refinement level + * \param[in] realspace_ba BoxArray in real space + * \param[in] dm DistributionMapping for the given BoxArray + * \param[in] norder_x spectral order along x + * \param[in] norder_y spectral order along y + * \param[in] norder_z spectral order along z + * \param[in] nodal whether the spectral solver is applied to a nodal or staggered grid + * \param[in] v_galilean three-dimensional vector containing the components of the Galilean + * velocity for the standard or averaged Galilean PSATD solvers + * \param[in] v_comoving three-dimensional vector containing the components of the comoving + * velocity for the comoving PSATD solver + * \param[in] dx AMREX_SPACEDIM-dimensional vector containing the cell sizes along each direction + * \param[in] dt time step for the analytical integration of Maxwell's equations + * \param[in] pml whether the boxes in the given BoxArray are PML boxes + * \param[in] periodic_single_box whether there is only one periodic single box + * (no domain decomposition) + * \param[in] update_with_rho whether rho is used in the field update equations + * \param[in] fft_do_time_averaging whether the time averaging algorithm is used + * \param[in] J_linear_in_time whether to use two currents computed at the beginning and + * the end of the time interval (instead of using one current + * compute at half time) + * \param[in] dive_cleaning whether to use div(E) cleaning to account for errors in Gauss law + * (new field F in the field update equations) + * \param[in] divb_cleaning whether to use div(B) cleaning to account for errors in magnetic Gauss law + * (new field G in the field update equations) + */ SpectralSolver (const int lev, const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, @@ -47,10 +75,11 @@ class SpectralSolver 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, - const bool update_with_rho=false, - const bool fft_do_time_averaging=false, + 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 J_linear_in_time = false, const bool dive_cleaning = false, const bool divb_cleaning = false); @@ -117,7 +146,46 @@ class SpectralSolver algorithm->VayDeposition(lev, field_data, current); } + /** + * \brief Copy spectral data from component \c src_comp to component \c dest_comp + * of \c field_data.fields. + * + * \param[in] src_comp component of the source FabArray from which the data are copied + * \param[in] dest_comp component of the destination FabArray where the data are copied + */ + void CopySpectralDataComp (const int src_comp, const int dest_comp) + { + // The last two arguments represent the number of components and + // the number of ghost cells to perform this operation + Copy(field_data.fields, field_data.fields, src_comp, dest_comp, 1, 0); + } + + /** + * \brief Set to zero the data on component \c icomp of \c field_data.fields + * + * \param[in] icomp component of the FabArray where the data are set to zero + */ + void ZeroOutDataComp (const int icomp) + { + // The last argument represents the number of components to perform this operation + field_data.fields.setVal(0., icomp, 1); + } + + /** + * \brief Scale the data on component \c icomp of \c field_data.fields + * by a given scale factor + * + * \param[in] icomp component of the FabArray where the data are scaled + * \param[in] scale_factor scale factor to use for scaling + */ + void ScaleDataComp (const int icomp, const amrex::Real scale_factor) + { + // The last argument represents the number of components to perform this operation + field_data.fields.mult(scale_factor, icomp, 1); + } + private: + void ReadParameters (); // Store field in spectral space and perform the Fourier transforms diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index d04961c5f..89a7ce1f5 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -17,21 +17,6 @@ #if WARPX_USE_PSATD -/* \brief Initialize the spectral Maxwell solver - * - * This function selects the spectral algorithm to be used, allocates the - * corresponding coefficients for the discretized field update equation, - * and prepares the structures that store the fields in spectral space. - * - * \param norder_x Order of accuracy of the spatial derivatives along x - * \param norder_y Order of accuracy of the spatial derivatives along y - * \param norder_z Order of accuracy of the spatial derivatives along z - * \param nodal Whether the solver is applied to a nodal or staggered grid - * \param dx Cell size along each dimension - * \param dt Time step - * \param pml Whether the boxes in which the solver is applied are PML boxes - * \param periodic_single_box Whether the full simulation domain consists of a single periodic box (i.e. the global domain is not MPI parallelized) - */ SpectralSolver::SpectralSolver( const int lev, const amrex::BoxArray& realspace_ba, @@ -44,9 +29,10 @@ SpectralSolver::SpectralSolver( const bool pml, const bool periodic_single_box, const bool update_with_rho, const bool fft_do_time_averaging, + const bool J_linear_in_time, const bool dive_cleaning, - const bool divb_cleaning) { - + const bool divb_cleaning) +{ // Initialize all structures using the same distribution mapping dm // - Initialize k space object (Contains info about the size of @@ -70,14 +56,13 @@ SpectralSolver::SpectralSolver( // PSATD algorithms: standard, Galilean, or averaged Galilean else { algorithm = std::make_unique<PsatdAlgorithm>( - k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt, update_with_rho, fft_do_time_averaging); + k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt, update_with_rho, fft_do_time_averaging, J_linear_in_time); } } // - Initialize arrays for fields in spectral space + FFT plans field_data = SpectralFieldData( lev, realspace_ba, k_space, dm, algorithm->getRequiredNumberOfFields(), periodic_single_box); - } void |