Forked from
agnumpde / dune-tectonic
657 commits behind the upstream repository.
-
Elias Pipping authoredElias Pipping authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
timestepping.org 12.55 KiB
#include <dune/fufem/arithmetic.hh>
\documentclass{scrartcl}
\usepackage{amsmath}
\begin{document}
Abstract TimeSteppingScheme
template <class VectorType, class MatrixType, class FunctionType, int dim>
class TimeSteppingScheme
{
public:
void virtual
setup(VectorType const &ell,
double _tau,
double time,
VectorType &problem_rhs,
VectorType &problem_iterate,
MatrixType &problem_A) = 0;
void virtual
postProcess(VectorType const &problem_iterate) = 0;
void virtual
extractDisplacement(VectorType &displacement) const = 0;
void virtual
extractVelocity(VectorType &velocity) const = 0;
};
TimeSteppingScheme: Implicit Euler
template <class VectorType, class MatrixType, class FunctionType, int dim>
class ImplicitEuler : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
{
public:
ImplicitEuler(MatrixType const &_A,
VectorType const &_u_initial,
VectorType const &_ud_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction)
: A(_A),
u(_u_initial),
ud(_ud_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction)
{}
void virtual
setup(VectorType const &ell,
double _tau,
double time,
VectorType &problem_rhs,
VectorType &problem_iterate,
MatrixType &problem_A)
{
postProcessCalled = false;
tau = _tau;
// This function is only called once for every timestep
ud_old = ud;
u_old = u;
problem_rhs = ell;
A.mmv(u_old, problem_rhs);
// For fixed tau, we'd only really have to do this once
problem_A = A;
problem_A *= tau;
// ud_old makes a good initial iterate; we could use anything, though
problem_iterate = ud_old;
for (size_t i=0; i < dirichletNodes.size(); ++i)
if (dirichletNodes[i].count() == dim) {
problem_iterate[i] = 0;
dirichletFunction.evaluate(time, problem_iterate[i][0]);
} else if (dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction prescribed
}
void virtual postProcess(VectorType const &problem_iterate)
{
postProcessCalled = true;
ud = problem_iterate;
u = u_old;
Arithmetic::addProduct(u, tau, ud);
}
void virtual extractDisplacement(VectorType &displacement) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
displacement = u;
}
void virtual extractVelocity(VectorType &velocity) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = ud;
}
private:
MatrixType const &A;
VectorType u;
VectorType ud;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
VectorType u_old;
VectorType ud_old;
double tau;
bool postProcessCalled = false;
};
TimeSteppingScheme: Implicit Two-Step
\documentclass{scrartcl}
\usepackage{amsmath}
\begin{document}
We start out with
\begin{align}
a(u_1, v - \dot u_1) + j(v) - j(\dot u_1) \ge \ell(v-\dot u_1)
\end{align}
With the two-step implicit scheme
\begin{equation*}
\tau \dot u_1 = \frac 32 u_1 - 2 u_0 + \frac 12 u_{-1}
\end{equation*}
or equivalently
\begin{equation*}
\frac 23 \left( \tau \dot u_1 + 2 u_0 - \frac 12 u_{-1} \right) = u_1
\end{equation*}
we obtain
\begin{equation*}
\frac 23 \tau a(\dot u_1, v - \dot u_1) + j(v) - j(\dot u_1) \ge \ell(v-\dot u_1) - a\left(\frac 43 u_0 - \frac 13 u_{-1}, v - \dot u_1\right)
\end{equation*}
\end{document}
template <class VectorType, class MatrixType, class FunctionType, int dim>
class ImplicitTwoStep : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
{
public:
ImplicitTwoStep(MatrixType const &_A,
VectorType const &_u_initial,
VectorType const &_ud_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction)
: A(_A),
u(_u_initial),
ud(_ud_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction)
{}
// FIXME: handle case run == 1
void virtual
setup(VectorType const &ell,
double _tau,
double time,
VectorType &problem_rhs,
VectorType &problem_iterate,
MatrixType &problem_A)
{
postProcessCalled = false;
tau = _tau;
// This function is only called once for every timestep
u_old_old = u_old;
ud_old = ud;
u_old = u;
switch (state) {
// Perform an implicit Euler step since we lack information
case NO_SETUP:
state = FIRST_SETUP;
problem_rhs = ell;
A.mmv(u_old, problem_rhs);
problem_A = A;
problem_A *= tau;
break;
case FIRST_SETUP:
state = SECOND_SETUP;
// FALLTHROUGH
case SECOND_SETUP:
problem_rhs = ell;
A.usmv(-4.0/3.0, u_old, problem_rhs);
A.usmv(+1.0/3.0, u_old_old, problem_rhs);
// For fixed tau, we'd only really have to do this once
problem_A = A;
problem_A *= 2.0/3.0 * tau;
break;
default:
assert(false);
}
// ud_old makes a good initial iterate; we could use anything, though
problem_iterate = ud_old;
for (size_t i=0; i < dirichletNodes.size(); ++i)
if (dirichletNodes[i].count() == dim) {
problem_iterate[i] = 0;
dirichletFunction.evaluate(time, problem_iterate[i][0]);
} else if (dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction prescribed
}
void virtual postProcess(VectorType const &problem_iterate)
{
postProcessCalled = true;
ud = problem_iterate;
switch (state) {
case FIRST_SETUP:
u = u_old;
Arithmetic::addProduct(u, tau, ud);
break;
case SECOND_SETUP:
u = 0.0;
Arithmetic::addProduct(u, tau, ud);
Arithmetic::addProduct(u, 2.0, u_old);
Arithmetic::addProduct(u, -.5, u_old_old);
u *= 2.0/3.0;
break;
default:
assert(false);
}
}
void virtual extractDisplacement(VectorType &displacement) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
displacement = u;
}
void virtual extractVelocity(VectorType &velocity) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = ud;
}
private:
MatrixType const &A;
VectorType u;
VectorType ud;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
VectorType u_old;
VectorType u_old_old;
VectorType ud_old;
double tau;
bool postProcessCalled = false;
// Handle a lack of information
enum state_type { NO_SETUP, FIRST_SETUP, SECOND_SETUP };
state_type state = NO_SETUP;
};
TimeSteppingScheme: Newmark
<<preamble>>
\noindent The Newmark scheme in its classical form with $\gamma = 1/2$
and $\beta = 1/4$ reads
\begin{align}
\label{eq:1} \dot u_1
&= \dot u_0 + \frac \tau 2 (\ddot u_0 + \ddot u_1 )\\
\label{eq:2} u_1
&= u_0 + \tau \dot u_0 + \frac {\tau^2}4 ( \ddot u_0 + \ddot u_1 )\text.
\intertext{We can also write \eqref{eq:1} as}
\label{eq:3} \ddot u_1
&= \frac 2\tau (\dot u_1 - \dot u_0) - \ddot u_0
\intertext{so that it yields}
\label{eq:4} u_1
&= u_0 + \tau \dot u_0 + \frac {\tau^2}4 \ddot u_0 + \frac {\tau^2}4 \left( \frac 2\tau (\dot u_1 - \dot u_0) - \ddot u_0\right)\\
&= u_0 + \tau \dot u_0 + \frac \tau 2 (\dot u_1 - \dot u_0)\nonumber\\
&= u_0 + \frac \tau 2 (\dot u_1 + \dot u_0)\nonumber
\end{align}
in conjunction with \eqref{eq:2}. If we now consider the EVI
\begin{align*}
b(\ddot u_1, v - \dot u_1) + a(u_1, v - \dot u_1) + j(v) - j(\dot u_1)
&\ge \ell (v - \dot u_1)
\intertext{with unknowns $u_1$, $\dot u_1$, and $\ddot u_1$, we first derive}
\frac 2\tau b(\dot u_1, v - \dot u_1) + a(u_1, v - \dot u_1) + j(v) - j(\dot u_1)
&\ge \ell (v - \dot u_1) + b\left(\frac 2\tau \dot u_0 + \ddot u_0, v - \dot u_1\right)
\intertext{from \eqref{eq:3} and then}
\frac 2\tau b(\dot u_1, v - \dot u_1) + \frac \tau 2 a(\dot u_1, v - \dot u_1) + j(v) - j(\dot u_1)
&\ge \ell (v - \dot u_1) + b\left(\frac 2\tau \dot u_0 + \ddot u_0, v - \dot u_1\right)\\
&\quad - a\left(u_0 + \frac \tau 2 \dot u_0, v - \dot u_1\right)
\end{align*}
from \eqref{eq:4}. The only unknown is now $\dot u_1$.
\end{document}
template <class VectorType, class MatrixType, class FunctionType, int dim>
class Newmark : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
{
public:
Newmark(MatrixType const &_A,
MatrixType const &_B,
VectorType const &_u_initial,
VectorType const &_ud_initial,
VectorType const &_udd_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction)
: A(_A),
B(_B),
u(_u_initial),
ud(_ud_initial),
udd(_udd_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction)
{}
void virtual
setup(VectorType const &ell,
double _tau,
double time,
VectorType &problem_rhs,
VectorType &problem_iterate,
MatrixType &problem_A)
{
postProcessCalled = false;
tau = _tau;
// This function is only called once for every timestep
udd_old = udd;
ud_old = ud;
u_old = u;
problem_rhs = ell;
B.usmv( 2.0/tau, ud_old, problem_rhs);
B.usmv( 1.0, udd_old, problem_rhs);
A.usmv( -1.0, u_old, problem_rhs);
A.usmv(-tau/2.0, ud_old, problem_rhs);
// For fixed tau, we'd only really have to do this once
problem_A = A;
problem_A *= tau/2.0;
Arithmetic::addProduct(problem_A, 2.0/tau, B);
// ud_old makes a good initial iterate; we could use anything, though
problem_iterate = ud_old;
for (size_t i=0; i < dirichletNodes.size(); ++i)
if (dirichletNodes[i].count() == dim) {
problem_iterate[i] = 0;
dirichletFunction.evaluate(time, problem_iterate[i][0]);
} else if (dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction prescribed
}
void virtual postProcess(VectorType const &problem_iterate)
{
postProcessCalled = true;
ud = problem_iterate;
u = u_old;
Arithmetic::addProduct(u, tau/2.0, ud);
Arithmetic::addProduct(u, tau/2.0, ud_old);
udd = 0;
Arithmetic::addProduct(udd, 2.0/tau, ud);
Arithmetic::addProduct(udd, -2.0/tau, ud_old);
Arithmetic::addProduct(udd, -1.0, udd_old);
}
void virtual extractDisplacement(VectorType &displacement) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
displacement = u;
}
void virtual extractVelocity(VectorType &velocity) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = ud;
}
private:
MatrixType const &A;
MatrixType const &B;
VectorType u;
VectorType ud;
VectorType udd;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
VectorType u_old;
VectorType ud_old;
VectorType udd_old;
double tau;
bool postProcessCalled = false;
};
C++ code generation
<<includes>>
<<TimeSteppingScheme>>
<<ImplicitEuler>>
<<ImplicitTwoStep>>
<<Newmark>>