dune-fem  2.6-git
automaticdifferenceoperator.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_AUTOMATICDIFFERENCEOPERATOR_HH
2 #define DUNE_FEM_AUTOMATICDIFFERENCEOPERATOR_HH
3 
4 #include <limits>
5 
8 
9 namespace Dune
10 {
11 
12  namespace Fem
13  {
14 
15  // Internal Forward Declarations
16  // -----------------------------
17 
18  template< class DomainFunction, class RangeFunction, class LinearOperator >
19  class AutomaticDifferenceOperator;
20 
21 
22 
23  // AutomaticDifferenceLinearOperator
24  // ---------------------------------
25 
26  template< class DomainFunction, class RangeFunction = DomainFunction >
28  : public Dune::Fem::Operator< DomainFunction, RangeFunction >
29  {
32 
33  friend class AutomaticDifferenceOperator< DomainFunction, RangeFunction, ThisType >;
34 
35  public:
37 
42  typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
43 
44  typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType;
45  typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType;
46 
47  AutomaticDifferenceLinearOperator ( const std::string &name, const DomainSpaceType &dSpace, const RangeSpaceType &rSpace )
48  : name_( name ),
49  op_( 0 ), // initial value for op_ is 'undefined'
50  u_( 0 ), // initial value for u_ is 'undefined'
51  b_( "AutomaticDifferenceOperator::b_", dSpace ),
52  op_u_( "AutomaticDifferenceOperator::op_u_", rSpace ),
53  norm_u_( 0 )
54  {}
55 
56  virtual void operator() ( const DomainFunctionType &arg, RangeFunctionType &dest ) const;
57 
58  protected:
59  void set ( const DomainFunctionType &u, const OperatorType &op, const RealType &eps );
60  const std::string name_;
61  const OperatorType *op_;
63 
66 
69  };
70 
71 
72 
73  // AutomaticDifferenceOperator
74  // ---------------------------
75 
82  template< class DomainFunction, class RangeFunction = DomainFunction,
85  : public Dune::Fem::DifferentiableOperator< LinearOperator >
86  {
88 
89  public:
94  typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
95 
97 
98  typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType;
99  typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType;
100 
102  : eps_( parameter.getValue< RangeFieldType >( "fem.differenceoperator.eps", RangeFieldType( 0 ) ) )
103  {}
104 
106  : eps_( eps )
107  {}
108 
109  virtual void jacobian ( const DomainFunctionType &u, JacobianOperatorType &jOp ) const
110  {
111  jOp.set( u, *this, eps_ );
112  }
113 
114  private:
115  const RealType eps_;
116  };
117 
118 
119 
120  // Implementation of AutomaticDifferenceLinearOperator
121  // ---------------------------------------------------
122 
123  template< class DomainFunction, class RangeFunction >
126  {
127  assert( op_ && u_ );
128 
129  // 'Normal' difference-quotient, i.e.
130  // dest = 1/eps (op_(*u +eps*arg) - op_(*u))
131  //
132  // eps is chosen dynamically, the same way as in
133  // dune-fem/dune/fem/solver/ode/quasi_exact_newton.cpp, see also
134  // http://www.freidok.uni-freiburg.de/volltexte/3762/pdf/diss.pdf, page 137
135  RealType eps = eps_;
136  if( eps <= RealType( 0. ) )
137  {
138  const RealType machine_eps = std::numeric_limits< RealType >::epsilon();
139  const RealType norm_p_sq = arg.normSquaredDofs( );
140  if( norm_p_sq > machine_eps )
141  eps = std::sqrt( (RealType( 1 ) + norm_u_) * machine_eps / norm_p_sq );
142  else
143  eps = std::sqrt( machine_eps );
144  }
145 
146  b_.assign( *u_ );
147  b_.axpy( eps, arg );
148  (*op_)( b_, dest );
149  dest -= op_u_;
150  dest *= RealType( 1 ) / eps;
151  }
152 
153 
154  template< class DomainFunction, class RangeFunction >
156  ::set ( const DomainFunctionType &u, const OperatorType &op, const RealType &eps )
157  {
158  u_ = &u;
159  op_ = &op;
160  (*op_)( *u_, op_u_ );
161 
162  eps_ = eps;
163  if( eps_ <= RealType( 0 ) )
164  norm_u_ = std::sqrt( u_->scalarProductDofs( *u_ ) );
165  }
166 
167  } // namespace Fem
168 
169 } // namespace Dune
170 
171 #endif // #ifndef DUNE_FEM_AUTOMATICDIFFERENCEOPERATOR_HH
Definition: bindguard.hh:11
static double sqrt(const Double &v)
Definition: double.hh:887
static DUNE_EXPORT ParameterContainer & container()
Definition: io/parameter.hh:192
operator providing a Jacobian through automatic differentiation
Definition: automaticdifferenceoperator.hh:86
Dune::FieldTraits< RangeFieldType >::real_type RealType
Definition: automaticdifferenceoperator.hh:94
BaseType::JacobianOperatorType JacobianOperatorType
Definition: automaticdifferenceoperator.hh:96
BaseType::DomainFieldType DomainFieldType
Definition: automaticdifferenceoperator.hh:93
BaseType::RangeFunctionType RangeFunctionType
Definition: automaticdifferenceoperator.hh:90
AutomaticDifferenceOperator(const RangeFieldType &eps)
Definition: automaticdifferenceoperator.hh:105
RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType
Definition: automaticdifferenceoperator.hh:98
DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType
Definition: automaticdifferenceoperator.hh:99
virtual void jacobian(const DomainFunctionType &u, JacobianOperatorType &jOp) const
Definition: automaticdifferenceoperator.hh:109
BaseType::DomainFunctionType DomainFunctionType
Definition: automaticdifferenceoperator.hh:91
AutomaticDifferenceOperator(const ParameterReader &parameter=Parameter::container())
Definition: automaticdifferenceoperator.hh:101
BaseType::RangeFieldType RangeFieldType
Definition: automaticdifferenceoperator.hh:92
Definition: automaticdifferenceoperator.hh:29
RangeFieldType norm_u_
Definition: automaticdifferenceoperator.hh:68
DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType
Definition: automaticdifferenceoperator.hh:45
Dune::Fem::Operator< DomainFunction, RangeFunction > OperatorType
Definition: automaticdifferenceoperator.hh:36
BaseType::RangeFunctionType RangeFunctionType
Definition: automaticdifferenceoperator.hh:38
RangeFieldType eps_
Definition: automaticdifferenceoperator.hh:67
Dune::FieldTraits< RangeFieldType >::real_type RealType
Definition: automaticdifferenceoperator.hh:42
const OperatorType * op_
Definition: automaticdifferenceoperator.hh:61
void set(const DomainFunctionType &u, const OperatorType &op, const RealType &eps)
Definition: automaticdifferenceoperator.hh:156
const DomainFunctionType * u_
Definition: automaticdifferenceoperator.hh:62
virtual void operator()(const DomainFunctionType &arg, RangeFunctionType &dest) const
Definition: automaticdifferenceoperator.hh:125
BaseType::DomainFieldType DomainFieldType
Definition: automaticdifferenceoperator.hh:41
AutomaticDifferenceLinearOperator(const std::string &name, const DomainSpaceType &dSpace, const RangeSpaceType &rSpace)
Definition: automaticdifferenceoperator.hh:47
DomainFunctionType b_
Definition: automaticdifferenceoperator.hh:64
RangeFunctionType op_u_
Definition: automaticdifferenceoperator.hh:65
const std::string name_
Definition: automaticdifferenceoperator.hh:60
BaseType::DomainFunctionType DomainFunctionType
Definition: automaticdifferenceoperator.hh:39
BaseType::RangeFieldType RangeFieldType
Definition: automaticdifferenceoperator.hh:40
RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType
Definition: automaticdifferenceoperator.hh:44
abstract differentiable operator
Definition: differentiableoperator.hh:29
BaseType::RangeFunctionType RangeFunctionType
type of discrete function in the operator's range
Definition: differentiableoperator.hh:40
BaseType::DomainFunctionType DomainFunctionType
type of discrete function in the operator's domain
Definition: differentiableoperator.hh:38
DomainFunctionType::RangeFieldType DomainFieldType
field type of the operator's domain
Definition: differentiableoperator.hh:43
RangeFunctionType::RangeFieldType RangeFieldType
field type of the operator's range
Definition: differentiableoperator.hh:45
abstract operator
Definition: operator.hh:26
DomainFunction DomainFunctionType
type of discrete function in the operator's domain
Definition: operator.hh:28
RangeFunction::RangeFieldType RangeFieldType
field type of the operator's range
Definition: operator.hh:35
DomainFunction::RangeFieldType DomainFieldType
field type of the operator's domain
Definition: operator.hh:33
RangeFunction RangeFunctionType
type of discrete function in the operator's range
Definition: operator.hh:30
abstract affine-linear operator
Definition: operator.hh:73