dune-fem  2.6-git
basicrow.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_BASICROW_HH
2 #define DUNE_FEM_SOLVER_RUNGEKUTTA_BASICROW_HH
3 
4 //- system includes
5 #include <algorithm>
6 #include <cmath>
7 #include <cassert>
8 #include <limits>
9 #include <sstream>
10 #include <string>
11 #include <vector>
12 
13 //- dune-common includes
14 #include <dune/common/dynmatrix.hh>
15 #include <dune/common/dynvector.hh>
16 
17 //- dune-fem includes
20 
21 namespace DuneODE
22 {
23 
24  // NoROWRungeKuttaSourceTerm
25  // ------------------------------
26 
28  {
29  template< class T >
30  bool operator() ( double time, double timeStepSize, int stage, const T &u, const std::vector< T * > &update, T &source )
31  {
32  return false;
33  }
34 
35  template< class T >
36  void limit( T& update, const double time ) {}
37 
38  template< class T >
39  double initialTimeStepEstimate ( double time, const T &u ) const
40  {
41  // return negative value to indicate that implicit time step should be used
42  return -1.0;
43  }
44 
45  double timeStepEstimate () const
46  {
48  }
49 
50  };
51 
53  : public Dune::Fem::LocalParameter< ROWSolverParameter, ROWSolverParameter >
54  {
55  protected:
56  // key prefix, default is fem.solver.row. (can be overloaded by user)
57  const std::string keyPrefix_;
58 
60 
61  public:
63  : keyPrefix_( keyPrefix ),
65  {}
66 
68  : keyPrefix_( "fem.solver.row." ),
70  {}
71 
72  const Dune::Fem::ParameterReader &parameter () const { return parameter_; }
73 
74  virtual double linAbsTolParameter ( ) const
75  {
76  return parameter().getValue< double >( keyPrefix_ + "linabstol", 1e-6 );
77  }
78 
79  virtual double linReductionParameter ( ) const
80  {
81  return parameter().getValue< double >( keyPrefix_ + "linreduction", 1e-4 );
82  }
83 
84  virtual bool linearSolverVerbose () const
85  {
86  return parameter().getValue< bool >( keyPrefix_ + "verbose", false );
87  }
88 
89  virtual int maxLinearIterationsParameter () const
90  {
91  return parameter().getValue< int >( keyPrefix_ + "maxlineariterations", std::numeric_limits< int >::max() );
92  }
93  };
94 
95 
96 
98  template< class HelmholtzOperator, class NonlinearSolver, class TimeStepControl, class SourceTerm = NoROWRungeKuttaSourceTerm >
100  : public OdeSolverInterface< typename HelmholtzOperator::DomainFunctionType >
101  {
104 
105  public:
108 
109  typedef HelmholtzOperator HelmholtzOperatorType;
110  typedef NonlinearSolver NonlinearSolverType;
111  typedef TimeStepControl TimeStepControlType;
112  typedef SourceTerm SourceTermType;
113 
115  typedef typename NonlinearSolver::ParametersType NonlinearSolverParametersType;
116 
117  typedef typename HelmholtzOperator::SpaceOperatorType::PreconditionOperatorType PreconditionOperatorType;
118 
120 
128  template< class ButcherTable >
130  TimeProviderType& timeProvider,
131  const ButcherTable &butcherTable,
132  const TimeStepControlType &timeStepControl,
133  const SourceTermType &sourceTerm,
134  const ParametersType& parameter,
135  const NonlinearSolverParametersType& nlsParam )
136  : helmholtzOp_( helmholtzOp ),
137  nonlinearSolver_( helmholtzOp_, nlsParam ),
138  timeStepControl_( timeStepControl ),
139  sourceTerm_( sourceTerm ),
140  stages_( butcherTable.stages() ),
141  alpha_( butcherTable.A() ),
142  alpha2_( butcherTable.B() ),
143  gamma_( stages() ),
144  beta_( stages() ),
145  c_( butcherTable.c() ),
146  rhs_( "RK rhs", helmholtzOp_.space() ),
147  temp_( "RK temp", helmholtzOp_.space() ),
148  update_( stages(), nullptr ),
149  linAbsTol_( parameter.linAbsTolParameter( ) ),
150  linReduction_( parameter.linReductionParameter( ) ),
151  linVerbose_( parameter.linearSolverVerbose() ),
152  maxLinearIterations_( parameter.maxLinearIterationsParameter() ),
153  preconditioner_(helmholtzOp.spaceOperator().preconditioner())
154  {
155  setup( butcherTable );
156  }
157 
158  template< class ButcherTable >
160  TimeProviderType& timeProvider,
161  const ButcherTable &butcherTable,
162  const TimeStepControlType &timeStepControl,
163  const SourceTermType &sourceTerm,
165  : BasicROWRungeKuttaSolver( helmholtzOp, timeProvider, butcherTable, timeStepControl, sourceTerm, ParametersType( parameter ),
166  NonlinearSolverParametersType( parameter ) )
167  {}
168 
175  template< class ButcherTable >
177  TimeProviderType& timeProvider,
178  const ButcherTable &butcherTable,
179  const TimeStepControlType &timeStepControl,
180  const ParametersType& parameter,
181  const NonlinearSolverParametersType& nlsParam )
182  : BasicROWRungeKuttaSolver( helmholtzOp, timeProvider, butcherTable, timeStepControl, SourceTermType(), parameter, nlsParam )
183  {}
184 
185  template< class ButcherTable >
187  TimeProviderType& timeProvider,
188  const ButcherTable &butcherTable,
189  const TimeStepControlType &timeStepControl,
191  : BasicROWRungeKuttaSolver( helmholtzOp, timeProvider, butcherTable, timeStepControl, SourceTermType(), ParametersType( parameter ),
192  NonlinearSolverParametersType( parameter ) )
193  {}
194 
200  template< class ButcherTable >
202  TimeProviderType& timeProvider,
203  const ButcherTable &butcherTable,
204  const ParametersType& parameter,
205  const NonlinearSolverParametersType& nlsParam )
206  : BasicROWRungeKuttaSolver( helmholtzOp, timeProvider, butcherTable, TimeStepControlType(), SourceTermType(), parameter, nlsParam )
207  {}
208 
209  template< class ButcherTable >
211  TimeProviderType& timeProvider,
212  const ButcherTable &butcherTable,
214  : BasicROWRungeKuttaSolver( helmholtzOp, timeProvider, butcherTable, TimeStepControlType(), SourceTermType(), ParametersType( parameter ),
215  NonlinearSolverParametersType( parameter ) )
216  {}
217 
218  template< class ButcherTable >
219  void setup( const ButcherTable& butcherTable )
220  {
222  {
223  std::cout << "ROW method of order=" << butcherTable.order() << " with " << stages_ << " stages" << std::endl;
224  }
225 
226  // create intermediate functions
227  for( int i = 0; i < stages(); ++i )
228  {
229  std::ostringstream name;
230  name << "RK stage " << i;
231  update_[ i ] = new DestinationType( name.str(), helmholtzOp_.space() );
232  }
233 
234  // compute coefficients
235  for( int i = 0; i < stages(); ++i )
236  gamma_[ i ] = alpha_[ i ][ i ];
237 
238  alpha_.invert();
239  alpha_.mtv( butcherTable.b(), beta_ );
240  alpha2_.rightmultiply( alpha_ );
241  }
242 
245  {
246  for( int i = 0; i < stages(); ++i )
247  delete update_[ i ];
248  }
249 
251  void initialize ( const DestinationType &U0 )
252  {
253  const double time = timeStepControl_.time();
254 
255  helmholtzOp_.setTime( time );
256  helmholtzOp_.initializeTimeStepSize( U0 );
257  const double helmholtzEstimate = helmholtzOp_.timeStepEstimate();
258 
259  double sourceTermEstimate = sourceTerm_.initialTimeStepEstimate( time, U0 );
260  // negative time step is given by the empty source term
261  if( sourceTermEstimate < 0.0 ) sourceTermEstimate = helmholtzEstimate ;
262 
263  timeStepControl_.initialTimeStepSize( helmholtzEstimate, sourceTermEstimate );
264  }
265 
266  using BaseType::solve;
267 
269  void solve ( DestinationType &U, MonitorType &monitor )
270  {
271  monitor.reset();
272 
273  typename HelmholtzOperatorType::JacobianOperatorType jOp( "jacobianOperator", U.space(), U.space() );
274 
275  const double time = timeStepControl_.time();
276  const double timeStepSize = timeStepControl_.timeStepSize();
277  assert( timeStepSize > 0.0 );
278  // std::cout << "solving... " << time << " : " << U << std::endl;
279  for( int s = 0; s < stages(); ++s )
280  {
281  // update for stage s
282  DestinationType& updateStage = *update_[ s ];
283 
284  rhs_.assign( U );
285  for( int k = 0; k < s; ++k )
286  rhs_.axpy( alpha2_[ s ][ k ], *update_[ k ] );
287  helmholtzOp_.spaceOperator()(rhs_,updateStage);
288  updateStage *= (gamma_[s]*timeStepSize);
289  for( int k = 0; k < s; ++k )
290  updateStage.axpy( -gamma_[s]*alpha_[ s ][ k ], *update_[ k ] );
291 
292  rhs_.assign( updateStage );
293 
294  // solve the system
295  const double stageTime = time + c_[ s ]*timeStepSize;
296  helmholtzOp_.setTime( stageTime );
297  // compute jacobian if the diagonal entry in the butcher tableau changes
298  // if ( s==0 || (gamma_[s-1] != gamma_[s]) )
299  {
300  // std::cout << "lambda=" << gamma_[ s ]*timeStepSize << std::endl;
301  helmholtzOp_.setLambda( gamma_[ s ]*timeStepSize );
302  helmholtzOp_.jacobian( U, jOp );
303  }
304  const int remLinearIts = maxLinearIterations_;
305  typename NonlinearSolverType::LinearInverseOperatorType jInv( linReduction_, linAbsTol_, remLinearIts, linVerbose_ );
306  if (preconditioner_)
307  {
308  jInv.bind( jOp, *preconditioner_ );
309  //typename NonlinearSolverType::LinearInverseOperatorType jInv( jOp, *preconditioner_, linReduction_, linAbsTol_, remLinearIts, linVerbose_ );
310  jInv( rhs_, updateStage );
311  monitor.linearSolverIterations_ += jInv.iterations();
312  }
313  else
314  {
315  jInv.bind( jOp );
316  jInv( rhs_, updateStage );
317  monitor.linearSolverIterations_ += jInv.iterations();
318  }
319 
320  jInv.unbind();
321  }
322 
323  double error = 0.0;
324  if(0 && timeStepControl_.computeError() )
325  {
326  // store U (to be revised)
327  DestinationType Uerr( U );
328 
329  // update solution
330  U *= delta_;
331  for( int s = 0; s < stages(); ++s )
332  U.axpy( beta_[ s ], *update_[ s ] );
333 
334  //error = infNorm( U, Uerr );
335  Uerr.axpy( -1.0, U );
336  const double errorU = Uerr.scalarProductDofs( Uerr );
337  const double normU = U.scalarProductDofs( U );
338 
339  if( normU > 0 && errorU > 0 )
340  {
341  error = std::sqrt( errorU / normU );
342  }
343  std::cout << std::scientific << "Error in RK = " << error << " norm " << errorU << " " << normU << std::endl;
344  //std::cout << std::scientific << "Error in RK = " << error << std::endl;
345  }
346  else
347  {
348  // update solution
349  for( int s = 0; s < stages(); ++s )
350  U.axpy( beta_[ s ], *update_[ s ] );
351  }
352  // set error to monitor
353  monitor.error_ = error;
354 
355  // update time step size
356  timeStepControl_.timeStepEstimate( helmholtzOp_.timeStepEstimate(), sourceTerm_.timeStepEstimate(), monitor );
357  }
358 
359  int stages () const { return stages_; }
360 
361  void description ( std::ostream &out ) const
362  {
363  out << "Generic " << stages() << "-stage implicit Runge-Kutta solver.\\\\" << std::endl;
364  }
365 
366  protected:
367 
368  double infNorm(const DestinationType& U, const DestinationType& Uerr ) const
369  {
370  typedef typename DestinationType :: ConstDofIteratorType ConstDofIteratorType ;
371  const ConstDofIteratorType uend = U.dend();
372  double res = 0;
373  for( ConstDofIteratorType u = U.dbegin(), uerr = Uerr.dbegin(); u != uend; ++u, ++uerr )
374  {
375  double uval = *u;
376  double uerrval = *uerr ;
377  double div = std::abs( std::max( uval, uerrval ) );
378 
379  double norm = std::abs( uval - uerrval );
380  if( std::abs(div) > 1e-12 )
381  norm /= div;
382  res = std::max( res, norm );
383  }
384  return res;
385  }
386 
389  TimeStepControl timeStepControl_;
390  SourceTerm sourceTerm_;
391 
392  int stages_;
393  double delta_;
394  Dune::DynamicMatrix< double > alpha_, alpha2_;
395  Dune::DynamicVector< double > gamma_, beta_, c_;
396 
398  std::vector< DestinationType * > update_;
399 
400  const double linAbsTol_, linReduction_;
401  const bool linVerbose_;
403 
405  };
406 
407 } // namespace DuneODE
408 
409 #endif // #ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_BASICROW_HH
Double abs(const Double &a)
Definition: double.hh:872
static double sqrt(const Double &v)
Definition: double.hh:887
static constexpr T max(T a)
Definition: utility.hh:77
Definition: multistep.hh:17
T getValue(const std::string &key) const
get mandatory parameter
Definition: reader.hh:149
static bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:437
static DUNE_EXPORT ParameterContainer & container()
Definition: io/parameter.hh:192
Definition: io/parameter.hh:544
Interface class for ODE Solver.
Definition: odesolverinterface.hh:21
Monitor MonitorType
monitor type
Definition: odesolverinterface.hh:59
virtual void solve(DestinationType &u)
solve where is the internal operator.
Definition: odesolverinterface.hh:75
DestinationImp DestinationType
type of destination
Definition: odesolverinterface.hh:62
Definition: odesolverinterface.hh:27
Definition: basicrow.hh:28
double timeStepEstimate() const
Definition: basicrow.hh:45
void limit(T &update, const double time)
Definition: basicrow.hh:36
double initialTimeStepEstimate(double time, const T &u) const
Definition: basicrow.hh:39
bool operator()(double time, double timeStepSize, int stage, const T &u, const std::vector< T * > &update, T &source)
Definition: basicrow.hh:30
Definition: basicrow.hh:54
virtual int maxLinearIterationsParameter() const
Definition: basicrow.hh:89
const Dune::Fem::ParameterReader & parameter() const
Definition: basicrow.hh:72
const std::string keyPrefix_
Definition: basicrow.hh:57
ROWSolverParameter(const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: basicrow.hh:67
virtual double linAbsTolParameter() const
Definition: basicrow.hh:74
Dune::Fem::ParameterReader parameter_
Definition: basicrow.hh:59
virtual bool linearSolverVerbose() const
Definition: basicrow.hh:84
ROWSolverParameter(const std::string keyPrefix, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: basicrow.hh:62
virtual double linReductionParameter() const
Definition: basicrow.hh:79
ROW RungeKutta ODE solver.
Definition: basicrow.hh:101
void solve(DestinationType &U, MonitorType &monitor)
solve the system
Definition: basicrow.hh:269
Dune::DynamicVector< double > beta_
Definition: basicrow.hh:395
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const ParametersType &parameter, const NonlinearSolverParametersType &nlsParam)
constructor
Definition: basicrow.hh:176
void description(std::ostream &out) const
print description of ODE solver to out stream
Definition: basicrow.hh:361
SourceTerm SourceTermType
Definition: basicrow.hh:112
int stages_
Definition: basicrow.hh:392
TimeStepControl TimeStepControlType
Definition: basicrow.hh:111
SourceTerm sourceTerm_
Definition: basicrow.hh:390
std::vector< DestinationType * > update_
Definition: basicrow.hh:398
NonlinearSolver::ParametersType NonlinearSolverParametersType
Definition: basicrow.hh:115
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const SourceTermType &sourceTerm, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: basicrow.hh:159
Dune::DynamicMatrix< double > alpha2_
Definition: basicrow.hh:394
HelmholtzOperatorType & helmholtzOp_
Definition: basicrow.hh:387
Dune::DynamicMatrix< double > alpha_
Definition: basicrow.hh:394
const PreconditionOperatorType * preconditioner_
Definition: basicrow.hh:404
BaseType::MonitorType MonitorType
Definition: basicrow.hh:106
BaseType::DestinationType DestinationType
Definition: basicrow.hh:107
NonlinearSolverType nonlinearSolver_
Definition: basicrow.hh:388
int stages() const
Definition: basicrow.hh:359
ROWSolverParameter ParametersType
Definition: basicrow.hh:114
const double linAbsTol_
Definition: basicrow.hh:400
const double linReduction_
Definition: basicrow.hh:400
TimeStepControl timeStepControl_
Definition: basicrow.hh:389
DestinationType temp_
Definition: basicrow.hh:397
NonlinearSolver NonlinearSolverType
Definition: basicrow.hh:110
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const SourceTermType &sourceTerm, const ParametersType &parameter, const NonlinearSolverParametersType &nlsParam)
constructor
Definition: basicrow.hh:129
double delta_
Definition: basicrow.hh:393
const bool linVerbose_
Definition: basicrow.hh:401
const int maxLinearIterations_
Definition: basicrow.hh:402
Dune::Fem::TimeProviderBase TimeProviderType
Definition: basicrow.hh:119
void initialize(const DestinationType &U0)
apply operator once to get dt estimate
Definition: basicrow.hh:251
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const ParametersType &parameter, const NonlinearSolverParametersType &nlsParam)
constructor
Definition: basicrow.hh:201
HelmholtzOperator::SpaceOperatorType::PreconditionOperatorType PreconditionOperatorType
Definition: basicrow.hh:117
HelmholtzOperator HelmholtzOperatorType
Definition: basicrow.hh:109
double infNorm(const DestinationType &U, const DestinationType &Uerr) const
Definition: basicrow.hh:368
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: basicrow.hh:210
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: basicrow.hh:186
Dune::DynamicVector< double > gamma_
Definition: basicrow.hh:395
DestinationType rhs_
Definition: basicrow.hh:397
Dune::DynamicVector< double > c_
Definition: basicrow.hh:395
~BasicROWRungeKuttaSolver()
destructor
Definition: basicrow.hh:244
void setup(const ButcherTable &butcherTable)
Definition: basicrow.hh:219
Definition: timestepcontrol.hh:166
general base for time providers
Definition: timeprovider.hh:36