Skip to content

Commit

Permalink
refactoring comments and how to Copy WarpXSolverVec.
Browse files Browse the repository at this point in the history
  • Loading branch information
@JustinRayAngus
JustinRayAngus committed Aug 24, 2024
1 parent 31da2a8 commit 2a55d34
Show file tree
Hide file tree
Showing 4 changed files with 56 additions and 44 deletions.
18 changes: 9 additions & 9 deletions Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,24 +59,24 @@ void SemiImplicitEM::OneStep ( amrex::Real a_time,
{
amrex::ignore_unused(a_step);

// Fields have E^{n}, B^{n-1/2}
// Fields have Eg^{n}, Bg^{n-1/2}
// Particles have p^{n} and x^{n}.

// Save the values at the start of the time step,
m_WarpX->SaveParticlesAtImplicitStepStart ( );

// Save the fields at the start of the step
m_Eold.Copy( m_WarpX->getMultiLevelField(FieldType::Efield_fp), FieldType::Efield_fp );
m_E.Copy(m_Eold); // initial guess for E
// Save electric field at the start of the step
m_Eold.Copy( FieldType::Efield_fp );

// Compute Bfield at time n+1/2
// Advance WarpX owned Bfield_fp to t_{n+1/2}
m_WarpX->EvolveB(a_dt, DtType::Full);
m_WarpX->ApplyMagneticFieldBCs();

const amrex::Real half_time = a_time + 0.5_rt*a_dt;

// Solve nonlinear system for E at t_{n+1/2}
// Solve nonlinear system for Eg at t_{n+1/2}
// Particles will be advanced to t_{n+1/2}
m_E.Copy( m_Eold ); // initial guess for Eg
m_nlsolver->Solve( m_E, m_Eold, half_time, a_dt );

// Update WarpX owned Efield_fp to t_{n+1/2}
Expand All @@ -85,8 +85,8 @@ void SemiImplicitEM::OneStep ( amrex::Real a_time,
// Advance particles from time n+1/2 to time n+1
m_WarpX->FinishImplicitParticleUpdate();

// Advance E fields from time n+1/2 to time n+1
// Eg^{n+1} = 2.0*E_g^{n+1/2} - E_g^n
// Advance Eg from time n+1/2 to time n+1
// Eg^{n+1} = 2.0*Eg^{n+1/2} - Eg^n
m_E.linComb( 2._rt, m_E, -1._rt, m_Eold );
m_WarpX->SetElectricFieldAndApplyBCs( m_E );

Expand All @@ -107,6 +107,6 @@ void SemiImplicitEM::ComputeRHS ( WarpXSolverVec& a_RHS,
// current state of the fields E and B. Deposit current density at time n+1/2.
m_WarpX->ImplicitPreRHSOp( a_time, a_dt, a_nl_iter, a_from_jacobian );

// RHS = cvac^2*0.5*dt*( curl(B^{n+1/2}) - mu0*J^{n+1/2} )
// RHS = cvac^2*0.5*dt*( curl(Bg^{n+1/2}) - mu0*Jg^{n+1/2} )
m_WarpX->ImplicitComputeRHSE(0.5_rt*a_dt, a_RHS);
}
18 changes: 9 additions & 9 deletions Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,15 +80,14 @@ void ThetaImplicitEM::OneStep ( const amrex::Real a_time,
{
amrex::ignore_unused(a_step);

// Fields have E^{n} and B^{n}
// Fields have Eg^{n} and Bg^{n}
// Particles have p^{n} and x^{n}.

// Save the values at the start of the time step,
m_WarpX->SaveParticlesAtImplicitStepStart ( );

// Save the fields at the start of the step
m_Eold.Copy( m_WarpX->getMultiLevelField(FieldType::Efield_fp), FieldType::Efield_fp );
m_E.Copy(m_Eold); // initial guess for E
// Save electric field at the start of the step
m_Eold.Copy( FieldType::Efield_fp );

const int num_levels = static_cast<int>(m_Bold.size());
for (int lev = 0; lev < num_levels; ++lev) {
Expand All @@ -101,8 +100,9 @@ void ThetaImplicitEM::OneStep ( const amrex::Real a_time,

const amrex::Real theta_time = a_time + m_theta*a_dt;

// Solve nonlinear system for E at t_{n+theta}
// Solve nonlinear system for Eg at t_{n+theta}
// Particles will be advanced to t_{n+1/2}
m_E.Copy( m_Eold ); // initial guess for Eg
m_nlsolver->Solve( m_E, m_Eold, theta_time, a_dt );

// Update WarpX owned Efield_fp and Bfield_fp to t_{n+theta}
Expand All @@ -111,7 +111,7 @@ void ThetaImplicitEM::OneStep ( const amrex::Real a_time,
// Advance particles from time n+1/2 to time n+1
m_WarpX->FinishImplicitParticleUpdate();

// Advance E and B fields from time n+theta to time n+1
// Advance Eg and Bg from time n+theta to time n+1
const amrex::Real new_time = a_time + a_dt;
FinishFieldUpdate( new_time );

Expand All @@ -132,7 +132,7 @@ void ThetaImplicitEM::ComputeRHS ( WarpXSolverVec& a_RHS,
// current state of the fields E and B. Deposit current density at time n+1/2.
m_WarpX->ImplicitPreRHSOp( a_time, a_dt, a_nl_iter, a_from_jacobian );

// RHS = cvac^2*m_theta*dt*( curl(B^{n+theta}) - mu0*J^{n+1/2} )
// RHS = cvac^2*m_theta*dt*( curl(Bg^{n+theta}) - mu0*Jg^{n+1/2} )
m_WarpX->ImplicitComputeRHSE(m_theta*a_dt, a_RHS);
}

Expand All @@ -154,8 +154,8 @@ void ThetaImplicitEM::FinishFieldUpdate ( amrex::Real a_new_time )
{
amrex::ignore_unused(a_new_time);

// Eg^{n+1} = (1/theta)*E_g^{n+theta} + (1-1/theta)*E_g^n
// Bg^{n+1} = (1/theta)*B_g^{n+theta} + (1-1/theta)*B_g^n
// Eg^{n+1} = (1/theta)*Eg^{n+theta} + (1-1/theta)*Eg^n
// Bg^{n+1} = (1/theta)*Bg^{n+theta} + (1-1/theta)*Bg^n

const amrex::Real c0 = 1._rt/m_theta;
const amrex::Real c1 = 1._rt - c0;
Expand Down
29 changes: 9 additions & 20 deletions Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,24 +75,7 @@ public:

[[nodiscard]] RT dotProduct( const WarpXSolverVec& a_X ) const;

inline
void Copy ( const amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > >& a_field_vec,
const warpx::fields::FieldType a_field_type )
{

WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
IsDefined(),
"WarpXSolverVec::Copy() called on undefined WarpXSolverVec");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
a_field_type==m_field_type,
"WarpXSolverVec::Copy() called with vecs of different types");
for (int lev=0; lev<m_num_amr_levels; ++lev) {
for (int n=0; n<3; ++n) {
amrex::MultiFab::Copy( *m_field_vec[lev][n], *a_field_vec[lev][n], 0, 0, m_ncomp,
amrex::IntVect::TheZeroVector() );
}
}
}
void Copy ( const warpx::fields::FieldType a_field_type );

inline
void Copy ( const WarpXSolverVec& a_solver_vec )
Expand All @@ -106,8 +89,14 @@ public:
"WarpXSolverVec::Copy(solver_vec) called with vecs of different types");
}
else { Define(a_solver_vec); }
const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& solver_vec = a_solver_vec.getVec();
Copy(solver_vec, a_solver_vec.getVecType());

const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& this_vec = a_solver_vec.getVec();
for (int lev=0; lev<m_num_amr_levels; ++lev) {
for (int n=0; n<3; ++n) {
amrex::MultiFab::Copy( *m_field_vec[lev][n], *this_vec[lev][n], 0, 0, m_ncomp,
amrex::IntVect::TheZeroVector() );
}
}
}

// Prohibit Copy assignment operator
Expand Down
35 changes: 29 additions & 6 deletions Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,20 +14,43 @@ void WarpXSolverVec::Define ( WarpX* a_WarpX,
!IsDefined(),
"WarpXSolverVec::Define(a_vec, a_type) called on already defined WarpXSolverVec");

const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& this_vec(a_WarpX->getMultiLevelField(a_field_type));
const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& this_vec(
a_WarpX->getMultiLevelField(a_field_type));
m_field_vec.resize(m_num_amr_levels);
const int lev = 0;
for (int n=0; n<3; n++) {
const amrex::MultiFab& mf_model = *this_vec[lev][n];
m_field_vec[lev][n] = std::make_unique<amrex::MultiFab>( mf_model.boxArray(), mf_model.DistributionMap(),
mf_model.nComp(), amrex::IntVect::TheZeroVector() );
for (int lev=0; lev<m_num_amr_levels; ++lev) {
for (int n=0; n<3; n++) {
const amrex::MultiFab& mf_model = *this_vec[lev][n];
m_field_vec[lev][n] = std::make_unique<amrex::MultiFab>( mf_model.boxArray(),
mf_model.DistributionMap(),
mf_model.nComp(),
amrex::IntVect::TheZeroVector() );
}
}
m_field_type = a_field_type;
m_is_defined = true;
SetWarpXPointer(a_WarpX);
SetDotMask();
}

void WarpXSolverVec::Copy ( const warpx::fields::FieldType a_field_type )
{
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
IsDefined(),
"WarpXSolverVec::Copy() called on undefined WarpXSolverVec");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
a_field_type==m_field_type,
"WarpXSolverVec::Copy() called with vecs of different types");

const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& this_vec(
m_WarpX->getMultiLevelField(a_field_type));
for (int lev=0; lev<m_num_amr_levels; ++lev) {
for (int n=0; n<3; ++n) {
amrex::MultiFab::Copy( *m_field_vec[lev][n], *this_vec[lev][n], 0, 0, m_ncomp,
amrex::IntVect::TheZeroVector() );
}
}
}

void WarpXSolverVec::SetWarpXPointer( WarpX* a_WarpX )
{
if (m_warpx_ptr_defined) { return; }
Expand Down

0 comments on commit 2a55d34

Please sign in to comment.