From baae6ebd5e7e20fd006928ea42e5bc3766ae2b3b Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Wed, 17 Jul 2013 23:24:04 +0200
Subject: [PATCH] [Problem] Damping
---
src/one-body-sample.cc | 10 ++++++++--
src/timestepping/newmark.cc | 18 +++++++++++-------
2 files changed, 19 insertions(+), 9 deletions(-)
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 153f5d3c..b78b4bfb 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -399,16 +399,22 @@ int main(int argc, char *argv[]) {
VectorType a_initial(finestSize);
a_initial = 0.0;
{
- // We solve Au + Ma + Psi(v) = ell, thus Ma = - (Au + Psi(v) - ell)
+ /* We solve Au + Cv + Ma + Psi(v) = ell, thus
+ Ma = - (Au + Cv + Psi(v) - ell)
+ */
+ MatrixType const &C = A;
+ double const wc = 0.0; // watch out -- do not choose 1.0 here!
VectorType problem_rhs_initial(finestSize);
{
problem_rhs_initial = 0.0;
Arithmetic::addProduct(problem_rhs_initial, A, u_initial);
+ Arithmetic::addProduct(problem_rhs_initial, wc, C, v_initial);
myGlobalNonlinearity->updateState(alpha_initial);
// NOTE: We assume differentiability of Psi at v0 here!
myGlobalNonlinearity->addGradient(v_initial, problem_rhs_initial);
problem_rhs_initial -= ell;
- problem_rhs_initial *= -1.0;
+ // instead of multiplying M by (1.0 - wc), we divide the RHS
+ problem_rhs_initial *= -1.0 / (1.0 - wc);
}
LinearFactoryType accelerationFactory(parset.sub("solver.tnnmg"), // FIXME
refinements, 1e-12, // FIXME,
diff --git a/src/timestepping/newmark.cc b/src/timestepping/newmark.cc
index 22d8dcdc..aa1fad7b 100644
--- a/src/timestepping/newmark.cc
+++ b/src/timestepping/newmark.cc
@@ -27,9 +27,11 @@ void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
tau = _tau;
+ MatrixType const &C = A;
+ double const wc = 0.0;
/* We start out with the formulation
- M a + A u = ell
+ M a + C v + A u = ell
Newmark means
@@ -38,13 +40,13 @@ void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
in summary, we get at time t=1
- M [2/tau ( u1 - u0 ) - a0]
+ M [2/tau ( u1 - u0 ) - a0] + C v1
+ A [tau/2 ( v1 + v0 ) + u0] = ell
or
- 2/tau M v1 + tau/2 A v1
- = [2/tau M + tau/2 A] v1
+ 2/tau M v1 + C v1 + tau/2 A v1
+ = [2/tau M + C + tau/2 A] v1
= ell + 2/tau M v0 + M a0
- tau/2 A v0 - A u0
*/
@@ -54,17 +56,19 @@ void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
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, 2.0 / tau, M);
+ Arithmetic::addProduct(problem_AM, (1.0 - wc) * 2.0 / tau, M);
+ Arithmetic::addProduct(problem_AM, wc, C);
Arithmetic::addProduct(problem_AM, tau / 2.0, A);
// set up RHS
{
problem_rhs = ell;
- Arithmetic::addProduct(problem_rhs, 2.0 / tau, M, v_o);
- Arithmetic::addProduct(problem_rhs, M, a_o);
+ Arithmetic::addProduct(problem_rhs, (1.0 - wc) * 2.0 / tau, M, v_o);
+ Arithmetic::addProduct(problem_rhs, 1.0 - wc, M, a_o);
Arithmetic::subtractProduct(problem_rhs, tau / 2.0, A, v_o);
Arithmetic::subtractProduct(problem_rhs, A, u_o);
}
--
GitLab