dune-fem  2.8-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  template< class HelmholtzOperator, class NonlinearSolver, class TimeStepControl, class SourceTerm = NoROWRungeKuttaSourceTerm >
55  : public OdeSolverInterface< typename HelmholtzOperator::DomainFunctionType >
56  {
59 
60  public:
63 
64  typedef HelmholtzOperator HelmholtzOperatorType;
65  typedef NonlinearSolver NonlinearSolverType;
66  typedef TimeStepControl TimeStepControlType;
67  typedef SourceTerm SourceTermType;
68 
69  typedef typename NonlinearSolver::ParameterType NonlinearSolverParameterType;
70  typedef typename NonlinearSolverType::LinearInverseOperatorType LinearInverseOperatorType;
71 
72  typedef typename HelmholtzOperator::SpaceOperatorType::PreconditionOperatorType PreconditionOperatorType;
73 
75 
83  template< class ButcherTable >
85  TimeProviderType& timeProvider,
86  const ButcherTable &butcherTable,
87  const TimeStepControlType &timeStepControl,
88  const SourceTermType &sourceTerm,
89  const NonlinearSolverParameterType& parameter )
90  : helmholtzOp_( helmholtzOp ),
91  linearSolver_( parameter.linear() ),
92  timeStepControl_( timeStepControl ),
93  sourceTerm_( sourceTerm ),
94  stages_( butcherTable.stages() ),
95  alpha_( butcherTable.A() ),
96  alpha2_( butcherTable.B() ),
97  gamma_( stages() ),
98  beta_( stages() ),
99  c_( butcherTable.c() ),
100  rhs_( "RK rhs", helmholtzOp_.space() ),
101  temp_( "RK temp", helmholtzOp_.space() ),
102  update_( stages(), nullptr ),
103  maxLinearIterations_( parameter.linear().maxIterations() ),
104  preconditioner_(helmholtzOp.spaceOperator().preconditioner())
105  {
106  setup( butcherTable );
107  }
108 
109  template< class ButcherTable >
111  TimeProviderType& timeProvider,
112  const ButcherTable &butcherTable,
113  const TimeStepControlType &timeStepControl,
114  const SourceTermType &sourceTerm,
116  : BasicROWRungeKuttaSolver( helmholtzOp, timeProvider, butcherTable, timeStepControl, sourceTerm, NonlinearSolverParameterType( parameter ) )
117  {}
118 
125  template< class ButcherTable >
127  TimeProviderType& timeProvider,
128  const ButcherTable &butcherTable,
129  const TimeStepControlType &timeStepControl,
130  const NonlinearSolverParameterType& parameter )
131  : BasicROWRungeKuttaSolver( helmholtzOp, timeProvider, butcherTable, timeStepControl, SourceTermType(), parameter )
132  {}
133 
134  template< class ButcherTable >
136  TimeProviderType& timeProvider,
137  const ButcherTable &butcherTable,
138  const TimeStepControlType &timeStepControl,
140  : BasicROWRungeKuttaSolver( helmholtzOp, timeProvider, butcherTable, timeStepControl, SourceTermType(), NonlinearSolverParameterType( parameter ) )
141  {}
142 
148  template< class ButcherTable >
150  TimeProviderType& timeProvider,
151  const ButcherTable &butcherTable,
152  const NonlinearSolverParameterType& parameter )
153  : BasicROWRungeKuttaSolver( helmholtzOp, timeProvider, butcherTable, TimeStepControlType(), SourceTermType(), parameter )
154  {}
155 
156  template< class ButcherTable >
158  TimeProviderType& timeProvider,
159  const ButcherTable &butcherTable,
161  : BasicROWRungeKuttaSolver( helmholtzOp, timeProvider, butcherTable, TimeStepControlType(), SourceTermType(), NonlinearSolverParameterType( parameter ) )
162  {}
163 
164  template< class ButcherTable >
165  void setup( const ButcherTable& butcherTable )
166  {
168  {
169  std::cout << "ROW method of order=" << butcherTable.order() << " with " << stages_ << " stages" << std::endl;
170  }
171 
172  // create intermediate functions
173  for( int i = 0; i < stages(); ++i )
174  {
175  std::ostringstream name;
176  name << "RK stage " << i;
177  update_[ i ] = new DestinationType( name.str(), helmholtzOp_.space() );
178  }
179 
180  // compute coefficients
181  for( int i = 0; i < stages(); ++i )
182  gamma_[ i ] = alpha_[ i ][ i ];
183 
184  alpha_.invert();
185  alpha_.mtv( butcherTable.b(), beta_ );
186  alpha2_.rightmultiply( alpha_ );
187  }
188 
191  {
192  for( int i = 0; i < stages(); ++i )
193  delete update_[ i ];
194  }
195 
197  void initialize ( const DestinationType &U0 )
198  {
199  const double time = timeStepControl_.time();
200 
201  helmholtzOp_.setTime( time );
202  helmholtzOp_.initializeTimeStepSize( U0 );
203  const double helmholtzEstimate = helmholtzOp_.timeStepEstimate();
204 
205  double sourceTermEstimate = sourceTerm_.initialTimeStepEstimate( time, U0 );
206  // negative time step is given by the empty source term
207  if( sourceTermEstimate < 0.0 ) sourceTermEstimate = helmholtzEstimate ;
208 
209  timeStepControl_.initialTimeStepSize( helmholtzEstimate, sourceTermEstimate );
210  }
211 
212  using BaseType::solve;
213 
215  void solve ( DestinationType &U, MonitorType &monitor )
216  {
217  monitor.reset();
218 
219  typename HelmholtzOperatorType::JacobianOperatorType jOp( "jacobianOperator", U.space(), U.space() );
220 
221  const double time = timeStepControl_.time();
222  const double timeStepSize = timeStepControl_.timeStepSize();
223  assert( timeStepSize > 0.0 );
224  // std::cout << "solving... " << time << " : " << U << std::endl;
225  for( int s = 0; s < stages(); ++s )
226  {
227  // update for stage s
228  DestinationType& updateStage = *update_[ s ];
229 
230  rhs_.assign( U );
231  for( int k = 0; k < s; ++k )
232  rhs_.axpy( alpha2_[ s ][ k ], *update_[ k ] );
233  helmholtzOp_.spaceOperator()(rhs_,updateStage);
234  updateStage *= (gamma_[s]*timeStepSize);
235  for( int k = 0; k < s; ++k )
236  updateStage.axpy( -gamma_[s]*alpha_[ s ][ k ], *update_[ k ] );
237 
238  rhs_.assign( updateStage );
239 
240  // solve the system
241  const double stageTime = time + c_[ s ]*timeStepSize;
242  helmholtzOp_.setTime( stageTime );
243  // compute jacobian if the diagonal entry in the butcher tableau changes
244  // if ( s==0 || (gamma_[s-1] != gamma_[s]) )
245  {
246  // std::cout << "lambda=" << gamma_[ s ]*timeStepSize << std::endl;
247  helmholtzOp_.setLambda( gamma_[ s ]*timeStepSize );
248  helmholtzOp_.jacobian( U, jOp );
249  }
250  const int remLinearIts = maxLinearIterations_;
251 
252  linearSolver_.parameter().setMaxIterations( remLinearIts );
253 
254  if (preconditioner_)
255  linearSolver_.bind( jOp, *preconditioner_ );
256  else
257  linearSolver_.bind( jOp );
258 
259  linearSolver_( rhs_, updateStage );
260  monitor.linearSolverIterations_ += linearSolver_.iterations();
261 
262  linearSolver_.unbind();
263  }
264 
265  double error = 0.0;
266  if(0 && timeStepControl_.computeError() )
267  {
268  // store U (to be revised)
269  DestinationType Uerr( U );
270 
271  // update solution
272  U *= delta_;
273  for( int s = 0; s < stages(); ++s )
274  U.axpy( beta_[ s ], *update_[ s ] );
275 
276  //error = infNorm( U, Uerr );
277  Uerr.axpy( -1.0, U );
278  const double errorU = Uerr.scalarProductDofs( Uerr );
279  const double normU = U.scalarProductDofs( U );
280 
281  if( normU > 0 && errorU > 0 )
282  {
283  error = std::sqrt( errorU / normU );
284  }
285  std::cout << std::scientific << "Error in RK = " << error << " norm " << errorU << " " << normU << std::endl;
286  //std::cout << std::scientific << "Error in RK = " << error << std::endl;
287  }
288  else
289  {
290  // update solution
291  for( int s = 0; s < stages(); ++s )
292  U.axpy( beta_[ s ], *update_[ s ] );
293  }
294  // set error to monitor
295  monitor.error_ = error;
296 
297  // update time step size
298  timeStepControl_.timeStepEstimate( helmholtzOp_.timeStepEstimate(), sourceTerm_.timeStepEstimate(), monitor );
299  }
300 
301  int stages () const { return stages_; }
302 
303  void description ( std::ostream &out ) const
304  {
305  out << "Generic " << stages() << "-stage implicit Runge-Kutta solver.\\\\" << std::endl;
306  }
307 
308  protected:
309 
310  double infNorm(const DestinationType& U, const DestinationType& Uerr ) const
311  {
312  typedef typename DestinationType :: ConstDofIteratorType ConstDofIteratorType ;
313  const ConstDofIteratorType uend = U.dend();
314  double res = 0;
315  for( ConstDofIteratorType u = U.dbegin(), uerr = Uerr.dbegin(); u != uend; ++u, ++uerr )
316  {
317  double uval = *u;
318  double uerrval = *uerr ;
319  double div = std::abs( std::max( uval, uerrval ) );
320 
321  double norm = std::abs( uval - uerrval );
322  if( std::abs(div) > 1e-12 )
323  norm /= div;
324  res = std::max( res, norm );
325  }
326  return res;
327  }
328 
331 
332  TimeStepControl timeStepControl_;
333  SourceTerm sourceTerm_;
334 
335  int stages_;
336  double delta_;
337  Dune::DynamicMatrix< double > alpha_, alpha2_;
338  Dune::DynamicVector< double > gamma_, beta_, c_;
339 
341  std::vector< DestinationType * > update_;
342 
344 
346  };
347 
348 } // namespace DuneODE
349 
350 #endif // #ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_BASICROW_HH
Double abs(const Double &a)
Definition: double.hh:871
static double sqrt(const Double &v)
Definition: double.hh:886
static constexpr T max(T a)
Definition: utility.hh:77
Definition: multistep.hh:17
static bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:445
static ParameterContainer & container()
Definition: io/parameter.hh:193
Interface class for ODE Solver.
Definition: odesolverinterface.hh:21
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
double error_
Definition: odesolverinterface.hh:30
void reset()
Definition: odesolverinterface.hh:43
int linearSolverIterations_
Definition: odesolverinterface.hh:35
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
ROW RungeKutta ODE solver.
Definition: basicrow.hh:56
void solve(DestinationType &U, MonitorType &monitor)
solve the system
Definition: basicrow.hh:215
Dune::DynamicVector< double > beta_
Definition: basicrow.hh:338
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const SourceTermType &sourceTerm, const NonlinearSolverParameterType &parameter)
constructor
Definition: basicrow.hh:84
void description(std::ostream &out) const
print description of ODE solver to out stream
Definition: basicrow.hh:303
NonlinearSolverType::LinearInverseOperatorType LinearInverseOperatorType
Definition: basicrow.hh:70
SourceTerm SourceTermType
Definition: basicrow.hh:67
int stages_
Definition: basicrow.hh:335
TimeStepControl TimeStepControlType
Definition: basicrow.hh:66
SourceTerm sourceTerm_
Definition: basicrow.hh:333
LinearInverseOperatorType linearSolver_
Definition: basicrow.hh:330
NonlinearSolver::ParameterType NonlinearSolverParameterType
Definition: basicrow.hh:69
std::vector< DestinationType * > update_
Definition: basicrow.hh:341
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:110
Dune::DynamicMatrix< double > alpha2_
Definition: basicrow.hh:337
HelmholtzOperatorType & helmholtzOp_
Definition: basicrow.hh:329
Dune::DynamicMatrix< double > alpha_
Definition: basicrow.hh:337
const PreconditionOperatorType * preconditioner_
Definition: basicrow.hh:345
BaseType::MonitorType MonitorType
Definition: basicrow.hh:61
BaseType::DestinationType DestinationType
Definition: basicrow.hh:62
int stages() const
Definition: basicrow.hh:301
TimeStepControl timeStepControl_
Definition: basicrow.hh:332
DestinationType temp_
Definition: basicrow.hh:340
NonlinearSolver NonlinearSolverType
Definition: basicrow.hh:65
double delta_
Definition: basicrow.hh:336
const int maxLinearIterations_
Definition: basicrow.hh:343
Dune::Fem::TimeProviderBase TimeProviderType
Definition: basicrow.hh:74
void initialize(const DestinationType &U0)
apply operator once to get dt estimate
Definition: basicrow.hh:197
HelmholtzOperator::SpaceOperatorType::PreconditionOperatorType PreconditionOperatorType
Definition: basicrow.hh:72
HelmholtzOperator HelmholtzOperatorType
Definition: basicrow.hh:64
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const NonlinearSolverParameterType &parameter)
constructor
Definition: basicrow.hh:149
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const NonlinearSolverParameterType &parameter)
constructor
Definition: basicrow.hh:126
double infNorm(const DestinationType &U, const DestinationType &Uerr) const
Definition: basicrow.hh:310
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: basicrow.hh:157
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: basicrow.hh:135
Dune::DynamicVector< double > gamma_
Definition: basicrow.hh:338
DestinationType rhs_
Definition: basicrow.hh:340
Dune::DynamicVector< double > c_
Definition: basicrow.hh:338
~BasicROWRungeKuttaSolver()
destructor
Definition: basicrow.hh:190
void setup(const ButcherTable &butcherTable)
Definition: basicrow.hh:165
general base for time providers
Definition: timeprovider.hh:36