[Problem] Damping

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,

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);
   }