Skip to content
Snippets Groups Projects
Commit f23df2c7 authored by Elias Pipping's avatar Elias Pipping Committed by Elias Pipping
Browse files

Convert timestepping.cc to org-babel

parent e42838d2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/.gitignore 0 → 100644
+ 1
− 0
View file @ f23df2c7
timestepping.cc
src/Makefile.am
+ 5
− 0
View file @ f23df2c7
......@@ -141,3 +141,8 @@ DISTCHECK_CONFIGURE_FLAGS = \
CXX="$(CXX)" CC="$(CC)"
include $(top_srcdir)/am/global-rules
$(srcdir)/timestepping.cc: $(srcdir)/timestepping.org
emacs -Q --batch --eval \
"(let (vc-handled-backends) \
(org-babel-tangle-file \"$<\" nil 'c++))"
src/timestepping.cc deleted 100644 → 0
+ 0
− 211
View file @ e42838d2
#include <dune/fufem/arithmetic.hh>
template <class VectorType, class MatrixType, class FunctionType, int dim>
class TimeSteppingScheme {
public:
TimeSteppingScheme(VectorType const &_ell, MatrixType const &_A,
VectorType const &_u_old,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction, double _time,
double _tau)
: ell(_ell),
A(_A),
u_old(_u_old),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction),
time(_time),
tau(_tau) {}
virtual ~TimeSteppingScheme() {}
void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
MatrixType &problem_A) const = 0;
void virtual extractSolution(VectorType const &problem_iterate,
VectorType &solution) const = 0;
void virtual extractVelocity(VectorType const &problem_iterate,
VectorType &velocity) const = 0;
protected:
VectorType const &ell;
MatrixType const &A;
VectorType const &u_old;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
double const time;
double const tau;
};
template <class VectorType, class MatrixType, class FunctionType, int dim>
class ImplicitEuler
: public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
public:
ImplicitEuler(
VectorType const &_ell, MatrixType const &_A, VectorType const &_u_old,
VectorType const *_u_old_old_ptr, // FIXME: this should not be necessary
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction, double _time, double _tau)
: TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(
_ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
u_old_old_ptr(_u_old_old_ptr) {}
void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
MatrixType &problem_A) const {
problem_rhs = this->ell;
this->A.mmv(this->u_old, problem_rhs);
problem_A = this->A;
problem_A *= this->tau;
// Use the old velocity as an initial iterate if possible
if (u_old_old_ptr) {
problem_iterate = this->u_old;
problem_iterate -= *u_old_old_ptr;
problem_iterate /= this->tau;
} else
problem_iterate = 0.0;
for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
if (this->dirichletNodes[i].count() == dim) {
problem_iterate[i] = 0;
this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
} else if (this->dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction prescribed
}
void virtual extractSolution(VectorType const &problem_iterate,
VectorType &solution) const {
solution = this->u_old;
solution.axpy(this->tau, problem_iterate);
}
void virtual extractVelocity(VectorType const &problem_iterate,
VectorType &velocity) const {
velocity = problem_iterate;
}
private:
VectorType const *u_old_old_ptr;
};
template <class VectorType, class MatrixType, class FunctionType, int dim>
class ImplicitTwoStep
: public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
public:
ImplicitTwoStep(VectorType const &_ell, MatrixType const &_A,
VectorType const &_u_old, VectorType const *_u_old_old_ptr,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction, double _time,
double _tau)
: TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(
_ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
u_old_old_ptr(_u_old_old_ptr) {}
void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
MatrixType &problem_A) const {
problem_rhs = this->ell;
this->A.usmv(-4.0 / 3.0, this->u_old, problem_rhs);
this->A.usmv(+1.0 / 3.0, *u_old_old_ptr, problem_rhs);
problem_A = this->A;
problem_A *= 2.0 / 3.0 * this->tau;
// Use the old velocity as an initial iterate
problem_iterate = this->u_old;
problem_iterate -= *u_old_old_ptr;
problem_iterate /= this->tau;
for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
if (this->dirichletNodes[i].count() == dim) {
problem_iterate[i] = 0;
this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
} else if (this->dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction prescribed
}
void virtual extractSolution(VectorType const &problem_iterate,
VectorType &solution) const {
solution = problem_iterate;
solution *= this->tau;
solution.axpy(2, this->u_old);
solution.axpy(-.5, *u_old_old_ptr);
solution *= 2.0 / 3.0;
// Check if we split correctly
{
VectorType test = problem_iterate;
test *= this->tau;
test.axpy(-1.5, solution);
test.axpy(+2, this->u_old);
test.axpy(-.5, *u_old_old_ptr);
assert(test.two_norm() < 1e-10);
}
}
void virtual extractVelocity(VectorType const &problem_iterate,
VectorType &velocity) const {
velocity = problem_iterate;
}
private:
VectorType const *u_old_old_ptr;
};
template <class VectorType, class MatrixType, class FunctionType, int dim>
class Newmark
: public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
public:
Newmark(VectorType const &_ell, MatrixType const &_A, MatrixType const &_B,
VectorType const &_u_old, VectorType const &_ud_old,
VectorType const &_udd_old,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction, double _time, double _tau)
: TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(
_ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
B(_B),
ud_old(_ud_old),
udd_old(_udd_old) {}
void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
MatrixType &problem_A) const {
problem_rhs = this->ell;
/* */ B.usmv(2.0 / this->tau, ud_old, problem_rhs);
/* */ B.usmv(1.0, udd_old, problem_rhs);
this->A.usmv(-1.0, this->u_old, problem_rhs);
this->A.usmv(-this->tau / 2.0, ud_old, problem_rhs);
// For fixed tau, we'd only really have to do this once
problem_A = this->A;
problem_A *= this->tau / 2.0;
Arithmetic::addProduct(problem_A, 2.0 / this->tau, B);
// ud_old makes a good initial iterate; we could use anything, though
problem_iterate = ud_old;
for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
if (this->dirichletNodes[i].count() == dim) {
problem_iterate[i] = 0;
this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
} else if (this->dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction prescribed
}
// u1 = u0 + tau/2 (du1 + du0)
void virtual extractSolution(VectorType const &problem_iterate,
VectorType &solution) const {
solution = this->u_old;
Arithmetic::addProduct(solution, this->tau / 2.0, problem_iterate); // ud
Arithmetic::addProduct(solution, this->tau / 2.0, ud_old);
}
void virtual extractVelocity(VectorType const &problem_iterate,
VectorType &velocity) const {
velocity = problem_iterate;
}
private:
MatrixType const &B;
VectorType const &ud_old;
VectorType const &udd_old;
};
src/timestepping.org 0 → 100644
+ 339
− 0
View file @ f23df2c7
#+name: includes
#+begin_src c++
#include <dune/fufem/arithmetic.hh>
#+end_src
#+name: preamble
#+begin_src latex
\documentclass{scrartcl}
\usepackage{amsmath}
\begin{document}
#+end_src
* Abstract TimeSteppingScheme
#+name: TimeSteppingScheme
#+begin_src c++
template <class VectorType, class MatrixType, class FunctionType, int dim>
class TimeSteppingScheme
{
public:
TimeSteppingScheme(VectorType const &_ell,
MatrixType const &_A,
VectorType const &_u_old,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction,
double _time,
double _tau)
: ell(_ell),
A(_A),
u_old(_u_old),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction),
time(_time),
tau(_tau)
{}
virtual ~TimeSteppingScheme()
{}
void virtual setup(VectorType &problem_rhs,
VectorType &problem_iterate,
MatrixType &problem_A) const = 0;
void virtual extractSolution(VectorType const &problem_iterate,
VectorType &solution) const = 0;
void virtual extractVelocity(VectorType const &problem_iterate,
VectorType &velocity) const = 0;
protected:
VectorType const &ell;
MatrixType const &A;
VectorType const &u_old;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
double const time;
double const tau;
};
#+end_src
* TimeSteppingScheme: Implicit Euler
#+name: ImplicitEuler
#+begin_src c++
template <class VectorType, class MatrixType, class FunctionType, int dim>
class ImplicitEuler : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
{
public:
ImplicitEuler(VectorType const &_ell,
MatrixType const &_A,
VectorType const &_u_old,
VectorType const *_u_old_old_ptr, // FIXME: this should not be necessary
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction,
double _time,
double _tau)
: TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
(_ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
u_old_old_ptr(_u_old_old_ptr)
{}
void virtual
setup(VectorType &problem_rhs,
VectorType &problem_iterate,
MatrixType &problem_A) const
{
problem_rhs = this->ell;
this->A.mmv(this->u_old, problem_rhs);
problem_A = this->A;
problem_A *= this->tau;
// Use the old velocity as an initial iterate if possible
if (u_old_old_ptr) {
problem_iterate = this->u_old;
problem_iterate -= *u_old_old_ptr;
problem_iterate /= this->tau;
} else
problem_iterate = 0.0;
for (size_t i=0; i < this->dirichletNodes.size(); ++i)
if (this->dirichletNodes[i].count() == dim) {
problem_iterate[i] = 0;
this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
} else if (this->dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction prescribed
}
void virtual extractSolution(VectorType const &problem_iterate,
VectorType &solution) const
{
solution = this->u_old;
solution.axpy(this->tau, problem_iterate);
}
void virtual extractVelocity(VectorType const &problem_iterate,
VectorType &velocity) const
{
velocity = problem_iterate;
}
private:
VectorType const *u_old_old_ptr;
};
#+end_src
* TimeSteppingScheme: Implicit Two-Step
#+begin_src latex :tangle twostep.tex :noweb yes
\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}
#+end_src
#+name: ImplicitTwoStep
#+begin_src c++
template <class VectorType, class MatrixType, class FunctionType, int dim>
class ImplicitTwoStep : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
{
public:
ImplicitTwoStep(VectorType const &_ell,
MatrixType const &_A,
VectorType const &_u_old,
VectorType const *_u_old_old_ptr,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction,
double _time,
double _tau)
: TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
(_ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
u_old_old_ptr(_u_old_old_ptr)
{}
void virtual setup(VectorType &problem_rhs,
VectorType &problem_iterate,
MatrixType &problem_A) const
{
problem_rhs = this->ell;
this->A.usmv(-4.0/3.0, this->u_old, problem_rhs);
this->A.usmv(+1.0/3.0, *u_old_old_ptr, problem_rhs);
problem_A = this->A;
problem_A *= 2.0/3.0 * this->tau;
// Use the old velocity as an initial iterate
problem_iterate = this->u_old;
problem_iterate -= *u_old_old_ptr;
problem_iterate /= this->tau;
for (size_t i=0; i < this->dirichletNodes.size(); ++i)
if (this->dirichletNodes[i].count() == dim) {
problem_iterate[i] = 0;
this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
} else if (this->dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction prescribed
}
void virtual extractSolution(VectorType const &problem_iterate,
VectorType &solution) const
{
solution = problem_iterate;
solution *= this->tau;
solution.axpy(2, this->u_old);
solution.axpy(-.5, *u_old_old_ptr);
solution *= 2.0/3.0;
// Check if we split correctly
{
VectorType test = problem_iterate;
test *= this->tau;
test.axpy(-1.5, solution);
test.axpy(+2, this->u_old);
test.axpy(-.5, *u_old_old_ptr);
assert(test.two_norm() < 1e-10);
}
}
void virtual extractVelocity(VectorType const &problem_iterate,
VectorType &velocity) const
{
velocity = problem_iterate;
}
private:
VectorType const *u_old_old_ptr;
};
#+end_src
* TimeSteppingScheme: Newmark
#+begin_src latex :tangle newmark.tex :noweb yes
<<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}
#+end_src
#+name: Newmark
#+begin_src c++
template <class VectorType, class MatrixType, class FunctionType, int dim>
class Newmark : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
{
public:
Newmark(VectorType const &_ell,
MatrixType const &_A,
MatrixType const &_B,
VectorType const &_u_old,
VectorType const &_ud_old,
VectorType const &_udd_old,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction,
double _time,
double _tau)
: TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
(_ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
B(_B),
ud_old(_ud_old),
udd_old(_udd_old)
{}
void virtual
setup(VectorType &problem_rhs,
VectorType &problem_iterate,
MatrixType &problem_A) const
{
problem_rhs = this->ell;
/* */ B.usmv( 2.0/this->tau, ud_old, problem_rhs);
/* */ B.usmv( 1.0, udd_old, problem_rhs);
this->A.usmv( -1.0, this->u_old, problem_rhs);
this->A.usmv(-this->tau/2.0, ud_old, problem_rhs);
// For fixed tau, we'd only really have to do this once
problem_A = this->A;
problem_A *= this->tau/2.0;
Arithmetic::addProduct(problem_A, 2.0/this->tau, B);
// ud_old makes a good initial iterate; we could use anything, though
problem_iterate = ud_old;
for (size_t i=0; i < this->dirichletNodes.size(); ++i)
if (this->dirichletNodes[i].count() == dim) {
problem_iterate[i] = 0;
this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
} else if (this->dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction prescribed
}
// u1 = u0 + tau/2 (du1 + du0)
void virtual extractSolution(VectorType const &problem_iterate,
VectorType &solution) const
{
solution = this->u_old;
Arithmetic::addProduct(solution, this->tau / 2.0, problem_iterate); // ud
Arithmetic::addProduct(solution, this->tau / 2.0, ud_old);
}
void virtual extractVelocity(VectorType const &problem_iterate,
VectorType &velocity) const
{
velocity = problem_iterate;
}
private:
MatrixType const &B;
VectorType const &ud_old;
VectorType const &udd_old;
};
#+end_src
* C++ code generation
#+name: timestepping
#+begin_src c++ :tangle timestepping.cc :noweb yes
<<includes>>
<<TimeSteppingScheme>>
<<ImplicitEuler>>
<<ImplicitTwoStep>>
<<Newmark>>
#+end_src
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment