Commit 9271f845 authored by lisa_julia.nebel_at_tu-dresden.de's avatar lisa_julia.nebel_at_tu-dresden.de
Browse files

Use Patricks PN-Solver

parent 1d21458e
Pipeline #51189 failed with stage
in 7 minutes and 32 seconds
......@@ -30,9 +30,11 @@
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/proximalnewtonsolver.hh>
#include <dune/solvers/solvers/cholmodsolver.hh>
#include <dune/solvers/norms/twonorm.hh>
#include <dune/elasticity/common/trustregionsolver.hh>
#include <dune/elasticity/common/pnsolver.hh>
#include <dune/elasticity/assemblers/localadolcstiffness.hh>
#include <dune/elasticity/assemblers/feassembler.hh>
#include <dune/elasticity/materials/localintegralenergy.hh>
......@@ -52,6 +54,34 @@ const int order = 2;
using namespace Dune;
template<class VectorType>
struct SecondOrdersolver
{
const static int blocksize = VectorType::value_type::dimension;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize, blocksize> > MatrixType;
template<class BitVector, class G>
void minimize( const MatrixType& hess, const VectorType& grad, const G& g, double reg, VectorType& dx, const BitVector& /*ignore*/) const
{
// add reg
auto HH = hess;
for ( size_t i=0; i<dx.size(); i++ )
for ( size_t j=0; j<blocksize; j++ )
HH[i][i][j][j] += reg;
// add g to grad
auto gg = grad;
gg *= -1.0;
for ( size_t i=0; i<dx.size(); i++ )
gg[i] += 1.0;
Dune::Solvers::CholmodSolver<MatrixType,VectorType> cholmod( HH, dx, gg );
cholmod.solve();
}
};
int main (int argc, char *argv[]) try
{
// initialize MPI, finalize is done automatically on exit
......@@ -374,24 +404,19 @@ int main (int argc, char *argv[]) try
// Create a proximal newton solver
// /////////////////////////////////////////////////
ProximalNewtonSolver<PowerBasis,SolutionType> solver;
solver.setup(*grid,
&assembler,
x,
dirichletDofs,
tolerance,
maxSolverSteps,
initialRegularization
);
Solvers::ProximalNewton::ZeroFunctional<SolutionType> zeroFunctional;
Solvers::ProximalNewton::SimpleRegUpdater regUpdater;
SecondOrdersolver<SolutionType> sos;
TwoNorm<SolutionType> norm;
double reg = initialRegularization;
Solvers::ProximalNewtonSolver proximalNewtonSolver(assembler, zeroFunctional, sos, x, norm, regUpdater, reg, maxSolverSteps, tolerance, Solver::FULL);
proximalNewtonSolver.setIgnore(dirichletDofs);
// /////////////////////////////////////////////////////
// Solve!
// /////////////////////////////////////////////////////
solver.setInitialIterate(x);
solver.solve();
x = solver.getSol();
proximalNewtonSolver.solve();
}
/////////////////////////////////
......
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