Newer
Older
#include <dune/common/timer.hh>
#include <dune/common/parallel/mpihelper.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>

Oliver Sander
committed
// 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>
#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
#include <dune/solvers/solvers/iterativesolver.hh>

Oliver Sander
committed
#include "maxnormtrustregion.hh"
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/norms/h1seminorm.hh>
#include <dune/gfe/parallel/matrixcommunicator.hh>
#include <dune/gfe/parallel/vectorcommunicator.hh>

Oliver Sander
committed
template <class GridType, class TargetSpace>
void RiemannianTrustRegionSolver<GridType,TargetSpace>::
setup(const GridType& grid,
const GeodesicFEAssembler<BasisType, TargetSpace>* assembler,

Oliver Sander
committed
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
int maxTrustRegionSteps,
double initialTrustRegionRadius,
int multigridIterations,
double mgTolerance,

Oliver Sander
committed
int nu1,
int nu2,
int baseIterations,
double baseTolerance,
bool instrumented)
int rank = grid.comm().rank();

Oliver Sander
committed
grid_ = &grid;
assembler_ = assembler;

Oliver Sander
committed
this->tolerance_ = tolerance;
maxTrustRegionSteps_ = maxTrustRegionSteps;
initialTrustRegionRadius_ = initialTrustRegionRadius;

Oliver Sander
committed
innerIterations_ = multigridIterations;
innerTolerance_ = mgTolerance;
instrumented_ = instrumented;
ignoreNodes_ = &dirichletNodes;
int numLevels = grid_->maxLevel()+1;
//////////////////////////////////////////////////////////////////
// Create global numbering for matrix and vector transfer
//////////////////////////////////////////////////////////////////
guIndex_ = std::unique_ptr<GUIndex>(new GUIndex(grid_->leafGridView()));
// ////////////////////////////////
// Create a multigrid solver
// ////////////////////////////////
#ifdef HAVE_IPOPT
// First create an IPOpt base solver
QuadraticIPOptSolver<MatrixType, CorrectionType>* baseSolver = new QuadraticIPOptSolver<MatrixType,CorrectionType>;
baseSolver->verbosity_ = NumProc::QUIET;
baseSolver->tolerance_ = baseTolerance;
#else
#warning IPOpt not installed -- falling back onto a Gauss-Seidel base solver
// First create a Gauss-seidel base solver
TrustRegionGSStep<MatrixType, CorrectionType>* baseSolverStep = new TrustRegionGSStep<MatrixType, CorrectionType>;
EnergyNorm<MatrixType, CorrectionType>* baseEnergyNorm = new EnergyNorm<MatrixType, CorrectionType>(*baseSolverStep,1e-3);
::LoopSolver<CorrectionType>* baseSolver = new ::LoopSolver<CorrectionType>(baseSolverStep,

Oliver Sander
committed
baseIterations,
baseTolerance,
baseEnergyNorm,
Solver::QUIET);
// Transfer all Dirichlet data to the master processor
VectorCommunicator<GUIndex, Dune::BitSetVector<blocksize> > vectorComm(*guIndex_, 0);
Dune::BitSetVector<blocksize>* globalDirichletNodes = NULL;
globalDirichletNodes = new Dune::BitSetVector<blocksize>(vectorComm.reduceCopy(dirichletNodes));
// 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>;

Oliver Sander
committed
mmgStep->setMGType(mu, nu1, nu2);
mmgStep->ignoreNodes_ = globalDirichletNodes;
mmgStep->basesolver_ = baseSolver;
mmgStep->setSmoother(presmoother, postsmoother);

Oliver Sander
committed
mmgStep->obstacleRestrictor_= new MandelObstacleRestrictor<CorrectionType>();
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
// //////////////////////////////////////////////////////////////////////////////////////

Oliver Sander
committed
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;

Oliver Sander
committed
ScalarMatrixType localA;

Oliver Sander
committed
operatorAssembler.assemble(laplaceStiffness, localA);
if (h1SemiNorm_)
delete h1SemiNorm_;

