Commit 43753f66 authored by lisa_julia.nebel_at_tu-dresden.de's avatar lisa_julia.nebel_at_tu-dresden.de
Browse files

HACK: Comment out PARMG in TR-Solver

parent f279a716
Pipeline #51164 failed with stage
in 12 minutes and 27 seconds
......@@ -18,7 +18,7 @@
// 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>
#include <dune/solvers/transferoperators/truncateddensemgtransfer.hh>
#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/twonorm.hh>
......@@ -64,7 +64,7 @@ setup(const typename BasisType::GridView::Grid& grid,
const auto dim = VectorType::value_type::dimension;
#if HAVE_DUNE_PARMG
#if 0 // HAVE_DUNE_PARMG
Dune::Timer setupTimer;
mgSetup_ = std::make_unique<MGSetup>(grid);
mgSetup_->overlap(false);
......@@ -249,7 +249,7 @@ setup(const typename BasisType::GridView::Grid& grid,
assembler_->getNeighborsPerVertex(indices);
indices.exportIdx(*hessianMatrix_);
}
#if !HAVE_DUNE_PARMG
#if 1 // !HAVE_DUNE_PARMG
innerSolver_ = std::make_shared<LoopSolver<CorrectionType> >(mmgStep,
innerIterations_,
innerTolerance_,
......@@ -290,8 +290,8 @@ setup(const typename BasisType::GridView::Grid& grid,
int numLevels = grid_->maxLevel()+1;
using TransferOperatorType = typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType;
std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<CorrectionType>>> transferOperators(isP1Basis ? numLevels-1 : numLevels);
using TransferOperatorType = typename TruncatedDenseMGTransfer<CorrectionType>::TransferOperatorType;
std::vector<std::shared_ptr<TruncatedDenseMGTransfer<CorrectionType>>> transferOperators(isP1Basis ? numLevels-1 : numLevels);
// Here we set up the restriction of the leaf grid space into the leaf grid P1/Q1 space
// TODO: IIRC that move to P1Basis was only implemented because I didn't have general
......@@ -318,17 +318,17 @@ setup(const typename BasisType::GridView::Grid& grid,
DuneFunctionsBasis<BasisType> >(pkToP1TransferMatrix,p1Basis,basis);
}
auto pkToP1 = std::make_shared< TruncatedCompressedMGTransfer<CorrectionType> >();
auto pkToP1 = std::make_shared< TruncatedDenseMGTransfer<CorrectionType> >();
transferOperators.back() = pkToP1;
std::shared_ptr<TransferOperatorType> topTransferOperator = std::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
pkToP1->setMatrix(topTransferOperator);
pkToP1->setMatrix(*topTransferOperator);
}
// Now the P1/Q1 restriction operators for the remaining levels
for (int i=0; i<numLevels-1; i++) {
// Construct the local multigrid transfer matrix
transferOperators[i] = std::make_shared<TruncatedCompressedMGTransfer<CorrectionType> >();
transferOperators[i] = std::make_shared<TruncatedDenseMGTransfer<CorrectionType> >();
transferOperators[i]->setup(*grid_,i,i+1);
}
......@@ -349,7 +349,7 @@ template <class BasisType, class VectorType>
void TrustRegionSolver<BasisType,VectorType>::solve()
{
int rank = grid_->comm().rank();
#if !HAVE_DUNE_PARMG
#if 1 // !HAVE_DUNE_PARMG
MonotoneMGStep<MatrixType,CorrectionType>* mgStep = nullptr;
// if the inner solver is a monotone multigrid set up a max-norm trust-region
......@@ -400,7 +400,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
*hessianMatrix_,
i==0 // assemble occupation pattern only for the first call
);
#if HAVE_DUNE_PARMG
#if 0 // HAVE_DUNE_PARMG
std::function<void(VectorType&)> accumulate = Dune::ParMG::makeAccumulate<VectorType>(*(mgSetup_->comms_.back()));
accumulate(rhs);
#endif
......@@ -423,7 +423,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
if (this->verbosity_ == Solver::FULL)
std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
#if HAVE_DUNE_PARMG
#if 0 // HAVE_DUNE_PARMG
if (!mgSetup_->overlap())
Dune::ParMG::collectDiagonal(*hessianMatrix_, *mgSetup_->comms_.back());
mgSetup_->matrix(hessianMatrix_);
......@@ -442,7 +442,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
//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
#if 1 // !HAVE_DUNE_PARMG
mgStep->setProblem(stiffnessMatrix, corr, rhs);
mgStep->setObstacles(trustRegionObstacles);
......@@ -533,7 +533,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
#if ! HAVE_DUNE_PARMG
#if 1 // !HAVE_DUNE_PARMG
if (mgStep)
corr = mgStep->getSol();
......@@ -549,19 +549,15 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
if (this->verbosity_ == NumProc::FULL)
std::cout << "Infinity norm of the correction: " << correctionInfinityNorm << std::endl;
if (correctionInfinityNorm < this->tolerance_) {
if (correctionInfinityNorm < this->tolerance_ && this->tolerance_*100000 < trustRegion.radius() ) {
if (this->verbosity_ == NumProc::FULL and rank==0)
{
if (correctionInfinityNorm < trustRegion.radius())
std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
else
std::cout << "TRUST-REGION UNDERFLOW!" << std::endl;
}
std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
if (this->verbosity_ != NumProc::QUIET and rank==0)
std::cout << i+1 << " trust-region steps were taken." << std::endl;
break;
}
} else if (correctionInfinityNorm < this->tolerance_ )
std::cout << "TRUST-REGION UNDERFLOW!" << std::endl;
// ////////////////////////////////////////////////////
// Check whether trust-region step can be accepted
......@@ -611,7 +607,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
printf("Richtung ist keine Abstiegsrichtung!\n");
}
if (energy >= oldEnergy && (std::fabs(relativeFunctionalDecrease) < 1e-9 || std::fabs(relativeModelDecrease) < 1e-9)) {
/*if (energy >= oldEnergy && (std::fabs(relativeFunctionalDecrease) < 1e-9 || std::fabs(relativeModelDecrease) < 1e-9)) {
if (this->verbosity_ == NumProc::FULL and rank==0)
std::cout << "Suspecting rounding problems" << std::endl;
......@@ -620,7 +616,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
x_ = newIterate;
break;
}
}*/
}
}
......@@ -656,6 +652,9 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
if (this->verbosity_ == NumProc::FULL and rank==0)
std::cout << "Unsuccessful iteration!" << std::endl;
if (trustRegion.radius() < 1e-9 && this->verbosity_ == NumProc::FULL and rank==0)
std::cout << "The radius is too small to continue with a meaningful calculation!" << std::endl;
}
std::cout << "iteration took " << totalTimer.elapsed() << " sec." << std::endl;
......
......@@ -25,7 +25,7 @@
#include <dune/elasticity/assemblers/feassembler.hh>
#if HAVE_DUNE_PARMG
#if 0 // HAVE_DUNE_PARMG
// enable a temporary hack to let parmg work with powerBases
// which is needed for the current implementation of elasticty
......@@ -65,7 +65,7 @@ class TrustRegionSolver
typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > CorrectionType;
typedef VectorType SolutionType;
#if HAVE_DUNE_PARMG
#if 0 // HAVE_DUNE_PARMG
using MGSetup = Dune::ParMG::ParallelMultiGridSetup<BasisType, MatrixType, VectorType>;
#endif
public:
......@@ -98,7 +98,7 @@ public:
void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
{
ignoreNodes_ = &ignoreNodes;
#if !HAVE_DUNE_PARMG
#if 1 // !HAVE_DUNE_PARMG
std::shared_ptr<LoopSolver<CorrectionType> > loopSolver = std::dynamic_pointer_cast<LoopSolver<CorrectionType> >(innerSolver_);
assert(loopSolver);
loopSolver->iterationStep_->ignoreNodes_ = ignoreNodes_;
......@@ -151,7 +151,7 @@ protected:
const Dune::Elasticity::FEAssembler<BasisType,VectorType>*,
const Dune::FEAssembler<DuneFunctionsBasis<BasisType>,VectorType>* > assembler_;
#if HAVE_DUNE_PARMG
#if 0 // HAVE_DUNE_PARMG
/** \brief The solver for the inner problem */
std::unique_ptr<MGSetup> mgSetup_;
......@@ -161,7 +161,7 @@ protected:
/** \brief The solver for the quadratic inner problems */
std::shared_ptr<Solver> innerSolver_;
#if ! HAVE_DUNE_PARMG
#if 1 // !HAVE_DUNE_PARMG
/** \brief Contains 'true' everywhere -- the trust-region is bounded */
Dune::BitSetVector<blocksize> hasObstacle_;
#endif
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment