Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Remove theta() function from ImplicitSolver base class. #5441

Merged
merged 3 commits into from
Nov 15, 2024
Merged

Remove theta() function from ImplicitSolver base class. #5441

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
570eeba
added time step (m_dt) as member to ImplicitSolver. Removed passing d…
JustinRayAngus Nov 4, 2024
8375975
dt passed to m_nlsolver->Solve() and eventually used by the PC is now…
JustinRayAngus Nov 4, 2024
72c56c2
incorporating changes into new strang implicit spectral time advance …
JustinRayAngus Nov 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 5 additions & 3 deletions Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,12 +83,9 @@ public:
virtual void ComputeRHS ( WarpXSolverVec& a_RHS,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I notice that the comments do not have the list of input parameters for the routines. This is not for this PR, but can be added in the future.

const WarpXSolverVec& a_E,
amrex::Real a_time,
amrex::Real a_dt,
int a_nl_iter,
bool a_from_jacobian ) = 0;

[[nodiscard]] virtual amrex::Real theta () const { return 1.0; }

[[nodiscard]] int numAMRLevels () const { return m_num_amr_levels; }

[[nodiscard]] const amrex::Geometry& GetGeometry (int) const;
Expand All @@ -111,6 +108,11 @@ protected:
*/
int m_num_amr_levels = 1;

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

/**
* \brief Nonlinear solver type and object
*/
Expand Down
1 change: 0 additions & 1 deletion Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H
View file
Open in desktop
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
Open in desktop
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, 0.5_rt*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/StrangImplicitSpectralEM.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,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 Expand Up @@ -93,8 +92,7 @@ private:
* \brief Update the E and B fields owned by WarpX
*/
void UpdateWarpXFields ( WarpXSolverVec const& 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
21 changes: 11 additions & 10 deletions Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,9 @@ void StrangImplicitSpectralEM::OneStep ( amrex::Real a_time,
// Fields have E^{n} and B^{n}
// Particles have p^{n} and x^{n}.

// Set the member time step
m_dt = a_dt;

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

Expand All @@ -73,20 +76,20 @@ void StrangImplicitSpectralEM::OneStep ( amrex::Real a_time,
m_Eold.Copy( FieldType::Efield_fp );
m_E.Copy(m_Eold); // initial guess for E

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

// Solve nonlinear system for E at t_{n+1/2}
// Particles will be advanced to t_{n+1/2}
m_nlsolver->Solve( m_E, m_Eold, half_time, a_dt );
m_nlsolver->Solve( m_E, m_Eold, half_time, 0.5_rt*m_dt );

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

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

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

// Advance the fields to time n+1 source free
Expand All @@ -97,29 +100,27 @@ void StrangImplicitSpectralEM::OneStep ( amrex::Real a_time,
void StrangImplicitSpectralEM::ComputeRHS ( WarpXSolverVec& a_RHS,
WarpXSolverVec const & 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
// E from the nonlinear solver at time n+1/2
UpdateWarpXFields( a_E, a_time, a_dt );
UpdateWarpXFields( a_E, a_time );

// Self consistently update particle positions and velocities using the
// 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 );
m_WarpX->ImplicitPreRHSOp( a_time, m_dt, a_nl_iter, a_from_jacobian );

// For Strang split implicit PSATD, the RHS = -dt*mu*c**2*J
bool const allow_type_mismatch = true;
a_RHS.Copy(FieldType::current_fp, warpx::fields::FieldType::None, allow_type_mismatch);
amrex::Real constexpr coeff = PhysConst::c * PhysConst::c * PhysConst::mu0;
a_RHS.scale(-coeff * 0.5_rt*a_dt);
a_RHS.scale(-coeff * 0.5_rt*m_dt);

}

void StrangImplicitSpectralEM::UpdateWarpXFields (WarpXSolverVec const & a_E,
amrex::Real /*a_time*/,
amrex::Real /*a_dt*/)
amrex::Real /*a_time*/ )
{

// Update Efield_fp owned by WarpX
Expand Down
6 changes: 1 addition & 5 deletions Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,12 +74,9 @@ 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;

[[nodiscard]] amrex::Real theta () const override { return m_theta; }

private:

/**
Expand All @@ -101,8 +98,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
Open in desktop
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_theta*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 );

}
5 changes: 3 additions & 2 deletions Source/NonlinearSolvers/CurlCurlMLMGPC.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -272,7 +272,8 @@ void CurlCurlMLMGPC<T,Ops>::Update (const T& a_U)
amrex::ignore_unused(a_U);

// set the coefficients alpha and beta for curl-curl op
const RT alpha = (m_ops->theta()*this->m_dt*PhysConst::c) * (m_ops->theta()*this->m_dt*PhysConst::c);
// (m_dt here is actually theta<=0.5 times simulation dt)
const RT alpha = (this->m_dt*PhysConst::c) * (this->m_dt*PhysConst::c);
const RT beta = RT(1.0);

// currently not implemented in 1D
Expand All @@ -282,7 +283,7 @@ void CurlCurlMLMGPC<T,Ops>::Update (const T& a_U)

if (m_verbose) {
amrex::Print() << "Updating " << amrex::getEnumNameString(PreconditionerType::pc_curl_curl_mlmg)
<< ": dt = " << this->m_dt << ", "
<< ": theta*dt = " << this->m_dt << ", "
<< " coefficients: "
<< "alpha = " << alpha << ", "
<< "beta = " << beta << "\n";
Expand Down
2 changes: 1 addition & 1 deletion Source/NonlinearSolvers/JacobianFunctionMF.H
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading