Chombo + EB  3.0
LevelRK4.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_H__
12 #define _LEVELRK4_H__
13 
14 /// Basic templated implementation of RK4 advance for a single AMRLevel
15 /**
16  This implements fourth-order Runge-Kutta time advance of an ODE.
17  ODE is d(Soln)/dt = RHS
18 
19  Template types:
20  \li TSoln is the datatype for the solution
21  \li TFR is a flux-register datatype
22  \li TOp is an object encapsulating the actual ODE operations
23 
24  TOp requirements: TOp must have the following functions:
25 \code
26  // evaluate d(soln)/dt at current time based on soln
27  void evalRHS(TSoln& rhs, // d(soln)/dt based on soln
28  TSoln& soln, // soln at current time
29  TFR& fineFR, // flux register w/ finer level
30  TFR& crseFR, // flux register w/ crse level
31  const TSoln& oldCrseSoln, // old-time crse solution
32  Real oldCrseTime, // old crse time
33  const TSoln& newCrseSoln, // new-time crse solution
34  Real newCrseTime, // new crse time
35  Real time, // current time centering of soln
36  Real fluxWeight // weight to apply to fluxRegister updates
37  )
38 
39  // implements soln += dt*rhs
40  void updateODE(TSoln& soln,
41  const TSoln& rhs,
42  Real dt)
43 
44  // define data holder newSoln based on existingSoln,
45  // including ghost cell specification
46  void defineSolnData(TSoln& newSoln,
47  const TSoln& existingSoln)
48 
49  // define data holder for RHS based on existingSoln
50  // including ghost cell specification
51  // (which in most cases is no ghost cells)
52  void defineRHSData(TSoln& newRHS, const TSoln& existingSoln)
53 
54  /// copy data in src into dest
55  void copySolnData(TSoln& dest, const TSoln& src)
56 \endcode
57 */
58 
59 template <class TSoln, class TFR, class TOp>
60 void RK4LevelAdvance(/// the cell-centered solution at the new time (a_time + a_dt)
61  TSoln& a_newSoln,
62  /// the cell-centered solution at the old time (a_time)
63  TSoln& a_oldSoln,
64  /// old-time coarse-level solution (if it exists)
65  const TSoln& a_oldCrseSoln,
66  /// time-centering of a_oldCrseSoln
67  Real a_oldCrseTime,
68  /// new-time coarse-level solution (if it exists)
69  const TSoln& a_newCrseSoln,
70  /// time-centering of a_newCrseSoln
71  Real a_newCrseTime,
72  /// the pointer to the flux register between this level and the next coarser level
73  TFR& a_crseFRPtr,
74  /// the pointer to the flux register between this level and the next finer level
75  TFR& a_fineFRPtr,
76  /// time centering of a_oldSoln
77  Real a_time,
78  /// time step
79  Real a_dt,
80  /// object which encapsulates the ODE functionality (must meet the requirements detailed above)
81  TOp& a_op,
82  /// whether a_newSoln is initialized with a_oldSoln (default true)
83  const bool a_initializeNewSoln = true)
84 {
85  // Advance by one time step using 4th-order Runge-Kutta.
86 
87  // allocate and define temp storage
88  TSoln RHSTmp;
89  TSoln solnTmp;
90  a_op.defineRHSData(RHSTmp,a_oldSoln);
91  a_op.defineSolnData(solnTmp,a_oldSoln);
92  a_op.copySolnData(solnTmp,a_oldSoln);
93 
94  // start by copying old solution into new solution
95  if (a_initializeNewSoln)
96  {
97  a_op.copySolnData(a_newSoln, a_oldSoln);
98  }
99 
100  // RK Stage 1: compute RHS.at old time
101  // RHSTmp gets RHS(tOld)
102  // flux registers are also incremented if appropriate
103  a_op.evalRHS(RHSTmp,a_oldSoln,
104  a_fineFRPtr,a_crseFRPtr,
105  a_oldCrseSoln, a_oldCrseTime,
106  a_newCrseSoln, a_newCrseTime,
107  a_time,a_dt/6);
108 
109  // RK Stage 1: compute update {phi}_1, partial update to new vals.
110 
111  // first part of update...
112  a_op.updateODE(a_newSoln,RHSTmp,a_dt/6);
113 
114 
115  // predict half-time value using LofPhiTmp
116  a_op.updateODE(solnTmp, RHSTmp,a_dt/2);
117 
118  // RK Stage 2: compute RHS based on first predicted value
119  // also increment flux registers
120  a_op.evalRHS(RHSTmp,solnTmp,
121  a_fineFRPtr,a_crseFRPtr,
122  a_oldCrseSoln, a_oldCrseTime,
123  a_newCrseSoln, a_newCrseTime,
124  a_time+a_dt/2,a_dt/3);
125 
126 
127  // RK Stage 2: compute update {phi}_2, partial update to new vals.
128 
129  a_op.updateODE(a_newSoln,RHSTmp,a_dt/3);
130 
131  // now predict half-time value using latest estimate of RHS
132  a_op.copySolnData(solnTmp,a_oldSoln);
133  a_op.updateODE(solnTmp,RHSTmp,a_dt/2);
134 
135  // RK Stage 3: compute new RHS.
136 
137  a_op.evalRHS(RHSTmp,solnTmp,
138  a_fineFRPtr,a_crseFRPtr,
139  a_oldCrseSoln, a_oldCrseTime,
140  a_newCrseSoln, a_newCrseTime,
141  a_time+a_dt/2,a_dt/3);
142 
143 
144  // RK Stage 3: compute update {phi}_3, partial update to new vals.
145 
146  a_op.updateODE(a_newSoln,RHSTmp,a_dt/3);
147 
148  // predict value at new time
149  a_op.copySolnData(solnTmp,a_oldSoln);
150  a_op.updateODE(solnTmp,RHSTmp,a_dt);
151 
152  // RK Stage 4: compute RHS.
153 
154  a_op.evalRHS(RHSTmp,solnTmp,
155  a_fineFRPtr,a_crseFRPtr,
156  a_oldCrseSoln, a_oldCrseTime,
157  a_newCrseSoln, a_newCrseTime,
158  a_time+a_dt,a_dt/6);
159 
160  // RK Stage 4: compute final update of solution.
161 
162  a_op.updateODE(a_newSoln,RHSTmp,a_dt/6);
163 }
164 
165 
166 /// Templated implementation of RK4 advance for a single AMRLevel, allowing interpolation in time
167 /**
168  This implements fourth-order Runge-Kutta time advance of an ODE,
169  and saves intermediate results for use in interpolation in time.
170  ODE is d(Soln)/dt = RHS
171 
172  Template types:
173  \li TSoln is the datatype for the solution
174  \li TInterp is a storage class for intermediate values
175  \li TFR is a flux-register datatype
176  \li TOp is an object encapsulating the actual ODE operations
177 
178  TInterp requirements: TInterp must have the following functions:
179 \code
180  // set time step
181  void setDt(Real dt)
182 
183  // save initial solution
184  void saveInitialSoln(TSoln& soln)
185 
186  // save this value of d(soln)/dt
187  void saveRHS(TSoln& rhs)
188 \endcode
189  TOp requirements: TOp must have the following functions:
190 \code
191  // evaluate d(soln)/dt at current time based on soln
192  void evalRHS(TSoln& rhs, // d(soln)/dt based on soln
193  TSoln& soln, // soln at current time
194  TFR& fineFR, // flux register w/ finer level
195  TFR& crseFR, // flux register w/ crse level
196  const TSoln& oldCrseSoln, // old-time crse solution
197  Real oldCrseTime, // old crse time
198  const TSoln& newCrseSoln, // new-time crse solution
199  Real newCrseTime, // new crse time
200  Real time, // current time centering of soln
201  Real fluxWeight // weight to apply to fluxRegister updates
202  )
203 
204  // implements soln += dt*rhs
205  void updateODE(TSoln& soln,
206  const TSoln& rhs,
207  Real dt)
208 
209  // define data holder newSoln based on existingSoln,
210  // including ghost cell specification
211  void defineSolnData(TSoln& newSoln,
212  const TSoln& existingSoln)
213 
214  // define data holder for RHS based on existingSoln
215  // including ghost cell specification
216  // (which in most cases is no ghost cells)
217  void defineRHSData(TSoln& newRHS, const TSoln& existingSoln)
218 
219  /// copy data in src into dest
220  void copySolnData(TSoln& dest, const TSoln& src)
221 \endcode
222 */
223 
224 template <class TSoln, class TInterp, class TFR, class TOp>
225 void RK4LevelAdvance(/// the cell-centered solution at the new time (a_time + a_dt)
226  TSoln& a_newSoln,
227  /// the cell-centered solution at the old time (a_time)
228  TSoln& a_oldSoln,
229  /// object that encapsulates time interpolation
230  TInterp& a_interpPtr,
231  /// old-time coarse-level solution (if it exists)
232  const TSoln& a_oldCrseSoln,
233  /// time-centering of a_oldCrseSoln
234  Real a_oldCrseTime,
235  /// new-time coarse-level solution (if it exists)
236  const TSoln& a_newCrseSoln,
237  /// time-centering of a_newCrseSoln
238  Real a_newCrseTime,
239  /// the pointer to the flux register between this level and the next coarser level
240  TFR& a_crseFRPtr,
241  /// the pointer to the flux register between this level and the next finer level
242  TFR& a_fineFRPtr,
243  /// time centering of a_oldSoln
244  Real a_time,
245  /// time step
246  Real a_dt,
247  /// object which encapsulates the ODE functionality (must meet the requirements detailed above)
248  TOp& a_op,
249  /// whether a_newSoln is initialized with a_oldSoln (default true)
250  const bool a_initializeNewSoln = true)
251 
252 
253 {
254 
255  // Advance by one time step using 4th-order Runge-Kutta.
256 
257  // allocate and define temp storage
258  TSoln RHSTmp;
259  TSoln solnTmp;
260  a_op.defineRHSData(RHSTmp,a_oldSoln);
261  a_op.defineSolnData(solnTmp,a_oldSoln);
262  a_op.copySolnData(solnTmp,a_oldSoln);
263 
264  // start by copying old solution into new solution
265  if (a_initializeNewSoln)
266  {
267  a_op.copySolnData(a_newSoln, a_oldSoln);
268  }
269  a_interpPtr.setDt(a_dt);
270  a_interpPtr.saveInitialSoln(a_oldSoln);
271 
272  // RK Stage 1: compute RHS.at old time
273  // RHSTmp gets RHS(tOld)
274  // flux registers are also incremented if appropriate
275  a_op.evalRHS(RHSTmp,a_oldSoln,
276  a_fineFRPtr,a_crseFRPtr,
277  a_oldCrseSoln, a_oldCrseTime,
278  a_newCrseSoln, a_newCrseTime,
279  a_time,a_dt/6);
280  a_interpPtr.saveRHS(RHSTmp);
281 
282  // RK Stage 1: compute update {phi}_1, partial update to new vals.
283 
284  // first part of update...
285  a_op.updateODE(a_newSoln,RHSTmp,a_dt/6);
286 
287 
288  // predict half-time value using LofPhiTmp
289  a_op.updateODE(solnTmp, RHSTmp,a_dt/2);
290 
291  // RK Stage 2: compute RHS based on first predicted value
292  // also increment flux registers
293  a_op.evalRHS(RHSTmp,solnTmp,
294  a_fineFRPtr,a_crseFRPtr,
295  a_oldCrseSoln, a_oldCrseTime,
296  a_newCrseSoln, a_newCrseTime,
297  a_time+a_dt/2,a_dt/3);
298  a_interpPtr.saveRHS(RHSTmp);
299 
300 
301  // RK Stage 2: compute update {phi}_2, partial update to new vals.
302  a_op.updateODE(a_newSoln,RHSTmp,a_dt/3);
303 
304  // now predict half-time value using latest estimate of RHS
305  a_op.copySolnData(solnTmp,a_oldSoln);
306  a_op.updateODE(solnTmp,RHSTmp,a_dt/2);
307 
308  // RK Stage 3: compute new RHS.
309 
310  a_op.evalRHS(RHSTmp,solnTmp,
311  a_fineFRPtr,a_crseFRPtr,
312  a_oldCrseSoln, a_oldCrseTime,
313  a_newCrseSoln, a_newCrseTime,
314  a_time+a_dt/2,a_dt/3);
315  a_interpPtr.saveRHS(RHSTmp);
316 
317 
318  // RK Stage 3: compute update {phi}_3, partial update to new vals.
319 
320  a_op.updateODE(a_newSoln,RHSTmp,a_dt/3);
321 
322  // predict value at new time
323  a_op.copySolnData(solnTmp,a_oldSoln);
324  a_op.updateODE(solnTmp,RHSTmp,a_dt);
325 
326  // RK Stage 4: compute RHS.
327 
328  a_op.evalRHS(RHSTmp,solnTmp,
329  a_fineFRPtr,a_crseFRPtr,
330  a_oldCrseSoln, a_oldCrseTime,
331  a_newCrseSoln, a_newCrseTime,
332  a_time+a_dt,a_dt/6);
333  a_interpPtr.saveRHS(RHSTmp);
334 
335  // RK Stage 4: compute final update of solution.
336 
337  a_op.updateODE(a_newSoln,RHSTmp,a_dt/6);
338 }
339 
340 #endif // end multiple-include preventer
double Real
Definition: REAL.H:33
void RK4LevelAdvance(TSoln &a_newSoln, TSoln &a_oldSoln, const TSoln &a_oldCrseSoln, Real a_oldCrseTime, const TSoln &a_newCrseSoln, Real a_newCrseTime, TFR &a_crseFRPtr, TFR &a_fineFRPtr, Real a_time, Real a_dt, TOp &a_op, const bool a_initializeNewSoln=true)
Basic templated implementation of RK4 advance for a single AMRLevel.
Definition: LevelRK4.H:60