Chombo + EB  3.2
ARK4DenseOutput.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 _ARK4DENSEOUTPUT_H_
12 #define _ARK4DENSEOUTPUT_H_
13 
14 
15 #include "NamespaceHeader.H"
16 
17 /// This is a more flexible, reduced memory version
18 /// 4th-order additive Runge-Kutta algorithm,
19 /// that also includes dense output coefficients
20 /**
21  This templated class encapsulates
22  the fourth-order additive Runge-Kutta method
23  "ARK4(3)6L[2]SA"
24  by Kennedy and Carpenter 2003 Appl. Numer. Math. 44: 139-181
25 
26  See also section 3 of Zhang, Johansen, and Colella,
27  SIAM J. Sci. Comput. 34, pp. B179-B201.
28 */
29 
30 // Soln is the type that operators are applied to
31 // Rhs is the type that operator produce
32 // IMEXOp wraps the explicit and implicit operators
33 template <class Soln, class Rhs, class IMEXOp>
35 {
36 public:
37 
39 
40  // This must be called first.
41  void define(const Soln& a_state, Real a_dt, bool a_denseOutput = false);
42 
43  // Advance one step.
44  void advance(Real a_time, Soln& a_state);
45 
46  // Return current dense output coefs, 0th power first, etc.
47  // NOTE: These are intended to be for <JU> conserved quantities
48  void denseOutputCoefs(Vector<Rhs*>& a_interpCoefs);
49 
50  // Reset the timestep.
51  void resetDt(Real a_dt);
52 
53  // Set whether the step starts at 0. and/or ends at 1.
54  void start0end1(bool a_start0, bool a_end1);
55 
56  // Access to the operators if we're defined already
57  // Caller is responsible for making sure they're in a good state
58  // (for example, like have resetDt() called, etc.
59  IMEXOp& getImExOp();
60 
61  bool isDefined() const { return m_isDefined; }
62 
63  bool hasDenseOutput() const { return m_hasDenseOutput; }
64 
65  /// Runge-Kutta coefficients
66  static const int s_nStages = 6;
67  static const Real s_aIdiag;
68  static const Real s_c[s_nStages];
69  static const Real s_aE[s_nStages][s_nStages];
70  static const Real s_aI[s_nStages][s_nStages];
71  static const Real s_b[s_nStages];
72  static const int s_nDenseCoefs = 3;
74 
75 protected:
78  bool m_hasDenseOutput; // If there is dense output data to interpolate
82  Rhs m_rhs;
84  Rhs m_kE;
85  Rhs m_kI;
86  IMEXOp m_opImEx;
87 
88 private:
89 
90 };
91 
92 //==============================================
93 
94 template <class Soln, class Rhs, class IMEXOp>
96 define(const Soln& a_state, Real a_dt, bool a_denseOutput)
97 {
98  CH_TIMERS("ARK4DenseOutput::define");
99  m_dt = a_dt;
100  m_denseOutput = a_denseOutput;
101  // define Soln types
102  for (int stage = 0; stage < s_nStages; stage++)
103  m_phi[stage].define(a_state);
104 
105  // define Rhs types
106  m_kE.define(a_state);
107  m_kI.define(a_state);
108  m_rhs.define(a_state);
109 
110  // define opImEx
111  m_opImEx.define(a_state, m_dt, s_aIdiag);
112 
113  // if dense output is requested, need more storage
114  for (int coef = 0; m_denseOutput && (coef < s_nDenseCoefs); coef++)
115  m_denseCoefs[coef].define(a_state);
116 
117  m_hasDenseOutput = false;
118  m_isDefined = true;
119 }
120 
121 /*
122  Get a reference to the implicit-explicit operator
123  */
124 template <class Soln, class Rhs, class IMEXOp>
127 {
128  return m_opImEx;
129 }
130 
131 /*
132  Reset the timestep.
133  */
134 template <class Soln, class Rhs, class IMEXOp>
137 {
138  CH_assert(isDefined());
139 
140  // Only update everything if dt has changed
141  Real reltol = 1e-14;
142  if (Abs(m_dt - a_dt) > m_dt*reltol)
143  {
144  m_dt = a_dt;
145  m_hasDenseOutput = false;
146  m_opImEx.resetDt(m_dt);
147  }
148 }
149 
150 /*
151  Set whether the step starts at 0. and/or ends at 1.
152  */
153 template <class Soln, class Rhs, class IMEXOp>
155 start0end1(bool a_start0, bool a_end1)
156 {
157  CH_assert(isDefined());
158 
159  int stage0 = -1; // if no stage is at time 0.
160  if (a_start0) stage0 = 0;
161 
162  int stage1 = -1; // if no stage is at time 1.
163  if (a_end1) stage1 = s_nStages-1;
164 
165  m_opImEx.stage0stage1(stage0, stage1);
166 }
167 
168 /*
169  Advance solution a_state in time, a_time to a_time + a_dt.
170  */
171 template <class Soln, class Rhs, class IMEXOp>
173 advance(Real a_time, Soln& a_state)
174 {
175  CH_TIMERS("ARK4DenseOutput::advance");
176  CH_assert(isDefined());
177 
178  CH_TIMER("ARK4DenseOutput::advance - explicit op", t1);
179  CH_TIMER("ARK4DenseOutput::advance - implicit op", t2);
180  CH_TIMER("ARK4DenseOutput::advance - implicit solve", t3);
181 
182  // Reset the dense output coefs
183  if (m_denseOutput)
184  {
185  m_hasDenseOutput = false;
186  for (int icoef=0; icoef < s_nDenseCoefs; ++icoef)
187  m_denseCoefs[icoef].zero();
188  }
189 
190  // Copy a_state into all stages
191  for (int stage = 0; stage < s_nStages; stage++)
192  m_phi[stage].copy(a_state);
193 
194  // Copy a_state values back to itself, to reset any internal state
195  // (like fine flux registers for this level)
196  a_state.copy(m_phi[0]);
197 
198  // For each stage
199  for (int stage = 0; stage < s_nStages; stage++)
200  {
201  Real t = a_time + s_c[stage]*m_dt;
202  if (stage > 0)
203  {
204  // Do the solve - copy rhs from phi[stage]
205  m_rhs.copy(m_phi[stage]);
206  // Solve for m_phi[stage] in
207  // (I - s_aIdiag * m_dt * FI) (m_phi[stage]) = m_rhs.
208  CH_START(t3);
209  m_opImEx.solve(m_phi[stage], m_rhs, t, stage);
210  CH_STOP(t3);
211  }
212 
213  // Calculate the operators for this stage
214  CH_START(t1);
215  m_opImEx.explicitOp(m_kE, t, m_phi[stage], stage);
216  CH_STOP(t1);
217  // Accumulate the explicit op into the stage rhs, final solution
218  // NOTE: increment any internal registers
219  a_state.increment(m_kE, m_dt*s_b[stage], true);
220  for (int k=stage+1; k < s_nStages; ++k)
221  m_phi[k].increment(m_kE, m_dt*s_aE[k][stage]);
222 
223  CH_START(t2);
224  m_opImEx.implicitOp(m_kI, t, m_phi[stage], stage);
225  CH_STOP(t2);
226 
227  // Now accumulate the implicit op into the stage rhs, final solution
228  a_state.increment(m_kI, m_dt*s_b[stage], true);
229  for (int k=stage+1; k < s_nStages; ++k)
230  m_phi[k].increment(m_kI, m_dt*s_aI[k][stage]);
231 
232  // This might provide an optimization
233  // Accumulate the final solution diff from last stage, explicit op only
234  // a_state.increment(m_kE, m_dt*(s_b[stage] - s_aE[s_nStages-1][stage]));
235 
236  if (m_denseOutput)
237  {
238  // pout() << "Stage: " << stage-1 << endl;
239  for (int icoef=0; icoef < s_nDenseCoefs; ++icoef)
240  {
241  m_denseCoefs[icoef].increment(m_kE, m_dt*s_bstar[icoef][stage]);
242  m_denseCoefs[icoef].increment(m_kI, m_dt*s_bstar[icoef][stage]);
243 
244  /*
245  LevelData<FArrayBox>& coef = m_denseCoefs[icoef].data();
246  DataIterator dit(coef.getBoxes());
247  for (dit.begin(); dit.ok(); ++dit)
248  pout() << " Coef[" << icoef+1 << "] = " << coef[dit].min() << endl;
249  */
250  }
251  }
252  }
253 
254  m_hasDenseOutput = m_denseOutput;
255  m_time = a_time;
256 }
257 
258 /*
259  Return the coefs to interpolate solution, in terms of power of the fraction
260  of time between t_old and t_new.
261  */
262 template <class Soln, class Rhs, class IMEXOp>
265 {
266  CH_TIMERS("ARK4DenseOutput::denseOutputCoefs");
267 
268  const int nCoef = s_nDenseCoefs+1;
269  CH_assert(m_hasDenseOutput);
270  CH_assert(a_interpCoefs.size() == nCoef);
271 
272  for (int icoef=0; icoef < nCoef; ++icoef)
273  CH_assert(a_interpCoefs[icoef] != NULL);
274 
275  // Copy over the dense coef values
276 
277  // First coeficient is just the old state
278  a_interpCoefs[0]->copy(m_phi[0]);
279 
280  // Next coefs are our dense output
281  for (int icoef = 1; icoef < nCoef ; ++icoef)
282  {
283  /*
284  LevelData<FArrayBox>& coef = m_denseCoefs[icoef].data();
285  for (dit.begin(); dit.ok(); ++dit)
286  pout() << " Coef[" << icoef+1 << "] = " << coef[dit].min() << endl;
287  */
288  a_interpCoefs[icoef]->copy(m_denseCoefs[icoef-1]);
289  }
290 }
291 
292 
293 /*
294  Static constants for ARK4
295  */
296 
297 template <class Soln, class Rhs, class IMEXOp>
299 s_aIdiag = 0.25;
300 
301 // Time coefficients for each stage
302 template <class Soln, class Rhs, class IMEXOp>
304 s_c[] = { 0.0, 0.5, 0.332, 0.62, 0.85, 1.0 };
305 
306 // Stage coefficients - each row is for that stage
307 template <class Soln, class Rhs, class IMEXOp>
310  {0., 0., 0., 0., 0., 0.},
311  {0.5, 0., 0., 0., 0., 0.},
312  {0.221776, 0.110224, 0., 0., 0., 0.},
313  {-0.04884659515311857, -0.17772065232640102, 0.8465672474795197, 0., 0., 0.},
314  {-0.15541685842491548, -0.3567050098221991, 1.0587258798684427, 0.30339598837867193, 0., 0.},
315  { 0.2014243506726763, 0.008742057842904185, 0.15993995707168115, 0.4038290605220775, 0.22606457389066084, 0.}
316 };
317 
318 // Implicit stage coefficients
319 template <class Soln, class Rhs, class IMEXOp>
322  {0., 0., 0., 0., 0., 0.},
323  {0.25, 0.25, 0., 0., 0., 0.},
324  {0.137776, -0.055776, 0.25, 0., 0., 0.},
325  {0.14463686602698217, -0.22393190761334475, 0.4492950415863626, 0.25, 0., 0.},
326  {0.09825878328356477, -0.5915442428196704, 0.8101210538282996, 0.283164405707806, 0.25, 0.},
327  {0.15791629516167136, 0., 0.18675894052400077, 0.6805652953093346, -0.27524053099500667, 0.25}
328 };
329 
330 // Final stage coefficients
331 template <class Soln, class Rhs, class IMEXOp>
333 s_b[] =
334  {0.15791629516167136, 0., 0.18675894052400077, 0.6805652953093346, -0.27524053099500667, 0.25};
335 
336 // Coefficients for dense ouput, 4th-order interpolation
337 template <class Soln, class Rhs, class IMEXOp>
340  {0.961753400252887, 0., 0.787405595186356, -2.74544192086633, 3.70351728061223, -1.70723435518514},
341  {-1.76418754019038, 0., -0.774504669155511, 9.64023584441292, -12.544886411271, 5.44334277620397},
342  {0.960350435099165, 0., 0.173858014493155, -6.21422862823726, 8.56612859966376, -3.48610842101883}
343 };
344 
345 #include "NamespaceFooter.H"
346 #endif
#define CH_TIMERS(name)
Definition: CH_Timer.H:133
#define CH_assert(cond)
Definition: CHArray.H:37
bool hasDenseOutput() const
Definition: ARK4DenseOutput.H:63
IMEXOp & getImExOp()
Definition: ARK4DenseOutput.H:126
one dimensional dynamic array
Definition: Vector.H:53
#define CH_START(tpointer)
Definition: CH_Timer.H:145
Real m_dt
Definition: ARK4DenseOutput.H:79
void resetDt(Real a_dt)
Definition: ARK4DenseOutput.H:136
static const Real s_c[s_nStages]
Definition: ARK4DenseOutput.H:68
Rhs m_rhs
Definition: ARK4DenseOutput.H:82
#define CH_TIMER(name, tpointer)
Definition: CH_Timer.H:63
static const Real s_bstar[s_nDenseCoefs][s_nStages]
Definition: ARK4DenseOutput.H:73
static const Real s_aI[s_nStages][s_nStages]
Definition: ARK4DenseOutput.H:70
bool m_isDefined
Definition: ARK4DenseOutput.H:76
bool m_hasDenseOutput
Definition: ARK4DenseOutput.H:78
static const Real s_b[s_nStages]
Definition: ARK4DenseOutput.H:71
Soln m_phi[s_nStages]
Definition: ARK4DenseOutput.H:81
void start0end1(bool a_start0, bool a_end1)
Definition: ARK4DenseOutput.H:155
IMEXOp m_opImEx
Definition: ARK4DenseOutput.H:86
double Real
Definition: REAL.H:33
bool m_denseOutput
Definition: ARK4DenseOutput.H:77
T Abs(const T &a_a)
Definition: Misc.H:53
Real m_time
Definition: ARK4DenseOutput.H:80
Rhs m_kI
Definition: ARK4DenseOutput.H:85
Definition: ARK4DenseOutput.H:34
static const Real s_aE[s_nStages][s_nStages]
Definition: ARK4DenseOutput.H:69
void advance(Real a_time, Soln &a_state)
Definition: ARK4DenseOutput.H:173
#define CH_STOP(tpointer)
Definition: CH_Timer.H:150
static const int s_nDenseCoefs
Definition: ARK4DenseOutput.H:72
size_t size() const
Definition: Vector.H:192
Rhs m_denseCoefs[s_nDenseCoefs]
Definition: ARK4DenseOutput.H:83
void define(const Soln &a_state, Real a_dt, bool a_denseOutput=false)
Definition: ARK4DenseOutput.H:96
void denseOutputCoefs(Vector< Rhs * > &a_interpCoefs)
Definition: ARK4DenseOutput.H:264
bool isDefined() const
Definition: ARK4DenseOutput.H:61
static const int s_nStages
Runge-Kutta coefficients.
Definition: ARK4DenseOutput.H:66
Rhs m_kE
Definition: ARK4DenseOutput.H:84
static const Real s_aIdiag
Definition: ARK4DenseOutput.H:67