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