Chombo + EB + MF  3.2
LevelRK4_v2.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 _LEVELRK4_V2_H__
12 #define _LEVELRK4_V2_H__
13 
14 /// Templated implementation of RK4 advance for a single AMRLevel, allowing interpolation in time
15 /**
16  The difference from LevelRK4.H is that this allows for any arguments to
17  the evalRHS operator. The arguments are perfectly forwarded.
18 
19  This implements fourth-order Runge-Kutta time advance of an ODE,
20  and saves intermediate results for use in interpolation in time.
21  ODE is d(Soln)/dt = RHS
22 
23  Template types:
24  \li TSoln is the datatype for the solution
25  \li TInterp is a storage class for intermediate values
26  \li TFR is a flux-register datatype
27  \li TExOp is an object encapsulating the actual ODE operations
28  \li TOpArgs are perfectly forwarded arguments to the operator
29 
30  TInterp requirements: TInterp must have the following functions:
31 \code
32  // set time step
33  void setDt(Real dt)
34 
35  // save initial solution
36  void saveInitialSoln(TSoln& soln)
37 
38  // save this value of d(soln)/dt
39  void saveRHS(TSoln& rhs)
40 \endcode
41  TOp requirements: TOp must have the following functions:
42 \code
43  // evaluate d(soln)/dt at current time based on soln
44  void evalRHS(TSoln& rhs, // d(soln)/dt based on soln
45  TSoln& soln, // soln at current time
46  TFR& fineFR, // flux register w/ finer level
47  TFR& crseFR, // flux register w/ crse level
48  const TSoln& oldCrseSoln, // old-time crse solution
49  Real oldCrseTime, // old crse time
50  const TSoln& newCrseSoln, // new-time crse solution
51  Real newCrseTime, // new crse time
52  Real time, // current time centering of soln
53  Real fluxWeight // weight to apply to fluxRegister updates
54  )
55 
56  // implements soln += dt*rhs
57  void updateODE(TSoln& soln,
58  const TSoln& rhs,
59  Real dt)
60 
61  // define data holder newSoln based on existingSoln,
62  // including ghost cell specification
63  void defineSolnData(TSoln& newSoln,
64  const TSoln& existingSoln)
65 
66  // define data holder for RHS based on existingSoln
67  // including ghost cell specification
68  // (which in most cases is no ghost cells)
69  void defineRHSData(TSoln& newRHS, const TSoln& existingSoln)
70 
71  /// copy data in src into dest
72  void copySolnData(TSoln& dest, const TSoln& src)
73 \endcode
74 */
75 
76 template <typename TSoln, typename TRhs, typename TInterp, typename TExOp, typename... TOpArgs>
77 void RK4LevelAdvance(/// the cell-centered solution at the new time (a_time + a_dt)
78  TSoln& a_newSoln,
79  /// the cell-centered solution at the old time (a_time)
80  TSoln& a_oldSoln,
81  /// object that encapsulates time interpolation
82  TInterp* a_timeInterpolator,
83  /// time centering of a_oldSoln
84  Real a_time,
85  /// time step
86  Real a_dt,
87  /// whether a_newSoln is initialized with a_oldSoln (default true)
88  const bool a_initializeNewSoln,
89  /// object which encapsulates the ODE functionality (must meet the requirements detailed above)
90  TExOp& a_op,
91  TOpArgs&&... a_opArgs)
92 {
93 
94  // Advance by one time step using 4th-order Runge-Kutta.
95 
96  // allocate and define temp storage
97  TRhs RHSTmp;
98  TSoln solnTmp;
99  a_op.defineRHSData(RHSTmp,a_oldSoln);
100  a_op.defineSolnData(solnTmp,a_oldSoln);
101  a_op.copySolnData(solnTmp,a_oldSoln);
102 
103  // start by copying old solution into new solution
104  if (a_initializeNewSoln)
105  {
106  a_op.copySolnData(a_newSoln, a_oldSoln);
107  }
108  if (a_timeInterpolator != nullptr)
109  {
110  a_timeInterpolator->setDt(a_dt);
111  a_timeInterpolator->saveInitialSoln(a_oldSoln);
112  }
113 
114  // RK Stage 1: compute RHS.at old time
115  // RHSTmp gets RHS(tOld)
116  // flux registers are also incremented if appropriate
117  a_op.evalRHS(RHSTmp,a_oldSoln,0,a_time,a_dt/6,std::forward<TOpArgs>(a_opArgs)...);
118  if (a_timeInterpolator != nullptr)
119  {
120  a_timeInterpolator->saveRHS(RHSTmp);
121  }
122 
123  // RK Stage 1: compute update {phi}_1, partial update to new vals.
124 
125  // first part of update...
126  a_op.updateODE(a_newSoln,RHSTmp,a_dt/6);
127 
128 
129  // predict half-time value using LofPhiTmp
130  a_op.updateODE(solnTmp, RHSTmp,a_dt/2);
131 
132  // RK Stage 2: compute RHS based on first predicted value
133  // also increment flux registers
134  a_op.evalRHS(RHSTmp,solnTmp,1,a_time+a_dt/2,a_dt/3,std::forward<TOpArgs>(a_opArgs)...);
135  if (a_timeInterpolator != nullptr)
136  {
137  a_timeInterpolator->saveRHS(RHSTmp);
138  }
139 
140 
141  // RK Stage 2: compute update {phi}_2, partial update to new vals.
142  a_op.updateODE(a_newSoln,RHSTmp,a_dt/3);
143 
144  // now predict half-time value using latest estimate of RHS
145  a_op.copySolnData(solnTmp,a_oldSoln);
146  a_op.updateODE(solnTmp,RHSTmp,a_dt/2);
147 
148  // RK Stage 3: compute new RHS.
149 
150  a_op.evalRHS(RHSTmp,solnTmp,2,a_time+a_dt/2,a_dt/3,std::forward<TOpArgs>(a_opArgs)...);
151  if (a_timeInterpolator != nullptr)
152  {
153  a_timeInterpolator->saveRHS(RHSTmp);
154  }
155 
156 
157  // RK Stage 3: compute update {phi}_3, partial update to new vals.
158 
159  a_op.updateODE(a_newSoln,RHSTmp,a_dt/3);
160 
161  // predict value at new time
162  a_op.copySolnData(solnTmp,a_oldSoln);
163  a_op.updateODE(solnTmp,RHSTmp,a_dt);
164 
165  // RK Stage 4: compute RHS.
166 
167  a_op.evalRHS(RHSTmp,solnTmp,3,a_time+a_dt,a_dt/6,std::forward<TOpArgs>(a_opArgs)...);
168  if (a_timeInterpolator != nullptr)
169  {
170  a_timeInterpolator->saveRHS(RHSTmp);
171  }
172 
173  // RK Stage 4: compute final update of solution.
174 
175  a_op.updateODE(a_newSoln,RHSTmp,a_dt/6);
176 }
177 
178 #endif /* ! defined _LEVELRK4_V2_H_ */
void RK4LevelAdvance(TSoln &a_newSoln, TSoln &a_oldSoln, TInterp *a_timeInterpolator, Real a_time, Real a_dt, const bool a_initializeNewSoln, TExOp &a_op, TOpArgs &&... a_opArgs)
Templated implementation of RK4 advance for a single AMRLevel, allowing interpolation in time...
Definition: LevelRK4_v2.H:77
double Real
Definition: REAL.H:33