Merge branch 'enhancements' into 'master'

Enhancements

See merge request !52

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;

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;

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

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();
 

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_;

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();
 

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_;