[Algorit] Carsten: Use eigenvalue approximation

dune/tectonic/ellipticenergy.hh

-#ifndef DUNE_TECTONIC_ELLIPTICENERGY_HH
-#define DUNE_TECTONIC_ELLIPTICENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/stdstreams.hh>
-
-#include <dune/fufem/arithmetic.hh>
-#include <dune/solvers/common/interval.hh>
-
-#include "localfriction.hh"
-
-template <size_t dim> class EllipticEnergy {
-public:
-  using LocalVector = Dune::FieldVector<double, dim>;
-  using LocalMatrix = Dune::FieldMatrix<double, dim, dim>;
-
-  using Nonlinearity = LocalFriction<dim>;
-
-  EllipticEnergy(LocalMatrix const &A, LocalVector const &b,
-                 std::shared_ptr<Nonlinearity const> phi,
-                 typename Dune::BitSetVector<dim>::const_reference ignore)
-      : A(A), b(b), phi(phi), ignore(ignore) {}
-
-  double operator()(LocalVector const &v) const {
-    return computeEnergy(A, v, b) + (*phi)(v);
-  }
-
-  LocalMatrix const &A;
-  LocalVector const &b;
-  std::shared_ptr<Nonlinearity const> const phi;
-  typename Dune::BitSetVector<dim>::const_reference const ignore;
-
-  void descentDirection(LocalVector const &x, LocalVector &y) const {
-    if (x.two_norm() >= 0.0) {
-      A.mv(x, y);
-      y -= b;
-      phi->addGradient(x, y);
-
-      for (size_t i = 0; i < dim; ++i)
-        if (ignore[i])
-          y[i] = 0;
-      y *= -1;
-    } else {
-      A.mv(x, y);
-      y -= b;
-
-      for (size_t i = 0; i < dim; ++i)
-        if (ignore[i])
-          y[i] = 0;
-      y *= -1;
-
-      // The contribution of phi to the directional derivative is the
-      // same for all directions here! Choose the one that is best of
-      // the smooth part and check if we have overall decline
-      Dune::Solvers::Interval<double> D;
-      phi->directionalSubDiff(x, y, D);
-      if (D[1] - y.two_norm2() >= 0.0)
-        y = 0.0;
-    }
-  }
-};
-#endif

dune/tectonic/minimisation.hh

@@ -11,57 +11,31 @@
 
 #include "mydirectionalconvexfunction.hh"
 
+// Warning: this exploits the property v*x = 0
 template <class Functional>
 double lineSearch(Functional const &J,
                   typename Functional::LocalVector const &x,
                   typename Functional::LocalVector const &v,
                   Bisection const &bisection) {
   MyDirectionalConvexFunction<typename Functional::Nonlinearity> const JRest(
-      computeDirectionalA(J.A, v), computeDirectionalb(J.A, J.b, x, v), *J.phi,
-      x, v);
-
+      J.alpha * v.two_norm2(), J.b * v, *J.phi, x, v);
   int count;
   return bisection.minimize(JRest, 0.0, 0.0, count);
 }
 
+/** Minimise a quadratic problem, for which both the quadratic and the
+    nonlinear term have gradients which point in the direction of
+    their arguments */
 template <class Functional>
 void minimise(Functional const &J, typename Functional::LocalVector &x,
-              size_t steps, Bisection const &bisection) {
-  using LocalVector = typename Functional::LocalVector;
-  auto const diff = [](LocalVector const &a, LocalVector const &b) {
-    LocalVector tmp = a;
-    tmp -= b;
-    return tmp.two_norm();
-  };
-
-  LocalVector x_initial = x;
-  LocalVector x_o = x;
-
-  for (size_t step = 0; step < steps; ++step) {
-    LocalVector v;
-    J.descentDirection(x, v);
-
-    double const vnorm = v.two_norm();
-    if (vnorm <= 0.0)
-      return;
-    v /= vnorm;
-
-    double const alpha = lineSearch(J, x, v, bisection);
-    Arithmetic::addProduct(x, alpha, v);
-
-    if (alpha < 1e-14) // TODO
-      break;
-
-    double const correction = diff(x, x_o);
-    double const overallCorrection = diff(x, x_initial);
-    if (overallCorrection <= 0.0)
-      return;
-
-    double const correctionQuotient = correction / overallCorrection;
-    if (correctionQuotient < 0.1) // enough descent; TODO
-      break;
-
-    x_o = x;
-  }
+              Bisection const &bisection) {
+  auto v = J.b;
+  double const vnorm = v.two_norm();
+  if (vnorm <= 0.0)
+    return;
+  v /= vnorm;
+
+  double const alpha = lineSearch(J, x, v, bisection);
+  Arithmetic::addProduct(x, alpha, v);
 }
 #endif

dune/tectonic/myblockproblem.hh

@@ -6,6 +6,7 @@
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/nullptr.hh>
 #include <dune/common/parametertree.hh>
+#include <dune/common/fmatrixev.hh>
 
 #include <dune/fufem/arithmetic.hh>
 #include <dune/solvers/common/interval.hh>