Oliver Sander
committed
MatrixCommunicator<GUIndex, ScalarMatrixType> matrixComm(*guIndex_, 0);
ScalarMatrixType* A = new ScalarMatrixType(matrixComm.reduceAdd(localA));
h1SemiNorm_ = new H1SemiNorm<CorrectionType>(*A);
innerSolver_ = std::shared_ptr<LoopSolver<CorrectionType> >(new ::LoopSolver<CorrectionType>(mmgStep,
innerIterations_,
innerTolerance_,
h1SemiNorm_,
// Write all intermediate solutions, if requested

Oliver Sander
committed
if (instrumented_
&& dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get()))
dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get())->historyBuffer_ = "tmp/mgHistory";
// ////////////////////////////////////////////////////////////
// Create Hessian matrix and its occupation structure
// ////////////////////////////////////////////////////////////
hessianMatrix_ = std::auto_ptr<MatrixType>(new MatrixType);

Oliver Sander
committed
Dune::MatrixIndexSet indices(grid_->size(1), grid_->size(1));

Oliver Sander
committed
assembler_->getNeighborsPerVertex(indices);
indices.exportIdx(*hessianMatrix_);
// ////////////////////////////////////
// Create the transfer operators
// ////////////////////////////////////
for (size_t k=0; k<mmgStep->mgTransfer_.size(); 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);

Oliver Sander
committed
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
// and construct a single global transfer operator there.
typedef GlobalUniqueIndex<typename GridType::LeafGridView, gridDim> LeafP1GUIndex;
LeafP1GUIndex p1Index(grid_->leafGridView());
typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
MatrixCommunicator<GUIndex, TransferOperatorType, LeafP1GUIndex> matrixComm(*guIndex_, p1Index, 0);
mmgStep->mgTransfer_.back() = new PKtoP1MGTransfer<CorrectionType>
(Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(topTransferOp->getMatrix())));
for (int i=0; i<mmgStep->mgTransfer_.size()-1; i++){

Oliver Sander
committed
// Construct the local multigrid transfer matrix
TruncatedCompressedMGTransfer<CorrectionType>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType>;
newTransferOp->setup(*grid_,i+1,i+2);
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
// and construct a single global transfer operator there.
typedef GlobalUniqueIndex<typename GridType::LevelGridView, gridDim> LevelGUIndex;
LevelGUIndex fineGUIndex(grid_->levelGridView(i+2));
LevelGUIndex coarseGUIndex(grid_->levelGridView(i+1));
typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
MatrixCommunicator<LevelGUIndex, TransferOperatorType> matrixComm(fineGUIndex, coarseGUIndex, 0);
mmgStep->mgTransfer_[i] = new TruncatedCompressedMGTransfer<CorrectionType>
(Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(newTransferOp->getMatrix())));
}
}
#else
for (size_t i=0; i<mmgStep->mgTransfer_.size(); i++){
// Construct the local multigrid transfer matrix

Oliver Sander
committed
TruncatedCompressedMGTransfer<CorrectionType>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType>;
newTransferOp->setup(*grid_,i,i+1);
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
// and construct a single global transfer operator there.
typedef GlobalUniqueIndex<typename GridType::LevelGridView, gridDim> LevelGUIndex;
LevelGUIndex fineGUIndex(grid_->levelGridView(i+1));
LevelGUIndex coarseGUIndex(grid_->levelGridView(i));
typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
MatrixCommunicator<LevelGUIndex, TransferOperatorType> matrixComm(fineGUIndex, coarseGUIndex, 0);
mmgStep->mgTransfer_[i] = new TruncatedCompressedMGTransfer<CorrectionType>
(Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(newTransferOp->getMatrix())));
// //////////////////////////////////////////////////////////
// Create obstacles
// //////////////////////////////////////////////////////////
hasObstacle_.resize(guIndex_->nGlobalEntity(), true);
mmgStep->hasObstacle_ = &hasObstacle_;

Oliver Sander
committed

