DUNE PDELab (2.7)

explicitonestep.hh
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>
11#include <dune/pdelab/common/logtag.hh>
12#include <dune/pdelab/instationary/onestepparameter.hh>
13
14namespace Dune {
15 namespace PDELab {
16
25 template<class R>
26 class TimeControllerInterface
27 {
28 public:
29 typedef R RealType;
30
33 virtual RealType suggestTimestep (RealType time, RealType givendt) = 0;
34
36 virtual ~TimeControllerInterface () {}
37 };
38
40
43 template<class R>
44 class SimpleTimeController : public TimeControllerInterface<R>
45 {
46 public:
47 typedef R RealType;
48
51 virtual RealType suggestTimestep (RealType time, RealType givendt) override
52 {
53 return givendt;
54 }
55 };
56
57
59
63 template<class R, class IGOS>
64 class CFLTimeController : public TimeControllerInterface<R>
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) override
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> >
109 class ExplicitOneStepMethod
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
157 ~ExplicitOneStepMethod ()
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
182 void setMethod (const TimeSteppingParameterInterface<T>& method_)
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;
629 TimeControllerInterface<T> *tc;
630 bool allocated;
631 };
632
634 } // end namespace PDELab
635} // end namespace Dune
636#endif // DUNE_PDELAB_INSTATIONARY_EXPLICITONESTEP_HH
Dune::PDELab::CFLTimeController
limit time step to maximum dt * CFL number
Definition: explicitonestep.hh:65
Dune::PDELab::CFLTimeController::suggestTimestep
virtual RealType suggestTimestep(RealType time, RealType givendt) override
Return name of the scheme.
Definition: explicitonestep.hh:82
Dune::PDELab::Exception
Base class for all PDELab exceptions.
Definition: exceptions.hh:19
Dune::PDELab::ExplicitOneStepMethod
Do one step of an explicit time-stepping scheme.
Definition: explicitonestep.hh:110
Dune::PDELab::ExplicitOneStepMethod::apply
T apply(T time, T dt, TrlV &xold, F &f, TrlV &xnew, Limiter &limiter)
do one step;
Definition: explicitonestep.hh:409
Dune::PDELab::ExplicitOneStepMethod::ExplicitOneStepMethod
ExplicitOneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, LS &ls_, TC &tc_)
construct a new one step scheme
Definition: explicitonestep.hh:149
Dune::PDELab::ExplicitOneStepMethod::apply
T apply(T time, T dt, TrlV &xold, F &f, TrlV &xnew)
do one step;
Definition: explicitonestep.hh:393
Dune::PDELab::ExplicitOneStepMethod::apply
T apply(T time, T dt, TrlV &xold, TrlV &xnew)
do one step;
Definition: explicitonestep.hh:196
Dune::PDELab::ExplicitOneStepMethod::apply
T apply(T time, T dt, TrlV &xold, TrlV &xnew, Limiter &limiter)
do one step;
Definition: explicitonestep.hh:211
Dune::PDELab::ExplicitOneStepMethod::setMethod
void setMethod(const TimeSteppingParameterInterface< T > &method_)
redefine the method to be used; can be done before every step
Definition: explicitonestep.hh:182
Dune::PDELab::ExplicitOneStepMethod::setVerbosityLevel
void setVerbosityLevel(int level)
change verbosity level; 0 means completely quiet
Definition: explicitonestep.hh:163
Dune::PDELab::ExplicitOneStepMethod::setStepNumber
void setStepNumber(int newstep)
change number of current step
Definition: explicitonestep.hh:172
Dune::PDELab::ExplicitOneStepMethod::ExplicitOneStepMethod
ExplicitOneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, LS &ls_)
construct a new one step scheme
Definition: explicitonestep.hh:127
Dune::PDELab::LocalTag
Insert standard boilerplate into log messages.
Definition: logtag.hh:172
Dune::PDELab::SimpleTimeController
Default time controller; just returns given dt.
Definition: explicitonestep.hh:45
Dune::PDELab::SimpleTimeController::suggestTimestep
virtual RealType suggestTimestep(RealType time, RealType givendt) override
Return name of the scheme.
Definition: explicitonestep.hh:51
Dune::PDELab::TimeControllerInterface
Controller interface for adaptive time stepping.
Definition: explicitonestep.hh:27
Dune::PDELab::TimeControllerInterface::~TimeControllerInterface
virtual ~TimeControllerInterface()
every abstract base class has a virtual destructor
Definition: explicitonestep.hh:36
Dune::PDELab::TimeControllerInterface::suggestTimestep
virtual RealType suggestTimestep(RealType time, RealType givendt)=0
Return name of the scheme.
Dune::ios_base_all_saver
Utility class for storing and resetting stream attributes.
Definition: ios_state.hh:32
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::PDELab::set_constrained_dofs
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:796
Dune::PDELab::TimeSteppingParameterInterface::implicit
virtual bool implicit() const =0
Return true if method is implicit.
Dune::PDELab::TimeSteppingParameterInterface::name
virtual std::string name() const =0
Return name of the scheme.
Dune::PDELab::TimeSteppingParameterInterface::d
virtual R d(int r) const =0
Return entries of the d Vector.
Dune::PDELab::TimeSteppingParameterInterface::s
virtual unsigned s() const =0
Return number of stages of the method.
Dune::Simd::Overloads::min
auto min(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::min()
Definition: defaults.hh:87
ios_state.hh
Utility class for storing and resetting stream attributes.
Dune
Dune namespace.
Definition: alignedallocator.hh:14
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)