trapezoid_rule.h
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 #ifndef OOMPH_TRAPEZOID_RULE_H
31 #define OOMPH_TRAPEZOID_RULE_H
32 
33 #include "problem.h"
34 #include "timesteppers.h"
35 
36 namespace oomph
37 {
38 
39  /// Trapezoid rule time stepping scheme.
40  ///
41  /// This method requires a value of dy/dt at the initial time. The
42  /// implementation of this calculation is exactly the same as is used for
43  /// explicit time stepping.
44  ///
45  /// The function setup_initial_derivative(Problem* problem_pt) should be
46  /// called after the initial conditions have been set, but before beginning
47  /// time stepping, to compute this initial value of dy/dt.
48  ///
49  /// Warning: moving nodes not implemented (I have no test case).
50  class TR : public TimeStepper
51  {
52  // The standard trapezoid rule is:
53 
54  // y_{n+1} = y_n + dt_n (f(t_n, y_n) + f(t_{n+1}, y_{n+1}))/2
55 
56  // where y is the vector of nodal values and f is the derivative function
57  // (as in standard ODE theory). However we want to calculate time steps
58  // using only one function evaluation (i.e. f) per time step because
59  // function evaluations are expensive. So we require f(t_0, y_0) as
60  // initial input and at each step store f(t_{n+1}, y_{n+1}) as a history
61  // value for use in the calculatio of the next time step.
62 
63 
64  public:
65  /// Constructor, storage for two history derivatives (one for TR and
66  /// one for the predictor step), one history value, present value and
67  /// predicted value.
68  TR(const bool& adaptive=false) : TimeStepper(2+2+1, 1)
69  {
70  //Set the weight for the zero-th derivative
71  Weight(0,0) = 1.0;
72 
73  Initial_derivative_set = false;
74 
75  Adaptive_Flag = adaptive;
76 
77  // Initialise adaptivity stuff
78  Predictor_weight.assign(4, 0.0);
79  Error_weight = 0.0;
80 
81  Shift_f = true;
82  }
83 
84  /// Virtual destructor
85  virtual ~TR() {}
86 
87  ///Return the actual order of the scheme
88  unsigned order() const {return 2;}
89 
90  /// Set the weights
91  void set_weights()
92  {
93  double dt = Time_pt->dt(0);
94  Weight(1,0) = 2.0/dt;
95  Weight(1,1) = -2.0/dt;
96  Weight(1,derivative_index(0)) = -1.0; // old derivative
97  }
98 
100  {
101  double dt=Time_pt->dt(0);
102  double dtprev=Time_pt->dt(1);
103  Error_weight = 1/(3*(1 + (dtprev/dt)));
104  }
105 
106  /// Function to set the predictor weights
108  {
109  // Read the value of the time steps
110  double dt=Time_pt->dt(0);
111  double dtprev=Time_pt->dt(1);
112  double dtr = dt/dtprev;
113 
114  // y weights
115  Predictor_weight[0] = 0.0;
116  Predictor_weight[1] = 1.0;
117 
118  // dy weights
119  Predictor_weight[derivative_index(0)] = (dt/2)*(2 + dtr);
120  Predictor_weight[derivative_index(1)] = -(dt/2)*dtr;
121  }
122 
123  /// Number of previous values available.
124  unsigned nprev_values() const {return 1;}
125 
126  /// Location in data of derivatives
127  unsigned derivative_index(const unsigned& t) const
128  {
129 #ifdef PARANOID
130  if(t >= 2)
131  {
132  std::string err = "Only storing two derivatives.";
133  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
134  OOMPH_CURRENT_FUNCTION);
135  }
136 #endif
137  return t + 2;
138  }
139 
140  /// Location of predicted value
141  unsigned predicted_value_index() const {return derivative_index(1)+1;}
142 
143  /// Number of timestep increments that need to be stored by the scheme
144  unsigned ndt() const {return 2;}
145 
146  /// \short Initialise the time-history for the Data values, corresponding
147  /// to an impulsive start.
148  void assign_initial_values_impulsive(Data* const &data_pt)
149  {
150  throw OomphLibError("Function not yet implemented",
151  OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
152  }
153 
154  /// \short Initialise the time-history for the nodal positions
155  /// corresponding to an impulsive start.
157  {
158  throw OomphLibError("Function not yet implemented",
159  OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
160  }
161 
162  void actions_after_timestep(Problem* problem_pt)
163  {
164  // only ever want this to be true for one step
165  Shift_f = true;
166  }
167 
169  {
170 #ifdef PARANOID
172  {
173  std::string err = "Initial derivative of TR has not been set";
174  throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
175  OOMPH_EXCEPTION_LOCATION);
176  }
177 #endif
178  }
179 
181  {
182  oomph_info << "Solving for derivative at initial time."
183  << " Warning: if residual is not in the correct form this may fail."
184  << std::endl;
185 
186  // Get the derivative at initial time
187  DoubleVector f0;
188  problem_pt->get_dvaluesdt(f0);
189 
190  // store in derivatives slot ready for use in timestepping.
191  problem_pt->set_dofs(this->derivative_index(0), f0);
192 
193  Initial_derivative_set = true;
194  Shift_f = false;
195  }
196 
197 
198  /// \short This function updates the Data's time history so that
199  /// we can advance to the next timestep.
200  void shift_time_values(Data* const &data_pt)
201  {
202  //Find number of values stored
203  unsigned n_value = data_pt->nvalue();
204 
205  // Get previous step dt, time_pt has already been shifted so it's in
206  // slot 1.
207  double dtn = time_pt()->dt(1);
208 
209  //Loop over the values
210  for(unsigned j=0; j<n_value; j++)
211  {
212  //Set previous values to the previous value, if not a copy
213  if(data_pt->is_a_copy(j) == false)
214  {
215  if(Shift_f)
216  {
217  // Calculate time derivative at step n by rearranging the TR
218  // formula.
219  double fnm1 = data_pt->value(derivative_index(0), j);
220  double ynm1 = data_pt->value(1, j);
221  double yn = data_pt->value(0, j);
222  double fn = (2/dtn)*(yn - ynm1) - fnm1;
223 
224  data_pt->set_value(derivative_index(0), j, fn);
225 
226  // Shift time derivative
227  data_pt->set_value(derivative_index(1), j, fnm1);
228  }
229  else
230  {
231  std::cout << "didn't shift derivatives" << std::endl;
232  }
233 
234  // Shift value
235  data_pt->set_value(1, j,
236  data_pt->value(0, j));
237  }
238  }
239 
240  }
241 
242 
244 
245  bool Shift_f;
246 
247  ///\short This function advances the time history of the positions at a
248  ///node.
249  void shift_time_positions(Node* const &node_pt)
250  {
251  // do nothing, not implemented for moving nodes
252  }
253 
254 
255 
256  /// Function to calculate predicted positions at a node
257  void calculate_predicted_positions(Node* const &node_pt)
258  {
259  // Loop over the dimensions
260  unsigned n_dim = node_pt->ndim();
261  for(unsigned j=0;j<n_dim;j++)
262  {
263  // If the node is not copied
264  if(!node_pt->position_is_a_copy(j))
265  {
266  // Initialise the predictor to zero
267  double predicted_value = 0.0;
268 
269  // Loop over all the stored data and add appropriate values
270  // to the predictor
271  for(unsigned i=1;i<4;i++)
272  {
273  predicted_value += node_pt->x(i,j)*Predictor_weight[i];
274  }
275 
276  // Store the predicted value
277  node_pt->x(predicted_value_index(), j) = predicted_value;
278  }
279  }
280  }
281 
282  /// Function to calculate predicted data values in a Data object
283  void calculate_predicted_values(Data* const &data_pt)
284  {
285  // Loop over the values
286  unsigned n_value = data_pt->nvalue();
287  for(unsigned j=0;j<n_value;j++)
288  {
289  // If the value is not copied
290  if(!data_pt->is_a_copy(j))
291  {
292  // Loop over all the stored data and add appropriate
293  // values to the predictor
294  double predicted_value = 0.0;
295  for(unsigned i=1;i<4;i++)
296  {
297  predicted_value += data_pt->value(i,j)*Predictor_weight[i];
298  }
299 
300  // Store the predicted value
301  data_pt->set_value(predicted_value_index(), j, predicted_value);
302  }
303  }
304  }
305 
306 
307  /// Compute the error in the position i at a node
308  double temporal_error_in_position(Node* const &node_pt,
309  const unsigned &i)
310  {
311  return Error_weight*(node_pt->x(i) -
312  node_pt->x(predicted_value_index(), i));
313  }
314 
315  /// Compute the error in the value i in a Data structure
316  double temporal_error_in_value(Data* const &data_pt,
317  const unsigned &i)
318  {
319  return Error_weight*(data_pt->value(i)
320  - data_pt->value(predicted_value_index(), i));
321  }
322 
323 
324 
325  private:
326 
327  ///Private data for the predictor weights
329 
330  /// Private data for the error weight
331  double Error_weight;
332 
333  /// Broken copy constructor
334  TR(const TR& dummy)
335  {BrokenCopy::broken_copy("TR");}
336 
337  /// Broken assignment operator
338  void operator=(const TR& dummy)
340 
341  };
342 
343 } // End of oomph namespace
344 
345 #endif
void shift_time_positions(Node *const &node_pt)
This function advances the time history of the positions at a node.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
cstr elem_len * i
Definition: cfortran.h:607
void actions_after_timestep(Problem *problem_pt)
The Problem class.
Definition: problem.h:152
char t
Definition: cfortran.h:572
double Error_weight
Private data for the error weight.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
OomphInfo oomph_info
void set_dofs(const DoubleVector &dofs)
Set the values of the dofs.
Definition: problem.cc:3361
unsigned order() const
Return the actual order of the scheme.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
void setup_initial_derivative(Problem *problem_pt)
void set_predictor_weights()
Function to set the predictor weights.
unsigned predicted_value_index() const
Location of predicted value.
bool Initial_derivative_set
void actions_before_timestep(Problem *problem_pt)
unsigned derivative_index(const unsigned &t) const
Location in data of derivatives.
DenseMatrix< double > Weight
Storage for the weights associated with the timestepper.
Definition: timesteppers.h:227
virtual void get_dvaluesdt(DoubleVector &f)
Get the time derivative of all values (using get_inverse_mass_matrix_times_residuals(..) with all time steppers set to steady) e.g. for use in explicit time steps. The approach used is slighty hacky, beware if you have a residual which is non-linear or implicit in the derivative or if you have overloaded get_jacobian(...).
Definition: problem.cc:3640
double temporal_error_in_value(Data *const &data_pt, const unsigned &i)
Compute the error in the value i in a Data structure.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
virtual bool is_a_copy() const
Return a boolean to indicate whether the Data objact contains any copied values. A base Data object c...
Definition: nodes.h:255
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:540
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
double & dt(const unsigned &t=0)
Definition: timesteppers.h:137
Time * Time_pt
Pointer to discrete time storage scheme.
Definition: timesteppers.h:224
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void assign_initial_positions_impulsive(Node *const &node_pt)
Initialise the time-history for the nodal positions corresponding to an impulsive start...
double temporal_error_in_position(Node *const &node_pt, const unsigned &i)
Compute the error in the position i at a node.
unsigned nprev_values() const
Number of previous values available.
void operator=(const TR &dummy)
Broken assignment operator.
TR(const TR &dummy)
Broken copy constructor.
Vector< double > Predictor_weight
Private data for the predictor weights.
bool Adaptive_Flag
Boolean variable to indicate whether the timestepping scheme can be adaptive.
Definition: timesteppers.h:235
void set_weights()
Set the weights.
virtual ~TR()
Virtual destructor.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
void assign_initial_values_impulsive(Data *const &data_pt)
Initialise the time-history for the Data values, corresponding to an impulsive start.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void set_error_weights()
Set the weights for the error computation, (currently empty – overwrite for specific scheme) ...
TR(const bool &adaptive=false)
void shift_time_values(Data *const &data_pt)
This function updates the Data's time history so that we can advance to the next timestep.
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
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
void calculate_predicted_positions(Node *const &node_pt)
Function to calculate predicted positions at a node.
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:291
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
void calculate_predicted_values(Data *const &data_pt)
Function to calculate predicted data values in a Data object.
virtual bool position_is_a_copy() const
Return whether any position coordinate has been copied (always false)
Definition: nodes.h:1048
unsigned ndt() const
Number of timestep increments that need to be stored by the scheme.