Move tau from nonlinearity to stiffness matrix

dune/tectonic/globalruinanonlinearity.hh

@@ -25,18 +25,17 @@ class GlobalRuinaNonlinearity
 
   GlobalRuinaNonlinearity(dataref nodalIntegrals, dataref a, dataref mu,
                           dataref eta, dataref normalStress, dataref b,
-                          dataref state, dataref L, double h)
+                          dataref state, dataref L)
       : restrictions(nodalIntegrals.size()) {
     auto trivialNonlinearity = make_shared<LocalNonlinearity<dim> const>(
         make_shared<TrivialFunction const>());
     for (size_t i = 0; i < restrictions.size(); ++i) {
-      restrictions[i] =
-          nodalIntegrals[i] == 0
-              ? trivialNonlinearity
-              : make_shared<LocalNonlinearity<dim> const>(
-                    make_shared<RuinaFunction const>(
-                        nodalIntegrals[i], a[i], mu[i], eta[i], normalStress[i],
-                        b[i], state[i], L[i], h));
+      restrictions[i] = nodalIntegrals[i] == 0
+                            ? trivialNonlinearity
+                            : make_shared<LocalNonlinearity<dim> const>(
+                                  make_shared<RuinaFunction const>(
+                                      nodalIntegrals[i], a[i], mu[i], eta[i],
+                                      normalStress[i], b[i], state[i], L[i]));
     }
   }
 

dune/tectonic/nicefunction.hh

@@ -41,18 +41,17 @@ class NiceFunction {
 class RuinaFunction : public NiceFunction {
 public:
   RuinaFunction(double coefficient, double a, double mu, double eta,
-                double normalStress, double b, double state, double L, double h)
+                double normalStress, double b, double state, double L)
       : NiceFunction(),
         a(a),
         eta(eta),
-        h(h),
         coefficientProduct(coefficient * normalStress),
         K(mu + b * (state - std::log(eta * L))), // state is assumed to be
                                                  // logarithmic
         rho(std::exp(-K / a)) {}
 
   double virtual evaluate(double x) const {
-    double const arg = x / h * eta;
+    double const arg = x * eta;
     if (arg <= rho)
       return 0;
 
@@ -61,7 +60,7 @@ class RuinaFunction : public NiceFunction {
   }
 
   double virtual leftDifferential(double x) const {
-    double const arg = x / h * eta;
+    double const arg = x * eta;
     if (arg <= rho)
       return 0;
 
@@ -74,7 +73,7 @@ class RuinaFunction : public NiceFunction {
   }
 
   double virtual second_deriv(double x) const {
-    double const arg = x / h * eta;
+    double const arg = x * eta;
     if (arg <= rho)
       return 0;
 
@@ -82,7 +81,7 @@ class RuinaFunction : public NiceFunction {
   }
 
   double virtual regularity(double x) const {
-    double const arg = x / h * eta;
+    double const arg = x * eta;
     // TODO: Make this controllable
     if (std::abs(arg - rho) < 1e-14)
       return std::numeric_limits<double>::infinity();
@@ -95,7 +94,6 @@ class RuinaFunction : public NiceFunction {
 private:
   double const a;
   double const eta;
-  double const h;
   double const coefficientProduct;
   double const K;
   double const rho;

src/assemblers.cc

@@ -84,7 +84,7 @@ assemble_nonlinearity(
 
       return Dune::make_shared<
           Dune::GlobalRuinaNonlinearity<MatrixType, VectorType> const>(
-          nodalIntegrals, a, mu, eta, normalStress, b, state, L, tau);
+          nodalIntegrals, a, mu, eta, normalStress, b, state, L);
     }
     case Config::Laursen:
       assert(false);

src/one-body-sample.cc

@@ -270,11 +270,11 @@ int main(int argc, char *argv[]) {
         if (run == 1 || !parset.get<bool>("implicitTwoStep"))
           timeSteppingScheme =
               Dune::make_shared<IE>(ell, stiffnessMatrix, u, u_old_old,
-                                    ignoreNodes, dirichletFunction, time);
+                                    ignoreNodes, dirichletFunction, time, tau);
         else
           timeSteppingScheme =
               Dune::make_shared<ITS>(ell, stiffnessMatrix, u, u_old_old,
-                                     ignoreNodes, dirichletFunction, time);
+                                     ignoreNodes, dirichletFunction, time, tau);
 
         timeSteppingScheme->setup(problem_rhs, problem_iterate, problem_A);
 

src/timestepping.cc

@@ -4,14 +4,16 @@ class TimeSteppingScheme {
   TimeSteppingScheme(VectorType const &_ell, MatrixType const &_A,
                      VectorType const &_u_old, VectorType const *_u_old_old,
                      Dune::BitSetVector<dim> const &_dirichletNodes,
-                     FunctionType const &_dirichletFunction, double _time)
+                     FunctionType const &_dirichletFunction, double _time,
+                     double _tau)
       : ell(_ell),
         A(_A),
         u_old(_u_old),
         u_old_old(_u_old_old),
         dirichletNodes(_dirichletNodes),
         dirichletFunction(_dirichletFunction),
-        time(_time) {}
+        time(_time),
+        tau(_tau) {}
 
   virtual ~TimeSteppingScheme() {}
 
@@ -31,7 +33,8 @@ class TimeSteppingScheme {
   VectorType const *u_old_old;
   Dune::BitSetVector<dim> const &dirichletNodes;
   FunctionType const &dirichletFunction;
-  double time;
+  double const time;
+  double const tau;
 };
 
 // Implicit Euler: Solve the problem
@@ -50,10 +53,13 @@ class ImplicitEuler
 
   void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
                      MatrixType &problem_A) const {
-    problem_A = this->A;
     problem_rhs = this->ell;
-    problem_A.mmv(this->u_old, problem_rhs);
+    this->A.mmv(this->u_old, problem_rhs);
+
+    problem_A = this->A;
+    problem_A *= this->tau;
 
+    // Treat as Delta u_n for now
     problem_iterate = this->u_old;
     if (this->u_old_old)
       problem_iterate -= *this->u_old_old;
@@ -67,17 +73,21 @@ class ImplicitEuler
             val - this->u_old[i][0]; // Time-dependent X direction
       } else if (this->dirichletNodes[i][1])
         problem_iterate[i][1] = 0; // Y direction described
+
+    // Make it a velocity
+    problem_iterate /= this->tau;
   }
 
   void virtual extractSolution(VectorType const &problem_iterate,
                                VectorType &solution) const {
     solution = this->u_old;
-    solution += problem_iterate;
+    solution.axpy(this->tau, problem_iterate);
   }
 
   void virtual extractDifference(VectorType const &problem_iterate,
                                  VectorType &difference) const {
     difference = problem_iterate;
+    difference *= this->tau;
   }
 };
 
@@ -94,13 +104,14 @@ class ImplicitTwoStep
 
   void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
                      MatrixType &problem_A) const {
-    problem_A = this->A;
-    problem_A /= 1.5;
     problem_rhs = this->ell;
-    problem_A.usmv(-2, this->u_old, problem_rhs);
-    problem_A.usmv(.5, *this->u_old_old, problem_rhs);
+    this->A.usmv(-4.0 / 3.0, this->u_old, problem_rhs);
+    this->A.usmv(+1.0 / 3.0, *this->u_old_old, problem_rhs);
+
+    problem_A = this->A;
+    problem_A *= 2.0 / 3.0 * this->tau;
 
-    // The finite difference makes a good start
+    // The finite difference makes a good start (treat it as such for now)
     problem_iterate = this->u_old;
     problem_iterate -= *this->u_old_old;
 
@@ -117,11 +128,14 @@ class ImplicitTwoStep
         // Y direction described
         problem_iterate[i][1] =
             -2 * this->u_old[i][1] + .5 * (*this->u_old_old)[i][1];
+    // Now treat it as a velocity
+    problem_iterate /= this->tau;
   }
 
   void virtual extractSolution(VectorType const &problem_iterate,
                                VectorType &solution) const {
     solution = problem_iterate;
+    solution *= this->tau;
     solution.axpy(2, this->u_old);
     solution.axpy(-.5, *this->u_old_old);
     solution *= 2.0 / 3.0;
@@ -129,6 +143,7 @@ class ImplicitTwoStep
     // 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, *this->u_old_old);
@@ -139,5 +154,6 @@ class ImplicitTwoStep
   void virtual extractDifference(VectorType const &problem_iterate,
                                  VectorType &difference) const {
     difference = problem_iterate;
+    difference *= this->tau;
   }
 };