1#ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_BASICROW_HH
2#define DUNE_FEM_SOLVER_RUNGEKUTTA_BASICROW_HH
14#include <dune/common/dynmatrix.hh>
15#include <dune/common/dynvector.hh>
30 bool operator() (
double time,
double timeStepSize,
int stage,
const T &u,
const std::vector< T * > &update, T &source )
36 void limit( T& update,
const double time ) {}
47 return std::numeric_limits< double >::max();
53 template<
class HelmholtzOperator,
class NonlinearSolver,
class TimeStepControl,
class SourceTerm = NoROWRungeKuttaSourceTerm >
83 template<
class ButcherTable >
86 const ButcherTable &butcherTable,
95 alpha_( butcherTable.A() ),
99 c_( butcherTable.c() ),
106 setup( butcherTable );
109 template<
class ButcherTable >
112 const ButcherTable &butcherTable,
125 template<
class ButcherTable >
128 const ButcherTable &butcherTable,
134 template<
class ButcherTable >
137 const ButcherTable &butcherTable,
148 template<
class ButcherTable >
151 const ButcherTable &butcherTable,
156 template<
class ButcherTable >
159 const ButcherTable &butcherTable,
164 template<
class ButcherTable >
165 void setup(
const ButcherTable& butcherTable )
169 std::cout <<
"ROW method of order=" << butcherTable.order() <<
" with " <<
stages_ <<
" stages" << std::endl;
173 for(
int i = 0; i <
stages(); ++i )
175 std::ostringstream name;
176 name <<
"RK stage " << i;
181 for(
int i = 0; i <
stages(); ++i )
192 for(
int i = 0; i <
stages(); ++i )
203 const double helmholtzEstimate =
helmholtzOp_.timeStepEstimate();
205 double sourceTermEstimate =
sourceTerm_.initialTimeStepEstimate( time, U0 );
207 if( sourceTermEstimate < 0.0 ) sourceTermEstimate = helmholtzEstimate ;
209 timeStepControl_.initialTimeStepSize( helmholtzEstimate, sourceTermEstimate );
219 typename HelmholtzOperatorType::JacobianOperatorType jOp(
"jacobianOperator", U.space(), U.space() );
223 assert( timeStepSize > 0.0 );
225 for(
int s = 0; s <
stages(); ++s )
231 for(
int k = 0; k < s; ++k )
234 updateStage *= (
gamma_[s]*timeStepSize);
235 for(
int k = 0; k < s; ++k )
238 rhs_.assign( updateStage );
241 const double stageTime = time +
c_[ s ]*timeStepSize;
273 for(
int s = 0; s <
stages(); ++s )
277 Uerr.axpy( -1.0, U );
278 const double errorU = Uerr.scalarProductDofs( Uerr );
279 const double normU = U.scalarProductDofs( U );
281 if( normU > 0 && errorU > 0 )
285 std::cout << std::scientific <<
"Error in RK = " << error <<
" norm " << errorU <<
" " << normU << std::endl;
291 for(
int s = 0; s <
stages(); ++s )
305 out <<
"Generic " <<
stages() <<
"-stage implicit Runge-Kutta solver.\\\\" << std::endl;
312 typedef typename DestinationType :: ConstDofIteratorType ConstDofIteratorType ;
313 const ConstDofIteratorType uend = U.dend();
315 for( ConstDofIteratorType u = U.dbegin(), uerr = Uerr.dbegin(); u != uend; ++u, ++uerr )
318 double uerrval = *uerr ;
321 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
static bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:445
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 ¶meter)
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 ¶meter=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 ¶meter)
constructor
Definition: basicrow.hh:149
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const NonlinearSolverParameterType ¶meter)
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 ¶meter=Dune::Fem::Parameter::container())
Definition: basicrow.hh:157
BasicROWRungeKuttaSolver(HelmholtzOperatorType &helmholtzOp, TimeProviderType &timeProvider, const ButcherTable &butcherTable, const TimeStepControlType &timeStepControl, const Dune::Fem::ParameterReader ¶meter=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