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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
31 #include "implicit_midpoint_rule.h"
32 #include "problem.h"
33 
34 // Needed for mipoint update...
35 #include "mesh.h"
36 #include "elements.h"
37 #include "nodes.h"
38 
39 namespace oomph
40 {
41 
42  /// \short This function advances the Data's time history so that
43  /// we can move on to the next timestep
44  void IMRBase::shift_time_values(Data* const &data_pt)
45  {
46  //Loop over the values, set previous values to the previous value, if
47  //not a copy.
48  for(unsigned j=0, nj=data_pt->nvalue(); j<nj; j++)
49  {
50  if(!data_pt->is_a_copy(j))
51  {
52  for(unsigned t=ndt(); t>0; t--)
53  {
54  data_pt->set_value(t,j,data_pt->value(t-1,j));
55  }
56  }
57  }
58  }
59 
60  ///\short This function advances the time history of the positions
61  ///at a node. ??ds Untested: I have no problems with moving nodes.
62  void IMRBase::shift_time_positions(Node* const &node_pt)
63  {
64  //Find the number of coordinates
65  unsigned n_dim = node_pt->ndim();
66  //Find the number of position types
67  unsigned n_position_type = node_pt->nposition_type();
68 
69  //Find number of stored timesteps
70  unsigned n_tstorage = ntstorage();
71 
72  //Storage for the velocity
73  double velocity[n_position_type][n_dim];
74 
75  //If adaptive, find the velocities
76  if(adaptive_flag())
77  {
78  //Loop over the variables
79  for(unsigned i=0;i<n_dim;i++)
80  {
81  for(unsigned k=0;k<n_position_type;k++)
82  {
83  //Initialise velocity to zero
84  velocity[k][i] =0.0;
85  //Loop over all history data
86  for(unsigned t=0;t<n_tstorage;t++)
87  {
88  velocity[k][i] += Weight(1,t)*node_pt->x_gen(t,k,i);
89  }
90  }
91  }
92  }
93 
94  //Loop over the positions
95  for(unsigned i=0;i<n_dim;i++)
96  {
97  //If the position is not a copy
98  if(node_pt->position_is_a_copy(i) == false)
99  {
100  //Loop over the position types
101  for(unsigned k=0;k<n_position_type;k++)
102  {
103  //Loop over stored times, and set values to previous values
104  for(unsigned t=ndt();t>0;t--)
105  {
106  node_pt->x_gen(t,k,i) = node_pt->x_gen(t-1,k,i);
107  }
108  }
109  }
110  }
111  }
112 
113 
114  /// Dummy - just check that the values that
115  /// problem::calculate_predicted_values() has been called right.
117  {
118  if(adaptive_flag())
119  {
120  // Can't do it here, but we can check that the predicted values have
121  // been updated.
123  }
124  }
125 
126 
127  double IMRBase::temporal_error_in_value(Data* const &data_pt,
128  const unsigned &i)
129  {
130  if(adaptive_flag())
131  {
132  // predicted value is more accurate so just compare with that
133  return data_pt->value(i) - data_pt->value(Predictor_storage_index, i);
134  }
135  else
136  {
137  std::string err("Tried to get the temporal error from a non-adaptive");
138  err += " time stepper.";
139  throw OomphLibError(err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
140  }
141  }
142 
143  /// Half the timestep before starting solve
145  {
146 
147  // Check that this is the only time stepper
148 #ifdef PARANOID
149  if(problem_pt->ntime_stepper() != 1)
150  {
151  std::string err = "This implementation of midpoint can only work with a ";
152  err += "single time stepper, try using IMR instead (but check ";
153  err += "your Jacobian and residual calculations very carefully for compatability).";
154  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
155  OOMPH_CURRENT_FUNCTION);
156  }
157 #endif
158 
159  time_pt()->dt() /= 2;
160  time_pt()->time() -= time_pt()->dt();
161 
162  // Set the weights again because we just changed the step size
163  set_weights();
164 
165 
166  if(problem_pt->use_predictor_values_as_initial_guess())
167  {
168  // Shift the initial guess to the midpoint so that it is an initial
169  // guess for the newton step that we actually take.
170 
171  // Optimisation possiblity: here we update all values three time
172  // (initial prediction, copy into initial guess, interpolate to
173  // midpoint), could probably avoid this with more fancy code if
174  // needed.
175 
176  // Get dofs at previous time step
177  DoubleVector dof_n;
178  problem_pt->get_dofs(1, dof_n);
179 
180  // Update dofs at current step to be the average of current and previous
181  for(unsigned i=0; i<problem_pt->ndof(); i++)
182  {
183  problem_pt->dof(i) = (problem_pt->dof(i) + dof_n[i])/2;
184  }
185  }
186  }
187 
188  /// Local (not exported in header) helper function to handle midpoint
189  /// update on a data object.
190  void post_midpoint_update(Data* dat_pt, const bool& update_pinned)
191  {
192  if(!dat_pt->is_a_copy())
193  {
194  for(unsigned j=0, nj=dat_pt->nvalue(); j<nj; j++)
195  {
196  int eqn = dat_pt->eqn_number(j);
197  if(update_pinned || eqn >= 0)
198  {
199  double ynp1 = 2*dat_pt->value(0, j) - dat_pt->value(1, j);
200  dat_pt->set_value(0, j, ynp1);
201  }
202  }
203  }
204  }
205 
206  /// Take problem from t={n+1/2} to t=n+1 by algebraic update and restore
207  /// time step.
209  {
210 #ifdef PARANOID
211  // Do it as dofs too to compare
212  const unsigned ndof = problem_pt->ndof();
213  DoubleVector dof_n, dof_np1;
214  problem_pt->get_dofs(1, dof_n);
215  problem_pt->get_dofs(dof_np1);
216 
217  for(unsigned i=0; i<ndof; i++)
218  {
219  dof_np1[i] += dof_np1[i] - dof_n[i];
220  }
221 #endif
222 
223  // First deal with global data
224  for(unsigned i=0, ni=problem_pt->nglobal_data(); i<ni; i++)
225  {
227  }
228 
229  // Next element internal data
230  for(unsigned i=0, ni=problem_pt->mesh_pt()->nelement(); i<ni; i++)
231  {
232  GeneralisedElement* ele_pt = problem_pt->mesh_pt()->element_pt(i);
233  for(unsigned j=0, nj=ele_pt->ninternal_data(); j<nj; j++)
234  {
236  }
237  }
238 
239  // Now the nodes
240  for(unsigned i=0, ni=problem_pt->mesh_pt()->nnode(); i<ni; i++)
241  {
242  post_midpoint_update(problem_pt->mesh_pt()->node_pt(i),
243  Update_pinned);
244  }
245 
246  // Update time
247  problem_pt->time_pt()->time() += problem_pt->time_pt()->dt();
248 
249  // Return step size to its full value
250  problem_pt->time_pt()->dt() *= 2;
251 
252 #ifdef PARANOID
253  using namespace StringConversion;
254  DoubleVector actual_dof_np1;
255  problem_pt->get_dofs(actual_dof_np1);
256 
257  for(unsigned j=0; j<actual_dof_np1.nrow(); j++)
258  {
259  if(std::abs(actual_dof_np1[j] - dof_np1[j]) > 1e-8)
260  {
261  std::string err = "Got different values doing midpoint update via extracted dofs than doing it in place!";
262  err += to_string(actual_dof_np1[j]) + " vs "
263  + to_string(dof_np1[j]);
264  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
265  OOMPH_CURRENT_FUNCTION);
266  }
267  }
268 
269 #endif
270 
271  }
272 
273 } // End of oomph namespace
A Generalised Element class.
Definition: elements.h:76
void shift_time_positions(Node *const &node_pt)
This function advances the time history of the positions at a node.
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:590
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...
cstr elem_len * i
Definition: cfortran.h:607
void actions_before_timestep(Problem *problem_pt)
Half the timestep before starting solve.
The Problem class.
Definition: problem.h:152
void post_midpoint_update(Data *dat_pt, const bool &update_pinned)
char t
Definition: cfortran.h:572
void check_predicted_values_up_to_date() const
Check that the predicted values are the ones we want.
Definition: timesteppers.h:406
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
int Predictor_storage_index
The time-index in each Data object where predicted values are stored. -1 if not set.
Definition: timesteppers.h:263
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1264
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
e
Definition: cfortran.h:575
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
bool Update_pinned
Should we update pinned variables after the half-step?
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:966
bool & use_predictor_values_as_initial_guess()
Definition: problem.h:1871
unsigned ntime_stepper() const
Return the number of time steppers.
Definition: problem.h:1578
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
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:2454
DenseMatrix< double > Weight
Storage for the weights associated with the timestepper.
Definition: timesteppers.h:227
double & dof(const unsigned &i)
i-th dof in the problem
Definition: problem.h:1696
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:255
void calculate_predicted_values(Data *const &data_pt)
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:540
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:267
double & dt(const unsigned &t=0)
Definition: timesteppers.h:137
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1574
void actions_after_timestep(Problem *problem_pt)
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1446
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i).
Definition: nodes.h:1055
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
unsigned nglobal_data() const
Return the number of global data values.
Definition: problem.h:1581
unsigned ndt() const
Number of timestep increments that are required by the scheme.
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:456
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:828
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
Data *& global_data_pt(const unsigned &i)
Return a pointer to the the i-th global data object.
Definition: problem.h:1557
double temporal_error_in_value(Data *const &data_pt, const unsigned &i)
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
bool adaptive_flag() const
Function to indicate whether the scheme is adaptive (false by default)
Definition: timesteppers.h:582
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
void set_weights()
Setup weights for time derivative calculations.
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:291
virtual bool position_is_a_copy() const
Return whether any position coordinate has been copied (always false)
Definition: nodes.h:1048
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)...