Skip to content

Commit

Permalink
[pre-commit.ci] auto fixes from pre-commit.com hooks
Browse files Browse the repository at this point in the history
for more information, see https://pre-commit.ci
  • Loading branch information
pre-commit-ci[bot] committed Feb 29, 2024
1 parent a0bbf12 commit 1826926
Show file tree
Hide file tree
Showing 16 changed files with 208 additions and 208 deletions.
2 changes: 1 addition & 1 deletion Docs/source/usage/parameters.rst
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@ Overall simulation parameters

* ``explicit``: Use an explicit solver, such as the standard FDTD or PSATD

* ``theta_implicit_em``: Use an fully implicit solver with a time-biasing parameter theta bound between 0.5 and 1.0. Exact energy conservation is achieved using theta = 0.5. Maximal damping of high-k modes is obtained using theta = 1.0. Choices for the nonlinear solver include a Picard iteration scheme and particle-suppressed (PS) JNFK.
* ``theta_implicit_em``: Use an fully implicit solver with a time-biasing parameter theta bound between 0.5 and 1.0. Exact energy conservation is achieved using theta = 0.5. Maximal damping of high-k modes is obtained using theta = 1.0. Choices for the nonlinear solver include a Picard iteration scheme and particle-suppressed (PS) JNFK.
The version implemented is an updated version that is relativistically correct, including the relativistic gamma factor for the particles.
The algorithm itself is numerical stable for large time steps. That is, it does not require time steps that resolve the plasma period or the CFL condition for light waves. However, the practicality of using a large time step depends on the nonlinear solver. Note that the Picard solver is for demonstration only. It is inefficient and will most like not converge when
:math:`\omega_{pe} \Delta t` is close to or greater than one or when the CFL condition for light waves is violated. The PS-JFNK method must be used in order to use large time steps. However, the current implementation of PS-JFNK is still inefficient because the JFNK solver is not preconditioned and there is no use of the mass matrices to minimize the cost of a linear iteration. The time step is limited by how many cells a particle can cross in a time step (MPI-related) and by the need to resolve the relavent physics.
Expand Down
8 changes: 4 additions & 4 deletions Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H
Original file line number Diff line number Diff line change
Expand Up @@ -26,16 +26,16 @@ public:
//

virtual void Define ( WarpX* const a_WarpX ) = 0;

virtual bool IsDefined () const = 0;

virtual void PrintParameters () const = 0;

virtual void Initialize () = 0;

virtual void GetParticleSolverParams (int& a_max_particle_iter,
amrex::ParticleReal& a_particle_tol ) = 0;

