Chombo + EB  3.2
RK4DenseOutput.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 _RK4DENSEOUTPUT_H_
12 #define _RK4DENSEOUTPUT_H_
13 
14 
15 #include "NamespaceHeader.H"
16 
17 /// This is a standard RK4 implementation
18 /// that also includes dense output coefficients
19 
20 // Soln is the type that operators are applied to
21 // Rhs is the type that operator produce
22 // EXOp wraps the explicit operator
23 template <class Soln, class Rhs, class EXOP>
25 {
26 public:
27 
29 
30  // This must be called first.
31  void define(const Soln& a_state, Real a_dt, bool a_denseOutput = false);
32 
33  // Advance one step.
34  void advance(Real a_time, Soln& a_state);
35 
36  // Return current dense output coefs, 0th power first, etc.
37  // NOTE: These are intended to be for <JU> conserved quantities
38  void denseOutputCoefs(Vector<Rhs*>& a_interpCoefs);
39 
40  // Reset the timestep.
41  void resetDt(Real a_dt);
42 
43  // Set whether the step starts at 0. and/or ends at 1.
44  void start0end1(bool a_start0, bool a_end1);
45 
46  // Access to the operators if we're defined already
47  // Caller is responsible for making sure they're in a good state
48  // (for example, like have resetDt() called, etc.
49  EXOP& getEXOP();
50 
51  bool isDefined() const { return m_isDefined; }
52 
53  bool hasDenseOutput() const { return m_hasDenseOutput; }
54 
55  /// Runge-Kutta coefficients
56  static const int s_nStages = 4;
57  static const Real s_c[s_nStages];
58  static const Real s_a[s_nStages][s_nStages];
59  static const Real s_b[s_nStages];
60  static const int s_nDenseCoefs = 3;
62 
63 protected:
66  bool m_hasDenseOutput; // If there is dense output data to interpolate
70  Rhs m_rhs;
72  Rhs m_kE;
73  EXOP m_opEx;
74 
75 private:
76 
77 };
78 
79 //==============================================
80 
81 template <class Soln, class Rhs, class EXOP>
83 define(const Soln& a_state, Real a_dt, bool a_denseOutput)
84 {
85  CH_TIMERS("RK4DenseOutput::define");
86  m_dt = a_dt;
87  m_denseOutput = a_denseOutput;
88  // define Soln types
89  for (int stage = 0; stage < s_nStages; stage++)
90  m_phi[stage].define(a_state);
91 
92  // define Rhs types
93  m_kE.define(a_state);
94  m_rhs.define(a_state);
95 
96  // define opEx
97  m_opEx.define(a_state, m_dt);
98 
99  // if dense output is requested, need more storage
100  for (int coef = 0; m_denseOutput && (coef < s_nDenseCoefs); coef++)
101  m_denseCoefs[coef].define(a_state);
102 
103  m_hasDenseOutput = false;
104  m_isDefined = true;
105 }
106 
107 /*
108  Get a reference to the implicit-explicit operator
109  */
110 template <class Soln, class Rhs, class EXOP>
113 {
114  return m_opEx;
115 }
116 
117 /*
118  Reset the timestep.
119  */
120 template <class Soln, class Rhs, class EXOP>
123 {
124  CH_assert(isDefined());
125 
126  // Only update everything if dt has changed
127  Real reltol = 1e-14;
128  if (Abs(m_dt - a_dt) > m_dt*reltol)
129  {
130  m_dt = a_dt;
131  m_hasDenseOutput = false;
132  m_opEx.resetDt(m_dt);
133  }
134 }
135 
136 /*
137  Set whether the step starts at 0. and/or ends at 1.
138  */
139 template <class Soln, class Rhs, class EXOP>
141 start0end1(bool a_start0, bool a_end1)
142 {
143  CH_assert(isDefined());
144 
145  int stage0 = -1; // if no stage is at time 0.
146  if (a_start0) stage0 = 0;
147 
148  int stage1 = -1; // if no stage is at time 1.
149  if (a_end1) stage1 = s_nStages-1;
150 
151  m_opEx.stage0stage1(stage0, stage1);
152 }
153 
154 /*
155  Advance solution a_state in time, a_time to a_time + a_dt.
156  */
157 template <class Soln, class Rhs, class EXOP>
159 advance(Real a_time, Soln& a_state)
160 {
161  CH_TIMERS("RK4DenseOutput::advance");
162  CH_assert(isDefined());
163 
164  CH_TIMER("RK4DenseOutput::advance - explicit op", t1);
165 
166  // Reset the dense output coefs
167  if (m_denseOutput)
168  {
169  m_hasDenseOutput = false;
170  for (int icoef=0; icoef < s_nDenseCoefs; ++icoef)
171  m_denseCoefs[icoef].zero();
172  }
173 
174  // Copy a_state into all stages
175  for (int stage = 0; stage < s_nStages; stage++)
176  m_phi[stage].copy(a_state);
177 
178  // Copy a_state values back to itself, to reset any internal state
179  // (like fine flux registers for this level)
180  a_state.copy(m_phi[0]);
181 
182  // For each stage
183  for (int stage = 0; stage < s_nStages; stage++)
184  {
185  Real t = a_time + s_c[stage]*m_dt;
186  if (stage > 0)
187  {
188  // Copy rhs from phi[stage]
189  m_rhs.copy(m_phi[stage]);
190  }
191 
192  // Calculate the operators for this stage
193  CH_START(t1);
194  m_opEx.explicitOp(m_kE, t, m_phi[stage], stage);
195  CH_STOP(t1);
196  // Accumulate the explicit op into the stage rhs, final solution
197  // NOTE: increment any internal registers
198  a_state.increment(m_kE, m_dt*s_b[stage], true);
199  for (int k=stage+1; k < s_nStages; ++k)
200  m_phi[k].increment(m_kE, m_dt*s_a[k][stage]);
201 
202  // This might provide an optimization
203  // Accumulate the final solution diff from last stage, explicit op only
204  // a_state.increment(m_kE, m_dt*(s_b[stage] - s_aE[s_nStages-1][stage]));
205 
206  if (m_denseOutput)
207  {
208  // pout() << "Stage: " << stage-1 << endl;
209  for (int icoef=0; icoef < s_nDenseCoefs; ++icoef)
210  {
211  m_denseCoefs[icoef].increment(m_kE, m_dt*s_bstar[icoef][stage]);
212 
213  /*
214  LevelData<FArrayBox>& coef = m_denseCoefs[icoef].data();
215  DataIterator dit(coef.getBoxes());
216  for (dit.begin(); dit.ok(); ++dit)
217  pout() << " Coef[" << icoef+1 << "] = " << coef[dit].min() << endl;
218  */
219  }
220  }
221  }
222 
223  m_hasDenseOutput = m_denseOutput;
224  m_time = a_time;
225 }
226 
227 /*
228  Return the coefs to interpolate solution, in terms of power of the fraction
229  of time between t_old and t_new.
230  */
231 template <class Soln, class Rhs, class EXOP>
234 {
235  CH_TIMERS("RK4DenseOutput::denseOutputCoefs");
236 
237  const int nCoef = s_nDenseCoefs+1;
238  CH_assert(m_hasDenseOutput);
239  CH_assert(a_interpCoefs.size() == nCoef);
240 
241  for (int icoef=0; icoef < nCoef; ++icoef)
242  CH_assert(a_interpCoefs[icoef] != NULL);
243 
244  // Copy over the dense coef values
245 
246  // First coeficient is just the old state
247  a_interpCoefs[0]->copy(m_phi[0]);
248 
249  // Next coefs are our dense output
250  for (int icoef = 1; icoef < nCoef ; ++icoef)
251  {
252  /*
253  LevelData<FArrayBox>& coef = m_denseCoefs[icoef].data();
254  for (dit.begin(); dit.ok(); ++dit)
255  pout() << " Coef[" << icoef+1 << "] = " << coef[dit].min() << endl;
256  */
257  a_interpCoefs[icoef]->copy(m_denseCoefs[icoef-1]);
258  }
259 }
260 
261 
262 /*
263  Static constants for RK4
264  */
265 
266 // Time coefficients for each stage
267 template <class Soln, class Rhs, class EXOP>
269 s_c[] = { 0.0, 0.5, 0.5, 1.0 };
270 
271 // Stage coefficients - each row is for that stage
272 template <class Soln, class Rhs, class EXOP>
275  {0., 0., 0., 0.},
276  {0.5, 0., 0., 0.},
277  {0., 0.5, 0., 0.},
278  {0., 0., 1.0, 0.},
279 };
280 
281 // Final stage coefficients
282 template <class Soln, class Rhs, class EXOP>
284 s_b[] =
285  {0.16666666666666667, 0.33333333333333333, 0.33333333333333333, 0.16666666666666667};
286 
287 // Coefficients for dense ouput, 4th-order interpolation
288 template <class Soln, class Rhs, class EXOP>
291  {1.0, 0., 0., 0.},
292  {-1.5, 1.0, 1.0, -0.5},
293  {0.66666666666666667, -0.66666666666666667, -0.66666666666666667, 0.66666666666666667}
294 };
295 
296 
297 #include "NamespaceFooter.H"
298 #endif
#define CH_TIMERS(name)
Definition: CH_Timer.H:133
EXOP & getEXOP()
Definition: RK4DenseOutput.H:112
bool m_isDefined
Definition: RK4DenseOutput.H:64
static const int s_nStages
Runge-Kutta coefficients.
Definition: RK4DenseOutput.H:56
#define CH_assert(cond)
Definition: CHArray.H:37
static const Real s_a[s_nStages][s_nStages]
Definition: RK4DenseOutput.H:58
one dimensional dynamic array
Definition: Vector.H:53
#define CH_START(tpointer)
Definition: CH_Timer.H:145
EXOP m_opEx
Definition: RK4DenseOutput.H:73
static const Real s_c[s_nStages]
Definition: RK4DenseOutput.H:57
#define CH_TIMER(name, tpointer)
Definition: CH_Timer.H:63
Rhs m_rhs
Definition: RK4DenseOutput.H:70
bool isDefined() const
Definition: RK4DenseOutput.H:51
Real m_dt
Definition: RK4DenseOutput.H:67
static const int s_nDenseCoefs
Definition: RK4DenseOutput.H:60
Definition: RK4DenseOutput.H:24
void start0end1(bool a_start0, bool a_end1)
Definition: RK4DenseOutput.H:141
double Real
Definition: REAL.H:33
void resetDt(Real a_dt)
Definition: RK4DenseOutput.H:122
T Abs(const T &a_a)
Definition: Misc.H:53
bool m_hasDenseOutput
Definition: RK4DenseOutput.H:66
void denseOutputCoefs(Vector< Rhs * > &a_interpCoefs)
Definition: RK4DenseOutput.H:233
Real m_time
Definition: RK4DenseOutput.H:68
void define(const Soln &a_state, Real a_dt, bool a_denseOutput=false)
Definition: RK4DenseOutput.H:83
Rhs m_denseCoefs[s_nDenseCoefs]
Definition: RK4DenseOutput.H:71
static const Real s_bstar[s_nDenseCoefs][s_nStages]
Definition: RK4DenseOutput.H:61
#define CH_STOP(tpointer)
Definition: CH_Timer.H:150
bool m_denseOutput
Definition: RK4DenseOutput.H:65
Rhs m_kE
Definition: RK4DenseOutput.H:72
bool hasDenseOutput() const
Definition: RK4DenseOutput.H:53
size_t size() const
Definition: Vector.H:192
Soln m_phi[s_nStages]
Definition: RK4DenseOutput.H:69
void advance(Real a_time, Soln &a_state)
Definition: RK4DenseOutput.H:159
static const Real s_b[s_nStages]
Definition: RK4DenseOutput.H:59