From 1094fd78a5d95fbc650b3196fe7aefa9cc5633d5 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Wed, 5 Sep 2012 16:11:54 +0200
Subject: [PATCH] Add a two-step implicit time stepping scheme
---
src/one-body-sample.cc | 174 +++++++++++++++++++++++++++++++++----
src/one-body-sample.parset | 2 +
src/one-body-sample.py | 12 +--
3 files changed, 162 insertions(+), 26 deletions(-)
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index a39f1f2b..aa9e475a 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -101,6 +101,119 @@ void setup_boundary(GridView const &gridView,
}
}
+// Implicit Euler: Solve the problem
+//
+// a(Delta u_new, v - Delta u_new) + j_tau(v) - j_tau(Delta u_new)
+// >= l(w - Delta u_new) - a(u_new, v - Delta u_new)
+//
+// Setup: Substract a(u_new, .) from rhs
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void implicitEulerSetup(VectorType const &ell, MatrixType const &A,
+ VectorType const &u_old, VectorType const *u_old_old,
+ VectorType &problem_rhs, VectorType &problem_iterate,
+ MatrixType &problem_A,
+ Dune::BitSetVector<dim> const &dirichletNodes,
+ FunctionType const &dirichletFunction, double time) {
+ problem_A = A;
+ problem_rhs = ell;
+ problem_A.mmv(u_old, problem_rhs);
+
+ problem_iterate = u_old;
+ if (u_old_old)
+ problem_iterate -= *u_old_old;
+
+ for (size_t i = 0; i < dirichletNodes.size(); ++i)
+ if (dirichletNodes[i].count() == dim) {
+ double val;
+ dirichletFunction.evaluate(time, val);
+ problem_iterate[i] = 0; // Everything prescribed
+ problem_iterate[i][0] = val - u_old[i][0]; // Time-dependent X direction
+ } else if (dirichletNodes[i][1])
+ problem_iterate[i][1] = 0; // Y direction described
+}
+
+// Extraction: Add previous solution
+template <class VectorType>
+void implicitEulerExtract(VectorType const &u_old, VectorType const *u_old_old,
+ VectorType const &problem_iterate,
+ VectorType &solution) {
+ solution = u_old;
+ solution += problem_iterate;
+}
+
+template <class VectorType>
+void implicitEulerExtractVelocity(VectorType const &u_old,
+ VectorType const *u_old_old,
+ VectorType const &problem_iterate,
+ VectorType &diff) {
+ diff = problem_iterate;
+}
+
+// two-Stage implicit algorithm: Solve the problem
+//
+// a(Delta u_new, v - Delta u_new) + j_tau(v) - j_tau(Delta u_new)
+// >= l(w - Delta u_new) - a(u_new, v - Delta u_new)
+//
+// Setup: Substract a(u_new, .) from rhs
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void twoStageImplicitSetup(VectorType const &ell, MatrixType const &A,
+ VectorType const &u_old, VectorType const *u_old_old,
+ VectorType &problem_rhs, VectorType &problem_iterate,
+ MatrixType &problem_A,
+ Dune::BitSetVector<dim> const &dirichletNodes,
+ FunctionType const &dirichletFunction, double time) {
+ problem_A = A;
+ problem_A /= 1.5;
+ problem_rhs = ell;
+ problem_A.usmv(-2, u_old, problem_rhs);
+ problem_A.usmv(.5, *u_old_old, problem_rhs);
+
+ // The finite difference makes a good start
+ problem_iterate = u_old;
+ problem_iterate -= *u_old_old;
+
+ for (size_t i = 0; i < dirichletNodes.size(); ++i)
+ if (dirichletNodes[i].count() == dim) {
+ double val;
+ dirichletFunction.evaluate(time, val);
+ problem_iterate[i] = 0;
+ problem_iterate[i].axpy(-2, u_old[i]);
+ problem_iterate[i].axpy(.5, (*u_old_old)[i]);
+ problem_iterate[i][0] =
+ 1.5 * val - 2 * u_old[i][0] + .5 * (*u_old_old)[i][0];
+ } else if (dirichletNodes[i][1])
+ // Y direction described
+ problem_iterate[i][1] = -2 * u_old[i][1] + .5 * (*u_old_old)[i][1];
+}
+
+template <class VectorType>
+void twoStageImplicitExtract(VectorType const &u_old,
+ VectorType const *u_old_old,
+ VectorType const &problem_iterate,
+ VectorType &solution) {
+ solution = problem_iterate;
+ solution.axpy(2, u_old);
+ solution.axpy(-.5, *u_old_old);
+ solution *= 2.0 / 3.0;
+
+ // Check if we split correctly
+ {
+ VectorType test = problem_iterate;
+ test.axpy(-1.5, solution);
+ test.axpy(+2, u_old);
+ test.axpy(-.5, *u_old_old);
+ assert(test.two_norm() < 1e-10);
+ }
+}
+
+template <class VectorType>
+void twoStageImplicitExtractVelocity(VectorType const &u_old,
+ VectorType const *u_old_old,
+ VectorType const &problem_iterate,
+ VectorType &diff) {
+ diff = problem_iterate;
+}
+
int main(int argc, char *argv[]) {
try {
typedef SharedPointerMap<std::string, Dune::VirtualFunction<double, double>>
@@ -191,15 +304,18 @@ int main(int argc, char *argv[]) {
// {{{ Initialise vectors
VectorType u(finestSize);
u = 0.0; // Has to be zero!
+ VectorType u_previous;
VectorType u_diff(finestSize);
u_diff = 0.0; // Has to be zero!
+ // temporary storage for u
+ VectorType solution(finestSize);
+
SingletonVectorType alpha_old(finestSize);
alpha_old = parset.get<double>("boundary.friction.state.initial");
SingletonVectorType alpha(alpha_old);
SingletonVectorType vonMisesStress;
- VectorType rhs;
// }}}
typedef MyConvexProblem<MatrixType, VectorType> MyConvexProblemType;
@@ -244,36 +360,61 @@ int main(int argc, char *argv[]) {
for (size_t run = 1; run <= timesteps; ++run) {
double const time = tau * run;
{
+ VectorType ell(finestSize);
assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
- leafView, p1Basis, neumannNodes, rhs, neumannFunction, time);
- stiffnessMatrix.mmv(u, rhs);
- // Apply Dirichlet condition
- for (size_t i = 0; i < finestSize; ++i)
- if (ignoreNodes[i].count() == dim) {
- dirichletFunction.evaluate(time, u_diff[i][0]);
- u_diff[i][0] /= timesteps;
- }
+ leafView, p1Basis, neumannNodes, ell, neumannFunction, time);
+
+ MatrixType problem_A;
+ VectorType problem_rhs(finestSize);
+ VectorType problem_iterate(finestSize);
+
+ auto setupFunc =
+ (run == 1 || !parset.get<bool>("twoStageImplicit"))
+ ? &implicitEulerSetup<VectorType, MatrixType,
+ decltype(dirichletFunction), dim>
+ : &twoStageImplicitSetup<VectorType, MatrixType,
+ decltype(dirichletFunction), dim>;
+ VectorType *u_old_old_ptr = (run == 1) ? nullptr : &u_previous;
+ setupFunc(ell, stiffnessMatrix, u, u_old_old_ptr, problem_rhs,
+ problem_iterate, problem_A, ignoreNodes, dirichletFunction,
+ time);
+
+ VectorType solution_saved = u;
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) {
auto myGlobalNonlinearity =
assemble_nonlinearity<MatrixType, VectorType>(
parset.sub("boundary.friction"), *nodalIntegrals, alpha, tau);
- MyConvexProblemType const myConvexProblem(stiffnessMatrix,
- *myGlobalNonlinearity, rhs);
+ MyConvexProblemType const myConvexProblem(
+ problem_A, *myGlobalNonlinearity, problem_rhs);
MyBlockProblemType myBlockProblem(parset, myConvexProblem);
auto multigridStep = mySolver.getSolver();
- multigridStep->setProblem(u_diff, myBlockProblem);
+ multigridStep->setProblem(problem_iterate, myBlockProblem);
- VectorType const u_diff_saved = u_diff;
LoopSolver<VectorType> overallSolver(
multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
solver_tolerance, &energyNorm, verbosity,
false); // absolute error
overallSolver.solve();
+ auto extractFunc = (run == 1 || !parset.get<bool>("twoStageImplicit"))
+ ? implicitEulerExtract<VectorType>
+ : twoStageImplicitExtract<VectorType>;
+
+ // Extract solution from solver
+ extractFunc(u, u_old_old_ptr, problem_iterate, solution);
+
+ auto extractDiffFunc =
+ (run == 1 || !parset.get<bool>("twoStageImplicit"))
+ ? implicitEulerExtractVelocity<VectorType>
+ : twoStageImplicitExtractVelocity<VectorType>;
+
+ // Extract difference from solver
+ extractDiffFunc(u, u_old_old_ptr, problem_iterate, u_diff);
+
// Update the state
for (size_t i = 0; i < frictionalNodes.size(); ++i) {
if (frictionalNodes[i][0]) {
@@ -299,9 +440,11 @@ int main(int argc, char *argv[]) {
std::cerr << '.';
std::cerr.flush();
}
- if (energyNorm.diff(u_diff_saved, u_diff) <
+ if (energyNorm.diff(solution_saved, solution) <
parset.get<double>("solver.tnnmg.fixed_point_tolerance"))
break;
+ else
+ solution_saved = solution;
if (state_fpi == state_fpi_max)
std::cerr << "[ref = " << refinements
@@ -363,7 +506,8 @@ int main(int argc, char *argv[]) {
std::cout << u_diff[i][0] * timesteps << " ";
std::cout << std::endl;
}
- u += u_diff;
+ u_previous = u;
+ u = solution;
alpha_old = alpha;
// Compute von Mises stress and write everything to a file
diff --git a/src/one-body-sample.parset b/src/one-body-sample.parset
index 63bfca90..3a623491 100644
--- a/src/one-body-sample.parset
+++ b/src/one-body-sample.parset
@@ -13,6 +13,8 @@ printVelocitySteppingComparison = false
enable_timer = false
+twoStageImplicit = false
+
[grid]
refinements = 4
diff --git a/src/one-body-sample.py b/src/one-body-sample.py
index e6db78f3..77b2f491 100644
--- a/src/one-body-sample.py
+++ b/src/one-body-sample.py
@@ -13,17 +13,7 @@ class neumannCondition:
class dirichletCondition:
def __call__(self, x):
- return .005
- # return 0
- fst = 3e-1
- snd = 5e-1
- trd = 3e-1
- if x < 1.0/5:
- return fst
- elif x < 3.0/5:
- return snd
- else:
- return trd
+ return x * .005
Functions = {
'neumannCondition' : neumannCondition(),
--
GitLab