dune-fem 2.8.0
Loading...
Searching...
No Matches
linearized.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_SCHEMES_LINEARIZED_HH
2#define DUNE_FEM_SCHEMES_LINEARIZED_HH
3
4#include <cstddef>
5
6#include <tuple>
7#include <type_traits>
8#include <utility>
9#include <vector>
10
13
14
15namespace Dune
16{
17
18 namespace Fem
19 {
20 // Given non linear operator L[u;(c)] where (c) are some coefficients
21 // we define a linear operator
22 // A[u;(c)] := L[ubar;(c)] + DL[ubar;(c)](u-ubar)
23 // = DL[ubar;(c)]u + L[ubar;(c)] - DL[ubar;(c)]ubar
24 // = DL[ubar;(c)]u + b
25 // Note: DA[ubar,(c)]u = DL[ubar,(c)]u
26 // and if L[u];(c)] = Au + c we have
27 // A[u;(c)] = Au + Aubar + c - Aubar = Au + c
28 //
29 // Use in Newton method:
30 // DL[ubar]d + L[ubar] = 0 ; d = -DL^{-1}L[ubar] ; u = u+d
31 // u = u-DL^{-1}L[ubar]
32 // A[u] = DL[ubar]u + L[ubar] - DL[ubar]ubar = 0
33 // u = -DL^{-1}(L[ubar]-DL ubar)
34 // u = ubar-DL^{-1}L[ubar]
35 template< class Scheme >
37 {
38 typedef Scheme SchemeType;
39 typedef typename SchemeType::DiscreteFunctionType DiscreteFunctionType;
40 typedef typename SchemeType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
41 typedef typename SchemeType::GridPartType GridPartType;
42 typedef typename SchemeType::LinearOperatorType LinearOperatorType;
43 typedef typename SchemeType::InverseOperatorType LinearInverseOperatorType;
44 typedef typename SchemeType::ModelType ModelType;
46 {
49 {}
50
53 };
54
57 : scheme_( scheme ),
58 linearOperator_( "linearized Op", scheme.space(), scheme.space() ),
59 linabstol_( 1e-12 ), // parameter.getValue< double >("linabstol", 1e-12 ) ),
60 linreduction_( 1e-12 ), // parameter_.getValue< double >("linreduction", 1e-12 ) ),
61 maxIter_( 1000 ),
62 inverseOperator_( nullptr ),
63 affineShift_( "affine shift", scheme.space() ),
64 parameter_( std::move( parameter ) ),
65 ubar_("ubar", scheme.space())
66 {
67 ubar_.clear();
68 setup(ubar_);
69 }
72 : scheme_( scheme ),
73 linearOperator_( "linearized Op", scheme.space(), scheme.space() ),
74 linabstol_( 1e-12 ), // parameter.getValue< double >("linabstol", 1e-12 ) ),
75 linreduction_( 1e-12 ), // parameter_.getValue< double >("linreduction", 1e-12 ) ),
76 maxIter_( 1000 ),
77 inverseOperator_( nullptr ),
78 affineShift_( "affine shift", scheme.space() ),
79 parameter_( std::move( parameter ) ),
80 ubar_("ubar", scheme.space())
81 {
82 setup(ubar);
83 }
84
85 void setup(const DiscreteFunctionType &ubar)
86 {
87 scheme_.fullOperator().jacobian(ubar, linearOperator_);
88 inverseOperator_ = std::make_shared<LinearInverseOperatorType>(linreduction_, linabstol_, maxIter_ );
90 sequence_ = scheme_.space().sequence();
91 ubar_.assign(ubar);
92 }
93
94 void constraint ( DiscreteFunctionType &u ) const { scheme_.constraint(u); }
95
97 {
98 assert( sequence_ == scheme_.space().sequence() );
100 linearOperator_(u,w);
101 w += affineShift_;
102 }
103
104 template <class GridFunction>
105 void operator() ( const GridFunction &arg, DiscreteFunctionType &dest ) const
106 {
107 assert( sequence_ == scheme_.space().sequence() );
108 DiscreteFunctionType tmp(dest);
109 Fem::interpolate(arg,tmp);
110 (*this)(tmp,dest);
111 }
112
114 {
115 int oldCount = (*inverseOperator_).iterations();
116 constraint(solution);
117 assert( sequence_ == scheme_.space().sequence() );
118 (*inverseOperator_)( rhs, solution );
119 return SolverInfo( true, (*inverseOperator_).iterations()-oldCount, 1 );
120 }
122 {
123 affineShift();
124 return solve(affineShift_,solution);
125 }
126
127 template< class GridFunction >
128 const LinearOperatorType &assemble ( const GridFunction &u )
129 {
130 assert( sequence_ == scheme_.space().sequence() );
131 return linearOperator_;
132 }
133
134 bool mark ( double tolerance ) { return scheme_.mark( tolerance ); }
135 double estimate ( const DiscreteFunctionType &solution ) { return scheme_.estimate( solution ); }
136
137 const GridPartType &gridPart () const { return scheme_.gridPart(); }
138 const DiscreteFunctionSpaceType &space() const { return scheme_.space(); }
139
140 protected:
141 void affineShift() const
142 {
144 tmp.clear();
145
146 affineShift_.clear();
148 scheme_.fullOperator()( ubar_, tmp );
149 affineShift_ -= tmp;
151 }
152
158 std::shared_ptr<LinearInverseOperatorType> inverseOperator_;
163 };
164
165 } // namespace Fem
166
167} // namespace Dune
168
169#endif // #ifndef DUNE_FEM_SCHEMES_LINEARIZED_HH
static void interpolate(const GridFunction &u, DiscreteFunction &v)
perform native interpolation of a discrete function space
Definition: common/interpolate.hh:54
STL namespace.
Definition: bindguard.hh:11
static ParameterContainer & container()
Definition: io/parameter.hh:193
Definition: linearized.hh:37
DiscreteFunctionType ubar_
Definition: linearized.hh:161
SchemeType::DiscreteFunctionType DiscreteFunctionType
Definition: linearized.hh:39
SchemeType & scheme_
Definition: linearized.hh:153
SolverInfo solve(const DiscreteFunctionType &rhs, DiscreteFunctionType &solution) const
Definition: linearized.hh:113
const LinearOperatorType & assemble(const GridFunction &u)
Definition: linearized.hh:128
SchemeType::ModelType ModelType
Definition: linearized.hh:44
bool mark(double tolerance)
Definition: linearized.hh:134
LinearizedScheme(SchemeType &scheme, const DiscreteFunctionType &ubar, Dune::Fem::ParameterReader parameter=Dune::Fem::Parameter::container())
Definition: linearized.hh:70
double linreduction_
Definition: linearized.hh:156
void setup(const DiscreteFunctionType &ubar)
Definition: linearized.hh:85
DiscreteFunctionType affineShift_
Definition: linearized.hh:159
const GridPartType & gridPart() const
Definition: linearized.hh:137
void operator()(const DiscreteFunctionType &u, DiscreteFunctionType &w) const
Definition: linearized.hh:96
Dune::Fem::ParameterReader parameter_
Definition: linearized.hh:160
SchemeType::GridPartType GridPartType
Definition: linearized.hh:41
Scheme SchemeType
Definition: linearized.hh:38
LinearOperatorType linearOperator_
Definition: linearized.hh:154
SchemeType::InverseOperatorType LinearInverseOperatorType
Definition: linearized.hh:43
int sequence_
Definition: linearized.hh:162
int maxIter_
Definition: linearized.hh:157
SchemeType::LinearOperatorType LinearOperatorType
Definition: linearized.hh:42
SolverInfo solve(DiscreteFunctionType &solution) const
Definition: linearized.hh:121
double estimate(const DiscreteFunctionType &solution)
Definition: linearized.hh:135
const DiscreteFunctionSpaceType & space() const
Definition: linearized.hh:138
SchemeType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType
Definition: linearized.hh:40
double linabstol_
Definition: linearized.hh:155
void affineShift() const
Definition: linearized.hh:141
LinearizedScheme(SchemeType &scheme, Dune::Fem::ParameterReader parameter=Dune::Fem::Parameter::container())
Definition: linearized.hh:55
std::shared_ptr< LinearInverseOperatorType > inverseOperator_
Definition: linearized.hh:158
void constraint(DiscreteFunctionType &u) const
Definition: linearized.hh:94
Definition: linearized.hh:46
int linearIterations
Definition: linearized.hh:52
int nonlinearIterations
Definition: linearized.hh:52
SolverInfo(bool converged, int linearIterations, int nonlinearIterations)
Definition: linearized.hh:47
bool converged
Definition: linearized.hh:51