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-2022 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"
35
36namespace 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)
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
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.
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.
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...