From d63d0b07e7676b9eff38d11a073687b1d8a087a5 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sun, 21 Jul 2013 13:16:05 +0200
Subject: [PATCH] [Algorit] Add damping to EulerPair
---
src/timestepping/eulerpair.cc | 61 ++++++++++++++++++++++++++---------
src/timestepping/eulerpair.hh | 4 +--
2 files changed, 47 insertions(+), 18 deletions(-)
diff --git a/src/timestepping/eulerpair.cc b/src/timestepping/eulerpair.cc
index 364e64fd..dfef35b8 100644
--- a/src/timestepping/eulerpair.cc
+++ b/src/timestepping/eulerpair.cc
@@ -1,12 +1,11 @@
-
template <class VectorType, class MatrixType, class FunctionType, int dim>
EulerPair<VectorType, MatrixType, FunctionType, dim>::EulerPair(
- MatrixType const &_A, MatrixType const &_B, VectorType const &_u_initial,
+ MatrixType const &_A, MatrixType const &_M, VectorType const &_u_initial,
VectorType const &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction)
: A(_A),
- B(_B),
+ M(_M),
u(_u_initial),
v(_v_initial),
dirichletNodes(_dirichletNodes),
@@ -21,26 +20,56 @@ void EulerPair<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
template <class VectorType, class MatrixType, class FunctionType, int dim>
void EulerPair<VectorType, MatrixType, FunctionType, dim>::setup(
VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
- VectorType &problem_iterate, MatrixType &problem_AB) {
+ VectorType &problem_iterate, MatrixType &problem_AM) {
postProcessCalled = false;
tau = _tau;
- problem_rhs = ell;
- Arithmetic::addProduct(problem_rhs, 1.0 / tau, B, v_o);
- Arithmetic::subtractProduct(problem_rhs, A, u_o);
+ MatrixType const &C = A;
+ double const wc = 0.0;
+ /* We start out with the formulation
+
+ M a + C v + A u = ell
+
+ Backward Euler means
+
+ a1 = 1.0/tau ( v1 - v0 )
+ u1 = tau v1 + u0
+
+ in summary, we get at time t=1
+
+ M [1.0/tau ( v1 - v0 )] + C v1
+ + A [tau v1 + u0] = ell
+
+ or
+
+ 1.0/tau M v1 + C v1 + tau A v1
+ = [1.0/tau M + C + tau A] v1
+ = ell + 1.0/tau M v0 - A u0
+ */
+
+ // set up LHS (for fixed tau, we'd only really have to do this once)
+ {
+ Dune::MatrixIndexSet indices(A.N(), A.M());
+ indices.import(A);
+ indices.import(M);
+ indices.import(C);
+ indices.exportIdx(problem_AM);
+ }
+ problem_AM = 0.0;
+ Arithmetic::addProduct(problem_AM, (1.0 - wc) / tau, M);
+ Arithmetic::addProduct(problem_AM, wc, C);
+ Arithmetic::addProduct(problem_AM, tau, A);
- // For fixed tau, we'd only really have to do this once
- Dune::MatrixIndexSet indices(A.N(), A.M());
- indices.import(A);
- indices.import(B);
- indices.exportIdx(problem_AB);
- problem_AB = 0.0;
- Arithmetic::addProduct(problem_AB, tau, A);
- Arithmetic::addProduct(problem_AB, 1.0 / tau, B);
+ // set up RHS
+ {
+ problem_rhs = ell;
+ Arithmetic::addProduct(problem_rhs, (1.0 - wc) / tau, M, v_o);
+ Arithmetic::subtractProduct(problem_rhs, A, u_o);
+ }
// v_o makes a good initial iterate; we could use anything, though
- problem_iterate = v_o;
+ problem_iterate = 0.0;
for (size_t i = 0; i < dirichletNodes.size(); ++i)
switch (dirichletNodes[i].count()) {
diff --git a/src/timestepping/eulerpair.hh b/src/timestepping/eulerpair.hh
index 289f98d2..105529e1 100644
--- a/src/timestepping/eulerpair.hh
+++ b/src/timestepping/eulerpair.hh
@@ -5,7 +5,7 @@ template <class VectorType, class MatrixType, class FunctionType, int dim>
class EulerPair
: public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
public:
- EulerPair(MatrixType const &_A, MatrixType const &_B,
+ EulerPair(MatrixType const &_A, MatrixType const &_M,
VectorType const &_u_initial, VectorType const &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction);
@@ -19,7 +19,7 @@ class EulerPair
private:
MatrixType const &A;
- MatrixType const &B;
+ MatrixType const &M;
VectorType u;
VectorType v;
Dune::BitSetVector<dim> const &dirichletNodes;
--
GitLab