dune-fem 2.8.0
Loading...
Searching...
No Matches
explicit.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_EXPLICIT_HH
2#define DUNE_FEM_SOLVER_RUNGEKUTTA_EXPLICIT_HH
3
4//- system includes
5#include <iostream>
6#include <cmath>
7#include <vector>
8#include <cassert>
9
10//- Dune includes
16
17namespace DuneODE
18{
19
20 using namespace Dune;
21 using namespace Fem;
22 using namespace std;
23
62 template<class DestinationImp>
64 : public OdeSolverInterface<DestinationImp>
65 {
66 public:
67 typedef DestinationImp DestinationType;
69 typedef typename DestinationType :: DiscreteFunctionSpaceType SpaceType;
70
72
73 using OdeSolverInterface<DestinationImp> :: solve ;
74 protected:
76 {
77 switch( order )
78 {
79 case 1: return explicitEulerButcherTable();
80 case 2: return tvd2ButcherTable();
81 case 3: return tvd3ButcherTable();
82 case 4: return rk4ButcherTable();
83 case 5:
84 case 6: return expl6ButcherTable();
85
86 default:
87 std::cerr<< "Warning: ExplicitRungeKutta of order "<< order << " not implemented, using order 6!" << std::endl;
88 return expl6ButcherTable();
89 }
90 }
91
92 public:
101 const SimpleButcherTable< double >& butcherTable,
102 bool verbose )
103 : A_( butcherTable.A() ),
104 b_( butcherTable.b() ),
105 c_( butcherTable.c() ),
106 Upd(),
107 op_(op),
108 tp_(tp),
109 ord_( butcherTable.order() ),
110 stages_( butcherTable.stages() ),
111 initialized_(false)
112 {
113 assert(ord_>0);
114
115 // create update memory
116 for (int i=0; i<stages_; ++i)
117 {
118 Upd.emplace_back( new DestinationType("URK",op_.space()) );
119 }
120 Upd.emplace_back(new DestinationType("Ustep",op_.space()) );
121
122 }
123
132 const int pord,
133 bool verbose )
134 : ExplicitRungeKuttaSolver( op, tp, defaultButcherTables( pord ), verbose )
135 {
136 }
137
140 const int pord,
142 : ExplicitRungeKuttaSolver( op, tp, pord, true )
143 {}
144
147 {
148 if( ! initialized_ )
149 {
150 // Compute Steps
151 op_(U0, *(Upd[0]));
152 initialized_ = true;
153
154 // provide operators time step estimate
156 }
157 }
158
160 void solve(DestinationType& U0, MonitorType& monitor )
161 {
162 // no information here
163 monitor.reset();
164
165 // initialize
166 if( ! initialized_ )
167 {
168 DUNE_THROW(InvalidStateException,"ExplicitRungeKuttaSolver wasn't initialized before first call!");
169 }
170
171 // get cfl * timeStepEstimate
172 const double dt = tp_.deltaT();
173 // get time
174 const double t = tp_.time();
175
176 // set new time
177 op_.setTime( t );
178
179 // Compute Steps
180 op_(U0, *(Upd[0]));
181
182 // provide operators time step estimate
184
185 for (int i=1; i<stages_; ++i)
186 {
187 (Upd[ord_])->assign(U0);
188 for (int j=0; j<i ; ++j)
189 {
190 (Upd[ord_])->axpy((A_[i][j]*dt), *(Upd[j]));
191 }
192
193 // set new time
194 op_.setTime( t + c_[i]*dt );
195
196 // apply operator
197 op_( *(Upd[ord_]), *(Upd[i]) );
198
199 // provide operators time step estimate
201 }
202
203 // Perform Update
204 for (int j=0; j<stages_; ++j)
205 {
206 U0.axpy((b_[j]*dt), *(Upd[j]));
207 }
208 }
209
210 void description(std::ostream& out) const
211 {
212 out << "ExplRungeKutta, steps: " << ord_
213 //<< ", cfl: " << this->tp_.factor()
214 << "\\\\" <<std::endl;
215 }
216
217 protected:
218 // Butcher table A,b,c
219 Dune::DynamicMatrix< double > A_;
220 Dune::DynamicVector< double > b_;
221 Dune::DynamicVector< double > c_;
222
223 // stages of Runge-Kutta solver
224 std::vector< std::unique_ptr< DestinationType > > Upd;
225
226 // operator to solve for
228 // time provider
230
231 // order of RK solver
232 const int ord_;
233 // number of stages
234 const int stages_;
235 // init flag
237 };
238
241} // namespace DuneODE
242
243#endif // #ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_EXPLICIT_HH
STL namespace.
Definition: bindguard.hh:11
void axpy(const T &a, const T &x, T &y)
Definition: space/basisfunctionset/functor.hh:38
Definition: multistep.hh:17
SimpleButcherTable< double > expl6ButcherTable()
Definition: butchertable.cc:76
SimpleButcherTable< double > explicitEulerButcherTable()
Definition: butchertable.cc:16
SimpleButcherTable< double > tvd3ButcherTable()
Definition: butchertable.cc:43
SimpleButcherTable< double > tvd2ButcherTable()
Definition: butchertable.cc:29
SimpleButcherTable< double > rk4ButcherTable()
Definition: butchertable.cc:58
static ParameterContainer & container()
Definition: io/parameter.hh:193
interface for time evolution operators
Definition: spaceoperatorif.hh:38
virtual void setTime(const double time)
set time for operators
Definition: spaceoperatorif.hh:76
virtual const DiscreteFunctionSpaceType & space() const =0
return reference to space (needed by ode solvers)
virtual double timeStepEstimate() const
estimate maximum time step
Definition: spaceoperatorif.hh:86
Interface class for ODE Solver.
Definition: odesolverinterface.hh:21
Definition: odesolverinterface.hh:27
void reset()
Definition: odesolverinterface.hh:43
Definition: butchertable.hh:17
Exlicit RungeKutta ODE solver.
Definition: explicit.hh:65
bool initialized_
Definition: explicit.hh:236
OperatorType & op_
Definition: explicit.hh:227
OdeSolverInterface< DestinationImp >::MonitorType MonitorType
Definition: explicit.hh:71
SimpleButcherTable< double > defaultButcherTables(const int order) const
Definition: explicit.hh:75
Dune::DynamicMatrix< double > A_
Definition: explicit.hh:219
DestinationType::DiscreteFunctionSpaceType SpaceType
Definition: explicit.hh:69
ExplicitRungeKuttaSolver(OperatorType &op, TimeProviderBase &tp, const int pord, bool verbose)
constructor
Definition: explicit.hh:130
DestinationImp DestinationType
Definition: explicit.hh:67
void description(std::ostream &out) const
print description of ODE solver to out stream
Definition: explicit.hh:210
std::vector< std::unique_ptr< DestinationType > > Upd
Definition: explicit.hh:224
const int stages_
Definition: explicit.hh:234
TimeProviderBase & tp_
Definition: explicit.hh:229
ExplicitRungeKuttaSolver(OperatorType &op, TimeProviderBase &tp, const SimpleButcherTable< double > &butcherTable, bool verbose)
constructor
Definition: explicit.hh:99
SpaceOperatorInterface< DestinationImp > OperatorType
Definition: explicit.hh:68
Dune::DynamicVector< double > b_
Definition: explicit.hh:220
void solve(DestinationType &U0, MonitorType &monitor)
solve the system
Definition: explicit.hh:160
ExplicitRungeKuttaSolver(OperatorType &op, TimeProviderBase &tp, const int pord, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: explicit.hh:138
Dune::DynamicVector< double > c_
Definition: explicit.hh:221
void initialize(const DestinationType &U0)
apply operator once to get dt estimate
Definition: explicit.hh:146
const int ord_
Definition: explicit.hh:232
general base for time providers
Definition: timeprovider.hh:36
double time() const
obtain the current time
Definition: timeprovider.hh:94
void provideTimeStepEstimate(const double dtEstimate)
set time step estimate to minimum of given value and internal time step estiamte
Definition: timeprovider.hh:142
double deltaT() const
obtain the size of the current time step
Definition: timeprovider.hh:113