Skip to content

Commit

Permalink
added time step (m_dt) as member to ImplicitSolver. Removed passing d…
Browse files Browse the repository at this point in the history 
…t to ComputeRHS().
  • Loading branch information
@JustinRayAngus
JustinRayAngus committed Nov 5, 2024
1 parent a4d5631 commit 259059b
Show file tree
Hide file tree
Showing 9 changed files with 33 additions and 30 deletions.
6 changes: 5 additions & 1 deletion Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,6 @@ public:
virtual void ComputeRHS ( WarpXSolverVec& a_RHS,
const WarpXSolverVec& a_E,
amrex::Real a_time,
amrex::Real a_dt,
int a_nl_iter,
bool a_from_jacobian ) = 0;

Expand Down Expand Up @@ -111,6 +110,11 @@ protected:
*/
int m_num_amr_levels = 1;

/**
* \brief Time step
*/
mutable amrex::Real m_dt = DBL_MAX;

/**
* \brief Nonlinear solver type and object
*/
Expand Down
1 change: 0 additions & 1 deletion Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,6 @@ public:
void ComputeRHS ( WarpXSolverVec& a_RHS,
const WarpXSolverVec& a_E,
amrex::Real a_time,
amrex::Real a_dt,
int a_nl_iter,
bool a_from_jacobian ) override;

