diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 0a9f21e6863..90f8c87df48 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -159,6 +159,7 @@ class FiniteDifferenceSolver std::unique_ptr const& rhofield, std::unique_ptr const& Pefield, std::array< std::unique_ptr, 3 > const& edge_lengths, + amrex::Real t, int lev, HybridPICModel const* hybrid_model, bool solve_for_Faraday ); @@ -244,6 +245,7 @@ class FiniteDifferenceSolver std::unique_ptr const& rhofield, std::unique_ptr const& Pefield, std::array< std::unique_ptr, 3 > const& edge_lengths, + amrex::Real t, int lev, HybridPICModel const* hybrid_model, bool solve_for_Faraday ); @@ -349,6 +351,7 @@ class FiniteDifferenceSolver std::unique_ptr const& rhofield, std::unique_ptr const& Pefield, std::array< std::unique_ptr, 3 > const& edge_lengths, + amrex::Real t, int lev, HybridPICModel const* hybrid_model, bool solve_for_Faraday ); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H index e17bb554275..ebd71be7e84 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H @@ -91,6 +91,7 @@ public: amrex::Vector, 3>> const& Bfield, amrex::Vector> const& rhofield, amrex::Vector, 3>> const& edge_lengths, + amrex::Real t, bool solve_for_Faraday); void HybridPICSolveE ( @@ -99,6 +100,7 @@ public: std::array< std::unique_ptr, 3> const& Bfield, std::unique_ptr const& rhofield, std::array< std::unique_ptr, 3> const& edge_lengths, + amrex::Real t, int lev, bool solve_for_Faraday); void HybridPICSolveE ( @@ -107,6 +109,7 @@ public: std::array< std::unique_ptr, 3> const& Bfield, std::unique_ptr const& rhofield, std::array< std::unique_ptr, 3> const& edge_lengths, + amrex::Real t, int lev, PatchType patch_type, bool solve_for_Faraday); void BfieldEvolveRK ( @@ -115,7 +118,7 @@ public: amrex::Vector, 3>> const& Jfield, amrex::Vector> const& rhofield, amrex::Vector, 3>> const& edge_lengths, - amrex::Real dt, DtType a_dt_type, + amrex::Real t, amrex::Real dt, DtType a_dt_type, amrex::IntVect ng, std::optional nodal_sync); void BfieldEvolveRK ( @@ -124,7 +127,7 @@ public: amrex::Vector, 3>> const& Jfield, amrex::Vector> const& rhofield, amrex::Vector, 3>> const& edge_lengths, - amrex::Real dt, int lev, DtType dt_type, + amrex::Real t, amrex::Real dt, int lev, DtType dt_type, amrex::IntVect ng, std::optional nodal_sync); void FieldPush ( @@ -133,7 +136,7 @@ public: amrex::Vector, 3>> const& Jfield, amrex::Vector> const& rhofield, amrex::Vector, 3>> const& edge_lengths, - amrex::Real dt, DtType dt_type, + amrex::Real t, amrex::Real dt, DtType dt_type, amrex::IntVect ng, std::optional nodal_sync); /** diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp index 621765d5a49..236dd708725 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp @@ -465,6 +465,7 @@ void HybridPICModel::HybridPICSolveE ( amrex::Vector, 3>> const& Bfield, amrex::Vector> const& rhofield, amrex::Vector, 3>> const& edge_lengths, + amrex::Real t, const bool solve_for_Faraday) { auto& warpx = WarpX::GetInstance(); @@ -472,7 +473,7 @@ void HybridPICModel::HybridPICSolveE ( { HybridPICSolveE( Efield[lev], Jfield[lev], Bfield[lev], rhofield[lev], - edge_lengths[lev], lev, solve_for_Faraday + edge_lengths[lev], t, lev, solve_for_Faraday ); } } @@ -483,12 +484,13 @@ void HybridPICModel::HybridPICSolveE ( std::array< std::unique_ptr, 3> const& Bfield, std::unique_ptr const& rhofield, std::array< std::unique_ptr, 3> const& edge_lengths, + amrex::Real t, const int lev, const bool solve_for_Faraday) { WARPX_PROFILE("WarpX::HybridPICSolveE()"); HybridPICSolveE( - Efield, Jfield, Bfield, rhofield, edge_lengths, lev, + Efield, Jfield, Bfield, rhofield, edge_lengths, t, lev, PatchType::fine, solve_for_Faraday ); if (lev > 0) @@ -504,6 +506,7 @@ void HybridPICModel::HybridPICSolveE ( std::array< std::unique_ptr, 3> const& Bfield, std::unique_ptr const& rhofield, std::array< std::unique_ptr, 3> const& edge_lengths, + amrex::Real t, const int lev, PatchType patch_type, const bool solve_for_Faraday) { @@ -514,7 +517,7 @@ void HybridPICModel::HybridPICSolveE ( Efield, current_fp_ampere[lev], Jfield, current_fp_external[lev], Bfield, rhofield, electron_pressure_fp[lev], - edge_lengths, lev, this, solve_for_Faraday + edge_lengths, t, lev, this, solve_for_Faraday ); warpx.ApplyEfieldBoundary(lev, patch_type); } @@ -576,14 +579,14 @@ void HybridPICModel::BfieldEvolveRK ( amrex::Vector, 3>> const& Jfield, amrex::Vector> const& rhofield, amrex::Vector, 3>> const& edge_lengths, - amrex::Real dt, DtType dt_type, + amrex::Real t, amrex::Real dt, DtType dt_type, IntVect ng, std::optional nodal_sync ) { auto& warpx = WarpX::GetInstance(); for (int lev = 0; lev <= warpx.finestLevel(); ++lev) { BfieldEvolveRK( - Bfield, Efield, Jfield, rhofield, edge_lengths, dt, lev, dt_type, + Bfield, Efield, Jfield, rhofield, edge_lengths, t, dt, lev, dt_type, ng, nodal_sync ); } @@ -595,7 +598,7 @@ void HybridPICModel::BfieldEvolveRK ( amrex::Vector, 3>> const& Jfield, amrex::Vector> const& rhofield, amrex::Vector, 3>> const& edge_lengths, - amrex::Real dt, int lev, DtType dt_type, + amrex::Real t, amrex::Real dt, int lev, DtType dt_type, IntVect ng, std::optional nodal_sync ) { // Make copies of the B-field multifabs at t = n and create multifabs for @@ -618,11 +621,13 @@ void HybridPICModel::BfieldEvolveRK ( K[ii].setVal(0.0); } + amrex::Real t_eval = t; + // The Runge-Kutta scheme begins here. // Step 1: FieldPush( Bfield, Efield, Jfield, rhofield, edge_lengths, - 0.5_rt*dt, dt_type, ng, nodal_sync + t_eval, 0.5_rt*dt, dt_type, ng, nodal_sync ); // The Bfield is now given by: @@ -636,9 +641,10 @@ void HybridPICModel::BfieldEvolveRK ( } // Step 2: + t_eval = t+0.5_rt*dt; FieldPush( Bfield, Efield, Jfield, rhofield, edge_lengths, - 0.5_rt*dt, dt_type, ng, nodal_sync + t_eval, 0.5_rt*dt, dt_type, ng, nodal_sync ); // The Bfield is now given by: @@ -658,7 +664,7 @@ void HybridPICModel::BfieldEvolveRK ( // Step 3: FieldPush( Bfield, Efield, Jfield, rhofield, edge_lengths, - dt, dt_type, ng, nodal_sync + t_eval, dt, dt_type, ng, nodal_sync ); // The Bfield is now given by: @@ -672,9 +678,10 @@ void HybridPICModel::BfieldEvolveRK ( } // Step 4: + t_eval = t + dt; FieldPush( Bfield, Efield, Jfield, rhofield, edge_lengths, - 0.5_rt*dt, dt_type, ng, nodal_sync + t_eval, 0.5_rt*dt, dt_type, ng, nodal_sync ); // The Bfield is now given by: @@ -708,7 +715,7 @@ void HybridPICModel::FieldPush ( amrex::Vector, 3>> const& Jfield, amrex::Vector> const& rhofield, amrex::Vector, 3>> const& edge_lengths, - amrex::Real dt, DtType dt_type, + amrex::Real t, amrex::Real dt, DtType dt_type, IntVect ng, std::optional nodal_sync ) { auto& warpx = WarpX::GetInstance(); @@ -716,7 +723,7 @@ void HybridPICModel::FieldPush ( // Calculate J = curl x B / mu0 CalculateCurrentAmpere(Bfield, edge_lengths); // Calculate the E-field from Ohm's law - HybridPICSolveE(Efield, Jfield, Bfield, rhofield, edge_lengths, true); + HybridPICSolveE(Efield, Jfield, Bfield, rhofield, edge_lengths, t, true); warpx.FillBoundaryE(ng, nodal_sync); // Push forward the B-field using Faraday's law warpx.EvolveB(dt, dt_type); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp index 2b36ea4c5f3..ed47fee8433 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp @@ -360,6 +360,7 @@ void FiniteDifferenceSolver::HybridPICSolveE ( std::unique_ptr const& rhofield, std::unique_ptr const& Pefield, std::array< std::unique_ptr, 3 > const& edge_lengths, + amrex::Real t, int lev, HybridPICModel const* hybrid_model, const bool solve_for_Faraday) { @@ -370,14 +371,14 @@ void FiniteDifferenceSolver::HybridPICSolveE ( HybridPICSolveECylindrical ( Efield, Jfield, Jifield, Jextfield, Bfield, rhofield, Pefield, - edge_lengths, lev, hybrid_model, solve_for_Faraday + edge_lengths, t, lev, hybrid_model, solve_for_Faraday ); #else HybridPICSolveECartesian ( Efield, Jfield, Jifield, Jextfield, Bfield, rhofield, Pefield, - edge_lengths, lev, hybrid_model, solve_for_Faraday + edge_lengths, t, lev, hybrid_model, solve_for_Faraday ); #endif @@ -398,6 +399,7 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( std::unique_ptr const& rhofield, std::unique_ptr const& Pefield, std::array< std::unique_ptr, 3 > const& edge_lengths, + amrex::Real t, int lev, HybridPICModel const* hybrid_model, const bool solve_for_Faraday ) { @@ -431,7 +433,6 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( const auto Ez_part = hybrid_model->m_E_external[2]; auto & warpx = WarpX::GetInstance(); - auto t = warpx.gett_new(lev); auto dx_lev = warpx.Geom(lev).CellSizeArray(); const RealBox& real_box = warpx.Geom(lev).ProbDomain(); @@ -776,6 +777,7 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( std::unique_ptr const& rhofield, std::unique_ptr const& Pefield, std::array< std::unique_ptr, 3 > const& edge_lengths, + amrex::Real t, int lev, HybridPICModel const* hybrid_model, const bool solve_for_Faraday ) { @@ -803,7 +805,6 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( const auto Ez_part = hybrid_model->m_E_external[2]; auto & warpx = WarpX::GetInstance(); - auto t = warpx.gett_new(lev); auto dx_lev = warpx.Geom(lev).CellSizeArray(); const RealBox& real_box = warpx.Geom(lev).ProbDomain(); diff --git a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp index c16f0193b8d..2f711f31456 100644 --- a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp +++ b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp @@ -88,6 +88,12 @@ void WarpX::HybridPICEvolveFields () } } + // Calculate the electron pressure at t=n using rho^n + m_hybrid_pic_model->CalculateElectronPressure(DtType::FirstHalf); + + amrex::Real t_start = gett_old(0); + amrex::Real sub_dt = 0.5_rt/sub_steps*dt[0]; + // Push the B field from t=n to t=n+1/2 using the current and density // at t=n, while updating the E field along with B using the electron // momentum equation @@ -95,7 +101,9 @@ void WarpX::HybridPICEvolveFields () { m_hybrid_pic_model->BfieldEvolveRK( Bfield_fp, Efield_fp, current_fp_temp, rho_fp_temp, - m_edge_lengths, 0.5_rt/sub_steps*dt[0], + m_edge_lengths, + t_start + static_cast(sub_step)*sub_dt, + sub_dt, DtType::FirstHalf, guard_cells.ng_FieldSolver, WarpX::sync_nodal_points ); @@ -113,12 +121,19 @@ void WarpX::HybridPICEvolveFields () ); } + // Calculate the electron pressure at t=n+1/2 + m_hybrid_pic_model->CalculateElectronPressure(DtType::SecondHalf); + + t_start += 0.5_rt*dt[0]; + // Now push the B field from t=n+1/2 to t=n+1 using the n+1/2 quantities for (int sub_step = 0; sub_step < sub_steps; sub_step++) { m_hybrid_pic_model->BfieldEvolveRK( Bfield_fp, Efield_fp, current_fp, rho_fp_temp, - m_edge_lengths, 0.5_rt/sub_steps*dt[0], + m_edge_lengths, + t_start + static_cast(sub_step)*sub_dt, + sub_dt, DtType::SecondHalf, guard_cells.ng_FieldSolver, WarpX::sync_nodal_points ); @@ -148,7 +163,7 @@ void WarpX::HybridPICEvolveFields () // Update the E field to t=n+1 using the extrapolated J_i^n+1 value m_hybrid_pic_model->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths); m_hybrid_pic_model->HybridPICSolveE( - Efield_fp, current_fp_temp, Bfield_fp, rho_fp, m_edge_lengths, false + Efield_fp, current_fp_temp, Bfield_fp, rho_fp, m_edge_lengths, gett_new(0), false ); FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);