aboutsummaryrefslogtreecommitdiff
path: root/Source/Evolve
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Evolve')
-rw-r--r--Source/Evolve/WarpXComputeDt.cpp61
-rw-r--r--Source/Evolve/WarpXEvolve.cpp99
2 files changed, 86 insertions, 74 deletions
diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp
index a58fa5e7b..6237f85c4 100644
--- a/Source/Evolve/WarpXComputeDt.cpp
+++ b/Source/Evolve/WarpXComputeDt.cpp
@@ -23,42 +23,37 @@ WarpX::ComputeDt ()
const amrex::Real* dx = geom[max_level].CellSize();
amrex::Real deltat = 0.;
-#if WARPX_USE_PSATD
- // Computation of dt for spectral algorithm
-
-# if (defined WARPX_DIM_RZ)
- // - In RZ geometry: dz/c
- deltat = cfl * dx[1]/PhysConst::c;
-# elif (defined WARPX_DIM_XZ)
- // - In Cartesian 2D geometry: determined by the minimum cell size in all direction
- deltat = cfl * std::min( dx[0], dx[1] )/PhysConst::c;
-# else
- // - In Cartesian 3D geometry: determined by the minimum cell size in all direction
- deltat = cfl * std::min( dx[0], std::min( dx[1], dx[2] ) )/PhysConst::c;
-# endif
-
-
+ if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+ // Computation of dt for spectral algorithm
+#if (defined WARPX_DIM_RZ)
+ // - In RZ geometry: dz/c
+ deltat = cfl * dx[1]/PhysConst::c;
+#elif (defined WARPX_DIM_XZ)
+ // - In Cartesian 2D geometry: determined by the minimum cell size in all direction
+ deltat = cfl * std::min( dx[0], dx[1] )/PhysConst::c;
#else
- // Computation of dt for FDTD algorithm
-
-# ifdef WARPX_DIM_RZ
- // - In RZ geometry
- if (maxwell_solver_id == MaxwellSolverAlgo::Yee) {
- deltat = cfl * CylindricalYeeAlgorithm::ComputeMaxDt(dx, n_rz_azimuthal_modes);
-# else
- // - In Cartesian geometry
- if (do_nodal) {
- deltat = cfl * CartesianNodalAlgorithm::ComputeMaxDt( dx );
- } else if (maxwell_solver_id == MaxwellSolverAlgo::Yee) {
- deltat = cfl * CartesianYeeAlgorithm::ComputeMaxDt( dx );
- } else if (maxwell_solver_id == MaxwellSolverAlgo::CKC) {
- deltat = cfl * CartesianCKCAlgorithm::ComputeMaxDt( dx );
-# endif
+ // - In Cartesian 3D geometry: determined by the minimum cell size in all direction
+ deltat = cfl * std::min(dx[0], std::min(dx[1], dx[2])) / PhysConst::c;
+#endif
} else {
- amrex::Abort("Unknown algorithm");
- }
-
+ // Computation of dt for FDTD algorithm
+#ifdef WARPX_DIM_RZ
+ // - In RZ geometry
+ if (maxwell_solver_id == MaxwellSolverAlgo::Yee) {
+ deltat = cfl * CylindricalYeeAlgorithm::ComputeMaxDt(dx, n_rz_azimuthal_modes);
+#else
+ // - In Cartesian geometry
+ if (do_nodal) {
+ deltat = cfl * CartesianNodalAlgorithm::ComputeMaxDt(dx);
+ } else if (maxwell_solver_id == MaxwellSolverAlgo::Yee) {
+ deltat = cfl * CartesianYeeAlgorithm::ComputeMaxDt(dx);
+ } else if (maxwell_solver_id == MaxwellSolverAlgo::CKC) {
+ deltat = cfl * CartesianCKCAlgorithm::ComputeMaxDt(dx);
#endif
+ } else {
+ amrex::Abort("ComputeDt: Unknown algorithm");
+ }
+ }
dt.resize(0);
dt.resize(max_level+1,deltat);
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
index bf25e4263..a037c973c 100644
--- a/Source/Evolve/WarpXEvolve.cpp
+++ b/Source/Evolve/WarpXEvolve.cpp
@@ -56,9 +56,8 @@ WarpX::Evolve (int numsteps)
amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(0);
if (cost) {
-#ifdef WARPX_USE_PSATD
- amrex::Abort("LoadBalance for PSATD: TODO");
-#endif
+ if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
+ amrex::Abort("LoadBalance for PSATD: TODO");
if (step > 0 && load_balance_intervals.contains(step+1))
{
LoadBalance();
@@ -114,10 +113,9 @@ WarpX::Evolve (int numsteps)
FillBoundaryE_avg(guard_cells.ng_FieldGather, guard_cells.ng_Extra);
FillBoundaryB_avg(guard_cells.ng_FieldGather, guard_cells.ng_Extra);
}
-#ifndef WARPX_USE_PSATD
// TODO Remove call to FillBoundaryAux before UpdateAuxilaryData?
- FillBoundaryAux(guard_cells.ng_UpdateAux);
-#endif
+ if (WarpX::maxwell_solver_id != MaxwellSolverAlgo::PSATD)
+ FillBoundaryAux(guard_cells.ng_UpdateAux);
UpdateAuxilaryData();
FillBoundaryAux(guard_cells.ng_UpdateAux);
}
@@ -288,14 +286,16 @@ WarpX::OneStep_nosub (Real cur_time)
// 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
- 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");
- if ( !fft_periodic_single_box && (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) )
- amrex::Abort("\nVay current deposition 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
+ if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+ 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");
+ if (!fft_periodic_single_box && (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay))
+ amrex::Abort(
+ "\nVay current deposition 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");
+ }
#ifdef WARPX_QED
doQEDEvents();
@@ -309,13 +309,13 @@ 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
- if ( fft_periodic_single_box && current_correction ) CurrentCorrection();
- if ( fft_periodic_single_box && (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) )
- VayDeposition();
-#endif
+ // Apply current correction in Fourier space: for periodic single-box global FFTs
+ // without guard cells, apply this after calling SyncCurrent
+ if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+ if (fft_periodic_single_box && current_correction) CurrentCorrection();
+ if (fft_periodic_single_box && (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay))
+ VayDeposition();
+ }
// At this point, J is up-to-date inside the domain, and E and B are
@@ -330,7 +330,7 @@ WarpX::OneStep_nosub (Real cur_time)
// Electromagnetic solver:
// Push E and B from {n} to {n+1}
// (And update guard cells immediately afterwards)
-#ifdef WARPX_USE_PSATD
+ if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
if (use_hybrid_QED)
{
WarpX::Hybrid_QED_Push(dt);
@@ -344,13 +344,12 @@ WarpX::OneStep_nosub (Real cur_time)
{
WarpX::Hybrid_QED_Push(dt);
FillBoundaryE(guard_cells.ng_alloc_EB, guard_cells.ng_Extra);
-
}
if (do_pml) DampPML();
-#else
- EvolveF(0.5*dt[0], DtType::FirstHalf);
+ } else {
+ EvolveF(0.5 * dt[0], DtType::FirstHalf);
FillBoundaryF(guard_cells.ng_FieldSolverF);
- EvolveB(0.5*dt[0]); // We now have B^{n+1/2}
+ EvolveB(0.5 * dt[0]); // We now have B^{n+1/2}
FillBoundaryB(guard_cells.ng_FieldSolver, IntVect::TheZeroVector());
if (WarpX::em_solver_medium == MediumForEM::Vacuum) {
@@ -364,8 +363,8 @@ WarpX::OneStep_nosub (Real cur_time)
}
FillBoundaryE(guard_cells.ng_FieldSolver, IntVect::TheZeroVector());
- EvolveF(0.5*dt[0], DtType::SecondHalf);
- EvolveB(0.5*dt[0]); // We now have B^{n+1}
+ EvolveF(0.5 * dt[0], DtType::SecondHalf);
+ EvolveB(0.5 * dt[0]); // We now have B^{n+1}
if (do_pml) {
FillBoundaryF(guard_cells.ng_alloc_F);
DampPML();
@@ -375,10 +374,10 @@ WarpX::OneStep_nosub (Real cur_time)
}
// E and B are up-to-date in the domain, but all guard cells are
// outdated.
- if ( safe_guard_cells )
+ if (safe_guard_cells)
FillBoundaryB(guard_cells.ng_alloc_EB, guard_cells.ng_Extra);
-#endif
- }
+ } // !PSATD
+ } // !do_electrostatic
}
/* /brief Perform one PIC iteration, with subcycling
@@ -735,27 +734,45 @@ WarpX::applyMirrors(Real time){
}
// Apply current correction in Fourier space
-#ifdef WARPX_USE_PSATD
void
WarpX::CurrentCorrection ()
{
- for ( int lev = 0; lev <= finest_level; ++lev )
+#ifdef WARPX_USE_PSATD
+ if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
{
- spectral_solver_fp[lev]->CurrentCorrection( current_fp[lev], rho_fp[lev] );
- if ( spectral_solver_cp[lev] ) spectral_solver_cp[lev]->CurrentCorrection( current_cp[lev], rho_cp[lev] );
+ for ( int lev = 0; lev <= finest_level; ++lev )
+ {
+ spectral_solver_fp[lev]->CurrentCorrection( current_fp[lev], rho_fp[lev] );
+ if ( spectral_solver_cp[lev] ) spectral_solver_cp[lev]->CurrentCorrection( current_cp[lev], rho_cp[lev] );
+ }
+ } else {
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false,
+ "WarpX::CurrentCorrection: only implemented for spectral solver.");
}
-}
+#else
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false,
+ "WarpX::CurrentCorrection: requires WarpX build with spectral solver support.");
#endif
+}
// Compute current from Vay deposition in Fourier space
-#ifdef WARPX_USE_PSATD
void
WarpX::VayDeposition ()
{
- for (int lev = 0; lev <= finest_level; ++lev)
+#ifdef WARPX_USE_PSATD
+ if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
{
- spectral_solver_fp[lev]->VayDeposition(current_fp[lev]);
- if (spectral_solver_cp[lev]) spectral_solver_cp[lev]->VayDeposition(current_cp[lev]);
+ for (int lev = 0; lev <= finest_level; ++lev)
+ {
+ spectral_solver_fp[lev]->VayDeposition(current_fp[lev]);
+ if (spectral_solver_cp[lev]) spectral_solver_cp[lev]->VayDeposition(current_cp[lev]);
+ }
+ } else {
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false,
+ "WarpX::VayDeposition: only implemented for spectral solver.");
}
-}
+#else
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false,
+ "WarpX::CurrentCorrection: requires WarpX build with spectral solver support.");
#endif
+}
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-11 22:22:28 +0000