timesteppers.h
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// This header contains classes and function prototypes for the Time,
27// TimeStepper and derived objects
28
29// Include guard to prevent multiple inclusions of the header
30#ifndef OOMPH_TIME_STEPPERS_HEADER
31#define OOMPH_TIME_STEPPERS_HEADER
32
33// Config header generated by autoconfig
34#ifdef HAVE_CONFIG_H
35#include <oomph-lib-config.h>
36#endif
37
38#ifdef OOMPH_HAS_MPI
39#include "mpi.h"
40#endif
41
42// oomph-lib headers
43#include "Vector.h"
44#include "nodes.h"
45#include "matrices.h"
46
47namespace oomph
48{
49 class Problem;
50 class ExplicitTimeStepper;
51
52 //====================================================================
53 /// Class to keep track of discrete/continous time. It is essential
54 /// to have a single Time object when using multiple time-stepping schemes;
55 /// e.g., in fluid-structure interaction problems, it is common to use
56 /// different schemes for the fluid and solid domains.
57 /// Storage is allocated for the current value
58 /// of the (continuous) time and a limited history of previous timesteps.
59 /// The number of previous timesteps must be equal to the number required
60 /// by the "highest order" scheme.
61 //====================================================================
62 class Time
63 {
64 private:
65 /// Pointer to the value of the continuous time.
67
68 /// Vector that stores the values of the current and previous timesteps.
70
71 public:
72 /// Constructor: Do not allocate any storage for previous timesteps,
73 /// but set the initial value of the time to zero
75
76 /// Constructor: Pass the number of timesteps to be stored
77 /// and set the initial value of time to zero.
78 Time(const unsigned& ndt) : Continuous_time(0.0)
79 {
80 // Allocate memory for the timestep storage and initialise values to one
81 // to avoid division by zero.
82 Dt.resize(ndt, 1.0);
83 }
84
85 /// Broken copy constructor
86 Time(const Time&) = delete;
87
88 /// Broken assignment operator
89 void operator=(const Time&) = delete;
90
91 /// Resize the vector holding the number of previous timesteps
92 /// and initialise the new values to zero.
93 void resize(const unsigned& n_dt)
94 {
95 Dt.resize(n_dt, 0.0);
96 }
97
98 /// Set all timesteps to the same value, dt.
99 void initialise_dt(const double& dt_)
100 {
101 unsigned ndt = Dt.size();
102 Dt.assign(ndt, dt_);
103 }
104
105 /// Set the value of the timesteps to be equal to the values passed
106 /// in a vector
108 {
109 // Assign the values from the vector dt_, but preserve the size of Dt and
110 // any Dt[i] s.t. i > dt_.size() (which is why we can't just use
111 // assignment operator).
112 unsigned n_dt = dt_.size();
113 for (unsigned i = 0; i < n_dt; i++)
114 {
115 Dt[i] = dt_[i];
116 }
117 }
118
119 /// Destructor: empty
120 ~Time() {}
121
122 /// Return the current value of the continuous time
123 double& time()
124 {
125 return Continuous_time;
126 }
127
128 /// Return the number of timesteps stored
129 unsigned ndt() const
130 {
131 return Dt.size();
132 }
133
134 /// Return the value of the t-th stored timestep (t=0: present; t>0:
135 /// previous).
136 double& dt(const unsigned& t = 0)
137 {
138 return Dt[t];
139 }
140
141 /// Return the value of the t-th stored timestep (t=0: present; t>0:
142 /// previous), const version.
143 double dt(const unsigned& t = 0) const
144 {
145 return Dt[t];
146 }
147
148 /// Return the value of the continuous time at the t-th previous
149 /// time level (t=0: current; t>0 previous).
150 double time(const unsigned& t = 0) const
151 {
152#ifdef PARANOID
153 if (t > ndt())
154 {
155 using namespace StringConversion;
156 std::string error_msg = "Timestep " + to_string(t) + " out of range!";
157 throw OomphLibError(
158 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
159 }
160#endif
161 // Load the current value of the time
162 double time_local = Continuous_time;
163 // Loop over the t previous timesteps and subtract each dt
164 for (unsigned i = 0; i < t; i++)
165 {
166 time_local -= Dt[i];
167 }
168 return time_local;
169 }
170
171 /// Update all stored values of dt by shifting each value along
172 /// the array. This function must be called before starting to solve at a
173 /// new time level.
174 void shift_dt()
175 {
176 unsigned n_dt = Dt.size();
177 // Return straight away if there are no stored timesteps
178 if (n_dt == 0)
179 {
180 return;
181 }
182 // Start from the end of the array and shift every entry back by one
183 for (unsigned i = (n_dt - 1); i > 0; i--)
184 {
185 Dt[i] = Dt[i - 1];
186 }
187 }
188 };
189
190
191 /// //////////////////////////////////////////////////////////////////////
192 /// //////////////////////////////////////////////////////////////////////
193 /// //////////////////////////////////////////////////////////////////////
194
195
196 //====================================================================
197 /// Base class for time-stepping schemes.
198 /// Timestepper provides an approximation of the temporal derivatives of
199 /// Data such that the i-th derivative of the j-th value in Data is
200 /// represented as
201 /// \code
202 /// unsigned n_value=data_pt->nvalue();
203 /// Vector<double> time_derivative(n_value);
204 /// for (unsigned j=0;j<n_value;j++)
205 /// {
206 /// time_derivative[j]=0.0;
207 /// for (unsigned t=0;t<ntstorage();t++)
208 /// {
209 /// time_derivative[j] += data_pt->value(t,j)*weight(i,t);
210 /// }
211 /// }
212 /// \endcode
213 /// where the time index \c t is such that
214 /// - \c t \c = \c 0 refers to the current value of the Data
215 /// - \c t \c > \c 0 refers to the `history' values, i.e. the doubles that
216 /// are stored in Data to represent the value's time derivatives. For BDF
217 /// schemes these doubles are the values at previous times;
218 /// for other timestepping schemes they can represent quantities
219 /// such as previous first and second derivatives, etc.
220 ///
221 /// TimeSteppers can evaluate all derivatives up to their \c order().
222 ///
223 /// The first \c nprev_values() `history values' represent the
224 /// values at the previous timesteps. For BDF schemes,
225 /// \c nprev_values() \c = \c ntstorage()-1 , i.e. all `history values'
226 /// represent previous values. For \c Newmark<NSTEPS>, only the first
227 /// NSTEPS `history values' represent previous values (NSTEPS=1
228 /// gives the normal Newmark scheme).
229 //====================================================================
231 {
232 protected:
233 /// Pointer to discrete time storage scheme
235
236 /// Storage for the weights associated with the timestepper
238
239 /// String that indicates the type of the timestepper
240 /// (e.g. "BDF", "Newmark", etc.)
242
243 /// Boolean variable to indicate whether the timestepping scheme can
244 /// be adaptive
246
247 /// Bool to indicate if the timestepper is steady, i.e. its
248 /// time-derivatives evaluate to zero. This status may be achieved
249 /// temporarily by calling make_steady(). It can be reset to the
250 /// appropriate default by the function undo_make_steady().
252
253 /// Boolean to indicate if the timestepper will output warnings when
254 /// setting possibly an incorrect number of initial data values from
255 /// function pointers
257
258 /// Flag: is adaptivity done by taking a separate step using an
259 /// ExplicitTimeStepper object?
261
262 /// Pointer to explicit time stepper to use as predictor if
263 /// Predict_by_explicit_step is set.
264 /// ??ds merge the two? predict by explicit if pointer is set?
266
267 /// Store the time that the predicted values currently stored are at, to
268 /// compare for paranoid checks.
270
271 /// The time-index in each Data object where predicted values are
272 /// stored. -1 if not set.
274
275 public:
276 /// Constructor. Pass the amount of storage required by
277 /// timestepper (present value + history values) and the
278 /// order of highest time-derivative.
279 TimeStepper(const unsigned& tstorage, const unsigned& max_deriv)
280 : Time_pt(0),
281 Adaptive_Flag(false),
282 Is_steady(false),
285 {
286 // Resize Weights matrix and initialise each weight to zero
287 Weight.resize(max_deriv + 1, tstorage, 0.0);
288
289 // Set the weight for zero-th derivative, which is always 1.0
290 Weight(0, 0) = 1.0;
291
292 // Set predictor storage index to negative value so that we can catch
293 // cases where it has not been set to a correct value by the inheriting
294 // constructor.
296
298 }
299
300 /// Broken empty constructor
302 {
303 throw OomphLibError("Don't call default constructor for TimeStepper!",
304 OOMPH_CURRENT_FUNCTION,
305 OOMPH_EXCEPTION_LOCATION);
306 }
307
308 /// Broken copy constructor
309 TimeStepper(const TimeStepper&) = delete;
310
311 /// Broken assignment operator
312 void operator=(const TimeStepper&) = delete;
313
314 /// virtual destructor
315 virtual ~TimeStepper();
316
317 /// Highest order derivative that the scheme can compute
318 unsigned highest_derivative() const
319 {
320 return Weight.nrow() - 1;
321 }
322
323 /// Actual order (accuracy) of the scheme
324 virtual unsigned order() const
325 {
326 return 0;
327 }
328
329 /// Return current value of continous time
330 // (can't have a paranoid test for null pointers because this could be
331 // used as a set function)
332 double& time()
333 {
334 return Time_pt->time();
335 }
336
337 /// Return current value of continous time.
338 double time() const
339 {
340#ifdef PARANOID
341 if (Time_pt == 0)
342 {
343 std::string error_msg = "Time pointer is null!";
344 throw OomphLibError(
345 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
346 }
347#endif
348 return Time_pt->time();
349 }
350
351 /// Number of timestep increments that are required by the scheme
352 virtual unsigned ndt() const = 0;
353
354 /// Number of previous values needed to calculate the value at the
355 /// current time. i.e. how many previous values must we loop over to
356 /// calculate the values at the evaluation time. For most methods this is
357 /// 1, i.e. just use the value at t_{n+1}. See midpoint method for a
358 /// counter-example.
360 {
361 return 1;
362 }
363
364 /// Number of previous values available: 0 for static, 1 for BDF<1>,...
365 virtual unsigned nprev_values() const = 0;
366
367 /// Function to set the weights for present timestep (don't
368 /// need to pass present timestep or previous timesteps as they
369 /// are available via Time_pt)
370 virtual void set_weights() = 0;
371
372 /// Function to make the time stepper temporarily steady. This
373 /// is trivially achieved by setting all the weights to zero
375 {
376 // Zero the matrix
378
379 // Weight of current step in calculation of current step is always 1 in
380 // steady state. All other entries left at zero.
381 Weight(0, 0) = 1.0;
382
383 // Update flag
384 Is_steady = true;
385 }
386
387 /// Flag to indicate if a timestepper has been made steady (possibly
388 /// temporarily to switch off time-dependence)
389 bool is_steady() const
390 {
391 return Is_steady;
392 }
393
394 /// Flag: is adaptivity done by taking a separate step using an
395 /// ExplicitTimeStepper object?
397 {
399 }
400
401 /// Get the pointer to the explicit timestepper to use as a predictor in
402 /// adaptivity if Predict_by_explicit_step is set.
404 {
406 }
407
408 /// Set the pointer to the explicit timestepper to use as a predictor in
409 /// adaptivity if Predict_by_explicit_step is set.
411 {
412 Explicit_predictor_pt = _pred_pt;
413 }
414
415 /// Set the time that the current predictions apply for, only needed for
416 /// paranoid checks when doing Predict_by_explicit_step.
417 void update_predicted_time(const double& new_time)
418 {
419 Predicted_time = new_time;
420 }
421
422 /// Check that the predicted values are the ones we want.
424 {
425#ifdef PARANOID
426 if (std::abs(time_pt()->time() - Predicted_time) > 1e-15)
427 {
428 throw OomphLibError(
429 "Predicted values are not from the correct time step",
430 OOMPH_EXCEPTION_LOCATION,
431 OOMPH_CURRENT_FUNCTION);
432 }
433#endif
434 }
435
436 /// Return the time-index in each Data where predicted values are
437 /// stored if the timestepper is adaptive.
438 unsigned predictor_storage_index() const
439 {
440 // Error if time stepper is non-adaptive (because then the index doesn't
441 // exist so using it would give a potentially difficult to find
442 // segault).
443#ifdef PARANOID
445 {
446 return unsigned(Predictor_storage_index);
447 }
448 else
449 {
450 std::string err = "Predictor storage index is negative, this probably";
451 err += " means it hasn't been set for this timestepper.";
452 throw OomphLibError(
453 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
454 }
455#else
456 return unsigned(Predictor_storage_index);
457#endif
458 }
459
460 /// Enable the output of warnings due to possible fct pointer vector
461 /// size mismatch in assign_initial_data_values (Default)
463 {
465 }
466
467 /// Disable the output of warnings due to possible fct pointer vector
468 /// size mismatch in assign_initial_data_values
470 {
472 }
473
474 /// Get a (const) pointer to the weights.
476 {
477 return &Weight;
478 }
479
480 /// Reset the is_steady status of a specific TimeStepper to its
481 /// default and re-assign the weights.
482 virtual void undo_make_steady()
483 {
484 Is_steady = false;
485 set_weights();
486 }
487
488 /// Return string that indicates the type of the timestepper
489 /// (e.g. "BDF", "Newmark", etc.)
491 {
492 return Type;
493 }
494
495 //??ds I think that the Data and Node cases below couldp robably be
496 // handled with a template argument. However they can't be handled by
497 // normal polymorphism because data.value is different to node.value and
498 // value is not a virtual function.
499
500 /// Evaluate i-th derivative of all values in Data and return in
501 /// Vector deriv[].
502 void time_derivative(const unsigned& i,
503 Data* const& data_pt,
504 Vector<double>& deriv)
505 {
506 // Number of values stored in the Data object
507 unsigned nvalue = data_pt->nvalue();
508 deriv.assign(nvalue, 0.0);
509
510 // Loop over all values
511 for (unsigned j = 0; j < nvalue; j++)
512 {
513 deriv[j] = time_derivative(i, data_pt, j);
514 }
515 }
516
517 /// Evaluate i-th derivative of j-th value in Data.
518 double time_derivative(const unsigned& i,
519 Data* const& data_pt,
520 const unsigned& j)
521 {
522 double deriv = 0.0;
523 unsigned n_tstorage = ntstorage();
524 // Loop over all history data and add the appropriate contribution
525 // to the derivative
526 for (unsigned t = 0; t < n_tstorage; t++)
527 {
528 deriv += Weight(i, t) * data_pt->value(t, j);
529 }
530 return deriv;
531 }
532
533 /// Evaluate i-th derivative of all values in Node and return in
534 /// Vector deriv[] (this can't be simply combined with time_derivative(..,
535 /// Data, ...) because of differences with haning nodes).
536 void time_derivative(const unsigned& i,
537 Node* const& node_pt,
538 Vector<double>& deriv)
539 {
540 // Number of values stored in the Data object
541 unsigned nvalue = node_pt->nvalue();
542 deriv.assign(nvalue, 0.0);
543
544 // Loop over all values
545 for (unsigned j = 0; j < nvalue; j++)
546 {
547 deriv[j] = time_derivative(i, node_pt, j);
548 }
549 }
550
551 /// Evaluate i-th derivative of j-th value in Node. Note the use of
552 /// the node's value() function so that hanging nodes are taken into
553 /// account (this is why the functions for Data and Node cannot be
554 /// combined through simple polymorphism: value is not virtual).
555 double time_derivative(const unsigned& i,
556 Node* const& node_pt,
557 const unsigned& j)
558 {
559 double deriv = 0.0;
560 unsigned n_tstorage = ntstorage();
561 // Loop over all history data and add the appropriate contribution
562 // to the derivative
563 for (unsigned t = 0; t < n_tstorage; t++)
564 {
565 deriv += weight(i, t) * node_pt->value(t, j);
566 }
567 return deriv;
568 }
569
570
571 /// Access function for the pointer to time (const version)
572 Time* const& time_pt() const
573 {
574#ifdef PARANOID
575 if (Time_pt == 0)
576 {
577 std::string error_msg(
578 "Time_pt is null, probably because it is a steady time stepper.");
579 throw OomphLibError(
580 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
581 }
582#endif
583 return Time_pt;
584 }
585
586 /// Access function for the pointer to time (can't have a paranoid test
587 /// for null pointers because this is used as a set function).
589 {
590 return Time_pt;
591 }
592
593 /// Access function for j-th weight for the i-th derivative
594 virtual double weight(const unsigned& i, const unsigned& j) const
595 {
596 return Weight(i, j);
597 }
598
599 /// Return the number of doubles required to represent history
600 /// (one for steady)
601 unsigned ntstorage() const
602 {
603 return (Weight.ncol());
604 }
605
606 /// Initialise the time-history for the Data values
607 /// corresponding to an impulsive start.
608 virtual void assign_initial_values_impulsive(Data* const& data_pt) = 0;
609
610 /// Initialiset the positions for the node corresponding to
611 /// an impulsive start
612 virtual void assign_initial_positions_impulsive(Node* const& node_pt) = 0;
613
614 /// This function advances the Data's time history so that
615 /// we can move on to the next timestep
616 virtual void shift_time_values(Data* const& data_pt) = 0;
617
618 /// This function advances the time history of the positions
619 /// at a node. The default should be OK, but would need to be overloaded
620 virtual void shift_time_positions(Node* const& node_pt) = 0;
621
622 /// Function to indicate whether the scheme is adaptive (false by default)
623 bool adaptive_flag() const
624 {
625 return Adaptive_Flag;
626 }
627
628 /// Set the weights for the predictor
629 /// previous timestep (currently empty -- overwrite for specific scheme)
630 virtual void set_predictor_weights() {}
631
632 /// Do the predictor step for data stored in a Data object
633 /// (currently empty -- overwrite for specific scheme)
634 virtual void calculate_predicted_values(Data* const& data_pt) {}
635
636 /// Do the predictor step for the positions at a node
637 /// (currently empty --- overwrite for a specific scheme)
638 virtual void calculate_predicted_positions(Node* const& node_pt) {}
639
640 /// Set the weights for the error computation,
641 /// (currently empty -- overwrite for specific scheme)
642 virtual void set_error_weights() {}
643
644 /// Compute the error in the position i at a node
645 /// zero here -- overwrite for specific scheme.
646 virtual double temporal_error_in_position(Node* const& node_pt,
647 const unsigned& i)
648 {
649 return 0.0;
650 }
651
652 /// Compute the error in the value i in a Data structure
653 /// zero here -- overwrite for specific scheme.
654 virtual double temporal_error_in_value(Data* const& data_pt,
655 const unsigned& i)
656 {
657 return 0.0;
658 }
659
660 /// Interface for any actions that need to be performed before a time
661 /// step.
662 virtual void actions_before_timestep(Problem* problem_pt) {}
663
664 /// Interface for any actions that need to be performed after a time
665 /// step.
666 virtual void actions_after_timestep(Problem* problem_pt) {}
667 };
668
669
670 /// //////////////////////////////////////////////////////////////////////
671 /// //////////////////////////////////////////////////////////////////////
672 /// //////////////////////////////////////////////////////////////////////
673
674
675 //====================================================================
676 /// Faux time-stepper for steady problems. Allows storage
677 /// for NSTEPS previous values.
678 //====================================================================
679 template<unsigned NSTEPS>
680 class Steady : virtual public TimeStepper
681 {
682 public:
683 /// Constructor: Creates storage for NSTEPS previous timesteps
684 /// and can evaluate up to 2nd derivatives (though it doesn't
685 /// actually do anything -- all time-derivatives
686 /// evaluate to zero)
687 Steady() : TimeStepper(NSTEPS + 1, 2)
688 {
689 Type = "Steady";
691
692 // Overwrite default assignment in base constructor -- this TimeStepper
693 // actually is steady all the time.
694 Is_steady = true;
695 }
696
697
698 /// Broken copy constructor
699 Steady(const Steady&) = delete;
700
701 /// Broken assignment operator
702 void operator=(const Steady&) = delete;
703
704 /// Return the actual order of the scheme. Returning zero here --
705 /// doesn't make much sense, though
706 unsigned order() const
707 {
708 return 0;
709 }
710
711 /// Initialise the time-history for the Data values,
712 /// corresponding to an impulsive start.
714 {
715 // Find number of values stored
716 unsigned n_value = data_pt->nvalue();
717 // Loop over values
718 for (unsigned j = 0; j < n_value; j++)
719 {
720 // Set previous values to the initial value, if not a copy
721 if (data_pt->is_a_copy(j) == false)
722 {
723 for (unsigned t = 1; t <= NSTEPS; t++)
724 {
725 data_pt->set_value(t, j, data_pt->value(j));
726 }
727 }
728 }
729 }
730
731 /// Initialise the time-history for the nodal positions
732 /// corresponding to an impulsive start.
734 {
735 // Find the number of coordinates
736 unsigned n_dim = node_pt->ndim();
737 // Find the number of positoin types
738 unsigned n_position_type = node_pt->nposition_type();
739
740 // Loop over values
741 for (unsigned i = 0; i < n_dim; i++)
742 {
743 // Set previous values to the initial value, if not a copy
744 if (node_pt->position_is_a_copy(i) == false)
745 {
746 for (unsigned k = 0; k < n_position_type; k++)
747 {
748 for (unsigned t = 1; t <= NSTEPS; t++)
749 {
750 node_pt->x_gen(t, k, i) = node_pt->x_gen(k, i);
751 }
752 }
753 }
754 }
755 }
756
757
758 /// Typedef for function that returns the (scalar) initial
759 /// value at a given value of the continuous time t.
760 typedef double (*InitialConditionFctPt)(const double& t);
761
762 /// Initialise the time-history for the Data values,
763 /// corresponding to given time history, specified by
764 /// Vector of function pointers.
766 Data* const& data_pt, Vector<InitialConditionFctPt> initial_value_fct)
767 {
768 // The time history stores the previous function values
769 unsigned n_time_value = ntstorage();
770
771 // Find number of values stored
772 unsigned n_value = data_pt->nvalue();
773
774 // Loop over current and stored timesteps
775 for (unsigned t = 0; t < n_time_value; t++)
776 {
777 // Get corresponding continous time
778 double time_local = Time_pt->time(t);
779
780 // Loop over values
781 for (unsigned j = 0; j < n_value; j++)
782 {
783 data_pt->set_value(t, j, initial_value_fct[j](time_local));
784 }
785 }
786 }
787
788 /// This function updates the Data's time history so that
789 /// we can advance to the next timestep. As for BDF schemes,
790 /// we simply push the values backwards...
791 void shift_time_values(Data* const& data_pt)
792 {
793 // Find number of values stored
794 unsigned n_value = data_pt->nvalue();
795
796 // Loop over the values
797 for (unsigned j = 0; j < n_value; j++)
798 {
799 // Set previous values to the previous value, if not a copy
800 if (data_pt->is_a_copy(j) == false)
801 {
802 // Loop over times, in reverse order
803 for (unsigned t = NSTEPS; t > 0; t--)
804 {
805 data_pt->set_value(t, j, data_pt->value(t - 1, j));
806 }
807 }
808 }
809 }
810
811 /// This function advances the time history of the positions
812 /// at a node.
813 void shift_time_positions(Node* const& node_pt)
814 {
815 // Find the number of coordinates
816 unsigned n_dim = node_pt->ndim();
817 // Find the number of position types
818 unsigned n_position_type = node_pt->nposition_type();
819
820 // Loop over the positions
821 for (unsigned i = 0; i < n_dim; i++)
822 {
823 // If the position is not a copy
824 if (node_pt->position_is_a_copy(i) == false)
825 {
826 for (unsigned k = 0; k < n_position_type; k++)
827 {
828 // Loop over stored times, and set values to previous values
829 for (unsigned t = NSTEPS; t > 0; t--)
830 {
831 node_pt->x_gen(t, k, i) = node_pt->x_gen(t - 1, k, i);
832 }
833 }
834 }
835 }
836 }
837
838 /// Set weights
840 {
841 // Loop over higher derivatives
842 unsigned max_deriv = highest_derivative();
843 for (unsigned i = 0; i < max_deriv; i++)
844 {
845 for (unsigned j = 0; j < NSTEPS; j++)
846 {
847 Weight(i, j) = 0.0;
848 }
849 }
850
851 // Zeroth derivative must return the value itself:
852 Weight(0, 0) = 1.0;
853 }
854
855 /// Number of previous values available.
856 unsigned nprev_values() const
857 {
858 return NSTEPS;
859 }
860
861 /// Number of timestep increments that need to be stored by the scheme
862 unsigned ndt() const
863 {
864 return NSTEPS;
865 }
866
867 /// Dummy: Access function for j-th weight for the i-th derivative
868 double weight(const unsigned& i, const unsigned& j) const
869 {
870 if ((i == 0) && (j == 0))
871 {
872 return One;
873 }
874 else
875 {
876 return Zero;
877 }
878 }
879
880 private:
881 /// Static variable to hold the value 1.0
882 static double One;
883
884 /// Static variable to hold the value 0.0
885 static double Zero;
886
887 // Default Time object
889 };
890
891
892 /// //////////////////////////////////////////////////////////////////////
893 /// //////////////////////////////////////////////////////////////////////
894 /// //////////////////////////////////////////////////////////////////////
895
896
897 //====================================================================
898 /// Newmark scheme for second time deriv. Stored data represents
899 /// - t=0: value at at present time, Time_pt->time()
900 /// - t=1: value at previous time, Time_pt->time()-dt
901 /// - ...
902 /// - t=NSTEPS: value at previous time, Time_pt->time()-NSTEPS*dt
903 /// - t=NSTEPS+1: 1st time deriv (= "velocity") at previous time,
904 /// Time_pt->time()-dt
905 /// - t=NSTEPS+2: 2nd time deriv (= "acceleration") at previous time,
906 /// Time_pt->time()-dt
907 ///
908 /// NSTEPS=1 gives normal Newmark.
909 //====================================================================
910 template<unsigned NSTEPS>
911 class Newmark : public TimeStepper
912 {
913 public:
914 /// Constructor: Pass pointer to global time. We set up a
915 /// timestepping scheme with NSTEPS+2 doubles to represent the history and
916 /// the highest deriv is 2.
917 Newmark() : TimeStepper(NSTEPS + 3, 2)
918 {
919 Type = "Newmark";
920
921 // Standard Newmark parameters
922 Beta1 = 0.5;
923 Beta2 = 0.5;
924 }
925
926 /// Broken copy constructor
927 Newmark(const Newmark&) = delete;
928
929 /// Broken assignment operator
930 void operator=(const Newmark&) = delete;
931
932
933 /// The actual order (accuracy of the scheme)
934 unsigned order() const
935 {
936 std::string error_message =
937 "Can't remember the order of the Newmark scheme";
938 error_message += " -- I think it's 2nd order...\n";
939
941 error_message, "Newmark::order()", OOMPH_EXCEPTION_LOCATION);
942 return 2;
943 }
944
945 /// Initialise the time-history for the values,
946 /// corresponding to an impulsive start.
947 void assign_initial_values_impulsive(Data* const& data_pt);
948
949 /// Initialise the time-history for the values,
950 /// corresponding to an impulsive start.
951 void assign_initial_positions_impulsive(Node* const& node_pt);
952
953 /// Typedef for function that returns the (scalar) initial
954 /// value at a given value of the continuous time t.
955 typedef double (*InitialConditionFctPt)(const double& t);
956
957 /// Initialise the time-history for the Data values,
958 /// so that the Newmark representations for current veloc and
959 /// acceleration are exact.
961 Data* const& data_pt,
962 Vector<InitialConditionFctPt> initial_value_fct,
963 Vector<InitialConditionFctPt> initial_veloc_fct,
964 Vector<InitialConditionFctPt> initial_accel_fct);
965
966
967 /// Typedef for function that returns the (scalar) initial
968 /// value at a given value of the continuous time t and the spatial
969 /// coordinate -- appropriate for assignement of initial conditions for
970 /// nodes
971 typedef double (*NodeInitialConditionFctPt)(const double& t,
972 const Vector<double>& x);
973
974 /// Initialise the time-history for the nodal values,
975 /// so that the Newmark representations for current veloc and
976 /// acceleration are exact.
978 Node* const& node_pt,
979 Vector<NodeInitialConditionFctPt> initial_value_fct,
980 Vector<NodeInitialConditionFctPt> initial_veloc_fct,
981 Vector<NodeInitialConditionFctPt> initial_accel_fct);
982
983 /// First step in a two-stage procedure to assign
984 /// the history values for the Newmark scheme so
985 /// that the veloc and accel that are computed by the scheme
986 /// are correct at the current time.
987 ///
988 /// Call this function for t_deriv=0,1,2,3.
989 /// When calling with
990 /// - t_deriv=0 : data_pt(0,*) should contain the values at the
991 /// previous timestep
992 /// - t_deriv=1 : data_pt(0,*) should contain the current values;
993 /// they get stored (temporarily) in data_pt(1,*).
994 /// - t_deriv=2 : data_pt(0,*) should contain the current velocities
995 /// (first time derivs); they get stored (temporarily)
996 /// in data_pt(NSTEPS+1,*).
997 /// - t_deriv=3 : data_pt(0,*) should contain the current accelerations
998 /// (second time derivs); they get stored (temporarily)
999 /// in data_pt(NSTEPS+2,*).
1000 /// .
1001 /// Follow this by calls to
1002 /// \code
1003 /// assign_initial_data_values_stage2(...)
1004 /// \endcode
1005 void assign_initial_data_values_stage1(const unsigned t_deriv,
1006 Data* const& data_pt);
1007
1008 /// Second step in a two-stage procedure to assign
1009 /// the history values for the Newmark scheme so
1010 /// that the veloc and accel that are computed by the scheme
1011 /// are correct at the current time.
1012 ///
1013 /// This assigns appropriate values for the "previous
1014 /// velocities and accelerations" so that their current
1015 /// values, which were defined in assign_initial_data_values_stage1(...),
1016 /// are represented exactly by the Newmark scheme.
1017 void assign_initial_data_values_stage2(Data* const& data_pt);
1018
1019
1020 /// This function updates the Data's time history so that
1021 /// we can advance to the next timestep.
1022 void shift_time_values(Data* const& data_pt);
1023
1024 /// This function updates a nodal time history so that
1025 /// we can advance to the next timestep.
1026 void shift_time_positions(Node* const& node_pt);
1027
1028 /// Set weights
1029 void set_weights();
1030
1031 /// Number of previous values available.
1032 unsigned nprev_values() const
1033 {
1034 return NSTEPS;
1035 }
1036
1037 /// Number of timestep increments that need to be stored by the scheme
1038 unsigned ndt() const
1039 {
1040 return NSTEPS;
1041 }
1042
1043
1044 protected:
1045 /// First Newmark parameter (usually 0.5)
1046 double Beta1;
1047
1048 /// Second Newmark parameter (usually 0.5)
1049 double Beta2;
1050 };
1051
1052
1053 /// //////////////////////////////////////////////////////////////////////
1054 /// //////////////////////////////////////////////////////////////////////
1055 /// //////////////////////////////////////////////////////////////////////
1056
1057
1058 //====================================================================
1059 /// Newmark scheme for second time deriv with first derivatives
1060 /// calculated using BDF. . Stored data represents
1061 /// - t=0: value at at present time, Time_pt->time()
1062 /// - t=1: value at previous time, Time_pt->time()-dt
1063 /// - ...
1064 /// - t=NSTEPS: value at previous time, Time_pt->time()-NSTEPS*dt
1065 /// - t=NSTEPS+1: 1st time deriv (= "velocity") at previous time,
1066 /// Time_pt->time()-dt
1067 /// - t=NSTEPS+2: 2nd time deriv (= "acceleration") at previous time,
1068 /// Time_pt->time()-dt
1069 ///
1070 /// NSTEPS=1 gives normal Newmark.
1071 //====================================================================
1072 template<unsigned NSTEPS>
1073 class NewmarkBDF : public Newmark<NSTEPS>
1074 {
1075 public:
1076 /// Constructor: Pass pointer to global time. We set up a
1077 /// timestepping scheme with NSTEPS+2 doubles to represent the history and
1078 /// the highest deriv is 2.
1080 {
1081 this->Type = "NewmarkBDF";
1083 Newmark_veloc_weight.resize(NSTEPS + 3);
1084 }
1085
1086 /// Broken copy constructor
1087 NewmarkBDF(const NewmarkBDF&) = delete;
1088
1089 /// Broken assignment operator
1090 void operator=(const NewmarkBDF&) = delete;
1091
1092 /// Set weights
1094
1095 /// This function updates the Data's time history so that
1096 /// we can advance to the next timestep.
1097 void shift_time_values(Data* const& data_pt);
1098
1099 /// This function updates a nodal time history so that
1100 /// we can advance to the next timestep.
1101 void shift_time_positions(Node* const& node_pt);
1102
1103 /// Degrade scheme to first order BDF (for first derivs/veloc);
1104 /// usually for start-up.
1106 {
1108 }
1109
1110
1111 /// Disable degradation to first order BDF.
1113 {
1115 }
1116
1117 private:
1118 /// Set original Newmark weights for velocities (needed when
1119 /// shifting history values -- they're used when updating the
1120 /// previous accelerations and doing this with bdf can make the
1121 /// scheme unstable...
1122 void set_newmark_veloc_weights(const double& dt)
1123 {
1124 Newmark_veloc_weight[0] = this->Beta1 * dt * this->Weight(2, 0);
1125 Newmark_veloc_weight[1] = this->Beta1 * dt * this->Weight(2, 1);
1126 for (unsigned t = 2; t <= NSTEPS; t++)
1127 {
1128 Newmark_veloc_weight[t] = 0.0;
1129 }
1130 Newmark_veloc_weight[NSTEPS + 1] =
1131 1.0 + this->Beta1 * dt * this->Weight(2, NSTEPS + 1);
1132 Newmark_veloc_weight[NSTEPS + 2] =
1133 dt * (1.0 - this->Beta1) +
1134 this->Beta1 * dt * this->Weight(2, NSTEPS + 2);
1135 }
1136
1137 /// Boolean flag to indicate degradation of scheme to first
1138 /// order BDF (for first derivs/veloc); usually for start-up.
1140
1141 /// Original Newmark weights for velocities (needed when
1142 /// shifting history values -- they're used when updating the
1143 /// previous accelerations and doing this with bdf can make the
1144 /// scheme unstable...
1146 };
1147
1148
1149 /// //////////////////////////////////////////////////////////////////////
1150 /// //////////////////////////////////////////////////////////////////////
1151 /// //////////////////////////////////////////////////////////////////////
1152
1153
1154 //====================================================================
1155 /// Templated class for BDF-type time-steppers with fixed or
1156 /// variable timestep.
1157 /// 1st time derivative recovered directly from the previous function
1158 /// values. Template parameter represents the number of previous timesteps
1159 /// stored, so that BDF<1> is the classical first order
1160 /// backward Euler scheme.
1161 /// Need to reset weights after every change in timestep.
1162 //====================================================================
1163 template<unsigned NSTEPS>
1164 class BDF : public TimeStepper
1165 {
1166 // A BDF<1> time data set consists of:
1167 // [y_np1, y_n, dy_n, y^P_np1]
1168 // Or in english:
1169 // * y-value at time being/just been solved for
1170 // * y-value at previous time
1171 // * approximation to y derivative at previous time (also refered to as
1172 // velocity in some places, presumably it corresponds to velocity in
1173 // solid mechanics).
1174 // * predicted y-value at time n+1
1175
1176 // A BDF<2> time data set consists of:
1177 // [y_np1, y_n, y_nm1, dy_n, y^P_np1]
1178 // i.e. the same thing but with one more previous time value. Also the
1179 // derivative approximation will be more accurate.
1180
1181 // If the adaptive flag is set to false then the final two pieces of data
1182 // in each are not stored (derivative and predictor value).
1183
1184 private:
1185 /// Private data for the predictor weights
1187
1188 /// Private data for the error weight
1190
1191 public:
1192 /// Constructor for the case when we allow adaptive timestepping
1193 BDF(const bool& adaptive = false) : TimeStepper(NSTEPS + 1, 1)
1194 {
1195 Type = "BDF";
1196
1197 // If it's adaptive, we need to allocate additional space to
1198 // carry along a prediction and an acceleration
1199 if (adaptive)
1200 {
1201 // Set the adaptive flag to be true
1202 Adaptive_Flag = true;
1203
1204 // Set the size of the Predictor_Weights Vector N.B. The size is
1205 // correct for BDF1 and 2, but may be wrong for others.
1206 Predictor_weight.resize(NSTEPS + 2);
1207
1208 // Resize the weights to the appropriate size
1209 Weight.resize(2, NSTEPS + 3, 0.0);
1210
1211 // Storing predicted values in slot after other information
1212 Predictor_storage_index = NSTEPS + 2;
1213 }
1214
1215 // Set the weight for the zero-th derivative
1216 Weight(0, 0) = 1.0;
1217 }
1218
1219
1220 /// Broken copy constructor
1221 BDF(const BDF&) = delete;
1222
1223 /// Broken assignment operator
1224 void operator=(const BDF&) = delete;
1225
1226 /// Return the actual order of the scheme
1227 unsigned order() const
1228 {
1229 return NSTEPS;
1230 }
1231
1232 /// Initialise the time-history for the Data values,
1233 /// corresponding to an impulsive start.
1235 {
1236 // Find number of values stored
1237 unsigned n_value = data_pt->nvalue();
1238 // Loop over values
1239 for (unsigned j = 0; j < n_value; j++)
1240 {
1241 // Set previous values to the initial value, if not a copy
1242 if (data_pt->is_a_copy(j) == false)
1243 {
1244 for (unsigned t = 1; t <= NSTEPS; t++)
1245 {
1246 data_pt->set_value(t, j, data_pt->value(j));
1247 }
1248
1249 // If it's adaptive
1250 if (adaptive_flag())
1251 {
1252 // Initial velocity is zero
1253 data_pt->set_value(NSTEPS + 1, j, 0.0);
1254 // Initial prediction is the value
1255 data_pt->set_value(NSTEPS + 2, j, data_pt->value(j));
1256 }
1257 }
1258 }
1259 }
1260
1261 /// Initialise the time-history for the nodal positions
1262 /// corresponding to an impulsive start.
1264 {
1265 // Find the dimension of the node
1266 unsigned n_dim = node_pt->ndim();
1267 // Find the number of position types at the node
1268 unsigned n_position_type = node_pt->nposition_type();
1269
1270 // Loop over the position variables
1271 for (unsigned i = 0; i < n_dim; i++)
1272 {
1273 // If the position is not copied
1274 // We copy entire coordinates at once
1275 if (node_pt->position_is_a_copy(i) == false)
1276 {
1277 // Loop over the position types
1278 for (unsigned k = 0; k < n_position_type; k++)
1279 {
1280 // Set previous values to the initial value, if not a copy
1281 for (unsigned t = 1; t <= NSTEPS; t++)
1282 {
1283 node_pt->x_gen(t, k, i) = node_pt->x_gen(k, i);
1284 }
1285
1286 // If it's adaptive
1287 if (adaptive_flag())
1288 {
1289 // Initial mesh velocity is zero
1290 node_pt->x_gen(NSTEPS + 1, k, i) = 0.0;
1291 // Initial prediction is the value
1292 node_pt->x_gen(NSTEPS + 2, k, i) = node_pt->x_gen(k, i);
1293 }
1294 }
1295 }
1296 }
1297 }
1298
1299
1300 /// Typedef for function that returns the (scalar) initial
1301 /// value at a given value of the continuous time t.
1302 typedef double (*InitialConditionFctPt)(const double& t);
1303
1304 /// Initialise the time-history for the Data values,
1305 /// corresponding to given time history, specified by
1306 /// Vector of function pointers.
1308 Data* const& data_pt, Vector<InitialConditionFctPt> initial_value_fct)
1309 {
1310 // The time history stores the previous function values
1311 unsigned n_time_value = ntstorage();
1312
1313 // Find number of values stored
1314 unsigned n_value = data_pt->nvalue();
1315
1316 // Loop over current and stored timesteps
1317 for (unsigned t = 0; t < n_time_value; t++)
1318 {
1319 // Get corresponding continous time
1320 double time_local = Time_pt->time(t);
1321
1322 // Loop over values
1323 for (unsigned j = 0; j < n_value; j++)
1324 {
1325 data_pt->set_value(t, j, initial_value_fct[j](time_local));
1326 }
1327 }
1328 }
1329
1330 /// This function updates the Data's time history so that
1331 /// we can advance to the next timestep. For BDF schemes,
1332 /// we simply push the values backwards...
1333 void shift_time_values(Data* const& data_pt)
1334 {
1335 // Find number of values stored
1336 unsigned n_value = data_pt->nvalue();
1337 // Storage for velocity need to be here to be in scope
1338 Vector<double> velocity(n_value);
1339
1340 // Find the number of history values that are stored
1341 const unsigned nt_value = nprev_values();
1342
1343 // If adaptive, find the velocities
1344 if (adaptive_flag())
1345 {
1346 time_derivative(1, data_pt, velocity);
1347 }
1348
1349 // Loop over the values
1350 for (unsigned j = 0; j < n_value; j++)
1351 {
1352 // Set previous values to the previous value, if not a copy
1353 if (data_pt->is_a_copy(j) == false)
1354 {
1355 // Loop over times, in reverse order
1356 for (unsigned t = nt_value; t > 0; t--)
1357 {
1358 data_pt->set_value(t, j, data_pt->value(t - 1, j));
1359 }
1360
1361 // If we are using the adaptive scheme
1362 if (adaptive_flag())
1363 {
1364 // Set the velocity
1365 data_pt->set_value(nt_value + 1, j, velocity[j]);
1366 }
1367 }
1368 }
1369 }
1370
1371 /// This function advances the time history of the positions
1372 /// at a node.
1373 void shift_time_positions(Node* const& node_pt)
1374 {
1375 // Find the number of coordinates
1376 unsigned n_dim = node_pt->ndim();
1377 // Find the number of position types
1378 unsigned n_position_type = node_pt->nposition_type();
1379
1380 // Find number of stored timesteps
1381 unsigned n_tstorage = ntstorage();
1382
1383 // Storage for the velocity
1384 double velocity[n_position_type][n_dim];
1385
1386 // If adaptive, find the velocities
1387 if (adaptive_flag())
1388 {
1389 // Loop over the variables
1390 for (unsigned i = 0; i < n_dim; i++)
1391 {
1392 for (unsigned k = 0; k < n_position_type; k++)
1393 {
1394 // Initialise velocity to zero
1395 velocity[k][i] = 0.0;
1396 // Loop over all history data
1397 for (unsigned t = 0; t < n_tstorage; t++)
1398 {
1399 velocity[k][i] += Weight(1, t) * node_pt->x_gen(t, k, i);
1400 }
1401 }
1402 }
1403 }
1404
1405 // Loop over the positions
1406 for (unsigned i = 0; i < n_dim; i++)
1407 {
1408 // If the position is not a copy
1409 if (node_pt->position_is_a_copy(i) == false)
1410 {
1411 // Loop over the position types
1412 for (unsigned k = 0; k < n_position_type; k++)
1413 {
1414 // Loop over stored times, and set values to previous values
1415 for (unsigned t = NSTEPS; t > 0; t--)
1416 {
1417 node_pt->x_gen(t, k, i) = node_pt->x_gen(t - 1, k, i);
1418 }
1419
1420 // If we are using the adaptive scheme, set the velocity
1421 if (adaptive_flag())
1422 {
1423 node_pt->x_gen(NSTEPS + 1, k, i) = velocity[k][i];
1424 }
1425 }
1426 }
1427 }
1428 }
1429
1430 /// Set the weights
1432
1433 /// Number of previous values available.
1434 unsigned nprev_values() const
1435 {
1436 return NSTEPS;
1437 }
1438
1439 /// Number of timestep increments that need to be stored by the scheme
1440 unsigned ndt() const
1441 {
1442 return NSTEPS;
1443 }
1444
1445 /// Function to set the predictor weights
1447
1448 /// Function to calculate predicted positions at a node
1450
1451 /// Function to calculate predicted data values in a Data object
1452 void calculate_predicted_values(Data* const& data_pt);
1453
1454 /// Function to set the error weights
1456
1457 /// Compute the error in the position i at a node
1458 double temporal_error_in_position(Node* const& node_pt, const unsigned& i);
1459
1460 /// Compute the error in the value i in a Data structure
1461 double temporal_error_in_value(Data* const& data_pt, const unsigned& i);
1462 };
1463
1464} // namespace oomph
1465
1466#endif
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
////////////////////////////////////////////////////////////////////// //////////////////////////////...
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 calculate_predicted_positions(Node *const &node_pt)
Function to calculate predicted positions at a node.
void assign_initial_positions_impulsive(Node *const &node_pt)
Initialise the time-history for the nodal positions corresponding to an impulsive start.
void calculate_predicted_values(Data *const &data_pt)
Function to calculate predicted data values in a Data object.
BDF(const bool &adaptive=false)
Constructor for the case when we allow adaptive timestepping.
void assign_initial_values_impulsive(Data *const &data_pt)
Initialise the time-history for the Data values, corresponding to an impulsive start.
void set_weights()
Set the weights.
double Error_weight
Private data for the error weight.
unsigned ndt() const
Number of timestep increments that need to be stored by the scheme.
unsigned order() const
Return the actual order of the scheme.
void operator=(const BDF &)=delete
Broken assignment operator.
double(* InitialConditionFctPt)(const double &t)
Typedef for function that returns the (scalar) initial value at a given value of the continuous time ...
Vector< double > Predictor_weight
Private data for the predictor 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.
void shift_time_positions(Node *const &node_pt)
This function advances the time history of the positions at a node.
BDF(const BDF &)=delete
Broken copy constructor.
void assign_initial_data_values(Data *const &data_pt, Vector< InitialConditionFctPt > initial_value_fct)
Initialise the time-history for the Data values, corresponding to given time history,...
double temporal_error_in_position(Node *const &node_pt, const unsigned &i)
Compute the error in the position i at a node.
unsigned nprev_values() const
Number of previous values available.
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
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:485
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:491
void initialise(const T &val)
Initialize all values in the matrix to val.
Definition: matrices.h:514
void resize(const unsigned long &n)
Resize to a square nxn matrix; any values already present will be transfered.
Definition: matrices.h:498
A Base class for explicit timesteppers.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
void set_weights()
Set weights.
void disable_degrade_first_derivatives_to_bdf1()
Disable degradation to first order BDF.
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.
NewmarkBDF(const NewmarkBDF &)=delete
Broken copy constructor.
bool Degrade_to_bdf1_for_first_derivs
Boolean flag to indicate degradation of scheme to first order BDF (for first derivs/veloc); usually f...
void shift_time_positions(Node *const &node_pt)
This function updates a nodal time history so that we can advance to the next timestep.
NewmarkBDF()
Constructor: Pass pointer to global time. We set up a timestepping scheme with NSTEPS+2 doubles to re...
void operator=(const NewmarkBDF &)=delete
Broken assignment operator.
Vector< double > Newmark_veloc_weight
Original Newmark weights for velocities (needed when shifting history values – they're used when upda...
void set_newmark_veloc_weights(const double &dt)
Set original Newmark weights for velocities (needed when shifting history values – they're used when ...
void enable_degrade_first_derivatives_to_bdf1()
Degrade scheme to first order BDF (for first derivs/veloc); usually for start-up.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:912
double Beta2
Second Newmark parameter (usually 0.5)
Newmark()
Constructor: Pass pointer to global time. We set up a timestepping scheme with NSTEPS+2 doubles to re...
Definition: timesteppers.h:917
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(* InitialConditionFctPt)(const double &t)
Typedef for function that returns the (scalar) initial value at a given value of the continuous time ...
Definition: timesteppers.h:955
unsigned nprev_values() const
Number of previous values available.
unsigned ndt() const
Number of timestep increments that need to be stored by the scheme.
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.
double Beta1
First Newmark parameter (usually 0.5)
void assign_initial_values_impulsive(Data *const &data_pt)
Initialise the time-history for the values, corresponding to an impulsive start.
unsigned order() const
The actual order (accuracy of the scheme)
Definition: timesteppers.h:934
void operator=(const Newmark &)=delete
Broken assignment operator.
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.
Newmark(const Newmark &)=delete
Broken copy constructor.
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(* NodeInitialConditionFctPt)(const double &t, const Vector< double > &x)
Typedef for function that returns the (scalar) initial value at a given value of the continuous time ...
Definition: timesteppers.h:971
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
virtual bool position_is_a_copy() const
Return whether any position coordinate has been copied (always false)
Definition: nodes.h:1113
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:1016
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i). ‘Type’: k; Coordinate direction: i.
Definition: nodes.h:1126
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation....
Definition: nodes.cc:2408
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: problem.h:151
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:681
double(* InitialConditionFctPt)(const double &t)
Typedef for function that returns the (scalar) initial value at a given value of the continuous time ...
Definition: timesteppers.h:760
void assign_initial_values_impulsive(Data *const &data_pt)
Initialise the time-history for the Data values, corresponding to an impulsive start.
Definition: timesteppers.h:713
void set_weights()
Set weights.
Definition: timesteppers.h:839
double weight(const unsigned &i, const unsigned &j) const
Dummy: Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:868
void operator=(const Steady &)=delete
Broken assignment operator.
void shift_time_positions(Node *const &node_pt)
This function advances the time history of the positions at a node.
Definition: timesteppers.h:813
static double One
Static variable to hold the value 1.0.
Definition: timesteppers.h:882
static Time Dummy_time
Definition: timesteppers.h:888
Steady(const Steady &)=delete
Broken copy constructor.
Steady()
Constructor: Creates storage for NSTEPS previous timesteps and can evaluate up to 2nd derivatives (th...
Definition: timesteppers.h:687
void assign_initial_data_values(Data *const &data_pt, Vector< InitialConditionFctPt > initial_value_fct)
Initialise the time-history for the Data values, corresponding to given time history,...
Definition: timesteppers.h:765
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....
Definition: timesteppers.h:791
static double Zero
Static variable to hold the value 0.0.
Definition: timesteppers.h:885
unsigned order() const
Return the actual order of the scheme. Returning zero here – doesn't make much sense,...
Definition: timesteppers.h:706
unsigned nprev_values() const
Number of previous values available.
Definition: timesteppers.h:856
void assign_initial_positions_impulsive(Node *const &node_pt)
Initialise the time-history for the nodal positions corresponding to an impulsive start.
Definition: timesteppers.h:733
unsigned ndt() const
Number of timestep increments that need to be stored by the scheme.
Definition: timesteppers.h:862
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
virtual unsigned ndt() const =0
Number of timestep increments that are required by the scheme.
virtual void shift_time_values(Data *const &data_pt)=0
This function advances the Data's time history so that we can move on to the next timestep.
virtual void set_weights()=0
Function to set the weights for present timestep (don't need to pass present timestep or previous tim...
TimeStepper(const TimeStepper &)=delete
Broken copy constructor.
DenseMatrix< double > Weight
Storage for the weights associated with the timestepper.
Definition: timesteppers.h:237
virtual unsigned nprev_values() const =0
Number of previous values available: 0 for static, 1 for BDF<1>,...
ExplicitTimeStepper * explicit_predictor_pt()
Get the pointer to the explicit timestepper to use as a predictor in adaptivity if Predict_by_explici...
Definition: timesteppers.h:403
virtual void calculate_predicted_values(Data *const &data_pt)
Do the predictor step for data stored in a Data object (currently empty – overwrite for specific sche...
Definition: timesteppers.h:634
virtual ~TimeStepper()
virtual destructor
Definition: timesteppers.cc:37
bool Shut_up_in_assign_initial_data_values
Boolean to indicate if the timestepper will output warnings when setting possibly an incorrect number...
Definition: timesteppers.h:256
virtual unsigned order() const
Actual order (accuracy) of the scheme.
Definition: timesteppers.h:324
virtual double temporal_error_in_value(Data *const &data_pt, const unsigned &i)
Compute the error in the value i in a Data structure zero here – overwrite for specific scheme.
Definition: timesteppers.h:654
virtual void set_predictor_weights()
Set the weights for the predictor previous timestep (currently empty – overwrite for specific scheme)
Definition: timesteppers.h:630
virtual void calculate_predicted_positions(Node *const &node_pt)
Do the predictor step for the positions at a node (currently empty — overwrite for a specific scheme)
Definition: timesteppers.h:638
bool Predict_by_explicit_step
Flag: is adaptivity done by taking a separate step using an ExplicitTimeStepper object?
Definition: timesteppers.h:260
unsigned predictor_storage_index() const
Return the time-index in each Data where predicted values are stored if the timestepper is adaptive.
Definition: timesteppers.h:438
void enable_warning_in_assign_initial_data_values()
Enable the output of warnings due to possible fct pointer vector size mismatch in assign_initial_data...
Definition: timesteppers.h:462
virtual void actions_before_timestep(Problem *problem_pt)
Interface for any actions that need to be performed before a time step.
Definition: timesteppers.h:662
void check_predicted_values_up_to_date() const
Check that the predicted values are the ones we want.
Definition: timesteppers.h:423
Time * Time_pt
Pointer to discrete time storage scheme.
Definition: timesteppers.h:234
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
std::string Type
String that indicates the type of the timestepper (e.g. "BDF", "Newmark", etc.)
Definition: timesteppers.h:241
virtual void actions_after_timestep(Problem *problem_pt)
Interface for any actions that need to be performed after a time step.
Definition: timesteppers.h:666
void time_derivative(const unsigned &i, Node *const &node_pt, Vector< double > &deriv)
Evaluate i-th derivative of all values in Node and return in Vector deriv[] (this can't be simply com...
Definition: timesteppers.h:536
ExplicitTimeStepper * Explicit_predictor_pt
Pointer to explicit time stepper to use as predictor if Predict_by_explicit_step is set....
Definition: timesteppers.h:265
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
virtual void shift_time_positions(Node *const &node_pt)=0
This function advances the time history of the positions at a node. The default should be OK,...
bool Is_steady
Bool to indicate if the timestepper is steady, i.e. its time-derivatives evaluate to zero....
Definition: timesteppers.h:251
double & time()
Return current value of continous time.
Definition: timesteppers.h:332
Time *& time_pt()
Access function for the pointer to time (can't have a paranoid test for null pointers because this is...
Definition: timesteppers.h:588
virtual unsigned nprev_values_for_value_at_evaluation_time() const
Number of previous values needed to calculate the value at the current time. i.e. how many previous v...
Definition: timesteppers.h:359
bool adaptive_flag() const
Function to indicate whether the scheme is adaptive (false by default)
Definition: timesteppers.h:623
TimeStepper()
Broken empty constructor.
Definition: timesteppers.h:301
void operator=(const TimeStepper &)=delete
Broken assignment operator.
virtual void assign_initial_values_impulsive(Data *const &data_pt)=0
Initialise the time-history for the Data values corresponding to an impulsive start.
void disable_warning_in_assign_initial_data_values()
Disable the output of warnings due to possible fct pointer vector size mismatch in assign_initial_dat...
Definition: timesteppers.h:469
void time_derivative(const unsigned &i, Data *const &data_pt, Vector< double > &deriv)
Evaluate i-th derivative of all values in Data and return in Vector deriv[].
Definition: timesteppers.h:502
void make_steady()
Function to make the time stepper temporarily steady. This is trivially achieved by setting all the w...
Definition: timesteppers.h:374
TimeStepper(const unsigned &tstorage, const unsigned &max_deriv)
Constructor. Pass the amount of storage required by timestepper (present value + history values) and ...
Definition: timesteppers.h:279
unsigned highest_derivative() const
Highest order derivative that the scheme can compute.
Definition: timesteppers.h:318
virtual void undo_make_steady()
Reset the is_steady status of a specific TimeStepper to its default and re-assign the weights.
Definition: timesteppers.h:482
double time() const
Return current value of continous time.
Definition: timesteppers.h:338
virtual void assign_initial_positions_impulsive(Node *const &node_pt)=0
Initialiset the positions for the node corresponding to an impulsive start.
double time_derivative(const unsigned &i, Data *const &data_pt, const unsigned &j)
Evaluate i-th derivative of j-th value in Data.
Definition: timesteppers.h:518
void update_predicted_time(const double &new_time)
Set the time that the current predictions apply for, only needed for paranoid checks when doing Predi...
Definition: timesteppers.h:417
bool predict_by_explicit_step() const
Flag: is adaptivity done by taking a separate step using an ExplicitTimeStepper object?
Definition: timesteppers.h:396
virtual double temporal_error_in_position(Node *const &node_pt, const unsigned &i)
Compute the error in the position i at a node zero here – overwrite for specific scheme.
Definition: timesteppers.h:646
const DenseMatrix< double > * weights_pt() const
Get a (const) pointer to the weights.
Definition: timesteppers.h:475
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
std::string type() const
Return string that indicates the type of the timestepper (e.g. "BDF", "Newmark", etc....
Definition: timesteppers.h:490
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:389
virtual void set_error_weights()
Set the weights for the error computation, (currently empty – overwrite for specific scheme)
Definition: timesteppers.h:642
double time_derivative(const unsigned &i, Node *const &node_pt, const unsigned &j)
Evaluate i-th derivative of j-th value in Node. Note the use of the node's value() function so that h...
Definition: timesteppers.h:555
double Predicted_time
Store the time that the predicted values currently stored are at, to compare for paranoid checks.
Definition: timesteppers.h:269
int Predictor_storage_index
The time-index in each Data object where predicted values are stored. -1 if not set.
Definition: timesteppers.h:273
bool Adaptive_Flag
Boolean variable to indicate whether the timestepping scheme can be adaptive.
Definition: timesteppers.h:245
void set_predictor_pt(ExplicitTimeStepper *_pred_pt)
Set the pointer to the explicit timestepper to use as a predictor in adaptivity if Predict_by_explici...
Definition: timesteppers.h:410
Class to keep track of discrete/continous time. It is essential to have a single Time object when usi...
Definition: timesteppers.h:63
double Continuous_time
Pointer to the value of the continuous time.
Definition: timesteppers.h:66
Time(const Time &)=delete
Broken copy constructor.
void operator=(const Time &)=delete
Broken assignment operator.
Vector< double > Dt
Vector that stores the values of the current and previous timesteps.
Definition: timesteppers.h:69
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
double time(const unsigned &t=0) const
Return the value of the continuous time at the t-th previous time level (t=0: current; t>0 previous).
Definition: timesteppers.h:150
double dt(const unsigned &t=0) const
Return the value of the t-th stored timestep (t=0: present; t>0: previous), const version.
Definition: timesteppers.h:143
void initialise_dt(const Vector< double > &dt_)
Set the value of the timesteps to be equal to the values passed in a vector.
Definition: timesteppers.h:107
void shift_dt()
Update all stored values of dt by shifting each value along the array. This function must be called b...
Definition: timesteppers.h:174
void initialise_dt(const double &dt_)
Set all timesteps to the same value, dt.
Definition: timesteppers.h:99
double & dt(const unsigned &t=0)
Return the value of the t-th stored timestep (t=0: present; t>0: previous).
Definition: timesteppers.h:136
unsigned ndt() const
Return the number of timesteps stored.
Definition: timesteppers.h:129
void resize(const unsigned &n_dt)
Resize the vector holding the number of previous timesteps and initialise the new values to zero.
Definition: timesteppers.h:93
Time()
Constructor: Do not allocate any storage for previous timesteps, but set the initial value of the tim...
Definition: timesteppers.h:74
~Time()
Destructor: empty.
Definition: timesteppers.h:120
Time(const unsigned &ndt)
Constructor: Pass the number of timesteps to be stored and set the initial value of time to zero.
Definition: timesteppers.h:78
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types).
//////////////////////////////////////////////////////////////////// ////////////////////////////////...