explicit_timesteppers.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Non-inline member functions for the explicit timesteppers
31 #include "explicit_timesteppers.h"
32 #include "timesteppers.h"
33 
34 namespace oomph
35 {
36  //===================================================================
37  ///Dummy value of time always set to zero
38  //==================================================================
40 
41  //==================================================================
42  ///\short A single virtual function that returns the residuals
43  ///vector multiplied by the inverse mass matrix
44  //=================================================================
47  {
48  std::ostringstream error_stream;
49  error_stream
50  << "Empty default function called.\n"
51  << "The function must return the solution x of the linear system\n"
52  << " M x = R\n"
53  << "in order for the object to be used by an ExplicitTimeStepper.\n"
54  << "NOTE: It is the responsibility of the object to set the size \n"
55  << " of the vector x\n";
56 
57 
58  throw OomphLibError(error_stream.str(),
59  OOMPH_CURRENT_FUNCTION,
60  OOMPH_EXCEPTION_LOCATION);
61  }
62 
63  //=======================================================================
64  /// Function that should get the values of the dofs in the object
65  //=======================================================================
67  {
68  std::ostringstream error_stream;
69  error_stream
70  << "Empty default function called.\n"
71  << "The function must return the current values of the degrees of \n"
72  << "freedom in the object.\n"
73  << "Note: It is the responsibility of the object to set the size\n"
74  << "of the vector\n";
75 
76  throw OomphLibError(error_stream.str(),
77  OOMPH_CURRENT_FUNCTION,
78  OOMPH_EXCEPTION_LOCATION);
79 
80  }
81 
82  //=======================================================================
83  /// Function that should get the values of the dofs in the object
84  //=======================================================================
86  DoubleVector &dofs) const
87  {
88  std::ostringstream error_stream;
89  error_stream
90  << "Empty default function called.\n"
91  << "The function must return the t'th history values of the degrees of \n"
92  << "freedom in the object.\n"
93  << "Note: It is the responsibility of the object to set the size\n"
94  << "of the vector\n";
95 
96  throw OomphLibError(error_stream.str(),
97  OOMPH_CURRENT_FUNCTION,
98  OOMPH_EXCEPTION_LOCATION);
99 
100  }
101 
102  //=====================================================================
103  /// Function that sets the values of the dofs in the object
104  //====================================================================
106  {
107  std::ostringstream error_stream;
108  error_stream
109  << "Empty default function called.\n"
110  << "The function must set the current values of the degrees of \n"
111  << "freedom in the object.\n";
112 
113  throw OomphLibError(error_stream.str(),
114  OOMPH_CURRENT_FUNCTION,
115  OOMPH_EXCEPTION_LOCATION);
116 
117  }
118 
119 
120  //====================================================================
121  ///Function that adds the lambda multiplied by the increment_dofs
122  ///vector to the dofs in the object
123  //====================================================================
125  const double &lambda,
126  const DoubleVector &increment_dofs)
127  {
128  std::ostringstream error_stream;
129  error_stream
130  << "Empty default function called.\n"
131  << "The function must add lambda multiplied by the contents of the\n"
132  << "input vector to the degrees of freedom in the object.\n"
133  << "Note: It is the responsibility of the object to ensure that the\n"
134  << " the degrees of freedom are in the same order as those \n"
135  << " returned by get_dvaluesdt()\n";
136 
137  throw OomphLibError(error_stream.str(),
138  OOMPH_CURRENT_FUNCTION,
139  OOMPH_EXCEPTION_LOCATION);
140 
141  }
142 
143  //==================================================================
144  ///\short Virtual function that should be overloaded to return access
145  ///to the local time in the object
146  //=================================================================
148  {
149  std::ostringstream error_stream;
150  error_stream
151  << "Empty default function called.\n"
152  << "The function must return a reference to the local time in the object\n";
153 
154  throw OomphLibError(error_stream.str(),
155  OOMPH_CURRENT_FUNCTION,
156  OOMPH_EXCEPTION_LOCATION);
157 
158  return Dummy_time_value;
159  }
160 
161  /// \short Virtual function that should be overloaded to return a pointer to a
162  /// Time object.
164  {
165  std::ostringstream error_stream;
166  error_stream
167  << "Empty default function called.\n"
168  << "The function must return a pointer to an oomph-lib Time object.\n";
169  throw OomphLibError(error_stream.str(),
170  OOMPH_CURRENT_FUNCTION,
171  OOMPH_EXCEPTION_LOCATION);
172 
173  return 0;
174  }
175 
176 
177  //================================================================
178  /// Euler timestepping x^{t+1} = x^{t} + dt M^{-1} L(x^{t})
179  //=================================================================
181  const double &dt)
182  {
184  object_pt->actions_before_explicit_stage();
185 
186  //Vector that will hold the inverse mass matrix multiplied by the
187  //residuals
188  DoubleVector minv_res;
189  //Get M^{-1} R
190  object_pt->get_dvaluesdt(minv_res);
191 
192  //Add the result to the unknowns
193  object_pt->add_to_dofs(dt,minv_res);
194  //Increment the time by the appropriate amount
195  object_pt->time() += dt;
196 
197  //Call any actions required after the change in the unknowns
198  object_pt->actions_after_explicit_stage();
199  object_pt->actions_after_explicit_timestep();
200  }
201 
202 
203 //====================================================================
204 //Broken default timestep function for RungeKutta schemes
205 //====================================================================
206 template<unsigned ORDER>
208  const double &dt)
209 {
210  std::ostringstream error_stream;
211  error_stream << "Timestep not implemented for order " << ORDER << "\n";
212  throw OomphLibError(error_stream.str(),
213  OOMPH_CURRENT_FUNCTION,
214  OOMPH_EXCEPTION_LOCATION);
215 }
216 
217 //===================================================================
218 ///Explicit specialisation for fourth-order RK scheme
219 //==================================================================
220 template<>
222  const double &dt)
223 {
225 
226  // Stage 1
227  // ============================================================
228  object_pt->actions_before_explicit_stage();
229 
230  //Store the initial values and initial time
231  DoubleVector u;
232  object_pt->get_dofs(u);
233 
234  //Now get the first unknowns
235  DoubleVector k1;
236  object_pt->get_dvaluesdt(k1);
237 
238  //Add to the residuals
239  object_pt->add_to_dofs(0.5*dt,k1);
240  //Increment the time
241  object_pt->time() += 0.5*dt;
242  object_pt->actions_after_explicit_stage();
243 
244 
245  // Stage 2
246  // ============================================================
247  object_pt->actions_before_explicit_stage();
248 
249  //Get the next unknowns
250  DoubleVector k2;
251  object_pt->get_dvaluesdt(k2);
252 
253  //Now reset the residuals
254  object_pt->set_dofs(u);
255  object_pt->add_to_dofs(0.5*dt,k2);
256  //Time remains the same
257  object_pt->actions_after_explicit_stage();
258 
259  // Stage 3
260  // ============================================================
261  object_pt->actions_before_explicit_stage();
262 
263  //Get the next unknowns
264  DoubleVector k3;
265  object_pt->get_dvaluesdt(k3);
266 
267  //Now reset the residuals
268  object_pt->set_dofs(u);
269  object_pt->add_to_dofs(dt,k3);
270  //Increment the time (now at initial_time + dt)
271  object_pt->time() += 0.5*dt;
272  object_pt->actions_after_explicit_stage();
273 
274  // Stage 4
275  // ============================================================
276  object_pt->actions_before_explicit_stage();
277 
278  //Get the final unknowns
279  DoubleVector k4;
280  object_pt->get_dvaluesdt(k4);
281 
282  //Set the final values of the unknowns
283  object_pt->set_dofs(u);
284  object_pt->add_to_dofs((dt/6.0),k1);
285  object_pt->add_to_dofs((dt/3.0),k2);
286  object_pt->add_to_dofs((dt/3.0),k3);
287  object_pt->add_to_dofs((dt/6.0),k4);
288  object_pt->actions_after_explicit_stage();
289 
290  object_pt->actions_after_explicit_timestep();
291 }
292 
293  //===================================================================
294  ///Explicit specialisation for second-order RK scheme
295  //==================================================================
296  template<>
298  const double &dt)
299  {
301 
302  // Stage 1
303  // ============================================================
304  object_pt->actions_before_explicit_stage();
305 
306  // Store the initial values
307  DoubleVector u;
308  object_pt->get_dofs(u);
309 
310  // Get f1 (time derivative at t0, y0) and add to dofs
311  DoubleVector f1;
312  object_pt->get_dvaluesdt(f1);
313  object_pt->add_to_dofs(dt, f1);
314 
315  // Advance time to t1 = t0 + dt
316  object_pt->time() += dt;
317 
318  object_pt->actions_after_explicit_stage();
319 
320 
321  // Stage 2
322  // ============================================================
323  object_pt->actions_before_explicit_stage();
324 
325  // get f2 (with t=t1, y = y0 + h f0)
326  DoubleVector f2;
327  object_pt->get_dvaluesdt(f2);
328 
329  // Final answer is starting dofs + h/2 * (f1 + f2)
330  object_pt->set_dofs(u);
331  object_pt->add_to_dofs(0.5*dt, f1);
332  object_pt->add_to_dofs(0.5*dt, f2);
333 
334  object_pt->actions_after_explicit_stage();
335 
336  // Done, do actions
337  object_pt->actions_after_explicit_timestep();
338  }
339 
340 
341 //=================================================================
342 //General constructor for LowOrder RK schemes
343 //==================================================================
344 template<unsigned ORDER>
346 {
347  Type = "LowStorageRungeKutta";
348 }
349 
350 //======================================================================
351 ///Broken default timestep for LowStorageRungeKutta
352 //======================================================================
353 template<unsigned ORDER>
355  ExplicitTimeSteppableObject* const &object_pt,
356  const double &dt)
357 {
358  std::ostringstream error_stream;
359  error_stream << "Timestep not implemented for order " << ORDER << "\n";
360  throw OomphLibError(error_stream.str(),
361  OOMPH_CURRENT_FUNCTION,
362  OOMPH_EXCEPTION_LOCATION);
363 }
364 
365 //================================================================
366 //Specialised constructor for fourth-order RK scheme
367 //================================================================
368 template<>
370 {
371  Type="LowStorageRungeKutta";
372 
373  A.resize(5);
374  A[0] = 0.0;
375  A[1] = - 567301805773.0/1357537059087.0;
376  A[2] = - 2404267990393.0/2016746695238.0;
377  A[3] = - 3550918686646.0/2091501179385.0;
378  A[4] = - 1275806237668.0/842570457699.0;
379 
380  B.resize(5);
381  B[0] = 1432997174477.0/9575080441755.0;
382  B[1] = 5161836677717.0/13612068292357.0;
383  B[2] = 1720146321549.0/2090206949498.0;
384  B[3] = 3134564353537.0/4481467310338.0;
385  B[4] = 2277821191437.0/14882151754819.0;
386 
387  C.resize(5);
388  C[0] = B[0];
389  C[1] = 2526269341429.0/6820363962896.0;
390  C[2] = 2006345519317.0/3224310063776.0;
391  C[3] = 2802321613138.0/2924317926251.0;
392  C[4] = 1.0;
393 }
394 
395 
396 
397 //Explicit specialisation for fourth-order RK scheme
398 template<>
400  ExplicitTimeSteppableObject* const &object_pt,
401  const double &dt)
402 {
404 
405  //Store the initial time
406  const double initial_time = object_pt->time();
407  //Temporary storage
408  DoubleVector k;
409 
410  //Temporary storage for the inverse mass matrix multiplied by the residuals
411  DoubleVector minv_res;
412 
413  //Loop over the number of steps in the scheme
414  for(unsigned i=0;i<5;i++)
415  {
416  object_pt->actions_before_explicit_stage();
417 
418  //Get the inverse mass matrix multiplied by the current value
419  //of the residuals
420  object_pt->get_dvaluesdt(minv_res);
421  //Get the values of k
422  const unsigned n_dof= minv_res.nrow();
423 
424  //First time round resize k and initialise to zero
425  if(i==0) {k.build(minv_res.distribution_pt(),0.0);}
426  //Now construct the next value of k
427  for(unsigned n=0;n<n_dof;n++)
428  {
429  k[n] *= A[i];
430  k[n] += dt*minv_res[n];
431  }
432 
433  //Add to the residuals
434  object_pt->add_to_dofs(B[i],k);
435  //Set the new time
436  object_pt->time() = initial_time + C[i]*dt;
437 
438  object_pt->actions_after_explicit_stage();
439  }
440 
441  object_pt->actions_after_explicit_timestep();
442 }
443 
444 // ??ds this could be heavily optimised if needed. Keeping it simple for
445 // now
447  const double &dt)
448 {
449  using namespace StringConversion;
450 
452  object_pt->actions_before_explicit_stage();
453 
454  // Storage indicies for the history values that we need
455  unsigned tn = 1;
456  unsigned tnm1 = tn+1;
457  unsigned tnm2 = tnm1+1;
458 
459  // Check dts are the same, this will need to be removed if ebdf3 is being
460  // used as something other than a predictor... But seeing as it isn't
461  // stable that isn't likely.
462 #ifdef PARANOID
463  if(std::abs(dt - object_pt->time_pt()->dt(0)) > 1e-15)
464  {
465  std::string err = "Inconsistent dts! Predictor is stepping by " + to_string(dt);
466  err += " but dt(0) = " + to_string(object_pt->time_pt()->dt(0));
467  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
468  OOMPH_CURRENT_FUNCTION);
469  }
470 #endif
471 
472  // Get older dt values
473  double dtnm1 = object_pt->time_pt()->dt(1);
474  double dtnm2 = object_pt->time_pt()->dt(2);
475 
476  // Calculate weights for these dts
477  set_weights(dt, dtnm1, dtnm2);
478 
479  // Get derivative value at step n (even though this uses values from t=0 =
480  // step n+1, it's ok because we haven't changed the values in that slot
481  // yet).
482  DoubleVector fn;
483  object_pt->get_dvaluesdt(fn);
484  fn *= Fn_weight;
485 
486  // Extract history values and multiply by their weights
487  DoubleVector ynp1, yn, ynm1, ynm2;
488  object_pt->get_dofs(tn, yn);
489  yn *= Yn_weight;
490  object_pt->get_dofs(tnm1, ynm1);
491  ynm1 *= Ynm1_weight;
492  object_pt->get_dofs(tnm2, ynm2);
493  ynm2 *= Ynm2_weight;
494 
495 
496  // Add everything together
497  ynp1 = yn;
498  ynp1 += ynm1;
499  ynp1 += ynm2;
500  ynp1 += fn;
501 
502  // Done, update things in the object
503  object_pt->set_dofs(ynp1);
504  object_pt->time() += dt;
505 
506  // Actions functions
507  object_pt->actions_after_explicit_stage();
508  object_pt->actions_after_explicit_timestep();
509 }
510 
511 /// Calculate the weights for this set of step sizes.
512 void EBDF3::set_weights(const double &dtn, const double &dtnm1,
513  const double &dtnm2)
514  {
515  using namespace std;
516 
517  // If this is slow we can probably optimise by doing direct
518  // multiplication instead of using pow.
519 
520  // Generated using sympy from my python ode code.
521 
522  double denom = pow(dtnm1,4)*dtnm2 + 2*pow(dtnm1,3)*pow(dtnm2,2)
523  + pow(dtnm1,2)*pow(dtnm2,3);
524 
525  Yn_weight = -(2*pow(dtn,3)*dtnm1*dtnm2 + pow(dtn,3)*pow(dtnm2,2)
526  + 3*pow(dtn,2)*pow(dtnm1,2)*dtnm2
527  + 3*pow(dtn,2)*dtnm1*pow(dtnm2,2) + pow(dtn,2)*pow(dtnm2,3)
528  - pow(dtnm1,4)*dtnm2 - 2*pow(dtnm1,3)*pow(dtnm2,2)
529  - pow(dtnm1,2)*pow(dtnm2,3))/denom;
530 
531  Ynm1_weight = -(-pow(dtn,3)*pow(dtnm1,2) - 2*pow(dtn,3)*dtnm1*dtnm2
532  - pow(dtn,3)*pow(dtnm2,2) - pow(dtn,2)*pow(dtnm1,3)
533  - 3*pow(dtn,2)*pow(dtnm1,2)*dtnm2
534  - 3*pow(dtn,2)*dtnm1*pow(dtnm2,2)
535  - pow(dtn,2)*pow(dtnm2,3))/denom;
536 
537  Ynm2_weight = -(pow(dtn,3)*pow(dtnm1,2) + pow(dtn,2)*pow(dtnm1,3))/denom;
538 
539  Fn_weight = -(-pow(dtn,3)*pow(dtnm1,2)*dtnm2 - pow(dtn,3)*dtnm1*pow(dtnm2,2)
540  - 2*pow(dtn,2)*pow(dtnm1,3)*dtnm2
541  - 3*pow(dtn,2)*pow(dtnm1,2)*pow(dtnm2,2)
542  - pow(dtn,2)*dtnm1*pow(dtnm2,3) - dtn*pow(dtnm1,4)*dtnm2
543  - 2*dtn*pow(dtnm1,3)*pow(dtnm2,2)
544  - dtn*pow(dtnm1,2)*pow(dtnm2,3))/denom;
545  }
546 
547 
548 
549 //Force build of templates
550 template class RungeKutta<4>;
551 template class LowStorageRungeKutta<4>;
552 }
void timestep(ExplicitTimeSteppableObject *const &object_pt, const double &dt)
Function that is used to advance time in the object.
LowStorageRungeKutta()
Constructor, set the type.
virtual void actions_after_explicit_timestep()
Empty virtual function that can be overloaded to do anything needed after an explicit step...
void timestep(ExplicitTimeSteppableObject *const &object_pt, const double &dt)
Function that is used to advance the solution by time dt.
cstr elem_len * i
Definition: cfortran.h:607
virtual double & time()
Broken virtual function that should be overloaded to return access to the local time in the object...
Class to keep track of discrete/continous time. It is essential to have a single Time object when usi...
Definition: timesteppers.h:67
virtual void set_dofs(const DoubleVector &dofs)
Function that sets the values of the dofs in the object.
unsigned nrow() const
access function to the number of global rows.
char t
Definition: cfortran.h:572
virtual void actions_after_explicit_stage()
Empty virtual function that should be overloaded to update and slaved data or boundary conditions tha...
virtual Time * time_pt() const
Virtual function that should be overloaded to return a pointer to a Time object.
e
Definition: cfortran.h:575
virtual void actions_before_explicit_stage()
Empty virtual function to do anything needed before a stage of an explicit time step (Runge-Kutta ste...
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void set_weights(const double &dtn, const double &dtnm1, const double &dtnm2)
Calculate the weights for this set of step sizes.
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void timestep(ExplicitTimeSteppableObject *const &object_pt, const double &dt)
Function that is used to advance the solution by time dt.
void timestep(ExplicitTimeSteppableObject *const &object_pt, const double &dt)
Overload function that is used to advance time in the object reference by object_pt by an amount dt...
virtual void actions_before_explicit_timestep()
Empty virtual function that can be overloaded to do anything needed before an explicit step...
double & dt(const unsigned &t=0)
Definition: timesteppers.h:137
virtual void add_to_dofs(const double &lambda, const DoubleVector &increment_dofs)
Function that adds the values to the dofs.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
virtual void get_dofs(DoubleVector &dofs) const
Function that gets the values of the dofs in the object.
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
virtual void get_dvaluesdt(DoubleVector &minv_res)
A single virtual function that returns the residuals vector multiplied by the inverse mass matrix...
static double Dummy_time_value
Dummy value of time always set to zero.
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types)...