Skip to content
Snippets Groups Projects
Commit 16374b9e authored by lisa_julia.nebel_at_tu-dresden.de's avatar lisa_julia.nebel_at_tu-dresden.de Committed by oliver.sander_at_tu-dresden.de
Browse files

Cleanup of the TrustRegionSolver class 

Mainly issues related to dune-parmg.
parent 91552bbb
No related branches found
No related tags found
1 merge request!23Various minor improvements by Lisa
Hide whitespace changes
Inline Side-by-side
dune/elasticity/common/trustregionsolver.cc
+ 49
− 94
View file @ 16374b9e
......@@ -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();
......
dune/elasticity/common/trustregionsolver.hh
+ 6
− 0
View file @ 16374b9e
......@@ -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_;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment