Skip to content
Snippets Groups Projects

Add a simple proximal newton solver

1 unresolved thread
4 files
+ 627
0
Compare changes
  • Side-by-side
  • Inline
Files
4
+ 275
0
 
#include <dune/common/bitsetvector.hh>
 
#include <dune/common/timer.hh>
 
 
#include <dune/istl/io.hh>
 
#include <dune/istl/scaledidmatrix.hh>
 
 
#include <dune/functions/functionspacebases/lagrangebasis.hh>
 
 
#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
 
#include <dune/fufem/assemblers/istlbackend.hh>
 
 
#if DUNE_VERSION_GTE(DUNE_SOLVERS, 2, 8)
 
// Using a cholmod solver as the inner solver, available only since 2.8
 
#include <dune/solvers/solvers/cholmodsolver.hh>
 
#else
 
// Using a umfpack solver as the inner solver
 
#include <dune/solvers/solvers/umfpacksolver.hh>
 
#endif
 
 
template <class BasisType, class VectorType>
 
void ProximalNewtonSolver<BasisType,VectorType>::
 
setup(const typename BasisType::GridView::Grid& grid,
 
std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
 
Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
 
const Dune::Elasticity::FEAssembler<BasisType,VectorType>*,
 
const Dune::FEAssembler<DuneFunctionsBasis<BasisType>,VectorType>* > assembler,
 
const SolutionType& x,
 
const Dune::BitSetVector<blocksize>& dirichletNodes,
 
double tolerance,
 
int maxProximalNewtonSteps,
 
double initialRegularization)
 
 
{
 
grid_ = &grid;
 
assembler_ = assembler;
 
x_ = x;
 
this->tolerance_ = tolerance;
 
maxProximalNewtonSteps_ = maxProximalNewtonSteps;
 
initialRegularization_ = initialRegularization;
 
ignoreNodes_ = std::make_shared<Dune::BitSetVector<blocksize>>(dirichletNodes);
 
 
const auto dim = VectorType::value_type::dimension;
 
 
//////////////////////////////////////////////////////////////////
 
// Create the inner solver using a cholmod (or umfpack) solver
 
//////////////////////////////////////////////////////////////////
 
#if DUNE_VERSION_GTE(DUNE_SOLVERS, 2, 8)
 
innerSolver_ = std::make_shared<Dune::Solvers::CholmodSolver<MatrixType,CorrectionType> >();
 
#else
 
std::cout << "using umfpacksolver" << std::endl;
 
innerSolver_ = std::make_shared<Dune::Solvers::UMFPackSolver<MatrixType,CorrectionType> >();
 
#endif
 
auto globalDirichletNodes = new Dune::BitSetVector<blocksize>(dirichletNodes);
 
innerSolver_->setIgnore(*globalDirichletNodes);
 
 
// ////////////////////////////////////////////////////////////
 
// Create pointer to Hessian matrix and set its occupation structure
 
// ////////////////////////////////////////////////////////////
 
std::conditional_t< // do we have a dune-functions basis?
 
Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
 
BasisType,
 
DuneFunctionsBasis<BasisType> > basis(assembler_->basis_);
 
std::conditional_t< // do we have a dune-functions basis?
 
Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
 
Dune::Fufem::DuneFunctionsOperatorAssembler<BasisType,BasisType>,
 
OperatorAssembler<DuneFunctionsBasis<BasisType>,DuneFunctionsBasis<BasisType>,Dune::Partitions::Interior> > operatorAssembler(basis,basis);
 
 
hessianMatrix_ = std::make_shared<MatrixType>();
 
if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
 
{
 
auto hessianBackend = Dune::Fufem::istlMatrixBackend(*hessianMatrix_);
 
auto hessianPatternBuilder = hessianBackend.patternBuilder();
 
operatorAssembler.assembleBulkPattern(hessianPatternBuilder);
 
hessianPatternBuilder.setupMatrix();
 
}
 
else
 
{
 
Dune::MatrixIndexSet indices(grid_->size(1), grid_->size(1));
 
assembler_->getNeighborsPerVertex(indices);
 
indices.exportIdx(*hessianMatrix_);
 
}
 
}
 
 
 
template <class BasisType, class VectorType>
 
void ProximalNewtonSolver<BasisType,VectorType>::solve()
 
