implicit_midpoint_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_IMPLICIT_MIDPOINT_RULE_H
27 #define OOMPH_IMPLICIT_MIDPOINT_RULE_H
28 
29 // oomph-lib headers
30 #include "Vector.h"
31 #include "nodes.h"
32 #include "matrices.h"
33 #include "timesteppers.h"
34 #include "explicit_timesteppers.h"
35 
36 namespace oomph
37 {
38  // Forward decl. so that we can have function of Problem*
39  class Problem;
40 
41 
42  /// Implicit midpoint rule base class for the two implementations.
43  class IMRBase : public TimeStepper
44  {
45  public:
46  /// Constructor with initialisation
47  IMRBase(const bool& adaptive = false)
48  : TimeStepper(2, 1) // initialise weights later
49  {
50  Adaptive_Flag = adaptive;
51  Is_steady = false;
52  Type = "Midpoint method";
54 
55  // If adaptive then we are storing predicted values in slot
56  // 4. Otherwise we aren't storing them so leave it as -1.
57  if (adaptive)
58  {
60  }
61 
62 
63  // Storage for weights needs to be 2x(2 + 0/2) (the +2 is extra history
64  // for ebdf3 predictor if adaptive). This means we provide ways to
65  // calculate the zeroth and first derivatives and in calculations we
66  // use 2 + 0/2 time values.
67 
68  // But actually (for some reason...) the amount of data storage
69  // allocated for the timestepper in each node is determined by the
70  // size of weight. So we need to store an extra dummy entry in order
71  // to have storage space for the predicted value at t_{n+1}.
72 
73  // Just leave space for predictor values anyway, it's not expensive.
74  Weight.resize(2, 5, 0.0);
75 
76  // Assume that set_weights() or make_steady() is called before use to
77  // set up the values stored in Weight.
78  }
79 
80  /// Destructor
81  virtual ~IMRBase() {}
82 
83  /// Setup weights for time derivative calculations.
84  virtual void set_weights() = 0;
85 
86  /// Number of history values to interpolate over to get the "current"
87  /// value.
88  virtual unsigned nprev_values_for_value_at_evaluation_time() const = 0;
89 
90  /// Actual order (accuracy) of the scheme
91  unsigned order() const
92  {
93  return 2;
94  }
95 
96  /// Number of timestep increments that are required by the scheme
97  unsigned ndt() const
98  {
99  return nprev_values();
100  }
101 
102  /// ??ds
103  unsigned nprev_values() const
104  {
105  return 4;
106  }
107 
108  /// This function advances the Data's time history so that
109  /// we can move on to the next timestep
110  void shift_time_values(Data* const& data_pt);
111 
112  /// This function advances the time history of the positions
113  /// at a node.
114  void shift_time_positions(Node* const& node_pt);
115 
116  /// Set the weights for the error computation. This is not used
117  /// by midpoint rule.
119 
120  /// Set the weights for the predictor previous timestep. This is not
121  /// used by midpint rule.
123 
124  /// not implemented (??ds TODO)
125  void assign_initial_values_impulsive(Data* const& data_pt)
126  {
127  std::string err = "Not implemented";
128  throw OomphLibError(
129  err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
130  }
132  {
133  std::string err = "Not implemented";
134  throw OomphLibError(
135  err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
136  }
137 
138 
139  void calculate_predicted_positions(Node* const& node_pt)
140  {
141  std::string err = "Not implemented";
142  throw OomphLibError(
143  err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
144  }
145 
146  double temporal_error_in_position(Node* const& node_pt, const unsigned& i)
147  {
148  std::string err = "Not implemented";
149  throw OomphLibError(
150  err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
151  }
152 
153  // Adaptivity
154  void calculate_predicted_values(Data* const& data_pt);
155  double temporal_error_in_value(Data* const& data_pt, const unsigned& i);
156  };
157 
158  /// The "real" implementation of the implicit midpoint rule. Implemented
159  /// by calculation of residuals etc. at half step. This requires
160  /// non-trivial modifications to the element's residual and Jacobian
161  /// calculation functions to interpolate values to the midpoint. As such
162  /// IMRByBDF should be preferred.
163  ///
164  /// However this class must be used when multiple different time steppers
165  /// are being used simultaneously for different parts of the problem.
166  class IMR : public IMRBase
167  {
168  public:
169  /// Common mistakes when using this implementation of midpoint:
170  /// * Didn't include the d_u_evaltime_by_du_np1 term in the Jacobian.
171  /// * Didn't properly interpolate time/values/x/derivatives in
172  /// jacobian/residual (all of these MUST be evaluated at the midpoint).
173 
174 
175  /// Constructor with initialisation
176  IMR(const bool& adaptive = false) : IMRBase(adaptive) {}
177 
178  /// Destructor, predictor_pt handled by base
179  virtual ~IMR() {}
180 
181  /// Setup weights for time derivative calculations.
182  void set_weights()
183  {
184  // Set the weights for zero-th derivative (i.e. the value to use in
185  // newton solver calculations, implicit midpoint method uses the
186  // average of previous and current values).
187  Weight(0, 0) = 0.5;
188  Weight(0, 1) = 0.5;
189 
190  // Set the weights for 1st time derivative
191  double dt = Time_pt->dt(0);
192  Weight(1, 0) = 1.0 / dt;
193  Weight(1, 1) = -1.0 / dt;
194  }
195 
196  /// Number of history values to interpolate over to get the
197  /// "current" value.
199  {
200  return 2;
201  }
202 
203  private:
204  /// Inaccessible copy constructor.
205  IMR(const IMR& dummy) {}
206 
207  /// Inaccessible assignment operator.
208  void operator=(const IMR& dummy) {}
209  };
210 
211 
212  /// Implementation of implicit midpoint rule by taking half a step of bdf1
213  /// then applying an update to all dofs. This implementation *should* work
214  /// with any existing problem for which the BDF methods work.
215  ///
216  /// The exception is when multiple different time steppers are being used
217  /// simultaneously for different parts of the problem. In this case the
218  /// IMR class should be used.
219  class IMRByBDF : public IMRBase
220  {
221  public:
222  /// Constructor with initialisation
223  IMRByBDF(const bool& adaptive = false) : IMRBase(adaptive)
224  {
225  Update_pinned = true;
226  }
227 
228  /// Destructor
229  virtual ~IMRByBDF() {}
230 
231  /// Setup weights for time derivative calculations.
232  void set_weights()
233  {
234  // Use weights from bdf1
235  double dt = Time_pt->dt(0);
236  Weight(1, 0) = 1.0 / dt;
237  Weight(1, 1) = -1.0 / dt;
238  }
239 
240  /// Number of history values to interpolate over to get the
241  /// "current" value. Evaluation time is the end of the bdf1 "half-step",
242  /// so only need one value as normal.
244  {
245  return 1;
246  }
247 
248  /// Half the timestep before starting solve
249  void actions_before_timestep(Problem* problem_pt);
250 
251  /// Take problem from t={n+1/2} to t=n+1 by algebraic update and restore
252  /// time step.
253  void actions_after_timestep(Problem* problem_pt);
254 
255  /// Should we update pinned variables after the half-step?
257 
258  private:
259  /// Inaccessible copy constructor.
260  IMRByBDF(const IMRByBDF& dummy) {}
261 
262  /// Inaccessible assignment operator.
263  void operator=(const IMRByBDF& dummy) {}
264  };
265 
266 } // namespace oomph
267 
268 #endif
cstr elem_len * i
Definition: cfortran.h:603
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
void resize(const unsigned long &n)
Resize to a square nxn matrix; any values already present will be transfered.
Definition: matrices.h:498
Implicit midpoint rule base class for the two implementations.
virtual ~IMRBase()
Destructor.
void calculate_predicted_values(Data *const &data_pt)
Dummy - just check that the values that problem::calculate_predicted_values() has been called right.
void assign_initial_values_impulsive(Data *const &data_pt)
not implemented (??ds TODO)
void shift_time_values(Data *const &data_pt)
This function advances the Data's time history so that we can move on to the next timestep.
virtual void set_weights()=0
Setup weights for time derivative calculations.
void shift_time_positions(Node *const &node_pt)
This function advances the time history of the positions at a node.
void set_predictor_weights()
Set the weights for the predictor previous timestep. This is not used by midpint rule.
double temporal_error_in_position(Node *const &node_pt, const unsigned &i)
Compute the error in the position i at a node zero here – overwrite for specific scheme.
unsigned ndt() const
Number of timestep increments that are required by the scheme.
void set_error_weights()
Set the weights for the error computation. This is not used by midpoint rule.
unsigned order() const
Actual order (accuracy) of the scheme.
virtual unsigned nprev_values_for_value_at_evaluation_time() const =0
Number of history values to interpolate over to get the "current" value.
unsigned nprev_values() const
??ds
double temporal_error_in_value(Data *const &data_pt, const unsigned &i)
Compute the error in the value i in a Data structure zero here – overwrite for specific scheme.
void calculate_predicted_positions(Node *const &node_pt)
Do the predictor step for the positions at a node (currently empty — overwrite for a specific scheme)
void assign_initial_positions_impulsive(Node *const &node_pt)
Initialiset the positions for the node corresponding to an impulsive start.
IMRBase(const bool &adaptive=false)
Constructor with initialisation.
Implementation of implicit midpoint rule by taking half a step of bdf1 then applying an update to all...
void operator=(const IMRByBDF &dummy)
Inaccessible assignment operator.
bool Update_pinned
Should we update pinned variables after the half-step?
void actions_after_timestep(Problem *problem_pt)
Take problem from t={n+1/2} to t=n+1 by algebraic update and restore time step.
unsigned nprev_values_for_value_at_evaluation_time() const
Number of history values to interpolate over to get the "current" value. Evaluation time is the end o...
void set_weights()
Setup weights for time derivative calculations.
IMRByBDF(const IMRByBDF &dummy)
Inaccessible copy constructor.
IMRByBDF(const bool &adaptive=false)
Constructor with initialisation.
virtual ~IMRByBDF()
Destructor.
void actions_before_timestep(Problem *problem_pt)
Half the timestep before starting solve.
The "real" implementation of the implicit midpoint rule. Implemented by calculation of residuals etc....
IMR(const IMR &dummy)
Inaccessible copy constructor.
virtual ~IMR()
Destructor, predictor_pt handled by base.
void operator=(const IMR &dummy)
Inaccessible assignment operator.
unsigned nprev_values_for_value_at_evaluation_time() const
Number of history values to interpolate over to get the "current" value.
IMR(const bool &adaptive=false)
Common mistakes when using this implementation of midpoint:
void set_weights()
Setup weights for time derivative calculations.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
DenseMatrix< double > Weight
Storage for the weights associated with the timestepper.
Definition: timesteppers.h:237
bool Predict_by_explicit_step
Flag: is adaptivity done by taking a separate step using an ExplicitTimeStepper object?
Definition: timesteppers.h:260
Time * Time_pt
Pointer to discrete time storage scheme.
Definition: timesteppers.h:234
std::string Type
String that indicates the type of the timestepper (e.g. "BDF", "Newmark", etc.)
Definition: timesteppers.h:241
bool Is_steady
Bool to indicate if the timestepper is steady, i.e. its time-derivatives evaluate to zero....
Definition: timesteppers.h:251
int Predictor_storage_index
The time-index in each Data object where predicted values are stored. -1 if not set.
Definition: timesteppers.h:273
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...