DUNE-FEM (unstable)

multistep.hh
1 #ifndef MULTISTEP_SOLVER_HH
2 #define MULTISTEP_SOLVER_HH
3 
4 // inlcude all used headers before, that they don not appear in DuneODE
5 
6 //- system includes
7 #include <iostream>
8 #include <cmath>
9 #include <vector>
10 #include <cassert>
11 
12 //- Dune includes
13 #include <dune/fem/solver/timeprovider.hh>
14 #include <dune/fem/solver/rungekutta/explicit.hh>
15 
16 namespace DuneODE
17 {
18 
19  using namespace Dune;
20  using namespace Fem;
21  using namespace std;
22 
30  template<class Operator>
31  class ExplMultiStepBase
32  {
33  public:
34  typedef typename Operator::DestinationType DestinationType;
35  typedef typename DestinationType :: DiscreteFunctionSpaceType SpaceType;
36  // typedef typename SpaceType :: GridType :: Traits :: Communication DuneCommunicatorType;
37  protected:
38  std::vector< std::vector<double> > a;
39  std::vector<double> b;
40  std::vector<double> c;
41  size_t steps_;
42  double gamma_;
43  std::vector<DestinationType*> Uj;
44  std::vector<DestinationType*> Fj;
45  std::vector<double> deltat_;
46  bool msInit,msFirst;
47  protected:
48  const int ord_;
49 
50  public:
57  ExplMultiStepBase(Operator& op, TimeProviderBase& tp,
58  int pord, bool verbose = true ) :
59  a(0),b(0),c(0)
60  , steps_(pord+1)
61  , gamma_(0.5) // 1./3. -> Shu
62  , Uj(0)
63  , Fj(0)
64  , deltat_(steps_)
65  , msInit(false)
66  , msFirst(true)
67  , ord_(pord)
68  , op_(op)
69  , tp_(tp)
70  , initialized_(false)
71  {
72  a.resize(ord_);
73  for (int i=0; i<ord_; ++i)
74  {
75  a[i].resize(ord_);
76  }
77  b.resize(ord_);
78  c.resize(ord_);
79 
80  switch (ord_) {
81  case 4 :
82  a[0][0]=0.; a[0][1]=0.; a[0][2]=0.; a[0][3]=0.;
83  a[1][0]=1.0; a[1][1]=0.; a[1][2]=0.; a[1][3]=0.;
84  a[2][0]=0.25; a[2][1]=0.25; a[2][2]=0.; a[2][3]=0.;
85  a[3][0]=1./6.; a[3][1]=1./6.; a[3][2]=2./3.; a[3][3]=0.;
86  b[0]=1./6.; b[1]=1./6.; b[2]=2./3.; b[3]=0.;
87  c[0]=0.; c[1]=1.0; c[2]=0.5; c[3]=1.0;
88  break;
89  case 3 :
90  a[0][0]=0.; a[0][1]=0.; a[0][2]=0.;
91  a[1][0]=1.0; a[1][1]=0.; a[1][2]=0.;
92  a[2][0]=0.25; a[2][1]=0.25; a[2][2]=0.;
93  b[0]=1./6.; b[1]=1./6.; b[2]=2./3.;
94  c[0]=0.; c[1]=1; c[2]=0.5;
95  break;
96  case 2 :
97  a[0][0]=0.; a[0][1]=0.;
98  a[1][0]=1.0; a[1][1]=0.;
99  b[0]=0.5; b[1]=0.5;
100  c[0]=0; c[1]=1;
101  break;
102  case 1:
103  a[0][0]=0.;
104  b[0]=1.;
105  c[0]=0.;
106  break;
107  default : std::cerr << "Runge-Kutta method of this order not implemented"
108  << std::endl;
109  abort();
110  }
111  // create memory for intermediate stages
112  for (size_t i=0; i<steps_; ++i)
113  {
114  Fj.push_back(new DestinationType("FMS",op_.space()) );
115  }
116  Uj.push_back(new DestinationType("UMS",op_.space()) );
117  }
118 
120  ~ExplMultiStepBase()
121  {
122  for(size_t i=0; i<Fj.size(); ++i) {
123  delete Fj[i];
124  }
125  for(size_t i=0; i<Uj.size(); ++i) {
126  delete Uj[i];
127  }
128  }
129 
131  void initialize(const DestinationType& U0)
132  {
133  if( ! initialized_ )
134  {
135  // Compute Steps
136  op_(U0, (*Fj[0]));
137  initialized_ = true;
138  }
139  }
140 
142  void solve(DestinationType& U0)
143  {
144  if (!msInit) {
145  // Startup phase using RungeKutta
146  Uj[Uj.size()-1]->assign(U0);
147  solveRK(U0);
148  deltat_[Uj.size()-1] = tp_.deltaT();
149  Uj.push_back(new DestinationType("UMS",op_.space()) );
150  if (Uj.size()==steps_) {
151  Uj[steps_-1]->assign(U0);
152  for (size_t i=0;i<steps_-1;i++) {
153  op_(*(Uj[i]),(*Fj[i]));
154  }
155  msInit=true;
156  }
157  return;
158  }
159  // now multistep
160 
161  // get cfl * timeStepEstimate
162  deltat_[steps_-1] = tp_.deltaT();
163  // get time
164 
165  // Compute Steps
166  Uj[steps_-1]->assign(U0);
167 
168  // set new time
169  op_.setTime( tp_.time() );
170 
171  op_(*(Uj[steps_-1]), *(Fj[steps_-1]));
172 
173  // provide time step estimate
174  tp_.provideTimeStepEstimate( op_.timeStepEstimate() );
175 
176  // Perform Update
177  double alpha[steps_];
178  double beta[steps_];
179  switch (steps_) {
180  case 2: {
181  std::cerr << "the two step scheme of order 2 has negative beta coeffs and requires adjoint space operator!" << std::endl;
182  double w = deltat_[steps_-1]/deltat_[steps_-2];
183  double alpha2 = (1.+2.*gamma_*w)/(1.+w);
184  alpha[1] = -((1.-2.*gamma_)*w-1.)/alpha2;
185  alpha[0] = -((2.*gamma_-1.)*w*w/(1.+w))/alpha2;
186  beta[1] = (1.+gamma_*w)/alpha2*deltat_[steps_-1];
187  beta[0] = -(gamma_*w)/alpha2*deltat_[steps_-1];
188  std::cout << "# MS Coeffs : "
189  << alpha[0] << " " << alpha[1] << " "
190  << alpha2 << " "
191  << beta[0] << " " << beta[1] << " "
192  << deltat_[0] << " " << deltat_[1]
193  << std::endl;
194  } break;
195  case 3: {
196  /* Shu
197  alpha[0] = 1./4.;
198  alpha[1] = 0.;
199  alpha[2] = 3./4.;
200  beta[0] = 0.*deltat_[steps_-1];
201  beta[1] = 0.*deltat_[steps_-1];
202  beta[2] = 3./2.*deltat_[steps_-1];
203  */
204  double w1 = deltat_[1]/deltat_[0];
205  double w2 = deltat_[2]/deltat_[1];
206  double g = w1*w2/(1+w1);
207  // CFL < 1-g -> TVD
208  // therefore need g<1!
209  assert(g<1);
210  alpha[0] = g*g;
211  alpha[1] = 0.;
212  alpha[2] = 1.-g*g; // > 0 due to CFL
213  beta[0] = 0.*deltat_[steps_-1];
214  beta[1] = 0.*deltat_[steps_-1];
215  beta[2] = (1+g)*deltat_[steps_-1]; // always > 0
216  /*
217  double w1 = deltat_[1]/deltat_[0];
218  double w2 = deltat_[2]/deltat_[1];
219  double w1_1=1.+w1;
220  alpha[0] = (w1*w2-gamma_*w1*w1)/(w1_1*w1_1);
221  alpha[1] = gamma_;
222  alpha[2] = 1.-alpha[1]-alpha[0];
223  beta[0] = 0.*deltat_[steps_-1];
224  beta[1] = 0.*deltat_[steps_-1];
225  beta[2] = (1+w1+w1*w2-gamma_/w2*(1.+2.*w1))/w1_1*
226  deltat_[steps_-1]; // always > 0
227  */
228  } break;
229  case 4: {
230 
231  alpha[0] = 11./27.; // 1./9.;
232  alpha[1] = 0.;
233  alpha[2] = 0.;
234  alpha[3] = 16./27.; // 8./9.;
235  beta[0] = 4./9.*deltat_[steps_-1]; // 4./9.
236  beta[1] = 0.*deltat_[steps_-1];
237  beta[2] = 0.*deltat_[steps_-1];
238  beta[3] = 16./9.*deltat_[steps_-1]; // 4./3.
239  /*
240  double k3 = deltat_[3];
241  double K2 = deltat_[0]+deltat_[1]+deltat_[2];
242  double K3 = K2 + k3;
243  double q1 = K3/K2;
244  double q2 = k3/K2;
245  double g = q2*q2*(3.+2.*q2);
246  alpha[0] = g;
247  alpha[1] = 0.;
248  alpha[2] = 0.;
249  alpha[3] = 1.-g;
250  beta[0] = q2*q1*deltat_[steps_-1]; // 0.
251  beta[1] = 0.*deltat_[steps_-1];
252  beta[2] = 0.*deltat_[steps_-1];
253  beta[3] = q1*q1*deltat_[steps_-1]; // 4./3.
254  */
255  } break;
256  }
257  #if 0 // !NDEBUG
258  { // Check order conditions
259  int p = 1;
260  double cond;
261  double M[steps_+1];
262  for (int j=0;j<=steps_;j++) {
263  M[j] = 0;
264  for (int i=0;i<j;i++)
265  M[j] += deltat_[i];
266  }
267  do {
268  cond = -1.;
269  for (int q=0;q<p;q++)
270  cond *= M[steps_];
271  for (int j=1;j<steps_;j++) {
272  double Mp=1;
273  for (int q=0;q<p-1;q++)
274  Mp *= M[j];
275  cond += (alpha[j]*M[j]+p*beta[j])*Mp;
276  }
277  if (p==1)
278  cond += beta[0];
279  std::cout << "# MS Cond for p= " << p
280  << " " << cond
281  << std::endl;
282  ++p;
283  } while (fabs(cond)<1e-7 && p<10);
284  std::cout << " # CFL and stability: ";
285  for (int i=0;i<steps_;i++) {
286  if (fabs(beta[i])>1e-10) {
287  std::cout << alpha[i]/beta[i]*
288  deltat_[steps_-1] << " ";
289  }
290  }
291  std::cout << std::endl;
292  }
293  #endif
294  // Update
295  U0 *= alpha[steps_-1];
296  U0.axpy(beta[steps_-1], *(Fj[steps_-1]) );
297  for (size_t i=0;i<steps_-1;i++) {
298  U0.axpy(alpha[i], *(Uj[i]));
299  U0.axpy(beta[i], *(Fj[i]));
300  }
301 
302  DestinationType* Utmp = Uj[0];
303  DestinationType* Ftmp = Fj[0];
304  for (size_t i=1;i<steps_;i++) {
305  deltat_[i-1] = deltat_[i];
306  Uj[i-1] = Uj[i];
307  Fj[i-1] = Fj[i];
308  }
309  Uj[steps_-1] = Utmp;
310  Fj[steps_-1] = Ftmp;
311  }
313  void solveRK(DestinationType& U0)
314  {
315  DestinationType Uval("Utmp",op_.space());
316 
317  // get cfl * timeStepEstimate
318  const double dt = tp_.deltaT();
319  // get time
320  const double t = tp_.time();
321 
322  // set new time
323  op_.setTime( t );
324 
325  // Compute Steps
326  op_(U0, *(Fj[0]));
327 
328  // provide time step estimate
329  tp_.provideTimeStepEstimate( op_.timeStepEstimate() );
330 
331  for (int i=1; i<ord_; ++i)
332  {
333  Uval.assign(U0);
334  for (int j=0; j<i ; ++j)
335  {
336  Uval.axpy((a[i][j]*dt), *(Fj[j]));
337  }
338 
339  // set new time
340  op_.setTime( t + c[i]*dt );
341 
342  // apply operator
343  op_(Uval,*(Fj[i]));
344 
345  // provide time step estimate
346  tp_.provideTimeStepEstimate( op_.timeStepEstimate() );
347  }
348 
349  // Perform Update
350  for (int j=0; j<ord_; ++j)
351  {
352  U0.axpy((b[j]*dt), *(Fj[j]));
353  }
354  }
355 
356  protected:
357  // operator to solve for
358  const Operator& op_;
359  // time provider
360  TimeProviderBase& tp_;
361  // init flag
362  bool initialized_;
363  };
364 
365 #if 0
368  template<class Operator>
369  class ExplMultiStep : public TimeProvider ,
370  public ExplMultiStepBase<Operator>
371  {
372  typedef ExplMultiStepBase<Operator> BaseType;
373  public:
374  typedef typename Operator :: DestinationType DestinationType;
375  typedef typename DestinationType :: DiscreteFunctionSpaceType SpaceType;
376  typedef typename SpaceType :: GridType :: Traits :: Communication DuneCommunicatorType;
377 
378  public:
385  ExplMultiStep (Operator& op,int pord,double cfl, bool verbose = true ) :
386  TimeProvider(0.0,(cfl/double(pord))*0.95),
387  BaseType(op,*this,pord,verbose),
388  tp_(op.space().gridPart().comm(),*this),
389  savetime_(0.0), savestep_(1)
390  {
391  op.timeProvider(this);
392  }
393 
394  void initialize(const DestinationType& U0)
395  {
396  if(! this->initialized_)
397  {
398  // initialize
399  BaseType :: initialize(U0);
400 
401  // global min of dt and reset of dtEstimate
402  this->tp_.syncTimeStep();
403  }
404  }
405 
406  double solve(typename Operator::DestinationType& U0)
407  {
408  initialize( U0 );
409 
410  // solve ode
411  BaseType :: solve (U0);
412 
413  // calls setTime ( t + dt );
414  this->tp_.augmentTime();
415 
416  // global min of dt and reset of dtEstimate
417  this->tp_.syncTimeStep();
418 
419  return this->tp_.time();
420  }
421 
422  void printGrid(int nr, const typename Operator::DestinationType& U)
423  {
424  if (time()>=savetime_) {
425  printSGrid(time(),savestep_*10+nr,this->op_.space(),U);
426  ++savestep_;
427  savetime_+=0.001;
428  }
429  }
430 
431  void printmyInfo(string filename) const {
432  std::ostringstream filestream;
433  filestream << filename;
434  std::ofstream ofs(filestream.str().c_str(), std::ios::app);
435  ofs << "ExplMultiStep, steps: " << this->ord_ << "\n\n";
436  ofs << " cfl: " << this->tp_.cfl() << "\\\\\n\n";
437  ofs.close();
438  this->op_.printmyInfo(filename);
439  }
440 
441  private:
442  // TimeProvider with communicator
443  ParallelTimeProvider<DuneCommunicatorType> tp_;
444  double savetime_;
445  int savestep_;
446  };
447 #endif
448 
454  template<class DestinationImp>
455  class ExplicitMultiStepSolver :
456  public OdeSolverInterface<DestinationImp> ,
457  public ExplMultiStepBase<SpaceOperatorInterface<DestinationImp> >
458  {
459  typedef DestinationImp DestinationType;
460  typedef SpaceOperatorInterface<DestinationImp> OperatorType;
461  typedef ExplMultiStepBase<OperatorType> BaseType;
462 
463  private:
464  TimeProviderBase &timeProvider_;
465  double cfl_;
466 
467  public:
474  ExplicitMultiStepSolver(OperatorType& op, TimeProviderBase& tp, int pord, bool verbose = false) :
475  BaseType(op,tp,pord,verbose),
476  timeProvider_(tp),
477  cfl_( std :: min( 0.45 / (2.0 * pord+1) / double(pord), 1.0 ) )
478  {
479  if(verbose)
480  std::cout << "ExplicitMultiStepSolver: cfl = " << cfl_ << "!" << std :: endl;
481  }
482 
484  virtual ~ExplicitMultiStepSolver() {}
485 
487  void initialize(const DestinationType& U0)
488  {
489  BaseType :: initialize(U0);
490  }
491 
493  void solve(DestinationType& U0)
494  {
495  // initialize
496  if( ! this->initialized_ )
497  {
498  DUNE_THROW(InvalidStateException,"ExplicitMultiStepSolver wasn't initialized before first call!");
499  }
500 
501  // solve ode
502  BaseType :: solve(U0);
503  }
504  };
505 
507 } // namespace DuneODE
508 
509 #endif // #ifndef MULTISTEP_SOLVER_HH
DuneODE::ExplMultiStepBase
Definition: multistep.hh:32
DuneODE::ExplMultiStepBase::~ExplMultiStepBase
~ExplMultiStepBase()
destructor
Definition: multistep.hh:120
DuneODE::ExplMultiStepBase::ExplMultiStepBase
ExplMultiStepBase(Operator &op, TimeProviderBase &tp, int pord, bool verbose=true)
constructor
Definition: multistep.hh:57
DuneODE::ExplMultiStepBase::solveRK
void solveRK(DestinationType &U0)
solve the system
Definition: multistep.hh:313
DuneODE::ExplMultiStepBase::initialize
void initialize(const DestinationType &U0)
apply operator once to get dt estimate
Definition: multistep.hh:131
DuneODE::ExplMultiStepBase::solve
void solve(DestinationType &U0)
solve the system
Definition: multistep.hh:142
DuneODE::ExplicitMultiStepSolver
Exlicit multi step ODE solver.
Definition: multistep.hh:458
DuneODE::ExplicitMultiStepSolver::solve
void solve(DestinationType &U0)
solve system
Definition: multistep.hh:493
DuneODE::ExplicitMultiStepSolver::ExplicitMultiStepSolver
ExplicitMultiStepSolver(OperatorType &op, TimeProviderBase &tp, int pord, bool verbose=false)
constructor
Definition: multistep.hh:474
DuneODE::ExplicitMultiStepSolver::initialize
void initialize(const DestinationType &U0)
apply operator once to get dt estimate
Definition: multistep.hh:487
DuneODE::ExplicitMultiStepSolver::~ExplicitMultiStepSolver
virtual ~ExplicitMultiStepSolver()
destructor
Definition: multistep.hh:484
DuneODE::OdeSolverInterface
Interface class for ODE Solver.
Definition: odesolverinterface.hh:21
Dune::Fem::FunctionSpace
A vector valued function space.
Definition: functionspace.hh:60
Dune::Fem::SpaceOperatorInterface
interface for time evolution operators
Definition: spaceoperatorif.hh:38
Dune::Fem::TimeProviderBase
general base for time providers
Definition: timeprovider.hh:36
Dune::InvalidStateException
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
Dune::Operator
An abstract operator Interface class for Operators. Operators are applied to Functions and the result...
Definition: operator.hh:229
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune::Hybrid::min
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
Dune::Simd::Overloads::cond
V cond(ADLTag< 0 >, const Mask< V > &mask, const V &ifTrue, const V &ifFalse)=delete
implements Simd::cond()
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::printGrid
void printGrid(const GridType &grid, const Dune::MPIHelper &helper, std::string output_file="printgrid", int size=2000, bool execute_plot=true, bool png=true, bool local_corner_indices=true, bool local_intersection_indices=true, bool outer_normals=true)
Print a grid as a gnuplot for testing and development.
Definition: printgrid.hh:73
Dune::Fem::Operator
abstract operator
Definition: operator.hh:34
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (May 16, 22:29, 2024)