virtual void OneStep ( const amrex::Real a_time,
const amrex::Real a_dt,
const int a_step ) = 0;
Expand Down
14 changes: 7 additions & 7 deletions Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ class SemiImplicitEM : public ImplicitSolver
{
public:

SemiImplicitEM()
SemiImplicitEM()
{
m_nlsolver_type = NonlinearSolverType::Picard;
m_max_particle_iterations = 21;
Expand All @@ -28,15 +28,15 @@ public:
}

virtual ~SemiImplicitEM() = default;

virtual void Define ( WarpX* const a_WarpX );

virtual bool IsDefined () const { return m_is_defined; }

virtual void PrintParameters () const;

virtual void Initialize () {;}

virtual void GetParticleSolverParams (int& a_max_particle_iterations,
amrex::ParticleReal& a_particle_tolerance )
{
Expand Down Expand Up @@ -65,15 +65,15 @@ private:
bool m_verbose;

WarpX* m_WarpX;

amrex::ParticleReal m_particle_tolerance;
int m_max_particle_iterations;

NonlinearSolverType m_nlsolver_type;
std::unique_ptr<NonlinearSolver<WarpXSolverVec,SemiImplicitEM>> m_nlsolver;

WarpXSolverVec m_E, m_Eold;

void UpdateWarpXState ( const WarpXSolverVec& a_E,
const amrex::Real a_time,
const amrex::Real a_dt );
Expand Down
30 changes: 15 additions & 15 deletions Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,13 +12,13 @@
* xp^{n+1} = xp^n + dt*up^{n+1/2}/(gammap^n + gammap^{n+1})
* up^{n+1} = up^n + dt*qp/mp*(Ep^{n+1/2} + up^{n+1/2}/gammap^{n+1/2} x Bp^{n+1/2})
* where f^{n+1/2} = (f^{n} + f^{n+1})/2.0, for all but Bg, which lives at half steps
*
*
* This algorithm is approximately energy conserving. The violation in energy conservation
* is typically negligible. The advantage of this method over the exactly energy-conserving
* theta-implicit EM method is that light wave dispersion is captured much better. However,
* is typically negligible. The advantage of this method over the exactly energy-conserving
* theta-implicit EM method is that light wave dispersion is captured much better. However,
* the CFL condition for light waves does have to be satisifed for numerical stability.
*
* See G. Chen, L. Chacon, L. Yin, B.J. Albright, D.J. Stark, R.F. Bird,
* See G. Chen, L. Chacon, L. Yin, B.J. Albright, D.J. Stark, R.F. Bird,
* "A semi-implicit energy- and charge-conserving particle-in-cell algorithm for the
* relativistic Vlasov-Maxwell equations.", JCP 407 (2020).
*/
Expand All @@ -28,19 +28,19 @@ void SemiImplicitEM::Define ( WarpX* const a_WarpX )
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
!m_is_defined,
"SemiImplicitEM object is already defined!");

// Retain a pointer back to main WarpX class
m_WarpX = a_WarpX;

// Define E vectors
m_E.Define( m_WarpX->getEfield_fp_vec() );
m_Eold.Define( m_WarpX->getEfield_fp_vec() );

// Need to define the WarpXSolverVec owned dot_mask to do dot
// product correctly for linear and nonlinear solvers
const amrex::Vector<amrex::Geometry>& Geom = m_WarpX->Geom();
m_E.SetDotMask(Geom);

// Parse implicit solver parameters
amrex::ParmParse pp("implicit_evolve");
pp.query("verbose", m_verbose);
Expand Down Expand Up @@ -93,7 +93,7 @@ void SemiImplicitEM::PrintParameters () const
void SemiImplicitEM::OneStep ( const amrex::Real a_old_time,
const amrex::Real a_dt,
const int a_step )
{
{
using namespace amrex::literals;
amrex::ignore_unused(a_step);

Expand All @@ -114,30 +114,30 @@ void SemiImplicitEM::OneStep ( const amrex::Real a_old_time,
// 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, a_old_time, a_dt );

// update WarpX owned Efield_fp and Bfield_fp to t_{n+1/2}
UpdateWarpXState( m_E, a_old_time, a_dt );

// Update field boundary probes prior to updating fields to t_{n+1}
//m_fields->updateBoundaryProbes( a_dt );

// Advance particles to step n+1
m_WarpX->FinishImplicitParticleUpdate();

// Advance fields to step n+1
m_WarpX->FinishImplicitField(m_E.getVec(), m_Eold.getVec(), 0.5);
m_WarpX->UpdateElectricField( m_E, false ); // JRA not sure about false here. is what DG had.

}

void SemiImplicitEM::PreRHSOp ( const WarpXSolverVec& a_E,
const amrex::Real a_time,
const amrex::Real a_dt,
const int a_nl_iter,
const bool a_from_jacobian )
{
{
amrex::ignore_unused(a_E);

// update derived variable B and then update WarpX owned Efield_fp and Bfield_fp
UpdateWarpXState( a_E, a_time, a_dt );

Expand All @@ -152,7 +152,7 @@ void SemiImplicitEM::ComputeRHS ( WarpXSolverVec& a_Erhs,
const WarpXSolverVec& a_E,
const amrex::Real a_time,
const amrex::Real a_dt )
{
{
amrex::ignore_unused(a_E, a_time);
m_WarpX->ComputeRHSE(0.5*a_dt, a_Erhs);
}
Expand Down
12 changes: 6 additions & 6 deletions Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ class ThetaImplicitEM : public ImplicitSolver
{
public:

ThetaImplicitEM()
ThetaImplicitEM()
{
m_nlsolver_type = NonlinearSolverType::Picard;

Expand All @@ -31,15 +31,15 @@ public:
}

virtual ~ThetaImplicitEM() = default;

virtual void Define ( WarpX* const a_WarpX );

virtual bool IsDefined () const { return m_is_defined; }

virtual void PrintParameters () const;

virtual void Initialize () {;}

virtual void GetParticleSolverParams (int& a_max_particle_iterations,
amrex::ParticleReal& a_particle_tolerance )
{
Expand Down Expand Up @@ -80,7 +80,7 @@ private:

WarpXSolverVec m_E, m_Eold;
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > m_Bold;

void UpdateWarpXState ( const WarpXSolverVec& a_E,
const amrex::Real a_time,
const amrex::Real a_dt );
Expand Down
36 changes: 18 additions & 18 deletions Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,11 +12,11 @@
* xp^{n+1} = xp^n + dt*up^{n+1/2}/(gammap^n + gammap^{n+1})
* up^{n+1} = up^n + dt*qp/mp*(Ep^{n+theta} + up^{n+1/2}/gammap^{n+1/2} x Bp^{n+theta})
* where f^{n+theta} = (1.0-theta)*f^{n} + theta*f^{n+1} with 0.5 <= theta <= 1.0
*
* The user-specified time-biasing parameter theta used for the fields on the RHS is bound
*
* The user-specified time-biasing parameter theta used for the fields on the RHS is bound
* between 0.5 and 1.0. The algorithm is exactly energy conserving for theta = 0.5.
* Signifcant damping of high-k modes will occur as theta approaches 1.0. The algorithm is
* numerially stable for any time step. I.e., the CFL condition for light waves does not
* numerially stable for any time step. I.e., the CFL condition for light waves does not
* have to be satisifed and the time step is not limited by the plasma period. However, how
* efficiently the algorithm can use large time steps depends strongly on the nonlinear solver.
* Furthermore, the time step should always be such that particles do not travel outside the
Expand All @@ -25,7 +25,7 @@
*
* See S. Markidis, G. Lapenta, "The energy conserving particle-in-cell method." JCP 230 (2011).
*
* See G. Chen, L. Chacon, D.C. Barnes, "An energy- and charge-conserving, implicit,
* See G. Chen, L. Chacon, D.C. Barnes, "An energy- and charge-conserving, implicit,
* elctrostatic particle-in-cell algorithm." JCP 230 (2011).
*
* See J.R. Angus, A. Link, A. Friedman, D. Ghosh, J. D. Johnson, "On numerical energy
Expand All @@ -42,30 +42,30 @@ void ThetaImplicitEM::Define ( WarpX* const a_WarpX )
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
!m_is_defined,
"ThetaImplicitEM object is already defined!");

// Retain a pointer back to main WarpX class
m_WarpX = a_WarpX;

// Define E vectors
m_E.Define( m_WarpX->getEfield_fp_vec() );
m_Eold.Define( m_WarpX->getEfield_fp_vec() );

// Need to define the WarpXSolverVec owned dot_mask to do dot
// product correctly for linear and nonlinear solvers
const amrex::Vector<amrex::Geometry>& Geom = m_WarpX->Geom();
m_E.SetDotMask(Geom);

// Define Bold vector
const int lev = 0;
m_Bold.resize(1); // size is number of levels
for (int n=0; n<3; n++) {
const amrex::MultiFab& Bfp = m_WarpX->getBfield_fp(lev,n);
m_Bold[lev][n] = std::make_unique<amrex::MultiFab>( Bfp.boxArray(),
m_Bold[lev][n] = std::make_unique<amrex::MultiFab>( Bfp.boxArray(),
Bfp.DistributionMap(),
Bfp.nComp(),
Bfp.nComp(),
Bfp.nGrowVect() );
}

// Parse implicit solver parameters
amrex::ParmParse pp("implicit_evolve");
pp.query("verbose", m_verbose);
Expand Down Expand Up @@ -128,7 +128,7 @@ void ThetaImplicitEM::PrintParameters () const
void ThetaImplicitEM::OneStep ( const amrex::Real a_old_time,
const amrex::Real a_dt,
const int a_step )
{
{
using namespace amrex::literals;
amrex::ignore_unused(a_step);

Expand All @@ -152,31 +152,31 @@ void ThetaImplicitEM::OneStep ( const amrex::Real a_old_time,
// Solve nonlinear system for E at t_{n+theta}
// Particles will be advanced to t_{n+1/2}
m_nlsolver->Solve( m_E, m_Eold, a_old_time, a_dt );

// update WarpX owned Efield_fp and Bfield_fp to t_{n+theta}
UpdateWarpXState( m_E, a_old_time, a_dt );

// Update field boundary probes prior to updating fields to t_{n+1}
//m_fields->updateBoundaryProbes( a_dt );

// Advance particles to step n+1
m_WarpX->FinishImplicitParticleUpdate();

// Advance fields to step n+1
m_WarpX->FinishImplicitField(m_E.getVec(), m_Eold.getVec(), m_theta);
m_WarpX->UpdateElectricField( m_E, false ); // JRA not sure about false here. is what DG had.
m_WarpX->FinishMagneticField( m_Bold, m_theta );

}

void ThetaImplicitEM::PreRHSOp ( const WarpXSolverVec& a_E,
const amrex::Real a_time,
const amrex::Real a_dt,
const int a_nl_iter,
const bool a_from_jacobian )
{
{
amrex::ignore_unused(a_E);

// update derived variable B and then update WarpX owned Efield_fp and Bfield_fp
UpdateWarpXState( a_E, a_time, a_dt );

Expand All @@ -191,7 +191,7 @@ void ThetaImplicitEM::ComputeRHS ( WarpXSolverVec& a_Erhs,
const WarpXSolverVec& a_E,
const amrex::Real a_time,
const amrex::Real a_dt )
{
{
amrex::ignore_unused(a_E, a_time);
m_WarpX->ComputeRHSE(m_theta*a_dt, a_Erhs);
}
Expand Down
Loading

0 comments on commit 1826926

Please sign in to comment.