@@ -13,10 +14,10 @@
 #include <dune/tnnmg/problem-classes/bisection.hh>
 #include <dune/tnnmg/problem-classes/blocknonlineargsproblem.hh>
 
-#include "ellipticenergy.hh"
 #include "globalnonlinearity.hh"
 #include "minimisation.hh"
 #include "mydirectionalconvexfunction.hh"
+#include "quadraticenergy.hh"
 
 /** \brief Base class for problems where each block can be solved with a
  * modified gradient method */
@@ -56,8 +57,16 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
 
   MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
       : BNGSP(parset, problem),
+        maxEigenvalues_(problem.f.size()),
         parset_(parset),
-        localBisection(0.0, 1.0, 1e-12, false) {}
+        localBisection(0.0, 1.0, 1e-12, false) {
+    for (size_t i = 0; i < problem.f.size(); ++i) {
+      LocalVectorType eigenvalues;
+      Dune::FMatrixHelp::eigenValues(problem.A[i][i], eigenvalues);
+      maxEigenvalues_[i] =
+          *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
+    }
+  }
 
   std::string getOutput(bool header = false) const {
     if (header) {
@@ -188,10 +197,12 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
 
   /** \brief Constructs and returns an iterate object */
   IterateObject getIterateObject() {
-    return IterateObject(parset_, localBisection, problem_);
+    return IterateObject(parset_, localBisection, problem_, maxEigenvalues_);
   }
 
 private:
+  std::vector<double> maxEigenvalues_;
+
   Dune::ParameterTree const &parset_;
 
   // problem data
@@ -213,10 +224,11 @@ class MyBlockProblem<ConvexProblem>::IterateObject {
    * \param problem The problem including quadratic part and nonlinear part
    */
   IterateObject(Dune::ParameterTree const &parset, Bisection const &bisection,
-                ConvexProblem const &problem)
+                ConvexProblem const &problem,
+                std::vector<double> const &maxEigenvalues)
       : problem(problem),
-        bisection_(bisection),
-        localsteps(parset.get<size_t>("localsolver.steps")) {}
+        maxEigenvalues_(maxEigenvalues),
+        bisection_(bisection) {}
 
 public:
   /** \brief Set the current iterate */
@@ -240,22 +252,24 @@ class MyBlockProblem<ConvexProblem>::IterateObject {
       LocalVectorType &ui, size_t m,
       typename Dune::BitSetVector<block_size>::const_reference ignore) {
     {
-      LocalMatrixType const *localA = nullptr;
-      LocalVectorType localb(problem.f[m]);
-
+      LocalVectorType localb = problem.f[m];
       auto const end = std::end(problem.A[m]);
       for (auto it = std::begin(problem.A[m]); it != end; ++it) {
         size_t const j = it.index();
-        if (j == m)
-          localA = &(*it); // localA = A[m][m]
-        else
-          Arithmetic::subtractProduct(localb, *it, u[j]);
+        Arithmetic::subtractProduct(localb, *it, u[j]); // also the diagonal!
       }
-      assert(localA != nullptr);
+      Arithmetic::addProduct(localb, maxEigenvalues_[m], u[m]);
 
-      auto const phi = problem.phi.restriction(m);
-      EllipticEnergy<block_size> localJ(*localA, localb, phi, ignore);
-      minimise(localJ, ui, localsteps, bisection_);
+      // We minimise over an affine subspace
+      for (size_t j = 0; j < block_size; ++j)
+        if (ignore[j])
+          localb[j] = 0;
+        else
+          ui[j] = 0;
+
+      QuadraticEnergy<typename ConvexProblem::NonlinearityType::Friction>
+      localJ(maxEigenvalues_[m], localb, problem.phi.restriction(m));
+      minimise(localJ, ui, bisection_);
     }
   }
 
@@ -263,12 +277,12 @@ class MyBlockProblem<ConvexProblem>::IterateObject {
   // problem data
   ConvexProblem const &problem;
 
+  std::vector<double> maxEigenvalues_;
+
   Bisection const bisection_;
 
   // state data for smoothing procedure used by:
   // setIterate, updateIterate, solveLocalProblem
   VectorType u;
-
-  size_t const localsteps;
 };
 #endif

dune/tectonic/quadraticenergy.hh

+#ifndef DUNE_TECTONIC_QUADRATICENERGY_HH
+#define DUNE_TECTONIC_QUADRATICENERGY_HH
+
+#include <memory>
+
+template <class NonlinearityTEMPLATE> class QuadraticEnergy {
+public:
+  using Nonlinearity = NonlinearityTEMPLATE;
+  using LocalVector = typename Nonlinearity::VectorType;
+
+  QuadraticEnergy(double alpha, LocalVector const &b,
+                  std::shared_ptr<Nonlinearity const> phi)
+      : alpha(alpha), b(b), phi(phi) {}
+
+  double const alpha;
+  LocalVector const &b;
+  std::shared_ptr<Nonlinearity const> const phi;
+};
+#endif

src/sand-wedge-data/parset.cfg

@@ -73,6 +73,3 @@ post               = 3
 pre   = 1
 multi = 5 # number of multigrid steps
 post  = 0
-
-[localsolver]
-steps = 100