Expand Down
14 changes: 8 additions & 6 deletions Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,9 @@ void SemiImplicitEM::OneStep ( amrex::Real a_time,
{
amrex::ignore_unused(a_step);

// Set the member time step
m_dt = a_dt;

// Fields have Eg^{n}, Bg^{n-1/2}
// Particles have up^{n} and xp^{n}.

Expand All @@ -68,15 +71,15 @@ void SemiImplicitEM::OneStep ( amrex::Real a_time,
m_Eold.Copy( FieldType::Efield_fp );

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

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

// 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^{n+1/2}
m_nlsolver->Solve( m_E, m_Eold, half_time, a_dt );
m_nlsolver->Solve( m_E, m_Eold, half_time, m_dt );

// Update WarpX owned Efield_fp to t_{n+1/2}
m_WarpX->SetElectricFieldAndApplyBCs( m_E );
Expand All @@ -94,7 +97,6 @@ void SemiImplicitEM::OneStep ( amrex::Real a_time,
void SemiImplicitEM::ComputeRHS ( WarpXSolverVec& a_RHS,
const WarpXSolverVec& a_E,
amrex::Real a_time,
amrex::Real a_dt,
int a_nl_iter,
bool a_from_jacobian )
{
Expand All @@ -104,8 +106,8 @@ void SemiImplicitEM::ComputeRHS ( WarpXSolverVec& a_RHS,

// Update particle positions and velocities using the current state
// of Eg and Bg. Deposit current density at time n+1/2
m_WarpX->ImplicitPreRHSOp( a_time, a_dt, a_nl_iter, a_from_jacobian );
m_WarpX->ImplicitPreRHSOp( a_time, m_dt, a_nl_iter, a_from_jacobian );

// 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);
m_WarpX->ImplicitComputeRHSE(0.5_rt*m_dt, a_RHS);
}
4 changes: 1 addition & 3 deletions Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,6 @@ public:
void ComputeRHS ( WarpXSolverVec& a_RHS,
const WarpXSolverVec& a_E,
amrex::Real a_time,
amrex::Real a_dt,
int a_nl_iter,
bool a_from_jacobian ) override;

Expand All @@ -101,8 +100,7 @@ private:
* \brief Update the E and B fields owned by WarpX
*/
void UpdateWarpXFields ( const WarpXSolverVec& a_E,
amrex::Real a_time,
amrex::Real a_dt );
amrex::Real a_time );

/**
* \brief Nonlinear solver is for the time-centered values of E. After
Expand Down
25 changes: 13 additions & 12 deletions Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,9 @@ void ThetaImplicitEM::OneStep ( const amrex::Real a_time,
// Fields have Eg^{n} and Bg^{n}
// Particles have up^{n} and xp^{n}.

// Set the member time step
m_dt = a_dt;

// Save up and xp at the start of the time step
m_WarpX->SaveParticlesAtImplicitStepStart ( );

Expand All @@ -99,47 +102,45 @@ void ThetaImplicitEM::OneStep ( const amrex::Real a_time,
}
}

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

// 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^{n+theta}
m_nlsolver->Solve( m_E, m_Eold, theta_time, a_dt );
m_nlsolver->Solve( m_E, m_Eold, theta_time, m_dt );

// Update WarpX owned Efield_fp and Bfield_fp to t_{n+theta}
UpdateWarpXFields( m_E, theta_time, a_dt );
UpdateWarpXFields( m_E, theta_time );

// Advance particles from time n+1/2 to time n+1
m_WarpX->FinishImplicitParticleUpdate();

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

}

void ThetaImplicitEM::ComputeRHS ( WarpXSolverVec& a_RHS,
const WarpXSolverVec& a_E,
amrex::Real a_time,
amrex::Real a_dt,
int a_nl_iter,
bool a_from_jacobian )
{
// Update WarpX-owned Efield_fp and Bfield_fp using current state of
// Eg from the nonlinear solver at time n+theta
UpdateWarpXFields( a_E, a_time, a_dt );
UpdateWarpXFields( a_E, a_time );

// Update particle positions and velocities using the current state
// of Eg and Bg. Deposit current density at time n+1/2
m_WarpX->ImplicitPreRHSOp( a_time, a_dt, a_nl_iter, a_from_jacobian );
m_WarpX->ImplicitPreRHSOp( a_time, m_dt, a_nl_iter, a_from_jacobian );

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

void ThetaImplicitEM::UpdateWarpXFields ( const WarpXSolverVec& a_E,
amrex::Real a_time,
amrex::Real a_dt )
amrex::Real a_time )
{
amrex::ignore_unused(a_time);

Expand All @@ -148,7 +149,7 @@ void ThetaImplicitEM::UpdateWarpXFields ( const WarpXSolverVec& a_E,

// Update Bfield_fp owned by WarpX
ablastr::fields::MultiLevelVectorField const& B_old = m_WarpX->m_fields.get_mr_levels_alldirs(FieldType::B_old, 0);
m_WarpX->UpdateMagneticFieldAndApplyBCs(B_old, m_theta * a_dt );
m_WarpX->UpdateMagneticFieldAndApplyBCs( B_old, m_theta*m_dt );

}

Expand All @@ -164,6 +165,6 @@ void ThetaImplicitEM::FinishFieldUpdate ( amrex::Real a_new_time )
m_E.linComb( c0, m_E, c1, m_Eold );
m_WarpX->SetElectricFieldAndApplyBCs( m_E );
ablastr::fields::MultiLevelVectorField const & B_old = m_WarpX->m_fields.get_mr_levels_alldirs(FieldType::B_old, 0);
m_WarpX->FinishMagneticFieldAndApplyBCs(B_old, m_theta );
m_WarpX->FinishMagneticFieldAndApplyBCs( B_old, m_theta );

}
2 changes: 1 addition & 1 deletion Source/NonlinearSolvers/JacobianFunctionMF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -253,7 +253,7 @@ void JacobianFunctionMF<T,Ops>::apply (T& a_dF, const T& a_dU)
const RT eps_inv = 1.0_rt/eps;

m_Z.linComb( 1.0, m_Y0, eps, a_dU ); // Z = Y0 + eps*dU
m_ops->ComputeRHS(m_R, m_Z, m_cur_time, m_dt, -1, true );
m_ops->ComputeRHS(m_R, m_Z, m_cur_time, -1, true );

// F(Y) = Y - b - R(Y) ==> dF = dF/dY*dU = [1 - dR/dY]*dU
// = dU - (R(Z)-R(Y0))/eps
Expand Down
6 changes: 2 additions & 4 deletions Source/NonlinearSolvers/NewtonSolver.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,6 @@ private:
const Vec& a_U,
const Vec& a_b,
amrex::Real a_time,
amrex::Real a_dt,
int a_iter ) const;

};
Expand Down Expand Up @@ -252,7 +251,7 @@ void NewtonSolver<Vec,Ops>::Solve ( Vec& a_U,
for (iter = 0; iter < m_maxits;) {

// Compute residual: F(U) = U - b - R(U)
EvalResidual(m_F, a_U, a_b, a_time, a_dt, iter);
EvalResidual(m_F, a_U, a_b, a_time, iter);

// Compute norm of the residual
norm_abs = m_F.norm2();
Expand Down Expand Up @@ -329,11 +328,10 @@ void NewtonSolver<Vec,Ops>::EvalResidual ( Vec& a_F,
const Vec& a_U,
const Vec& a_b,
amrex::Real a_time,
amrex::Real a_dt,
int a_iter ) const
{

m_ops->ComputeRHS( m_R, a_U, a_time, a_dt, a_iter, false );
m_ops->ComputeRHS( m_R, a_U, a_time, a_iter, false );

// set base U and R(U) for matrix-free Jacobian action calculation
m_linear_function->setBaseSolution(a_U);
Expand Down
2 changes: 1 addition & 1 deletion Source/NonlinearSolvers/NonlinearSolver.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
* This class is templated on a vector class Vec, and an operator class Ops.
*
* The Ops class must have the following function:
* ComputeRHS( R_vec, U_vec, time, dt, nl_iter, from_jacobian ),
* ComputeRHS( R_vec, U_vec, time, nl_iter, from_jacobian ),
* where U_vec and R_vec are of type Vec.
*
* The Vec class must have basic math operators, such as Copy, +=, -=,
Expand Down
3 changes: 2 additions & 1 deletion Source/NonlinearSolvers/PicardSolver.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,7 @@ void PicardSolver<Vec,Ops>::Solve ( Vec& a_U,
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
this->m_is_defined,
"PicardSolver::Solve() called on undefined object");
amrex::ignore_unused(a_dt);
using namespace amrex::literals;

//
Expand All @@ -156,7 +157,7 @@ void PicardSolver<Vec,Ops>::Solve ( Vec& a_U,
m_Usave.Copy(a_U);

// Update the solver state (a_U = a_b + m_R)
m_ops->ComputeRHS( m_R, a_U, a_time, a_dt, iter, false );
m_ops->ComputeRHS( m_R, a_U, a_time, iter, false );
a_U.Copy(a_b);
a_U += m_R;

Expand Down

0 comments on commit 259059b

Please sign in to comment.