From ec072594fb1bddb4631c55fb3018050cbf461243 Mon Sep 17 00:00:00 2001 From: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> Date: Wed, 2 Feb 2022 16:29:52 -0800 Subject: Rename PSATD Classes (#2805) * Rename PSATD Classes * Rename PsatdAlgorithmJLinear as PsatdAlgorithmJLinearInTime --- .../SpectralAlgorithms/CMakeLists.txt | 10 +- .../SpectralAlgorithms/ComovingPsatdAlgorithm.H | 105 ---- .../SpectralAlgorithms/ComovingPsatdAlgorithm.cpp | 526 --------------------- .../SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H | 78 --- .../GalileanPsatdAlgorithmRZ.cpp | 385 --------------- .../SpectralSolver/SpectralAlgorithms/Make.package | 10 +- .../SpectralAlgorithms/PMLPsatdAlgorithm.H | 91 ---- .../SpectralAlgorithms/PMLPsatdAlgorithm.cpp | 421 ----------------- .../SpectralAlgorithms/PMLPsatdAlgorithmRZ.H | 72 --- .../SpectralAlgorithms/PMLPsatdAlgorithmRZ.cpp | 177 ------- .../SpectralAlgorithms/PsatdAlgorithmComoving.H | 105 ++++ .../SpectralAlgorithms/PsatdAlgorithmComoving.cpp | 526 +++++++++++++++++++++ .../SpectralAlgorithms/PsatdAlgorithmGalileanRZ.H | 78 +++ .../PsatdAlgorithmGalileanRZ.cpp | 385 +++++++++++++++ .../PsatdAlgorithmJLinearInTime.H | 147 ++++++ .../PsatdAlgorithmJLinearInTime.cpp | 453 ++++++++++++++++++ .../SpectralAlgorithms/PsatdAlgorithmMultiJ.H | 147 ------ .../SpectralAlgorithms/PsatdAlgorithmMultiJ.cpp | 453 ------------------ .../SpectralAlgorithms/PsatdAlgorithmPml.H | 91 ++++ .../SpectralAlgorithms/PsatdAlgorithmPml.cpp | 421 +++++++++++++++++ .../SpectralAlgorithms/PsatdAlgorithmPmlRZ.H | 72 +++ .../SpectralAlgorithms/PsatdAlgorithmPmlRZ.cpp | 177 +++++++ 22 files changed, 2465 insertions(+), 2465 deletions(-) delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.cpp delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.H delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.cpp create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.H create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.H create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmMultiJ.H delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmMultiJ.cpp create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.H create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.H create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.cpp (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms') diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt index 3fbe99765..eb88b1f4e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt @@ -1,10 +1,10 @@ target_sources(WarpX PRIVATE PsatdAlgorithm.cpp - PsatdAlgorithmMultiJ.cpp - PMLPsatdAlgorithm.cpp + PsatdAlgorithmJLinearInTime.cpp + PsatdAlgorithmPml.cpp SpectralBaseAlgorithm.cpp - ComovingPsatdAlgorithm.cpp + PsatdAlgorithmComoving.cpp ) if(WarpX_DIMS STREQUAL RZ) @@ -12,7 +12,7 @@ if(WarpX_DIMS STREQUAL RZ) PRIVATE SpectralBaseAlgorithmRZ.cpp PsatdAlgorithmRZ.cpp - GalileanPsatdAlgorithmRZ.cpp - PMLPsatdAlgorithmRZ.cpp + PsatdAlgorithmGalileanRZ.cpp + PsatdAlgorithmPmlRZ.cpp ) endif() diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H deleted file mode 100644 index a1e4a4412..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H +++ /dev/null @@ -1,105 +0,0 @@ -#ifndef WARPX_COMOVING_PSATD_ALGORITHM_H_ -#define WARPX_COMOVING_PSATD_ALGORITHM_H_ - -#include "FieldSolver/SpectralSolver/SpectralFieldData.H" -#include "FieldSolver/SpectralSolver/SpectralKSpace.H" -#include "SpectralBaseAlgorithm.H" - -#include -#include -#include - -#include - -#include -#include - -#if WARPX_USE_PSATD - -/* \brief Class that updates the field in spectral space and stores the coefficients - * of the corresponding update equation, according to the comoving spectral scheme. - */ -class ComovingPsatdAlgorithm : public SpectralBaseAlgorithm -{ - public: - - /** - * \brief Class constructor - */ - ComovingPsatdAlgorithm (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::Vector& v_comoving, - const amrex::Real dt, - const bool update_with_rho); - - /** - * \brief Override the update equations in Fourier space - */ - virtual void pushSpectralFields (SpectralFieldData& f) const override final; - - /* \brief Initialize the coefficients needed in the update equations - */ - void InitializeSpectralCoefficients (const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, - const amrex::Real dt); - - /** - * \brief Virtual function for current correction in Fourier space. - * 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,3>& current, - const std::unique_ptr& rho) override final; - - /** - * \brief Virtual function for Vay current deposition in Fourier space. - * 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,3>& current) override final; - - private: - - // Real and complex spectral coefficients - SpectralRealCoefficients C_coef, S_ck_coef; - SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef; - - SpectralFieldIndex m_spectral_index; - - // k vectors - KVectorComponent kx_vec; -#if defined(WARPX_DIM_3D) - KVectorComponent ky_vec; -#endif - KVectorComponent kz_vec; - - // Additional member variables - amrex::Vector m_v_comoving; - amrex::Real m_dt; -}; - -#endif // WARPX_USE_PSATD -#endif // WARPX_COMOVING_PSATD_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp deleted file mode 100644 index 0e2953624..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp +++ /dev/null @@ -1,526 +0,0 @@ -#include "ComovingPsatdAlgorithm.H" - -#include "Utils/WarpXConst.H" -#include "Utils/WarpX_Complex.H" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#if WARPX_USE_PSATD - -using namespace amrex; - -ComovingPsatdAlgorithm::ComovingPsatdAlgorithm (const SpectralKSpace& spectral_kspace, - const 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::Vector& v_comoving, - const amrex::Real dt, - const bool update_with_rho) - // Members initialization - : SpectralBaseAlgorithm(spectral_kspace, dm, spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards), - m_spectral_index(spectral_index), - // Initialize the infinite-order k vectors (the argument n_order = -1 selects - // the infinite order option, the argument nodal = false is then irrelevant) - kx_vec(spectral_kspace.getModifiedKComponent(dm, 0, -1, false)), -#if defined(WARPX_DIM_3D) - ky_vec(spectral_kspace.getModifiedKComponent(dm, 1, -1, false)), - kz_vec(spectral_kspace.getModifiedKComponent(dm, 2, -1, false)), -#else - kz_vec(spectral_kspace.getModifiedKComponent(dm, 1, -1, false)), -#endif - m_v_comoving(v_comoving), - m_dt(dt) -{ - amrex::ignore_unused(update_with_rho); - - const BoxArray& ba = spectral_kspace.spectralspace_ba; - - // Allocate arrays of real spectral coefficients - C_coef = SpectralRealCoefficients(ba, dm, 1, 0); - S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0); - - // Allocate arrays of complex spectral coefficients - X1_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - X2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - X3_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - X4_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - Theta2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - - // Initialize real and complex spectral coefficients - InitializeSpectralCoefficients(spectral_kspace, dm, dt); -} - -void -ComovingPsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const -{ - 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 fields = f.fields[mfi].array(); - - // Extract arrays for the coefficients - amrex::Array4 C_arr = C_coef [mfi].array(); - amrex::Array4 S_ck_arr = S_ck_coef[mfi].array(); - amrex::Array4 X1_arr = X1_coef [mfi].array(); - amrex::Array4 X2_arr = X2_coef [mfi].array(); - amrex::Array4 X3_arr = X3_coef [mfi].array(); - amrex::Array4 X4_arr = X4_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 = fields(i,j,k,Idx.Jx); - const Complex Jy = fields(i,j,k,Idx.Jy); - const Complex Jz = fields(i,j,k,Idx.Jz); - const Complex rho_old = fields(i,j,k,Idx.rho_old); - const Complex rho_new = fields(i,j,k,Idx.rho_new); - - // k vector values - const amrex::Real kx_mod = modified_kx_arr[i]; -#if defined(WARPX_DIM_3D) - const amrex::Real ky_mod = modified_ky_arr[j]; - const amrex::Real kz_mod = modified_kz_arr[k]; -#else - constexpr amrex::Real ky_mod = 0._rt; - const amrex::Real kz_mod = modified_kz_arr[j]; -#endif - - // Physical constant c**2 and imaginary unit - constexpr amrex::Real c2 = PhysConst::c*PhysConst::c; - constexpr Complex I = Complex{0._rt,1._rt}; - - // The definition of these coefficients is explained in more detail - // in the function InitializeSpectralCoefficients below - const amrex::Real C = C_arr(i,j,k); - const amrex::Real S_ck = S_ck_arr(i,j,k); - const Complex X1 = X1_arr(i,j,k); - const Complex X2 = X2_arr(i,j,k); - const Complex X3 = X3_arr(i,j,k); - const Complex X4 = X4_arr(i,j,k); - - // Update E - fields(i,j,k,Idx.Ex) = C*Ex_old + S_ck*c2*I*(ky_mod*Bz_old - kz_mod*By_old) - + X4*Jx - I*(X2*rho_new - X3*rho_old)*kx_mod; - - fields(i,j,k,Idx.Ey) = C*Ey_old + S_ck*c2*I*(kz_mod*Bx_old - kx_mod*Bz_old) - + X4*Jy - I*(X2*rho_new - X3*rho_old)*ky_mod; - - fields(i,j,k,Idx.Ez) = C*Ez_old + S_ck*c2*I*(kx_mod*By_old - ky_mod*Bx_old) - + X4*Jz - I*(X2*rho_new - X3*rho_old)*kz_mod; - - // Update B - fields(i,j,k,Idx.Bx) = C*Bx_old - S_ck*I*(ky_mod*Ez_old - kz_mod*Ey_old) - + X1*I*(ky_mod*Jz - kz_mod*Jy); - - fields(i,j,k,Idx.By) = C*By_old - S_ck*I*(kz_mod*Ex_old - kx_mod*Ez_old) - + X1*I*(kz_mod*Jx - kx_mod*Jz); - - fields(i,j,k,Idx.Bz) = C*Bz_old - S_ck*I*(kx_mod*Ey_old - ky_mod*Ex_old) - + X1*I*(kx_mod*Jy - ky_mod*Jx); - }); - } -} - -void ComovingPsatdAlgorithm::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_mod = modified_kx_vec[mfi].dataPtr(); - const amrex::Real* kx = kx_vec[mfi].dataPtr(); -#if defined(WARPX_DIM_3D) - const amrex::Real* ky_mod = modified_ky_vec[mfi].dataPtr(); - const amrex::Real* ky = ky_vec[mfi].dataPtr(); -#endif - const amrex::Real* kz_mod = modified_kz_vec[mfi].dataPtr(); - const amrex::Real* kz = kz_vec[mfi].dataPtr(); - - // Extract arrays for the coefficients - amrex::Array4 C = C_coef [mfi].array(); - amrex::Array4 S_ck = S_ck_coef [mfi].array(); - amrex::Array4 X1 = X1_coef [mfi].array(); - amrex::Array4 X2 = X2_coef [mfi].array(); - amrex::Array4 X3 = X3_coef [mfi].array(); - amrex::Array4 X4 = X4_coef [mfi].array(); - amrex::Array4 T2 = Theta2_coef[mfi].array(); - - // Store comoving velocity - const amrex::Real vx = m_v_comoving[0]; -#if defined(WARPX_DIM_3D) - const amrex::Real vy = m_v_comoving[1]; -#endif - const amrex::Real vz = m_v_comoving[2]; - - // Loop over indices within one box - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - // Calculate norm of finite-order k vector - const amrex::Real knorm_mod = std::sqrt( - std::pow(kx_mod[i], 2) + -#if defined(WARPX_DIM_3D) - std::pow(ky_mod[j], 2) + - std::pow(kz_mod[k], 2)); -#else - std::pow(kz_mod[j], 2)); -#endif - // Calculate norm of infinite-order k vector - const amrex::Real knorm = std::sqrt( - std::pow(kx[i], 2) + -#if defined(WARPX_DIM_3D) - std::pow(ky[j], 2) + - std::pow(kz[k], 2)); -#else - std::pow(kz[j], 2)); -#endif - // Physical constants c, c**2, and epsilon_0, and imaginary unit - constexpr amrex::Real c = PhysConst::c; - constexpr amrex::Real c2 = c*c; - constexpr amrex::Real ep0 = PhysConst::ep0; - constexpr Complex I = Complex{0._rt, 1._rt}; - - // Auxiliary coefficients used when update_with_rho=false - const amrex::Real dt2 = dt * dt; - - // Calculate dot product of k vector with comoving velocity - const amrex::Real kv = kx[i]*vx + -#if defined(WARPX_DIM_3D) - ky[j]*vy + kz[k]*vz; -#else - kz[j]*vz; -#endif - - if (knorm_mod != 0. && knorm != 0.) { - - // Auxiliary coefficients - const amrex::Real om_mod = c * knorm_mod; - const amrex::Real om2_mod = om_mod * om_mod; - const amrex::Real om = c * knorm; - const amrex::Real om2 = om * om; - const Complex tmp1 = amrex::exp( I * om_mod * dt); - const Complex tmp2 = amrex::exp(- I * om_mod * dt); - const Complex tmp1_sqrt = amrex::exp( I * om_mod * dt * 0.5_rt); - const Complex tmp2_sqrt = amrex::exp(- I * om_mod * dt * 0.5_rt); - - C (i,j,k) = std::cos(om_mod * dt); - S_ck(i,j,k) = std::sin(om_mod * dt) / om_mod; - - const amrex::Real nu = - kv / om; - const Complex theta = amrex::exp( I * nu * om * dt * 0.5_rt); - const Complex theta_star = amrex::exp(- I * nu * om * dt * 0.5_rt); - - T2(i,j,k) = theta * theta; - - if ( (nu != om_mod/om) && (nu != -om_mod/om) && (nu != 0.) ) { - - Complex x1 = om2 / (om2_mod - nu * nu * om2) - * (theta_star - theta * C(i,j,k) + I * nu * om * theta * S_ck(i,j,k)); - - // X1 multiplies i*(k \times J) in the update equation for B - X1(i,j,k) = x1 / (ep0 * om2); - - // X2 multiplies rho_new in the update equation for E - // X3 multiplies rho_old in the update equation for E - X2(i,j,k) = c2 * (x1 * om2_mod - theta * (1._rt - C(i,j,k)) * om2) - / (theta_star - theta) / (ep0 * om2 * om2_mod); - X3(i,j,k) = c2 * (x1 * om2_mod - theta_star * (1._rt - C(i,j,k)) * om2) - / (theta_star - theta) / (ep0 * om2 * om2_mod); - - // X4 multiplies J in the update equation for E - X4(i,j,k) = I * nu * om * X1(i,j,k) - theta * S_ck(i,j,k) / ep0; - } - - // Limits for nu = 0 - if (nu == 0.) { - - // X1 multiplies i*(k \times J) in the update equation for B - X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2_mod); - - // X2 multiplies rho_new in the update equation for E - // X3 multiplies rho_old in the update equation for E - X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2_mod); - X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2_mod); - - // Coefficient multiplying J in update equation for E - X4(i,j,k) = - S_ck(i,j,k) / ep0; - } - - // Limits for nu = omega_mod/omega - if (nu == om_mod/om) { - - // X1 multiplies i*(k \times J) in the update equation for B - X1(i,j,k) = tmp1_sqrt * (1._rt - tmp2 * tmp2 - 2._rt * I * om_mod * dt) / (4._rt * ep0 * om2_mod); - - // X2 multiplies rho_new in the update equation for E - // X3 multiplies rho_old in the update equation for E - X2(i,j,k) = c2 * (- 4._rt + 3._rt * tmp1 + tmp2 - 2._rt * I * om_mod * dt * tmp1) - / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt)); - X3(i,j,k) = c2 * (2._rt - tmp2 - 3._rt * tmp1 + 2._rt * tmp1 * tmp1 - 2._rt * I * om_mod * dt * tmp1) - / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt)); - - // Coefficient multiplying J in update equation for E - X4(i,j,k) = tmp1_sqrt * (I - I * tmp2 * tmp2 - 2._rt * om_mod * dt) / (4._rt * ep0 * om_mod); - } - - // Limits for nu = -omega_mod/omega - if (nu == -om_mod/om) { - - // X1 multiplies i*(k \times J) in the update equation for B - X1(i,j,k) = tmp2_sqrt * (1._rt - tmp1 * tmp1 + 2._rt * I * om_mod * dt) / (4._rt * ep0 * om2_mod); - - // X2 multiplies rho_new in the update equation for E - // X3 multiplies rho_old in the update equation for E - X2(i,j,k) = c2 * (- 3._rt + 4._rt * tmp1 - tmp1 * tmp1 - 2._rt * I * om_mod * dt) - / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt)); - X3(i,j,k) = c2 * (3._rt - 2._rt * tmp2 - 2._rt * tmp1 + tmp1 * tmp1 - 2._rt * I * om_mod * dt) - / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt)); - - // Coefficient multiplying J in update equation for E - X4(i,j,k) = tmp2_sqrt * (- I + I * tmp1 * tmp1 - 2._rt * om_mod * dt) / (4._rt * ep0 * om_mod); - } - } - - // Limits for omega = 0 only - else if (knorm_mod != 0. && knorm == 0.) { - - const amrex::Real om_mod = c * knorm_mod; - const amrex::Real om2_mod = om_mod * om_mod; - - C (i,j,k) = std::cos(om_mod * dt); - S_ck(i,j,k) = std::sin(om_mod * dt) / om_mod; - T2(i,j,k) = 1._rt; - - // X1 multiplies i*(k \times J) in the update equation for B - X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2_mod); - - // X2 multiplies rho_new in the update equation for E - // X3 multiplies rho_old in the update equation for E - X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2_mod); - X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2_mod); - - // Coefficient multiplying J in update equation for E - X4(i,j,k) = - S_ck(i,j,k) / ep0; - - } - - // Limits for omega_mod = 0 only - else if (knorm_mod == 0. && knorm != 0.) { - - const amrex::Real om = c * knorm; - const amrex::Real om2 = om * om; - const amrex::Real nu = - kv / om; - const Complex theta = amrex::exp(I * nu * om * dt * 0.5_rt); - const Complex theta_star = amrex::exp(- I * nu * om * dt * 0.5_rt); - - C(i,j,k) = 1._rt; - S_ck(i,j,k) = dt; - T2(i,j,k) = theta * theta; - - if (nu != 0.) { - // X1 multiplies i*(k \times J) in the update equation for B - X1(i,j,k) = (-theta_star + theta - I * nu * om * dt * theta) - / (ep0 * nu * nu * om2); - - // X2 multiplies rho_new in the update equation for E - // X3 multiplies rho_old in the update equation for E - X2(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om * dt * T2(i,j,k) - + 0.5_rt * nu * nu * om2 * dt * dt * T2(i,j,k)) - / (ep0 * nu * nu * om2 * (T2(i,j,k) - 1._rt)); - X3(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om * dt * T2(i,j,k) - + 0.5_rt * nu * nu * om2 * dt * dt) - / (ep0 * nu * nu * om2 * (T2(i,j,k) - 1._rt)); - - // Coefficient multiplying J in update equation for E - X4(i,j,k) = I * (theta - theta_star) / (ep0 * nu * om); - } - - else { - // X1 multiplies i*(k \times J) in the update equation for B - X1(i,j,k) = dt2 / (2._rt * ep0); - - // X2 multiplies rho_new in the update equation for E - // X3 multiplies rho_old in the update equation for E - X2(i,j,k) = c2 * dt2 / (6._rt * ep0); - X3(i,j,k) = - c2 * dt2 / (3._rt * ep0); - - // Coefficient multiplying J in update equation for E - X4(i,j,k) = -dt / ep0; - } - } - - // Limits for omega_mod = 0 and omega = 0 - else if (knorm_mod == 0. && knorm == 0.) { - - C(i,j,k) = 1._rt; - S_ck(i,j,k) = dt; - T2(i,j,k) = 1._rt; - - // X1 multiplies i*(k \times J) in the update equation for B - X1(i,j,k) = dt2 / (2._rt * ep0); - - // X2 multiplies rho_new in the update equation for E - // X3 multiplies rho_old in the update equation for E - X2(i,j,k) = c2 * dt2 / (6._rt * ep0); - X3(i,j,k) = - c2 * dt2 / (3._rt * ep0); - - // Coefficient multiplying J in update equation for E - X4(i,j,k) = -dt / ep0; - } - }); - } -} - -void -ComovingPsatdAlgorithm::CurrentCorrection (const int lev, - SpectralFieldData& field_data, - std::array,3>& current, - const std::unique_ptr& rho) -{ - // Profiling - BL_PROFILE("ComovingAlgorithm::CurrentCorrection"); - - const SpectralFieldIndex& Idx = m_spectral_index; - - // Forward Fourier transform of J and rho - field_data.ForwardTransform(lev, *current[0], Idx.Jx, 0); - field_data.ForwardTransform(lev, *current[1], Idx.Jy, 0); - field_data.ForwardTransform(lev, *current[2], Idx.Jz, 0); - field_data.ForwardTransform(lev, *rho, Idx.rho_old, 0); - field_data.ForwardTransform(lev, *rho, Idx.rho_new, 1); - - const amrex::IntVect& fill_guards = m_fill_guards; - - // Loop over boxes - for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){ - - const amrex::Box& bx = field_data.fields[mfi].box(); - - // Extract arrays for the fields to be updated - amrex::Array4 fields = field_data.fields[mfi].array(); - - // Extract pointers for the k vectors - const amrex::Real* const modified_kx_arr = modified_kx_vec[mfi].dataPtr(); - const amrex::Real* const kx_arr = kx_vec[mfi].dataPtr(); -#if defined(WARPX_DIM_3D) - const amrex::Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr(); - const amrex::Real* const ky_arr = ky_vec[mfi].dataPtr(); -#endif - const amrex::Real* const modified_kz_arr = modified_kz_vec[mfi].dataPtr(); - const amrex::Real* const kz_arr = kz_vec[mfi].dataPtr(); - - // Local copy of member variables before GPU loop - const amrex::Real dt = m_dt; - - // Comoving velocity - const amrex::Real vx = m_v_comoving[0]; - const amrex::Real vy = m_v_comoving[1]; - const amrex::Real vz = m_v_comoving[2]; - - // Loop over indices within one box - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - // Shortcuts for the values of J and rho - const Complex Jx = fields(i,j,k,Idx.Jx); - const Complex Jy = fields(i,j,k,Idx.Jy); - const Complex Jz = fields(i,j,k,Idx.Jz); - const Complex rho_old = fields(i,j,k,Idx.rho_old); - const Complex rho_new = fields(i,j,k,Idx.rho_new); - - // k vector values, and coefficients - const amrex::Real kx_mod = modified_kx_arr[i]; - const amrex::Real kx = kx_arr[i]; -#if defined(WARPX_DIM_3D) - const amrex::Real ky_mod = modified_ky_arr[j]; - const amrex::Real kz_mod = modified_kz_arr[k]; - const amrex::Real ky = ky_arr[j]; - const amrex::Real kz = kz_arr[k]; -#else - constexpr amrex::Real ky_mod = 0._rt; - const amrex::Real kz_mod = modified_kz_arr[j]; - constexpr amrex::Real ky = 0._rt; - const amrex::Real kz = kz_arr[j]; -#endif - constexpr Complex I = Complex{0._rt,1._rt}; - - const amrex::Real knorm_mod = std::sqrt(kx_mod * kx_mod + ky_mod * ky_mod + kz_mod * kz_mod); - - // Correct J - if (knorm_mod != 0._rt) - { - const Complex kmod_dot_J = kx_mod * Jx + ky_mod * Jy + kz_mod * Jz; - const amrex::Real k_dot_v = kx * vx + ky * vy + kz * vz; - - if ( k_dot_v != 0._rt ) { - - const Complex theta = amrex::exp(- I * k_dot_v * dt * 0.5_rt); - const Complex den = 1._rt - theta * theta; - - fields(i,j,k,Idx.Jx) = Jx - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kx_mod / (knorm_mod * knorm_mod); - fields(i,j,k,Idx.Jy) = Jy - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * ky_mod / (knorm_mod * knorm_mod); - fields(i,j,k,Idx.Jz) = Jz - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kz_mod / (knorm_mod * knorm_mod); - - } else { - - fields(i,j,k,Idx.Jx) = Jx - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kx_mod / (knorm_mod * knorm_mod); - fields(i,j,k,Idx.Jy) = Jy - (kmod_dot_J - I * (rho_new - rho_old) / dt) * ky_mod / (knorm_mod * knorm_mod); - fields(i,j,k,Idx.Jz) = Jz - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kz_mod / (knorm_mod * knorm_mod); - } - } - }); - } - - // Backward Fourier transform of J - field_data.BackwardTransform(lev, *current[0], Idx.Jx, 0, fill_guards); - field_data.BackwardTransform(lev, *current[1], Idx.Jy, 0, fill_guards); - field_data.BackwardTransform(lev, *current[2], Idx.Jz, 0, fill_guards); -} - -void -ComovingPsatdAlgorithm::VayDeposition (const int /*lev*/, - SpectralFieldData& /*field_data*/, - std::array,3>& /*current*/) -{ - amrex::Abort("Vay deposition not implemented for comoving PSATD"); -} - -#endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H deleted file mode 100644 index 10ad6dbc5..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H +++ /dev/null @@ -1,78 +0,0 @@ -/* Copyright 2019 David Grote - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#ifndef WARPX_GALILEAN_PSATD_ALGORITHM_RZ_H_ -#define WARPX_GALILEAN_PSATD_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 GalileanPsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ -{ - - public: - GalileanPsatdAlgorithmRZ (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, - const amrex::Vector& v_galilean, - amrex::Real const dt_step, - bool const update_with_rho); - // 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 - * This function overrides the virtual function \c CurrentCorrection in the - * base class \c SpectralBaseAlgorithmRZ (qualifier \c override) and cannot be - * overridden by further derived classes (qualifier \c final). - * - * \param[in] lev mesh-refinement level - * \param[in,out] field_data all fields in Fourier space - * \param[in,out] current two-dimensional array of unique pointers to MultiFab - * storing the three components of the current density - * \param[in] rho unique pointer to MultiFab storing the charge density - */ - virtual void CurrentCorrection ( const int lev, - SpectralFieldDataRZ& field_data, - std::array,3>& current, - const std::unique_ptr& rho ) override final; - - /** - * \brief Virtual function for Vay current deposition in Fourier space - * This function overrides the virtual function \c VayDeposition in the - * base class \c SpectralBaseAlgorithmRZ and cannot be overridden by further - * derived classes. - * - * \param[in] lev 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, SpectralFieldDataRZ& field_data, - std::array,3>& current) override final; - - private: - - SpectralFieldIndex m_spectral_index; - - bool coefficients_initialized; - // Note that dt and v_galilean are saved to use in InitializeSpectralCoefficients - amrex::Real const m_dt; - amrex::Vector m_v_galilean; - bool m_update_with_rho; - - SpectralRealCoefficients C_coef, S_ck_coef; - SpectralComplexCoefficients Theta2_coef, T_rho_coef, X1_coef, X2_coef, X3_coef, X4_coef; - -}; - -#endif // WARPX_GALILEAN_PSATD_ALGORITHM_RZ_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.cpp deleted file mode 100644 index 4fa68b681..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.cpp +++ /dev/null @@ -1,385 +0,0 @@ -/* Copyright 2019-2020 David Grote - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#include "GalileanPsatdAlgorithmRZ.H" -#include "Utils/WarpXConst.H" -#include "Utils/WarpXProfilerWrapper.H" -#include "WarpX.H" - -#include - -using namespace amrex::literals; - - -/* \brief Initialize coefficients for the update equation */ -GalileanPsatdAlgorithmRZ::GalileanPsatdAlgorithmRZ (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, - const amrex::Vector& v_galilean, - amrex::Real const dt, - bool const update_with_rho) - // Initialize members of base class - : SpectralBaseAlgorithmRZ(spectral_kspace, dm, spectral_index, norder_z, nodal), - m_spectral_index(spectral_index), - m_dt(dt), - m_v_galilean(v_galilean), - m_update_with_rho(update_with_rho) -{ - - // 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); - X1_coef = SpectralComplexCoefficients(ba, dm, n_rz_azimuthal_modes, 0); - X2_coef = SpectralComplexCoefficients(ba, dm, n_rz_azimuthal_modes, 0); - X3_coef = SpectralComplexCoefficients(ba, dm, n_rz_azimuthal_modes, 0); - X4_coef = SpectralComplexCoefficients(ba, dm, n_rz_azimuthal_modes, 0); - Theta2_coef = SpectralComplexCoefficients(ba, dm, n_rz_azimuthal_modes, 0); - T_rho_coef = SpectralComplexCoefficients(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 - * The algorithm is described in https://doi.org/10.1103/PhysRevE.94.053305 - * */ -void -GalileanPsatdAlgorithmRZ::pushSpectralFields (SpectralFieldDataRZ & f) -{ - - bool const update_with_rho = m_update_with_rho; - - 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 const& fields = f.fields[mfi].array(); - // Extract arrays for the coefficients - amrex::Array4 const& C_arr = C_coef[mfi].array(); - amrex::Array4 const& S_ck_arr = S_ck_coef[mfi].array(); - amrex::Array4 const& X1_arr = X1_coef[mfi].array(); - amrex::Array4 const& X2_arr = X2_coef[mfi].array(); - amrex::Array4 const& X3_arr = X3_coef[mfi].array(); - amrex::Array4 const& X4_arr = X4_coef[mfi].array(); - amrex::Array4 const& Theta2_arr = Theta2_coef[mfi].array(); - amrex::Array4 const& T_rho_arr = T_rho_coef[mfi].array(); - - // Extract pointers for the k vectors - auto const & kr_modes = f.getKrArray(mfi); - amrex::Real const* kr_arr = kr_modes.dataPtr(); - amrex::Real const* modified_kz_arr = modified_kz_vec[mfi].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 - auto const Ep_m = Idx.Ex + Idx.n_fields*mode; - auto const Em_m = Idx.Ey + Idx.n_fields*mode; - auto const Ez_m = Idx.Ez + Idx.n_fields*mode; - auto const Bp_m = Idx.Bx + Idx.n_fields*mode; - auto const Bm_m = Idx.By + Idx.n_fields*mode; - auto const Bz_m = Idx.Bz + Idx.n_fields*mode; - auto const Jp_m = Idx.Jx + Idx.n_fields*mode; - auto const Jm_m = Idx.Jy + Idx.n_fields*mode; - auto const Jz_m = Idx.Jz + Idx.n_fields*mode; - auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode; - auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode; - - // Record old values of the fields to be updated - Complex const Ep_old = fields(i,j,k,Ep_m); - Complex const Em_old = fields(i,j,k,Em_m); - Complex const Ez_old = fields(i,j,k,Ez_m); - Complex const Bp_old = fields(i,j,k,Bp_m); - Complex const Bm_old = fields(i,j,k,Bm_m); - Complex const Bz_old = fields(i,j,k,Bz_m); - // Shortcut for the values of J and rho - Complex const Jp = fields(i,j,k,Jp_m); - Complex const Jm = fields(i,j,k,Jm_m); - Complex const Jz = fields(i,j,k,Jz_m); - Complex const rho_old = fields(i,j,k,rho_old_m); - Complex const rho_new = fields(i,j,k,rho_new_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]; - amrex::Real const kz = modified_kz_arr[j]; - - 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); - Complex const X1 = X1_arr(i,j,k,mode); - Complex const X2 = X2_arr(i,j,k,mode); - Complex const X3 = X3_arr(i,j,k,mode); - Complex const X4 = X4_arr(i,j,k,mode); - Complex const T2 = Theta2_arr(i,j,k,mode); - Complex const T_rho = T_rho_arr(i,j,k,mode); - - Complex rho_diff; - if (update_with_rho) { - rho_diff = X2*rho_new - T2*X3*rho_old; - } else { - Complex const divE = kr*(Ep_old - Em_old) + I*kz*Ez_old; - Complex const divJ = kr*(Jp - Jm) + I*kz*Jz; - - auto const myeps0 = PhysConst::ep0; // temporary for NVCC - rho_diff = T2*(X2 - X3)*myeps0*divE + T_rho*X2*divJ; - } - - // Update E (see WarpX online documentation: theory section) - fields(i,j,k,Ep_m) = T2*C*Ep_old - + T2*S_ck*(-c2*I*kr/2._rt*Bz_old + c2*kz*Bp_old) - + X4*Jp + 0.5_rt*kr*rho_diff; - fields(i,j,k,Em_m) = T2*C*Em_old - + T2*S_ck*(-c2*I*kr/2._rt*Bz_old - c2*kz*Bm_old) - + X4*Jm - 0.5_rt*kr*rho_diff; - fields(i,j,k,Ez_m) = T2*C*Ez_old - + T2*S_ck*(c2*I*kr*Bp_old + c2*I*kr*Bm_old) - + X4*Jz - I*kz*rho_diff; - // Update B (see WarpX online documentation: theory section) - // Note: here X1 is T2*x1/(ep0*c*c*k_norm*k_norm), where - // x1 has the same definition as in the original paper - fields(i,j,k,Bp_m) = T2*C*Bp_old - - T2*S_ck*(-I*kr/2._rt*Ez_old + kz*Ep_old) - + X1*(-I*kr/2._rt*Jz + kz*Jp); - fields(i,j,k,Bm_m) = T2*C*Bm_old - - T2*S_ck*(-I*kr/2._rt*Ez_old - kz*Em_old) - + X1*(-I*kr/2._rt*Jz - kz*Jm); - fields(i,j,k,Bz_m) = T2*C*Bz_old - - T2*S_ck*I*(kr*Ep_old + kr*Em_old) - + X1*I*(kr*Jp + kr*Jm); - }); - } -} - -void GalileanPsatdAlgorithmRZ::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 const& C = C_coef[mfi].array(); - amrex::Array4 const& S_ck = S_ck_coef[mfi].array(); - amrex::Array4 const& X1 = X1_coef[mfi].array(); - amrex::Array4 const& X2 = X2_coef[mfi].array(); - amrex::Array4 const& X3 = X3_coef[mfi].array(); - amrex::Array4 const& X4 = X4_coef[mfi].array(); - amrex::Array4 const& Theta2 = Theta2_coef[mfi].array(); - amrex::Array4 const& T_rho = T_rho_coef[mfi].array(); - - // Extract real (for portability on GPU) - amrex::Real vz = m_v_galilean[2]; - - auto 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 - { - constexpr amrex::Real c = PhysConst::c; - constexpr amrex::Real ep0 = PhysConst::ep0; - Complex const I = Complex{0._rt,1._rt}; - - // 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 - 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); - - // Calculate dot product with galilean velocity - amrex::Real const kv = kz*vz; - - amrex::Real const nu = kv/(k_norm*c); - Complex const theta = amrex::exp( 0.5_rt*I*kv*dt ); - Complex const theta_star = amrex::exp( -0.5_rt*I*kv*dt ); - Complex const e_theta = amrex::exp( I*c*k_norm*dt ); - - Theta2(i,j,k,mode) = theta*theta; - - if (kz == 0._rt) { - T_rho(i,j,k,mode) = -dt; - } else { - T_rho(i,j,k,mode) = (1._rt - theta*theta)/(I*kz*vz); - } - - if ( (nu != 1._rt) && (nu != 0._rt) ) { - - // Note: the coefficients X1, X2, X do not correspond - // exactly to the original Galilean paper, but the - // update equation have been modified accordingly so that - // the expressions below (with the update equations) - // are mathematically equivalent to those of the paper. - Complex x1 = 1._rt/(1._rt-nu*nu) * - (theta_star - C(i,j,k,mode)*theta + I*kv*S_ck(i,j,k,mode)*theta); - // x1, above, is identical to the original paper - X1(i,j,k,mode) = theta*x1/(ep0*c*c*k_norm*k_norm); - // The difference betwen X2 and X3 below, and those - // from the original paper is the factor ep0*k_norm*k_norm - X2(i,j,k,mode) = (x1 - theta*(1._rt - C(i,j,k,mode))) - /(theta_star-theta)/(ep0*k_norm*k_norm); - X3(i,j,k,mode) = (x1 - theta_star*(1._rt - C(i,j,k,mode))) - /(theta_star-theta)/(ep0*k_norm*k_norm); - X4(i,j,k,mode) = I*kv*X1(i,j,k,mode) - theta*theta*S_ck(i,j,k,mode)/ep0; - - } else if (nu == 0._rt) { - - X1(i,j,k,mode) = (1._rt - C(i,j,k,mode))/(ep0 * c*c * k_norm*k_norm); - X2(i,j,k,mode) = (1._rt - S_ck(i,j,k,mode)/dt)/(ep0 * k_norm*k_norm); - X3(i,j,k,mode) = (C(i,j,k,mode) - S_ck(i,j,k,mode)/dt)/(ep0 * k_norm*k_norm); - X4(i,j,k,mode) = -S_ck(i,j,k,mode)/ep0; - - } else if ( nu == 1._rt) { - X1(i,j,k,mode) = (1._rt - e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*c*c*ep0*k_norm*k_norm); - X2(i,j,k,mode) = (3._rt - 4._rt*e_theta + e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*ep0*k_norm*k_norm*(1._rt - e_theta)); - X3(i,j,k,mode) = (3._rt - 2._rt/e_theta - 2._rt*e_theta + e_theta*e_theta - 2._rt*I*c*k_norm*dt) / (4._rt*ep0*(e_theta - 1._rt)*k_norm*k_norm); - X4(i,j,k,mode) = I*(-1._rt + e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*ep0*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; - X1(i,j,k,mode) = 0.5_rt * dt*dt / ep0; - X2(i,j,k,mode) = c*c * dt*dt / (6._rt*ep0); - X3(i,j,k,mode) = - c*c * dt*dt / (3._rt*ep0); - X4(i,j,k,mode) = -dt/ep0; - Theta2(i,j,k,mode) = 1._rt; - } - }); - } -} - -void -GalileanPsatdAlgorithmRZ::CurrentCorrection (const int lev, - SpectralFieldDataRZ& field_data, - std::array,3>& current, - const std::unique_ptr& rho ) -{ - // Profiling - WARPX_PROFILE( "GalileanPsatdAlgorithmRZ::CurrentCorrection" ); - - const SpectralFieldIndex& Idx = m_spectral_index; - - // Forward Fourier transform of J and rho - field_data.ForwardTransform( lev, *current[0], Idx.Jx, - *current[1], Idx.Jy); - field_data.ForwardTransform( lev, *current[2], Idx.Jz, 0); - field_data.ForwardTransform( lev, *rho, Idx.rho_old, 0 ); - field_data.ForwardTransform( lev, *rho, Idx.rho_new, 1 ); - - // Loop over boxes - for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){ - - amrex::Box const & bx = field_data.fields[mfi].box(); - - // Extract arrays for the fields to be updated - amrex::Array4 fields = field_data.fields[mfi].array(); - - // Extract pointers for the k vectors - auto const & kr_modes = field_data.getKrArray(mfi); - amrex::Real const* kr_arr = kr_modes.dataPtr(); - amrex::Real const* modified_kz_arr = modified_kz_vec[mfi].dataPtr(); - int const nr = bx.length(0); - - // Local copy of member variables before GPU loop - amrex::Real vz = m_v_galilean[2]; - amrex::Real const dt = m_dt; - - // Loop over indices within one box - int const modes = field_data.n_rz_azimuthal_modes; - 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 - auto const Jp_m = Idx.Jx + Idx.n_fields*mode; - auto const Jm_m = Idx.Jy + Idx.n_fields*mode; - auto const Jz_m = Idx.Jz + Idx.n_fields*mode; - auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode; - auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode; - - // Shortcuts for the values of J and rho - Complex const Jp = fields(i,j,k,Jp_m); - Complex const Jm = fields(i,j,k,Jm_m); - Complex const Jz = fields(i,j,k,Jz_m); - Complex const rho_old = fields(i,j,k,rho_old_m); - Complex const rho_new = fields(i,j,k,rho_new_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]; - amrex::Real const kz = modified_kz_arr[j]; - amrex::Real const k_norm2 = kr*kr + kz*kz; - - constexpr Complex I = Complex{0._rt,1._rt}; - - // Correct J - if ( k_norm2 != 0._rt ) - { - Complex const theta2 = amrex::exp(I*kz*vz*dt); - Complex const inv_1_T2 = 1._rt/(kz*vz == 0._rt ? 1._rt : 1._rt - theta2); - Complex const j_corr_coef = (kz == 0._rt ? 1._rt/dt : -I*kz*vz*inv_1_T2); - Complex const F = - (j_corr_coef*(rho_new - rho_old*theta2) + I*kz*Jz + kr*(Jp - Jm))/k_norm2; - - fields(i,j,k,Jp_m) += +0.5_rt*kr*F; - fields(i,j,k,Jm_m) += -0.5_rt*kr*F; - fields(i,j,k,Jz_m) += -I*kz*F; - } - }); - } - - // Backward Fourier transform of J - field_data.BackwardTransform( lev, - *current[0], Idx.Jx, - *current[1], Idx.Jy); - field_data.BackwardTransform( lev, *current[2], Idx.Jz, 0 ); - -} - -void -GalileanPsatdAlgorithmRZ::VayDeposition (const int /*lev*/, - SpectralFieldDataRZ& /*field_data*/, - std::array,3>& /*current*/) -{ - amrex::Abort("Vay deposition not implemented in RZ geometry"); -} diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package index e4c32ea2b..7c2bf4281 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -1,14 +1,14 @@ CEXE_sources += SpectralBaseAlgorithm.cpp CEXE_sources += PsatdAlgorithm.cpp -CEXE_sources += PsatdAlgorithmMultiJ.cpp -CEXE_sources += PMLPsatdAlgorithm.cpp -CEXE_sources += ComovingPsatdAlgorithm.cpp +CEXE_sources += PsatdAlgorithmJLinearInTime.cpp +CEXE_sources += PsatdAlgorithmPml.cpp +CEXE_sources += PsatdAlgorithmComoving.cpp ifeq ($(USE_RZ),TRUE) CEXE_sources += SpectralBaseAlgorithmRZ.cpp CEXE_sources += PsatdAlgorithmRZ.cpp - CEXE_sources += GalileanPsatdAlgorithmRZ.cpp - CEXE_sources += PMLPsatdAlgorithmRZ.cpp + CEXE_sources += PsatdAlgorithmGalileanRZ.cpp + CEXE_sources += PsatdAlgorithmPmlRZ.cpp endif VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H deleted file mode 100644 index d6d7e074e..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H +++ /dev/null @@ -1,91 +0,0 @@ -/* Copyright 2019 Axel Huebl, Remi Lehe - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#ifndef WARPX_PML_PSATD_ALGORITHM_H_ -#define WARPX_PML_PSATD_ALGORITHM_H_ - -#include "SpectralBaseAlgorithm.H" - -#include "FieldSolver/SpectralSolver/SpectralFieldData_fwd.H" -#include "FieldSolver/SpectralSolver/SpectralKSpace_fwd.H" - -#include - -#include - -#include -#include - -#if WARPX_USE_PSATD - -/* \brief Class that updates the field in spectral space - * and stores the coefficients of the corresponding update equation. - */ -class PMLPsatdAlgorithm : public SpectralBaseAlgorithm -{ - public: - PMLPsatdAlgorithm(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 dive_cleaning, - const bool divb_cleaning); - - void InitializeSpectralCoefficients( - const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, - const amrex::Real dt); - - // Redefine functions from base class - virtual void pushSpectralFields(SpectralFieldData& f) const override final; - - /** - * \brief Virtual function for current correction in Fourier space - * ( Vay et al, 2013). - * 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,3>& current, - const std::unique_ptr& rho) override final; - - /** - * \brief Virtual function for Vay current deposition in Fourier space - * ( Vay et al, 2013). - * 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,3>& current) override final; - - private: - SpectralFieldIndex m_spectral_index; - 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 -#endif // WARPX_PML_PSATD_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp deleted file mode 100644 index bfe02a238..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp +++ /dev/null @@ -1,421 +0,0 @@ -/* Copyright 2019 Remi Lehe - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#include "PMLPsatdAlgorithm.H" - -#include "FieldSolver/SpectralSolver/SpectralFieldData.H" -#include "FieldSolver/SpectralSolver/SpectralKSpace.H" -#include "Utils/WarpXConst.H" -#include "Utils/WarpX_Complex.H" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#if WARPX_USE_PSATD - -using namespace amrex; - -/* \brief Initialize coefficients for the update equation */ -PMLPsatdAlgorithm::PMLPsatdAlgorithm(const SpectralKSpace& spectral_kspace, - const 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 Real dt, - const bool dive_cleaning, const bool divb_cleaning) - // Initialize members of base class - : SpectralBaseAlgorithm(spectral_kspace, dm, spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards), - m_spectral_index(spectral_index), - m_dt(dt), - m_dive_cleaning(dive_cleaning), - m_divb_cleaning(divb_cleaning) -{ - const BoxArray& ba = spectral_kspace.spectralspace_ba; - - // Allocate the arrays of coefficients - C_coef = SpectralRealCoefficients(ba, dm, 1, 0); - S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0); - inv_k2_coef = SpectralRealCoefficients(ba, dm, 1, 0); - - InitializeSpectralCoefficients(spectral_kspace, dm, dt); -} - -/* Advance the E and B field in spectral space (stored in `f`) - * over one time step */ -void -PMLPsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const { - - const bool dive_cleaning = m_dive_cleaning; - const bool divb_cleaning = m_divb_cleaning; - - const SpectralFieldIndex& Idx = m_spectral_index; - - // Loop over boxes - for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){ - - const Box& bx = f.fields[mfi].box(); - - // Extract arrays for the fields to be updated - Array4 fields = f.fields[mfi].array(); - - // Extract arrays for the coefficients - Array4 C_arr = C_coef[mfi].array(); - Array4 S_ck_arr = S_ck_coef[mfi].array(); - Array4 inv_k2_arr = inv_k2_coef[mfi].array(); - - // Extract pointers for the k vectors - const Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr(); -#if defined(WARPX_DIM_3D) - const Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr(); -#endif - const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr(); - - const amrex::Real dt = m_dt; - - // Loop over indices within one box - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - // Record old values of the fields to be updated - const Complex Exy = fields(i,j,k,Idx.Exy); - const Complex Exz = fields(i,j,k,Idx.Exz); - const Complex Eyx = fields(i,j,k,Idx.Eyx); - const Complex Eyz = fields(i,j,k,Idx.Eyz); - const Complex Ezx = fields(i,j,k,Idx.Ezx); - const Complex Ezy = fields(i,j,k,Idx.Ezy); - - const Complex Bxy = fields(i,j,k,Idx.Bxy); - const Complex Bxz = fields(i,j,k,Idx.Bxz); - const Complex Byx = fields(i,j,k,Idx.Byx); - const Complex Byz = fields(i,j,k,Idx.Byz); - 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 defined(WARPX_DIM_3D) - const Real ky = modified_ky_arr[j]; - const Real kz = modified_kz_arr[k]; -#else - constexpr Real ky = 0._rt; - const Real kz = modified_kz_arr[j]; -#endif - constexpr Real c2 = PhysConst::c * PhysConst::c; - - const Complex I = Complex{0._rt, 1._rt}; - - const Real kx2 = kx*kx; - const Real ky2 = ky*ky; - const Real kz2 = kz*kz; - const Real k_norm = std::sqrt(kx2 + ky2 + kz2); - - if (k_norm != 0._rt) { - - 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; - } - } - }); - } -} - -void PMLPsatdAlgorithm::InitializeSpectralCoefficients ( - const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, - const amrex::Real dt) -{ - const BoxArray& ba = spectral_kspace.spectralspace_ba; - - // Fill them with the right values: - // Loop over boxes and allocate the corresponding coefficients - // for each box owned by the local MPI proc - for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi) { - - const Box& bx = ba[mfi]; - - // Extract pointers for the k vectors - const Real* modified_kx = modified_kx_vec[mfi].dataPtr(); -#if defined(WARPX_DIM_3D) - const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); -#endif - const Real* modified_kz = modified_kz_vec[mfi].dataPtr(); - - // Extract arrays for the coefficients - Array4 C = C_coef[mfi].array(); - Array4 S_ck = S_ck_coef[mfi].array(); - Array4 inv_k2 = inv_k2_coef[mfi].array(); - - // Loop over indices within one box - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - const Real kx = modified_kx[i]; -#if defined(WARPX_DIM_3D) - const Real ky = modified_ky[j]; - const Real kz = modified_kz[k]; -#else - constexpr Real ky = 0._rt; - const Real kz = modified_kz[j]; -#endif - - // Calculate norm of vector - const Real k_norm = std::sqrt(kx*kx + ky*ky + kz*kz); - const Real k2 = k_norm * k_norm; - - // Calculate coefficients - constexpr Real c = PhysConst::c; - - // Coefficients for k_norm = 0 do not need to be set - if (k_norm != 0._rt) { - C(i,j,k) = std::cos(c * k_norm * dt); - S_ck(i,j,k) = std::sin(c * k_norm * dt) / (c * k_norm); - inv_k2(i,j,k) = 1._rt / k2; - } - }); - } -} - -void -PMLPsatdAlgorithm::CurrentCorrection (const int /*lev*/, - SpectralFieldData& /*field_data*/, - std::array,3>& /*current*/, - const std::unique_ptr& /*rho*/) -{ - amrex::Abort("Current correction not implemented for PML PSATD"); -} - -void -PMLPsatdAlgorithm::VayDeposition (const int /*lev*/, - SpectralFieldData& /*field_data*/, - std::array,3>& /*current*/) -{ - amrex::Abort("Vay deposition not implemented for PML PSATD"); -} - -#endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.H deleted file mode 100644 index e0956a514..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.H +++ /dev/null @@ -1,72 +0,0 @@ -/* 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 - * ( Vay et al, 2013). - * 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,3>& current, - const std::unique_ptr& rho) override final; - - /** - * \brief Virtual function for Vay current deposition in Fourier space - * ( Vay et al, 2013). - * 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,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 deleted file mode 100644 index 892c542c7..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.cpp +++ /dev/null @@ -1,177 +0,0 @@ -/* 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 - -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 const& fields = f.fields[mfi].array(); - // Extract arrays for the coefficients - amrex::Array4 const& C_arr = C_coef[mfi].array(); - amrex::Array4 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 const& C = C_coef[mfi].array(); - amrex::Array4 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,3>& /* current */, - const std::unique_ptr& /* rho */) -{ - amrex::Abort("Current correction not implemented in RZ geometry PML"); -} - -void -PMLPsatdAlgorithmRZ::VayDeposition (const int /* lev */, - SpectralFieldDataRZ& /*field_data*/, - std::array,3>& /*current*/) -{ - amrex::Abort("Vay deposition not implemented in RZ geometry PML"); -} diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.H new file mode 100644 index 000000000..a01719a9a --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.H @@ -0,0 +1,105 @@ +#ifndef WARPX_PSATD_ALGORITHM_COMOVING_H_ +#define WARPX_PSATD_ALGORITHM_COMOVING_H_ + +#include "FieldSolver/SpectralSolver/SpectralFieldData.H" +#include "FieldSolver/SpectralSolver/SpectralKSpace.H" +#include "SpectralBaseAlgorithm.H" + +#include +#include +#include + +#include + +#include +#include + +#if WARPX_USE_PSATD + +/* \brief Class that updates the field in spectral space and stores the coefficients + * of the corresponding update equation, according to the comoving spectral scheme. + */ +class PsatdAlgorithmComoving : public SpectralBaseAlgorithm +{ + public: + + /** + * \brief Class constructor + */ + PsatdAlgorithmComoving (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::Vector& v_comoving, + const amrex::Real dt, + const bool update_with_rho); + + /** + * \brief Override the update equations in Fourier space + */ + virtual void pushSpectralFields (SpectralFieldData& f) const override final; + + /* \brief Initialize the coefficients needed in the update equations + */ + void InitializeSpectralCoefficients (const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const amrex::Real dt); + + /** + * \brief Virtual function for current correction in Fourier space. + * 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,3>& current, + const std::unique_ptr& rho) override final; + + /** + * \brief Virtual function for Vay current deposition in Fourier space. + * 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,3>& current) override final; + + private: + + // Real and complex spectral coefficients + SpectralRealCoefficients C_coef, S_ck_coef; + SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef; + + SpectralFieldIndex m_spectral_index; + + // k vectors + KVectorComponent kx_vec; +#if defined(WARPX_DIM_3D) + KVectorComponent ky_vec; +#endif + KVectorComponent kz_vec; + + // Additional member variables + amrex::Vector m_v_comoving; + amrex::Real m_dt; +}; + +#endif // WARPX_USE_PSATD +#endif // WARPX_PSATD_ALGORITHM_COMOVING_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp new file mode 100644 index 000000000..b14ded092 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp @@ -0,0 +1,526 @@ +#include "PsatdAlgorithmComoving.H" + +#include "Utils/WarpXConst.H" +#include "Utils/WarpX_Complex.H" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#if WARPX_USE_PSATD + +using namespace amrex; + +PsatdAlgorithmComoving::PsatdAlgorithmComoving (const SpectralKSpace& spectral_kspace, + const 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::Vector& v_comoving, + const amrex::Real dt, + const bool update_with_rho) + // Members initialization + : SpectralBaseAlgorithm(spectral_kspace, dm, spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards), + m_spectral_index(spectral_index), + // Initialize the infinite-order k vectors (the argument n_order = -1 selects + // the infinite order option, the argument nodal = false is then irrelevant) + kx_vec(spectral_kspace.getModifiedKComponent(dm, 0, -1, false)), +#if defined(WARPX_DIM_3D) + ky_vec(spectral_kspace.getModifiedKComponent(dm, 1, -1, false)), + kz_vec(spectral_kspace.getModifiedKComponent(dm, 2, -1, false)), +#else + kz_vec(spectral_kspace.getModifiedKComponent(dm, 1, -1, false)), +#endif + m_v_comoving(v_comoving), + m_dt(dt) +{ + amrex::ignore_unused(update_with_rho); + + const BoxArray& ba = spectral_kspace.spectralspace_ba; + + // Allocate arrays of real spectral coefficients + C_coef = SpectralRealCoefficients(ba, dm, 1, 0); + S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0); + + // Allocate arrays of complex spectral coefficients + X1_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + X2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + X3_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + X4_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + Theta2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + + // Initialize real and complex spectral coefficients + InitializeSpectralCoefficients(spectral_kspace, dm, dt); +} + +void +PsatdAlgorithmComoving::pushSpectralFields (SpectralFieldData& f) const +{ + 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 fields = f.fields[mfi].array(); + + // Extract arrays for the coefficients + amrex::Array4 C_arr = C_coef [mfi].array(); + amrex::Array4 S_ck_arr = S_ck_coef[mfi].array(); + amrex::Array4 X1_arr = X1_coef [mfi].array(); + amrex::Array4 X2_arr = X2_coef [mfi].array(); + amrex::Array4 X3_arr = X3_coef [mfi].array(); + amrex::Array4 X4_arr = X4_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 = fields(i,j,k,Idx.Jx); + const Complex Jy = fields(i,j,k,Idx.Jy); + const Complex Jz = fields(i,j,k,Idx.Jz); + const Complex rho_old = fields(i,j,k,Idx.rho_old); + const Complex rho_new = fields(i,j,k,Idx.rho_new); + + // k vector values + const amrex::Real kx_mod = modified_kx_arr[i]; +#if defined(WARPX_DIM_3D) + const amrex::Real ky_mod = modified_ky_arr[j]; + const amrex::Real kz_mod = modified_kz_arr[k]; +#else + constexpr amrex::Real ky_mod = 0._rt; + const amrex::Real kz_mod = modified_kz_arr[j]; +#endif + + // Physical constant c**2 and imaginary unit + constexpr amrex::Real c2 = PhysConst::c*PhysConst::c; + constexpr Complex I = Complex{0._rt,1._rt}; + + // The definition of these coefficients is explained in more detail + // in the function InitializeSpectralCoefficients below + const amrex::Real C = C_arr(i,j,k); + const amrex::Real S_ck = S_ck_arr(i,j,k); + const Complex X1 = X1_arr(i,j,k); + const Complex X2 = X2_arr(i,j,k); + const Complex X3 = X3_arr(i,j,k); + const Complex X4 = X4_arr(i,j,k); + + // Update E + fields(i,j,k,Idx.Ex) = C*Ex_old + S_ck*c2*I*(ky_mod*Bz_old - kz_mod*By_old) + + X4*Jx - I*(X2*rho_new - X3*rho_old)*kx_mod; + + fields(i,j,k,Idx.Ey) = C*Ey_old + S_ck*c2*I*(kz_mod*Bx_old - kx_mod*Bz_old) + + X4*Jy - I*(X2*rho_new - X3*rho_old)*ky_mod; + + fields(i,j,k,Idx.Ez) = C*Ez_old + S_ck*c2*I*(kx_mod*By_old - ky_mod*Bx_old) + + X4*Jz - I*(X2*rho_new - X3*rho_old)*kz_mod; + + // Update B + fields(i,j,k,Idx.Bx) = C*Bx_old - S_ck*I*(ky_mod*Ez_old - kz_mod*Ey_old) + + X1*I*(ky_mod*Jz - kz_mod*Jy); + + fields(i,j,k,Idx.By) = C*By_old - S_ck*I*(kz_mod*Ex_old - kx_mod*Ez_old) + + X1*I*(kz_mod*Jx - kx_mod*Jz); + + fields(i,j,k,Idx.Bz) = C*Bz_old - S_ck*I*(kx_mod*Ey_old - ky_mod*Ex_old) + + X1*I*(kx_mod*Jy - ky_mod*Jx); + }); + } +} + +void PsatdAlgorithmComoving::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_mod = modified_kx_vec[mfi].dataPtr(); + const amrex::Real* kx = kx_vec[mfi].dataPtr(); +#if defined(WARPX_DIM_3D) + const amrex::Real* ky_mod = modified_ky_vec[mfi].dataPtr(); + const amrex::Real* ky = ky_vec[mfi].dataPtr(); +#endif + const amrex::Real* kz_mod = modified_kz_vec[mfi].dataPtr(); + const amrex::Real* kz = kz_vec[mfi].dataPtr(); + + // Extract arrays for the coefficients + amrex::Array4 C = C_coef [mfi].array(); + amrex::Array4 S_ck = S_ck_coef [mfi].array(); + amrex::Array4 X1 = X1_coef [mfi].array(); + amrex::Array4 X2 = X2_coef [mfi].array(); + amrex::Array4 X3 = X3_coef [mfi].array(); + amrex::Array4 X4 = X4_coef [mfi].array(); + amrex::Array4 T2 = Theta2_coef[mfi].array(); + + // Store comoving velocity + const amrex::Real vx = m_v_comoving[0]; +#if defined(WARPX_DIM_3D) + const amrex::Real vy = m_v_comoving[1]; +#endif + const amrex::Real vz = m_v_comoving[2]; + + // Loop over indices within one box + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Calculate norm of finite-order k vector + const amrex::Real knorm_mod = std::sqrt( + std::pow(kx_mod[i], 2) + +#if defined(WARPX_DIM_3D) + std::pow(ky_mod[j], 2) + + std::pow(kz_mod[k], 2)); +#else + std::pow(kz_mod[j], 2)); +#endif + // Calculate norm of infinite-order k vector + const amrex::Real knorm = std::sqrt( + std::pow(kx[i], 2) + +#if defined(WARPX_DIM_3D) + std::pow(ky[j], 2) + + std::pow(kz[k], 2)); +#else + std::pow(kz[j], 2)); +#endif + // Physical constants c, c**2, and epsilon_0, and imaginary unit + constexpr amrex::Real c = PhysConst::c; + constexpr amrex::Real c2 = c*c; + constexpr amrex::Real ep0 = PhysConst::ep0; + constexpr Complex I = Complex{0._rt, 1._rt}; + + // Auxiliary coefficients used when update_with_rho=false + const amrex::Real dt2 = dt * dt; + + // Calculate dot product of k vector with comoving velocity + const amrex::Real kv = kx[i]*vx + +#if defined(WARPX_DIM_3D) + ky[j]*vy + kz[k]*vz; +#else + kz[j]*vz; +#endif + + if (knorm_mod != 0. && knorm != 0.) { + + // Auxiliary coefficients + const amrex::Real om_mod = c * knorm_mod; + const amrex::Real om2_mod = om_mod * om_mod; + const amrex::Real om = c * knorm; + const amrex::Real om2 = om * om; + const Complex tmp1 = amrex::exp( I * om_mod * dt); + const Complex tmp2 = amrex::exp(- I * om_mod * dt); + const Complex tmp1_sqrt = amrex::exp( I * om_mod * dt * 0.5_rt); + const Complex tmp2_sqrt = amrex::exp(- I * om_mod * dt * 0.5_rt); + + C (i,j,k) = std::cos(om_mod * dt); + S_ck(i,j,k) = std::sin(om_mod * dt) / om_mod; + + const amrex::Real nu = - kv / om; + const Complex theta = amrex::exp( I * nu * om * dt * 0.5_rt); + const Complex theta_star = amrex::exp(- I * nu * om * dt * 0.5_rt); + + T2(i,j,k) = theta * theta; + + if ( (nu != om_mod/om) && (nu != -om_mod/om) && (nu != 0.) ) { + + Complex x1 = om2 / (om2_mod - nu * nu * om2) + * (theta_star - theta * C(i,j,k) + I * nu * om * theta * S_ck(i,j,k)); + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = x1 / (ep0 * om2); + + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = c2 * (x1 * om2_mod - theta * (1._rt - C(i,j,k)) * om2) + / (theta_star - theta) / (ep0 * om2 * om2_mod); + X3(i,j,k) = c2 * (x1 * om2_mod - theta_star * (1._rt - C(i,j,k)) * om2) + / (theta_star - theta) / (ep0 * om2 * om2_mod); + + // X4 multiplies J in the update equation for E + X4(i,j,k) = I * nu * om * X1(i,j,k) - theta * S_ck(i,j,k) / ep0; + } + + // Limits for nu = 0 + if (nu == 0.) { + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2_mod); + + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2_mod); + X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2_mod); + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = - S_ck(i,j,k) / ep0; + } + + // Limits for nu = omega_mod/omega + if (nu == om_mod/om) { + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = tmp1_sqrt * (1._rt - tmp2 * tmp2 - 2._rt * I * om_mod * dt) / (4._rt * ep0 * om2_mod); + + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = c2 * (- 4._rt + 3._rt * tmp1 + tmp2 - 2._rt * I * om_mod * dt * tmp1) + / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt)); + X3(i,j,k) = c2 * (2._rt - tmp2 - 3._rt * tmp1 + 2._rt * tmp1 * tmp1 - 2._rt * I * om_mod * dt * tmp1) + / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt)); + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = tmp1_sqrt * (I - I * tmp2 * tmp2 - 2._rt * om_mod * dt) / (4._rt * ep0 * om_mod); + } + + // Limits for nu = -omega_mod/omega + if (nu == -om_mod/om) { + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = tmp2_sqrt * (1._rt - tmp1 * tmp1 + 2._rt * I * om_mod * dt) / (4._rt * ep0 * om2_mod); + + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = c2 * (- 3._rt + 4._rt * tmp1 - tmp1 * tmp1 - 2._rt * I * om_mod * dt) + / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt)); + X3(i,j,k) = c2 * (3._rt - 2._rt * tmp2 - 2._rt * tmp1 + tmp1 * tmp1 - 2._rt * I * om_mod * dt) + / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt)); + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = tmp2_sqrt * (- I + I * tmp1 * tmp1 - 2._rt * om_mod * dt) / (4._rt * ep0 * om_mod); + } + } + + // Limits for omega = 0 only + else if (knorm_mod != 0. && knorm == 0.) { + + const amrex::Real om_mod = c * knorm_mod; + const amrex::Real om2_mod = om_mod * om_mod; + + C (i,j,k) = std::cos(om_mod * dt); + S_ck(i,j,k) = std::sin(om_mod * dt) / om_mod; + T2(i,j,k) = 1._rt; + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2_mod); + + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2_mod); + X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2_mod); + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = - S_ck(i,j,k) / ep0; + + } + + // Limits for omega_mod = 0 only + else if (knorm_mod == 0. && knorm != 0.) { + + const amrex::Real om = c * knorm; + const amrex::Real om2 = om * om; + const amrex::Real nu = - kv / om; + const Complex theta = amrex::exp(I * nu * om * dt * 0.5_rt); + const Complex theta_star = amrex::exp(- I * nu * om * dt * 0.5_rt); + + C(i,j,k) = 1._rt; + S_ck(i,j,k) = dt; + T2(i,j,k) = theta * theta; + + if (nu != 0.) { + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = (-theta_star + theta - I * nu * om * dt * theta) + / (ep0 * nu * nu * om2); + + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om * dt * T2(i,j,k) + + 0.5_rt * nu * nu * om2 * dt * dt * T2(i,j,k)) + / (ep0 * nu * nu * om2 * (T2(i,j,k) - 1._rt)); + X3(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om * dt * T2(i,j,k) + + 0.5_rt * nu * nu * om2 * dt * dt) + / (ep0 * nu * nu * om2 * (T2(i,j,k) - 1._rt)); + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = I * (theta - theta_star) / (ep0 * nu * om); + } + + else { + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = dt2 / (2._rt * ep0); + + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = c2 * dt2 / (6._rt * ep0); + X3(i,j,k) = - c2 * dt2 / (3._rt * ep0); + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = -dt / ep0; + } + } + + // Limits for omega_mod = 0 and omega = 0 + else if (knorm_mod == 0. && knorm == 0.) { + + C(i,j,k) = 1._rt; + S_ck(i,j,k) = dt; + T2(i,j,k) = 1._rt; + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = dt2 / (2._rt * ep0); + + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = c2 * dt2 / (6._rt * ep0); + X3(i,j,k) = - c2 * dt2 / (3._rt * ep0); + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = -dt / ep0; + } + }); + } +} + +void +PsatdAlgorithmComoving::CurrentCorrection (const int lev, + SpectralFieldData& field_data, + std::array,3>& current, + const std::unique_ptr& rho) +{ + // Profiling + BL_PROFILE("PsatdAlgorithmComoving::CurrentCorrection"); + + const SpectralFieldIndex& Idx = m_spectral_index; + + // Forward Fourier transform of J and rho + field_data.ForwardTransform(lev, *current[0], Idx.Jx, 0); + field_data.ForwardTransform(lev, *current[1], Idx.Jy, 0); + field_data.ForwardTransform(lev, *current[2], Idx.Jz, 0); + field_data.ForwardTransform(lev, *rho, Idx.rho_old, 0); + field_data.ForwardTransform(lev, *rho, Idx.rho_new, 1); + + const amrex::IntVect& fill_guards = m_fill_guards; + + // Loop over boxes + for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){ + + const amrex::Box& bx = field_data.fields[mfi].box(); + + // Extract arrays for the fields to be updated + amrex::Array4 fields = field_data.fields[mfi].array(); + + // Extract pointers for the k vectors + const amrex::Real* const modified_kx_arr = modified_kx_vec[mfi].dataPtr(); + const amrex::Real* const kx_arr = kx_vec[mfi].dataPtr(); +#if defined(WARPX_DIM_3D) + const amrex::Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr(); + const amrex::Real* const ky_arr = ky_vec[mfi].dataPtr(); +#endif + const amrex::Real* const modified_kz_arr = modified_kz_vec[mfi].dataPtr(); + const amrex::Real* const kz_arr = kz_vec[mfi].dataPtr(); + + // Local copy of member variables before GPU loop + const amrex::Real dt = m_dt; + + // Comoving velocity + const amrex::Real vx = m_v_comoving[0]; + const amrex::Real vy = m_v_comoving[1]; + const amrex::Real vz = m_v_comoving[2]; + + // Loop over indices within one box + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Shortcuts for the values of J and rho + const Complex Jx = fields(i,j,k,Idx.Jx); + const Complex Jy = fields(i,j,k,Idx.Jy); + const Complex Jz = fields(i,j,k,Idx.Jz); + const Complex rho_old = fields(i,j,k,Idx.rho_old); + const Complex rho_new = fields(i,j,k,Idx.rho_new); + + // k vector values, and coefficients + const amrex::Real kx_mod = modified_kx_arr[i]; + const amrex::Real kx = kx_arr[i]; +#if defined(WARPX_DIM_3D) + const amrex::Real ky_mod = modified_ky_arr[j]; + const amrex::Real kz_mod = modified_kz_arr[k]; + const amrex::Real ky = ky_arr[j]; + const amrex::Real kz = kz_arr[k]; +#else + constexpr amrex::Real ky_mod = 0._rt; + const amrex::Real kz_mod = modified_kz_arr[j]; + constexpr amrex::Real ky = 0._rt; + const amrex::Real kz = kz_arr[j]; +#endif + constexpr Complex I = Complex{0._rt,1._rt}; + + const amrex::Real knorm_mod = std::sqrt(kx_mod * kx_mod + ky_mod * ky_mod + kz_mod * kz_mod); + + // Correct J + if (knorm_mod != 0._rt) + { + const Complex kmod_dot_J = kx_mod * Jx + ky_mod * Jy + kz_mod * Jz; + const amrex::Real k_dot_v = kx * vx + ky * vy + kz * vz; + + if ( k_dot_v != 0._rt ) { + + const Complex theta = amrex::exp(- I * k_dot_v * dt * 0.5_rt); + const Complex den = 1._rt - theta * theta; + + fields(i,j,k,Idx.Jx) = Jx - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kx_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jy) = Jy - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * ky_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jz) = Jz - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kz_mod / (knorm_mod * knorm_mod); + + } else { + + fields(i,j,k,Idx.Jx) = Jx - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kx_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jy) = Jy - (kmod_dot_J - I * (rho_new - rho_old) / dt) * ky_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jz) = Jz - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kz_mod / (knorm_mod * knorm_mod); + } + } + }); + } + + // Backward Fourier transform of J + field_data.BackwardTransform(lev, *current[0], Idx.Jx, 0, fill_guards); + field_data.BackwardTransform(lev, *current[1], Idx.Jy, 0, fill_guards); + field_data.BackwardTransform(lev, *current[2], Idx.Jz, 0, fill_guards); +} + +void +PsatdAlgorithmComoving::VayDeposition (const int /*lev*/, + SpectralFieldData& /*field_data*/, + std::array,3>& /*current*/) +{ + amrex::Abort("Vay deposition not implemented for comoving PSATD"); +} + +#endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.H new file mode 100644 index 000000000..cee84f99d --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.H @@ -0,0 +1,78 @@ +/* Copyright 2019 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PSATD_ALGORITHM_GALILEAN_RZ_H_ +#define WARPX_PSATD_ALGORITHM_GALILEAN_RZ_H_ + +#include "SpectralBaseAlgorithmRZ.H" + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + */ +class PsatdAlgorithmGalileanRZ : public SpectralBaseAlgorithmRZ +{ + + public: + PsatdAlgorithmGalileanRZ (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, + const amrex::Vector& v_galilean, + amrex::Real const dt_step, + bool const update_with_rho); + // 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 + * This function overrides the virtual function \c CurrentCorrection in the + * base class \c SpectralBaseAlgorithmRZ (qualifier \c override) and cannot be + * overridden by further derived classes (qualifier \c final). + * + * \param[in] lev mesh-refinement level + * \param[in,out] field_data all fields in Fourier space + * \param[in,out] current two-dimensional array of unique pointers to MultiFab + * storing the three components of the current density + * \param[in] rho unique pointer to MultiFab storing the charge density + */ + virtual void CurrentCorrection ( const int lev, + SpectralFieldDataRZ& field_data, + std::array,3>& current, + const std::unique_ptr& rho ) override final; + + /** + * \brief Virtual function for Vay current deposition in Fourier space + * This function overrides the virtual function \c VayDeposition in the + * base class \c SpectralBaseAlgorithmRZ and cannot be overridden by further + * derived classes. + * + * \param[in] lev 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, SpectralFieldDataRZ& field_data, + std::array,3>& current) override final; + + private: + + SpectralFieldIndex m_spectral_index; + + bool coefficients_initialized; + // Note that dt and v_galilean are saved to use in InitializeSpectralCoefficients + amrex::Real const m_dt; + amrex::Vector m_v_galilean; + bool m_update_with_rho; + + SpectralRealCoefficients C_coef, S_ck_coef; + SpectralComplexCoefficients Theta2_coef, T_rho_coef, X1_coef, X2_coef, X3_coef, X4_coef; + +}; + +#endif // WARPX_PSATD_ALGORITHM_GALILEAN_RZ_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp new file mode 100644 index 000000000..2faf48c04 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp @@ -0,0 +1,385 @@ +/* Copyright 2019-2020 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "PsatdAlgorithmGalileanRZ.H" +#include "Utils/WarpXConst.H" +#include "Utils/WarpXProfilerWrapper.H" +#include "WarpX.H" + +#include + +using namespace amrex::literals; + + +/* \brief Initialize coefficients for the update equation */ +PsatdAlgorithmGalileanRZ::PsatdAlgorithmGalileanRZ (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, + const amrex::Vector& v_galilean, + amrex::Real const dt, + bool const update_with_rho) + // Initialize members of base class + : SpectralBaseAlgorithmRZ(spectral_kspace, dm, spectral_index, norder_z, nodal), + m_spectral_index(spectral_index), + m_dt(dt), + m_v_galilean(v_galilean), + m_update_with_rho(update_with_rho) +{ + + // 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); + X1_coef = SpectralComplexCoefficients(ba, dm, n_rz_azimuthal_modes, 0); + X2_coef = SpectralComplexCoefficients(ba, dm, n_rz_azimuthal_modes, 0); + X3_coef = SpectralComplexCoefficients(ba, dm, n_rz_azimuthal_modes, 0); + X4_coef = SpectralComplexCoefficients(ba, dm, n_rz_azimuthal_modes, 0); + Theta2_coef = SpectralComplexCoefficients(ba, dm, n_rz_azimuthal_modes, 0); + T_rho_coef = SpectralComplexCoefficients(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 + * The algorithm is described in https://doi.org/10.1103/PhysRevE.94.053305 + * */ +void +PsatdAlgorithmGalileanRZ::pushSpectralFields (SpectralFieldDataRZ & f) +{ + + bool const update_with_rho = m_update_with_rho; + + 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 const& fields = f.fields[mfi].array(); + // Extract arrays for the coefficients + amrex::Array4 const& C_arr = C_coef[mfi].array(); + amrex::Array4 const& S_ck_arr = S_ck_coef[mfi].array(); + amrex::Array4 const& X1_arr = X1_coef[mfi].array(); + amrex::Array4 const& X2_arr = X2_coef[mfi].array(); + amrex::Array4 const& X3_arr = X3_coef[mfi].array(); + amrex::Array4 const& X4_arr = X4_coef[mfi].array(); + amrex::Array4 const& Theta2_arr = Theta2_coef[mfi].array(); + amrex::Array4 const& T_rho_arr = T_rho_coef[mfi].array(); + + // Extract pointers for the k vectors + auto const & kr_modes = f.getKrArray(mfi); + amrex::Real const* kr_arr = kr_modes.dataPtr(); + amrex::Real const* modified_kz_arr = modified_kz_vec[mfi].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 + auto const Ep_m = Idx.Ex + Idx.n_fields*mode; + auto const Em_m = Idx.Ey + Idx.n_fields*mode; + auto const Ez_m = Idx.Ez + Idx.n_fields*mode; + auto const Bp_m = Idx.Bx + Idx.n_fields*mode; + auto const Bm_m = Idx.By + Idx.n_fields*mode; + auto const Bz_m = Idx.Bz + Idx.n_fields*mode; + auto const Jp_m = Idx.Jx + Idx.n_fields*mode; + auto const Jm_m = Idx.Jy + Idx.n_fields*mode; + auto const Jz_m = Idx.Jz + Idx.n_fields*mode; + auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode; + auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode; + + // Record old values of the fields to be updated + Complex const Ep_old = fields(i,j,k,Ep_m); + Complex const Em_old = fields(i,j,k,Em_m); + Complex const Ez_old = fields(i,j,k,Ez_m); + Complex const Bp_old = fields(i,j,k,Bp_m); + Complex const Bm_old = fields(i,j,k,Bm_m); + Complex const Bz_old = fields(i,j,k,Bz_m); + // Shortcut for the values of J and rho + Complex const Jp = fields(i,j,k,Jp_m); + Complex const Jm = fields(i,j,k,Jm_m); + Complex const Jz = fields(i,j,k,Jz_m); + Complex const rho_old = fields(i,j,k,rho_old_m); + Complex const rho_new = fields(i,j,k,rho_new_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]; + amrex::Real const kz = modified_kz_arr[j]; + + 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); + Complex const X1 = X1_arr(i,j,k,mode); + Complex const X2 = X2_arr(i,j,k,mode); + Complex const X3 = X3_arr(i,j,k,mode); + Complex const X4 = X4_arr(i,j,k,mode); + Complex const T2 = Theta2_arr(i,j,k,mode); + Complex const T_rho = T_rho_arr(i,j,k,mode); + + Complex rho_diff; + if (update_with_rho) { + rho_diff = X2*rho_new - T2*X3*rho_old; + } else { + Complex const divE = kr*(Ep_old - Em_old) + I*kz*Ez_old; + Complex const divJ = kr*(Jp - Jm) + I*kz*Jz; + + auto const myeps0 = PhysConst::ep0; // temporary for NVCC + rho_diff = T2*(X2 - X3)*myeps0*divE + T_rho*X2*divJ; + } + + // Update E (see WarpX online documentation: theory section) + fields(i,j,k,Ep_m) = T2*C*Ep_old + + T2*S_ck*(-c2*I*kr/2._rt*Bz_old + c2*kz*Bp_old) + + X4*Jp + 0.5_rt*kr*rho_diff; + fields(i,j,k,Em_m) = T2*C*Em_old + + T2*S_ck*(-c2*I*kr/2._rt*Bz_old - c2*kz*Bm_old) + + X4*Jm - 0.5_rt*kr*rho_diff; + fields(i,j,k,Ez_m) = T2*C*Ez_old + + T2*S_ck*(c2*I*kr*Bp_old + c2*I*kr*Bm_old) + + X4*Jz - I*kz*rho_diff; + // Update B (see WarpX online documentation: theory section) + // Note: here X1 is T2*x1/(ep0*c*c*k_norm*k_norm), where + // x1 has the same definition as in the original paper + fields(i,j,k,Bp_m) = T2*C*Bp_old + - T2*S_ck*(-I*kr/2._rt*Ez_old + kz*Ep_old) + + X1*(-I*kr/2._rt*Jz + kz*Jp); + fields(i,j,k,Bm_m) = T2*C*Bm_old + - T2*S_ck*(-I*kr/2._rt*Ez_old - kz*Em_old) + + X1*(-I*kr/2._rt*Jz - kz*Jm); + fields(i,j,k,Bz_m) = T2*C*Bz_old + - T2*S_ck*I*(kr*Ep_old + kr*Em_old) + + X1*I*(kr*Jp + kr*Jm); + }); + } +} + +void PsatdAlgorithmGalileanRZ::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 const& C = C_coef[mfi].array(); + amrex::Array4 const& S_ck = S_ck_coef[mfi].array(); + amrex::Array4 const& X1 = X1_coef[mfi].array(); + amrex::Array4 const& X2 = X2_coef[mfi].array(); + amrex::Array4 const& X3 = X3_coef[mfi].array(); + amrex::Array4 const& X4 = X4_coef[mfi].array(); + amrex::Array4 const& Theta2 = Theta2_coef[mfi].array(); + amrex::Array4 const& T_rho = T_rho_coef[mfi].array(); + + // Extract real (for portability on GPU) + amrex::Real vz = m_v_galilean[2]; + + auto 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 + { + constexpr amrex::Real c = PhysConst::c; + constexpr amrex::Real ep0 = PhysConst::ep0; + Complex const I = Complex{0._rt,1._rt}; + + // 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 + 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); + + // Calculate dot product with galilean velocity + amrex::Real const kv = kz*vz; + + amrex::Real const nu = kv/(k_norm*c); + Complex const theta = amrex::exp( 0.5_rt*I*kv*dt ); + Complex const theta_star = amrex::exp( -0.5_rt*I*kv*dt ); + Complex const e_theta = amrex::exp( I*c*k_norm*dt ); + + Theta2(i,j,k,mode) = theta*theta; + + if (kz == 0._rt) { + T_rho(i,j,k,mode) = -dt; + } else { + T_rho(i,j,k,mode) = (1._rt - theta*theta)/(I*kz*vz); + } + + if ( (nu != 1._rt) && (nu != 0._rt) ) { + + // Note: the coefficients X1, X2, X do not correspond + // exactly to the original Galilean paper, but the + // update equation have been modified accordingly so that + // the expressions below (with the update equations) + // are mathematically equivalent to those of the paper. + Complex x1 = 1._rt/(1._rt-nu*nu) * + (theta_star - C(i,j,k,mode)*theta + I*kv*S_ck(i,j,k,mode)*theta); + // x1, above, is identical to the original paper + X1(i,j,k,mode) = theta*x1/(ep0*c*c*k_norm*k_norm); + // The difference betwen X2 and X3 below, and those + // from the original paper is the factor ep0*k_norm*k_norm + X2(i,j,k,mode) = (x1 - theta*(1._rt - C(i,j,k,mode))) + /(theta_star-theta)/(ep0*k_norm*k_norm); + X3(i,j,k,mode) = (x1 - theta_star*(1._rt - C(i,j,k,mode))) + /(theta_star-theta)/(ep0*k_norm*k_norm); + X4(i,j,k,mode) = I*kv*X1(i,j,k,mode) - theta*theta*S_ck(i,j,k,mode)/ep0; + + } else if (nu == 0._rt) { + + X1(i,j,k,mode) = (1._rt - C(i,j,k,mode))/(ep0 * c*c * k_norm*k_norm); + X2(i,j,k,mode) = (1._rt - S_ck(i,j,k,mode)/dt)/(ep0 * k_norm*k_norm); + X3(i,j,k,mode) = (C(i,j,k,mode) - S_ck(i,j,k,mode)/dt)/(ep0 * k_norm*k_norm); + X4(i,j,k,mode) = -S_ck(i,j,k,mode)/ep0; + + } else if ( nu == 1._rt) { + X1(i,j,k,mode) = (1._rt - e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*c*c*ep0*k_norm*k_norm); + X2(i,j,k,mode) = (3._rt - 4._rt*e_theta + e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*ep0*k_norm*k_norm*(1._rt - e_theta)); + X3(i,j,k,mode) = (3._rt - 2._rt/e_theta - 2._rt*e_theta + e_theta*e_theta - 2._rt*I*c*k_norm*dt) / (4._rt*ep0*(e_theta - 1._rt)*k_norm*k_norm); + X4(i,j,k,mode) = I*(-1._rt + e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*ep0*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; + X1(i,j,k,mode) = 0.5_rt * dt*dt / ep0; + X2(i,j,k,mode) = c*c * dt*dt / (6._rt*ep0); + X3(i,j,k,mode) = - c*c * dt*dt / (3._rt*ep0); + X4(i,j,k,mode) = -dt/ep0; + Theta2(i,j,k,mode) = 1._rt; + } + }); + } +} + +void +PsatdAlgorithmGalileanRZ::CurrentCorrection (const int lev, + SpectralFieldDataRZ& field_data, + std::array,3>& current, + const std::unique_ptr& rho ) +{ + // Profiling + WARPX_PROFILE( "PsatdAlgorithmGalileanRZ::CurrentCorrection" ); + + const SpectralFieldIndex& Idx = m_spectral_index; + + // Forward Fourier transform of J and rho + field_data.ForwardTransform( lev, *current[0], Idx.Jx, + *current[1], Idx.Jy); + field_data.ForwardTransform( lev, *current[2], Idx.Jz, 0); + field_data.ForwardTransform( lev, *rho, Idx.rho_old, 0 ); + field_data.ForwardTransform( lev, *rho, Idx.rho_new, 1 ); + + // Loop over boxes + for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){ + + amrex::Box const & bx = field_data.fields[mfi].box(); + + // Extract arrays for the fields to be updated + amrex::Array4 fields = field_data.fields[mfi].array(); + + // Extract pointers for the k vectors + auto const & kr_modes = field_data.getKrArray(mfi); + amrex::Real const* kr_arr = kr_modes.dataPtr(); + amrex::Real const* modified_kz_arr = modified_kz_vec[mfi].dataPtr(); + int const nr = bx.length(0); + + // Local copy of member variables before GPU loop + amrex::Real vz = m_v_galilean[2]; + amrex::Real const dt = m_dt; + + // Loop over indices within one box + int const modes = field_data.n_rz_azimuthal_modes; + 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 + auto const Jp_m = Idx.Jx + Idx.n_fields*mode; + auto const Jm_m = Idx.Jy + Idx.n_fields*mode; + auto const Jz_m = Idx.Jz + Idx.n_fields*mode; + auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode; + auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode; + + // Shortcuts for the values of J and rho + Complex const Jp = fields(i,j,k,Jp_m); + Complex const Jm = fields(i,j,k,Jm_m); + Complex const Jz = fields(i,j,k,Jz_m); + Complex const rho_old = fields(i,j,k,rho_old_m); + Complex const rho_new = fields(i,j,k,rho_new_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]; + amrex::Real const kz = modified_kz_arr[j]; + amrex::Real const k_norm2 = kr*kr + kz*kz; + + constexpr Complex I = Complex{0._rt,1._rt}; + + // Correct J + if ( k_norm2 != 0._rt ) + { + Complex const theta2 = amrex::exp(I*kz*vz*dt); + Complex const inv_1_T2 = 1._rt/(kz*vz == 0._rt ? 1._rt : 1._rt - theta2); + Complex const j_corr_coef = (kz == 0._rt ? 1._rt/dt : -I*kz*vz*inv_1_T2); + Complex const F = - (j_corr_coef*(rho_new - rho_old*theta2) + I*kz*Jz + kr*(Jp - Jm))/k_norm2; + + fields(i,j,k,Jp_m) += +0.5_rt*kr*F; + fields(i,j,k,Jm_m) += -0.5_rt*kr*F; + fields(i,j,k,Jz_m) += -I*kz*F; + } + }); + } + + // Backward Fourier transform of J + field_data.BackwardTransform( lev, + *current[0], Idx.Jx, + *current[1], Idx.Jy); + field_data.BackwardTransform( lev, *current[2], Idx.Jz, 0 ); + +} + +void +PsatdAlgorithmGalileanRZ::VayDeposition (const int /*lev*/, + SpectralFieldDataRZ& /*field_data*/, + std::array,3>& /*current*/) +{ + amrex::Abort("Vay deposition not implemented in RZ geometry"); +} diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H new file mode 100644 index 000000000..e9636fc32 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H @@ -0,0 +1,147 @@ +/* Copyright 2019 + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PSATD_ALGORITHM_J_LINEAR_IN_TIME_H_ +#define WARPX_PSATD_ALGORITHM_J_LINEAR_IN_TIME_H_ + +#include "FieldSolver/SpectralSolver/SpectralFieldData.H" +#include "FieldSolver/SpectralSolver/SpectralKSpace.H" +#include "SpectralBaseAlgorithm.H" + +#include +#include +#include + +#include + +#include +#include + +#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 PsatdAlgorithmJLinearInTime : public SpectralBaseAlgorithm +{ + public: + + /** + * \brief Constructor of the class PsatdAlgorithmJLinearInTime + * + * \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 + */ + PsatdAlgorithmJLinearInTime ( + 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 + * ( Vay et al, 2013). + * 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,3>& current, + const std::unique_ptr& rho) override final; + + /** + * \brief Virtual function for Vay current deposition in Fourier space + * ( Vay et al, 2013). + * 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,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_J_LINEAR_IN_TIME_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp new file mode 100644 index 000000000..5db997612 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp @@ -0,0 +1,453 @@ +/* Copyright 2019 + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "PsatdAlgorithmJLinearInTime.H" + +#include "Utils/WarpXConst.H" +#include "Utils/WarpX_Complex.H" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#if WARPX_USE_PSATD + +using namespace amrex::literals; + +PsatdAlgorithmJLinearInTime::PsatdAlgorithmJLinearInTime( + 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 +PsatdAlgorithmJLinearInTime::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 fields = f.fields[mfi].array(); + + // These coefficients are always allocated + amrex::Array4 C_arr = C_coef[mfi].array(); + amrex::Array4 S_ck_arr = S_ck_coef[mfi].array(); + amrex::Array4 X1_arr = X1_coef[mfi].array(); + amrex::Array4 X2_arr = X2_coef[mfi].array(); + amrex::Array4 X3_arr = X3_coef[mfi].array(); + + amrex::Array4 X5_arr; + amrex::Array4 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 PsatdAlgorithmJLinearInTime::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 C = C_coef[mfi].array(); + amrex::Array4 S_ck = S_ck_coef[mfi].array(); + amrex::Array4 X1 = X1_coef[mfi].array(); + amrex::Array4 X2 = X2_coef[mfi].array(); + amrex::Array4 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 PsatdAlgorithmJLinearInTime::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 C = C_coef[mfi].array(); + amrex::Array4 S_ck = S_ck_coef[mfi].array(); + + amrex::Array4 X5 = X5_coef[mfi].array(); + amrex::Array4 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 +PsatdAlgorithmJLinearInTime::CurrentCorrection ( + const int lev, + SpectralFieldData& field_data, + std::array,3>& current, + const std::unique_ptr& rho) +{ + // Profiling + BL_PROFILE("PsatdAlgorithmJLinearInTime::CurrentCorrection"); + + amrex::ignore_unused(lev, field_data, current, rho); + amrex::Abort("Current correction not implemented for multi-J PSATD algorithm"); +} + +void +PsatdAlgorithmJLinearInTime::VayDeposition ( + const int lev, + SpectralFieldData& field_data, + std::array,3>& current) +{ + // Profiling + BL_PROFILE("PsatdAlgorithmJLinearInTime::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/SpectralAlgorithms/PsatdAlgorithmMultiJ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmMultiJ.H deleted file mode 100644 index 215e35676..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmMultiJ.H +++ /dev/null @@ -1,147 +0,0 @@ -/* 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 -#include -#include - -#include - -#include -#include - -#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 - * ( Vay et al, 2013). - * 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,3>& current, - const std::unique_ptr& rho) override final; - - /** - * \brief Virtual function for Vay current deposition in Fourier space - * ( Vay et al, 2013). - * 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,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 deleted file mode 100644 index 7ecde8d79..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmMultiJ.cpp +++ /dev/null @@ -1,453 +0,0 @@ -/* 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 -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#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 fields = f.fields[mfi].array(); - - // These coefficients are always allocated - amrex::Array4 C_arr = C_coef[mfi].array(); - amrex::Array4 S_ck_arr = S_ck_coef[mfi].array(); - amrex::Array4 X1_arr = X1_coef[mfi].array(); - amrex::Array4 X2_arr = X2_coef[mfi].array(); - amrex::Array4 X3_arr = X3_coef[mfi].array(); - - amrex::Array4 X5_arr; - amrex::Array4 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 C = C_coef[mfi].array(); - amrex::Array4 S_ck = S_ck_coef[mfi].array(); - amrex::Array4 X1 = X1_coef[mfi].array(); - amrex::Array4 X2 = X2_coef[mfi].array(); - amrex::Array4 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 C = C_coef[mfi].array(); - amrex::Array4 S_ck = S_ck_coef[mfi].array(); - - amrex::Array4 X5 = X5_coef[mfi].array(); - amrex::Array4 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,3>& current, - const std::unique_ptr& 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,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/SpectralAlgorithms/PsatdAlgorithmPml.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.H new file mode 100644 index 000000000..39e99012f --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.H @@ -0,0 +1,91 @@ +/* Copyright 2019 Axel Huebl, Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PSATD_ALGORITHM_PML_H_ +#define WARPX_PSATD_ALGORITHM_PML_H_ + +#include "SpectralBaseAlgorithm.H" + +#include "FieldSolver/SpectralSolver/SpectralFieldData_fwd.H" +#include "FieldSolver/SpectralSolver/SpectralKSpace_fwd.H" + +#include + +#include + +#include +#include + +#if WARPX_USE_PSATD + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + */ +class PsatdAlgorithmPml : public SpectralBaseAlgorithm +{ + public: + PsatdAlgorithmPml(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 dive_cleaning, + const bool divb_cleaning); + + void InitializeSpectralCoefficients( + const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const amrex::Real dt); + + // Redefine functions from base class + virtual void pushSpectralFields(SpectralFieldData& f) const override final; + + /** + * \brief Virtual function for current correction in Fourier space + * ( Vay et al, 2013). + * 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,3>& current, + const std::unique_ptr& rho) override final; + + /** + * \brief Virtual function for Vay current deposition in Fourier space + * ( Vay et al, 2013). + * 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,3>& current) override final; + + private: + SpectralFieldIndex m_spectral_index; + 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 +#endif // WARPX_PSATD_ALGORITHM_PML_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp new file mode 100644 index 000000000..0780b8c52 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp @@ -0,0 +1,421 @@ +/* Copyright 2019 Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "PsatdAlgorithmPml.H" + +#include "FieldSolver/SpectralSolver/SpectralFieldData.H" +#include "FieldSolver/SpectralSolver/SpectralKSpace.H" +#include "Utils/WarpXConst.H" +#include "Utils/WarpX_Complex.H" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#if WARPX_USE_PSATD + +using namespace amrex; + +/* \brief Initialize coefficients for the update equation */ +PsatdAlgorithmPml::PsatdAlgorithmPml(const SpectralKSpace& spectral_kspace, + const 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 Real dt, + const bool dive_cleaning, const bool divb_cleaning) + // Initialize members of base class + : SpectralBaseAlgorithm(spectral_kspace, dm, spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards), + m_spectral_index(spectral_index), + m_dt(dt), + m_dive_cleaning(dive_cleaning), + m_divb_cleaning(divb_cleaning) +{ + const BoxArray& ba = spectral_kspace.spectralspace_ba; + + // Allocate the arrays of coefficients + C_coef = SpectralRealCoefficients(ba, dm, 1, 0); + S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0); + inv_k2_coef = SpectralRealCoefficients(ba, dm, 1, 0); + + InitializeSpectralCoefficients(spectral_kspace, dm, dt); +} + +/* Advance the E and B field in spectral space (stored in `f`) + * over one time step */ +void +PsatdAlgorithmPml::pushSpectralFields(SpectralFieldData& f) const { + + const bool dive_cleaning = m_dive_cleaning; + const bool divb_cleaning = m_divb_cleaning; + + const SpectralFieldIndex& Idx = m_spectral_index; + + // Loop over boxes + for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){ + + const Box& bx = f.fields[mfi].box(); + + // Extract arrays for the fields to be updated + Array4 fields = f.fields[mfi].array(); + + // Extract arrays for the coefficients + Array4 C_arr = C_coef[mfi].array(); + Array4 S_ck_arr = S_ck_coef[mfi].array(); + Array4 inv_k2_arr = inv_k2_coef[mfi].array(); + + // Extract pointers for the k vectors + const Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr(); +#if defined(WARPX_DIM_3D) + const Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr(); +#endif + const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr(); + + const amrex::Real dt = m_dt; + + // Loop over indices within one box + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Record old values of the fields to be updated + const Complex Exy = fields(i,j,k,Idx.Exy); + const Complex Exz = fields(i,j,k,Idx.Exz); + const Complex Eyx = fields(i,j,k,Idx.Eyx); + const Complex Eyz = fields(i,j,k,Idx.Eyz); + const Complex Ezx = fields(i,j,k,Idx.Ezx); + const Complex Ezy = fields(i,j,k,Idx.Ezy); + + const Complex Bxy = fields(i,j,k,Idx.Bxy); + const Complex Bxz = fields(i,j,k,Idx.Bxz); + const Complex Byx = fields(i,j,k,Idx.Byx); + const Complex Byz = fields(i,j,k,Idx.Byz); + 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 defined(WARPX_DIM_3D) + const Real ky = modified_ky_arr[j]; + const Real kz = modified_kz_arr[k]; +#else + constexpr Real ky = 0._rt; + const Real kz = modified_kz_arr[j]; +#endif + constexpr Real c2 = PhysConst::c * PhysConst::c; + + const Complex I = Complex{0._rt, 1._rt}; + + const Real kx2 = kx*kx; + const Real ky2 = ky*ky; + const Real kz2 = kz*kz; + const Real k_norm = std::sqrt(kx2 + ky2 + kz2); + + if (k_norm != 0._rt) { + + 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; + } + } + }); + } +} + +void PsatdAlgorithmPml::InitializeSpectralCoefficients ( + const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const amrex::Real dt) +{ + const BoxArray& ba = spectral_kspace.spectralspace_ba; + + // Fill them with the right values: + // Loop over boxes and allocate the corresponding coefficients + // for each box owned by the local MPI proc + for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi) { + + const Box& bx = ba[mfi]; + + // Extract pointers for the k vectors + const Real* modified_kx = modified_kx_vec[mfi].dataPtr(); +#if defined(WARPX_DIM_3D) + const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); +#endif + const Real* modified_kz = modified_kz_vec[mfi].dataPtr(); + + // Extract arrays for the coefficients + Array4 C = C_coef[mfi].array(); + Array4 S_ck = S_ck_coef[mfi].array(); + Array4 inv_k2 = inv_k2_coef[mfi].array(); + + // Loop over indices within one box + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const Real kx = modified_kx[i]; +#if defined(WARPX_DIM_3D) + const Real ky = modified_ky[j]; + const Real kz = modified_kz[k]; +#else + constexpr Real ky = 0._rt; + const Real kz = modified_kz[j]; +#endif + + // Calculate norm of vector + const Real k_norm = std::sqrt(kx*kx + ky*ky + kz*kz); + const Real k2 = k_norm * k_norm; + + // Calculate coefficients + constexpr Real c = PhysConst::c; + + // Coefficients for k_norm = 0 do not need to be set + if (k_norm != 0._rt) { + C(i,j,k) = std::cos(c * k_norm * dt); + S_ck(i,j,k) = std::sin(c * k_norm * dt) / (c * k_norm); + inv_k2(i,j,k) = 1._rt / k2; + } + }); + } +} + +void +PsatdAlgorithmPml::CurrentCorrection (const int /*lev*/, + SpectralFieldData& /*field_data*/, + std::array,3>& /*current*/, + const std::unique_ptr& /*rho*/) +{ + amrex::Abort("Current correction not implemented for PML PSATD"); +} + +void +PsatdAlgorithmPml::VayDeposition (const int /*lev*/, + SpectralFieldData& /*field_data*/, + std::array,3>& /*current*/) +{ + amrex::Abort("Vay deposition not implemented for PML PSATD"); +} + +#endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.H new file mode 100644 index 000000000..fac84b90a --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.H @@ -0,0 +1,72 @@ +/* Copyright 2021 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PSATD_ALGORITHM_PML_RZ_H_ +#define WARPX_PSATD_ALGORITHM_PML_RZ_H_ + +#include "SpectralBaseAlgorithmRZ.H" + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + */ +class PsatdAlgorithmPmlRZ : public SpectralBaseAlgorithmRZ +{ + + public: + PsatdAlgorithmPmlRZ (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 + * ( Vay et al, 2013). + * 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,3>& current, + const std::unique_ptr& rho) override final; + + /** + * \brief Virtual function for Vay current deposition in Fourier space + * ( Vay et al, 2013). + * 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,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_PSATD_ALGORITHM_PML_RZ_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.cpp new file mode 100644 index 000000000..16f4f1e37 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.cpp @@ -0,0 +1,177 @@ +/* Copyright 2021 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "PsatdAlgorithmPmlRZ.H" +#include "FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.H" +#include "Utils/WarpXConst.H" +#include "Utils/WarpXProfilerWrapper.H" +#include "WarpX.H" + +#include + +using amrex::operator""_rt; + + +/* \brief Initialize coefficients for the update equation */ +PsatdAlgorithmPmlRZ::PsatdAlgorithmPmlRZ (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 +PsatdAlgorithmPmlRZ::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 const& fields = f.fields[mfi].array(); + // Extract arrays for the coefficients + amrex::Array4 const& C_arr = C_coef[mfi].array(); + amrex::Array4 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 PsatdAlgorithmPmlRZ::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 const& C = C_coef[mfi].array(); + amrex::Array4 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 +PsatdAlgorithmPmlRZ::CurrentCorrection (const int /* lev */, + SpectralFieldDataRZ& /* field_data */, + std::array,3>& /* current */, + const std::unique_ptr& /* rho */) +{ + amrex::Abort("Current correction not implemented in RZ geometry PML"); +} + +void +PsatdAlgorithmPmlRZ::VayDeposition (const int /* lev */, + SpectralFieldDataRZ& /*field_data*/, + std::array,3>& /*current*/) +{ + amrex::Abort("Vay deposition not implemented in RZ geometry PML"); +} -- cgit v1.2.3