Oliver Sander
committed
template <class GridType, class TargetSpace>
void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
int argc = 0;
char** argv;
Dune::MPIHelper& mpiHelper = Dune::MPIHelper::instance(argc,argv);
int rank = mpiHelper.rank();

Oliver Sander
committed
MonotoneMGStep<MatrixType,CorrectionType>* mgStep = NULL;

Oliver Sander
committed
// 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_);

Oliver Sander
committed
MaxNormTrustRegion<blocksize> trustRegion(guIndex_->nGlobalEntity(), initialTrustRegionRadius_);

Oliver Sander
committed
std::vector<std::vector<BoxConstraint<field_type,blocksize> > > trustRegionObstacles((mgStep)

Oliver Sander
committed
? mgStep->numLevels()
: 0);

Oliver Sander
committed
// /////////////////////////////////////////////////////
// Set up the log file, if requested
// /////////////////////////////////////////////////////
if (instrumented_) {
fp = fopen("statistics", "w");
if (!fp)
DUNE_THROW(Dune::IOError, "Couldn't open statistics file for writing!");
// /////////////////////////////////////////////////////
// Trust-Region Solver
// /////////////////////////////////////////////////////

Oliver Sander
committed
double oldEnergy = assembler_->computeEnergy(x_);
oldEnergy = mpiHelper.getCollectiveCommunication().sum(oldEnergy);
bool recomputeGradientHessian = true;
CorrectionType rhs;
MatrixType stiffnessMatrix;
CorrectionType rhs_global;
VectorCommunicator<GUIndex, CorrectionType> vectorComm(*guIndex_, 0);
MatrixCommunicator<GUIndex, MatrixType> matrixComm(*guIndex_, 0);
for (int i=0; i<maxTrustRegionSteps_; i++) {
/* std::cout << "current iterate:\n";
std::cout << x_[j] << std::endl;*/
if (this->verbosity_ == Solver::FULL and rank==0) {
std::cout << "----------------------------------------------------" << std::endl;
std::cout << " Trust-Region Step Number: " << i
<< ", radius: " << trustRegion.radius()
<< ", energy: " << oldEnergy << std::endl;
std::cout << "----------------------------------------------------" << std::endl;
}
CorrectionType corr(x_.size());
corr = 0;
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
if (this->verbosity_ == Solver::FULL)
std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
// Transfer matrix data
stiffnessMatrix = matrixComm.reduceAdd(*hessianMatrix_);
// Transfer vector data
rhs_global = vectorComm.reduceAdd(rhs);
recomputeGradientHessian = false;
CorrectionType corr_global(rhs_global.size());
corr_global = 0;
if (rank==0)
{
mgStep->setProblem(stiffnessMatrix, corr_global, rhs_global);
trustRegionObstacles.back() = trustRegion.obstacles();
mgStep->obstacles_ = &trustRegionObstacles.back();
innerSolver_->preprocess();
///////////////////////////////
// Solve !
///////////////////////////////
std::cout << "Solve quadratic problem..." << std::endl;
innerSolver_->solve();
if (mgStep)
corr_global = mgStep->getSol();
//std::cout << "Correction: " << std::endl << corr_global << std::endl;
}
// Distribute solution
std::cout << "Transfer solution back to root process ..." << std::endl;
corr = vectorComm.scatter(corr_global);
if (instrumented_) {
fprintf(fp, "Trust-region step: %d, trust-region radius: %g\n",
i, trustRegion.radius());
// ///////////////////////////////////////////////////////////////
// Compute and measure progress against the exact solution
// for each trust region step
// ///////////////////////////////////////////////////////////////
CorrectionType exactSolution = corr;
// Start from 0
double oldError = 0;
double totalConvRate = 1;
double convRate = 1;
// Write statistics of the initial solution
// Compute the energy norm
oldError = h1SemiNorm_->operator()(exactSolution);

Oliver Sander
committed
for (int j=0; j<innerIterations_; j++) {
// read iteration from file
CorrectionType intermediateSol(grid_->size(gridDim));
intermediateSol = 0;
char iSolFilename[100];
sprintf(iSolFilename, "tmp/mgHistory/intermediatesolution_%04d", j);
FILE* fpInt = fopen(iSolFilename, "rb");
if (!fpInt)
DUNE_THROW(Dune::IOError, "Couldn't open intermediate solution");
for (size_t k=0; k<intermediateSol.size(); k++)
for (int l=0; l<blocksize; l++)
fread(&intermediateSol[k][l], sizeof(double), 1, fpInt);
fclose(fpInt);
//std::cout << "intermediateSol\n" << intermediateSol << std::endl;
// Compute errors
intermediateSol -= exactSolution;
//std::cout << "error\n" << intermediateSol << std::endl;
// Compute the H1 norm
double error = h1SemiNorm_->operator()(intermediateSol);
convRate = error / oldError;
totalConvRate *= convRate;
if (error < 1e-12)
break;
std::cout << "Iteration: " << j << " ";
std::cout << "Errors: error " << error << ", convergence rate: " << convRate
<< ", total conv rate " << pow(totalConvRate, 1/((double)j+1)) << std::endl;
fprintf(fp, "%d %g %g %g\n", j+1, error, convRate, pow(totalConvRate, 1/((double)j+1)));
oldError = error;
if (this->verbosity_ == NumProc::FULL)
std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;

Oliver Sander
committed
if (corr.infinity_norm() < this->tolerance_) {
if (this->verbosity_ == NumProc::FULL and rank==0)
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;
// ////////////////////////////////////////////////////
// Check whether trust-region step can be accepted
// ////////////////////////////////////////////////////
for (size_t j=0; j<newIterate.size(); j++)

Oliver Sander
committed
newIterate[j] = TargetSpace::exp(newIterate[j], corr[j]);
double energy = assembler_->computeEnergy(newIterate);
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.size());
tmp = 0;
double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
modelDecrease = mpiHelper.getCollectiveCommunication().sum(modelDecrease);
double relativeModelDecrease = modelDecrease / std::fabs(energy);
if (this->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 (energy >= oldEnergy and rank==0) {

Oliver Sander
committed
if (this->verbosity_ == NumProc::FULL)
printf("Richtung ist keine Abstiegsrichtung!\n");

Oliver Sander
committed
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;

Oliver Sander
committed
x_ = newIterate;
break;
}
// //////////////////////////////////////////////
// Check for acceptance of the step
// //////////////////////////////////////////////
if ( (oldEnergy-energy) / modelDecrease > 0.9) {
// very successful iteration
trustRegion.scale(2);

Oliver Sander
committed
// 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;

Oliver Sander
committed
// current energy becomes 'oldEnergy' for the next iteration
oldEnergy = energy;
recomputeGradientHessian = true;
// Decrease the trust-region radius
trustRegion.scale(0.5);
if (this->verbosity_ == NumProc::FULL and rank==0)
std::cout << "Unsuccessful iteration!" << std::endl;
// /////////////////////////////////////////////////////////////////////
// Write the iterate to disk for later convergence rate measurement
// /////////////////////////////////////////////////////////////////////
if (instrumented_) {

Oliver Sander
committed
char iFilename[100];
sprintf(iFilename, "tmp/intermediateSolution_%04d", i);

Oliver Sander
committed
FILE* fpIterate = fopen(iFilename, "wb");
if (!fpIterate)
DUNE_THROW(SolverError, "Couldn't open file " << iFilename << " for writing");
for (size_t j=0; j<x_.size(); j++)

Oliver Sander
committed
fwrite(&x_[j], sizeof(TargetSpace), 1, fpIterate);

Oliver Sander
committed
fclose(fpIterate);
if (rank==0)
std::cout << "iteration took " << totalTimer.elapsed() << " sec." << std::endl;
// //////////////////////////////////////////////
// Close logfile
// //////////////////////////////////////////////
if (instrumented_)
fclose(fp);