diff options
Diffstat (limited to 'Source/FieldSolver/SpectralSolver')
8 files changed, 285 insertions, 10 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt index f46fd9eb9..1ce24961d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt @@ -12,5 +12,6 @@ if(WarpX_DIMS STREQUAL RZ) SpectralBaseAlgorithmRZ.cpp PsatdAlgorithmRZ.cpp GalileanPsatdAlgorithmRZ.cpp + PMLPsatdAlgorithmRZ.cpp ) endif() diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package index 11a10bb91..ab3256452 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -7,6 +7,7 @@ ifeq ($(USE_RZ),TRUE) CEXE_sources += SpectralBaseAlgorithmRZ.cpp CEXE_sources += PsatdAlgorithmRZ.cpp CEXE_sources += GalileanPsatdAlgorithmRZ.cpp + CEXE_sources += PMLPsatdAlgorithmRZ.cpp endif VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.H new file mode 100644 index 000000000..e0956a514 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.H @@ -0,0 +1,72 @@ +/* Copyright 2021 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PMLPSATD_ALGORITHM_RZ_H_ +#define WARPX_PMLPSATD_ALGORITHM_RZ_H_ + +#include "SpectralBaseAlgorithmRZ.H" + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + */ +class PMLPsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ +{ + + public: + PMLPsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, + amrex::DistributionMapping const & dm, + const SpectralFieldIndex& spectral_index, + int const n_rz_azimuthal_modes, int const norder_z, + bool const nodal, amrex::Real const dt_step); + + // Redefine functions from base class + virtual void pushSpectralFields (SpectralFieldDataRZ & f) override final; + + void InitializeSpectralCoefficients (SpectralFieldDataRZ const & f); + + /** + * \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 SpectralBaseAlgorithmRZ and cannot be overridden by further + * derived classes. + * + * \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, + SpectralFieldDataRZ& 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 SpectralBaseAlgorithmRZ and cannot be overridden by further + * derived classes. + * + * \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, + SpectralFieldDataRZ& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + + private: + + SpectralFieldIndex m_spectral_index; + + bool coefficients_initialized; + // Note that dt is saved to use in InitializeSpectralCoefficients + amrex::Real m_dt; + SpectralRealCoefficients C_coef, S_ck_coef; +}; + +#endif // WARPX_PMLPSATD_ALGORITHM_RZ_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.cpp new file mode 100644 index 000000000..892c542c7 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.cpp @@ -0,0 +1,177 @@ +/* Copyright 2021 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "PMLPsatdAlgorithmRZ.H" +#include "FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.H" +#include "Utils/WarpXConst.H" +#include "Utils/WarpXProfilerWrapper.H" +#include "WarpX.H" + +#include <cmath> + +using amrex::operator""_rt; + + +/* \brief Initialize coefficients for the update equation */ +PMLPsatdAlgorithmRZ::PMLPsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, + amrex::DistributionMapping const & dm, + const SpectralFieldIndex& spectral_index, + int const n_rz_azimuthal_modes, int const norder_z, + bool const nodal, amrex::Real const dt) + // Initialize members of base class + : SpectralBaseAlgorithmRZ(spectral_kspace, dm, spectral_index, norder_z, nodal), + m_spectral_index(spectral_index), + m_dt(dt) +{ + // Allocate the arrays of coefficients + amrex::BoxArray const & ba = spectral_kspace.spectralspace_ba; + C_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0); + S_ck_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0); + + coefficients_initialized = false; +} + +/* Advance the E and B field in spectral space (stored in `f`) + * over one time step */ +void +PMLPsatdAlgorithmRZ::pushSpectralFields (SpectralFieldDataRZ & f) +{ + + if (not coefficients_initialized) { + // This is called from here since it needs the kr values + // which can be obtained from the SpectralFieldDataRZ + InitializeSpectralCoefficients(f); + coefficients_initialized = true; + } + + const SpectralFieldIndex& Idx = m_spectral_index; + + // Loop over boxes + for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi){ + + amrex::Box const & bx = f.fields[mfi].box(); + + // Extract arrays for the fields to be updated + amrex::Array4<Complex> const& fields = f.fields[mfi].array(); + // Extract arrays for the coefficients + amrex::Array4<const amrex::Real> const& C_arr = C_coef[mfi].array(); + amrex::Array4<const amrex::Real> const& S_ck_arr = S_ck_coef[mfi].array(); + + // Extract pointers for the k vectors + HankelTransform::RealVector const & kr_modes = f.getKrArray(mfi); + amrex::Real const* kr_arr = kr_modes.dataPtr(); + int const nr = bx.length(0); + + // Loop over indices within one box + // Note that k = 0 + int const modes = f.n_rz_azimuthal_modes; + amrex::ParallelFor(bx, modes, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept + { + + // All of the fields of each mode are grouped together + int const Ep_m_pml = Idx.Er_pml + Idx.n_fields*mode; + int const Em_m_pml = Idx.Et_pml + Idx.n_fields*mode; + int const Bp_m_pml = Idx.Br_pml + Idx.n_fields*mode; + int const Bm_m_pml = Idx.Bt_pml + Idx.n_fields*mode; + int const Ez_m = Idx.Ez + Idx.n_fields*mode; + int const Bz_m = Idx.Bz + Idx.n_fields*mode; + + // Record old values of the fields to be updated + Complex const Ep_old_pml = fields(i,j,k,Ep_m_pml); + Complex const Em_old_pml = fields(i,j,k,Em_m_pml); + Complex const Bp_old_pml = fields(i,j,k,Bp_m_pml); + Complex const Bm_old_pml = fields(i,j,k,Bm_m_pml); + Complex const Ez_old = fields(i,j,k,Ez_m); + Complex const Bz_old = fields(i,j,k,Bz_m); + + // k vector values, and coefficients + // The k values for each mode are grouped together + int const ir = i + nr*mode; + amrex::Real const kr = kr_arr[ir]; + + constexpr amrex::Real c2 = PhysConst::c*PhysConst::c; + Complex const I = Complex{0._rt,1._rt}; + amrex::Real const C = C_arr(i,j,k,mode); + amrex::Real const S_ck = S_ck_arr(i,j,k,mode); + + // Update E (see WarpX online documentation: theory section) + fields(i,j,k,Ep_m_pml) = C*Ep_old_pml + + S_ck*(-c2*I*kr/2._rt*Bz_old); + fields(i,j,k,Em_m_pml) = C*Em_old_pml + + S_ck*(-c2*I*kr/2._rt*Bz_old); + // Update B (see WarpX online documentation: theory section) + fields(i,j,k,Bp_m_pml) = C*Bp_old_pml + - S_ck*(-I*kr/2._rt*Ez_old); + fields(i,j,k,Bm_m_pml) = C*Bm_old_pml + - S_ck*(-I*kr/2._rt*Ez_old); + + }); + } +} + +void PMLPsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const & f) +{ + + // Fill them with the right values: + // Loop over boxes and allocate the corresponding coefficients + // for each box owned by the local MPI proc + for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi){ + + amrex::Box const & bx = f.fields[mfi].box(); + + // Extract pointers for the k vectors + amrex::Real const* const modified_kz = modified_kz_vec[mfi].dataPtr(); + + // Extract arrays for the coefficients + amrex::Array4<amrex::Real> const& C = C_coef[mfi].array(); + amrex::Array4<amrex::Real> const& S_ck = S_ck_coef[mfi].array(); + + HankelTransform::RealVector const & kr_modes = f.getKrArray(mfi); + amrex::Real const* kr_arr = kr_modes.dataPtr(); + int const nr = bx.length(0); + amrex::Real const dt = m_dt; + + // Loop over indices within one box + int const modes = f.n_rz_azimuthal_modes; + amrex::ParallelFor(bx, modes, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept + { + // Calculate norm of vector + int const ir = i + nr*mode; + amrex::Real const kr = kr_arr[ir]; + amrex::Real const kz = modified_kz[j]; + amrex::Real const k_norm = std::sqrt(kr*kr + kz*kz); + + // Calculate coefficients + constexpr amrex::Real c = PhysConst::c; + if (k_norm != 0._rt){ + C(i,j,k,mode) = std::cos(c*k_norm*dt); + S_ck(i,j,k,mode) = std::sin(c*k_norm*dt)/(c*k_norm); + } else { // Handle k_norm = 0, by using the analytical limit + C(i,j,k,mode) = 1._rt; + S_ck(i,j,k,mode) = dt; + } + }); + } +} + +void +PMLPsatdAlgorithmRZ::CurrentCorrection (const int /* lev */, + SpectralFieldDataRZ& /* field_data */, + std::array<std::unique_ptr<amrex::MultiFab>,3>& /* current */, + const std::unique_ptr<amrex::MultiFab>& /* rho */) +{ + amrex::Abort("Current correction not implemented in RZ geometry PML"); +} + +void +PMLPsatdAlgorithmRZ::VayDeposition (const int /* lev */, + SpectralFieldDataRZ& /*field_data*/, + std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/) +{ + amrex::Abort("Vay deposition not implemented in RZ geometry PML"); +} diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 9b748a048..4f0ad1c20 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -50,13 +50,16 @@ class SpectralFieldIndex * div(B) = 0 law (new field G in the update equations) * \param[in] pml whether the indices are used to access spectral data * for the PML spectral solver + * \param[in] pml_rz whether the indices are used to access spectral data + * for the RZ PML spectral solver */ SpectralFieldIndex (const bool update_with_rho, const bool time_averaging, const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning, - const bool pml); + const bool pml, + const bool pml_rz = false); /** * \brief Default constructor @@ -97,6 +100,9 @@ class SpectralFieldIndex // PML with div(E) and/or div(B) cleaning int Exx = -1, Eyy = -1, Ezz = -1, Bxx = -1, Byy = -1, Bzz = -1; int Fx = -1, Fy = -1, Fz = -1, Gx = -1, Gy = -1, Gz = -1; + + // PML RZ + int Er_pml = -1, Et_pml = -1, Br_pml = -1, Bt_pml = -1; }; /** \brief Class that stores the fields in spectral space, and performs the diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 6d7d18b5f..3926df286 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -38,7 +38,8 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning, - const bool pml) + const bool pml, + const bool pml_rz) { // TODO Use these to allocate rho_old, rho_new, F, and G only when needed amrex::ignore_unused(update_with_rho); @@ -73,6 +74,13 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, Jy_new = c++; Jz_new = c++; } + + if (pml_rz) + { + Er_pml = c++; Et_pml = c++; + Br_pml = c++; Bt_pml = c++; + } + } else // PML { diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H index 827d226c2..9f9b8dbef 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H @@ -34,6 +34,7 @@ class SpectralSolverRZ int const norder_z, bool const nodal, const amrex::Vector<amrex::Real>& v_galilean, amrex::RealVect const dx, amrex::Real const dt, + bool const with_pml, bool const update_with_rho, const bool fft_do_time_averaging, const bool do_multi_J, @@ -63,7 +64,7 @@ class SpectralSolverRZ amrex::MultiFab& field_mf2, int const field_index2); /* \brief Update the fields in spectral space, over one timestep */ - void pushSpectralFields (); + void pushSpectralFields (const bool doing_pml=false); /* \brief Initialize K space filtering arrays */ void InitFilter (amrex::IntVect const & filter_npass_each_dir, @@ -160,6 +161,7 @@ class SpectralSolverRZ SpectralFieldDataRZ field_data; // Store field in spectral space // and perform the Fourier transforms std::unique_ptr<SpectralBaseAlgorithmRZ> algorithm; + std::unique_ptr<SpectralBaseAlgorithmRZ> PML_algorithm; // Defines field update equation in spectral space, // and the associated coefficients. // SpectralBaseAlgorithmRZ is a base class ; this pointer is meant diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp index e187b9d27..558d4a78a 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp @@ -6,6 +6,7 @@ */ #include "SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H" #include "SpectralAlgorithms/PsatdAlgorithmRZ.H" +#include "SpectralAlgorithms/PMLPsatdAlgorithmRZ.H" #include "SpectralKSpaceRZ.H" #include "SpectralSolverRZ.H" #include "Utils/WarpXProfilerWrapper.H" @@ -22,8 +23,7 @@ * \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 - * PML is not supported. + * \param with_pml Whether PML boundary will be used */ SpectralSolverRZ::SpectralSolverRZ (const int lev, amrex::BoxArray const & realspace_ba, @@ -32,6 +32,7 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, int const norder_z, bool const nodal, const amrex::Vector<amrex::Real>& v_galilean, amrex::RealVect const dx, amrex::Real const dt, + bool const with_pml, bool const update_with_rho, const bool fft_do_time_averaging, const bool do_multi_J, @@ -45,13 +46,16 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, // the spectral space corresponding to each box in `realspace_ba`, // as well as the value of the corresponding k coordinates. - const bool pml = false; + const bool is_pml = false; m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging, - do_multi_J, dive_cleaning, divb_cleaning, pml); + do_multi_J, dive_cleaning, divb_cleaning, is_pml, with_pml); // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space - // PML is not supported. + if (with_pml) { + PML_algorithm = std::make_unique<PMLPsatdAlgorithmRZ>( + k_space, dm, m_spectral_index, n_rz_azimuthal_modes, norder_z, nodal, dt); + } if (v_galilean[2] == 0) { // v_galilean is 0: use standard PSATD algorithm algorithm = std::make_unique<PsatdAlgorithmRZ>( @@ -117,12 +121,16 @@ SpectralSolverRZ::BackwardTransform (const int lev, /* \brief Update the fields in spectral space, over one timestep */ void -SpectralSolverRZ::pushSpectralFields () { +SpectralSolverRZ::pushSpectralFields (const bool doing_pml) { WARPX_PROFILE("SpectralSolverRZ::pushSpectralFields"); // Virtual function: the actual function used here depends // on the sub-class of `SpectralBaseAlgorithm` that was // initialized in the constructor of `SpectralSolverRZ` - algorithm->pushSpectralFields(field_data); + if (doing_pml) { + PML_algorithm->pushSpectralFields(field_data); + } else { + algorithm->pushSpectralFields(field_data); + } } /** |