Skip to content

Commit 1826926

Browse files
[pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci
1 parent a0bbf12 commit 1826926

16 files changed

+208
-208
lines changed

Docs/source/usage/parameters.rst

+1-1
Original file line numberDiff line numberDiff line change
@@ -86,7 +86,7 @@ Overall simulation parameters
8686

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

89-
* ``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.
89+
* ``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.
9090
The version implemented is an updated version that is relativistically correct, including the relativistic gamma factor for the particles.
9191
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
9292
: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.

Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H

+4-4
Original file line numberDiff line numberDiff line change
@@ -26,16 +26,16 @@ public:
2626
//
2727

2828
virtual void Define ( WarpX* const a_WarpX ) = 0;
29-
29+
3030
virtual bool IsDefined () const = 0;
31-
31+
3232
virtual void PrintParameters () const = 0;
3333

3434
virtual void Initialize () = 0;
35-
35+
3636
virtual void GetParticleSolverParams (int& a_max_particle_iter,
3737
amrex::ParticleReal& a_particle_tol ) = 0;
38-
38+
3939
virtual void OneStep ( const amrex::Real a_time,
4040
const amrex::Real a_dt,
4141
const int a_step ) = 0;

Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H

+7-7
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ class SemiImplicitEM : public ImplicitSolver
1818
{
1919
public:
2020

21-
SemiImplicitEM()
21+
SemiImplicitEM()
2222
{
2323
m_nlsolver_type = NonlinearSolverType::Picard;
2424
m_max_particle_iterations = 21;
@@ -28,15 +28,15 @@ public:
2828
}
2929

3030
virtual ~SemiImplicitEM() = default;
31-
31+
3232
virtual void Define ( WarpX* const a_WarpX );
33-
33+
3434
virtual bool IsDefined () const { return m_is_defined; }
35-
35+
3636
virtual void PrintParameters () const;
3737

3838
virtual void Initialize () {;}
39-
39+
4040
virtual void GetParticleSolverParams (int& a_max_particle_iterations,
4141
amrex::ParticleReal& a_particle_tolerance )
4242
{
@@ -65,15 +65,15 @@ private:
6565
bool m_verbose;
6666

6767
WarpX* m_WarpX;
68-
68+
6969
amrex::ParticleReal m_particle_tolerance;
7070
int m_max_particle_iterations;
7171

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

7575
WarpXSolverVec m_E, m_Eold;
76-
76+
7777
void UpdateWarpXState ( const WarpXSolverVec& a_E,
7878
const amrex::Real a_time,
7979
const amrex::Real a_dt );

Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp

+15-15
Original file line numberDiff line numberDiff line change
@@ -12,13 +12,13 @@
1212
* xp^{n+1} = xp^n + dt*up^{n+1/2}/(gammap^n + gammap^{n+1})
1313
* 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})
1414
* where f^{n+1/2} = (f^{n} + f^{n+1})/2.0, for all but Bg, which lives at half steps
15-
*
15+
*
1616
* This algorithm is approximately energy conserving. The violation in energy conservation
17-
* is typically negligible. The advantage of this method over the exactly energy-conserving
18-
* theta-implicit EM method is that light wave dispersion is captured much better. However,
17+
* is typically negligible. The advantage of this method over the exactly energy-conserving
18+
* theta-implicit EM method is that light wave dispersion is captured much better. However,
1919
* the CFL condition for light waves does have to be satisifed for numerical stability.
2020
*
21-
* See G. Chen, L. Chacon, L. Yin, B.J. Albright, D.J. Stark, R.F. Bird,
21+
* See G. Chen, L. Chacon, L. Yin, B.J. Albright, D.J. Stark, R.F. Bird,
2222
* "A semi-implicit energy- and charge-conserving particle-in-cell algorithm for the
2323
* relativistic Vlasov-Maxwell equations.", JCP 407 (2020).
2424
*/
@@ -28,19 +28,19 @@ void SemiImplicitEM::Define ( WarpX* const a_WarpX )
2828
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
2929
!m_is_defined,
3030
"SemiImplicitEM object is already defined!");
31-
31+
3232
// Retain a pointer back to main WarpX class
3333
m_WarpX = a_WarpX;
34-
34+
3535
// Define E vectors
3636
m_E.Define( m_WarpX->getEfield_fp_vec() );
3737
m_Eold.Define( m_WarpX->getEfield_fp_vec() );
38-
38+
3939
// Need to define the WarpXSolverVec owned dot_mask to do dot
4040
// product correctly for linear and nonlinear solvers
4141
const amrex::Vector<amrex::Geometry>& Geom = m_WarpX->Geom();
4242
m_E.SetDotMask(Geom);
43-
43+
4444
// Parse implicit solver parameters
4545
amrex::ParmParse pp("implicit_evolve");
4646
pp.query("verbose", m_verbose);
@@ -93,7 +93,7 @@ void SemiImplicitEM::PrintParameters () const
9393
void SemiImplicitEM::OneStep ( const amrex::Real a_old_time,
9494
const amrex::Real a_dt,
9595
const int a_step )
96-
{
96+
{
9797
using namespace amrex::literals;
9898
amrex::ignore_unused(a_step);
9999

@@ -114,30 +114,30 @@ void SemiImplicitEM::OneStep ( const amrex::Real a_old_time,
114114
// Solve nonlinear system for E at t_{n+1/2}
115115
// Particles will be advanced to t_{n+1/2}
116116
m_nlsolver->Solve( m_E, m_Eold, a_old_time, a_dt );
117-
117+
118118
// update WarpX owned Efield_fp and Bfield_fp to t_{n+1/2}
119119
UpdateWarpXState( m_E, a_old_time, a_dt );
120120

121121
// Update field boundary probes prior to updating fields to t_{n+1}
122122
//m_fields->updateBoundaryProbes( a_dt );
123-
123+
124124
// Advance particles to step n+1
125125
m_WarpX->FinishImplicitParticleUpdate();
126126

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

133133
void SemiImplicitEM::PreRHSOp ( const WarpXSolverVec& a_E,
134134
const amrex::Real a_time,
135135
const amrex::Real a_dt,
136136
const int a_nl_iter,
137137
const bool a_from_jacobian )
138-
{
138+
{
139139
amrex::ignore_unused(a_E);
140-
140+
141141
// update derived variable B and then update WarpX owned Efield_fp and Bfield_fp
142142
UpdateWarpXState( a_E, a_time, a_dt );
143143

@@ -152,7 +152,7 @@ void SemiImplicitEM::ComputeRHS ( WarpXSolverVec& a_Erhs,
152152
const WarpXSolverVec& a_E,
153153
const amrex::Real a_time,
154154
const amrex::Real a_dt )
155-
{
155+
{
156156
amrex::ignore_unused(a_E, a_time);
157157
m_WarpX->ComputeRHSE(0.5*a_dt, a_Erhs);
158158
}

Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H

+6-6
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ class ThetaImplicitEM : public ImplicitSolver
1818
{
1919
public:
2020

21-
ThetaImplicitEM()
21+
ThetaImplicitEM()
2222
{
2323
m_nlsolver_type = NonlinearSolverType::Picard;
2424

@@ -31,15 +31,15 @@ public:
3131
}
3232

3333
virtual ~ThetaImplicitEM() = default;
34-
34+
3535
virtual void Define ( WarpX* const a_WarpX );
36-
36+
3737
virtual bool IsDefined () const { return m_is_defined; }
38-
38+
3939
virtual void PrintParameters () const;
4040

4141
virtual void Initialize () {;}
42-
42+
4343
virtual void GetParticleSolverParams (int& a_max_particle_iterations,
4444
amrex::ParticleReal& a_particle_tolerance )
4545
{
@@ -80,7 +80,7 @@ private:
8080

8181
WarpXSolverVec m_E, m_Eold;
8282
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > m_Bold;
83-
83+
8484
void UpdateWarpXState ( const WarpXSolverVec& a_E,
8585
const amrex::Real a_time,
8686
const amrex::Real a_dt );

Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp

+18-18
Original file line numberDiff line numberDiff line change
@@ -12,11 +12,11 @@
1212
* xp^{n+1} = xp^n + dt*up^{n+1/2}/(gammap^n + gammap^{n+1})
1313
* up^{n+1} = up^n + dt*qp/mp*(Ep^{n+theta} + up^{n+1/2}/gammap^{n+1/2} x Bp^{n+theta})
1414
* where f^{n+theta} = (1.0-theta)*f^{n} + theta*f^{n+1} with 0.5 <= theta <= 1.0
15-
*
16-
* The user-specified time-biasing parameter theta used for the fields on the RHS is bound
15+
*
16+
* The user-specified time-biasing parameter theta used for the fields on the RHS is bound
1717
* between 0.5 and 1.0. The algorithm is exactly energy conserving for theta = 0.5.
1818
* Signifcant damping of high-k modes will occur as theta approaches 1.0. The algorithm is
19-
* numerially stable for any time step. I.e., the CFL condition for light waves does not
19+
* numerially stable for any time step. I.e., the CFL condition for light waves does not
2020
* have to be satisifed and the time step is not limited by the plasma period. However, how
2121
* efficiently the algorithm can use large time steps depends strongly on the nonlinear solver.
2222
* Furthermore, the time step should always be such that particles do not travel outside the
@@ -25,7 +25,7 @@
2525
*
2626
* See S. Markidis, G. Lapenta, "The energy conserving particle-in-cell method." JCP 230 (2011).
2727
*
28-
* See G. Chen, L. Chacon, D.C. Barnes, "An energy- and charge-conserving, implicit,
28+
* See G. Chen, L. Chacon, D.C. Barnes, "An energy- and charge-conserving, implicit,
2929
* elctrostatic particle-in-cell algorithm." JCP 230 (2011).
3030
*
3131
* See J.R. Angus, A. Link, A. Friedman, D. Ghosh, J. D. Johnson, "On numerical energy
@@ -42,30 +42,30 @@ void ThetaImplicitEM::Define ( WarpX* const a_WarpX )
4242
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
4343
!m_is_defined,
4444
"ThetaImplicitEM object is already defined!");
45-
45+
4646
// Retain a pointer back to main WarpX class
4747
m_WarpX = a_WarpX;
48-
48+
4949
// Define E vectors
5050
m_E.Define( m_WarpX->getEfield_fp_vec() );
5151
m_Eold.Define( m_WarpX->getEfield_fp_vec() );
52-
52+
5353
// Need to define the WarpXSolverVec owned dot_mask to do dot
5454
// product correctly for linear and nonlinear solvers
5555
const amrex::Vector<amrex::Geometry>& Geom = m_WarpX->Geom();
5656
m_E.SetDotMask(Geom);
57-
57+
5858
// Define Bold vector
5959
const int lev = 0;
6060
m_Bold.resize(1); // size is number of levels
6161
for (int n=0; n<3; n++) {
6262
const amrex::MultiFab& Bfp = m_WarpX->getBfield_fp(lev,n);
63-
m_Bold[lev][n] = std::make_unique<amrex::MultiFab>( Bfp.boxArray(),
63+
m_Bold[lev][n] = std::make_unique<amrex::MultiFab>( Bfp.boxArray(),
6464
Bfp.DistributionMap(),
65-
Bfp.nComp(),
65+
Bfp.nComp(),
6666
Bfp.nGrowVect() );
6767
}
68-
68+
6969
// Parse implicit solver parameters
7070
amrex::ParmParse pp("implicit_evolve");
7171
pp.query("verbose", m_verbose);
@@ -128,7 +128,7 @@ void ThetaImplicitEM::PrintParameters () const
128128
void ThetaImplicitEM::OneStep ( const amrex::Real a_old_time,
129129
const amrex::Real a_dt,
130130
const int a_step )
131-
{
131+
{
132132
using namespace amrex::literals;
133133
amrex::ignore_unused(a_step);
134134

@@ -152,31 +152,31 @@ void ThetaImplicitEM::OneStep ( const amrex::Real a_old_time,
152152
// Solve nonlinear system for E at t_{n+theta}
153153
// Particles will be advanced to t_{n+1/2}
154154
m_nlsolver->Solve( m_E, m_Eold, a_old_time, a_dt );
155-
155+
156156
// update WarpX owned Efield_fp and Bfield_fp to t_{n+theta}
157157
UpdateWarpXState( m_E, a_old_time, a_dt );
158158

159159
// Update field boundary probes prior to updating fields to t_{n+1}
160160
//m_fields->updateBoundaryProbes( a_dt );
161-
161+
162162
// Advance particles to step n+1
163163
m_WarpX->FinishImplicitParticleUpdate();
164164

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

172172
void ThetaImplicitEM::PreRHSOp ( const WarpXSolverVec& a_E,
173173
const amrex::Real a_time,
174174
const amrex::Real a_dt,
175175
const int a_nl_iter,
176176
const bool a_from_jacobian )
177-
{
177+
{
178178
amrex::ignore_unused(a_E);
179-
179+
180180
// update derived variable B and then update WarpX owned Efield_fp and Bfield_fp
181181
UpdateWarpXState( a_E, a_time, a_dt );
182182

@@ -191,7 +191,7 @@ void ThetaImplicitEM::ComputeRHS ( WarpXSolverVec& a_Erhs,
191191
const WarpXSolverVec& a_E,
192192
const amrex::Real a_time,
193193
const amrex::Real a_dt )
194-
{
194+
{
195195
amrex::ignore_unused(a_E, a_time);
196196
m_WarpX->ComputeRHSE(m_theta*a_dt, a_Erhs);
197197
}

0 commit comments

Comments
 (0)