dune-pdelab  2.4-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 
197  T apply (T time, T dt, TrlV& xold, TrlV& xnew)
198  {
199  DefaultLimiter limiter;
200  return apply(time,dt,xold,xnew,limiter);
201  }
202 
203  template<typename Limiter>
204  T apply (T time, T dt, TrlV& xold, TrlV& xnew, Limiter& limiter)
205  {
206  // save formatting attributes
207  ios_base_all_saver format_attribute_saver(std::cout);
208  LocalTag mytag;
209  mytag << "ExplicitOneStepMethod::apply(): ";
210 
211  std::vector<TrlV*> x(1); // vector of pointers to all steps
212  x[0] = &xold; // initially we have only one
213  if(verbosityLevel>=4)
214  std::cout << mytag << "Creating residual vectors alpha and beta..."
215  << std::endl;
216  TstV alpha(igos.testGridFunctionSpace()), beta(igos.testGridFunctionSpace()); // split residual vectors
217  if(verbosityLevel>=4)
218  std::cout << mytag
219  << "Creating residual vectors alpha and beta... done."
220  << std::endl;
221 
222  if (verbosityLevel>=1){
223  std::ios_base::fmtflags oldflags = std::cout.flags();
224  std::cout << "TIME STEP [" << method->name() << "] "
225  << std::setw(6) << step
226  << " time (from): "
227  << std::setw(12) << std::setprecision(4) << std::scientific
228  << time
229  << " dt: "
230  << std::setw(12) << std::setprecision(4) << std::scientific
231  << dt
232  << " time (to): "
233  << std::setw(12) << std::setprecision(4) << std::scientific
234  << time+dt
235  << std::endl;
236  std::cout.flags(oldflags);
237  }
238 
239  // prepare assembler
240  if(verbosityLevel>=4)
241  std::cout << mytag << "Preparing assembler..." << std::endl;
242  igos.preStep(*method,time,dt);
243  if(verbosityLevel>=4)
244  std::cout << mytag << "Preparing assembler... done." << std::endl;
245 
246  // loop over all stages
247  for(unsigned r=1; r<=method->s(); ++r)
248  {
249  LocalTag stagetag(mytag);
250  stagetag << "stage " << r << ": ";
251  if (verbosityLevel>=4)
252  std::cout << stagetag << "Start." << std::endl;
253 
254  if (verbosityLevel>=2){
255  std::ios_base::fmtflags oldflags = std::cout.flags();
256  std::cout << "STAGE "
257  << r
258  << " time (to): "
259  << std::setw(12) << std::setprecision(4) << std::scientific
260  << time+method->d(r)*dt
261  << "." << std::endl;
262  std::cout.flags(oldflags);
263  }
264 
265  // get vector for current stage
266  if (r==method->s())
267  {
268  // last stage
269  x.push_back(&xnew);
270  if (r>1) xnew = *(x[r-1]); // if r=1 then xnew has already initial guess
271  }
272  else
273  {
274  // intermediate step
275  x.push_back(new TrlV(igos.trialGridFunctionSpace()));
276  if (r>1)
277  *(x[r]) = *(x[r-1]); // use result of last stage as initial guess
278  else
279  *(x[r]) = xnew;
280  }
281 
282  // compute residuals and jacobian
283  if (verbosityLevel>=4) std::cout << "assembling D, alpha, beta ..." << std::endl;
284  D = Real(0.0);
285  alpha = 0.0;
286  beta = 0.0;
287 
288  //apply slope limiter to old solution (e.g for finite volume reconstruction scheme)
289  limiter.prestage(*x[r-1]);
290 
291  if(verbosityLevel>=4)
292  std::cout << stagetag << "Assembling residual..." << std::endl;
293  igos.explicit_jacobian_residual(r,x,D,alpha,beta);
294  if(verbosityLevel>=4)
295  std::cout << stagetag << "Assembling residual... done."
296  << std::endl;
297 
298  // let time controller compute the optimal dt in first stage
299  if (r==1)
300  {
301  T newdt = tc->suggestTimestep(time,dt);
302  newdt = std::min(newdt, dt);
303 
304  if (verbosityLevel>=4){
305  std::ios_base::fmtflags oldflags = std::cout.flags();
306  std::cout << stagetag
307  << "current dt: "
308  << std::setw(12) << std::setprecision(4) << std::scientific
309  << dt
310  << " suggested dt: "
311  << std::setw(12) << std::setprecision(4) << std::scientific
312  << newdt
313  << std::endl;
314  std::cout.flags(oldflags);
315  }
316 
317  if (verbosityLevel>=2 && newdt!=dt)
318  {
319  std::ios_base::fmtflags oldflags = std::cout.flags();
320  std::cout << "changed dt to "
321  << std::setw(12) << std::setprecision(4) << std::scientific
322  << newdt
323  << std::endl;
324  std::cout.flags(oldflags);
325  }
326  dt = newdt;
327  }
328 
329  // combine residual with selected dt
330  if (verbosityLevel>=4)
331  std::cout << stagetag
332  << "Combining residuals with selected dt..."
333  << std::endl;
334  alpha.axpy(dt,beta);
335  if (verbosityLevel>=4)
336  std::cout << stagetag
337  << "Combining residuals with selected dt... done."
338  << std::endl;
339 
340  // solve diagonal system
341  if (verbosityLevel>=4)
342  std::cout << stagetag << "Solving diagonal system..."
343  << std::endl;
344  ls.apply(D,*x[r],alpha,0.99); // dummy reduction
345  if (verbosityLevel>=4)
346  std::cout << stagetag << "Solving diagonal system... done."
347  << std::endl;
348 
349  // apply slope limiter to new solution (e.g DG scheme)
350  limiter.poststage(*x[r]);
351 
352  // stage cleanup
353  if (verbosityLevel>=4)
354  std::cout << stagetag << "Cleanup..." << std::endl;
355  igos.postStage();
356  if (verbosityLevel>=4)
357  std::cout << stagetag << "Cleanup... done" << std::endl;
358 
359  if (verbosityLevel>=4)
360  std::cout << stagetag << "Finished." << std::endl;
361  }
362 
363  // delete intermediate steps
364  for(unsigned i=1; i<method->s(); ++i) delete x[i];
365 
366  // step cleanup
367  if (verbosityLevel>=4)
368  std::cout << mytag << "Cleanup..." << std::endl;
369  igos.postStep();
370  if (verbosityLevel>=4)
371  std::cout << mytag << "Cleanup... done." << std::endl;
372 
373  step++;
374  return dt;
375  }
376 
377  private:
378 
380  class DefaultLimiter
381  {
382  public:
383  template<typename V>
384  void prestage(V& v)
385  {}
386 
387  template<typename V>
388  void poststage(V& v)
389  {}
390  };
391 
392  const TimeSteppingParameterInterface<T> *method;
393  IGOS& igos;
394  LS& ls;
395  int verbosityLevel;
396  int step;
397  M D;
398  TimeControllerInterface<T> *tc;
399  bool allocated;
400  };
401 
403  } // end namespace PDELab
404 } // end namespace Dune
405 #endif
CFLTimeController(R cfl_, const IGOS &igos_)
Definition: explicitonestep.hh:69
Default time controller; just returns given dt.
Definition: explicitonestep.hh:44
ExplicitOneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, LS &ls_, TC &tc_)
construct a new one step scheme
Definition: explicitonestep.hh:149
virtual std::string name() const =0
Return name of the scheme.
virtual RealType suggestTimestep(RealType time, RealType givendt)
Return name of the scheme.
Definition: explicitonestep.hh:82
void setTarget(R target_)
Definition: explicitonestep.hh:75
virtual unsigned s() const =0
Return number of stages of the method.
R RealType
Definition: explicitonestep.hh:67
limit time step to maximum dt * CFL number
Definition: explicitonestep.hh:64
R RealType
Definition: explicitonestep.hh:29
R RealType
Definition: explicitonestep.hh:47
void setMethod(const TimeSteppingParameterInterface< T > &method_)
redefine the method to be used; can be done before every step
Definition: explicitonestep.hh:182
virtual ~TimeControllerInterface()
every abstract base class has a virtual destructor
Definition: explicitonestep.hh:36
Do one step of an explicit time-stepping scheme.
Definition: explicitonestep.hh:109
ExplicitOneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, LS &ls_)
construct a new one step scheme
Definition: explicitonestep.hh:127
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
virtual RealType suggestTimestep(RealType time, RealType givendt)=0
Return name of the scheme.
void setVerbosityLevel(int level)
change verbosity level; 0 means completely quiet
Definition: explicitonestep.hh:163
Definition: adaptivity.hh:27
void setStepNumber(int newstep)
change number of current step
Definition: explicitonestep.hh:172
virtual bool implicit() const =0
Return true if method is implicit.
Controller interface for adaptive time stepping.
Definition: explicitonestep.hh:26
virtual RealType suggestTimestep(RealType time, RealType givendt)
Return name of the scheme.
Definition: explicitonestep.hh:51
CFLTimeController(R cfl_, R target_, const IGOS &igos_)
Definition: explicitonestep.hh:72
virtual R d(int r) const =0
Return entries of the d Vector.
T apply(T time, T dt, TrlV &xold, TrlV &xnew, Limiter &limiter)
Definition: explicitonestep.hh:204
~ExplicitOneStepMethod()
Definition: explicitonestep.hh:157
T apply(T time, T dt, TrlV &xold, TrlV &xnew)
do one step;
Definition: explicitonestep.hh:197
Insert standard boilerplate into log messages.
Definition: logtag.hh:180