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// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #ifndef OOMPH_TRAPEZOID_RULE_H
27 #define OOMPH_TRAPEZOID_RULE_H
28 
29 #include "problem.h"
30 #include "timesteppers.h"
31 
32 namespace oomph
33 {
34  /// Trapezoid rule time stepping scheme.
35  ///
36  /// This method requires a value of dy/dt at the initial time. The
37  /// implementation of this calculation is exactly the same as is used for
38  /// explicit time stepping.
39  ///
40  /// The function setup_initial_derivative(Problem* problem_pt) should be
41  /// called after the initial conditions have been set, but before beginning
42  /// time stepping, to compute this initial value of dy/dt.
43  ///
44  /// Warning: moving nodes not implemented (I have no test case).
45  class TR : public TimeStepper
46  {
47  // The standard trapezoid rule is:
48 
49  // y_{n+1} = y_n + dt_n (f(t_n, y_n) + f(t_{n+1}, y_{n+1}))/2
50 
51  // where y is the vector of nodal values and f is the derivative function
52  // (as in standard ODE theory). However we want to calculate time steps
53  // using only one function evaluation (i.e. f) per time step because
54  // function evaluations are expensive. So we require f(t_0, y_0) as
55  // initial input and at each step store f(t_{n+1}, y_{n+1}) as a history
56  // value for use in the calculatio of the next time step.
57 
58 
59  public:
60  /// Constructor, storage for two history derivatives (one for TR and
61  /// one for the predictor step), one history value, present value and
62  /// predicted value.
63  TR(const bool& adaptive = false) : TimeStepper(2 + 2 + 1, 1)
64  {
65  // Set the weight for the zero-th derivative
66  Weight(0, 0) = 1.0;
67 
68  Initial_derivative_set = false;
69 
70  Adaptive_Flag = adaptive;
71 
72  // Initialise adaptivity stuff
73  Predictor_weight.assign(4, 0.0);
74  Error_weight = 0.0;
75 
76  Shift_f = true;
77  }
78 
79  /// Virtual destructor
80  virtual ~TR() {}
81 
82  /// Return the actual order of the scheme
83  unsigned order() const
84  {
85  return 2;
86  }
87 
88  /// Set the weights
89  void set_weights()
90  {
91  double dt = Time_pt->dt(0);
92  Weight(1, 0) = 2.0 / dt;
93  Weight(1, 1) = -2.0 / dt;
94  Weight(1, derivative_index(0)) = -1.0; // old derivative
95  }
96 
98  {
99  double dt = Time_pt->dt(0);
100  double dtprev = Time_pt->dt(1);
101  Error_weight = 1 / (3 * (1 + (dtprev / dt)));
102  }
103 
104  /// Function to set the predictor weights
106  {
107  // Read the value of the time steps
108  double dt = Time_pt->dt(0);
109  double dtprev = Time_pt->dt(1);
110  double dtr = dt / dtprev;
111 
112  // y weights
113  Predictor_weight[0] = 0.0;
114  Predictor_weight[1] = 1.0;
115 
116  // dy weights
117  Predictor_weight[derivative_index(0)] = (dt / 2) * (2 + dtr);
118  Predictor_weight[derivative_index(1)] = -(dt / 2) * dtr;
119  }
120 
121  /// Number of previous values available.
122  unsigned nprev_values() const
123  {
124  return 1;
125  }
126 
127  /// Location in data of derivatives
128  unsigned derivative_index(const unsigned& t) const
129  {
130 #ifdef PARANOID
131  if (t >= 2)
132  {
133  std::string err = "Only storing two derivatives.";
134  throw OomphLibError(
135  err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
136  }
137 #endif
138  return t + 2;
139  }
140 
141  /// Location of predicted value
142  unsigned predicted_value_index() const
143  {
144  return derivative_index(1) + 1;
145  }
146 
147  /// Number of timestep increments that need to be stored by the scheme
148  unsigned ndt() const
149  {
150  return 2;
151  }
152 
153  /// Initialise the time-history for the Data values, corresponding
154  /// to an impulsive start.
155  void assign_initial_values_impulsive(Data* const& data_pt)
156  {
157  throw OomphLibError("Function not yet implemented",
158  OOMPH_EXCEPTION_LOCATION,
159  OOMPH_CURRENT_FUNCTION);
160  }
161 
162  /// Initialise the time-history for the nodal positions
163  /// corresponding to an impulsive start.
165  {
166  throw OomphLibError("Function not yet implemented",
167  OOMPH_EXCEPTION_LOCATION,
168  OOMPH_CURRENT_FUNCTION);
169  }
170 
171  void actions_after_timestep(Problem* problem_pt)
172  {
173  // only ever want this to be true for one step
174  Shift_f = true;
175  }
176 
178  {
179 #ifdef PARANOID
181  {
182  std::string err = "Initial derivative of TR has not been set";
183  throw OomphLibError(
184  err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
185  }
186 #endif
187  }
188 
190  {
191  oomph_info
192  << "Solving for derivative at initial time."
193  << " Warning: if residual is not in the correct form this may fail."
194  << std::endl;
195 
196  // Get the derivative at initial time
197  DoubleVector f0;
198  problem_pt->get_dvaluesdt(f0);
199 
200  // store in derivatives slot ready for use in timestepping.
201  problem_pt->set_dofs(this->derivative_index(0), f0);
202 
203  Initial_derivative_set = true;
204  Shift_f = false;
205  }
206 
207 
208  /// This function updates the Data's time history so that
209  /// we can advance to the next timestep.
210  void shift_time_values(Data* const& data_pt)
211  {
212  // Find number of values stored
213  unsigned n_value = data_pt->nvalue();
214 
215  // Get previous step dt, time_pt has already been shifted so it's in
216  // slot 1.
217  double dtn = time_pt()->dt(1);
218 
219  // Loop over the values
220  for (unsigned j = 0; j < n_value; j++)
221  {
222  // Set previous values to the previous value, if not a copy
223  if (data_pt->is_a_copy(j) == false)
224  {
225  if (Shift_f)
226  {
227  // Calculate time derivative at step n by rearranging the TR
228  // formula.
229  double fnm1 = data_pt->value(derivative_index(0), j);
230  double ynm1 = data_pt->value(1, j);
231  double yn = data_pt->value(0, j);
232  double fn = (2 / dtn) * (yn - ynm1) - fnm1;
233 
234  data_pt->set_value(derivative_index(0), j, fn);
235 
236  // Shift time derivative
237  data_pt->set_value(derivative_index(1), j, fnm1);
238  }
239  else
240  {
241  std::cout << "didn't shift derivatives" << std::endl;
242  }
243 
244  // Shift value
245  data_pt->set_value(1, j, data_pt->value(0, j));
246  }
247  }
248  }
249 
250 
252 
253  bool Shift_f;
254 
255  /// This function advances the time history of the positions at a
256  /// node.
257  void shift_time_positions(Node* const& node_pt)
258  {
259  // do nothing, not implemented for moving nodes
260  }
261 
262 
263  /// Function to calculate predicted positions at a node
264  void calculate_predicted_positions(Node* const& node_pt)
265  {
266  // Loop over the dimensions
267  unsigned n_dim = node_pt->ndim();
268  for (unsigned j = 0; j < n_dim; j++)
269  {
270  // If the node is not copied
271  if (!node_pt->position_is_a_copy(j))
272  {
273  // Initialise the predictor to zero
274  double predicted_value = 0.0;
275 
276  // Loop over all the stored data and add appropriate values
277  // to the predictor
278  for (unsigned i = 1; i < 4; i++)
279  {
280  predicted_value += node_pt->x(i, j) * Predictor_weight[i];
281  }
282 
283  // Store the predicted value
284  node_pt->x(predicted_value_index(), j) = predicted_value;
285  }
286  }
287  }
288 
289  /// Function to calculate predicted data values in a Data object
290  void calculate_predicted_values(Data* const& data_pt)
291  {
292  // Loop over the values
293  unsigned n_value = data_pt->nvalue();
294  for (unsigned j = 0; j < n_value; j++)
295  {
296  // If the value is not copied
297  if (!data_pt->is_a_copy(j))
298  {
299  // Loop over all the stored data and add appropriate
300  // values to the predictor
301  double predicted_value = 0.0;
302  for (unsigned i = 1; i < 4; i++)
303  {
304  predicted_value += data_pt->value(i, j) * Predictor_weight[i];
305  }
306 
307  // Store the predicted value
308  data_pt->set_value(predicted_value_index(), j, predicted_value);
309  }
310  }
311  }
312 
313 
314  /// Compute the error in the position i at a node
315  double temporal_error_in_position(Node* const& node_pt, const unsigned& i)
316  {
317  return Error_weight *
318  (node_pt->x(i) - node_pt->x(predicted_value_index(), i));
319  }
320 
321  /// Compute the error in the value i in a Data structure
322  double temporal_error_in_value(Data* const& data_pt, const unsigned& i)
323  {
324  return Error_weight *
325  (data_pt->value(i) - data_pt->value(predicted_value_index(), i));
326  }
327 
328 
329  private:
330  /// Private data for the predictor weights
332 
333  /// Private data for the error weight
334  double Error_weight;
335 
336  /// Broken copy constructor
337  TR(const TR& dummy) = delete;
338 
339  /// Broken assignment operator
340  void operator=(const TR& dummy) = delete;
341  };
342 
343 } // namespace oomph
344 
345 #endif
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
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:271
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:253
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
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:293
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
virtual bool position_is_a_copy() const
Return whether any position coordinate has been copied (always false)
Definition: nodes.h:1113
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
virtual void get_dvaluesdt(DoubleVector &f)
Get the time derivative of all values (using get_inverse_mass_matrix_times_residuals(....
Definition: problem.cc:3770
void set_dofs(const DoubleVector &dofs)
Set the values of the dofs.
Definition: problem.cc:3497
Trapezoid rule time stepping scheme.
unsigned order() const
Return the actual order of the scheme.
void shift_time_positions(Node *const &node_pt)
This function advances the time history of the positions at a node.
Vector< double > Predictor_weight
Private data for the predictor weights.
void calculate_predicted_values(Data *const &data_pt)
Function to calculate predicted data values in a Data object.
virtual ~TR()
Virtual destructor.
unsigned predicted_value_index() const
Location of predicted value.
void assign_initial_positions_impulsive(Node *const &node_pt)
Initialise the time-history for the nodal positions corresponding to an impulsive start.
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.
void operator=(const TR &dummy)=delete
Broken assignment operator.
double Error_weight
Private data for the error weight.
void calculate_predicted_positions(Node *const &node_pt)
Function to calculate predicted positions at a node.
void assign_initial_values_impulsive(Data *const &data_pt)
Initialise the time-history for the Data values, corresponding to an impulsive start.
void actions_before_timestep(Problem *problem_pt)
Interface for any actions that need to be performed before a time step.
void set_error_weights()
Set the weights for the error computation, (currently empty – overwrite for specific scheme)
void setup_initial_derivative(Problem *problem_pt)
unsigned derivative_index(const unsigned &t) const
Location in data of derivatives.
unsigned ndt() const
Number of timestep increments that need to be stored by the scheme.
TR(const bool &adaptive=false)
Constructor, storage for two history derivatives (one for TR and one for the predictor step),...
unsigned nprev_values() const
Number of previous values available.
bool Initial_derivative_set
void actions_after_timestep(Problem *problem_pt)
Interface for any actions that need to be performed after a time step.
double temporal_error_in_position(Node *const &node_pt, const unsigned &i)
Compute the error in the position i at a node.
void set_predictor_weights()
Function to set the predictor weights.
TR(const TR &dummy)=delete
Broken copy constructor.
void set_weights()
Set the weights.
double temporal_error_in_value(Data *const &data_pt, const unsigned &i)
Compute the error in the value i in a Data structure.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
DenseMatrix< double > Weight
Storage for the weights associated with the timestepper.
Definition: timesteppers.h:237
Time * Time_pt
Pointer to discrete time storage scheme.
Definition: timesteppers.h:234
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
bool Adaptive_Flag
Boolean variable to indicate whether the timestepping scheme can be adaptive.
Definition: timesteppers.h:245
double & dt(const unsigned &t=0)
Return the value of the t-th stored timestep (t=0: present; t>0: previous).
Definition: timesteppers.h:136
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...