diff options
Diffstat (limited to 'Source/FieldSolver/SpectralSolver')
7 files changed, 632 insertions, 197 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt index 1ce24961d..3fbe99765 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt @@ -1,6 +1,7 @@ target_sources(WarpX PRIVATE PsatdAlgorithm.cpp + PsatdAlgorithmMultiJ.cpp PMLPsatdAlgorithm.cpp SpectralBaseAlgorithm.cpp ComovingPsatdAlgorithm.cpp diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package index ab3256452..e4c32ea2b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -1,5 +1,6 @@ CEXE_sources += SpectralBaseAlgorithm.cpp CEXE_sources += PsatdAlgorithm.cpp +CEXE_sources += PsatdAlgorithmMultiJ.cpp CEXE_sources += PMLPsatdAlgorithm.cpp CEXE_sources += ComovingPsatdAlgorithm.cpp diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index 5549eda07..f9a155742 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -43,9 +43,6 @@ 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] do_multi_J whether the multi-J algorithm is used (hence two currents - * computed at the beginning and the end of the time interval - * instead of one current computed at half time) * \param[in] dive_cleaning Update F as part of the field update, so that errors in divE=rho propagate away at the speed of light * \param[in] divb_cleaning Update G as part of the field update, so that errors in divB=0 propagate away at the speed of light */ @@ -62,7 +59,6 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm 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); @@ -86,19 +82,6 @@ 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 * @@ -153,8 +136,6 @@ 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; @@ -172,7 +153,6 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm amrex::Real m_dt; bool m_update_with_rho; bool m_time_averaging; - bool m_do_multi_J; bool m_dive_cleaning; bool m_divb_cleaning; bool m_is_galilean; 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, diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmMultiJ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmMultiJ.H new file mode 100644 index 000000000..215e35676 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmMultiJ.H @@ -0,0 +1,147 @@ +/* Copyright 2019 + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PSATD_ALGORITHM_MULTIJ_H_ +#define WARPX_PSATD_ALGORITHM_MULTIJ_H_ + +#include "FieldSolver/SpectralSolver/SpectralFieldData.H" +#include "FieldSolver/SpectralSolver/SpectralKSpace.H" +#include "SpectralBaseAlgorithm.H" + +#include <AMReX_Array.H> +#include <AMReX_Config.H> +#include <AMReX_REAL.H> + +#include <AMReX_BaseFwd.H> + +#include <array> +#include <memory> + +#if WARPX_USE_PSATD +/* + * \brief Class that updates the fields in spectral space according to the multi-J algorithm + * and stores the coefficients of the corresponding update equations. J is assumed to be + * linear in time and two currents, deposited at the beginning and the end of the time step, + * are used for the PSATD update equations, instead of only one current deposited at half time. + */ +class PsatdAlgorithmMultiJ : public SpectralBaseAlgorithm +{ + public: + + /** + * \brief Constructor of the class PsatdAlgorithmMultiJ + * + * \param[in] spectral_kspace spectral space + * \param[in] dm distribution mapping + * \param[in] spectral_index object containing indices to access data in spectral space + * \param[in] norder_x order of the spectral solver along x + * \param[in] norder_y order of the spectral solver along y + * \param[in] norder_z order of the spectral solver along z + * \param[in] nodal whether the E and B fields are defined on a fully nodal grid or a Yee grid + * \param[in] fill_guards Update the guard cells (in addition to the valid cells) when pushing the fields in time + * \param[in] dt time step of the simulation + * \param[in] time_averaging whether to use time averaging for large time steps + * \param[in] dive_cleaning Update F as part of the field update, so that errors in divE=rho propagate away at the speed of light + * \param[in] divb_cleaning Update G as part of the field update, so that errors in divB=0 propagate away at the speed of light + */ + PsatdAlgorithmMultiJ ( + const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const SpectralFieldIndex& spectral_index, + const int norder_x, + const int norder_y, + const int norder_z, + const bool nodal, + const amrex::IntVect& fill_guards, + const amrex::Real dt, + const bool time_averaging, + const bool dive_cleaning, + const bool divb_cleaning); + + /** + * \brief Updates the E and B fields in spectral space, according to the multi-J PSATD equations + * + * \param[in,out] f all the fields in spectral space + */ + virtual void pushSpectralFields (SpectralFieldData& f) const override final; + + /** + * \brief Initializes the coefficients used in \c pushSpectralFields to update the E and B fields + * + * \param[in] spectral_kspace spectral space + * \param[in] dm distribution mapping + * \param[in] dt time step of the simulation + */ + void InitializeSpectralCoefficients ( + const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + 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 InitializeSpectralCoefficientsAveraging ( + const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const amrex::Real dt); + + /** + * \brief Virtual function for current correction in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This function overrides the virtual function \c CurrentCorrection in the + * base class \c SpectralBaseAlgorithm and cannot be overridden by further + * derived classes. + * + * \param[in] lev The mesh-refinement level + * \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 ( + const int lev, + 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 + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This function overrides the virtual function \c VayDeposition in the + * base class \c SpectralBaseAlgorithm and cannot be overridden by further + * derived classes. + * + * \param[in] lev The mesh-refinement level + * \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 ( + const int lev, + SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + + private: + + // These real and complex coefficients are always allocated + SpectralRealCoefficients C_coef, S_ck_coef; + SpectralRealCoefficients X1_coef, X2_coef, X3_coef, X5_coef, X6_coef; + + SpectralFieldIndex m_spectral_index; + + // Other member variables + amrex::Real m_dt; + bool m_time_averaging; + bool m_dive_cleaning; + bool m_divb_cleaning; +}; +#endif // WARPX_USE_PSATD +#endif // WARPX_PSATD_ALGORITHM_MULTIJ_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmMultiJ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmMultiJ.cpp new file mode 100644 index 000000000..7ecde8d79 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmMultiJ.cpp @@ -0,0 +1,453 @@ +/* Copyright 2019 + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "PsatdAlgorithmMultiJ.H" + +#include "Utils/WarpXConst.H" +#include "Utils/WarpX_Complex.H" + +#include <AMReX_Array4.H> +#include <AMReX_BLProfiler.H> +#include <AMReX_BaseFab.H> +#include <AMReX_BoxArray.H> +#include <AMReX_GpuComplex.H> +#include <AMReX_GpuLaunch.H> +#include <AMReX_GpuQualifiers.H> +#include <AMReX_IntVect.H> +#include <AMReX_MFIter.H> +#include <AMReX_PODVector.H> + +#include <cmath> + +#if WARPX_USE_PSATD + +using namespace amrex::literals; + +PsatdAlgorithmMultiJ::PsatdAlgorithmMultiJ( + const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const SpectralFieldIndex& spectral_index, + const int norder_x, + const int norder_y, + const int norder_z, + const bool nodal, + const amrex::IntVect& fill_guards, + const amrex::Real dt, + const bool time_averaging, + const bool dive_cleaning, + const bool divb_cleaning) + // Initializer list + : SpectralBaseAlgorithm(spectral_kspace, dm, spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards), + m_spectral_index(spectral_index), + m_dt(dt), + m_time_averaging(time_averaging), + m_dive_cleaning(dive_cleaning), + m_divb_cleaning(divb_cleaning) +{ + const amrex::BoxArray& ba = spectral_kspace.spectralspace_ba; + + // Always allocate these coefficients + C_coef = SpectralRealCoefficients(ba, dm, 1, 0); + S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0); + X1_coef = SpectralRealCoefficients(ba, dm, 1, 0); + X2_coef = SpectralRealCoefficients(ba, dm, 1, 0); + X3_coef = SpectralRealCoefficients(ba, dm, 1, 0); + + InitializeSpectralCoefficients(spectral_kspace, dm, dt); + + // Allocate these coefficients only with time averaging + if (time_averaging) + { + X5_coef = SpectralRealCoefficients(ba, dm, 1, 0); + X6_coef = SpectralRealCoefficients(ba, dm, 1, 0); + InitializeSpectralCoefficientsAveraging(spectral_kspace, dm, dt); + } +} + +void +PsatdAlgorithmMultiJ::pushSpectralFields (SpectralFieldData& f) const +{ + const bool time_averaging = m_time_averaging; + const bool dive_cleaning = m_dive_cleaning; + const bool divb_cleaning = m_divb_cleaning; + + const amrex::Real dt = m_dt; + + const SpectralFieldIndex& Idx = m_spectral_index; + + // 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(); + + // These coefficients are always allocated + 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 amrex::Real> X1_arr = X1_coef[mfi].array(); + amrex::Array4<const amrex::Real> X2_arr = X2_coef[mfi].array(); + amrex::Array4<const amrex::Real> X3_arr = X3_coef[mfi].array(); + + amrex::Array4<const amrex::Real> X5_arr; + amrex::Array4<const amrex::Real> X6_arr; + if (time_averaging) + { + 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) + 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 + 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_old = fields(i,j,k,Idx.Jx); + const Complex Jy_old = fields(i,j,k,Idx.Jy); + const Complex Jz_old = fields(i,j,k,Idx.Jz); + 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); + const Complex rho_old = fields(i,j,k,Idx.rho_old); + const Complex rho_new = fields(i,j,k,Idx.rho_new); + + Complex F_old, G_old; + if (dive_cleaning) F_old = fields(i,j,k,Idx.F); + if (divb_cleaning) G_old = fields(i,j,k,Idx.G); + + // k vector values + const amrex::Real kx = modified_kx_arr[i]; +#if defined(WARPX_DIM_3D) + const amrex::Real ky = modified_ky_arr[j]; + const amrex::Real kz = modified_kz_arr[k]; +#else + constexpr amrex::Real ky = 0._rt; + const amrex::Real kz = modified_kz_arr[j]; +#endif + // Physical constants and imaginary unit + constexpr amrex::Real c2 = PhysConst::c * PhysConst::c; + constexpr amrex::Real ep0 = PhysConst::ep0; + constexpr amrex::Real inv_ep0 = 1._rt / PhysConst::ep0; + constexpr Complex I = Complex{0._rt, 1._rt}; + + // These coefficients are initialized in the function InitializeSpectralCoefficients + const amrex::Real C = C_arr(i,j,k); + const amrex::Real S_ck = S_ck_arr(i,j,k); + const amrex::Real X1 = X1_arr(i,j,k); + const amrex::Real X2 = X2_arr(i,j,k); + const amrex::Real X3 = X3_arr(i,j,k); + const amrex::Real X4 = - S_ck / PhysConst::ep0; + + // Update equations for E in the formulation with rho + + fields(i,j,k,Idx.Ex) = C * Ex_old + + I * c2 * S_ck * (ky * Bz_old - kz * By_old) + + X4 * Jx_old - I * (X2 * rho_new - X3 * rho_old) * kx - X1 * (Jx_new - Jx_old) / dt; + + fields(i,j,k,Idx.Ey) = C * Ey_old + + I * c2 * S_ck * (kz * Bx_old - kx * Bz_old) + + X4 * Jy_old - I * (X2 * rho_new - X3 * rho_old) * ky - X1 * (Jy_new - Jy_old) / dt; + + fields(i,j,k,Idx.Ez) = C * Ez_old + + I * c2 * S_ck * (kx * By_old - ky * Bx_old) + + X4 * Jz_old - I * (X2 * rho_new - X3 * rho_old) * kz - X1 * (Jz_new - Jz_old) / dt; + + // Update equations for B + + fields(i,j,k,Idx.Bx) = C * Bx_old + - I * S_ck * (ky * Ez_old - kz * Ey_old) + I * X1 * (ky * Jz_old - kz * Jy_old) + + I * X2/c2 * (ky * (Jz_new - Jz_old) - kz * (Jy_new - Jy_old)); + + fields(i,j,k,Idx.By) = C * By_old + - I * S_ck * (kz * Ex_old - kx * Ez_old) + I * X1 * (kz * Jx_old - kx * Jz_old) + + I * X2/c2 * (kz * (Jx_new - Jx_old) - kx * (Jz_new - Jz_old)); + + fields(i,j,k,Idx.Bz) = C * Bz_old + - I * S_ck * (kx * Ey_old - ky * Ex_old) + I * X1 * (kx * Jy_old - ky * Jx_old) + + I * X2/c2 * (kx * (Jy_new - Jy_old) - ky * (Jx_new - Jx_old)); + + if (dive_cleaning) + { + const Complex k_dot_E = kx * Ex_old + ky * Ey_old + kz * Ez_old; + const Complex k_dot_J = kx * Jx_old + ky * Jy_old + kz * Jz_old; + const Complex k_dot_dJ = kx * (Jx_new - Jx_old) + ky * (Jy_new - Jy_old) + kz * (Jz_new - Jz_old); + + fields(i,j,k,Idx.Ex) += I * c2 * S_ck * F_old * kx; + fields(i,j,k,Idx.Ey) += I * c2 * S_ck * F_old * ky; + fields(i,j,k,Idx.Ez) += I * c2 * S_ck * F_old * kz; + + fields(i,j,k,Idx.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; + } + + if (divb_cleaning) + { + const Complex k_dot_B = kx * Bx_old + ky * By_old + kz * Bz_old; + + fields(i,j,k,Idx.Bx) += I * S_ck * G_old * kx; + fields(i,j,k,Idx.By) += I * S_ck * G_old * ky; + fields(i,j,k,Idx.Bz) += I * S_ck * G_old * kz; + + fields(i,j,k,Idx.G) = C * G_old + I * c2 * S_ck * k_dot_B; + } + + if (time_averaging) + { + const amrex::Real X5 = X5_arr(i,j,k); + const amrex::Real 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_old - 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_old - 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_old - 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_old - kz * Jy_old) - 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_old - kx * Jz_old) - 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_old - ky * Jx_old) - 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; + } + } + }); + } +} + +void PsatdAlgorithmMultiJ::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_s = modified_kx_vec[mfi].dataPtr(); +#if defined(WARPX_DIM_3D) + const amrex::Real* ky_s = modified_ky_vec[mfi].dataPtr(); +#endif + const amrex::Real* kz_s = modified_kz_vec[mfi].dataPtr(); + + // Coefficients always allocated + amrex::Array4<amrex::Real> C = C_coef[mfi].array(); + amrex::Array4<amrex::Real> S_ck = S_ck_coef[mfi].array(); + amrex::Array4<amrex::Real> X1 = X1_coef[mfi].array(); + amrex::Array4<amrex::Real> X2 = X2_coef[mfi].array(); + amrex::Array4<amrex::Real> X3 = X3_coef[mfi].array(); + + // Loop over indices within one box + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Calculate norm of k vector + const amrex::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 amrex::Real c = PhysConst::c; + constexpr amrex::Real ep0 = PhysConst::ep0; + + const amrex::Real c2 = std::pow(c, 2); + const amrex::Real dt2 = std::pow(dt, 2); + + const amrex::Real om_s = c * knorm_s; + const amrex::Real om2_s = std::pow(om_s, 2); + + // C + C(i,j,k) = std::cos(om_s * dt); + + // S_ck + if (om_s != 0.) + { + S_ck(i,j,k) = std::sin(om_s * dt) / om_s; + } + else // om_s = 0 + { + S_ck(i,j,k) = dt; + } + + // X1 (multiplies i*([k] \times J) in the update equation for update B) + if (om_s != 0.) + { + X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2_s); + } + else // om_s = 0 + { + X1(i,j,k) = 0.5_rt * dt2 / ep0; + } + + // X2 (multiplies rho_new in the update equation for E) + if (om_s != 0.) + { + X2(i,j,k) = c2 * (dt - S_ck(i,j,k)) / (ep0 * dt * om2_s); + } + else // om_s = 0 + { + X2(i,j,k) = c2 * dt2 / (6._rt * ep0); + } + + // X3 (multiplies rho_old in the update equation for E) + if (om_s != 0.) + { + X3(i,j,k) = c2 * (dt * C(i,j,k) - S_ck(i,j,k)) / (ep0 * dt * om2_s); + } + else // om_s = 0 + { + X3(i,j,k) = - c2 * dt2 / (3._rt * ep0); + } + }); + } +} + +void PsatdAlgorithmMultiJ::InitializeSpectralCoefficientsAveraging ( + 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_s = modified_kx_vec[mfi].dataPtr(); +#if defined(WARPX_DIM_3D) + const amrex::Real* ky_s = modified_ky_vec[mfi].dataPtr(); +#endif + const amrex::Real* kz_s = modified_kz_vec[mfi].dataPtr(); + + amrex::Array4<amrex::Real const> C = C_coef[mfi].array(); + amrex::Array4<amrex::Real const> S_ck = S_ck_coef[mfi].array(); + + amrex::Array4<amrex::Real> X5 = X5_coef[mfi].array(); + amrex::Array4<amrex::Real> X6 = X6_coef[mfi].array(); + + // Loop over indices within one box + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Calculate norm of k vector + const amrex::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 amrex::Real c = PhysConst::c; + constexpr amrex::Real c2 = c*c; + constexpr amrex::Real ep0 = PhysConst::ep0; + + // Auxiliary coefficients + const amrex::Real dt3 = dt * dt * dt; + + const amrex::Real om_s = c * knorm_s; + const amrex::Real om2_s = om_s * om_s; + const amrex::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 +PsatdAlgorithmMultiJ::CurrentCorrection ( + const int lev, + SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current, + const std::unique_ptr<amrex::MultiFab>& rho) +{ + // Profiling + BL_PROFILE("PsatdAlgorithmMultiJ::CurrentCorrection"); + + amrex::ignore_unused(lev, field_data, current, rho); + amrex::Abort("Current correction not implemented for multi-J PSATD algorithm"); +} + +void +PsatdAlgorithmMultiJ::VayDeposition ( + const int lev, + SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) +{ + // Profiling + BL_PROFILE("PsatdAlgorithmMultiJ::VayDeposition()"); + + amrex::ignore_unused(lev, field_data, current); + amrex::Abort("Vay deposition not implemented for multi-J PSATD algorithm"); +} + +#endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index 70800f732..ee169398c 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -9,6 +9,7 @@ #include "SpectralAlgorithms/ComovingPsatdAlgorithm.H" #include "SpectralAlgorithms/PMLPsatdAlgorithm.H" #include "SpectralAlgorithms/PsatdAlgorithm.H" +#include "SpectralAlgorithms/PsatdAlgorithmMultiJ.H" #include "SpectralKSpace.H" #include "SpectralSolver.H" #include "Utils/WarpXProfilerWrapper.H" @@ -47,24 +48,36 @@ SpectralSolver::SpectralSolver( // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space - if (pml) { + if (pml) // PSATD equations in the PML grids + { algorithm = std::make_unique<PMLPsatdAlgorithm>( k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards, dt, dive_cleaning, divb_cleaning); } - else { + else // PSATD equations in the regulard grids + { // Comoving PSATD algorithm - if (v_comoving[0] != 0. || v_comoving[1] != 0. || v_comoving[2] != 0.) { + if (v_comoving[0] != 0. || v_comoving[1] != 0. || v_comoving[2] != 0.) + { algorithm = std::make_unique<ComovingPsatdAlgorithm>( k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards, v_comoving, dt, update_with_rho); } - // PSATD algorithms: standard, Galilean, or averaged Galilean - else { - algorithm = std::make_unique<PsatdAlgorithm>( - k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards, - v_galilean, dt, update_with_rho, fft_do_time_averaging, do_multi_J, - dive_cleaning, divb_cleaning); + else // PSATD algorithms: standard, Galilean, averaged Galilean, multi-J + { + if (do_multi_J) + { + algorithm = std::make_unique<PsatdAlgorithmMultiJ>( + k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, + fill_guards, dt, fft_do_time_averaging, dive_cleaning, divb_cleaning); + } + else // standard, Galilean, averaged Galilean + { + algorithm = std::make_unique<PsatdAlgorithm>( + k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, + fill_guards, v_galilean, dt, update_with_rho, fft_do_time_averaging, + dive_cleaning, divb_cleaning); + } } } |