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-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 
27 #include "implicit_midpoint_rule.h"
28 #include "problem.h"
29 
30 // Needed for mipoint update...
31 #include "mesh.h"
32 #include "elements.h"
33 #include "nodes.h"
34 
35 namespace 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 
122  double IMRBase::temporal_error_in_value(Data* const& data_pt,
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 
162  if (problem_pt->use_predictor_values_as_initial_guess())
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
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
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
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.
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
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
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
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
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:1678
bool & use_predictor_values_as_initial_guess()
Definition: problem.h:2123
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1504
unsigned nglobal_data() const
Return the number of global data values.
Definition: problem.h:1690
double & dof(const unsigned &i)
i-th dof in the problem
Definition: problem.h:1817
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
unsigned ntime_stepper() const
Return the number of time steppers.
Definition: problem.h:1684
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:2565
Data *& global_data_pt(const unsigned &i)
Return a pointer to the the i-th global data object.
Definition: problem.h:1647
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 & dt(const unsigned &t=0)
Return the value of the t-th stored timestep (t=0: present; t>0: previous).
Definition: timesteppers.h:136
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
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.