diff options
| -rw-r--r-- | Docs/source/running_cpp/parameters.rst | 5 | ||||
| -rwxr-xr-x | Examples/Tests/Langmuir/analysis_langmuir_multi.py | 14 | ||||
| -rwxr-xr-x | Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py | 14 | ||||
| -rw-r--r-- | Regression/WarpX-tests.ini | 37 | ||||
| -rw-r--r-- | Source/Evolve/WarpXEvolve.cpp | 43 | ||||
| -rw-r--r-- | Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H | 19 | ||||
| -rw-r--r-- | Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp | 96 | ||||
| -rw-r--r-- | Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H | 16 | ||||
| -rw-r--r-- | Source/FieldSolver/SpectralSolver/SpectralSolver.H | 15 | ||||
| -rw-r--r-- | Source/WarpX.H | 14 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 1 |
11 files changed, 268 insertions, 6 deletions
diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 2afb42772..b96e9c471 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -1028,6 +1028,11 @@ Numerics and algorithms See `this section of the FFTW documentation <http://www.fftw.org/fftw3_doc/Planner-Flags.html>`__ for more information. +* ``psatd.do_current_correction`` (`0` or `1`; default: `0`) + If true, the current correction defined by equation (19) of + `(Vay et al, JCP 243, 2013) <https://doi.org/10.1016/j.jcp.2013.03.010>`_ is applied. + Only used when compiled and running with the PSATD solver. + * ``pstad.v_galilean`` (`3 floats`, in units of the speed of light; default `0. 0. 0.`) Defines the galilean velocity. Non-zero `v_galilean` activates Galilean algorithm, which suppresses the Numerical Cherenkov instability diff --git a/Examples/Tests/Langmuir/analysis_langmuir_multi.py b/Examples/Tests/Langmuir/analysis_langmuir_multi.py index 9a2335bfe..e19f855cf 100755 --- a/Examples/Tests/Langmuir/analysis_langmuir_multi.py +++ b/Examples/Tests/Langmuir/analysis_langmuir_multi.py @@ -15,6 +15,7 @@ # $$ E_y = \epsilon \,\frac{m_e c^2 k_y}{q_e}\cos(k_x x)\sin(k_y y)\cos(k_z z)\sin( \omega_p t)$$ # $$ E_z = \epsilon \,\frac{m_e c^2 k_z}{q_e}\cos(k_x x)\cos(k_y y)\sin(k_z z)\sin( \omega_p t)$$ import sys +import re import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt @@ -26,6 +27,9 @@ from scipy.constants import e, m_e, epsilon_0, c # this will be the name of the plot file fn = sys.argv[1] +# Parse test name and check if current correction (psatd.do_current_correction=1) is applied +current_correction = True if re.search( 'current_correction', fn ) else False + # Parameters (these parameters must match the parameters in `inputs.multi.rt`) epsilon = 0.01 n = 4.e24 @@ -122,3 +126,13 @@ print("error_rel : " + str(error_rel)) print("tolerance_rel: " + str(tolerance_rel)) assert( error_rel < tolerance_rel ) + +# Check relative L-infinity spatial norm of rho/epsilon_0 - div(E) when +# current correction (psatd.do_current_correction=1) is applied +if current_correction: + rho = data['rho' ].to_ndarray() + divE = data['divE'].to_ndarray() + Linf_norm = np.amax( np.abs( rho/epsilon_0 - divE ) ) / np.amax( np.abs( rho/epsilon_0 ) ) + print("error: " + str(Linf_norm)) + print("tolerance: 1.e-9") + assert( Linf_norm < 1.e-9 ) diff --git a/Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py b/Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py index 3a225fe0f..1e417b9b8 100755 --- a/Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py +++ b/Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py @@ -15,6 +15,7 @@ # $$ E_y = \epsilon \,\frac{m_e c^2 k_y}{q_e}\cos(k_x x)\sin(k_y y)\cos(k_z z)\sin( \omega_p t)$$ # $$ E_z = \epsilon \,\frac{m_e c^2 k_z}{q_e}\cos(k_x x)\cos(k_y y)\sin(k_z z)\sin( \omega_p t)$$ import sys +import re import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt @@ -26,6 +27,9 @@ from scipy.constants import e, m_e, epsilon_0, c # this will be the name of the plot file fn = sys.argv[1] +# Parse test name and check if current correction (psatd.do_current_correction=1) is applied +current_correction = True if re.search( 'current_correction', fn ) else False + # Parameters (these parameters must match the parameters in `inputs.multi.rt`) epsilon = 0.01 n = 4.e24 @@ -95,3 +99,13 @@ print("error_rel : " + str(error_rel)) print("tolerance_rel: " + str(tolerance_rel)) assert( error_rel < tolerance_rel ) + +# Check relative L-infinity spatial norm of rho/epsilon_0 - div(E) when +# current correction (psatd.do_current_correction=1) is applied +if current_correction: + rho = data['rho' ].to_ndarray() + divE = data['divE'].to_ndarray() + Linf_norm = np.amax( np.abs( rho/epsilon_0 - divE ) ) / np.amax( np.abs( rho/epsilon_0 ) ) + print("error: " + str(Linf_norm)) + print("tolerance: 1.e-9") + assert( Linf_norm < 1.e-9 ) diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index f2807124b..553dbaeb4 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -379,6 +379,25 @@ analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi.py analysisOutputImage = langmuir_multi_analysis.png tolerance = 5.e-11 +[Langmuir_multi_psatd_current_correction] +buildDir = . +inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt +runtime_params = psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 psatd.do_current_correction=1 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz part_per_cell rho divE +dim = 3 +addToCompileString = USE_PSATD=TRUE +restartTest = 0 +useMPI = 1 +numprocs = 1 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +tolerance = 5.e-11 +particleTypes = electrons positrons +analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi.py +analysisOutputImage = langmuir_multi_analysis.png + [Langmuir_multi_psatd_nodal] buildDir = . inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt @@ -455,6 +474,24 @@ analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py analysisOutputImage = langmuir_multi_2d_analysis.png tolerance = 1.e-14 +[Langmuir_multi_2d_psatd_current_correction] +buildDir = . +inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt +runtime_params = amr.max_grid_size=128 psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 psatd.do_current_correction=1 diag1.electrons.variables=w ux uy uz Ex Ey Ez diag1.positrons.variables=w ux uy uz Ex Ey Ez diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE +dim = 2 +addToCompileString = USE_PSATD=TRUE +restartTest = 0 +useMPI = 1 +numprocs = 1 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = electrons positrons +analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py +analysisOutputImage = langmuir_multi_2d_analysis.png + [Langmuir_multi_2d_psatd_nodal] buildDir = . inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt 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; |