diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/Evolve/WarpXEvolve.cpp | 23 | ||||
| -rw-r--r-- | Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp | 18 | ||||
| -rw-r--r-- | Source/WarpX.H | 9 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 2 |
4 files changed, 24 insertions, 28 deletions
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 72621e52c..d53c5ea5b 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -288,14 +288,13 @@ WarpX::OneStep_nosub (Real cur_time) if (warpx_py_afterdeposition) warpx_py_afterdeposition(); #endif +// TODO +// Apply current correction in Fourier space: for domain decomposition with local +// FFTs over guard cells, apply this before calling SyncCurrent #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(); - } + if ( !fft_periodic_single_box && current_correction ) + amrex::Abort("\nCurrent correction does not guarantee charge conservation with local FFTs over guard cells:\n" + "set psatd.periodic_single_box_fft=1 too, in order to guarantee charge conservation"); #endif #ifdef WARPX_QED @@ -306,14 +305,10 @@ WarpX::OneStep_nosub (Real cur_time) SyncCurrent(); SyncRho(); +// Apply current correction in Fourier space: for periodic single-box global FFTs +// without guard cells, apply this after calling SyncCurrent #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(); - } + if ( fft_periodic_single_box && current_correction ) CurrentCorrection(); #endif diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index 9684fde27..c75e8b8c7 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -23,8 +23,8 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const int norder_x, const int norder_y, const int norder_z, const bool nodal, const Real dt) // Initialize members of base class - : SpectralBaseAlgorithm( spectral_kspace, dm, - norder_x, norder_y, norder_z, nodal ) + : m_dt( dt ), + SpectralBaseAlgorithm( spectral_kspace, dm, norder_x, norder_y, norder_z, nodal ) { const BoxArray& ba = spectral_kspace.spectralspace_ba; @@ -37,8 +37,6 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, // Initialize coefficients for update equations InitializeSpectralCoefficients(spectral_kspace, dm, dt); - - m_dt = dt; } /** @@ -214,6 +212,7 @@ PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data, // 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) @@ -225,17 +224,18 @@ PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data, const Real dt = m_dt; // Loop over indices within one box - ParallelFor(bx, - [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + 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) @@ -245,10 +245,10 @@ PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data, 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 ); + constexpr Complex I = Complex{0,1}; + // div(J) in Fourier space const Complex k_dot_J = kx*Jx + ky*Jy + kz*Jz; @@ -259,7 +259,7 @@ PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data, 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 diff --git a/Source/WarpX.H b/Source/WarpX.H index 715e55a4f..ed373d2f0 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -136,10 +136,11 @@ public: static int em_solver_medium; static int macroscopic_solver_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; +#ifdef WARPX_USE_PSATD + // If true (overwritten by the user in the input file), the current correction + // defined in equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010 is applied + bool current_correction = false; +#endif // div E cleaning static int do_dive_cleaning; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 1c5002076..9999dd40b 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -605,7 +605,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("current_correction", 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; |