diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/Evolve/WarpXEvolve.cpp | 56 | ||||
| -rw-r--r-- | Source/FieldSolver/WarpXPushFieldsEM.cpp | 133 | ||||
| -rw-r--r-- | Source/WarpX.H | 16 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 18 |
4 files changed, 176 insertions, 47 deletions
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index b7a8c64f1..cc71e78d5 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -398,19 +398,9 @@ WarpX::OneStep_nosub (Real cur_time) ExecutePythonCallback("afterdeposition"); - // Synchronize J and rho: filter, exchange boundary, interpolate across levels. - // With Vay current deposition, the current deposited at this point is not yet - // the actual current J. This is computed later in WarpX::PushPSATD, by calling - // WarpX::PSATDVayDeposition. The function SyncCurrent is called after that, - // instead of here, so that we synchronize the correct current. - // With current centering, the nodal current is deposited in 'current_fp_nodal': - // SyncCurrent stores the result of its centering into 'current_fp' and then - // performs both filtering, if used, and exchange of guard cells. - if (WarpX::current_deposition_algo != CurrentDepositionAlgo::Vay) - { - SyncCurrent(current_fp, current_cp); - } - SyncRho(); + // Synchronize J and rho: + // filter (if used), exchange guard cells, interpolate across MR levels + SyncCurrentAndRho(); // 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 @@ -495,6 +485,46 @@ WarpX::OneStep_nosub (Real cur_time) ExecutePythonCallback("afterEsolve"); } +void WarpX::SyncCurrentAndRho () +{ + if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) + { + if (fft_periodic_single_box) + { + // With periodic single box, synchronize J and rho here, + // even with current correction or Vay deposition + if (current_deposition_algo == CurrentDepositionAlgo::Vay) + { + // TODO Replace current_cp with current_cp_vay once Vay deposition is implemented with MR + SyncCurrent(current_fp_vay, current_cp); + SyncRho(); + } + else + { + SyncCurrent(current_fp, current_cp); + SyncRho(); + } + } + else // no periodic single box + { + // Without periodic single box, synchronize J and rho here, + // except with current correction or Vay deposition: + // in these cases, synchronize later (in WarpX::PushPSATD) + if (current_correction == false && + current_deposition_algo != CurrentDepositionAlgo::Vay) + { + SyncCurrent(current_fp, current_cp); + SyncRho(); + } + } + } + else // FDTD + { + SyncCurrent(current_fp, current_cp); + SyncRho(); + } +} + void WarpX::OneStep_multiJ (const amrex::Real cur_time) { diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 6d604bc15..7499ee140 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -270,7 +270,8 @@ WarpX::PSATDBackwardTransformG () void WarpX::PSATDForwardTransformJ ( const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_fp, - const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_cp) + const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_cp, + const bool apply_kspace_filter) { SpectralFieldIndex Idx; int idx_jx, idx_jy, idx_jz; @@ -299,7 +300,7 @@ void WarpX::PSATDForwardTransformJ ( #ifdef WARPX_DIM_RZ // Apply filter in k space if needed - if (WarpX::use_kspace_filter) + if (use_kspace_filter && apply_kspace_filter) { for (int lev = 0; lev <= finest_level; ++lev) { @@ -311,6 +312,8 @@ void WarpX::PSATDForwardTransformJ ( } } } +#else + amrex::ignore_unused(apply_kspace_filter); #endif } @@ -349,8 +352,10 @@ void WarpX::PSATDBackwardTransformJ ( void WarpX::PSATDForwardTransformRho ( const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& charge_fp, const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& charge_cp, - const int icomp, const int dcomp) + const int icomp, const int dcomp, const bool apply_kspace_filter) { + if (charge_fp[0] == nullptr) return; + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; // Select index in k space @@ -368,7 +373,7 @@ void WarpX::PSATDForwardTransformRho ( #ifdef WARPX_DIM_RZ // Apply filter in k space if needed - if (WarpX::use_kspace_filter) + if (use_kspace_filter && apply_kspace_filter) { for (int lev = 0; lev <= finest_level; ++lev) { @@ -380,6 +385,8 @@ void WarpX::PSATDForwardTransformRho ( } } } +#else + amrex::ignore_unused(apply_kspace_filter); #endif } @@ -635,46 +642,114 @@ WarpX::PushPSATD () "PushFieldsEM: PSATD solver selected but not built")); #else - PSATDForwardTransformEB(Efield_fp, Bfield_fp, Efield_cp, Bfield_cp); - - amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_fp = - (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ? current_fp_vay : current_fp; + if (fft_periodic_single_box) + { + if (current_correction) + { + // FFT of J and rho + PSATDForwardTransformJ(current_fp, current_cp); + PSATDForwardTransformRho(rho_fp, rho_cp, 0, 0); // rho old + PSATDForwardTransformRho(rho_fp, rho_cp, 1, 1); // rho new - PSATDForwardTransformJ(J_fp, current_cp); + // Correct J in k-space + PSATDCurrentCorrection(); - // Do rho FFTs only if needed - if (WarpX::update_with_rho || WarpX::current_correction || WarpX::do_dive_cleaning) - { - PSATDForwardTransformRho(rho_fp, rho_cp, 0,0); // rho old - PSATDForwardTransformRho(rho_fp, rho_cp, 1,1); // rho new + // Inverse FFT of J + PSATDBackwardTransformJ(current_fp, current_cp); + } + else if (current_deposition_algo == CurrentDepositionAlgo::Vay) + { + // FFT of D and rho (if used) + // TODO Replace current_cp with current_cp_vay once Vay deposition is implemented with MR + PSATDForwardTransformJ(current_fp_vay, current_cp); + PSATDForwardTransformRho(rho_fp, rho_cp, 0, 0); // rho old + PSATDForwardTransformRho(rho_fp, rho_cp, 1, 1); // rho new + + // Compute J from D in k-space + PSATDVayDeposition(); + + // Inverse FFT of J, subtract cumulative sums of D + PSATDBackwardTransformJ(current_fp, current_cp); + // TODO Cumulative sums need to be fixed with periodic single box + PSATDSubtractCurrentPartialSumsAvg(); + + // FFT of J after subtraction of cumulative sums + PSATDForwardTransformJ(current_fp, current_cp); + } + else // no current correction, no Vay deposition + { + // FFT of J and rho (if used) + PSATDForwardTransformJ(current_fp, current_cp); + PSATDForwardTransformRho(rho_fp, rho_cp, 0, 0); // rho old + PSATDForwardTransformRho(rho_fp, rho_cp, 1, 1); // rho new + } } - - // Correct the current in Fourier space so that the continuity equation is satisfied, and - // transform back to real space so that the current correction is reflected in the diagnostics - if (WarpX::current_correction) + else // no periodic single box { - PSATDCurrentCorrection(); - PSATDBackwardTransformJ(current_fp, current_cp); - } + if (current_correction) + { + // FFT of J and rho +#ifdef WARPX_DIM_RZ + // In RZ geometry, do not apply filtering here, since it is + // applied in the subsequent calls to these functions (below) + const bool apply_kspace_filter = false; + PSATDForwardTransformJ(current_fp, current_cp, apply_kspace_filter); + PSATDForwardTransformRho(rho_fp, rho_cp, 0, 0, apply_kspace_filter); // rho old + PSATDForwardTransformRho(rho_fp, rho_cp, 1, 1, apply_kspace_filter); // rho new +#else + PSATDForwardTransformJ(current_fp, current_cp); + PSATDForwardTransformRho(rho_fp, rho_cp, 0, 0); // rho old + PSATDForwardTransformRho(rho_fp, rho_cp, 1, 1); // rho new +#endif - // Compute the current in Fourier space according to the Vay deposition scheme, and - // transform back to real space so that the Vay deposition is reflected in the diagnostics - if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) - { - PSATDVayDeposition(); - PSATDBackwardTransformJ(current_fp, current_cp); - PSATDSubtractCurrentPartialSumsAvg(); - SyncCurrent(current_fp, current_cp); + // Correct J in k-space + PSATDCurrentCorrection(); + + // Inverse FFT of J + PSATDBackwardTransformJ(current_fp, current_cp); + + // Synchronize J and rho + SyncCurrent(current_fp, current_cp); + SyncRho(); + } + else if (current_deposition_algo == CurrentDepositionAlgo::Vay) + { + // FFT of D + PSATDForwardTransformJ(current_fp_vay, current_cp); + + // Compute J from D in k-space + PSATDVayDeposition(); + + // Inverse FFT of J, subtract cumulative sums of D + PSATDBackwardTransformJ(current_fp, current_cp); + PSATDSubtractCurrentPartialSumsAvg(); + + // Synchronize J and rho (if used) + SyncCurrent(current_fp, current_cp); + SyncRho(); + } + + // FFT of J and rho (if used) PSATDForwardTransformJ(current_fp, current_cp); + PSATDForwardTransformRho(rho_fp, rho_cp, 0, 0); // rho old + PSATDForwardTransformRho(rho_fp, rho_cp, 1, 1); // rho new } + // FFT of E and B + PSATDForwardTransformEB(Efield_fp, Bfield_fp, Efield_cp, Bfield_cp); + #ifdef WARPX_DIM_RZ if (pml_rz[0]) pml_rz[0]->PushPSATD(0); #endif + // FFT of F and G if (WarpX::do_dive_cleaning) PSATDForwardTransformF(); if (WarpX::do_divb_cleaning) PSATDForwardTransformG(); + + // Update E, B, F, and G in k-space PSATDPushSpectralFields(); + + // Inverse FFT of E, B, F, and G PSATDBackwardTransformEB(Efield_fp, Bfield_fp, Efield_cp, Bfield_cp); if (WarpX::fft_do_time_averaging) PSATDBackwardTransformEBavg(Efield_avg_fp, Bfield_avg_fp, Efield_avg_cp, Bfield_avg_cp); diff --git a/Source/WarpX.H b/Source/WarpX.H index c16f30408..7a051e30d 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -647,6 +647,13 @@ public: void FillBoundaryAux (int lev, amrex::IntVect ng); /** + * \brief Synchronize J and rho: + * filter (if used), exchange guard cells, interpolate across MR levels. + * Contains separate calls to WarpX::SyncCurrent and WarpX::SyncRho. + */ + void SyncCurrentAndRho (); + + /** * \brief Apply filter and sum guard cells across MR levels. * If current centering is used, center the current from a nodal grid * to a staggered grid. For each MR level beyond level 0, interpolate @@ -1506,10 +1513,13 @@ private: * storing the fine patch current to be transformed * \param J_cp Vector of three-dimensional arrays (for each level) * storing the coarse patch current to be transformed + * \param[in] apply_kspace_filter Control whether to apply filtering + * (only used in RZ geometry to avoid double filtering) */ void PSATDForwardTransformJ ( const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_fp, - const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_cp); + const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_cp, + const bool apply_kspace_filter=true); /** * \brief Backward FFT of J on all mesh refinement levels @@ -1531,11 +1541,13 @@ private: * \param charge_cp Vector (for each level) storing the coarse patch charge to be transformed * \param[in] icomp index of fourth component (0 for rho_old, 1 for rho_new) * \param[in] dcomp index of spectral component (0 for rho_old, 1 for rho_new) + * \param[in] apply_kspace_filter Control whether to apply filtering + * (only used in RZ geometry to avoid double filtering) */ void PSATDForwardTransformRho ( const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& charge_fp, const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& charge_cp, - const int icomp, const int dcomp); + const int icomp, const int dcomp, const bool apply_kspace_filter=true); /** * \brief Copy rho_new to rho_old in spectral space diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 12747b7b9..86961dcaf 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1194,6 +1194,13 @@ WarpX::ReadParameters () "Option algo.current_deposition=vay must be used with psatd.periodic_single_box_fft=0."); } + if (current_deposition_algo == CurrentDepositionAlgo::Vay) + { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + current_correction == false, + "Options algo.current_deposition=vay and psatd.current_correction=1 cannot be combined together."); + } + // Auxiliary: boosted_frame = true if warpx.gamma_boost is set in the inputs amrex::ParmParse pp_warpx("warpx"); const bool boosted_frame = pp_warpx.query("gamma_boost", gamma_boost); @@ -1328,10 +1335,15 @@ WarpX::ReadParameters () } } - // Fill guard cells with backward FFTs if Vay current deposition is used - if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) + // Without periodic single box, fill guard cells with backward FFTs, + // with current correction or Vay deposition + if (fft_periodic_single_box == false) { - WarpX::m_fill_guards_current = amrex::IntVect(1); + if (current_correction || + current_deposition_algo == CurrentDepositionAlgo::Vay) + { + WarpX::m_fill_guards_current = amrex::IntVect(1); + } } } |