From 79e02a3c32cc1441584a4e5e6ae5e81c1757cd99 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sat, 31 May 2014 01:48:47 +0200
Subject: [PATCH] [Algorit] Carsten: Use eigenvalue approximation
---
dune/tectonic/ellipticenergy.hh | 63 --------------------------------
dune/tectonic/minimisation.hh | 54 +++++++--------------------
dune/tectonic/myblockproblem.hh | 52 ++++++++++++++++----------
dune/tectonic/quadraticenergy.hh | 19 ++++++++++
src/sand-wedge-data/parset.cfg | 3 --
5 files changed, 66 insertions(+), 125 deletions(-)
delete mode 100644 dune/tectonic/ellipticenergy.hh
create mode 100644 dune/tectonic/quadraticenergy.hh
diff --git a/dune/tectonic/ellipticenergy.hh b/dune/tectonic/ellipticenergy.hh
deleted file mode 100644
index e3405a1c..00000000
--- a/dune/tectonic/ellipticenergy.hh
+++ /dev/null
@@ -1,63 +0,0 @@
-#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
diff --git a/dune/tectonic/minimisation.hh b/dune/tectonic/minimisation.hh
index 3ae1464a..4d4b172b 100644
--- a/dune/tectonic/minimisation.hh
+++ b/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
diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index 9770d244..0f43db0b 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/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
diff --git a/dune/tectonic/quadraticenergy.hh b/dune/tectonic/quadraticenergy.hh
new file mode 100644
index 00000000..f724a1e8
--- /dev/null
+++ b/dune/tectonic/quadraticenergy.hh
@@ -0,0 +1,19 @@
+#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
diff --git a/src/sand-wedge-data/parset.cfg b/src/sand-wedge-data/parset.cfg
index 4287864a..f6ad4f34 100644
--- a/src/sand-wedge-data/parset.cfg
+++ b/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
--
GitLab