diff --git a/dune/gfe/mixedlocalgfeadolcstiffness.hh b/dune/gfe/mixedlocalgfeadolcstiffness.hh
index ddf89f9637cfe2bb06a2f1b28dd8f9cbadff10b8..75a34a0ae6ac4151be5c64ae9be782e43fb7e16b 100644
--- a/dune/gfe/mixedlocalgfeadolcstiffness.hh
+++ b/dune/gfe/mixedlocalgfeadolcstiffness.hh
@@ -86,12 +86,13 @@ energy(const typename Basis::LocalView& localView,
const std::vector<TargetSpace0>& localConfiguration0,
const std::vector<TargetSpace1>& localConfiguration1) const
{
+ int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
double pureEnergy;
std::vector<ATargetSpace0> localAConfiguration0(localConfiguration0.size());
std::vector<ATargetSpace1> localAConfiguration1(localConfiguration1.size());
- trace_on(1);
+ trace_on(rank);
adouble energy = 0;
@@ -122,10 +123,14 @@ energy(const typename Basis::LocalView& localView,
}
using namespace Dune::TypeTree::Indices;
- energy = localEnergy_->energy(localView,
+ try {
+ energy = localEnergy_->energy(localView,
localAConfiguration0,
localAConfiguration1);
-
+ } catch (Dune::Exception &e) {
+ trace_off();
+ throw e;
+ }
energy >>= pureEnergy;
trace_off();
@@ -197,6 +202,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
std::vector<typename TargetSpace0::TangentVector>& localGradient0,
std::vector<typename TargetSpace1::TangentVector>& localGradient1)
{
+ int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localConfiguration0, localConfiguration1);
@@ -222,7 +228,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
// Compute gradient
std::vector<double> g(nDoubles);
- gradient(1,nDoubles,xp.data(),g.data()); // gradient evaluation
+ gradient(rank,nDoubles,xp.data(),g.data()); // gradient evaluation
// Copy into Dune type
std::vector<typename TargetSpace0::EmbeddedTangentVector> localEmbeddedGradient0(localConfiguration0.size());
@@ -338,7 +344,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
for (int i=0; i<embeddedBlocksize1; i++)
tangent[nDofs0*embeddedBlocksize0 + (j/blocksize1)*embeddedBlocksize1+i][nDofs0*blocksize0 + j] = orthonormalFrame1[j/blocksize1][j%blocksize1][i];
- hess_mat(1,nDoubles,nDirections,xp.data(),tangent,rawHessian);
+ hess_mat(rank,nDoubles,nDirections,xp.data(),tangent,rawHessian);
// Copy Hessian into Dune data type
size_t offset0 = nDofs0*embeddedBlocksize0;
diff --git a/dune/gfe/mixedriemanniantrsolver.cc b/dune/gfe/mixedriemanniantrsolver.cc
index c09c7bf0e2fef75761b54e8d6be06391031f7640..2cd74a7ee6a9956e54588af84fdc364137a4cad5 100644
--- a/dune/gfe/mixedriemanniantrsolver.cc
+++ b/dune/gfe/mixedriemanniantrsolver.cc
@@ -385,6 +385,7 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,
corr_global[_0].resize(rhs_global[_0].size());
corr_global[_1].resize(rhs_global[_1].size());
corr_global = 0;
+ bool solvedByInnerSolver = true;
if (rank==0)
{
@@ -405,35 +406,46 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,
mmgStep1->preprocess();
+ ///////////////////////////////
+ // Solve !
+ ///////////////////////////////
std::cout << "Solve quadratic problem..." << std::endl;
double oldEnergy = 0;
Dune::Timer solutionTimer;
- for (int ii=0; ii<innerIterations_; ii++)
- {
- residual[_0] = rhs_global[_0];
- stiffnessMatrix[_0][_1].mmv(corr_global[_1], residual[_0]);
- mmgStep0->setRhs(residual[_0]);
- mmgStep0->iterate();
-
- residual[_1] = rhs_global[_1];
- stiffnessMatrix[_1][_0].mmv(corr_global[_0], residual[_1]);
- mmgStep1->setRhs(residual[_1]);
- mmgStep1->iterate();
-
- // Compute energy
- CorrectionType tmp(corr_global);
- stiffnessMatrix.mv(corr_global,tmp);
- double energy = 0.5 * (tmp*corr_global) - (rhs_global*corr_global);
-
- if (energy > oldEnergy)
- //DUNE_THROW(Dune::Exception, "energy increase!");
- std::cout << "Warning: energy increase!" << std::endl;
-
- oldEnergy = energy;
+ int ii = 0;
+ try {
+ for (; ii<innerIterations_; ii++) {
+ residual[_0] = rhs_global[_0];
+ stiffnessMatrix[_0][_1].mmv(corr_global[_1], residual[_0]);
+ mmgStep0->setRhs(residual[_0]);
+ mmgStep0->iterate();
+
+ residual[_1] = rhs_global[_1];
+ stiffnessMatrix[_1][_0].mmv(corr_global[_0], residual[_1]);
+ mmgStep1->setRhs(residual[_1]);
+ mmgStep1->iterate();
+
+ // Compute energy
+ CorrectionType tmp(corr_global);
+ stiffnessMatrix.mv(corr_global,tmp);
+ double energy = 0.5 * (tmp*corr_global) - (rhs_global*corr_global);
+
+ if (energy > oldEnergy)
+ //DUNE_THROW(Dune::Exception, "energy increase!");
+ std::cout << "Warning: energy increase!" << std::endl;
+
+ oldEnergy = energy;
+ if (corr_global.infinity_norm() < innerTolerance_)
+ break;
+ }
+ } catch (Dune::Exception &e) {
+ std::cerr << "Error while solving: " << e << std::endl;
+ solvedByInnerSolver = false;
+ corr_global = 0;
}
- std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
+ std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds and " << ii << " steps." << std::endl;
//std::cout << "Correction: " << std::endl << corr_global << std::endl;
@@ -445,73 +457,88 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,
//corr = vectorComm.scatter(corr_global);
corr = corr_global;
-
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
-
- if (corr_global.infinity_norm() < tolerance_) {
- if (verbosity_ == NumProc::FULL and rank==0)
- std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
-
- if (verbosity_ != NumProc::QUIET and rank==0)
- std::cout << i+1 << " trust-region steps were taken." << std::endl;
- break;
- }
-
- // ////////////////////////////////////////////////////
- // Check whether trust-region step can be accepted
- // ////////////////////////////////////////////////////
-
+ double energy = 0;
+ double modelDecrease = 0;
SolutionType newIterate = x_;
- for (size_t j=0; j<newIterate[_0].size(); j++)
- newIterate[_0][j] = TargetSpace0::exp(newIterate[_0][j], corr[_0][j]);
-
- for (size_t j=0; j<newIterate[_1].size(); j++)
- newIterate[_1][j] = TargetSpace1::exp(newIterate[_1][j], corr[_1][j]);
-
- double energy = assembler_->computeEnergy(newIterate[_0],newIterate[_1]);
- energy = mpiHelper.getCollectiveCommunication().sum(energy);
-
- // compute the model decrease
- // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
- // Note that rhs = -g
- CorrectionType tmp(corr);
- hessianMatrix_->mv(corr,tmp);
- double modelDecrease = rhs*corr - 0.5 * (corr*tmp);
- modelDecrease = mpiHelper.getCollectiveCommunication().sum(modelDecrease);
-
- double relativeModelDecrease = modelDecrease / std::fabs(energy);
-
- if (verbosity_ == NumProc::FULL and rank==0) {
- std::cout << "Absolute model decrease: " << modelDecrease
- << ", functional decrease: " << oldEnergy - energy << std::endl;
- std::cout << "Relative model decrease: " << relativeModelDecrease
- << ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
- }
+ if (solvedByInnerSolver) {
+
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
+
+ if (corr_global.infinity_norm() < tolerance_) {
+ if (verbosity_ == NumProc::FULL and rank==0)
+ std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
+
+ if (verbosity_ != NumProc::QUIET and rank==0)
+ std::cout << i+1 << " trust-region steps were taken." << std::endl;
+ break;
+ }
+
+ // ////////////////////////////////////////////////////
+ // Check whether trust-region step can be accepted
+ // ////////////////////////////////////////////////////
+
+ for (size_t j=0; j<newIterate[_0].size(); j++)
+ newIterate[_0][j] = TargetSpace0::exp(newIterate[_0][j], corr[_0][j]);
+
+ for (size_t j=0; j<newIterate[_1].size(); j++)
+ newIterate[_1][j] = TargetSpace1::exp(newIterate[_1][j], corr[_1][j]);
+
+ try {
+ energy = assembler_->computeEnergy(newIterate[_0],newIterate[_1]);
+ } catch (Dune::Exception &e) {
+ std::cerr << "Error while computing the energy of the new Iterate: " << e << std::endl;
+ std::cerr << "Redoing trust region step with smaller radius..." << std::endl;
+ newIterate = x_;
+ solvedByInnerSolver = false;
+ energy = oldEnergy;
+ }
+
+ if (solvedByInnerSolver) {
+ energy = mpiHelper.getCollectiveCommunication().sum(energy);
+
+ // compute the model decrease
+ // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
+ // Note that rhs = -g
+ CorrectionType tmp(corr);
+ hessianMatrix_->mv(corr,tmp);
+ modelDecrease = rhs*corr - 0.5 * (corr*tmp);
+ modelDecrease = mpiHelper.getCollectiveCommunication().sum(modelDecrease);
+
+ double relativeModelDecrease = modelDecrease / std::fabs(energy);
+
+ if (verbosity_ == NumProc::FULL and rank==0) {
+ std::cout << "Absolute model decrease: " << modelDecrease
+ << ", functional decrease: " << oldEnergy - energy << std::endl;
+ std::cout << "Relative model decrease: " << relativeModelDecrease
+ << ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
+ }
- assert(modelDecrease >= 0);
+ assert(modelDecrease >= 0);
- if (energy >= oldEnergy and rank==0) {
- if (this->verbosity_ == NumProc::FULL)
- printf("Richtung ist keine Abstiegsrichtung!\n");
- }
+ if (energy >= oldEnergy and rank==0) {
+ if (this->verbosity_ == NumProc::FULL)
+ printf("Richtung ist keine Abstiegsrichtung!\n");
+ }
- if (energy >= oldEnergy &&
- (std::abs((oldEnergy-energy)/energy) < 1e-9 || relativeModelDecrease < 1e-9)) {
- if (this->verbosity_ == NumProc::FULL and rank==0)
- std::cout << "Suspecting rounding problems" << std::endl;
+ if (energy >= oldEnergy &&
+ (std::abs((oldEnergy-energy)/energy) < 1e-9 || relativeModelDecrease < 1e-9)) {
+ if (this->verbosity_ == NumProc::FULL and rank==0)
+ std::cout << "Suspecting rounding problems" << std::endl;
- if (this->verbosity_ != NumProc::QUIET and rank==0)
- std::cout << i+1 << " trust-region steps were taken." << std::endl;
+ if (this->verbosity_ != NumProc::QUIET and rank==0)
+ std::cout << i+1 << " trust-region steps were taken." << std::endl;
- x_ = newIterate;
- break;
- }
+ x_ = newIterate;
+ break;
+ }
+ }
+ }
// //////////////////////////////////////////////
// Check for acceptance of the step
// //////////////////////////////////////////////
- if ( (oldEnergy-energy) / modelDecrease > 0.9) {
+ if ( solvedByInnerSolver && (oldEnergy-energy) / modelDecrease > 0.9) {
// very successful iteration
x_ = newIterate;
@@ -523,7 +550,7 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,
recomputeGradientHessian = true;
- } else if ( (oldEnergy-energy) / modelDecrease > 0.01
+ } else if (solvedByInnerSolver && (oldEnergy-energy) / modelDecrease > 0.01
|| std::abs(oldEnergy-energy) < 1e-12) {
// successful iteration
x_ = newIterate;
diff --git a/dune/gfe/realtuple.hh b/dune/gfe/realtuple.hh
index 5d47b0d178b959b9b8465589e8b66b5a66042ff0..c1f52db2e40acabdc6345ae94e88ed3beddf0616 100644
--- a/dune/gfe/realtuple.hh
+++ b/dune/gfe/realtuple.hh
@@ -66,6 +66,17 @@ public:
return *this;
}
+ RealTuple& operator-=(const Dune::FieldVector<T,N>& other) {
+ data_ -= other;
+ return *this;
+ }
+
+ template <class T2>
+ RealTuple& operator-=(const RealTuple<T2,N>& other) {
+ data_ -= other.data_;
+ return *this;
+ }
+
/** \brief Assigment from RealTuple with different type -- used for automatic differentiation with ADOL-C */
template <class T2>
RealTuple& operator <<= (const RealTuple<T2,N>& other) {
@@ -74,6 +85,11 @@ public:
return *this;
}
+ /** \brief Const random-access operator*/
+ T operator[] (const size_t indexVariable ) const {
+ return data_[indexVariable];
+ }
+
/** \brief Rebind the RealTuple to another coordinate type */
template<class U>
struct rebind
diff --git a/dune/gfe/riemannianpnsolver.cc b/dune/gfe/riemannianpnsolver.cc
index 0324629b6f6603a7d0c6d56932b0f645e80dcd55..4937720c4adabd0abb1ce46a66eed4ab1a2605be 100644
--- a/dune/gfe/riemannianpnsolver.cc
+++ b/dune/gfe/riemannianpnsolver.cc
@@ -22,10 +22,10 @@
#include <dune/gfe/parallel/vectorcommunicator.hh>
#endif
-template <class Basis, class TargetSpace>
-void RiemannianProximalNewtonSolver<Basis,TargetSpace>::
+template <class Basis, class TargetSpace, class Assembler>
+void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::
setup(const GridType& grid,
- const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
+ const Assembler* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -150,8 +150,8 @@ setup(const GridType& grid,
}
-template <class Basis, class TargetSpace>
-void RiemannianProximalNewtonSolver<Basis,TargetSpace>::solve()
+template <class Basis, class TargetSpace, class Assembler>
+void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
{
int rank = grid_->comm().rank();
diff --git a/dune/gfe/riemannianpnsolver.hh b/dune/gfe/riemannianpnsolver.hh
index a1c1aca5605b2e1c6fd8a579f562295ecfc39f9a..1024a72a7aaa1fe36268c2a84a02f98af279c8b4 100644
--- a/dune/gfe/riemannianpnsolver.hh
+++ b/dune/gfe/riemannianpnsolver.hh
@@ -17,7 +17,6 @@
#include <dune/gfe/periodic1dpq1nodalbasis.hh>
-#include "geodesicfeassembler.hh"
#include "riemanniantrsolver.hh"
#include <dune/grid/utility/globalindexset.hh>
#include <dune/gfe/parallel/globalmapper.hh>
@@ -26,7 +25,7 @@
#include <dune/gfe/parallel/p2mapper.hh>
/** \brief Riemannian proximal-newton solver for geodesic finite-element problems */
-template <class Basis, class TargetSpace>
+template <class Basis, class TargetSpace, class Assembler = GeodesicFEAssembler<Basis,TargetSpace>>
class RiemannianProximalNewtonSolver
: public IterativeSolver<std::vector<TargetSpace>,
Dune::BitSetVector<TargetSpace::TangentVector::dimension> >
@@ -70,7 +69,7 @@ public:
/** \brief Set up the solver using a choldmod solver as the inner solver */
void setup(const GridType& grid,
- const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
+ const Assembler* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -121,7 +120,7 @@ protected:
std::unique_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */
- const GeodesicFEAssembler<Basis, TargetSpace>* assembler_;
+ const Assembler* assembler_;
/** \brief The solver for the quadratic inner problems */
std::shared_ptr<typename Dune::Solvers::CholmodSolver<MatrixType,CorrectionType>> innerSolver_;
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index b8604dfc6c94f033342cbc61addc33454fd577fa..b45d502d01d5cf19495d1430d9e31c572112b456 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -28,10 +28,10 @@
#include <dune/gfe/parallel/vectorcommunicator.hh>
#endif
-template <class Basis, class TargetSpace>
-void RiemannianTrustRegionSolver<Basis,TargetSpace>::
+template <class Basis, class TargetSpace, class Assembler>
+void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::
setup(const GridType& grid,
- const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
+ const Assembler* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -306,8 +306,8 @@ setup(const GridType& grid,
}
-template <class Basis, class TargetSpace>
-void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
+template <class Basis, class TargetSpace, class Assembler>
+void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
{
int rank = grid_->comm().rank();
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index 70c69ffa7a518e134910d8a9a23c439205e7bec1..3184c86f91dfbcbbd56177908ae5366f53a4fb03 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -69,7 +69,7 @@ struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,3> >
};
/** \brief Riemannian trust-region solver for geodesic finite-element problems */
-template <class Basis, class TargetSpace>
+template <class Basis, class TargetSpace, class Assembler = GeodesicFEAssembler<Basis,TargetSpace>>
class RiemannianTrustRegionSolver
: public IterativeSolver<std::vector<TargetSpace>,
Dune::BitSetVector<TargetSpace::TangentVector::dimension> >
@@ -115,7 +115,7 @@ public:
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
- const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
+ const Assembler* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -188,7 +188,7 @@ protected:
std::unique_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */
- const GeodesicFEAssembler<Basis, TargetSpace>* assembler_;
+ const Assembler* assembler_;
/** \brief The solver for the quadratic inner problems */
std::shared_ptr<Solver> innerSolver_;