timesteppers.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 //Include the header
31 #include "timesteppers.h"
32 
33 // Required so that we can delete the predictor.
34 #include "explicit_timesteppers.h"
35 
36 
37 namespace oomph
38 {
39 
40  /// Destructor for time stepper, in .cc file so that we don't have to
41  /// include explicit_timesteppers.h in header.
43  {
44  // Make sure explicit predictor gets deleted if it was set
46  }
47 
48 
49 ////////////////////////////////////////////////////////////////////////
50 ////////////////////////////////////////////////////////////////////////
51 // Static variables for steady timestepper
52 ////////////////////////////////////////////////////////////////////////
53 ////////////////////////////////////////////////////////////////////////
54 
55 template<unsigned NSTEPS>
56 double Steady<NSTEPS>::One=1.0;
57 
58 template<unsigned NSTEPS>
59 double Steady<NSTEPS>::Zero=0.0;
60 
61 template<unsigned NSTEPS>
63 
64 // Force the instantiation for all NSTEPS values that may be
65 // required in other timesteppers
66  template class Steady<0>;
67  template class Steady<1>;
68  template class Steady<2>;
69  template class Steady<3>;
70  template class Steady<4>;
71 
72 
73 ////////////////////////////////////////////////////////////////////////
74 ////////////////////////////////////////////////////////////////////////
75 // Weights for specific bdf schemes
76 ////////////////////////////////////////////////////////////////////////
77 ////////////////////////////////////////////////////////////////////////
78 
79 
80 //=======================================================================
81  ///Assign the values of the weights
82 //=======================================================================
83 template<>
85  {
86  double dt=Time_pt->dt(0);
87  Weight(1,0) = 1.0/dt;
88  Weight(1,1) = -1.0/dt;
89  }
90 
91 
92 //======================================================================
93 /// Set the predictor weights (Gresho and Sani pg. 270)
94 //======================================================================
95 template<>
97 {
98  //If it's adaptive set the predictor weights
99  if(adaptive_flag())
100  {
101  double dt=Time_pt->dt(0);
102 
103  //Set the predictor weights
104  Predictor_weight[0] = 0.0;
105  Predictor_weight[1] = 1.0;
106 
107  // Velocity weight
108  Predictor_weight[2] = dt;
109  }
110 }
111 
112 
113 //=======================================================================
114 ///Calculate the predicted values and store them at the appropriate
115 ///location in the data structure
116 ///This function must be called after the time-values have been shifted!
117 //=======================================================================
118 template<>
120 {
121  //If it's adaptive calculate the values
122  if(adaptive_flag())
123  {
124  //Find number of values
125  unsigned n_value = data_pt->nvalue();
126  //Loop over the values
127  for(unsigned j=0;j<n_value;j++)
128  {
129  //If the value is not copied
130  if(data_pt->is_a_copy(j) == false)
131  {
132  //Now Initialise the predictor to zero
133  double predicted_value = 0.0;
134  //Now loop over all the stored data and add appropriate values
135  //to the predictor
136  for(unsigned i=1;i<3;i++)
137  {
138  predicted_value += data_pt->value(i,j)*Predictor_weight[i];
139  }
140  //Store the predicted value
141  data_pt->set_value(Predictor_storage_index, j, predicted_value);
142  }
143  }
144  }
145 }
146 
147 
148 //=======================================================================
149 ///Calculate predictions for the positions
150 //=======================================================================
151 template<>
153  {
154  //Only do this if adaptive
155  if(adaptive_flag())
156  {
157  //Find number of dimensions of the problem
158  unsigned n_dim = node_pt->ndim();
159  //Loop over the dimensions
160  for(unsigned j=0;j<n_dim;j++)
161  {
162  //If the node is not copied
163  if(node_pt->position_is_a_copy(j) == false)
164  {
165  //Initialise the predictor to zero
166  double predicted_value = 0.0;
167  //Now loop over all the stored data and add appropriate values
168  //to the predictor
169  for(unsigned i=1;i<3;i++)
170  {
171  predicted_value += node_pt->x(i,j)*Predictor_weight[i];
172  }
173  //Store the predicted value
174  node_pt->x(Predictor_storage_index, j) = predicted_value;
175  }
176  }
177  }
178  }
179 
180 
181 //=======================================================================
182 /// Set the error weights (Gresho and Sani pg. 270)
183 //=======================================================================
184 template<>
186 {
187  Error_weight = 0.5;
188 }
189 
190 //===================================================================
191 ///Function to compute the error in position i at node
192 //===================================================================
193 template<>
194 double BDF<1>::temporal_error_in_position(Node* const &node_pt,
195  const unsigned &i)
196 {
197  if(adaptive_flag())
198  {
199  //Just return the error
200  return Error_weight*(node_pt->x(i) - node_pt->x(Predictor_storage_index,i));
201  }
202  else
203  {
204  return 0.0;
205  }
206 }
207 
208 
209 //=========================================================================
210 ///Function to calculate the error in the data value i
211 //=========================================================================
212 template<>
213 double BDF<1>::temporal_error_in_value(Data* const &data_pt, const unsigned &i)
214 {
215  if(adaptive_flag())
216  {
217  //Just return the error
218  return Error_weight*(data_pt->value(i)
219  - data_pt->value(Predictor_storage_index,i));
220  }
221  else
222  {
223  return 0.0;
224  }
225 }
226 
227 //=======================================================================
228 /// Assign the values of the weights. The scheme used is from Gresho &
229 /// Sani, Incompressible Flow and the Finite Element Method
230 /// (advection-diffusion and isothermal laminar flow), 1998, pg. 715 (with
231 /// some algebraic rearrangement).
232 //=======================================================================
233 template<>
235  {
236  double dt=Time_pt->dt(0);
237  double dtprev=Time_pt->dt(1);
238  Weight(1,0) = 1.0/dt + 1.0/(dt + dtprev);
239  Weight(1,1) = -(dt + dtprev)/(dt*dtprev);
240  Weight(1,2) = dt/((dt+dtprev)*dtprev);
241 
242  if(adaptive_flag())
243  {
244  Weight(1,3) = 0.0;
245  Weight(1,4) = 0.0;
246  }
247  }
248 
249 //========================================================================
250 /// Calculate the predictor weights. The scheme used is from Gresho &
251 /// Sani, Incompressible Flow and the Finite Element Method
252 /// (advection-diffusion and isothermal laminar flow), 1998, pg. 715.
253 //========================================================================
254 template<>
256 {
257  //If it's adaptive set the predictor weights
258  if(adaptive_flag())
259  {
260  //Read the value of the previous timesteps
261  double dt=Time_pt->dt(0);
262  double dtprev=Time_pt->dt(1);
263 
264  //Set the predictor weights
265  Predictor_weight[0] = 0.0;
266  Predictor_weight[1] = 1.0 - (dt*dt)/(dtprev*dtprev);
267  Predictor_weight[2] = (dt*dt)/(dtprev*dtprev);
268  //Acceleration term
269  Predictor_weight[3] = (1.0 + dt/dtprev)*dt;
270  }
271 }
272 
273 //=======================================================================
274 ///Calculate the predicted values and store them at the appropriate
275 ///location in the data structure
276 ///This function must be called after the time-values have been shifted!
277 //=======================================================================
278 template<>
280 {
281  //If it's adaptive calculate the values
282  if(adaptive_flag())
283  {
284  //Find number of values
285  unsigned n_value = data_pt->nvalue();
286  //Loop over the values
287  for(unsigned j=0;j<n_value;j++)
288  {
289  //If the value is not copied
290  if(data_pt->is_a_copy(j) == false)
291  {
292  //Now Initialise the predictor to zero
293  double predicted_value = 0.0;
294  //Now loop over all the stored data and add appropriate values
295  //to the predictor
296  for(unsigned i=1;i<4;i++)
297  {
298  predicted_value += data_pt->value(i,j)*Predictor_weight[i];
299  }
300  //Store the predicted value
301  //Note that predictor is stored as the FIFTH entry in this scheme
302  data_pt->set_value(Predictor_storage_index,j,predicted_value);
303  }
304  }
305  }
306 }
307 
308 //=======================================================================
309 ///Calculate predictions for the positions
310 //=======================================================================
311 template<>
313  {
314  //Only do this if adaptive
315  if(adaptive_flag())
316  {
317  //Find number of dimensions of the problem
318  unsigned n_dim = node_pt->ndim();
319  //Loop over the dimensions
320  for(unsigned j=0;j<n_dim;j++)
321  {
322  //If the node is not copied
323  if(node_pt->position_is_a_copy(j) == false)
324  {
325  //Initialise the predictor to zero
326  double predicted_value = 0.0;
327  //Now loop over all the stored data and add appropriate values
328  //to the predictor
329  for(unsigned i=1;i<4;i++)
330  {
331  predicted_value += node_pt->x(i,j)*Predictor_weight[i];
332  }
333  //Store the predicted value
334  node_pt->x(Predictor_storage_index, j) = predicted_value;
335  }
336  }
337  }
338  }
339 
340 
341 //=======================================================================
342 ///Function that sets the error weights
343 //=======================================================================
344 template<>
346 {
347  if(adaptive_flag())
348  {
349  double dt=Time_pt->dt(0);
350  double dtprev=Time_pt->dt(1);
351  //Calculate the error weight
352  Error_weight = pow((1.0 + dtprev/dt),2.0)/
353  (1.0 + 3.0*(dtprev/dt) + 4.0*pow((dtprev/dt),2.0) +
354  2.0*pow((dtprev/dt),3.0));
355  }
356 }
357 
358 
359 //===================================================================
360 ///Function to compute the error in position i at node
361 //===================================================================
362 template<>
364  const unsigned &i)
365 {
366  if(adaptive_flag())
367  {
368  //Just return the error
369  return Error_weight*(node_pt->x(i) - node_pt->x(Predictor_storage_index,i));
370  }
371  else
372  {
373  return 0.0;
374  }
375 }
376 
377 
378 
379 //=========================================================================
380 ///Function to calculate the error in the data value i
381 //=========================================================================
382 template<>
383 double BDF<2>::temporal_error_in_value(Data* const &data_pt, const unsigned &i)
384 {
385  if(adaptive_flag())
386  {
387  //Just return the error
388  return Error_weight*(data_pt->value(i)
389  - data_pt->value(Predictor_storage_index,i));
390  }
391  else
392  {
393  return 0.0;
394  }
395 }
396 
397 
398 //====================================================================
399  ///Assign the values of the weights; pass the value of the timestep
400 //====================================================================
401  template<>
403  {
404 #ifdef PARANOID
405  double dt0=Time_pt->dt(0);
406  for (unsigned i=0;i<Time_pt->ndt();i++)
407  {
408  if (dt0!=Time_pt->dt(i))
409  {
410  throw OomphLibError(
411  "BDF4 currently only works for fixed timesteps \n",
412  OOMPH_CURRENT_FUNCTION,
413  OOMPH_EXCEPTION_LOCATION);
414  }
415  }
416 #endif
417  double dt=Time_pt->dt(0);
418  Weight(1,0) = 25.0/12.0/dt;
419  Weight(1,1) = -48.0/12.0/dt;
420  Weight(1,2) = 36.0/12.0/dt;
421  Weight(1,3) = -16.0/12.0/dt;
422  Weight(1,4) = 3.0/12.0/dt;
423  }
424 
425 
426 //======================================================================
427 ///Calculate the predictor weights
428 //======================================================================
429 template<>
431 {
432  //Read the value of the previous timesteps
433  //double dt=Time_pt->dt(0);
434  //double dtprev=Time_pt->dt(1);
435 
436  throw OomphLibError("Not implemented yet",
437  OOMPH_CURRENT_FUNCTION,
438  OOMPH_EXCEPTION_LOCATION);
439 }
440 
441 //=======================================================================
442 ///Calculate the predicted values and store them at the appropriate
443 ///location in the data structure
444 ///This function must be called after the time-values have been shifted!
445 //=======================================================================
446 template<>
448 {
449  throw OomphLibError("Not implemented yet",
450  OOMPH_CURRENT_FUNCTION,
451  OOMPH_EXCEPTION_LOCATION);
452 
453 }
454 
455 ///Calculate predictions for the positions
456 template<>
458  {
459  throw OomphLibError("Not implemented yet",
460  OOMPH_CURRENT_FUNCTION,
461  OOMPH_EXCEPTION_LOCATION);
462  }
463 
464 ///Function that sets the error weights
465 template<>
467 {
468  throw OomphLibError("Not implemented yet",
469  OOMPH_CURRENT_FUNCTION,
470  OOMPH_EXCEPTION_LOCATION);
471 }
472 
473 //===================================================================
474 ///Function to compute the error in position i at node
475 //===================================================================
476 template<>
477 double BDF<4>::temporal_error_in_position(Node* const &node_pt,
478  const unsigned &i)
479 {
480  throw OomphLibError("Not implemented yet",
481  OOMPH_CURRENT_FUNCTION,
482  OOMPH_EXCEPTION_LOCATION);
483  return 0.0;
484 }
485 
486 
487 //=========================================================================
488 ///Function to calculate the error in the data value i
489 //=========================================================================
490 template<>
491 double BDF<4>::temporal_error_in_value(Data* const &data_pt, const unsigned &i)
492 {
493  throw OomphLibError("Not implemented yet",
494  OOMPH_CURRENT_FUNCTION,
495  OOMPH_EXCEPTION_LOCATION);
496  return 0.0;
497 }
498 
499 
500 
501 /////////////////////////////////////////////////////////////////////
502 /////////////////////////////////////////////////////////////////////
503 /////////////////////////////////////////////////////////////////////
504 
505 
506 
507 
508 //=========================================================================
509 /// \short Initialise the time-history for the values,
510 /// corresponding to an impulsive start.
511 //=========================================================================
512 template<unsigned NSTEPS>
514 {
515  //Find number of values stored in Data object
516  unsigned n_value = data_pt->nvalue();
517 
518  //Loop over values
519  for(unsigned j=0;j<n_value;j++)
520  {
521  //Set previous values to the initial value, if not a copy
522  if(data_pt->is_a_copy(j) == false)
523  {
524  for (unsigned t=1;t<=NSTEPS;t++)
525  {
526  data_pt->set_value(t,j,data_pt->value(j));
527  }
528  }
529 
530  //Initial velocity and initial acceleration are zero
531  data_pt->set_value(NSTEPS+1,j,0.0);
532  data_pt->set_value(NSTEPS+2,j,0.0);
533  }
534 }
535 
536 
537 //=========================================================================
538 /// \short Initialise the time-history for the values,
539 /// corresponding to an impulsive start.
540 //=========================================================================
541 template<unsigned NSTEPS>
543 {
544  //Find the number of coordinates
545  unsigned n_dim = node_pt->ndim();
546  //Find the number of position types
547  unsigned n_position_type = node_pt->nposition_type();
548 
549  //Loop over values
550  for(unsigned i=0;i<n_dim;i++)
551  {
552  //If not a copy
553  if(node_pt->position_is_a_copy(i) == false)
554  {
555  //Loop over the position types
556  for(unsigned k=0;k<n_position_type;k++)
557  {
558  //Set previous value to the initial value, if not a copy
559  for (unsigned t=1;t<=NSTEPS;t++)
560  {
561  node_pt->x_gen(t,k,i) = node_pt->x_gen(k,i);
562  }
563 
564  //Initial velocity and initial acceleration are zero
565  node_pt->x_gen(NSTEPS+1,k,i) = 0.0;
566  node_pt->x_gen(NSTEPS+2,k,i) = 0.0;
567  }
568  }
569  }
570 }
571 
572 
573 
574 //=========================================================================
575 /// \short Initialise the time-history for the Data values,
576 /// so that the Newmark representations for current veloc and
577 /// acceleration are exact.
578 //=========================================================================
579 template<unsigned NSTEPS>
582  initial_value_fct,
584  initial_veloc_fct,
586  initial_accel_fct)
587 {
588  // Set weights
589  set_weights();
590 
591 #ifdef PARANOID
592  //Check if the 3 vectors of functions have the same size
593  if(initial_value_fct.size() != initial_veloc_fct.size() ||
594  initial_value_fct.size() != initial_accel_fct.size())
595  {
596  throw OomphLibError(
597  "The Vectors of fcts for values, velocities and acceleration must be the same size",
598  OOMPH_CURRENT_FUNCTION,
599  OOMPH_EXCEPTION_LOCATION);
600  }
601 #endif
602 
603  //The number of functions in each vector (they should be the same)
604  unsigned n_fcts = initial_value_fct.size();
605 
606 #ifdef PARANOID
607 
608  //If there are more data values at the node than functions, issue a warning
609  unsigned n_value = data_pt->nvalue();
610  if(n_value>n_fcts && !Shut_up_in_assign_initial_data_values)
611  {
612  std::stringstream message;
613  message << "There are more data values than initial value fcts.\n";
614  message << "Only the first " << n_fcts << " data values will be set\n";
615  OomphLibWarning(message.str(),
616  OOMPH_CURRENT_FUNCTION,
617  OOMPH_EXCEPTION_LOCATION);
618  }
619 #endif
620 
621  //Loop over elements of the function vectors
622  for(unsigned j=0;j<n_fcts;j++)
623  {
624  if (initial_value_fct[j]==0)
625  {
626 #ifdef PARANOID
627  if(!Shut_up_in_assign_initial_data_values)
628  {
629  std::stringstream message;
630  message << "Ignoring value " << j << " in assignment of ic.\n";
631  OomphLibWarning(message.str(),
632  "Newmark<NSTEPS>::assign_initial_data_values",
633  OOMPH_EXCEPTION_LOCATION);
634  }
635 #endif
636  }
637  else
638  {
639  // Value itself at current and previous times
640  for (unsigned t=0;t<=NSTEPS;t++)
641  {
642  double time_local=Time_pt->time(t);
643  data_pt->set_value(t,j,initial_value_fct[j](time_local));
644  }
645 
646  // Now, rather than assigning the values for veloc and accel
647  // in the Newmark scheme directly, we solve a linear system
648  // to determine the values required to make the Newmark
649  // representations of the veloc and accel at the current (!)
650  // time are exact.
651 
652  // Initial time: The value itself
653  double time_=Time_pt->time();
654  double U = initial_value_fct[j](time_);
655 
656  // Value itself at previous time
657  time_=Time_pt->time(1);
658  double U0 = initial_value_fct[j](time_);
659 
660  // Veloc (time deriv) at present time -- this is what the
661  // Newmark scheme should return!
662  time_=Time_pt->time(0);
663  double Udot = initial_veloc_fct[j](time_);
664 
665  // Acccel (2nd time deriv) at present time -- this is what the
666  // Newmark scheme should return!
667  time_=Time_pt->time(0);
668  double Udotdot = initial_accel_fct[j](time_);
669 
670  Vector<double> vect(2);
671  vect[0]=Udotdot-Weight(2,0)*U-Weight(2,1)*U0;
672  vect[1]=Udot -Weight(1,0)*U-Weight(1,1)*U0;
673 
674  DenseDoubleMatrix matrix(2,2);
675  matrix(0,0)=Weight(2,NSTEPS+1);
676  matrix(0,1)=Weight(2,NSTEPS+2);
677  matrix(1,0)=Weight(1,NSTEPS+1);
678  matrix(1,1)=Weight(1,NSTEPS+2);
679 
680  matrix.solve(vect);
681 
682  // Discrete veloc (time deriv) at previous time , to be entered into the
683  // Newmark scheme so that its prediction for the *current* veloc
684  // and accel is correct.
685  data_pt->set_value(NSTEPS+1,j,vect[0]);
686 
687  // Discrete veloc (2nd time deriv) at previous time, to be entered into
688  // the Newmark scheme so that its prediction for the *current* veloc
689  // and accel is correct.
690  data_pt->set_value(NSTEPS+2,j,vect[1]);
691  }
692  }
693 }
694 
695 
696 
697 //=========================================================================
698 /// \short Initialise the time-history for the Data values,
699 /// so that the Newmark representations for current veloc and
700 /// acceleration are exact.
701 //=========================================================================
702 template<unsigned NSTEPS>
704  Node* const & node_pt,
705  Vector<NodeInitialConditionFctPt> initial_value_fct,
706  Vector<NodeInitialConditionFctPt> initial_veloc_fct,
707  Vector<NodeInitialConditionFctPt> initial_accel_fct)
708 {
709  // Set weights
710  set_weights();
711 
712 
713 #ifdef PARANOID
714  //Check if the 3 vectors of functions have the same size
715  if(initial_value_fct.size() != initial_veloc_fct.size() ||
716  initial_value_fct.size() != initial_accel_fct.size())
717  {
718  throw OomphLibError("The Vectors of fcts for values, velocities and acceleration must be the same size",
719  "Newmark<NSTEPS>::assign_initial_data_values",
720  OOMPH_EXCEPTION_LOCATION);
721  }
722 #endif
723 
724  //The number of functions in each vector (they should be the same)
725  unsigned n_fcts = initial_value_fct.size();
726 
727 #ifdef PARANOID
728 
729  //If there are more data values at the node than functions, issue a warning
730  unsigned n_value = node_pt->nvalue();
731  if(n_value>n_fcts && !Shut_up_in_assign_initial_data_values)
732  {
733  std::stringstream message;
734  message << "There are more nodal values than initial value fcts.\n";
735  message << "Only the first " << n_fcts << " nodal values will be set\n";
736  OomphLibWarning(message.str(),
737  OOMPH_CURRENT_FUNCTION,
738  OOMPH_EXCEPTION_LOCATION);
739  }
740 #endif
741 
742  // Get current nodal coordinates
743  unsigned n_dim=node_pt->ndim();
744  Vector<double> x(n_dim);
745  for (unsigned i=0;i<n_dim;i++) x[i]=node_pt->x(i);
746 
747  //Loop over values
748  for(unsigned j=0;j<n_fcts;j++)
749  {
750  if (initial_value_fct[j]==0)
751  {
752 #ifdef PARANOID
753  if(!Shut_up_in_assign_initial_data_values)
754  {
755  std::stringstream message;
756  message << "Ignoring value " << j << " in assignment of ic.\n";
757  OomphLibWarning(message.str(),
758  OOMPH_CURRENT_FUNCTION,
759  OOMPH_EXCEPTION_LOCATION);
760  }
761 #endif
762  }
763  else
764  {
765  // Value itself at current and previous times
766  for (unsigned t=0;t<=NSTEPS;t++)
767  {
768  double time_local=Time_pt->time(t);
769  node_pt->set_value(t,j,initial_value_fct[j](time_local,x));
770  }
771 
772  // Now, rather than assigning the values for veloc and accel
773  // in the Newmark scheme directly, we solve a linear system
774  // to determine the values required to make the Newmark
775  // representations of the veloc and accel at the current (!)
776  // time are exact.
777 
778  // Initial time: The value itself
779  double time_=Time_pt->time();
780  double U = initial_value_fct[j](time_,x);
781 
782  // Value itself at previous time
783  time_=Time_pt->time(1);
784  double U0 = initial_value_fct[j](time_,x);
785 
786  // Veloc (time deriv) at present time -- this is what the
787  // Newmark scheme should return!
788  time_=Time_pt->time(0);
789  double Udot = initial_veloc_fct[j](time_,x);
790 
791  // Acccel (2nd time deriv) at present time -- this is what the
792  // Newmark scheme should return!
793  time_=Time_pt->time(0);
794  double Udotdot = initial_accel_fct[j](time_,x);
795 
796  Vector<double> vect(2);
797  vect[0]=Udotdot-Weight(2,0)*U-Weight(2,1)*U0;
798  vect[1]=Udot -Weight(1,0)*U-Weight(1,1)*U0;
799 
800  DenseDoubleMatrix matrix(2,2);
801  matrix(0,0)=Weight(2,NSTEPS+1);
802  matrix(0,1)=Weight(2,NSTEPS+2);
803  matrix(1,0)=Weight(1,NSTEPS+1);
804  matrix(1,1)=Weight(1,NSTEPS+2);
805 
806  matrix.solve(vect);
807 
808  // Discrete veloc (time deriv) at previous time , to be entered into the
809  // Newmark scheme so that its prediction for the *current* veloc
810  // and accel is correct.
811  node_pt->set_value(NSTEPS+1,j,vect[0]);
812 
813  // Discrete veloc (2nd time deriv) at previous time, to be entered into
814  // the Newmark scheme so that its prediction for the *current* veloc
815  // and accel is correct.
816  node_pt->set_value(NSTEPS+2,j,vect[1]);
817  }
818  }
819 }
820 
821 
822 //=========================================================================
823 /// \short First step in a two-stage procedure to assign
824 /// the history values for the Newmark scheme so
825 /// that the veloc and accel that are computed by the scheme
826 /// are correct at the current time.
827 ///
828 /// Call this function for t_deriv=0,1,2,3.
829 /// When calling with
830 /// - t_deriv=0 : data_pt(0,*) should contain the values at the
831 /// previous timestep
832 /// - t_deriv=1 : data_pt(0,*) should contain the current values;
833 /// they get stored (temporarily) in data_pt(1,*).
834 /// - t_deriv=2 : data_pt(0,*) should contain the current velocities
835 /// (first time derivs); they get stored (temporarily)
836 /// in data_pt(NSTEPS+1,*).
837 /// - t_deriv=3 : data_pt(0,*) should contain the current accelerations
838 /// (second time derivs); they get stored (temporarily)
839 /// in data_pt(NSTEPS+2,*).
840 /// .
841 /// Follow this by calls to
842 /// \code
843 /// assign_initial_data_values_stage2(...)
844 /// \endcode
845 //=========================================================================
846 template<unsigned NSTEPS>
848  Data* const &data_pt)
849 {
850  //Find number of values stored
851  unsigned n_value = data_pt->nvalue();
852 
853  //Loop over values
854  for(unsigned j=0;j<n_value;j++)
855  {
856  // Which deriv are we dealing with?
857  switch (t_deriv)
858  {
859 
860  // 0-th deriv: the value itself at present time
861  case 0:
862 
863  data_pt->set_value(1,j,data_pt->value(j));
864  break;
865 
866  // 1st deriv: the veloc at present time
867  case 1:
868 
869  data_pt->set_value(NSTEPS+1,j,data_pt->value(j));
870  break;
871 
872  // 2nd deriv: the accel at present time
873  case 2:
874 
875  data_pt->set_value(NSTEPS+2,j,data_pt->value(j));
876  break;
877 
878 
879  // None other are possible!
880  default:
881 
882  std::ostringstream error_message_stream;
883  error_message_stream << "t_deriv " << t_deriv
884  << " is not possible with a Newmark scheme "
885  << std::endl;
886 
887  throw OomphLibError(
888  error_message_stream.str(),
889  OOMPH_CURRENT_FUNCTION,
890  OOMPH_EXCEPTION_LOCATION);
891  }
892  }
893 }
894 
895 //=========================================================================
896 /// \short Second step in a two-stage procedure to assign
897 /// the history values for the Newmark scheme so
898 /// that the veloc and accel that are computed by the scheme
899 /// are correct at the current time.
900 ///
901 /// This assigns appropriate values for the "previous
902 /// velocities and accelerations" so that their current
903 /// values, which were defined in assign_initial_data_values_stage1(...),
904 /// are represented exactly by the Newmark scheme.
905 //=========================================================================
906 template<unsigned NSTEPS>
908 {
909  //Find number of values stored
910  unsigned n_value = data_pt->nvalue();
911 
912  //Loop over values
913  for(unsigned j=0;j<n_value;j++)
914  {
915 
916  // Rather than assigning the previous veloc and accel
917  // directly, solve linear system to determine their values
918  // so that the Newmark representations are exact.
919 
920  // The exact value at present time (stored in stage 1)
921  double U = data_pt->value(1,j);
922 
923  // Exact value at previous time
924  double U0 = data_pt->value(0,j);
925 
926  // Veloc (1st time deriv) at present time (stored in stage 1)
927  double Udot = data_pt->value(NSTEPS+1,j);
928 
929  // Acccel (2nd time deriv) at present time (stored in stage 1)
930  double Udotdot = data_pt->value(NSTEPS+2,j);
931 
932  Vector<double> vect(2);
933  vect[0]=Udotdot-Weight(2,0)*U-Weight(2,1)*U0;
934  vect[1]=Udot -Weight(1,0)*U-Weight(1,1)*U0;
935 
936  DenseDoubleMatrix matrix(2,2);
937  matrix(0,0)=Weight(2,NSTEPS+1);
938  matrix(0,1)=Weight(2,NSTEPS+2);
939  matrix(1,0)=Weight(1,NSTEPS+1);
940  matrix(1,1)=Weight(1,NSTEPS+2);
941 
942  matrix.solve(vect);
943 
944 
945  // Now assign the discrete coefficients
946  // so that the Newmark scheme returns the correct
947  // results for the present values, and 1st and 2nd derivatives:
948 
949  // Value at present time
950  data_pt->set_value(0,j,U);
951 
952  // Value at previous time
953  data_pt->set_value(1,j,U0);
954 
955  // Veloc (time deriv) at previous time
956  data_pt->set_value(NSTEPS+1,j,vect[0]);
957 
958  // Acccel (2nd time deriv) at previous time
959  data_pt->set_value(NSTEPS+2,j,vect[1]);
960  }
961 
962 }
963 
964 
965 //=========================================================================
966 /// \short This function updates the Data's time history so that
967 /// we can advance to the next timestep.
968 //=========================================================================
969 template<unsigned NSTEPS>
971 {
972  //Find number of values stored
973  unsigned n_value = data_pt->nvalue();
974 
975  Vector<double> veloc(n_value);
976  time_derivative(1,data_pt,veloc);
977 
978  Vector<double> accel(n_value);
979  time_derivative(2,data_pt,accel);
980 
981  //Loop over values
982  for(unsigned j=0;j<n_value;j++)
983  {
984  //Set previous values/veloc/accel to present values/veloc/accel,
985  // if not a copy
986  if(data_pt->is_a_copy(j) == false)
987  {
988  for (unsigned t=NSTEPS;t>0;t--)
989  {
990  data_pt->set_value(t,j,data_pt->value(t-1,j));
991  }
992  data_pt->set_value(NSTEPS+1,j,veloc[j]);
993  data_pt->set_value(NSTEPS+2,j,accel[j]);
994  }
995  }
996 }
997 
998 //=========================================================================
999 /// \short This function updates a nodal time history so that
1000 /// we can advance to the next timestep.
1001 //=========================================================================
1002 template<unsigned NSTEPS>
1004 {
1005  //Find the number of coordinates
1006  unsigned n_dim = node_pt->ndim();
1007  //Find the number of position types
1008  unsigned n_position_type = node_pt->nposition_type();
1009 
1010  //Storage for the velocity and acceleration
1011  double veloc[n_position_type][n_dim];
1012  double accel[n_position_type][n_dim];
1013 
1014  //Find number of stored time values
1015  unsigned n_tstorage = ntstorage();
1016 
1017  //Loop over the variables
1018  for(unsigned i=0;i<n_dim;i++)
1019  {
1020  for(unsigned k=0;k<n_position_type;k++)
1021  {
1022  veloc[k][i] = accel[k][i] = 0.0;
1023  //Loop over all the history data
1024  for(unsigned t=0;t<n_tstorage;t++)
1025  {
1026  veloc[k][i] += weight(1,t)*node_pt->x_gen(t,k,i);
1027  accel[k][i] += weight(2,t)*node_pt->x_gen(t,k,i);
1028  }
1029  }
1030  }
1031 
1032  //Loop over the position variables
1033  for(unsigned i=0;i<n_dim;i++)
1034  {
1035  //If not a copy
1036  if(node_pt->position_is_a_copy(i) == false)
1037  {
1038  //Loop over position types
1039  for(unsigned k=0;k<n_position_type;k++)
1040  {
1041  //Set previous values/veloc/accel to present values/veloc/accel,
1042  // if not a copy
1043  for(unsigned t=NSTEPS;t>0;t--)
1044  {
1045  node_pt->x_gen(t,k,i) = node_pt->x_gen(t-1,k,i);
1046  }
1047 
1048  node_pt->x_gen(NSTEPS+1,k,i) = veloc[k][i];
1049  node_pt->x_gen(NSTEPS+2,k,i) = accel[k][i];
1050  }
1051  }
1052  }
1053 }
1054 
1055 //=========================================================================
1056 ///Set weights
1057 //=========================================================================
1058 template<unsigned NSTEPS>
1060 {
1061 
1062  double dt=Time_pt->dt(0);
1063 
1064  // Weights for second derivs
1065  Weight(2,0)= 2.0/(Beta2*dt*dt);
1066  Weight(2,1)=-2.0/(Beta2*dt*dt);
1067  for (unsigned t=2;t<=NSTEPS;t++)
1068  {
1069  Weight(2,t)=0.0;
1070  }
1071  Weight(2,NSTEPS+1)=-2.0/(dt*Beta2);
1072  Weight(2,NSTEPS+2)=(Beta2-1.0)/Beta2;
1073 
1074  // Weights for first derivs.
1075  Weight(1,0)=Beta1*dt*Weight(2,0);
1076  Weight(1,1)=Beta1*dt*Weight(2,1);
1077  for (unsigned t=2;t<=NSTEPS;t++)
1078  {
1079  Weight(1,t)=0.0;
1080  }
1081  Weight(1,NSTEPS+1)=1.0+Beta1*dt*Weight(2,NSTEPS+1);
1082  Weight(1,NSTEPS+2)=dt*(1.0-Beta1)+Beta1*dt*Weight(2,NSTEPS+2);
1083 }
1084 
1085 
1086 //===================================================================
1087 //Force build of templates: These are all the ones that might be
1088 //needed if Newmark is used together with BDF schemes (they require
1089 //up to 4 previous timesteps)
1090 //===================================================================
1091 template class Newmark<1>;
1092 template class Newmark<2>;
1093 template class Newmark<3>;
1094 template class Newmark<4>;
1095 
1096 
1097 
1098 //=========================================================================
1099 /// \short This function updates the Data's time history so that
1100 /// we can advance to the next timestep.
1101 //=========================================================================
1102 template<unsigned NSTEPS>
1104 {
1105  //Find number of values stored
1106  unsigned n_value = data_pt->nvalue();
1107 
1108  //Storage for the velocity and acceleration
1109  Vector<double> veloc(n_value,0.0);
1110  Vector<double> accel(n_value,0.0);
1111 
1112  //Find number of stored time values
1113  unsigned n_tstorage = this->ntstorage();
1114 
1115  //Loop over the variables
1116  for(unsigned i=0;i<n_value;i++)
1117  {
1118  //Loop over all the history data
1119  for(unsigned t=0;t<n_tstorage;t++)
1120  {
1121  veloc[i] += Newmark_veloc_weight[t]*data_pt->value(t,i);
1122  accel[i] += this->weight(2,t)*data_pt->value(t,i);
1123  }
1124  }
1125 
1126 
1127  //Loop over values
1128  for(unsigned j=0;j<n_value;j++)
1129  {
1130  //Set previous values/veloc/accel to present values/veloc/accel,
1131  // if not a copy
1132  if(data_pt->is_a_copy(j) == false)
1133  {
1134  for (unsigned t=NSTEPS;t>0;t--)
1135  {
1136  data_pt->set_value(t,j,data_pt->value(t-1,j));
1137  }
1138  data_pt->set_value(NSTEPS+1,j,veloc[j]);
1139  data_pt->set_value(NSTEPS+2,j,accel[j]);
1140  }
1141  }
1142 }
1143 
1144 //=========================================================================
1145 /// \short This function updates a nodal time history so that
1146 /// we can advance to the next timestep.
1147 //=========================================================================
1148 template<unsigned NSTEPS>
1150 {
1151  //Find the number of coordinates
1152  unsigned n_dim = node_pt->ndim();
1153  //Find the number of position types
1154  unsigned n_position_type = node_pt->nposition_type();
1155 
1156  //Storage for the velocity and acceleration
1157  double veloc[n_position_type][n_dim];
1158  double accel[n_position_type][n_dim];
1159 
1160  //Find number of stored time values
1161  unsigned n_tstorage =this->ntstorage();
1162 
1163  //Loop over the variables
1164  for(unsigned i=0;i<n_dim;i++)
1165  {
1166  for(unsigned k=0;k<n_position_type;k++)
1167  {
1168  veloc[k][i] = accel[k][i] = 0.0;
1169  //Loop over all the history data
1170  for(unsigned t=0;t<n_tstorage;t++)
1171  {
1172  veloc[k][i] += Newmark_veloc_weight[t]*node_pt->x_gen(t,k,i);
1173  accel[k][i] += this->weight(2,t)*node_pt->x_gen(t,k,i);
1174  }
1175  }
1176  }
1177 
1178  //Loop over the position variables
1179  for(unsigned i=0;i<n_dim;i++)
1180  {
1181  //If not a copy
1182  if(node_pt->position_is_a_copy(i) == false)
1183  {
1184  //Loop over position types
1185  for(unsigned k=0;k<n_position_type;k++)
1186  {
1187  //Set previous values/veloc/accel to present values/veloc/accel,
1188  // if not a copy
1189  for(unsigned t=NSTEPS;t>0;t--)
1190  {
1191  node_pt->x_gen(t,k,i) = node_pt->x_gen(t-1,k,i);
1192  }
1193 
1194  node_pt->x_gen(NSTEPS+1,k,i) = veloc[k][i];
1195  node_pt->x_gen(NSTEPS+2,k,i) = accel[k][i];
1196  }
1197  }
1198  }
1199 }
1200 
1201 
1202 //=========================================================================
1203 ///Set weights
1204 //=========================================================================
1205 template<>
1207 {
1208  // Set Newmark weights for second derivatives
1209  double dt=Time_pt->dt(0);
1210  Weight(2,0)= 2.0/(Beta2*dt*dt);
1211  Weight(2,1)=-2.0/(Beta2*dt*dt);
1212  Weight(2,2)=-2.0/(dt*Beta2);
1213  Weight(2,3)=(Beta2-1.0)/Beta2;
1214 
1215  // Set BDF weights for first derivatives
1216  Weight(1,0) = 1.0/dt;
1217  Weight(1,1) = -1.0/dt;
1218  Weight(1,2)=0.0;
1219  Weight(1,3)=0.0;
1220 
1221  // Orig Newmark weights for first derivs.
1222  set_newmark_veloc_weights(dt);
1223 }
1224 
1225 //=========================================================================
1226 ///Set weights
1227 //=========================================================================
1228 template<>
1230 {
1231  // Set Newmark weights for second derivatives
1232  double dt=Time_pt->dt(0);
1233  Weight(2,0)= 2.0/(Beta2*dt*dt);
1234  Weight(2,1)=-2.0/(Beta2*dt*dt);
1235  Weight(2,2)=0.0;
1236  Weight(2,3)=-2.0/(dt*Beta2);
1237  Weight(2,4)=(Beta2-1.0)/Beta2;
1238 
1239  // Set BDF weights for first derivatives
1240  if (Degrade_to_bdf1_for_first_derivs)
1241  {
1242  this->Weight(1,0) = 1.0/dt;
1243  this->Weight(1,1) = -1.0/dt;
1244  unsigned nweights=this->Weight.ncol();
1245  for (unsigned i=2;i<nweights;i++)
1246  {
1247  this->Weight(1,i)=0.0;
1248  }
1249  }
1250  else
1251  {
1252  double dtprev=Time_pt->dt(1);
1253  Weight(1,0) = 1.0/dt + 1.0/(dt + dtprev);
1254  Weight(1,1) = -(dt + dtprev)/(dt*dtprev);
1255  Weight(1,2) = dt/((dt+dtprev)*dtprev);
1256  Weight(1,3)=0.0;
1257  Weight(1,4)=0.0;
1258  }
1259 
1260  // Orig Newmark weights for first derivs.
1261  set_newmark_veloc_weights(dt);
1262 }
1263 
1264 //=========================================================================
1265 ///Set weights
1266 //=========================================================================
1267 template<>
1269 {
1270  // Set Newmark weights for second derivatives
1271  double dt=Time_pt->dt(0);
1272  Weight(2,0)= 2.0/(Beta2*dt*dt);
1273  Weight(2,1)=-2.0/(Beta2*dt*dt);
1274  Weight(2,2)=0.0;
1275  Weight(2,3)=0.0;
1276  Weight(2,4)=0.0;
1277  Weight(2,5)=-2.0/(dt*Beta2);
1278  Weight(2,6)=(Beta2-1.0)/Beta2;
1279 
1280  // Set BDF weights for first derivatives
1281 #ifdef PARANOID
1282  for (unsigned i=0;i<Time_pt->ndt();i++)
1283  {
1284  if (dt!=Time_pt->dt(i))
1285  {
1286  throw OomphLibError(
1287  "BDF4 currently only works for fixed timesteps \n",
1288  OOMPH_CURRENT_FUNCTION,
1289  OOMPH_EXCEPTION_LOCATION);
1290  }
1291  }
1292 #endif
1293 
1294  // Set BDF weights for first derivatives
1295  if (Degrade_to_bdf1_for_first_derivs)
1296  {
1297  this->Weight(1,0) = 1.0/dt;
1298  this->Weight(1,1) = -1.0/dt;
1299  unsigned nweights=this->Weight.ncol();
1300  for (unsigned i=2;i<nweights;i++)
1301  {
1302  this->Weight(1,i)=0.0;
1303  }
1304  }
1305  else
1306  {
1307  Weight(1,0) = 25.0/12.0/dt;
1308  Weight(1,1) = -48.0/12.0/dt;
1309  Weight(1,2) = 36.0/12.0/dt;
1310  Weight(1,3) = -16.0/12.0/dt;
1311  Weight(1,4) = 3.0/12.0/dt;
1312  Weight(1,5) = 0.0;
1313  Weight(1,6) = 0.0;
1314  }
1315 
1316  // Orig Newmark weights for first derivs.
1317  set_newmark_veloc_weights(dt);
1318 
1319 }
1320 
1321 //===================================================================
1322 //Force build of templates.
1323 //===================================================================
1324 template class NewmarkBDF<1>;
1325 template class NewmarkBDF<2>;
1326 template class NewmarkBDF<4>;
1327 
1328 
1329 }
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1234
void set_weights()
Set weights.
cstr elem_len * i
Definition: cfortran.h:607
Class to keep track of discrete/continous time. It is essential to have a single Time object when usi...
Definition: timesteppers.h:67
void assign_initial_positions_impulsive(Node *const &node_pt)
Initialise the time-history for the values, corresponding to an impulsive start.
virtual ~TimeStepper()
virtual destructor
Definition: timesteppers.cc:42
ExplicitTimeStepper * Explicit_predictor_pt
Definition: timesteppers.h:255
void set_predictor_weights()
Function to set the predictor weights.
void shift_time_positions(Node *const &node_pt)
This function updates a nodal time history so that we can advance to the next timestep.
char t
Definition: cfortran.h:572
void calculate_predicted_positions(Node *const &node_pt)
Function to calculate predicted positions at a node.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void assign_initial_data_values(Data *const &data_pt, Vector< InitialConditionFctPt > initial_value_fct, Vector< InitialConditionFctPt > initial_veloc_fct, Vector< InitialConditionFctPt > initial_accel_fct)
Initialise the time-history for the Data values, so that the Newmark representations for current velo...
static Time Dummy_time
Definition: timesteppers.h:842
void set_error_weights()
Function to set the error weights.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
Newmark scheme for second time deriv. Stored data represents.
Definition: timesteppers.h:868
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:966
void solve(DoubleVector &rhs)
Complete LU solve (replaces matrix by its LU decomposition and overwrites RHS with solution)...
Definition: matrices.cc:55
void assign_initial_values_impulsive(Data *const &data_pt)
Initialise the time-history for the values, corresponding to an impulsive start.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
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 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
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void set_weights()
Set the weights.
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i).
Definition: nodes.h:1055
double temporal_error_in_position(Node *const &node_pt, const unsigned &i)
Compute the error in the position i at a node.
void shift_time_values(Data *const &data_pt)
This function updates the Data's time history so that we can advance to the next timestep.
double temporal_error_in_value(Data *const &data_pt, const unsigned &i)
Compute the error in the value i in a Data structure.
void shift_time_values(Data *const &data_pt)
This function updates the Data's time history so that we can advance to the next timestep.
void assign_initial_data_values_stage1(const unsigned t_deriv, Data *const &data_pt)
First step in a two-stage procedure to assign the history values for the Newmark scheme so that the v...
void assign_initial_data_values_stage2(Data *const &data_pt)
Second step in a two-stage procedure to assign the history values for the Newmark scheme so that the ...
Newmark scheme for second time deriv with first derivatives calculated using BDF. ...
void calculate_predicted_values(Data *const &data_pt)
Function to calculate predicted data values in a Data object.
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
void shift_time_positions(Node *const &node_pt)
This function updates a nodal time history so that we can advance to the next timestep.
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
void set_weights()
Set weights.