dune-fem  2.8-git
basicimplicit.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_BASICIMPLICIT_HH
2 #define DUNE_FEM_SOLVER_RUNGEKUTTA_BASICIMPLICIT_HH
3 
4 //- system includes
5 #include <cassert>
6 #include <cmath>
7 #include <limits>
8 #include <sstream>
9 #include <vector>
10 
11 //- dune-common includes
12 #include <dune/common/dynmatrix.hh>
13 #include <dune/common/dynvector.hh>
14 
15 //- dune-fem includes
17 
18 namespace DuneODE
19 {
20 
21  // NoImplicitRungeKuttaSourceTerm
22  // ------------------------------
23 
25  {
26  template< class T >
27  bool operator() ( double time, double timeStepSize, int stage, const T &u, const std::vector< T * > &update, T &source )
28  {
29  return false;
30  }
31 
32  template< class T >
33  void limit( T& update, const double time ) {}
34 
35  template< class T >
36  double initialTimeStepEstimate ( double time, const T &u ) const
37  {
38  // return negative value to indicate that implicit time step should be used
39  return -1.0;
40  }
41 
42  double timeStepEstimate () const
43  {
45  }
46  };
47 
48 
49 
51  template< class HelmholtzOperator, class NonlinearSolver, class TimeStepControl, class SourceTerm = NoImplicitRungeKuttaSourceTerm >
53  : public OdeSolverInterface< typename HelmholtzOperator::DomainFunctionType >
54  {
57 
58  public:
61 
62  typedef HelmholtzOperator HelmholtzOperatorType;
63  typedef NonlinearSolver NonlinearSolverType;
64  typedef TimeStepControl TimeStepControlType;
65  typedef SourceTerm SourceTermType;
66 
68 
69  typedef typename TimeStepControlType::ParameterType ParameterType;
70  typedef typename NonlinearSolver::ParameterType NonlinearSolverParameterType;
71 
79  template< class ButcherTable >
81  const ButcherTable &butcherTable,
82  const TimeStepControlType &timeStepControl,
83  const SourceTermType &sourceTerm,
84  const NonlinearSolverParameterType& parameters )
85  : helmholtzOp_( helmholtzOp ),
86  nonlinearSolver_( parameters ),
87  timeStepControl_( timeStepControl ),
88  sourceTerm_( sourceTerm ),
89  stages_( butcherTable.stages() ),
90  alpha_( butcherTable.A() ),
91  gamma_( stages() ),
92  beta_( stages() ),
93  c_( butcherTable.c() ),
94  rhs_( "RK rhs", helmholtzOp_.space() ),
96  update_( stages(), nullptr )
97  {
98  setup( butcherTable );
99  }
100 
107  template< class ButcherTable >
109  const ButcherTable &butcherTable,
110  const TimeStepControlType &timeStepControl,
111  const NonlinearSolverParameterType& parameters )
112  : helmholtzOp_( helmholtzOp ),
113  nonlinearSolver_( parameters ),
114  timeStepControl_( timeStepControl ),
115  stages_( butcherTable.stages() ),
116  alpha_( butcherTable.A() ),
117  gamma_( stages() ),
118  beta_( stages() ),
119  c_( butcherTable.c() ),
120  rhs_( "RK rhs", helmholtzOp_.space() ),
121  updateStorage_(),
122  update_( stages(), nullptr )
123  {
124  setup( butcherTable );
125  }
126 
127  template< class ButcherTable >
129  const ButcherTable &butcherTable,
130  const TimeStepControlType &timeStepControl,
132  : BasicImplicitRungeKuttaSolver( helmholtzOp, butcherTable, timeStepControl,
133  NonlinearSolverParameterType( parameter ) )
134  {
135  }
136 
137  template< class ButcherTable >
139  const ButcherTable &butcherTable,
141  : BasicImplicitRungeKuttaSolver( helmholtzOp, butcherTable,
142  TimeStepControlType( parameter ), NonlinearSolverParameterType( parameter ) )
143  {
144  }
145 
146  template< class ButcherTable >
147  void setup( const ButcherTable& butcherTable )
148  {
149  update_.clear();
150  update_.resize( stages(), nullptr );
151  updateStorage_.resize( stages() );
152 
153  // create intermediate functions
154  for( int i = 0; i < stages(); ++i )
155  {
156  std::ostringstream name;
157  name << "RK stage " << i;
158  updateStorage_[ i ].reset( new DestinationType( name.str(), helmholtzOp_.space() ) );
159  update_[ i ] = updateStorage_[ i ].operator ->();
160  }
161 
162  // compute coefficients
163  Dune::DynamicMatrix< double > AL( alpha_ );
164  for( int i = 0; i < stages(); ++i )
165  {
166  gamma_[ i ] = AL[ i ][ i ];
167  AL[ i ][ i ] = 0.0;
168  }
169 
170  alpha_.invert();
171  alpha_.mtv( butcherTable.b(), beta_ );
172 
173  alpha_.leftmultiply( AL );
174  for( int i = 0; i < stages(); ++i )
175  alpha_[ i ][ i ] = gamma_[ i ];
176 
177  for( int i = 0; i < stages(); ++i )
178  {
179  gamma_[ i ] = 1.0;
180  for( int j = 0; j < i; ++j )
181  gamma_[ i ] -= alpha_[ i ][ j ];
182  }
183 
184  delta_ = 1.0;
185  for( int i = 0; i < stages(); ++i )
186  delta_ -= beta_[ i ];
187 
188  }
189 
191  void initialize ( const DestinationType &U0 )
192  {
193  const double time = timeStepControl_.time();
194 
195  helmholtzOp_.setTime( time );
196  helmholtzOp_.initializeTimeStepSize( U0 );
197  const double helmholtzEstimate = helmholtzOp_.timeStepEstimate();
198 
199  double sourceTermEstimate = sourceTerm_.initialTimeStepEstimate( time, U0 );
200  // negative time step is given by the empty source term
201  if( sourceTermEstimate < 0.0 ) sourceTermEstimate = helmholtzEstimate ;
202 
203  timeStepControl_.initialTimeStepSize( helmholtzEstimate, sourceTermEstimate );
204  }
205 
206  using BaseType::solve;
207 
209  void solve ( DestinationType &U, MonitorType &monitor )
210  {
211  monitor.reset();
212 
213  const double time = timeStepControl_.time();
214  const double timeStepSize = timeStepControl_.timeStepSize();
215  assert( timeStepSize > 0.0 );
216  for( int s = 0; s < stages(); ++s )
217  {
218  assert( update_[ s ] );
219  // update for stage s
220  DestinationType& updateStage = *update_[ s ];
221 
222  // assemble rhs of nonlinear equation
223  updateStage.assign( U );
224  updateStage *= gamma_[ s ];
225  for( int k = 0; k < s; ++k )
226  updateStage.axpy( alpha_[ s ][ k ], *update_[ k ] );
227 
228  const double stageTime = time + c_[ s ]*timeStepSize;
229  if( sourceTerm_( time, timeStepSize, s, U, update_, rhs_ ) )
230  {
231  updateStage.axpy( alpha_[ s ][ s ]*timeStepSize, rhs_ );
232  sourceTerm_.limit( updateStage, stageTime );
233  }
234 
235  // apply Helmholtz operator to right hand side
236  helmholtzOp_.setTime( stageTime );
237  helmholtzOp_.setLambda( 0 );
238  helmholtzOp_( updateStage, rhs_ );
239 
240  // solve the system
241  helmholtzOp_.setLambda( alpha_[ s ][ s ]*timeStepSize );
243  nonlinearSolver_( rhs_, updateStage );
244  nonlinearSolver_.unbind();
245 
246  // update monitor
247  monitor.newtonIterations_ += nonlinearSolver_.iterations();
248  monitor.linearSolverIterations_ += nonlinearSolver_.linearIterations();
249 
250  // on failure break solving
251  if( !nonlinearSolver_.converged() )
252  return timeStepControl_.reduceTimeStep( helmholtzOp_.timeStepEstimate(), sourceTerm_.timeStepEstimate(), monitor );
253  }
254 
255  double error = 0.0;
256  if( timeStepControl_.computeError() )
257  {
258  // store U (to be revised)
259  DestinationType Uerr( U );
260 
261  // update solution
262  U *= delta_;
263  for( int s = 0; s < stages(); ++s )
264  U.axpy( beta_[ s ], *update_[ s ] );
265 
266  //error = infNorm( U, Uerr );
267  Uerr.axpy( -1.0, U );
268  const double errorU = Uerr.scalarProductDofs( Uerr );
269  const double normU = U.scalarProductDofs( U );
270 
271  if( normU > 0 && errorU > 0 )
272  {
273  error = std::sqrt( errorU / normU );
274  }
275  std::cout << std::scientific << "Error in RK = " << error << " norm " << errorU << " " << normU << std::endl;
276  //std::cout << std::scientific << "Error in RK = " << error << std::endl;
277  }
278  else
279  {
280  // update solution
281  U *= delta_;
282  for( int s = 0; s < stages(); ++s )
283  U.axpy( beta_[ s ], *update_[ s ] );
284  }
285  // set error to monitor
286  monitor.error_ = error;
287 
288  // update time step size
289  timeStepControl_.timeStepEstimate( helmholtzOp_.timeStepEstimate(), sourceTerm_.timeStepEstimate(), monitor );
290  }
291 
292  int stages () const { return stages_; }
293 
294  void description ( std::ostream &out ) const
295  {
296  out << "Generic " << stages() << "-stage implicit Runge-Kutta solver.\\\\" << std::endl;
297  }
298 
299  protected:
300  double infNorm(const DestinationType& U, const DestinationType& Uerr ) const
301  {
302  typedef typename DestinationType :: ConstDofIteratorType ConstDofIteratorType ;
303  const ConstDofIteratorType uend = U.dend();
304  double res = 0;
305  for( ConstDofIteratorType u = U.dbegin(), uerr = Uerr.dbegin(); u != uend; ++u, ++uerr )
306  {
307  double uval = *u;
308  double uerrval = *uerr ;
309  double div = std::abs( std::max( uval, uerrval ) );
310 
311  double norm = std::abs( uval - uerrval );
312  if( std::abs(div) > 1e-12 )
313  norm /= div;
314  res = std::max( res, norm );
315  }
316  return res;
317  }
318 
321  TimeStepControl timeStepControl_;
322  SourceTerm sourceTerm_;
323 
324  int stages_;
325  double delta_;
326  Dune::DynamicMatrix< double > alpha_;
327  Dune::DynamicVector< double > gamma_, beta_, c_;
328 
330  std::vector< std::unique_ptr< DestinationType > > updateStorage_;
331  std::vector< DestinationType* > update_;
332  };
333 
334 } // namespace DuneODE
335 
336 #endif // #ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_BASICIMPLICIT_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 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
int newtonIterations_
Definition: odesolverinterface.hh:34
double error_
Definition: odesolverinterface.hh:30
void reset()
Definition: odesolverinterface.hh:43
int linearSolverIterations_
Definition: odesolverinterface.hh:35
Definition: basicimplicit.hh:25
void limit(T &update, const double time)
Definition: basicimplicit.hh:33
double initialTimeStepEstimate(double time, const T &u) const
Definition: basicimplicit.hh:36
double timeStepEstimate() const
Definition: basicimplicit.hh:42
bool operator()(double time, double timeStepSize, int stage, const T &u, const std::vector< T * > &update, T &source)
Definition: basicimplicit.hh:27
Implicit RungeKutta ODE solver.
Definition: basicimplicit.hh:54
void description(std::ostream &out) const
print description of ODE solver to out stream
Definition: basicimplicit.hh:294
std::vector< std::unique_ptr< DestinationType > > updateStorage_
Definition: basicimplicit.hh:330
BasicImplicitRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, const ButcherTable &butcherTable, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: basicimplicit.hh:138
double delta_
Definition: basicimplicit.hh:325
DestinationType rhs_
Definition: basicimplicit.hh:329
NonlinearSolver::ParameterType NonlinearSolverParameterType
Definition: basicimplicit.hh:70
Dune::Fem::TimeProviderBase TimeProviderType
Definition: basicimplicit.hh:67
Dune::DynamicMatrix< double > alpha_
Definition: basicimplicit.hh:326
double infNorm(const DestinationType &U, const DestinationType &Uerr) const
Definition: basicimplicit.hh:300
int stages() const
Definition: basicimplicit.hh:292
BasicImplicitRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const SourceTermType &sourceTerm, const NonlinearSolverParameterType &parameters)
constructor
Definition: basicimplicit.hh:80
int stages_
Definition: basicimplicit.hh:324
Dune::DynamicVector< double > beta_
Definition: basicimplicit.hh:327
Dune::DynamicVector< double > gamma_
Definition: basicimplicit.hh:327
BaseType::DestinationType DestinationType
Definition: basicimplicit.hh:60
HelmholtzOperatorType & helmholtzOp_
Definition: basicimplicit.hh:319
SourceTerm sourceTerm_
Definition: basicimplicit.hh:322
void initialize(const DestinationType &U0)
apply operator once to get dt estimate
Definition: basicimplicit.hh:191
NonlinearSolverType nonlinearSolver_
Definition: basicimplicit.hh:320
SourceTerm SourceTermType
Definition: basicimplicit.hh:65
TimeStepControl TimeStepControlType
Definition: basicimplicit.hh:64
BaseType::MonitorType MonitorType
Definition: basicimplicit.hh:59
Dune::DynamicVector< double > c_
Definition: basicimplicit.hh:327
std::vector< DestinationType * > update_
Definition: basicimplicit.hh:331
NonlinearSolver NonlinearSolverType
Definition: basicimplicit.hh:63
TimeStepControlType::ParameterType ParameterType
Definition: basicimplicit.hh:69
BasicImplicitRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const NonlinearSolverParameterType &parameters)
constructor
Definition: basicimplicit.hh:108
BasicImplicitRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: basicimplicit.hh:128
void solve(DestinationType &U, MonitorType &monitor)
solve the system
Definition: basicimplicit.hh:209
TimeStepControl timeStepControl_
Definition: basicimplicit.hh:321
HelmholtzOperator HelmholtzOperatorType
Definition: basicimplicit.hh:62
void setup(const ButcherTable &butcherTable)
Definition: basicimplicit.hh:147
general base for time providers
Definition: timeprovider.hh:36