Chombo + EB  3.2
ARK4.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 #ifndef _ARK4_H_
12 #define _ARK4_H_
13 
14 /// 4th-order additive Runge-Kutta algorithm
15 /**
16  This templated class encapsulates
17  the fourth-order additive Runge-Kutta method
18  "ARK4(3)6L[2]SA"
19  by Kennedy and Carpenter 2003 Appl. Numer. Math. 44: 139-181
20 
21  See also section 3 of Zhang, Johansen, and Colella,
22  SIAM J. Sci. Comput. 34, pp. B179-B201.
23 */
24 
25 // Forward declaration
26 class dXConcept;
27 
28 // \class XConcept
29 // Concept of class X, the type of your state data
30 class XConcept
31 {
32  // allocate space for this
33  void define(const XConcept& a_state);
34 
35  // copy data to this
36  void copy(const XConcept& a_state);
37 
38  // increment this by a_factor * a_increment
39  void increment(const dXConcept& a_increment, Real a_factor = 1.);
40 
41  ~XConcept() {};
42 };
43 
44 // \class dXConcept
45 // Concept of class dX, same type as the operator evaluations
46 class dXConcept
47 {
48  // allocate space for this
49  void define(const XConcept& a_state);
50 
51  // initialize with given X
52  void init(const XConcept& a_state);
53 
54  // increment this by a_factor * a_increment
55  void increment(const dXConcept& a_increment, Real a_factor = 1.);
56 
57  ~dXConcept() {};
58 };
59 
60 // \class FEConcept
61 // Concept of class FE, which defines your explicit operator
62 class FEConcept
63 {
64  // allocate space and define infrastructure for this
65  void define(const XConcept& a_state, Real a_dt);
66 
67  // reset the timestep
68  void resetDt(Real a_dt);
69 
70  // apply the operator, including multiplication by timestep
71  void operator()(dXConcept& a_result, Real a_time, const XConcept& a_state);
72 
73  ~FEConcept() {};
74 };
75 
76 // \class FIConcept
77 // Concept of class FI, which defines your implicit operator
78 class FIConcept
79 {
80  // allocate space and define infrastructure for this
81  void define(const XConcept& a_state, Real a_dt, Real a_gamma);
82 
83  // reset the timestep
84  void resetDt(Real a_dt);
85 
86  // apply the operator, including multiplication by timestep
87  void operator()(dXConcept& a_result, Real a_time, const XConcept& a_state);
88 
89  // solve for a_soln in (I - s_gamma * m_dt * FI(a_time))(a_soln) = a_rhs
90  void solve(XConcept& a_soln, const dXConcept& a_rhs, Real a_time);
91 
92  ~FIConcept() {};
93 };
94 
95 template <class X, class FI, class FE, class dX>
96 class ARK4
97 {
98 public:
99 
100  ARK4<X, FI, FE, dX>() { m_isDefined = false; }
101 
102  // This must be called first.
103  void define(const X& a_state, Real a_dt);
104 
105  // Advance one step.
106  void advance(Real a_time, X& a_state);
107 
108  // Reset the timestep.
109  void resetDt(Real a_dt);
110 
111  bool isDefined() const { return m_isDefined; }
112 
113  /// Runge-Kutta coefficients
114  static const int s_nStages = 6;
115  static const Real s_gamma;
116  static const Real s_aE[s_nStages][s_nStages];
117  static const Real s_aI[s_nStages][s_nStages];
118  static const Real s_b[s_nStages];
119  static const Real s_c[s_nStages];
120 
121 protected:
122 
125  X m_phi[s_nStages];
126  dX m_kE[s_nStages-1];
127  dX m_kI[s_nStages-1];
129  dX m_rhs;
130  FI m_fI;
131  FE m_fE;
132 
133 private:
134 
135  // Make sure the concepts match the API
137 };
138 
139 //==============================================
140 
141 template <class X, class FI, class FE, class dX>
142 void ARK4<X, FI, FE, dX>::define(const X& a_state, Real a_dt)
143 {
144  m_dt = a_dt;
145  // define X
146  for (int stage = 0; stage < s_nStages; stage++)
147  {
148  m_phi[stage].define(a_state);
149  }
150  // define dX
151  for (int stage = 0; stage < s_nStages-1; stage++)
152  {
153  m_kE[stage].define(a_state);
154  m_kI[stage].define(a_state);
155  }
156  m_kEfinal.define(a_state);
157  m_rhs.define(a_state);
158  // define fI
159  m_fI.define(a_state, m_dt, s_gamma);
160  // define fE
161  m_fE.define(a_state, m_dt);
162  m_isDefined = true;
163 }
164 
165 /*
166  Reset the timestep.
167  */
168 template <class X, class FI, class FE, class dX>
170 {
171  CH_assert(isDefined());
172 
173  // Only update everything if dt has changed
174  Real reltol = 1e-10;
175  if (Abs(m_dt - a_dt) > m_dt*reltol)
176  {
177  m_dt = a_dt;
178  m_fI.resetDt(m_dt);
179  m_fE.resetDt(m_dt);
180  }
181 }
182 
183 /*
184  Advance solution a_state in time, a_time to a_time + a_dt.
185  */
186 template <class X, class FI, class FE, class dX>
187 void ARK4<X, FI, FE, dX>::advance(Real a_time, X& a_state)
188 {
189  CH_assert(isDefined());
190  // Set intermediate times t from a_time and m_dt.
191  Real t[s_nStages];
192  for (int stage = 0; stage < s_nStages; stage++)
193  {
194  t[stage] = a_time + s_c[stage] * m_dt;
195  }
196 
197  // Set m_phi[0] := a_state.
198  m_phi[0].copy(a_state);
199  for (int stage = 1; stage < s_nStages; stage++)
200  {
201  m_fE(m_kE[stage-1], t[stage-1], m_phi[stage-1]);
202  m_fI(m_kI[stage-1], t[stage-1], m_phi[stage-1]);
203 
204  // Set m_rhs := a_state +
205  // sum_{j=0:stage-1} ( s_aE[stage][j] * FE(m_phi[j], t[j]) +
206  // s_aI[stage][j] * FI(m_phi[j]) ).
207  m_rhs.init(a_state);
208  for (int j = 0; j < stage; j++)
209  {
210  m_rhs.increment(m_kE[j], s_aE[stage][j]);
211  m_rhs.increment(m_kI[j], s_aI[stage][j]);
212  }
213  // Solve for m_phi[stage] in
214  // (I - s_gamma * m_dt * FI) (m_phi[stage]) = m_rhs.
215  m_phi[stage].define(a_state);
216  m_fI.solve(m_phi[stage], m_rhs, t[stage]);
217  }
218 
219  // Set a_state := m_phi[s_nStages-1] +
220  // sum_{j=0:s_nStages-1} ( (s_b[j] - s_aE[s_nStages-1][j]) *
221  // FE(m_phi[j], t[s_nStages-1]) ).
222  // Note these m_fE evaluations are for m_phi[j] at time t[s_nStages-1],
223  // not time t[j], so you cannot reuse m_kE[j].
224  const int stageFinal = s_nStages-1;
225  Real tFinal = t[stageFinal];
226  // Calculate the final stage explicit operator
227  m_fE(m_kEfinal, tFinal, m_phi[stageFinal]);
228  a_state.copy(m_phi[stageFinal]);
229  for (int j = 0; j < s_nStages-1; j++)
230  {
231  a_state.increment(m_kE[j], s_b[j] - s_aE[stageFinal][j]);
232  }
233  // Add in the final stage explicit operator
234  a_state.increment(m_kEfinal, s_b[stageFinal] /* - s_aE, which is 0 */);
235 }
236 
237 template<class X, class FI, class FE, class dX>
239 
240 // Time coefficients for each stage
241 template<class X, class FI, class FE, class dX>
242 const Real ARK4<X, FI, FE, dX>::s_c[] = { 0.0, 0.5, 0.332, 0.62, 0.85, 1.0 };
243 
244 // Stage coefficients - each row is for that stage
245 template<class X, class FI, class FE, class dX>
247  {0., 0., 0., 0., 0., 0.},
248  {0.5, 0., 0., 0., 0., 0.},
249  {0.221776, 0.110224, 0., 0., 0., 0.},
250  {-0.04884659515311857, -0.17772065232640102, 0.8465672474795197, 0., 0., 0.},
251  {-0.15541685842491548, -0.3567050098221991, 1.0587258798684427, 0.30339598837867193, 0., 0.},
252  { 0.2014243506726763, 0.008742057842904185, 0.15993995707168115, 0.4038290605220775, 0.22606457389066084, 0.}
253 };
254 
255 template<class X, class FI, class FE, class dX>
257  {0., 0., 0., 0., 0., 0.},
258  {0.25, 0.25, 0., 0., 0., 0.},
259  {0.137776, -0.055776, 0.25, 0., 0., 0.},
260  {0.14463686602698217, -0.22393190761334475, 0.4492950415863626, 0.25, 0., 0.},
261  {0.09825878328356477, -0.5915442428196704, 0.8101210538282996, 0.283164405707806, 0.25, 0.},
262  {0.15791629516167136, 0., 0.18675894052400077, 0.6805652953093346, -0.27524053099500667, 0.25}
263 };
264 
265 template<class X, class FI, class FE, class dX>
267  {0.15791629516167136, 0., 0.18675894052400077, 0.6805652953093346, -0.27524053099500667, 0.25};
268 
269 #endif
void copy(const XConcept &a_state)
#define CH_assert(cond)
Definition: CHArray.H:37
void define(const XConcept &a_state)
bool m_isDefined
Definition: ARK4.H:123
void increment(const dXConcept &a_increment, Real a_factor=1.)
dX m_kEfinal
Definition: ARK4.H:128
Definition: ARK4.H:78
dX m_rhs
Definition: ARK4.H:129
FI m_fI
Definition: ARK4.H:130
~XConcept()
Definition: ARK4.H:41
Definition: ARK4.H:62
void advance(Real a_time, X &a_state)
Definition: ARK4.H:187
void define(const X &a_state, Real a_dt)
Definition: ARK4.H:142
Real m_dt
Definition: ARK4.H:124
double Real
Definition: REAL.H:33
static ARK4< XConcept, FIConcept, FEConcept, dXConcept > testConcept
Definition: ARK4.H:136
~FIConcept()
Definition: ARK4.H:92
T Abs(const T &a_a)
Definition: Misc.H:53
~dXConcept()
Definition: ARK4.H:57
bool isDefined() const
Definition: ARK4.H:111
Definition: ARK4.H:30
FE m_fE
Definition: ARK4.H:131
static const Real s_gamma
Definition: ARK4.H:115
~FEConcept()
Definition: ARK4.H:73
Definition: ARK4.H:96
void resetDt(Real a_dt)
Definition: ARK4.H:169
Definition: ARK4.H:46