aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/Evolve/WarpXEvolve.cpp43
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H19
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp96
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H16
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralSolver.H15
-rw-r--r--Source/WarpX.H14
-rw-r--r--Source/WarpX.cpp1
7 files changed, 198 insertions, 6 deletions
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
index 8ca059184..7cd1abd8a 100644
--- a/Source/Evolve/WarpXEvolve.cpp
+++ b/Source/Evolve/WarpXEvolve.cpp
@@ -16,6 +16,9 @@
#ifdef WARPX_USE_PY
# include "Python/WarpX_py.H"
#endif
+#ifdef WARPX_USE_PSATD
+#include "FieldSolver/SpectralSolver/SpectralSolver.H"
+#endif
#ifdef BL_USE_SENSEI_INSITU
# include <AMReX_AmrMeshInSituBridge.H>
@@ -341,15 +344,36 @@ WarpX::OneStep_nosub (Real cur_time)
if (warpx_py_afterdeposition) warpx_py_afterdeposition();
#endif
+#ifdef WARPX_USE_PSATD
+ // Apply current correction in Fourier space
+ // (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010)
+ if ( fft_periodic_single_box == false ) {
+ // For domain decomposition with local FFT over guard cells,
+ // apply this before `SyncCurrent`, i.e. before exchanging guard cells for J
+ if ( do_current_correction ) CurrentCorrection();
+ }
+#endif
+
#ifdef WARPX_QED
//Do QED processes
mypc->doQedEvents();
#endif
+ // Synchronize J and rho
SyncCurrent();
-
SyncRho();
+#ifdef WARPX_USE_PSATD
+ // Apply current correction in Fourier space
+ // (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010)
+ if ( fft_periodic_single_box == true ) {
+ // For periodic, single-box FFT (FFT without guard cells)
+ // apply this after `SyncCurrent`, i.e. after exchanging guard cells for J
+ if ( do_current_correction ) CurrentCorrection();
+ }
+#endif
+
+
// At this point, J is up-to-date inside the domain, and E and B are
// up-to-date including enough guard cells for first step of the field
// solve.
@@ -785,3 +809,20 @@ WarpX::applyMirrors(Real time){
}
}
}
+
+#ifdef WARPX_USE_PSATD
+void
+WarpX::CurrentCorrection ()
+{
+ for ( int lev = 0; lev <= finest_level; ++lev )
+ {
+ // Apply correction on fine patch
+ spectral_solver_fp[lev]->CurrentCorrection( current_fp[lev], rho_fp[lev] );
+ if ( spectral_solver_cp[lev] )
+ {
+ // Apply correction on coarse patch
+ spectral_solver_cp[lev]->CurrentCorrection( current_cp[lev], rho_cp[lev] );
+ }
+ }
+}
+#endif
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
index d71c0ab18..eb440d118 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
@@ -1,4 +1,4 @@
-/* Copyright 2019 Maxence Thevenet, Remi Lehe, Revathi Jambunathan
+/* Copyright 2019 Maxence Thevenet, Remi Lehe, Revathi Jambunathan, Edoardo Zoni
*
*
* This file is part of WarpX.
@@ -36,8 +36,25 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
const amrex::DistributionMapping& dm,
const amrex::Real dt);
+ /**
+ * \brief Virtual function for current correction in Fourier space
+ * (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010).
+ * This function overrides the virtual function \c CurrentCorrection in the
+ * base class \c SpectralBaseAlgorithm (qualifier \c override) and cannot be
+ * overridden by further derived classes (qualifier \c final).
+ *
+ * \param[in,out] field_data all fields in Fourier space
+ * \param[in,out] current three-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 ( SpectralFieldData& field_data,
+ std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+ const std::unique_ptr<amrex::MultiFab>& rho ) override final;
+
private:
SpectralRealCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef;
+ amrex::Real m_dt;
};
#endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
index 8fee0967d..9684fde27 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
@@ -1,11 +1,13 @@
-/* Copyright 2019 Remi Lehe, Revathi Jambunathan
+/* Copyright 2019 Remi Lehe, Revathi Jambunathan, Edoardo Zoni
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
+#include "WarpX.H"
#include "PsatdAlgorithm.H"
#include "Utils/WarpXConst.H"
+#include "Utils/WarpXProfilerWrapper.H"
#include <cmath>
@@ -13,7 +15,9 @@
#if WARPX_USE_PSATD
using namespace amrex;
-/* \brief Initialize coefficients for the update equation */
+/**
+ * \brief Constructor
+ */
PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace,
const DistributionMapping& dm,
const int norder_x, const int norder_y,
@@ -31,11 +35,15 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace,
X2_coef = SpectralRealCoefficients(ba, dm, 1, 0);
X3_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+ // Initialize coefficients for update equations
InitializeSpectralCoefficients(spectral_kspace, dm, dt);
+
+ m_dt = dt;
}
-/* Advance the E and B field in spectral space (stored in `f`)
- * over one time step */
+/**
+ * \brief Advance E and B fields in spectral space (stored in `f`) over one time step
+ */
void
PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
@@ -120,6 +128,9 @@ PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
}
};
+/**
+ * \brief Initialize coefficients for update equations
+ */
void PsatdAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace,
const amrex::DistributionMapping& dm,
const amrex::Real dt)
@@ -179,4 +190,81 @@ void PsatdAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spectr
});
}
}
+
+void
+PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data,
+ std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+ const std::unique_ptr<amrex::MultiFab>& rho ) {
+ // Profiling
+ WARPX_PROFILE( "PsatdAlgorithm::CurrentCorrection" );
+
+ using Idx = SpectralFieldIndex;
+
+ // Forward Fourier transform of J and rho
+ field_data.ForwardTransform( *current[0], Idx::Jx, 0 );
+ field_data.ForwardTransform( *current[1], Idx::Jy, 0 );
+ field_data.ForwardTransform( *current[2], Idx::Jz, 0 );
+ field_data.ForwardTransform( *rho, Idx::rho_old, 0 );
+ field_data.ForwardTransform( *rho, Idx::rho_new, 1 );
+
+ // Loop over boxes
+ for (MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
+
+ const Box& bx = field_data.fields[mfi].box();
+
+ // Extract arrays for the fields to be updated
+ Array4<Complex> fields = field_data.fields[mfi].array();
+ // Extract pointers for the k vectors
+ const Real* const modified_kx_arr = modified_kx_vec[mfi].dataPtr();
+#if (AMREX_SPACEDIM==3)
+ const Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr();
+#endif
+ const Real* const modified_kz_arr = modified_kz_vec[mfi].dataPtr();
+
+ // Local copy of member variables before GPU loop
+ const 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
+ using Idx = SpectralFieldIndex;
+ // 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 Real kx = modified_kx_arr[i];
+#if (AMREX_SPACEDIM==3)
+ const Real ky = modified_ky_arr[j];
+ const Real kz = modified_kz_arr[k];
+#else
+ constexpr Real ky = 0;
+ const Real kz = modified_kz_arr[j];
+#endif
+ constexpr Complex I = Complex{0,1};
+
+ const Real k_norm = std::sqrt( kx*kx + ky*ky + kz*kz );
+
+ // div(J) in Fourier space
+ const Complex k_dot_J = kx*Jx + ky*Jy + kz*Jz;
+
+ // Correct J
+ if ( k_norm != 0 )
+ {
+ fields(i,j,k,Idx::Jx) = Jx - (k_dot_J-I*(rho_new-rho_old)/dt)*kx/(k_norm*k_norm);
+ fields(i,j,k,Idx::Jy) = Jy - (k_dot_J-I*(rho_new-rho_old)/dt)*ky/(k_norm*k_norm);
+ fields(i,j,k,Idx::Jz) = Jz - (k_dot_J-I*(rho_new-rho_old)/dt)*kz/(k_norm*k_norm);
+ }
+ });
+ }
+
+ // Backward Fourier transform of J
+ field_data.BackwardTransform( *current[0], Idx::Jx, 0 );
+ field_data.BackwardTransform( *current[1], Idx::Jy, 0 );
+ field_data.BackwardTransform( *current[2], Idx::Jz, 0 );
+}
#endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
index 7ce774ab2..7143b7aa9 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
@@ -32,6 +32,22 @@ class SpectralBaseAlgorithm
virtual ~SpectralBaseAlgorithm() {};
/**
+ * \brief Virtual function for current correction in Fourier space
+ * (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010).
+ * This virtual function is not pure: it can be overridden in derived
+ * classes (e.g. PsatdAlgorithm, PMLPsatdAlgorithm), but a base-class
+ * implementation is provided (empty implementation in this case).
+ *
+ * \param[in,out] field_data all fields in Fourier space
+ * \param[in,out] current three-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 ( SpectralFieldData& field_data,
+ std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+ const std::unique_ptr<amrex::MultiFab>& rho ) {};
+
+ /**
* \brief Compute spectral divergence of E
*/
void ComputeSpectralDivE ( SpectralFieldData& field_data,
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
index 08eef19ad..6a6410d49 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
@@ -67,6 +67,21 @@ class SpectralSolver
algorithm->ComputeSpectralDivE( field_data, Efield, divE );
};
+ /**
+ * \brief Public interface to call the virtual function \c CurrentCorrection,
+ * defined in the base class SpectralBaseAlgorithm and possibly overridden
+ * by its derived classes (e.g. PsatdAlgorithm, PMLPsatdAlgorithm), from
+ * objects of class SpectralSolver through the private unique pointer \c algorithm
+ *
+ * \param[in,out] current three-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
+ */
+ void CurrentCorrection ( std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+ const std::unique_ptr<amrex::MultiFab>& rho ) {
+ algorithm->CurrentCorrection( field_data, current, rho );
+ };
+
private:
// Store field in spectral space and perform the Fourier transforms
diff --git a/Source/WarpX.H b/Source/WarpX.H
index a19d9bf81..26b96ecfd 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -140,6 +140,11 @@ public:
static int maxwell_fdtd_solver_id;
static long load_balance_costs_update_algo;
+ // Static public variable to switch on current correction in Fourier space
+ // (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010): can be set
+ // true by the user as an input parameter
+ bool do_current_correction = false;
+
// div E cleaning
static int do_dive_cleaning;
@@ -552,6 +557,15 @@ private:
void AddRhoFromFineLevelandSumBoundary (int lev, int icomp, int ncomp);
void NodalSyncRho (int lev, PatchType patch_type, int icomp, int ncomp);
+ /**
+ * \brief Private function for current correction in Fourier space
+ * (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010):
+ * loops over the MR levels and applies the correction on the fine and coarse
+ * patches (calls the virtual method \c CurrentCorrection of the spectral
+ * algorithm in use, via the public interface defined in the class SpectralSolver).
+ */
+ void CurrentCorrection ();
+
void ReadParameters ();
void InitFromScratch ();
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 17a47caaa..5962394e6 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -635,6 +635,7 @@ WarpX::ReadParameters ()
pp.query("nox", nox_fft);
pp.query("noy", noy_fft);
pp.query("noz", noz_fft);
+ pp.query("do_current_correction", do_current_correction);
pp.query("v_galilean", v_galilean);
// Scale the velocity by the speed of light
for (int i=0; i<3; i++) v_galilean[i] *= PhysConst::c;