diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index e2072e521bb8ec7eada87580c954cd4caeb22d1c..0e0ce17a09ce6971d47c9280077ff922005083f1 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -242,19 +242,15 @@ int main(int argc, char *argv[]) {
// {{{ Initialise vectors
VectorType u(finestSize);
u = 0.0;
- VectorType u_old(finestSize);
- u_old = 0.0; // Has to be zero!
- VectorType u_old_old(finestSize);
-
VectorType ud(finestSize);
ud = 0.0;
- VectorType ud_old(finestSize);
- ud_old = 0.0;
- VectorType udd(finestSize);
- udd = 0.0;
- VectorType udd_old(finestSize);
- udd_old = 0.0;
+ VectorType u_initial(finestSize);
+ u_initial = 0.0;
+ VectorType ud_initial(finestSize);
+ ud_initial = 0.0;
+ VectorType udd_initial(finestSize);
+ udd_initial = 0.0;
SingletonVectorType alpha_old(finestSize);
alpha_old = parset.get<double>("boundary.friction.state.initial");
@@ -302,6 +298,30 @@ int main(int argc, char *argv[]) {
first_frictional_node < frictionalNodes.size())
++first_frictional_node;
+ Dune::shared_ptr<TimeSteppingScheme<VectorType, MatrixType,
+ decltype(dirichletFunction), dim>>
+ timeSteppingScheme;
+ switch (parset.get<Config::scheme>("timeSteppingScheme")) {
+ case Config::ImplicitTwoStep:
+ timeSteppingScheme = Dune::make_shared<ImplicitTwoStep<
+ VectorType, MatrixType, decltype(dirichletFunction), dim>>(
+ stiffnessMatrix, u_initial, ud_initial, ignoreNodes,
+ dirichletFunction);
+ break;
+ case Config::ImplicitEuler:
+ timeSteppingScheme = Dune::make_shared<ImplicitEuler<
+ VectorType, MatrixType, decltype(dirichletFunction), dim>>(
+ stiffnessMatrix, u_initial, ud_initial, ignoreNodes,
+ dirichletFunction);
+ break;
+ case Config::Newmark:
+ timeSteppingScheme = Dune::make_shared<
+ Newmark<VectorType, MatrixType, decltype(dirichletFunction), dim>>(
+ stiffnessMatrix, massMatrix, u_initial, ud_initial, udd_initial,
+ ignoreNodes, dirichletFunction);
+ break;
+ }
+
for (size_t run = 1; run <= timesteps; ++run) {
double const time = tau * run;
{
@@ -310,43 +330,16 @@ int main(int argc, char *argv[]) {
leafView, p1Basis, neumannNodes, ell, neumannFunction, time);
ell += gravityFunctional;
+ // Copy everything for now
MatrixType problem_A;
VectorType problem_rhs(finestSize);
VectorType problem_iterate(finestSize);
- VectorType *u_old_old_ptr = (run == 1) ? nullptr : &u_old_old;
-
- Dune::shared_ptr<TimeSteppingScheme<VectorType, MatrixType,
- decltype(dirichletFunction), dim>>
- timeSteppingScheme;
- {
- switch (parset.get<Config::scheme>("timeSteppingScheme")) {
- case Config::ImplicitTwoStep:
- if (run != 1) {
- timeSteppingScheme = Dune::make_shared<ImplicitTwoStep<
- VectorType, MatrixType, decltype(dirichletFunction), dim>>(
- ell, stiffnessMatrix, u_old, u_old_old_ptr, ignoreNodes,
- dirichletFunction, time, tau);
- break;
- }
- // Fall through
- case Config::ImplicitEuler:
- timeSteppingScheme = Dune::make_shared<ImplicitEuler<
- VectorType, MatrixType, decltype(dirichletFunction), dim>>(
- ell, stiffnessMatrix, u_old, u_old_old_ptr, ignoreNodes,
- dirichletFunction, time, tau);
- break;
- case Config::Newmark:
- timeSteppingScheme = Dune::make_shared<Newmark<
- VectorType, MatrixType, decltype(dirichletFunction), dim>>(
- ell, stiffnessMatrix, massMatrix, u_old, ud_old, udd_old,
- ignoreNodes, dirichletFunction, time, tau);
- break;
- }
- }
- timeSteppingScheme->setup(problem_rhs, problem_iterate, problem_A);
+ timeSteppingScheme->setup(ell, tau, time, problem_rhs, problem_iterate,
+ problem_A);
- VectorType ud_saved = ud_old;
+ VectorType ud_saved(finestSize);
+ ud_saved = 0.0; // FIXME
auto const state_fpi_max =
parset.get<size_t>("solver.tnnmg.fixed_point_iterations");
for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) {
@@ -366,13 +359,9 @@ int main(int argc, char *argv[]) {
false); // absolute error
overallSolver.solve();
- timeSteppingScheme->extractDisplacement(problem_iterate, u);
- timeSteppingScheme->extractVelocity(problem_iterate, ud);
-
- udd = 0;
- Arithmetic::addProduct(udd, 2.0 / tau, ud);
- Arithmetic::addProduct(udd, -2.0 / tau, ud_old);
- Arithmetic::addProduct(udd, -1.0, udd_old);
+ timeSteppingScheme->postProcess(problem_iterate);
+ timeSteppingScheme->extractDisplacement(u);
+ timeSteppingScheme->extractVelocity(ud);
// Update the state
for (size_t i = 0; i < frictionalNodes.size(); ++i) {
@@ -465,10 +454,6 @@ int main(int argc, char *argv[]) {
std::cout << ud[i][0] << " ";
std::cout << std::endl;
}
- u_old_old = u_old;
- u_old = u;
- ud_old = ud;
- udd_old = udd;
alpha_old = alpha;
// Compute von Mises stress and write everything to a file
diff --git a/src/timestepping.org b/src/timestepping.org
index 92644bd67a63b681aea0f56cf0e76db7e3da1c31..911af73eb77591f5a6d485d83e307d875d8a9011 100644
--- a/src/timestepping.org
+++ b/src/timestepping.org
@@ -17,43 +17,22 @@
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
+ setup(VectorType const &ell,
+ double _tau,
+ double time,
+ VectorType &problem_rhs,
+ VectorType &problem_iterate,
+ MatrixType &problem_A) = 0;
- void virtual extractDisplacement(VectorType const &problem_iterate,
- VectorType &solution) const = 0;
+ void virtual
+ postProcess(VectorType const &problem_iterate) = 0;
- void virtual extractVelocity(VectorType const &problem_iterate,
- VectorType &velocity) const = 0;
+ void virtual
+ extractDisplacement(VectorType &displacement) const = 0;
- protected:
- VectorType const ℓ
- MatrixType const &A;
- VectorType const &u_old;
- Dune::BitSetVector<dim> const &dirichletNodes;
- FunctionType const &dirichletFunction;
- double const time;
- double const tau;
+ void virtual
+ extractVelocity(VectorType &velocity) const = 0;
};
#+end_src
@@ -64,61 +43,91 @@
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
+ ImplicitEuler(MatrixType const &_A,
+ VectorType const &_u_initial,
+ VectorType const &_ud_initial,
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)
+ FunctionType const &_dirichletFunction)
+ : A(_A),
+ u(_u_initial),
+ ud(_ud_initial),
+ dirichletNodes(_dirichletNodes),
+ dirichletFunction(_dirichletFunction)
{}
void virtual
- setup(VectorType &problem_rhs,
+ setup(VectorType const &ell,
+ double _tau,
+ double time,
+ VectorType &problem_rhs,
VectorType &problem_iterate,
- MatrixType &problem_A) const
+ MatrixType &problem_A)
{
- 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) {
+ 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;
- this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
- } else if (this->dirichletNodes[i][1])
+ dirichletFunction.evaluate(time, problem_iterate[i][0]);
+ } else if (dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction prescribed
}
- void virtual extractDisplacement(VectorType const &problem_iterate,
- VectorType &solution) const
+ 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
{
- solution = this->u_old;
- solution.axpy(this->tau, problem_iterate);
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ displacement = u;
}
- void virtual extractVelocity(VectorType const &problem_iterate,
- VectorType &velocity) const
+ void virtual extractVelocity(VectorType &velocity) const
{
- velocity = problem_iterate;
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ velocity = ud;
}
private:
- VectorType const *u_old_old_ptr;
+ 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;
};
#+end_src
@@ -152,71 +161,131 @@
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,
+ ImplicitTwoStep(MatrixType const &_A,
+ VectorType const &_u_initial,
+ VectorType const &_ud_initial,
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)
+ FunctionType const &_dirichletFunction)
+ : A(_A),
+ u(_u_initial),
+ ud(_ud_initial),
+ dirichletNodes(_dirichletNodes),
+ dirichletFunction(_dirichletFunction)
{}
- void virtual setup(VectorType &problem_rhs,
- VectorType &problem_iterate,
- MatrixType &problem_A) const
+ // FIXME: handle case run == 1
+ void virtual
+ setup(VectorType const &ell,
+ double _tau,
+ double time,
+ VectorType &problem_rhs,
+ VectorType &problem_iterate,
+ MatrixType &problem_A)
{
- 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);
+ postProcessCalled = false;
+
+ tau = _tau;
- problem_A = this->A;
- problem_A *= 2.0/3.0 * this->tau;
+ // This function is only called once for every timestep
+ u_old_old = u_old;
+ ud_old = ud;
+ u_old = u;
- // Use the old velocity as an initial iterate
- problem_iterate = this->u_old;
- problem_iterate -= *u_old_old_ptr;
- problem_iterate /= this->tau;
+ switch (state) {
+ // Perform an implicit Euler step since we lack information
+ case NO_SETUP:
+ state = FIRST_SETUP;
- for (size_t i=0; i < this->dirichletNodes.size(); ++i)
- if (this->dirichletNodes[i].count() == dim) {
+ 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;
- this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
- } else if (this->dirichletNodes[i][1])
+ dirichletFunction.evaluate(time, problem_iterate[i][0]);
+ } else if (dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction prescribed
}
- void virtual extractDisplacement(VectorType const &problem_iterate,
- VectorType &solution) const
+ void virtual postProcess(VectorType const &problem_iterate)
{
- 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);
+ 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 extractVelocity(VectorType const &problem_iterate,
- VectorType &velocity) const
+ 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
{
- velocity = problem_iterate;
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ velocity = ud;
}
private:
- VectorType const *u_old_old_ptr;
+ 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;
};
#+end_src
@@ -261,70 +330,111 @@
class Newmark : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
{
public:
- Newmark(VectorType const &_ell,
- MatrixType const &_A,
+ Newmark(MatrixType const &_A,
MatrixType const &_B,
- VectorType const &_u_old,
- VectorType const &_ud_old,
- VectorType const &_udd_old,
+ VectorType const &_u_initial,
+ VectorType const &_ud_initial,
+ VectorType const &_udd_initial,
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)
+ FunctionType const &_dirichletFunction)
+ : A(_A),
+ B(_B),
+ u(_u_initial),
+ ud(_ud_initial),
+ udd(_udd_initial),
+ dirichletNodes(_dirichletNodes),
+ dirichletFunction(_dirichletFunction)
{}
void virtual
- setup(VectorType &problem_rhs,
+ setup(VectorType const &ell,
+ double _tau,
+ double time,
+ VectorType &problem_rhs,
VectorType &problem_iterate,
- MatrixType &problem_A) const
+ MatrixType &problem_A)
{
- 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);
+ 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 = this->A;
- problem_A *= this->tau/2.0;
- Arithmetic::addProduct(problem_A, 2.0/this->tau, B);
+ 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 < this->dirichletNodes.size(); ++i)
- if (this->dirichletNodes[i].count() == dim) {
+ for (size_t i=0; i < dirichletNodes.size(); ++i)
+ if (dirichletNodes[i].count() == dim) {
problem_iterate[i] = 0;
- this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
- } else if (this->dirichletNodes[i][1])
+ dirichletFunction.evaluate(time, problem_iterate[i][0]);
+ } else if (dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction prescribed
}
- // u1 = u0 + tau/2 (du1 + du0)
- void virtual extractDisplacement(VectorType const &problem_iterate,
- VectorType &solution) const
+ 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
{
- solution = this->u_old;
- Arithmetic::addProduct(solution, this->tau / 2.0, problem_iterate); // ud
- Arithmetic::addProduct(solution, this->tau / 2.0, ud_old);
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ displacement = u;
}
- void virtual extractVelocity(VectorType const &problem_iterate,
- VectorType &velocity) const
+ void virtual extractVelocity(VectorType &velocity) const
{
- velocity = problem_iterate;
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ velocity = ud;
}
private:
+ MatrixType const &A;
MatrixType const &B;
- VectorType const &ud_old;
- VectorType const &udd_old;
+ 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;
};
+
#+end_src
* C++ code generation