dune-pdelab  2.5-dev
explicitonestep.hh
Go to the documentation of this file.
1 // -*- tab-width: 2; indent-tabs-mode: nil -*-
2 // vi: set et ts=2 sw=2 sts=2:
3 
4 #ifndef DUNE_PDELAB_INSTATIONARY_EXPLICITONESTEP_HH
5 #define DUNE_PDELAB_INSTATIONARY_EXPLICITONESTEP_HH
6 
7 #include <iostream>
8 #include <iomanip>
9 
10 #include <dune/common/ios_state.hh>
13 
14 namespace Dune {
15  namespace PDELab {
16 
25  template<class R>
27  {
28  public:
29  typedef R RealType;
30 
33  virtual RealType suggestTimestep (RealType time, RealType givendt) = 0;
34 
37  };
38 
40 
43  template<class R>
45  {
46  public:
47  typedef R RealType;
48 
51  virtual RealType suggestTimestep (RealType time, RealType givendt)
52  {
53  return givendt;
54  }
55  };
56 
57 
59 
63  template<class R, class IGOS>
65  {
66  public:
67  typedef R RealType;
68 
69  CFLTimeController (R cfl_, const IGOS& igos_) : cfl(cfl_), target(1e100), igos(igos_)
70  {}
71 
72  CFLTimeController (R cfl_, R target_, const IGOS& igos_) : cfl(cfl_), target(target_), igos(igos_)
73  {}
74 
75  void setTarget (R target_)
76  {
77  target = target_;
78  }
79 
82  virtual RealType suggestTimestep (RealType time, RealType givendt)
83  {
84  RealType suggested = cfl*igos.suggestTimestep(givendt);
85  if (time+2.0*suggested<target)
86  return suggested;
87  if (time+suggested<target)
88  return 0.5*(target-time);
89  return target-time;
90  }
91 
92  private:
93  R cfl;
94  R target;
95  const IGOS& igos;
96  };
97 
98 
100 
108  template<class T, class IGOS, class LS, class TrlV, class TstV = TrlV, class TC = SimpleTimeController<T> >
110  {
111  typedef typename TrlV::ElementType Real;
112  typedef typename IGOS::template MatrixContainer<Real>::Type M;
113 
114  public:
116 
127  ExplicitOneStepMethod(const TimeSteppingParameterInterface<T>& method_, IGOS& igos_, LS& ls_)
128  : method(&method_), igos(igos_), ls(ls_), verbosityLevel(1), step(1), D(igos),
129  tc(new SimpleTimeController<T>()), allocated(true)
130  {
131  if (method->implicit())
132  DUNE_THROW(Exception,"explicit one step method called with implicit scheme");
133  if (igos.trialGridFunctionSpace().gridView().comm().rank()>0)
134  verbosityLevel = 0;
135  }
136 
138 
149  ExplicitOneStepMethod(const TimeSteppingParameterInterface<T>& method_, IGOS& igos_, LS& ls_, TC& tc_)
150  : method(&method_), igos(igos_), ls(ls_), verbosityLevel(1), step(1), D(igos),
151  tc(&tc_), allocated(false)
152  {
153  if (method->implicit())
154  DUNE_THROW(Exception,"explicit one step method called with implicit scheme");
155  }
156 
158  {
159  if (allocated) delete tc;
160  }
161 
163  void setVerbosityLevel (int level)
164  {
165  if (igos.trialGridFunctionSpace().gridView().comm().rank()>0)
166  verbosityLevel = 0;
167  else
168  verbosityLevel = level;
169  }
170 
172  void setStepNumber(int newstep) { step = newstep; }
173 
175 
183  {
184  method = &method_;
185  if (method->implicit())
186  DUNE_THROW(Exception,"explicit one step method called with implicit scheme");
187  }
188 
196  T apply (T time, T dt, TrlV& xold, TrlV& xnew)
197  {
198  DefaultLimiter limiter;
199  return apply(time,dt,xold,xnew,limiter);
200  }
201 
210  template<typename Limiter>
211  T apply (T time, T dt, TrlV& xold, TrlV& xnew, Limiter& limiter)
212  {
213  // save formatting attributes
214  ios_base_all_saver format_attribute_saver(std::cout);
215  LocalTag mytag;
216  mytag << "ExplicitOneStepMethod::apply(): ";
217 
218  std::vector<TrlV*> x(1); // vector of pointers to all steps
219  x[0] = &xold; // initially we have only one
220  if(verbosityLevel>=4)
221  std::cout << mytag << "Creating residual vectors alpha and beta..."
222  << std::endl;
223  TstV alpha(igos.testGridFunctionSpace()), beta(igos.testGridFunctionSpace()); // split residual vectors
224  if(verbosityLevel>=4)
225  std::cout << mytag
226  << "Creating residual vectors alpha and beta... done."
227  << std::endl;
228 
229  if (verbosityLevel>=1){
230  std::ios_base::fmtflags oldflags = std::cout.flags();
231  std::cout << "TIME STEP [" << method->name() << "] "
232  << std::setw(6) << step
233  << " time (from): "
234  << std::setw(12) << std::setprecision(4) << std::scientific
235  << time
236  << " dt: "
237  << std::setw(12) << std::setprecision(4) << std::scientific
238  << dt
239  << " time (to): "
240  << std::setw(12) << std::setprecision(4) << std::scientific
241  << time+dt
242  << std::endl;
243  std::cout.flags(oldflags);
244  }
245 
246  // prepare assembler
247  if(verbosityLevel>=4)
248  std::cout << mytag << "Preparing assembler..." << std::endl;
249  igos.preStep(*method,time,dt);
250  if(verbosityLevel>=4)
251  std::cout << mytag << "Preparing assembler... done." << std::endl;
252 
253  // loop over all stages
254  for(unsigned r=1; r<=method->s(); ++r)
255  {
256  LocalTag stagetag(mytag);
257  stagetag << "stage " << r << ": ";
258  if (verbosityLevel>=4)
259  std::cout << stagetag << "Start." << std::endl;
260 
261  if (verbosityLevel>=2){
262  std::ios_base::fmtflags oldflags = std::cout.flags();
263  std::cout << "STAGE "
264  << r
265  << " time (to): "
266  << std::setw(12) << std::setprecision(4) << std::scientific
267  << time+method->d(r)*dt
268  << "." << std::endl;
269  std::cout.flags(oldflags);
270  }
271 
272  // get vector for current stage
273  if (r==method->s())
274  {
275  // last stage
276  x.push_back(&xnew);
277  if (r>1) xnew = *(x[r-1]); // if r=1 then xnew has already initial guess
278  }
279  else
280  {
281  // intermediate step
282  x.push_back(new TrlV(igos.trialGridFunctionSpace()));
283  if (r>1)
284  *(x[r]) = *(x[r-1]); // use result of last stage as initial guess
285  else
286  *(x[r]) = xnew;
287  }
288 
289  // compute residuals and jacobian
290  if (verbosityLevel>=4) std::cout << "assembling D, alpha, beta ..." << std::endl;
291  D = Real(0.0);
292  alpha = 0.0;
293  beta = 0.0;
294 
295  //apply slope limiter to old solution (e.g for finite volume reconstruction scheme)
296  limiter.prestage(*x[r-1]);
297 
298  if(verbosityLevel>=4)
299  std::cout << stagetag << "Assembling residual..." << std::endl;
300  igos.explicit_jacobian_residual(r,x,D,alpha,beta);
301  if(verbosityLevel>=4)
302  std::cout << stagetag << "Assembling residual... done."
303  << std::endl;
304 
305  // let time controller compute the optimal dt in first stage
306  if (r==1)
307  {
308  T newdt = tc->suggestTimestep(time,dt);
309  newdt = std::min(newdt, dt);
310 
311  if (verbosityLevel>=4){
312  std::ios_base::fmtflags oldflags = std::cout.flags();
313  std::cout << stagetag
314  << "current dt: "
315  << std::setw(12) << std::setprecision(4) << std::scientific
316  << dt
317  << " suggested dt: "
318  << std::setw(12) << std::setprecision(4) << std::scientific
319  << newdt
320  << std::endl;
321  std::cout.flags(oldflags);
322  }
323 
324  if (verbosityLevel>=2 && newdt!=dt)
325  {
326  std::ios_base::fmtflags oldflags = std::cout.flags();
327  std::cout << "changed dt to "
328  << std::setw(12) << std::setprecision(4) << std::scientific
329  << newdt
330  << std::endl;
331  std::cout.flags(oldflags);
332  }
333  dt = newdt;
334  }
335 
336  // combine residual with selected dt
337  if (verbosityLevel>=4)
338  std::cout << stagetag
339  << "Combining residuals with selected dt..."
340  << std::endl;
341  alpha.axpy(dt,beta);
342  if (verbosityLevel>=4)
343  std::cout << stagetag
344  << "Combining residuals with selected dt... done."
345  << std::endl;
346 
347  // solve diagonal system
348  if (verbosityLevel>=4)
349  std::cout << stagetag << "Solving diagonal system..."
350  << std::endl;
351  ls.apply(D,*x[r],alpha,0.99); // dummy reduction
352  if (verbosityLevel>=4)
353  std::cout << stagetag << "Solving diagonal system... done."
354  << std::endl;
355 
356  // apply slope limiter to new solution (e.g DG scheme)
357  limiter.poststage(*x[r]);
358 
359  // stage cleanup
360  if (verbosityLevel>=4)
361  std::cout << stagetag << "Cleanup..." << std::endl;
362  igos.postStage();
363  if (verbosityLevel>=4)
364  std::cout << stagetag << "Cleanup... done" << std::endl;
365 
366  if (verbosityLevel>=4)
367  std::cout << stagetag << "Finished." << std::endl;
368  }
369 
370  // delete intermediate steps
371  for(unsigned i=1; i<method->s(); ++i) delete x[i];
372 
373  // step cleanup
374  if (verbosityLevel>=4)
375  std::cout << mytag << "Cleanup..." << std::endl;
376  igos.postStep();
377  if (verbosityLevel>=4)
378  std::cout << mytag << "Cleanup... done." << std::endl;
379 
380  step++;
381  return dt;
382  }
383 
392  template <typename F>
393  T apply (T time, T dt, TrlV& xold, F& f, TrlV& xnew)
394  {
395  DefaultLimiter limiter;
396  return apply(time,dt,xold,f,xnew,limiter);
397  }
398 
408  template<typename F, typename Limiter>
409  T apply (T time, T dt, TrlV& xold, F&f, TrlV& xnew, Limiter& limiter)
410  {
411  // save formatting attributes
412  ios_base_all_saver format_attribute_saver(std::cout);
413  LocalTag mytag;
414  mytag << "ExplicitOneStepMethod::apply(): ";
415 
416  std::vector<TrlV*> x(1); // vector of pointers to all steps
417  x[0] = &xold; // initially we have only one
418  if(verbosityLevel>=4)
419  std::cout << mytag << "Creating residual vectors alpha and beta..."
420  << std::endl;
421  TstV alpha(igos.testGridFunctionSpace()), beta(igos.testGridFunctionSpace()); // split residual vectors
422  if(verbosityLevel>=4)
423  std::cout << mytag
424  << "Creating residual vectors alpha and beta... done."
425  << std::endl;
426 
427  // In order to apply boundary constraints correctly we need to
428  // solve the linear system from each stage in residual form
429  // with appropriate constraints. 'residual' and 'update' are
430  // two vectors needed for the residual formulation.
431  if(verbosityLevel>=4)
432  std::cout << mytag << "Creating residual vector and update for residual formulation of linear problem per stage"
433  << std::endl;
434  TrlV residual(igos.testGridFunctionSpace());
435  TrlV update(igos.testGridFunctionSpace());
436  if(verbosityLevel>=4)
437  std::cout << mytag << "Creating residual vector and update for residual... done."
438  << std::endl;
439 
440 
441  if (verbosityLevel>=1){
442  std::ios_base::fmtflags oldflags = std::cout.flags();
443  std::cout << "TIME STEP [" << method->name() << "] "
444  << std::setw(6) << step
445  << " time (from): "
446  << std::setw(12) << std::setprecision(4) << std::scientific
447  << time
448  << " dt: "
449  << std::setw(12) << std::setprecision(4) << std::scientific
450  << dt
451  << " time (to): "
452  << std::setw(12) << std::setprecision(4) << std::scientific
453  << time+dt
454  << std::endl;
455  std::cout.flags(oldflags);
456  }
457 
458  // prepare assembler
459  if(verbosityLevel>=4)
460  std::cout << mytag << "Preparing assembler..." << std::endl;
461  igos.preStep(*method,time,dt);
462  if(verbosityLevel>=4)
463  std::cout << mytag << "Preparing assembler... done." << std::endl;
464 
465  // loop over all stages
466  for(unsigned r=1; r<=method->s(); ++r)
467  {
468  LocalTag stagetag(mytag);
469  stagetag << "stage " << r << ": ";
470  if (verbosityLevel>=4)
471  std::cout << stagetag << "Start." << std::endl;
472 
473  if (verbosityLevel>=2){
474  std::ios_base::fmtflags oldflags = std::cout.flags();
475  std::cout << "STAGE "
476  << r
477  << " time (to): "
478  << std::setw(12) << std::setprecision(4) << std::scientific
479  << time+method->d(r)*dt
480  << "." << std::endl;
481  std::cout.flags(oldflags);
482  }
483 
484  // get vector for current stage
485  if (r==method->s())
486  {
487  // last stage
488  x.push_back(&xnew);
489  if (r>1) xnew = *(x[r-1]); // if r=1 then xnew has already initial guess
490  }
491  else
492  {
493  // intermediate step
494  x.push_back(new TrlV(igos.trialGridFunctionSpace()));
495  if (r>1)
496  *(x[r]) = *(x[r-1]); // use result of last stage as initial guess
497  else
498  *(x[r]) = xnew;
499  }
500 
501  // compute residuals and jacobian
502  if (verbosityLevel>=4) std::cout << "assembling D, alpha, beta ..." << std::endl;
503  D = Real(0.0);
504  alpha = 0.0;
505  beta = 0.0;
506 
507  // set x[r] to x[r-1] with boundary conditions interpolated from f
508  igos.interpolate(r,*x[r-1],f,*x[r]);
509 
510  // apply slope limiter to old solution (e.g for finite volume reconstruction scheme)
511  limiter.prestage(*x[r-1]);
512 
513  if(verbosityLevel>=4)
514  std::cout << stagetag << "Assembling residual..." << std::endl;
515  igos.explicit_jacobian_residual(r,x,D,alpha,beta);
516  if(verbosityLevel>=4)
517  std::cout << stagetag << "Assembling residual... done."
518  << std::endl;
519 
520  // let time controller compute the optimal dt in first stage
521  if (r==1)
522  {
523  T newdt = tc->suggestTimestep(time,dt);
524  newdt = std::min(newdt, dt);
525 
526  if (verbosityLevel>=4){
527  std::ios_base::fmtflags oldflags = std::cout.flags();
528  std::cout << stagetag
529  << "current dt: "
530  << std::setw(12) << std::setprecision(4) << std::scientific
531  << dt
532  << " suggested dt: "
533  << std::setw(12) << std::setprecision(4) << std::scientific
534  << newdt
535  << std::endl;
536  std::cout.flags(oldflags);
537  }
538 
539  if (verbosityLevel>=2 && newdt!=dt)
540  {
541  std::ios_base::fmtflags oldflags = std::cout.flags();
542  std::cout << "changed dt to "
543  << std::setw(12) << std::setprecision(4) << std::scientific
544  << newdt
545  << std::endl;
546  std::cout.flags(oldflags);
547  }
548  dt = newdt;
549  }
550 
551  // combine residual with selected dt
552  if (verbosityLevel>=4)
553  std::cout << stagetag
554  << "Combining residuals with selected dt..."
555  << std::endl;
556  alpha.axpy(dt,beta);
557  if (verbosityLevel>=4)
558  std::cout << stagetag
559  << "Combining residuals with selected dt... done."
560  << std::endl;
561 
562 
563 
564  // Set up residual formulation (for Dx[r]=alpha) and
565  // compute update by solving diagonal system
566  using Backend::native;
567  native(D).mv(native(*x[r]), native(residual));
568  residual -= alpha;
569  auto cc = igos.trialConstraints();
570  Dune::PDELab::set_constrained_dofs(cc, 0.0, residual);
571  if (verbosityLevel>=4)
572  std::cout << stagetag << "Solving diagonal system..."
573  << std::endl;
574  ls.apply(D, update, residual, 0.99); // dummy reduction
575  if (verbosityLevel>=4)
576  std::cout << stagetag << "Solving diagonal system... done."
577  << std::endl;
578  *x[r] -= update;
579 
580  // apply slope limiter to new solution (e.g DG scheme)
581  limiter.poststage(*x[r]);
582 
583  // stage cleanup
584  if (verbosityLevel>=4)
585  std::cout << stagetag << "Cleanup..." << std::endl;
586  igos.postStage();
587  if (verbosityLevel>=4)
588  std::cout << stagetag << "Cleanup... done" << std::endl;
589 
590  if (verbosityLevel>=4)
591  std::cout << stagetag << "Finished." << std::endl;
592  }
593 
594  // delete intermediate steps
595  for(unsigned i=1; i<method->s(); ++i) delete x[i];
596 
597  // step cleanup
598  if (verbosityLevel>=4)
599  std::cout << mytag << "Cleanup..." << std::endl;
600  igos.postStep();
601  if (verbosityLevel>=4)
602  std::cout << mytag << "Cleanup... done." << std::endl;
603 
604  step++;
605  return dt;
606  }
607 
608  private:
609 
611  class DefaultLimiter
612  {
613  public:
614  template<typename V>
615  void prestage(V& v)
616  {}
617 
618  template<typename V>
619  void poststage(V& v)
620  {}
621  };
622 
623  const TimeSteppingParameterInterface<T> *method;
624  IGOS& igos;
625  LS& ls;
626  int verbosityLevel;
627  int step;
628  M D;
630  bool allocated;
631  };
632 
634  } // end namespace PDELab
635 } // end namespace Dune
636 #endif // DUNE_PDELAB_INSTATIONARY_EXPLICITONESTEP_HH
R RealType
Definition: explicitonestep.hh:67
Default time controller; just returns given dt.
Definition: explicitonestep.hh:44
R RealType
Definition: explicitonestep.hh:47
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
Do one step of an explicit time-stepping scheme.
Definition: explicitonestep.hh:109
T apply(T time, T dt, TrlV &xold, TrlV &xnew)
do one step;
Definition: explicitonestep.hh:196
~ExplicitOneStepMethod()
Definition: explicitonestep.hh:157
Insert standard boilerplate into log messages.
Definition: logtag.hh:180
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:798
ExplicitOneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, LS &ls_, TC &tc_)
construct a new one step scheme
Definition: explicitonestep.hh:149
T apply(T time, T dt, TrlV &xold, F &f, TrlV &xnew)
do one step;
Definition: explicitonestep.hh:393
void setMethod(const TimeSteppingParameterInterface< T > &method_)
redefine the method to be used; can be done before every step
Definition: explicitonestep.hh:182
CFLTimeController(R cfl_, R target_, const IGOS &igos_)
Definition: explicitonestep.hh:72
virtual RealType suggestTimestep(RealType time, RealType givendt)
Return name of the scheme.
Definition: explicitonestep.hh:51
void setStepNumber(int newstep)
change number of current step
Definition: explicitonestep.hh:172
Controller interface for adaptive time stepping.
Definition: explicitonestep.hh:26
ExplicitOneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, LS &ls_)
construct a new one step scheme
Definition: explicitonestep.hh:127
limit time step to maximum dt * CFL number
Definition: explicitonestep.hh:64
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
void setVerbosityLevel(int level)
change verbosity level; 0 means completely quiet
Definition: explicitonestep.hh:163
T apply(T time, T dt, TrlV &xold, F &f, TrlV &xnew, Limiter &limiter)
do one step;
Definition: explicitonestep.hh:409
virtual ~TimeControllerInterface()
every abstract base class has a virtual destructor
Definition: explicitonestep.hh:36
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > &>::type native(T &t)
Definition: backend/interface.hh:192
void setTarget(R target_)
Definition: explicitonestep.hh:75
virtual RealType suggestTimestep(RealType time, RealType givendt)=0
Return name of the scheme.
virtual RealType suggestTimestep(RealType time, RealType givendt)
Return name of the scheme.
Definition: explicitonestep.hh:82
T apply(T time, T dt, TrlV &xold, TrlV &xnew, Limiter &limiter)
do one step;
Definition: explicitonestep.hh:211
R RealType
Definition: explicitonestep.hh:29
CFLTimeController(R cfl_, const IGOS &igos_)
Definition: explicitonestep.hh:69