implicit_midpoint_rule.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// 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
28#include "problem.h"
29
30// Needed for mipoint update...
31#include "mesh.h"
32#include "elements.h"
33#include "nodes.h"
34
35namespace oomph
36{
37 /// This function advances the Data's time history so that
38 /// we can move on to the next timestep
39 void IMRBase::shift_time_values(Data* const& data_pt)
40 {
41 // Loop over the values, set previous values to the previous value, if
42 // not a copy.
43 for (unsigned j = 0, nj = data_pt->nvalue(); j < nj; j++)
44 {
45 if (!data_pt->is_a_copy(j))
46 {
47 for (unsigned t = ndt(); t > 0; t--)
48 {
49 data_pt->set_value(t, j, data_pt->value(t - 1, j));
50 }
51 }
52 }
53 }
54
55 /// This function advances the time history of the positions
56 /// at a node. ??ds Untested: I have no problems with moving nodes.
57 void IMRBase::shift_time_positions(Node* const& node_pt)
58 {
59 // Find the number of coordinates
60 unsigned n_dim = node_pt->ndim();
61 // Find the number of position types
62 unsigned n_position_type = node_pt->nposition_type();
63
64 // Find number of stored timesteps
65 unsigned n_tstorage = ntstorage();
66
67 // Storage for the velocity
68 double velocity[n_position_type][n_dim];
69
70 // If adaptive, find the velocities
71 if (adaptive_flag())
72 {
73 // Loop over the variables
74 for (unsigned i = 0; i < n_dim; i++)
75 {
76 for (unsigned k = 0; k < n_position_type; k++)
77 {
78 // Initialise velocity to zero
79 velocity[k][i] = 0.0;
80 // Loop over all history data
81 for (unsigned t = 0; t < n_tstorage; t++)
82 {
83 velocity[k][i] += Weight(1, t) * node_pt->x_gen(t, k, i);
84 }
85 }
86 }
87 }
88
89 // Loop over the positions
90 for (unsigned i = 0; i < n_dim; i++)
91 {
92 // If the position is not a copy
93 if (node_pt->position_is_a_copy(i) == false)
94 {
95 // Loop over the position types
96 for (unsigned k = 0; k < n_position_type; k++)
97 {
98 // Loop over stored times, and set values to previous values
99 for (unsigned t = ndt(); t > 0; t--)
100 {
101 node_pt->x_gen(t, k, i) = node_pt->x_gen(t - 1, k, i);
102 }
103 }
104 }
105 }
106 }
107
108
109 /// Dummy - just check that the values that
110 /// problem::calculate_predicted_values() has been called right.
112 {
113 if (adaptive_flag())
114 {
115 // Can't do it here, but we can check that the predicted values have
116 // been updated.
118 }
119 }
120
121
123 const unsigned& i)
124 {
125 if (adaptive_flag())
126 {
127 // predicted value is more accurate so just compare with that
128 return data_pt->value(i) - data_pt->value(Predictor_storage_index, i);
129 }
130 else
131 {
132 std::string err("Tried to get the temporal error from a non-adaptive");
133 err += " time stepper.";
134 throw OomphLibError(
135 err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
136 }
137 }
138
139 /// Half the timestep before starting solve
141 {
142 // Check that this is the only time stepper
143#ifdef PARANOID
144 if (problem_pt->ntime_stepper() != 1)
145 {
146 std::string err = "This implementation of midpoint can only work with a ";
147 err += "single time stepper, try using IMR instead (but check ";
148 err += "your Jacobian and residual calculations very carefully for "
149 "compatability).";
150 throw OomphLibError(
151 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
152 }
153#endif
154
155 time_pt()->dt() /= 2;
156 time_pt()->time() -= time_pt()->dt();
157
158 // Set the weights again because we just changed the step size
159 set_weights();
160
161
163 {
164 // Shift the initial guess to the midpoint so that it is an initial
165 // guess for the newton step that we actually take.
166
167 // Optimisation possiblity: here we update all values three time
168 // (initial prediction, copy into initial guess, interpolate to
169 // midpoint), could probably avoid this with more fancy code if
170 // needed.
171
172 // Get dofs at previous time step
173 DoubleVector dof_n;
174 problem_pt->get_dofs(1, dof_n);
175
176 // Update dofs at current step to be the average of current and previous
177 for (unsigned i = 0; i < problem_pt->ndof(); i++)
178 {
179 problem_pt->dof(i) = (problem_pt->dof(i) + dof_n[i]) / 2;
180 }
181 }
182 }
183
184 /// Local (not exported in header) helper function to handle midpoint
185 /// update on a data object.
186 void post_midpoint_update(Data* dat_pt, const bool& update_pinned)
187 {
188 if (!dat_pt->is_a_copy())
189 {
190 for (unsigned j = 0, nj = dat_pt->nvalue(); j < nj; j++)
191 {
192 int eqn = dat_pt->eqn_number(j);
193 if (update_pinned || eqn >= 0)
194 {
195 double ynp1 = 2 * dat_pt->value(0, j) - dat_pt->value(1, j);
196 dat_pt->set_value(0, j, ynp1);
197 }
198 }
199 }
200 }
201
202 /// Take problem from t={n+1/2} to t=n+1 by algebraic update and restore
203 /// time step.
205 {
206#ifdef PARANOID
207 // Do it as dofs too to compare
208 const unsigned ndof = problem_pt->ndof();
209 DoubleVector dof_n, dof_np1;
210 problem_pt->get_dofs(1, dof_n);
211 problem_pt->get_dofs(dof_np1);
212
213 for (unsigned i = 0; i < ndof; i++)
214 {
215 dof_np1[i] += dof_np1[i] - dof_n[i];
216 }
217#endif
218
219 // First deal with global data
220 for (unsigned i = 0, ni = problem_pt->nglobal_data(); i < ni; i++)
221 {
223 }
224
225 // Next element internal data
226 for (unsigned i = 0, ni = problem_pt->mesh_pt()->nelement(); i < ni; i++)
227 {
228 GeneralisedElement* ele_pt = problem_pt->mesh_pt()->element_pt(i);
229 for (unsigned j = 0, nj = ele_pt->ninternal_data(); j < nj; j++)
230 {
232 }
233 }
234
235 // Now the nodes
236 for (unsigned i = 0, ni = problem_pt->mesh_pt()->nnode(); i < ni; i++)
237 {
239 }
240
241 // Update time
242 problem_pt->time_pt()->time() += problem_pt->time_pt()->dt();
243
244 // Return step size to its full value
245 problem_pt->time_pt()->dt() *= 2;
246
247#ifdef PARANOID
248 using namespace StringConversion;
249 DoubleVector actual_dof_np1;
250 problem_pt->get_dofs(actual_dof_np1);
251
252 for (unsigned j = 0; j < actual_dof_np1.nrow(); j++)
253 {
254 if (std::abs(actual_dof_np1[j] - dof_np1[j]) > 1e-8)
255 {
256 std::string err = "Got different values doing midpoint update via "
257 "extracted dofs than doing it in place!";
258 err += to_string(actual_dof_np1[j]) + " vs " + to_string(dof_np1[j]);
259 throw OomphLibError(
260 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
261 }
262 }
263
264#endif
265 }
266
267} // namespace oomph
e
Definition: cfortran.h:571
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
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
unsigned nrow() const
access function to the number of global rows.
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
A Generalised Element class.
Definition: elements.h:73
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:823
void calculate_predicted_values(Data *const &data_pt)
Dummy - just check that the values that problem::calculate_predicted_values() has been called right.
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.
void shift_time_positions(Node *const &node_pt)
This function advances the time history of the positions at a node.
unsigned ndt() const
Number of timestep increments that are required by the scheme.
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.
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.
void set_weights()
Setup weights for time derivative calculations.
void actions_before_timestep(Problem *problem_pt)
Half the timestep before starting solve.
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
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
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:1016
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i). ‘Type’: k; Coordinate direction: i.
Definition: nodes.h:1126
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1674
double & dof(const unsigned &i)
i-th dof in the problem
Definition: problem.h:1813
bool & use_predictor_values_as_initial_guess()
Definition: problem.h:1994
unsigned nglobal_data() const
Return the number of global data values.
Definition: problem.h:1686
Data *& global_data_pt(const unsigned &i)
Return a pointer to the the i-th global data object.
Definition: problem.h:1647
unsigned ntime_stepper() const
Return the number of time steppers.
Definition: problem.h:1680
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
void get_dofs(DoubleVector &dofs) const
Return the vector of dofs, i.e. a vector containing the current values of all unknowns.
Definition: problem.cc:2479
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1504
DenseMatrix< double > Weight
Storage for the weights associated with the timestepper.
Definition: timesteppers.h:237
void check_predicted_values_up_to_date() const
Check that the predicted values are the ones we want.
Definition: timesteppers.h:423
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
bool adaptive_flag() const
Function to indicate whether the scheme is adaptive (false by default)
Definition: timesteppers.h:623
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
int Predictor_storage_index
The time-index in each Data object where predicted values are stored. -1 if not set.
Definition: timesteppers.h:273
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
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.
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).
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void post_midpoint_update(Data *dat_pt, const bool &update_pinned)
Local (not exported in header) helper function to handle midpoint update on a data object.