ode_elements.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#ifndef OOMPH_ODE_ELEMENTS_H
27#define OOMPH_ODE_ELEMENTS_H
28
29#include "../generic/oomph_definitions.h"
30#include "../generic/oomph_utilities.h"
31
32#include "../generic/matrices.h"
33#include "../generic/Vector.h"
34#include "../generic/elements.h"
35#include "../generic/timesteppers.h"
36
37namespace oomph
38{
39 /// Element for integrating an initial value ODE
41 {
42 public:
43 /// Default constructor: null any pointers
45 {
47
48 Use_fd_jacobian = false;
49 }
50
51 /// Constructor: Pass time stepper and a solution function pointer, then
52 /// build the element.
53 ODEElement(TimeStepper* time_stepper_pt,
54 SolutionFunctorBase* exact_solution_pt)
55 {
56 build(time_stepper_pt, exact_solution_pt);
57 }
58
59 /// Store pointers, create internal data.
60 void build(TimeStepper* time_stepper_pt,
61 SolutionFunctorBase* exact_solution_pt)
62 {
63 Exact_solution_pt = exact_solution_pt;
64
65 Vector<double> exact = this->exact_solution(0);
66 unsigned nvalue = exact.size();
67
68 add_internal_data(new Data(time_stepper_pt, nvalue));
69
70 Use_fd_jacobian = false;
71 }
72
73
74 virtual ~ODEElement() {}
75
76 unsigned nvalue() const
77 {
78 return internal_data_pt(0)->nvalue();
79 }
80
81 /// Get residuals
83 {
84 // Get pointer to one-and-only internal data object
85 Data* dat_pt = internal_data_pt(0);
86
87 // Get it's values
88 Vector<double> u(nvalue(), 0.0);
89 dat_pt->value(u);
90
91 // Get time stepper
92 TimeStepper* time_stepper_pt = dat_pt->time_stepper_pt();
93
94 // Get continuous time
95 double t = time_stepper_pt->time();
96
98 for (unsigned j = 0, nj = deriv.size(); j < nj; j++)
99 {
100 // Get dudt approximation from time stepper
101 double dudt = time_stepper_pt->time_derivative(1, dat_pt, j);
102
103 // Residual is difference between the exact derivative and our
104 // time stepper's derivative estimate.
105 residuals[j] = deriv[j] - dudt;
106 }
107 }
108
110 DenseMatrix<double>& jacobian)
111 {
112 // Get residuals
114
116 {
117 // get df/du jacobian
118 double t = internal_data_pt(0)->time_stepper_pt()->time();
119 Vector<double> dummy, u(nvalue(), 0.0);
120 internal_data_pt(0)->value(u);
121 Exact_solution_pt->jacobian(t, dummy, u, jacobian);
122
123 // We need jacobian of residual = f - dudt so subtract diagonal
124 // (dudt)/du term.
125 const double a = internal_data_pt(0)->time_stepper_pt()->weight(1, 0);
126 const unsigned n = nvalue();
127 for (unsigned i = 0; i < n; i++)
128 {
129 jacobian(i, i) -= a;
130 }
131 }
132 else
133 {
134 // Use FD for jacobian
136 residuals, jacobian, true);
137 }
138 }
139
142 {
144 for (unsigned j = 0, nj = nvalue(); j < nj; j++)
145 {
146 mm(j, j) = 1;
147 }
148 }
149
150 /// Exact solution
151 Vector<double> exact_solution(const double& t) const
152 {
153#ifdef PARANOID
154 if (Exact_solution_pt == 0)
155 {
156 throw OomphLibError("No exact solution function",
157 OOMPH_EXCEPTION_LOCATION,
158 OOMPH_CURRENT_FUNCTION);
159 }
160#endif
161 Vector<double> dummy_x;
162 return (*Exact_solution_pt)(t, dummy_x);
163 }
164
165 /// Exact solution
167 {
168#ifdef PARANOID
169 if (Exact_solution_pt == 0)
170 {
171 throw OomphLibError("No derivative function",
172 OOMPH_EXCEPTION_LOCATION,
173 OOMPH_CURRENT_FUNCTION);
174 }
175#endif
176 Vector<double> dummy_x;
177 return Exact_solution_pt->derivative(t, dummy_x, u);
178 }
179
181
183 };
184
185} // namespace oomph
186
187#endif
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
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
A Generalised Element class.
Definition: elements.h:73
void fill_in_jacobian_from_internal_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
Calculate the contributions to the jacobian from the internal degrees of freedom using finite differe...
Definition: elements.cc:1102
unsigned add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
Definition: elements.cc:62
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
Element for integrating an initial value ODE.
Definition: ode_elements.h:41
ODEElement(TimeStepper *time_stepper_pt, SolutionFunctorBase *exact_solution_pt)
Constructor: Pass time stepper and a solution function pointer, then build the element.
Definition: ode_elements.h:53
Vector< double > derivative_function(const double &t, const Vector< double > &u)
Exact solution.
Definition: ode_elements.h:166
Vector< double > exact_solution(const double &t) const
Exact solution.
Definition: ode_elements.h:151
virtual void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the elemental contribution to the jacobian matrix. and the residuals vector. Note that this funct...
Definition: ode_elements.h:109
SolutionFunctorBase * Exact_solution_pt
Definition: ode_elements.h:180
ODEElement()
Default constructor: null any pointers.
Definition: ode_elements.h:44
virtual void fill_in_contribution_to_residuals(Vector< double > &residuals)
Get residuals.
Definition: ode_elements.h:82
void build(TimeStepper *time_stepper_pt, SolutionFunctorBase *exact_solution_pt)
Store pointers, create internal data.
Definition: ode_elements.h:60
virtual ~ODEElement()
Definition: ode_elements.h:74
unsigned nvalue() const
Definition: ode_elements.h:76
virtual void fill_in_contribution_to_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &mm)
Add the elemental contribution to the mass matrix matrix. and the residuals vector....
Definition: ode_elements.h:140
An OomphLibError object which should be thrown when an run-time error is encountered....
Function base class for exact solutions/initial conditions/boundary conditions. This is needed so tha...
virtual void jacobian(const double &t, const Vector< double > &x, const Vector< double > &u, DenseMatrix< double > &jacobian) const
The derivatives of the derivative function with respect to u (note that this is not quite the jacobia...
virtual Vector< double > derivative(const double &t, const Vector< double > &x, const Vector< double > &u) const =0
Call the derivative function.
virtual bool have_jacobian() const
Is a jacobian function implemented?
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
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
double & time()
Return current value of continous time.
Definition: timesteppers.h:332
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...