Newer
Older
Lisa Julia Nebel
committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
#ifndef RIEMANNIAN_PROXIMAL_NEWTON_SOLVER_HH
#define RIEMANNIAN_PROXIMAL_NEWTON_SOLVER_HH
#include <vector>
#include <dune/common/bitsetvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/solvers/common/boxconstraint.hh>
#include <dune/solvers/norms/h1seminorm.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/solvers/cholmodsolver.hh>
#include <dune/gfe/periodic1dpq1nodalbasis.hh>
#include "geodesicfeassembler.hh"
#include "riemanniantrsolver.hh"
#include <dune/grid/utility/globalindexset.hh>
#include <dune/gfe/parallel/globalmapper.hh>
#include <dune/gfe/parallel/globalp1mapper.hh>
#include <dune/gfe/parallel/globalp2mapper.hh>
#include <dune/gfe/parallel/p2mapper.hh>
/** \brief Riemannian proximal-newton solver for geodesic finite-element problems */
template <class Basis, class TargetSpace>
class RiemannianProximalNewtonSolver
: public IterativeSolver<std::vector<TargetSpace>,
Dune::BitSetVector<TargetSpace::TangentVector::dimension> >
{
typedef typename Basis::GridView::Grid GridType;
const static int blocksize = TargetSpace::TangentVector::dimension;
const static int gridDim = GridType::dimension;
// Centralize the field type here
typedef double field_type;
// Some types that I need
typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize, blocksize> > MatrixType;
typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > CorrectionType;
typedef std::vector<TargetSpace> SolutionType;
#if HAVE_MPI
typedef typename MapperFactory<Basis>::GlobalMapper GlobalMapper;
typedef typename MapperFactory<Basis>::LocalMapper LocalMapper;
#endif
/** \brief Records information about the last run of the RiemannianProximalNewtonSolver
*
* This is used primarily for unit testing.
*/
struct Statistics
{
std::size_t finalIteration;
field_type finalEnergy;
};
public:
RiemannianProximalNewtonSolver()
: IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
hessianMatrix_(nullptr), h1SemiNorm_(NULL)
{}
/** \brief Set up the solver using a choldmod solver as the inner solver */
void setup(const GridType& grid,
const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
int maxProximalNewtonSteps,
double initialRegularization,
bool instrumented);
void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
{
ignoreNodes_ = &ignoreNodes;
innerSolver_->ignoreNodes_ = ignoreNodes_;
}
void solve();
void setInitialSolution(const SolutionType& x) DUNE_DEPRECATED {
x_ = x;
}
void setInitialIterate(const SolutionType& x) {
x_ = x;
}
SolutionType getSol() const {return x_;}
const Statistics& getStatistics() const {return statistics_;}
protected:
#if HAVE_MPI
std::unique_ptr<GlobalMapper> globalMapper_;
#endif
/** \brief The grid */
const GridType* grid_;
/** \brief The solution vector */
SolutionType x_;
/** \brief The initial regularization parameter for the proximal newton step */
double initialRegularization_;
double tolerance_;
/** \brief Maximum number of proximal-newton steps */
std::size_t maxProximalNewtonSteps_;
/** \brief Hessian matrix */
std::unique_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */
const GeodesicFEAssembler<Basis, TargetSpace>* assembler_;
/** \brief The solver for the quadratic inner problems */
std::shared_ptr<typename Dune::Solvers::CholmodSolver<MatrixType,CorrectionType>> innerSolver_;
/** \brief The Dirichlet nodes */
const Dune::BitSetVector<blocksize>* ignoreNodes_;
/** \brief The norm used to measure convergence for statistics*/
std::shared_ptr<H1SemiNorm<CorrectionType> > h1SemiNorm_;
/** \brief An L2-norm, really. The H1SemiNorm class is badly named */
std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_;
/** \brief If set to true we log convergence speed and other stuff */
bool instrumented_;
/** \brief Store information about solver runs for unit testing */
Statistics statistics_;
#if DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
std::shared_ptr<Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > > A;
std::shared_ptr<Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > > massMatrix;
#endif
};
#include "riemannianpnsolver.cc"
#endif