{
 
double oldEnergy = assembler_->computeEnergy(x_);
 
 
bool recomputeGradientHessian = true;
 
CorrectionType rhs;
 
MatrixType stiffnessMatrix;
 
 
Dune::Timer problemTimer;
 
 
double regularization = initialRegularization_;
 
 
for (int i=0; i<maxProximalNewtonSteps_; i++) {
 
 
Dune::Timer totalTimer;
 
if (this->verbosity_ == Solver::FULL) {
 
std::cout << "----------------------------------------------------" << std::endl;
 
std::cout << " ProximalNewton Step Number: " << i
 
<< ", regularization parameter: " << regularization
 
<< ", 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 (size_t k=0; k<gradient[j].size(); k++)
 
if ((*ignoreNodes_)[j][k])
 
gradient[j][k] = 0;
 
 
if (this->verbosity_ == Solver::FULL)
 
std::cout << "Gradient infinity norm: " << gradient.infinity_norm() << 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;
 
bool solvedByInnerSolver = true;
 
 
// Add the regularization - Identity Matrix for now
 
for (std::size_t i=0; i<stiffnessMatrix.N(); i++)
 
for(int j=0; j<blocksize; j++)
 
stiffnessMatrix[i][i][j][j] += regularization;
 
 
innerSolver_->setProblem(stiffnessMatrix,corr,rhs);
 
innerSolver_->preprocess();
 
 
///////////////////////////////
 
// Solve !
 
///////////////////////////////
 
 
std::cout << "Solve quadratic problem..." << std::endl;
 
 
Dune::Timer solutionTimer;
 
 
try {
 
innerSolver_->solve();
 
} catch (Dune::Exception &e) {
 
std::cerr << "Error while solving: " << e << std::endl;
 
solvedByInnerSolver = false;
 
corr = 0;
 
}
 
double energy = 0;
 
double totalModelDecrease = 0;
 
 
double correctionInfinityNorm = corr.infinity_norm();
 
if (std::isnan(correctionInfinityNorm))
 
solvedByInnerSolver = false;
 
 
SolutionType newIterate = x_;
 
 
if (this->verbosity_ == NumProc::FULL)
 
std::cout << "Infinity norm of the correction: " << correctionInfinityNorm << std::endl;
 
 
if (correctionInfinityNorm < this->tolerance_ && this->tolerance_*100000 < 1.0/regularization ) {
 
if (this->verbosity_ == NumProc::FULL)
 
std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
 
 
if (this->verbosity_ != NumProc::QUIET)
 
std::cout << i+1 << " proximal-newton steps were taken." << std::endl;
 
break;
 
} else if (correctionInfinityNorm < this->tolerance_ )
 
std::cout << "regularization overflow!" << std::endl;
 
 
// ////////////////////////////////////////////////////
 
// Check whether proximal-newton step can be accepted
 
// ////////////////////////////////////////////////////
 
 
for (size_t j=0; j<newIterate.size(); j++)
 
newIterate[j] += corr[j];
 
 
try {
 
energy = assembler_->computeEnergy(newIterate);
 
} catch (Dune::Exception &e) {
 
std::cerr << "Error while computing the energy of the new iterate: " << e << std::endl;
 
std::cerr << "Redoing proximal-newton step with higher regularization parameter..." << std::endl;
 
newIterate = x_;
 
solvedByInnerSolver = false;
 
energy = oldEnergy;
 
}
 
 
if (!solvedByInnerSolver) {
 
energy = oldEnergy;
 
newIterate = x_;
 
} else {
 
// 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);
 
 
assert(modelDecrease >= 0);
 
 
double relativeModelDecrease = modelDecrease / std::fabs(energy);
 
double relativeFunctionalDecrease = (oldEnergy - energy)/std::fabs(energy);
 
 
if (this->verbosity_ == NumProc::FULL) {
 
std::cout << "Absolute model decrease: " << modelDecrease << std::endl;
 
std::cout << "Functional decrease: " << oldEnergy - energy << std::endl;
 
std::cout << "oldEnergy = " << oldEnergy << " and new energy = " << energy << std::endl;
 
std::cout << "Relative model decrease: " << relativeModelDecrease << std::endl;
 
std::cout << "Relative functional decrease: " << relativeFunctionalDecrease << std::endl;
 
}
 
 
 
if (energy >= oldEnergy) {
 
if (this->verbosity_ == NumProc::FULL)
 
printf("Correction caused an energy increase!\n");
 
}
 
}
 
 
// //////////////////////////////////////////////
 
// Check for acceptance of the step
 
// //////////////////////////////////////////////
 
if (solvedByInnerSolver && energy < oldEnergy && (oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.9) {
 
// very successful iteration
 
 
x_ = newIterate;
 
regularization *= 0.5;
 
 
// current energy becomes 'oldEnergy' for the next iteration
 
oldEnergy = energy;
 
 
recomputeGradientHessian = true;
 
 
} else if (solvedByInnerSolver && ((oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.01 || std::fabs(oldEnergy-energy) < 1e-12)) {
 
// successful iteration
 
x_ = newIterate;
 
 
// current energy becomes 'oldEnergy' for the next iteration
 
oldEnergy = energy;
 
 
recomputeGradientHessian = true;
 
 
} else {
 
 
// unsuccessful iteration
 
 
// Increase the regularization parameter
 
regularization *= 2;
 
 
if (this->verbosity_ == NumProc::FULL)
 
std::cout << "Unsuccessful iteration!" << std::endl;
 
}
 
 
std::cout << "iteration took " << totalTimer.elapsed() << " sec." << std::endl;
 
}
 
std::cout << "The whole proximal-newton step took " << problemTimer.elapsed() << " sec." << std::endl;
 
 
}
Loading