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

Compile time stepping separately

parent cf03bdab
No related branches found
No related tags found
No related merge requests found
timestepping.cc
timestepping.hh
timestepping_tmpl.cc
......@@ -10,6 +10,7 @@ SOURCES = \
compute_state_ruina.cc \
mysolver.cc \
one-body-sample.cc \
timestepping.cc \
vtk.cc
## 2D
......@@ -101,7 +102,7 @@ DISTCHECK_CONFIGURE_FLAGS = \
include $(top_srcdir)/am/global-rules
BUILT_SOURCES = timestepping.cc
BUILT_SOURCES = timestepping.hh timestepping.cc
$(srcdir)/timestepping.cc: $(srcdir)/timestepping.org
emacs -Q --batch --eval \
"(let (vc-handled-backends) \
......
......@@ -75,7 +75,7 @@
#include "enum_state_model.cc"
#include "enum_scheme.cc"
#include "timestepping.cc"
#include "timestepping.hh"
int const dim = DIM;
......
#+name: includes
* Preamble
#+name: includes.hh
#+begin_src c++
#include <dune/common/bitsetvector.hh>
#+end_src
#+name: preamble.cc
#+begin_src c++
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/fufem/arithmetic.hh>
#include "timestepping.hh"
#+end_src
#+name: preamble
#+name: preamble.tex
#+begin_src latex
\documentclass{scrartcl}
\usepackage{amsmath}
\begin{document}
#+end_src
#+name: virtual_function_declaration
#+begin_src c++
void virtual setup(VectorType const &, double, double, VectorType &,
VectorType &, MatrixType &);
void virtual postProcess(VectorType const &);
void virtual extractDisplacement(VectorType &) const;
void virtual extractVelocity(VectorType &) const;
#+end_src
* Abstract TimeSteppingScheme
#+name: TimeSteppingScheme
#+name: TimeSteppingScheme.hh
#+begin_src c++
template <class VectorType, class MatrixType, class FunctionType, int dim>
class TimeSteppingScheme
......@@ -25,20 +41,14 @@
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;
void virtual postProcess(VectorType const &problem_iterate) = 0;
void virtual extractDisplacement(VectorType &displacement) const = 0;
void virtual extractVelocity(VectorType &velocity) const = 0;
};
#+end_src
* TimeSteppingScheme: Implicit Euler
#+name: ImplicitEuler
#+begin_src c++
#+name: ImplicitEuler.hh
#+begin_src c++ noweb: yes
template <class VectorType, class MatrixType, class FunctionType, int dim>
class ImplicitEuler : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
{
......@@ -47,90 +57,109 @@
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)
{}
FunctionType const &_dirichletFunction);
<<virtual_function_declaration>>
void virtual
setup(VectorType const &ell,
double _tau,
double time,
VectorType &problem_rhs,
VectorType &problem_iterate,
MatrixType &problem_A)
{
postProcessCalled = false;
private:
MatrixType const &A;
VectorType u;
VectorType ud;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
tau = _tau;
VectorType u_old;
VectorType ud_old;
// This function is only called once for every timestep
ud_old = ud;
u_old = u;
double tau;
problem_rhs = ell;
A.mmv(u_old, problem_rhs);
bool postProcessCalled = false;
};
#+end_src
#+name: ImplicitEuler.cc
#+begin_src c++
template <class VectorType, class MatrixType, class FunctionType, int dim>
ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
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)
{}
// For fixed tau, we'd only really have to do this once
problem_A = A;
problem_A *= tau;
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
setup(VectorType const &ell,
double _tau,
double time,
VectorType &problem_rhs,
VectorType &problem_iterate,
MatrixType &problem_A)
{
postProcessCalled = false;
// ud_old makes a good initial iterate; we could use anything, though
problem_iterate = ud_old;
tau = _tau;
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
}
// This function is only called once for every timestep
ud_old = ud;
u_old = u;
void virtual postProcess(VectorType const &problem_iterate)
{
postProcessCalled = true;
problem_rhs = ell;
A.mmv(u_old, problem_rhs);
ud = problem_iterate;
// For fixed tau, we'd only really have to do this once
problem_A = A;
problem_A *= tau;
u = u_old;
Arithmetic::addProduct(u, tau, ud);
}
// ud_old makes a good initial iterate; we could use anything, though
problem_iterate = ud_old;
void virtual extractDisplacement(VectorType &displacement) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
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
}
displacement = u;
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
postProcess(VectorType const &problem_iterate)
{
postProcessCalled = true;
void virtual extractVelocity(VectorType &velocity) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
ud = problem_iterate;
velocity = ud;
}
u = u_old;
Arithmetic::addProduct(u, tau, ud);
}
private:
MatrixType const &A;
VectorType u;
VectorType ud;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
extractDisplacement(VectorType &displacement) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
VectorType u_old;
VectorType ud_old;
displacement = u;
}
double tau;
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
extractVelocity(VectorType &velocity) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = ud;
}
bool postProcessCalled = false;
};
#+end_src
* TimeSteppingScheme: Implicit Two-Step
#+begin_src latex :tangle twostep.tex :noweb yes
\documentclass{scrartcl}
......@@ -154,9 +183,8 @@
\end{equation*}
\end{document}
#+end_src
#+name: ImplicitTwoStep
#+begin_src c++
#+name: ImplicitTwoStep.hh
#+begin_src c++ noweb: yes
template <class VectorType, class MatrixType, class FunctionType, int dim>
class ImplicitTwoStep : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
{
......@@ -165,108 +193,9 @@
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;
}
FunctionType const &_dirichletFunction);
<<virtual_function_declaration>>
private:
MatrixType const &A;
......@@ -288,10 +217,126 @@
state_type state = NO_SETUP;
};
#+end_src
#+name: ImplicitTwoStep.cc
#+begin_src c++
template <class VectorType, class MatrixType, class FunctionType, int dim>
ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
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)
{}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
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
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
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);
}
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
extractDisplacement(VectorType &displacement) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
displacement = u;
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
extractVelocity(VectorType &velocity) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = ud;
}
#+end_src
* TimeSteppingScheme: Newmark
#+begin_src latex :tangle newmark.tex :noweb yes
<<preamble>>
<<preamble.tex>>
\noindent The Newmark scheme in its classical form with $\gamma = 1/2$
and $\beta = 1/4$ reads
\begin{align}
......@@ -323,9 +368,8 @@
from \eqref{eq:4}. The only unknown is now $\dot u_1$.
\end{document}
#+end_src
#+name: Newmark
#+begin_src c++
#+name: Newmark.hh
#+begin_src c++ noweb: yes
template <class VectorType, class MatrixType, class FunctionType, int dim>
class Newmark : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
{
......@@ -336,86 +380,9 @@
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)
{}
FunctionType const &_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;
}
<<virtual_function_declaration>>
private:
MatrixType const &A;
......@@ -434,16 +401,148 @@
bool postProcessCalled = false;
};
#+end_src
#+name: Newmark.cc
#+begin_src c++
template <class VectorType, class MatrixType, class FunctionType, int dim>
Newmark<VectorType, MatrixType, FunctionType, dim>::
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)
{}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void Newmark<VectorType, MatrixType, FunctionType, dim>::
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
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void Newmark<VectorType, MatrixType, FunctionType, dim>::
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);
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void Newmark<VectorType, MatrixType, FunctionType, dim>::
extractDisplacement(VectorType &displacement) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
displacement = u;
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void Newmark<VectorType, MatrixType, FunctionType, dim>::
extractVelocity(VectorType &velocity) const
{
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = ud;
}
#+end_src
* C++ code generation
#+name: timestepping
#+begin_src c++ :tangle timestepping.hh :noweb yes
// GENERATED -- DO NOT MODIFY
#ifndef DUNE_TECTONIC_TIMESTEPPING_HH
#define DUNE_TECTONIC_TIMESTEPPING_HH
<<includes.hh>>
<<TimeSteppingScheme.hh>>
<<ImplicitEuler.hh>>
<<ImplicitTwoStep.hh>>
<<Newmark.hh>>
#endif
#+end_src
#+begin_src c++ :tangle timestepping.cc :noweb yes
<<includes>>
// GENERATED -- DO NOT MODIFY
<<preamble.cc>>
<<TimeSteppingScheme>>
<<ImplicitEuler>>
<<ImplicitTwoStep>>
<<Newmark>>
<<ImplicitEuler.cc>>
<<ImplicitTwoStep.cc>>
<<Newmark.cc>>
#include "timestepping_tmpl.cc"
#+end_src
#+begin_src c++ :tangle timestepping_tmpl.cc :noweb yes
// GENERATED -- DO NOT MODIFY
#ifndef DIM
#error DIM unset
#endif
#include <dune/common/fmatrix.hh>
#include <dune/common/function.hh>
#include <dune/common/fvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
typedef Dune::FieldVector<double, DIM> SmallVector;
typedef Dune::FieldMatrix<double, DIM, DIM> SmallMatrix;
typedef Dune::BCRSMatrix<SmallMatrix> MatrixType;
typedef Dune::BlockVector<SmallVector> VectorType;
typedef Dune::VirtualFunction<double, double> FunctionType;
template class ImplicitEuler<VectorType, MatrixType, FunctionType, DIM>;
template class ImplicitTwoStep<VectorType, MatrixType, FunctionType, DIM>;
template class Newmark<VectorType, MatrixType, FunctionType, DIM>;
#+end_src
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment