Rename: h -> tau

src/assemblers.cc

@@ -57,7 +57,7 @@ Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType> const>
 assemble_nonlinearity(
     Dune::ParameterTree const &parset,
     Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double h) {
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double tau) {
   auto const size = nodalIntegrals.size();
 
   typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
@@ -84,7 +84,7 @@ assemble_nonlinearity(
 
       return Dune::make_shared<
           Dune::GlobalRuinaNonlinearity<MatrixType, VectorType> const>(
-          nodalIntegrals, a, mu, eta, normalStress, b, state, L, h);
+          nodalIntegrals, a, mu, eta, normalStress, b, state, L, tau);
     }
     case Config::Laursen:
       assert(false);

src/assemblers_tmpl.cc

@@ -34,4 +34,4 @@ template Dune::shared_ptr<
 assemble_nonlinearity<MatrixType, VectorType>(
     Dune::ParameterTree const &parset,
     Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double h);
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double tau);

src/compute_state.cc

@@ -17,11 +17,11 @@ class DecayingExponential {
   typedef Dune::FieldVector<double, 1> VectorType;
   typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
 
-  DecayingExponential(double h) : h(h) {}
+  DecayingExponential(double tau) : tau(tau) {}
 
   void directionalSubDiff(VectorType const &u, VectorType const &v,
                           Interval<double> &D) const {
-    D[0] = D[1] = v[0] * (-h * std::exp(-u[0]));
+    D[0] = D[1] = v[0] * (-tau * std::exp(-u[0]));
   }
 
   void directionalDomain(VectorType const &, VectorType const &,
@@ -31,14 +31,14 @@ class DecayingExponential {
   }
 
 private:
-  double const h;
+  double const tau;
 };
 
-double state_update_dieterich_bisection(double h, double uol,
+double state_update_dieterich_bisection(double tau, double uol,
                                         double old_state) {
   DecayingExponential::VectorType const start(0);
   DecayingExponential::VectorType const direction(1);
-  DecayingExponential const phi(h);
+  DecayingExponential const phi(tau);
   MyDirectionalConvexFunction<DecayingExponential> const J(
       1.0, old_state - uol, phi, start, direction);
 
@@ -47,18 +47,18 @@ double state_update_dieterich_bisection(double h, double uol,
   return bisection.minimize(J, 0.0, 0.0, bisectionsteps); // TODO
 }
 
-double state_update_dieterich(double h, double uol, double old_state) {
-  double const ret = state_update_dieterich_bisection(h, uol, old_state);
+double state_update_dieterich(double tau, double uol, double old_state) {
+  double const ret = state_update_dieterich_bisection(tau, uol, old_state);
   /* We have
 
-     ret - old_state + uol = h*exp(-ret)
+     ret - old_state + uol = tau*exp(-ret)
 
      or
 
-     log((ret - old_state + uol)/h) = -ret
+     log((ret - old_state + uol)/tau) = -ret
   */
-  assert(std::min(std::abs(ret - old_state + uol - h * std::exp(-ret)),
-                  std::abs(std::log((ret - old_state + uol) / h) + ret)) <
+  assert(std::min(std::abs(ret - old_state + uol - tau * std::exp(-ret)),
+                  std::abs(std::log((ret - old_state + uol) / tau) + ret)) <
          1e-12);
   return ret;
 }

src/compute_state_ruina.cc

@@ -4,10 +4,10 @@
 
 #include "compute_state_ruina.hh"
 
-double state_update_ruina(double h, double uol, double old_state) {
+double state_update_ruina(double tau, double uol, double old_state) {
   if (uol == 0)
     return old_state;
 
-  double const ret = old_state - uol * std::log(uol / h);
+  double const ret = old_state - uol * std::log(uol / tau);
   return ret / (1 + uol);
 }

src/one-body-sample.cc

@@ -222,7 +222,7 @@ int main(int argc, char *argv[]) {
     auto const eta = parset.get<double>("boundary.friction.eta");
     auto const mu = parset.get<double>("boundary.friction.mu");
     auto const timesteps = parset.get<size_t>("timesteps");
-    double const h = 1.0 / timesteps;
+    double const tau = 1.0 / timesteps;
 
     octave_writer << "# name: A" << std::endl << "# type: matrix" << std::endl
                   << "# rows: " << timesteps << std::endl << "# columns: 3"
@@ -242,7 +242,7 @@ int main(int argc, char *argv[]) {
       ++first_frictional_node;
 
     for (size_t run = 1; run <= timesteps; ++run) {
-      double const time = h * run;
+      double const time = tau * run;
       if (parset.get<bool>("printProgress")) {
         std::cout << '*';
         std::cout.flush();
@@ -263,7 +263,7 @@ int main(int argc, char *argv[]) {
              ++state_fpi) {
           auto myGlobalNonlinearity =
               assemble_nonlinearity<MatrixType, VectorType>(
-                  parset.sub("boundary.friction"), *nodalIntegrals, alpha, h);
+                  parset.sub("boundary.friction"), *nodalIntegrals, alpha, tau);
           MyConvexProblemType const myConvexProblem(stiffnessMatrix,
                                                     *myGlobalNonlinearity, rhs);
 
@@ -285,15 +285,16 @@ int main(int argc, char *argv[]) {
 
               // // the (logarithmic) steady state corresponding to the
               // // current velocity
-              // std::cout << std::log(L/unorm * h) << std::endl;
+              // std::cout << std::log(L/unorm * tau) << std::endl;
 
               switch (parset.get<Config::state_model>(
                   "boundary.friction.state.model")) {
                 case Config::Dieterich:
-                  alpha[i] = state_update_dieterich(h, unorm / L, alpha_old[i]);
+                  alpha[i] =
+                      state_update_dieterich(tau, unorm / L, alpha_old[i]);
                   break;
                 case Config::Ruina:
-                  alpha[i] = state_update_ruina(h, unorm / L, alpha_old[i]);
+                  alpha[i] = state_update_ruina(tau, unorm / L, alpha_old[i]);
                   break;
               }
             }
@@ -321,7 +322,7 @@ int main(int argc, char *argv[]) {
         // with the Ruina state evolution law.
         // Jumps at 120 and 360 timesteps; v1 = .6 * v2;
         if (parset.get<bool>("printVelocitySteppingComparison")) {
-          double const v = u_diff[first_frictional_node].two_norm() / h / L;
+          double const v = u_diff[first_frictional_node].two_norm() / tau / L;
 
           double const euler = alpha[first_frictional_node];
           double direct;
@@ -330,15 +331,15 @@ int main(int argc, char *argv[]) {
           } else if (run < 360) {
             double const v2 = v;
             double const v1 = 0.6 * v2;
-            direct =
-                std::log(1.0 / v2 *
-                         std::pow((v2 / v1), std::exp(-v2 * (run - 120) * h)));
+            direct = std::log(
+                1.0 / v2 *
+                std::pow((v2 / v1), std::exp(-v2 * (run - 120) * tau)));
           } else {
             double const v1 = v;
             double const v2 = v1 / 0.6;
-            direct =
-                std::log(1.0 / v1 *
-                         std::pow((v1 / v2), std::exp(-v1 * (run - 360) * h)));
+            direct = std::log(
+                1.0 / v1 *
+                std::pow((v1 / v2), std::exp(-v1 * (run - 360) * tau)));
           }
           velocity_stepping_writer << euler << " " << direct << std::endl;
         }