-
Notifications
You must be signed in to change notification settings - Fork 200
/
Copy pathImplicitSolver.H
178 lines (143 loc) · 6.04 KB
/
ImplicitSolver.H
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
/* Copyright 2024 Justin Angus, Debojyoti Ghosh
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#ifndef Implicit_Solver_H_
#define Implicit_Solver_H_
#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H"
#include "NonlinearSolvers/NonlinearSolverLibrary.H"
#include "Utils/WarpXAlgorithmSelection.H"
#include <AMReX_Array.H>
#include <AMReX_REAL.H>
#include <AMReX_LO_BCTYPES.H>
/**
* \brief Base class for implicit time solvers. The base functions are those
* needed by an implicit solver to be used through WarpX and those needed
* to use the nonlinear solvers, such as Picard or Newton (i.e., JFNK).
*/
class WarpX;
class ImplicitSolver
{
public:
ImplicitSolver() = default;
virtual ~ImplicitSolver() = default;
// Prohibit Move and Copy operations
ImplicitSolver(const ImplicitSolver&) = delete;
ImplicitSolver& operator=(const ImplicitSolver&) = delete;
ImplicitSolver(ImplicitSolver&&) = delete;
ImplicitSolver& operator=(ImplicitSolver&&) = delete;
//
// the following routines are called by WarpX
//
/**
* \brief Read user-provided parameters that control the implicit solver.
* Allocate internal arrays for intermediate field values needed by the solver.
*/
virtual void Define ( WarpX* a_WarpX ) = 0;
[[nodiscard]] bool IsDefined () const { return m_is_defined; }
virtual void PrintParameters () const = 0;
void GetParticleSolverParams (int& a_max_particle_iter,
amrex::ParticleReal& a_particle_tol ) const
{
a_max_particle_iter = m_max_particle_iterations;
a_particle_tol = m_particle_tolerance;
}
void CreateParticleAttributes () const;
/**
* \brief Advance fields and particles by one time step using the specified implicit algorithm
*/
virtual void OneStep ( amrex::Real a_time,
amrex::Real a_dt,
int a_step ) = 0;
//
// the following routines are called by the linear and nonlinear solvers
//
/**
* \brief Computes the RHS of the equation corresponding to the specified implicit algorithm.
* The discrete equations corresponding to numerical integration of ODEs are often
* written in the form U = b + RHS(U), where U is the variable to be solved for (e.g.,
* the solution at the next time step), b is a constant (i.e., the solution from the
* previous time step), and RHS(U) is the right-hand-side of the equation. Iterative
* solvers, such as Picard and Newton, and higher-order Runge-Kutta methods, need to
* compute RHS(U) multiple times for different values of U. Thus, a routine that
* returns this value is needed.
* e.g., Ebar - E^n = cvac^2*0.5*dt*(curl(Bbar) - mu0*Jbar(Ebar,Bbar))
* Here, U = Ebar, b = E^n, and the expression on the right is RHS(U).
*/
virtual void ComputeRHS ( WarpXSolverVec& a_RHS,
const WarpXSolverVec& a_E,
amrex::Real a_time,
int a_nl_iter,
bool a_from_jacobian ) = 0;
[[nodiscard]] int numAMRLevels () const { return m_num_amr_levels; }
[[nodiscard]] const amrex::Geometry& GetGeometry (int) const;
[[nodiscard]] const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& GetFieldBoundaryLo () const;
[[nodiscard]] const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& GetFieldBoundaryHi () const;
[[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> GetLinOpBCLo () const;
[[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> GetLinOpBCHi () const;
[[nodiscard]] amrex::Real GetTheta () const { return m_theta; }
protected:
/**
* \brief Pointer back to main WarpX class
*/
WarpX* m_WarpX;
bool m_is_defined = false;
/**
* \brief Number of AMR levels
*/
int m_num_amr_levels = 1;
/**
* \brief Time step
*/
mutable amrex::Real m_dt = 0.0;
/**
* \brief Time-biasing parameter for fields used on RHS to advance system.
*/
amrex::Real m_theta = 0.5;
/**
* \brief Nonlinear solver type and object
*/
NonlinearSolverType m_nlsolver_type;
std::unique_ptr<NonlinearSolver<WarpXSolverVec,ImplicitSolver>> m_nlsolver;
/**
* \brief tolerance used by the iterative method used to obtain a self-consistent
* update of the particle positions and velocities for given E and B on the grid
*/
amrex::ParticleReal m_particle_tolerance = 1.0e-10;
/**
* \brief maximum iterations for the iterative method used to obtain a self-consistent
* update of the particle positions and velocities for given E and B on the grid
*/
int m_max_particle_iterations = 21;
/**
* \brief parse nonlinear solver parameters (if one is used)
*/
void parseNonlinearSolverParams( const amrex::ParmParse& pp )
{
std::string nlsolver_type_str;
pp.get("nonlinear_solver", nlsolver_type_str);
if (nlsolver_type_str=="picard") {
m_nlsolver_type = NonlinearSolverType::Picard;
m_nlsolver = std::make_unique<PicardSolver<WarpXSolverVec,ImplicitSolver>>();
m_max_particle_iterations = 1;
m_particle_tolerance = 0.0;
}
else if (nlsolver_type_str=="newton") {
m_nlsolver_type = NonlinearSolverType::Newton;
m_nlsolver = std::make_unique<NewtonSolver<WarpXSolverVec,ImplicitSolver>>();
pp.query("max_particle_iterations", m_max_particle_iterations);
pp.query("particle_tolerance", m_particle_tolerance);
}
else {
WARPX_ABORT_WITH_MESSAGE(
"invalid nonlinear_solver specified. Valid options are picard and newton.");
}
}
/**
* \brief Convert from WarpX FieldBoundaryType to amrex::LinOpBCType
*/
[[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> convertFieldBCToLinOpBC ( const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& ) const;
};
#endif