[Noise] Clean up initial conditions

src/sand-wedge-data/boundaryconditions.py

@@ -6,6 +6,7 @@ class neumannCondition:
 
 class velocityDirichletCondition:
     def __call__(self, relativeTime):
+        # Assumed to vanish at time zero
         finalVelocity = -5e-5
         if relativeTime <= 0.1:
             return finalVelocity * ( 1.0 - math.cos(relativeTime * math.pi / 0.1) ) / 2.0

src/sand-wedge.cc

@@ -194,7 +194,6 @@ int main(int argc, char *argv[]) {
     myAssembler.assembleViscosity(body.getShearViscosityField(),
                                   body.getBulkViscosityField(), C);
     myAssembler.assembleMass(body.getDensityField(), M);
-    EnergyNorm<Matrix, Vector> const ANorm(A), MNorm(M);
 
     ScalarMatrix frictionalBoundaryMass;
     myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
@@ -213,21 +212,18 @@ int main(int argc, char *argv[]) {
                                       _relativeTime);
           _ell += gravityFunctional;
         };
-    Vector ell0(leafVertexCount);
-    computeExternalForces(0.0, ell0);
-
-    // {{{ Initial conditions
-    using LinearFactory =
-        SolverFactory<dims,
-                      BlockNonlinearTNNMGProblem<ConvexProblem<
-                          ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
-                      Grid>;
-    ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
 
+    // helper
     auto const solveLinearProblem = [&](
         Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
-        Vector const &_rhs, Vector &_x, EnergyNorm<Matrix, Vector> const &_norm,
+        Vector const &_rhs, Vector &_x,
         Dune::ParameterTree const &_localParset) {
+
+      using LinearFactory = SolverFactory<
+          dims, BlockNonlinearTNNMGProblem<ConvexProblem<
+                    ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
+          Grid>;
+      ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
       LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
                             refinements, *grid, _dirichletNodes);
 
@@ -237,9 +233,10 @@ int main(int argc, char *argv[]) {
 
       auto multigridStep = factory.getStep();
       multigridStep->setProblem(_x, problem);
+      EnergyNorm<Matrix, Vector> const norm(_matrix);
       LoopSolver<Vector> solver(
           multigridStep.get(), _localParset.get<size_t>("maximumIterations"),
-          _localParset.get<double>("tolerance"), &_norm,
+          _localParset.get<double>("tolerance"), &norm,
           _localParset.get<Solver::VerbosityMode>("verbosity"),
           false); // absolute error
 
@@ -247,52 +244,45 @@ int main(int argc, char *argv[]) {
       solver.solve();
     };
 
-    // Solve the stationary problem
+    // {{{ Initial conditions
+    Vector ell0(leafVertexCount);
+    computeExternalForces(0.0, ell0);
+
+    // Start from a situation of minimal stress, which is
+    // automatically attained in the case [v = 0 = a]. Assuming
+    // dPhi(v = 0) = 0, we thus only have to solve Au = ell0
     Vector u_initial(leafVertexCount);
-    u_initial = 0.0;
-    solveLinearProblem(dirichletNodes, A, ell0, u_initial, ANorm,
+    solveLinearProblem(dirichletNodes, A, ell0, u_initial,
                        parset.sub("u0.solver"));
 
-    ScalarVector alpha_initial(leafVertexCount);
-    alpha_initial = parset.get<double>("boundary.friction.initialAlpha");
     ScalarVector normalStress(leafVertexCount);
     myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
                                      body.getYoungModulus(),
                                      body.getPoissonRatio(), u_initial);
-    MyGlobalFrictionData<LocalVector> frictionInfo(
-        parset.sub("boundary.friction"), weakPatch);
-    auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
-        parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
-        frictionalBoundary, frictionInfo, normalStress);
-    myGlobalFriction->updateAlpha(alpha_initial);
 
+    ScalarVector alpha_initial(leafVertexCount);
+    alpha_initial = parset.get<double>("boundary.friction.initialAlpha");
+
+    // Start from zero velocity
     Vector v_initial(leafVertexCount);
     v_initial = 0.0;
-    {
-      double v_initial_const;
-      velocityDirichletFunction.evaluate(0.0, v_initial_const);
-      assert(std::abs(v_initial_const) < 1e-14);
-    }
 
+    // Compute the acceleration from Ma = ell0 - Au (without Dirichlet
+    // constraints), again assuming dPhi(v = 0) = 0
     Vector a_initial(leafVertexCount);
-    a_initial = 0.0;
-    {
-      // We solve Ma = ell0 - [Au + Cv + DPsi(v)]
-      Vector accelerationRHS(leafVertexCount);
-      {
-        accelerationRHS = 0.0;
-        Arithmetic::addProduct(accelerationRHS, A, u_initial);
-        Arithmetic::addProduct(accelerationRHS, C, v_initial);
-        // NOTE: We assume differentiability of Psi at 0 here!
-        myGlobalFriction->addGradient(v_initial, accelerationRHS);
-        accelerationRHS *= -1.0;
-        accelerationRHS += ell0;
-      }
-      solveLinearProblem(noNodes, M, accelerationRHS, a_initial, MNorm,
-                         parset.sub("a0.solver"));
-    }
+    Vector accelerationRHS = ell0;
+    Arithmetic::subtractProduct(accelerationRHS, A, u_initial);
+    solveLinearProblem(noNodes, M, accelerationRHS, a_initial,
+                       parset.sub("a0.solver"));
     // }}}
 
+    MyGlobalFrictionData<LocalVector> frictionInfo(
+        parset.sub("boundary.friction"), weakPatch);
+    auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
+        parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
+        frictionalBoundary, frictionInfo, normalStress);
+    myGlobalFriction->updateAlpha(alpha_initial);
+
     Vector vertexCoordinates(leafVertexCount);
     {
       Dune::MultipleCodimMultipleGeomTypeMapper<