dune-fem  2.8-git
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 
12 #include <dune/fem/io/parameter.hh>
13 
14 
15 namespace 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;
45  struct SolverInfo
46  {
49  {}
50 
51  bool converged;
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() );
99  affineShift();
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 
155  double linabstol_;
157  int maxIter_;
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
Definition: bindguard.hh:11
static ParameterContainer & container()
Definition: io/parameter.hh:193
Definition: linearized.hh:37
const GridPartType & gridPart() const
Definition: linearized.hh:137
DiscreteFunctionType ubar_
Definition: linearized.hh:161
SchemeType::DiscreteFunctionType DiscreteFunctionType
Definition: linearized.hh:39
SchemeType & scheme_
Definition: linearized.hh:153
const DiscreteFunctionSpaceType & space() const
Definition: linearized.hh:138
SolverInfo solve(const DiscreteFunctionType &rhs, DiscreteFunctionType &solution) const
Definition: linearized.hh:113
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
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
const LinearOperatorType & assemble(const GridFunction &u)
Definition: linearized.hh:128
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
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