dune-fem  2.8-git
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
11 #include <dune/fem/io/parameter.hh>
16 
17 namespace 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:
100  TimeProviderBase& tp,
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 
131  TimeProviderBase& tp,
132  const int pord,
133  bool verbose )
134  : ExplicitRungeKuttaSolver( op, tp, defaultButcherTables( pord ), verbose )
135  {
136  }
137 
139  TimeProviderBase& tp,
140  const int pord,
142  : ExplicitRungeKuttaSolver( op, tp, pord, true )
143  {}
144 
146  void initialize(const DestinationType& U0)
147  {
148  if( ! initialized_ )
149  {
150  // Compute Steps
151  op_(U0, *(Upd[0]));
152  initialized_ = true;
153 
154  // provide operators time step estimate
155  tp_.provideTimeStepEstimate( op_.timeStepEstimate() );
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
183  tp_.provideTimeStepEstimate( op_.timeStepEstimate() );
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
200  tp_.provideTimeStepEstimate( op_.timeStepEstimate() );
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
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
Interface class for ODE Solver.
Definition: odesolverinterface.hh:21
DestinationImp DestinationType
type of destination
Definition: odesolverinterface.hh:62
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
SimpleButcherTable< double > defaultButcherTables(const int order) const
Definition: explicit.hh:75
Dune::DynamicMatrix< double > A_
Definition: explicit.hh:219
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
OdeSolverInterface< DestinationImp >::MonitorType MonitorType
Definition: explicit.hh:71
DestinationType ::DiscreteFunctionSpaceType SpaceType
Definition: explicit.hh:69
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