1#ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_BASICIMPLICIT_HH
2#define DUNE_FEM_SOLVER_RUNGEKUTTA_BASICIMPLICIT_HH
12#include <dune/common/dynmatrix.hh>
13#include <dune/common/dynvector.hh>
27 bool operator() (
double time,
double timeStepSize,
int stage,
const T &u,
const std::vector< T * > &update, T &source )
33 void limit( T& update,
const double time ) {}
44 return std::numeric_limits< double >::max();
51 template<
class HelmholtzOperator,
class NonlinearSolver,
class TimeStepControl,
class SourceTerm = NoImplicitRungeKuttaSourceTerm >
79 template<
class ButcherTable >
81 const ButcherTable &butcherTable,
90 alpha_( butcherTable.A() ),
93 c_( butcherTable.c() ),
98 setup( butcherTable );
107 template<
class ButcherTable >
109 const ButcherTable &butcherTable,
116 alpha_( butcherTable.A() ),
119 c_( butcherTable.c() ),
124 setup( butcherTable );
127 template<
class ButcherTable >
129 const ButcherTable &butcherTable,
137 template<
class ButcherTable >
139 const ButcherTable &butcherTable,
146 template<
class ButcherTable >
147 void setup(
const ButcherTable& butcherTable )
154 for(
int i = 0; i <
stages(); ++i )
156 std::ostringstream name;
157 name <<
"RK stage " << i;
163 Dune::DynamicMatrix< double > AL(
alpha_ );
164 for(
int i = 0; i <
stages(); ++i )
166 gamma_[ i ] = AL[ i ][ i ];
173 alpha_.leftmultiply( AL );
174 for(
int i = 0; i <
stages(); ++i )
177 for(
int i = 0; i <
stages(); ++i )
180 for(
int j = 0; j < i; ++j )
185 for(
int i = 0; i <
stages(); ++i )
197 const double helmholtzEstimate =
helmholtzOp_.timeStepEstimate();
199 double sourceTermEstimate =
sourceTerm_.initialTimeStepEstimate( time, U0 );
201 if( sourceTermEstimate < 0.0 ) sourceTermEstimate = helmholtzEstimate ;
203 timeStepControl_.initialTimeStepSize( helmholtzEstimate, sourceTermEstimate );
215 assert( timeStepSize > 0.0 );
216 for(
int s = 0; s <
stages(); ++s )
223 updateStage.assign( U );
224 updateStage *=
gamma_[ s ];
225 for(
int k = 0; k < s; ++k )
228 const double stageTime = time +
c_[ s ]*timeStepSize;
231 updateStage.axpy(
alpha_[ s ][ s ]*timeStepSize,
rhs_ );
263 for(
int s = 0; s <
stages(); ++s )
267 Uerr.axpy( -1.0, U );
268 const double errorU = Uerr.scalarProductDofs( Uerr );
269 const double normU = U.scalarProductDofs( U );
271 if( normU > 0 && errorU > 0 )
275 std::cout << std::scientific <<
"Error in RK = " << error <<
" norm " << errorU <<
" " << normU << std::endl;
282 for(
int s = 0; s <
stages(); ++s )
296 out <<
"Generic " <<
stages() <<
"-stage implicit Runge-Kutta solver.\\\\" << std::endl;
302 typedef typename DestinationType :: ConstDofIteratorType ConstDofIteratorType ;
303 const ConstDofIteratorType uend = U.dend();
305 for( ConstDofIteratorType u = U.dbegin(), uerr = Uerr.dbegin(); u != uend; ++u, ++uerr )
308 double uerrval = *uerr ;
311 double norm =
std::abs( uval - uerrval );
double sqrt(const Dune::Fem::Double &v)
Definition: double.hh:977
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
double max(const Dune::Fem::Double &v, const double p)
Definition: double.hh:965
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 ¶meter=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 ¶meters)
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 ¶meters)
constructor
Definition: basicimplicit.hh:108
BasicImplicitRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const Dune::Fem::ParameterReader ¶meter=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