From 6934ce2f8ad1ee25de4dae0338bb8948f21c755e Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 23 Oct 2014 01:15:13 +0000
Subject: [PATCH] A Euclidean trust region solver
[[Imported from SVN: r9927]]
---
dune/gfe/Makefile.am | 2 +
dune/gfe/trustregionsolver.cc | 402 ++++++++++++++++++++++++++++++++++
dune/gfe/trustregionsolver.hh | 137 ++++++++++++
3 files changed, 541 insertions(+)
create mode 100644 dune/gfe/trustregionsolver.cc
create mode 100644 dune/gfe/trustregionsolver.hh
diff --git a/dune/gfe/Makefile.am b/dune/gfe/Makefile.am
index cb75f032..58ee4494 100644
--- a/dune/gfe/Makefile.am
+++ b/dune/gfe/Makefile.am
@@ -41,6 +41,8 @@ srcinclude_HEADERS = adolcnamespaceinjections.hh \
targetspacertrsolver.hh \
tensorssd.hh \
tensor3.hh \
+ trustregionsolver.hh \
+ trustregionsolver.cc \
unitvector.hh \
vtkfile.hh
diff --git a/dune/gfe/trustregionsolver.cc b/dune/gfe/trustregionsolver.cc
new file mode 100644
index 00000000..b5894182
--- /dev/null
+++ b/dune/gfe/trustregionsolver.cc
@@ -0,0 +1,402 @@
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/timer.hh>
+
+#include <dune/istl/io.hh>
+
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/assemblers/operatorassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
+
+// Using a monotone multigrid as the inner solver
+#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
+#include <dune/solvers/iterationsteps/mmgstep.hh>
+#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
+#if defined THIRD_ORDER || defined SECOND_ORDER
+#include <dune/gfe/pktop1mgtransfer.hh>
+#endif
+#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include "maxnormtrustregion.hh"
+
+#include <dune/solvers/norms/twonorm.hh>
+#include <dune/solvers/norms/h1seminorm.hh>
+
+template <class GridType, class VectorType>
+void TrustRegionSolver<GridType,VectorType>::
+setup(const GridType& grid,
+ const FEAssembler<BasisType, VectorType>* assembler,
+ const SolutionType& x,
+ const Dune::BitSetVector<blocksize>& dirichletNodes,
+ double tolerance,
+ int maxTrustRegionSteps,
+ double initialTrustRegionRadius,
+ int multigridIterations,
+ double mgTolerance,
+ int mu,
+ int nu1,
+ int nu2,
+ int baseIterations,
+ double baseTolerance)
+{
+ grid_ = &grid;
+ assembler_ = assembler;
+ x_ = x;
+ this->tolerance_ = tolerance;
+ maxTrustRegionSteps_ = maxTrustRegionSteps;
+ initialTrustRegionRadius_ = initialTrustRegionRadius;
+ innerIterations_ = multigridIterations;
+ innerTolerance_ = mgTolerance;
+ ignoreNodes_ = &dirichletNodes;
+
+ int numLevels = grid_->maxLevel()+1;
+
+ // ////////////////////////////////
+ // Create a multigrid solver
+ // ////////////////////////////////
+
+#if 0//def HAVE_IPOPT
+ // First create an IPOpt base solver
+ QuadraticIPOptSolver<MatrixType, CorrectionType>* baseSolver = new QuadraticIPOptSolver<MatrixType,CorrectionType>;
+ baseSolver->verbosity_ = NumProc::QUIET;
+ baseSolver->tolerance_ = baseTolerance;
+#else
+ // First create a Gauss-seidel base solver
+ TrustRegionGSStep<MatrixType, CorrectionType>* baseSolverStep = new TrustRegionGSStep<MatrixType, CorrectionType>;
+
+ // Hack: the two-norm may not scale all that well, but it is fast!
+ TwoNorm<CorrectionType>* baseNorm = new TwoNorm<CorrectionType>;
+
+ ::LoopSolver<CorrectionType>* baseSolver = new ::LoopSolver<CorrectionType>(baseSolverStep,
+ baseIterations,
+ baseTolerance,
+ baseNorm,
+ Solver::QUIET);
+#endif
+
+ // Make pre and postsmoothers
+ TrustRegionGSStep<MatrixType, CorrectionType>* presmoother = new TrustRegionGSStep<MatrixType, CorrectionType>;
+ TrustRegionGSStep<MatrixType, CorrectionType>* postsmoother = new TrustRegionGSStep<MatrixType, CorrectionType>;
+
+ MonotoneMGStep<MatrixType, CorrectionType>* mmgStep = new MonotoneMGStep<MatrixType, CorrectionType>;
+
+ mmgStep->setMGType(mu, nu1, nu2);
+ mmgStep->ignoreNodes_ = &dirichletNodes;
+ mmgStep->basesolver_ = baseSolver;
+ mmgStep->setSmoother(presmoother, postsmoother);
+ mmgStep->obstacleRestrictor_= new MandelObstacleRestrictor<CorrectionType>();
+ mmgStep->verbosity_ = Solver::QUIET;
+
+ // //////////////////////////////////////////////////////////////////////////////////////
+ // Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
+ // //////////////////////////////////////////////////////////////////////////////////////
+
+ BasisType basis(grid.leafGridView());
+ OperatorAssembler<BasisType,BasisType> operatorAssembler(basis, basis);
+
+ LaplaceAssembler<GridType, typename BasisType::LocalFiniteElement, typename BasisType::LocalFiniteElement> laplaceStiffness;
+ typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
+ ScalarMatrixType localA;
+
+ operatorAssembler.assemble(laplaceStiffness, localA);
+
+ if (h1SemiNorm_)
+ delete h1SemiNorm_;
+
+ ScalarMatrixType* A = new ScalarMatrixType(localA);
+
+ h1SemiNorm_ = new H1SemiNorm<CorrectionType>(*A);
+
+ innerSolver_ = std::shared_ptr<LoopSolver<CorrectionType> >(new ::LoopSolver<CorrectionType>(mmgStep,
+ innerIterations_,
+ innerTolerance_,
+ h1SemiNorm_,
+ Solver::REDUCED));
+
+ // //////////////////////////////////////////////////////////////////////////////////////
+ // Assemble a mass matrix to create a norm that's equivalent to the L2-norm
+ // This will be used to monitor the gradient
+ // //////////////////////////////////////////////////////////////////////////////////////
+
+ MassAssembler<GridType, typename BasisType::LocalFiniteElement, typename BasisType::LocalFiniteElement> massStiffness;
+ ScalarMatrixType localMassMatrix;
+
+ operatorAssembler.assemble(massStiffness, localMassMatrix);
+
+ ScalarMatrixType* massMatrix = new ScalarMatrixType(localMassMatrix);
+ l2Norm_ = std::make_shared<H1SemiNorm<CorrectionType> >(*massMatrix);
+
+ // ////////////////////////////////////////////////////////////
+ // Create Hessian matrix and its occupation structure
+ // ////////////////////////////////////////////////////////////
+
+ hessianMatrix_ = std::auto_ptr<MatrixType>(new MatrixType);
+ Dune::MatrixIndexSet indices(grid_->size(1), grid_->size(1));
+ assembler_->getNeighborsPerVertex(indices);
+ indices.exportIdx(*hessianMatrix_);
+
+ // ////////////////////////////////////
+ // Create the transfer operators
+ // ////////////////////////////////////
+
+ for (size_t k=0; k<mmgStep->mgTransfer_.size(); k++)
+ delete(mmgStep->mgTransfer_[k]);
+
+ mmgStep->mgTransfer_.resize(numLevels-1);
+
+#if defined THIRD_ORDER || defined SECOND_ORDER
+ if (numLevels>1) {
+ P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafGridView());
+
+ PKtoP1MGTransfer<CorrectionType>* topTransferOp = new PKtoP1MGTransfer<CorrectionType>;
+ topTransferOp->setup(basis,p1Basis);
+
+ mmgStep->mgTransfer_.back() = topTransferOp;
+
+ for (int i=0; i<mmgStep->mgTransfer_.size()-1; i++){
+ // Construct the local multigrid transfer matrix
+ TruncatedCompressedMGTransfer<CorrectionType>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType>;
+ newTransferOp->setup(*grid_,i+1,i+2);
+
+ mmgStep->mgTransfer_[i] = newTransferOp;
+ }
+
+ }
+
+#else
+ for (size_t i=0; i<mmgStep->mgTransfer_.size(); i++){
+
+ // Construct the local multigrid transfer matrix
+ TruncatedCompressedMGTransfer<CorrectionType>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType>;
+ newTransferOp->setup(*grid_,i,i+1);
+
+ mmgStep->mgTransfer_[i] = newTransferOp;;
+ }
+#endif
+
+ // //////////////////////////////////////////////////////////
+ // Create obstacles
+ // //////////////////////////////////////////////////////////
+
+ hasObstacle_.resize(basis.size(), true);
+ mmgStep->hasObstacle_ = &hasObstacle_;
+
+}
+
+
+template <class GridType, class VectorType>
+void TrustRegionSolver<GridType,VectorType>::solve()
+{
+ MonotoneMGStep<MatrixType,CorrectionType>* mgStep = NULL;
+
+ // if the inner solver is a monotone multigrid set up a max-norm trust-region
+ if (dynamic_cast<LoopSolver<CorrectionType>*>(innerSolver_.get())) {
+ mgStep = dynamic_cast<MonotoneMGStep<MatrixType,CorrectionType>*>(dynamic_cast<LoopSolver<CorrectionType>*>(innerSolver_.get())->iterationStep_);
+
+ }
+
+ BasisType basis(grid_->leafGridView());
+ MaxNormTrustRegion<blocksize> trustRegion(basis.size(), initialTrustRegionRadius_);
+
+ std::vector<BoxConstraint<field_type,blocksize> > trustRegionObstacles;
+
+ // /////////////////////////////////////////////////////
+ // Trust-Region Solver
+ // /////////////////////////////////////////////////////
+
+ double oldEnergy = assembler_->computeEnergy(x_);
+
+ bool recomputeGradientHessian = true;
+ CorrectionType rhs;
+ MatrixType stiffnessMatrix;
+
+ for (int i=0; i<maxTrustRegionSteps_; i++) {
+
+ Dune::Timer totalTimer;
+ if (this->verbosity_ == Solver::FULL) {
+ std::cout << "----------------------------------------------------" << std::endl;
+ std::cout << " Trust-Region Step Number: " << i
+ << ", radius: " << trustRegion.radius()
+ << ", energy: " << oldEnergy << std::endl;
+ std::cout << "----------------------------------------------------" << std::endl;
+ }
+
+ Dune::Timer gradientTimer;
+
+ if (recomputeGradientHessian) {
+
+ assembler_->assembleGradientAndHessian(x_,
+ rhs,
+ *hessianMatrix_,
+ i==0 // assemble occupation pattern only for the first call
+ );
+
+ rhs *= -1; // The right hand side is the _negative_ gradient
+
+ // Compute gradient norm to monitor convergence
+ CorrectionType gradient = rhs;
+ for (size_t j=0; j<gradient.size(); j++)
+ for (int k=0; k<gradient[j].size(); k++)
+ if ((*ignoreNodes_)[j][k])
+ gradient[j][k] = 0;
+
+ if (this->verbosity_ == Solver::FULL)
+ std::cout << "Gradient norm: " << l2Norm_->operator()(gradient) << std::endl;
+
+ if (this->verbosity_ == Solver::FULL)
+ std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
+
+ // Transfer matrix data
+ stiffnessMatrix = *hessianMatrix_;
+
+ recomputeGradientHessian = false;
+
+ }
+
+ CorrectionType corr(rhs.size());
+ corr = 0;
+
+ mgStep->setProblem(stiffnessMatrix, corr, rhs);
+
+ trustRegionObstacles = trustRegion.obstacles();
+ mgStep->obstacles_ = &trustRegionObstacles;
+
+ innerSolver_->preprocess();
+
+ ///////////////////////////////
+ // Solve !
+ ///////////////////////////////
+
+ std::cout << "Solve quadratic problem..." << std::endl;
+
+ Dune::Timer solutionTimer;
+ innerSolver_->solve();
+ std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
+
+ if (mgStep)
+ corr = mgStep->getSol();
+
+ //std::cout << "Correction: " << std::endl << corr_global << std::endl;
+
+ // Output correction for debugging
+ Dune::VTKWriter<typename GridType::LeafGridView> vtkWriter(grid_->leafGridView());
+
+ Dune::BlockVector<Dune::FieldVector<double,3> > displacement(x_.size());
+ for (size_t j=0; j<x_.size(); j++)
+ displacement[j] = x_[j] - identity_[j];
+
+ BasisType basis(grid_->leafGridView());
+ Dune::shared_ptr<VTKBasisGridFunction<BasisType,Dune::BlockVector<Dune::FieldVector<double,3> > > > vtkDisplacement
+ = Dune::make_shared<VTKBasisGridFunction<BasisType,Dune::BlockVector<Dune::FieldVector<double,3> > > >
+ (basis, displacement, "Displacement");
+
+ Dune::shared_ptr<VTKBasisGridFunction<BasisType,Dune::BlockVector<Dune::FieldVector<double,3> > > > vtkCorrection
+ = Dune::make_shared<VTKBasisGridFunction<BasisType,Dune::BlockVector<Dune::FieldVector<double,3> > > >
+ (basis, corr, "Correction");
+
+ Dune::shared_ptr<VTKBasisGridFunction<BasisType,Dune::BlockVector<Dune::FieldVector<double,3> > > > vtkGradient
+ = Dune::make_shared<VTKBasisGridFunction<BasisType,Dune::BlockVector<Dune::FieldVector<double,3> > > >
+ (basis, rhs, "Gradient");
+
+ vtkWriter.addVertexData(vtkDisplacement);
+ vtkWriter.addVertexData(vtkCorrection);
+ vtkWriter.addVertexData(vtkGradient);
+ vtkWriter.write("hencky_correction_" + std::to_string(i+1));
+
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
+
+ if (corr.infinity_norm() < this->tolerance_) {
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
+
+ if (this->verbosity_ != NumProc::QUIET)
+ std::cout << i+1 << " trust-region steps were taken." << std::endl;
+ break;
+ }
+
+ // ////////////////////////////////////////////////////
+ // Check whether trust-region step can be accepted
+ // ////////////////////////////////////////////////////
+
+ SolutionType newIterate = x_;
+ for (size_t j=0; j<newIterate.size(); j++)
+ newIterate[j] += corr[j];
+
+ double energy = assembler_->computeEnergy(newIterate);
+
+ // compute the model decrease
+ // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
+ // Note that rhs = -g
+ CorrectionType tmp(corr.size());
+ tmp = 0;
+ hessianMatrix_->umv(corr, tmp);
+ double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
+
+ double relativeModelDecrease = modelDecrease / std::fabs(energy);
+
+ if (this->verbosity_ == NumProc::FULL) {
+ 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);
+
+ if (energy >= oldEnergy) {
+ 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)
+ std::cout << "Suspecting rounding problems" << std::endl;
+
+ if (this->verbosity_ != NumProc::QUIET)
+ std::cout << i+1 << " trust-region steps were taken." << std::endl;
+
+ x_ = newIterate;
+ break;
+ }
+
+ // //////////////////////////////////////////////
+ // Check for acceptance of the step
+ // //////////////////////////////////////////////
+ if ( (oldEnergy-energy) / modelDecrease > 0.9) {
+ // very successful iteration
+
+ x_ = newIterate;
+ trustRegion.scale(2);
+
+ // current energy becomes 'oldEnergy' for the next iteration
+ oldEnergy = energy;
+
+ recomputeGradientHessian = true;
+
+ } else if ( (oldEnergy-energy) / modelDecrease > 0.01
+ || std::abs(oldEnergy-energy) < 1e-12) {
+ // successful iteration
+ x_ = newIterate;
+
+ // current energy becomes 'oldEnergy' for the next iteration
+ oldEnergy = energy;
+
+ recomputeGradientHessian = true;
+
+ } else {
+
+ // unsuccessful iteration
+
+ // Decrease the trust-region radius
+ trustRegion.scale(0.5);
+
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Unsuccessful iteration!" << std::endl;
+ }
+
+ std::cout << "iteration took " << totalTimer.elapsed() << " sec." << std::endl;
+ }
+
+}
diff --git a/dune/gfe/trustregionsolver.hh b/dune/gfe/trustregionsolver.hh
new file mode 100644
index 00000000..101d1ddf
--- /dev/null
+++ b/dune/gfe/trustregionsolver.hh
@@ -0,0 +1,137 @@
+#ifndef DUNE_GFE_TRUST_REGION_SOLVER_HH
+#define DUNE_GFE_TRUST_REGION_SOLVER_HH
+
+#include <vector>
+
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/solvers/common/boxconstraint.hh>
+#include <dune/solvers/norms/h1seminorm.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/functionspacebases/p2nodalbasis.hh>
+#include <dune/fufem/functionspacebases/p3nodalbasis.hh>
+
+#include <dune/gfe/feassembler.hh>
+
+/** \brief Trust-region solver */
+template <class GridType, class VectorType>
+class TrustRegionSolver
+ : public IterativeSolver<VectorType,
+ Dune::BitSetVector<VectorType::value_type::dimension> >
+{
+ const static int blocksize = VectorType::value_type::dimension;
+
+ const static int gridDim = GridType::dimension;
+
+ // Centralize the field type here
+ typedef double field_type;
+
+ // Some types that I need
+ typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize, blocksize> > MatrixType;
+ typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > CorrectionType;
+ typedef VectorType SolutionType;
+
+#ifdef THIRD_ORDER
+ typedef P3NodalBasis<typename GridType::LeafGridView,double> BasisType;
+#elif defined SECOND_ORDER
+ typedef P2NodalBasis<typename GridType::LeafGridView,double> BasisType;
+#else
+ typedef P1NodalBasis<typename GridType::LeafGridView,double> BasisType;
+#endif
+
+public:
+
+ TrustRegionSolver()
+ : IterativeSolver<VectorType, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
+ hessianMatrix_(std::auto_ptr<MatrixType>(NULL)), h1SemiNorm_(NULL)
+ {}
+
+ /** \brief Set up the solver using a monotone multigrid method as the inner solver */
+ void setup(const GridType& grid,
+ const FEAssembler<BasisType,VectorType>* assembler,
+ const SolutionType& x,
+ const Dune::BitSetVector<blocksize>& dirichletNodes,
+ double tolerance,
+ int maxTrustRegionSteps,
+ double initialTrustRegionRadius,
+ int multigridIterations,
+ double mgTolerance,
+ int mu,
+ int nu1,
+ int nu2,
+ int baseIterations,
+ double baseTolerance);
+
+ void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
+ {
+ ignoreNodes_ = &ignoreNodes;
+ Dune::shared_ptr<LoopSolver<CorrectionType> > loopSolver = std::dynamic_pointer_cast<LoopSolver<CorrectionType> >(innerSolver_);
+ assert(loopSolver);
+ loopSolver->iterationStep_->ignoreNodes_ = ignoreNodes_;
+ }
+
+ void solve();
+
+ void setInitialSolution(const SolutionType& x) DUNE_DEPRECATED {
+ x_ = x;
+ }
+
+ void setInitialIterate(const SolutionType& x) {
+ x_ = x;
+ }
+
+ SolutionType getSol() const {return x_;}
+
+protected:
+
+ /** \brief The grid */
+ const GridType* grid_;
+
+ /** \brief The solution vector */
+ SolutionType x_;
+
+ /** \brief The initial trust-region radius in the maximum-norm */
+ double initialTrustRegionRadius_;
+
+ /** \brief Maximum number of trust-region steps */
+ int maxTrustRegionSteps_;
+
+ /** \brief Maximum number of multigrid iterations */
+ int innerIterations_;
+
+ /** \brief Error tolerance of the multigrid QP solver */
+ double innerTolerance_;
+
+ /** \brief Hessian matrix */
+ std::auto_ptr<MatrixType> hessianMatrix_;
+
+ /** \brief The assembler for the material law */
+ const FEAssembler<BasisType, VectorType>* assembler_;
+
+ /** \brief The solver for the quadratic inner problems */
+ std::shared_ptr<Solver> innerSolver_;
+
+ /** \brief Contains 'true' everywhere -- the trust-region is bounded */
+ Dune::BitSetVector<blocksize> hasObstacle_;
+
+ /** \brief The Dirichlet nodes */
+ const Dune::BitSetVector<blocksize>* ignoreNodes_;
+
+ /** \brief The norm used to measure multigrid convergence */
+ H1SemiNorm<CorrectionType>* h1SemiNorm_;
+
+ /** \brief An L2-norm, really. The H1SemiNorm class is badly named */
+ std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_;
+public:
+ VectorType identity_;
+};
+
+#include "trustregionsolver.cc"
+
+#endif
--
GitLab