From 16374b9e2a4666d6344ed2bce4aeb87d1a2b075c Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Thu, 16 Apr 2020 15:22:54 +0200
Subject: [PATCH] Cleanup of the TrustRegionSolver class
Mainly issues related to dune-parmg.
---
dune/elasticity/common/trustregionsolver.cc | 143 +++++++-------------
dune/elasticity/common/trustregionsolver.hh | 6 +
2 files changed, 55 insertions(+), 94 deletions(-)
diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index 5d16bf1..0d1da51 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/dune/elasticity/common/trustregionsolver.cc
@@ -40,16 +40,8 @@ setup(const typename BasisType::GridView::Grid& grid,
int nu2,
int baseIterations,
double baseTolerance)
-{
-#if HAVE_DUNE_PARMG
- Dune::Timer setupTimer;
- mgSetup_ = std::make_unique<MGSetup>(grid);
- mgSetup_->overlap(true);
- mgSetup_->setupSmoother(1.0);
- if (grid.comm().rank()==0)
- std::cout << "Parallel multigrid setup took " << setupTimer.stop() << " seconds." << std::endl;
-#endif
+{
grid_ = &grid;
assembler_ = assembler;
x_ = x;
@@ -59,10 +51,30 @@ setup(const typename BasisType::GridView::Grid& grid,
innerIterations_ = multigridIterations;
innerTolerance_ = mgTolerance;
ignoreNodes_ = std::make_shared<Dune::BitSetVector<blocksize>>(dirichletNodes);
+ baseIterations_ = baseIterations;
+ baseTolerance_ = baseTolerance;
int numLevels = grid_->maxLevel()+1;
-#if !HAVE_DUNE_PARMG
+#if HAVE_DUNE_PARMG
+ Dune::Timer setupTimer;
+ mgSetup_ = std::make_unique<MGSetup>(grid);
+ mgSetup_->overlap(false);
+ mgSetup_->ignore(std::const_pointer_cast<Dune::BitSetVector<blocksize> >(ignoreNodes_));
+
+ BasisType feBasis(grid.levelGridView(grid.maxLevel()));
+ DuneFunctionsBasis<BasisType> basis(feBasis);
+
+ innerMultigridStep_ = std::make_unique<Dune::ParMG::Multigrid<VectorType> >();
+ innerMultigridStep_->mu_ = mu;
+ innerMultigridStep_->preSmootherSteps_ = nu1;
+ innerMultigridStep_->postSmootherSteps_ = nu2;
+
+ if (grid.comm().rank()==0)
+ std::cout << "Parallel multigrid setup took " << setupTimer.stop() << " seconds." << std::endl;
+#else
+ DuneFunctionsBasis<BasisType> basis(grid.leafGridView());
+
// ////////////////////////////////
// Create a multigrid solver
// ////////////////////////////////
@@ -87,11 +99,6 @@ setup(const typename BasisType::GridView::Grid& grid,
baseNorm,
Solver::QUIET);
#endif
-#endif
-
-#if HAVE_DUNE_PARMG
- mgSetup_->ignore(std::const_pointer_cast<Dune::BitSetVector<blocksize> >(ignoreNodes_));
-#else
// Make pre and postsmoothers
auto presmoother = std::make_shared< TrustRegionGSStep<MatrixType, CorrectionType> >();
auto postsmoother = std::make_shared< TrustRegionGSStep<MatrixType, CorrectionType> >();
@@ -108,13 +115,6 @@ setup(const typename BasisType::GridView::Grid& grid,
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
// //////////////////////////////////////////////////////////////////////////////////////
-
-#if HAVE_DUNE_PARMG
- BasisType feBasis(grid.levelGridView(grid.maxLevel()));
- DuneFunctionsBasis<BasisType> basis(feBasis);
-#else
- DuneFunctionsBasis<BasisType> basis(grid.leafGridView());
-#endif
OperatorAssembler<DuneFunctionsBasis<BasisType>,DuneFunctionsBasis<BasisType>,Dune::Partitions::Interior> operatorAssembler(basis, basis);
LaplaceAssembler<GridType,
@@ -129,19 +129,6 @@ setup(const typename BasisType::GridView::Grid& grid,
h1SemiNorm_ = std::make_shared< H1SemiNorm<CorrectionType> >(*A);
-#if HAVE_DUNE_PARMG
- innerMultigridStep_ = std::make_unique<Dune::ParMG::Multigrid<VectorType> >();
- innerMultigridStep_->mu_ = mu;
- innerMultigridStep_->preSmootherSteps_ = nu1;
- innerMultigridStep_->postSmootherSteps_ = nu2;
-#else
- innerSolver_ = std::make_shared<LoopSolver<CorrectionType> >(mmgStep,
- innerIterations_,
- innerTolerance_,
- h1SemiNorm_,
- Solver::REDUCED);
-#endif
-
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a mass matrix to create a norm that's equivalent to the L2-norm
// This will be used to monitor the gradient
@@ -167,6 +154,11 @@ setup(const typename BasisType::GridView::Grid& grid,
indices.exportIdx(*hessianMatrix_);
#if !HAVE_DUNE_PARMG
+ innerSolver_ = std::make_shared<LoopSolver<CorrectionType> >(mmgStep,
+ innerIterations_,
+ innerTolerance_,
+ h1SemiNorm_,
+ Solver::REDUCED);
// ////////////////////////////////////
// Create the transfer operators
// ////////////////////////////////////
@@ -248,9 +240,6 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
BasisType basis(grid_->levelGridView(grid_->maxLevel()));
BasisType coarseBasis(grid_->levelGridView(0));
std::vector<BoxConstraint<typename VectorType::field_type, blocksize>> coarseTrustRegionObstacles(coarseBasis.size());
- int numLevels = grid_->maxLevel()+1;
-
- auto& levelOp = mgSetup_->levelOps_;
#endif
MaxNormTrustRegion<blocksize> trustRegion(basis.size(), initialTrustRegionRadius_);
@@ -287,6 +276,10 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
*hessianMatrix_,
i==0 // assemble occupation pattern only for the first call
);
+#if HAVE_DUNE_PARMG
+ std::function<void(VectorType&)> accumulate = Dune::ParMG::makeAccumulate<VectorType>(*(mgSetup_->comms_.back()));
+ accumulate(rhs);
+#endif
rhs *= -1; // The right hand side is the _negative_ gradient
@@ -306,15 +299,13 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
if (this->verbosity_ == Solver::FULL)
std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
- // Transfer matrix data
- stiffnessMatrix = *hessianMatrix_;
-
#if HAVE_DUNE_PARMG
- Dune::ParMG::collectDiagonal(stiffnessMatrix, *mgSetup_->comms_.back());
+ if (!mgSetup_->overlap())
+ Dune::ParMG::collectDiagonal(*hessianMatrix_, *mgSetup_->comms_.back());
mgSetup_->matrix(hessianMatrix_);
-
- levelOp.back().maybeRestrictToMaster(rhs);
#endif
+ // Transfer matrix data
+ stiffnessMatrix = *hessianMatrix_;
recomputeGradientHessian = false;
@@ -323,17 +314,22 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
CorrectionType corr(rhs.size());
corr = 0;
+ //Take the obstacles on the finest grid and give them to the multigrid solver, it will create a hierarchy for all coarser grids
+ trustRegionObstacles = trustRegion.obstacles();
+
#if ! HAVE_DUNE_PARMG
mgStep->setProblem(stiffnessMatrix, corr, rhs);
- trustRegionObstacles = trustRegion.obstacles();
mgStep->setObstacles(trustRegionObstacles);
#else
- mgSetup_->setupLevelOps();
+ mgSetup_->setupObstacles(std::make_shared< std::vector<BoxConstraint<typename VectorType::field_type, blocksize>> >(trustRegionObstacles));
+ mgSetup_->setupSmoother(damping_);
+
+ auto& levelOp = mgSetup_->levelOps_;
bool enableCoarseCorrection = true;
if (enableCoarseCorrection)
- mgSetup_->setupCoarseIPOPTSolver();
+ mgSetup_->setupCoarseIPOPTSolver(baseTolerance_, baseIterations_);
else
mgSetup_->setupCoarseNullSolver();
@@ -343,29 +339,15 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
using ProjGS = Dune::ParMG::ParallelProjectedGS<MatrixType, VectorType>;
using namespace Dune::ParMG;
- auto& levelOp = mgSetup_->levelOps_;
- trustRegionObstacles = trustRegion.obstacles();
- for (int i = 0; i < coarseTrustRegionObstacles.size(); i++)
- coarseTrustRegionObstacles[i] = trustRegionObstacles[i];
- mgSetup_->setCoarseObstacles(coarseTrustRegionObstacles);
ProjGS projGs(&stiffnessMatrix, &trustRegionObstacles, mgSetup_->ignores_.back().get());
projGs.accumulate([&](VectorType& x) {
levelOp.back().maybeCopyFromMaster(x);
});
projGs.dampening(1.0);
- std::function<void(VectorType&)> collect = Dune::ParMG::makeCollect<VectorType>(*(mgSetup_->comms_.back()));
std::function<void(VectorType&)> restrictToMaster = [op=levelOp.back()](VectorType& x) { op.maybeRestrictToMaster(x); };
-
- auto energyFunctional = Dune::ParMG::makeParallelEnergyFunctional(
- stiffnessMatrix,
- rhs,
- grid_->comm(),
- //collect
- restrictToMaster
- );
- // matrix, b, dofmap, master);
+ restrictToMaster(rhs);
//Distributed energy norm
auto energyNorm = Dune::ParMG::parallelEnergyNorm<VectorType>(stiffnessMatrix, restrictToMaster, grid_->comm());
@@ -377,30 +359,16 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
levelOp.back().maybeCopyFromMaster(corr);
- double J = 0.0;
VectorType b;
unsigned step = 0;
- int activeObstacles;
MPI_Barrier(grid_->comm());
- auto realIterationStep = [&](VectorType& correction) {
- auto& project = projGs.project_;
- projGs.apply(correction, rhs);
- // prepare truncation
- auto currentIgnore = std::make_shared<Dune::BitSetVector<VectorType::block_type::dimension> >(*ignoreNodes_);
- auto extraIgnore = std::make_shared<Dune::BitSetVector<VectorType::block_type::dimension> >(ignoreNodes_->size(), false);
- project.ignoreActive(correction, *currentIgnore);
- project.ignoreActive(correction, *extraIgnore);
- mgSetup_->ignore(currentIgnore);
- assert(currentIgnore->count() - ignoreNodes_->count() == extraIgnore->count());
- activeObstacles = extraIgnore->count();
-
+ auto realIterationStep = [&](VectorType& innerIterate) {
// Create vector b to keep rhs untouched!
b = rhs;
- innerMultigridStep_->apply(correction, b);
-
- project(correction);
- J = energyFunctional(correction);
+ // The variable innerIterate is the correction that will be added
+ // to the current iterate in the trust region step.
+ innerMultigridStep_->apply(innerIterate, b);
++step;
};
@@ -410,22 +378,9 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
auto innerSolver_ = std::make_shared< LoopSolver<CorrectionType> >(iterationStep,
innerIterations_,
innerTolerance_,
- // h1SemiNorm_,
+ // h1SemiNorm_, //TODO: test this!
solverNorm,
NumProc::QUIET);
- innerSolver_->addCriterion(
- [&]() {
- return Dune::formatString(" % 12.10e", J);
- },
- " energy "
- );
-
- innerSolver_->addCriterion(
- [&]() {
- return Dune::formatString(" % 9d", activeObstacles);
- },
- " activeObs"
- );
#endif
innerSolver_->preprocess();
diff --git a/dune/elasticity/common/trustregionsolver.hh b/dune/elasticity/common/trustregionsolver.hh
index 1b87165..d819a2a 100644
--- a/dune/elasticity/common/trustregionsolver.hh
+++ b/dune/elasticity/common/trustregionsolver.hh
@@ -119,6 +119,12 @@ protected:
/** \brief Error tolerance of the multigrid QP solver */
double innerTolerance_;
+ /** \brief Maximum number of base solver iterations */
+ int baseIterations_;
+
+ /** \brief Error tolerance of the base solver inside the multigrid QP solver */
+ double baseTolerance_;
+
/** \brief Hessian matrix */
std::shared_ptr<MatrixType> hessianMatrix_;
--
GitLab