DUNE PDELab (git)

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
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>
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) override
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) 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> >
110 {
111 typedef typename TrlV::ElementType Real;
112 typedef typename IGOS::template MatrixContainer<Real>::Type M;
113
114 public:
116
127 [[deprecated(
128 "Constructor with default assigned reduction is no longer supported. "
129 "Use the new constructor by passing a suitable reduction to the linear "
130 "system. In particular, use a dummy reduction (i.e. 0.99) on systems "
131 "known to have diagonal mass matrices (e.g. FV or DG) and where the "
132 "linear solver performs a block inversion on the first iteration "
133 "(essentially any). Otherwise, specify a proper reduction suitable "
134 "for your needs."
135 )]]
136 ExplicitOneStepMethod(const TimeSteppingParameterInterface<T>& method_, IGOS& igos_, LS& ls_)
137 : method(&method_), igos(igos_), ls(ls_), verbosityLevel(1), step(1), D(igos),
138 tc(new SimpleTimeController<T>()), allocated(true), ls_reduct{0.99}
139 {
140 if (method->implicit())
141 DUNE_THROW(Exception,"explicit one step method called with implicit scheme");
142 if (igos.trialGridFunctionSpace().gridView().comm().rank()>0)
143 verbosityLevel = 0;
144 }
145
159 IGOS& igos_,
160 LS& ls_,
161 double ls_reduct_)
162 : method(&method_)
163 , igos(igos_)
164 , ls(ls_)
165 , verbosityLevel(1)
166 , step(1)
167 , D(igos)
168 , tc(new SimpleTimeController<T>())
169 , allocated(true)
170 , ls_reduct{ ls_reduct_ }
171 {
172 if (method->implicit())
174 "explicit one step method called with implicit scheme");
175 if (igos.trialGridFunctionSpace().gridView().comm().rank() > 0)
176 verbosityLevel = 0;
177 }
178
180
191 [[deprecated(
192 "Constructor with default assigned reduction is no longer supported. "
193 "Use the new constructor by passing a suitable reduction to the linear "
194 "system. In particular, use a dummy reduction (i.e. 0.99) on systems \n"
195 "known to have diagonal mass matrices (e.g. FV or DG) and where the \n"
196 "linear solver performs a block inversion on the first iteration \n"
197 "(essentially any). Otherwise, specify a proper reduction suitable \n"
198 "for your needs.\n"
199 )]]
200 ExplicitOneStepMethod(const TimeSteppingParameterInterface<T>& method_, IGOS& igos_, LS& ls_, TC& tc_)
201 : method(&method_), igos(igos_), ls(ls_), verbosityLevel(1), step(1), D(igos),
202 tc(&tc_), allocated(false), ls_reduct{0.99}
203 {
204 if (method->implicit())
205 DUNE_THROW(Exception,"explicit one step method called with implicit scheme");
206 }
207
221 IGOS& igos_,
222 LS& ls_,
223 TC& tc_,
224 double ls_reduct_)
225 : method(&method_)
226 , igos(igos_)
227 , ls(ls_)
228 , verbosityLevel(1)
229 , step(1)
230 , D(igos)
231 , tc(&tc_)
232 , allocated(false)
233 , ls_reduct{ ls_reduct_ }
234 {
235 if (method->implicit())
237 "explicit one step method called with implicit scheme");
238 }
239
241 {
242 if (allocated) delete tc;
243 }
244
246 void setVerbosityLevel (int level)
247 {
248 if (igos.trialGridFunctionSpace().gridView().comm().rank()>0)
249 verbosityLevel = 0;
250 else
251 verbosityLevel = level;
252 }
253
255 void setStepNumber(int newstep) { step = newstep; }
256
268 void setReduction(const double& ls_reduction_) { ls_reduct = ls_reduction_; }
269
271
279 {
280 method = &method_;
281 if (method->implicit())
282 DUNE_THROW(Exception,"explicit one step method called with implicit scheme");
283 }
284
292 T apply (T time, T dt, TrlV& xold, TrlV& xnew)
293 {
294 DefaultLimiter limiter;
295 return apply(time,dt,xold,xnew,limiter);
296 }
297
306 template<typename Limiter>
307 T apply (T time, T dt, TrlV& xold, TrlV& xnew, Limiter& limiter)
308 {
309 // save formatting attributes
310 ios_base_all_saver format_attribute_saver(std::cout);
311 LocalTag mytag;
312 mytag << "ExplicitOneStepMethod::apply(): ";
313
314 std::vector<TrlV*> x(1); // vector of pointers to all steps
315 x[0] = &xold; // initially we have only one
316 if(verbosityLevel>=4)
317 std::cout << mytag << "Creating residual vectors alpha and beta..."
318 << std::endl;
319 TstV alpha(igos.testGridFunctionSpace()), beta(igos.testGridFunctionSpace()); // split residual vectors
320 if(verbosityLevel>=4)
321 std::cout << mytag
322 << "Creating residual vectors alpha and beta... done."
323 << std::endl;
324
325 if (verbosityLevel>=1){
326 std::ios_base::fmtflags oldflags = std::cout.flags();
327 std::cout << "TIME STEP [" << method->name() << "] "
328 << std::setw(6) << step
329 << " time (from): "
330 << std::setw(12) << std::setprecision(4) << std::scientific
331 << time
332 << " dt: "
333 << std::setw(12) << std::setprecision(4) << std::scientific
334 << dt
335 << " time (to): "
336 << std::setw(12) << std::setprecision(4) << std::scientific
337 << time+dt
338 << std::endl;
339 std::cout.flags(oldflags);
340 }
341
342 // prepare assembler
343 if(verbosityLevel>=4)
344 std::cout << mytag << "Preparing assembler..." << std::endl;
345 igos.preStep(*method,time,dt);
346 if(verbosityLevel>=4)
347 std::cout << mytag << "Preparing assembler... done." << std::endl;
348
349 // loop over all stages
350 for(unsigned r=1; r<=method->s(); ++r)
351 {
352 LocalTag stagetag(mytag);
353 stagetag << "stage " << r << ": ";
354 if (verbosityLevel>=4)
355 std::cout << stagetag << "Start." << std::endl;
356
357 if (verbosityLevel>=2){
358 std::ios_base::fmtflags oldflags = std::cout.flags();
359 std::cout << "STAGE "
360 << r
361 << " time (to): "
362 << std::setw(12) << std::setprecision(4) << std::scientific
363 << time+method->d(r)*dt
364 << "." << std::endl;
365 std::cout.flags(oldflags);
366 }
367
368 // get vector for current stage
369 if (r==method->s())
370 {
371 // last stage
372 x.push_back(&xnew);
373 if (r>1) xnew = *(x[r-1]); // if r=1 then xnew has already initial guess
374 }
375 else
376 {
377 // intermediate step
378 x.push_back(new TrlV(igos.trialGridFunctionSpace()));
379 if (r>1)
380 *(x[r]) = *(x[r-1]); // use result of last stage as initial guess
381 else
382 *(x[r]) = xnew;
383 }
384
385 // compute residuals and jacobian
386 if (verbosityLevel>=4) std::cout << "assembling D, alpha, beta ..." << std::endl;
387 D = Real(0.0);
388 alpha = 0.0;
389 beta = 0.0;
390
391 //apply slope limiter to old solution (e.g for finite volume reconstruction scheme)
392 limiter.prestage(*x[r-1]);
393
394 if(verbosityLevel>=4)
395 std::cout << stagetag << "Assembling residual..." << std::endl;
396 igos.explicit_jacobian_residual(r,x,D,alpha,beta);
397 if(verbosityLevel>=4)
398 std::cout << stagetag << "Assembling residual... done."
399 << std::endl;
400
401 // let time controller compute the optimal dt in first stage
402 if (r==1)
403 {
404 T newdt = tc->suggestTimestep(time,dt);
405 newdt = std::min(newdt, dt);
406
407 if (verbosityLevel>=4){
408 std::ios_base::fmtflags oldflags = std::cout.flags();
409 std::cout << stagetag
410 << "current dt: "
411 << std::setw(12) << std::setprecision(4) << std::scientific
412 << dt
413 << " suggested dt: "
414 << std::setw(12) << std::setprecision(4) << std::scientific
415 << newdt
416 << std::endl;
417 std::cout.flags(oldflags);
418 }
419
420 if (verbosityLevel>=2 && newdt!=dt)
421 {
422 std::ios_base::fmtflags oldflags = std::cout.flags();
423 std::cout << "changed dt to "
424 << std::setw(12) << std::setprecision(4) << std::scientific
425 << newdt
426 << std::endl;
427 std::cout.flags(oldflags);
428 }
429 dt = newdt;
430 }
431
432 // combine residual with selected dt
433 if (verbosityLevel>=4)
434 std::cout << stagetag
435 << "Combining residuals with selected dt..."
436 << std::endl;
437 alpha.axpy(dt,beta);
438 if (verbosityLevel>=4)
439 std::cout << stagetag
440 << "Combining residuals with selected dt... done."
441 << std::endl;
442
443 // solve linear system
444 if (verbosityLevel>=4)
445 std::cout << stagetag << "Solving linear system..."
446 << std::endl;
447 ls.apply(D,*x[r],alpha,ls_reduct);
448 if (verbosityLevel>=4)
449 std::cout << stagetag << "Solving linear system... done."
450 << std::endl;
451
452 // apply slope limiter to new solution (e.g DG scheme)
453 limiter.poststage(*x[r]);
454
455 // stage cleanup
456 if (verbosityLevel>=4)
457 std::cout << stagetag << "Cleanup..." << std::endl;
458 igos.postStage();
459 if (verbosityLevel>=4)
460 std::cout << stagetag << "Cleanup... done" << std::endl;
461
462 if (verbosityLevel>=4)
463 std::cout << stagetag << "Finished." << std::endl;
464 }
465
466 // delete intermediate steps
467 for(unsigned i=1; i<method->s(); ++i) delete x[i];
468
469 // step cleanup
470 if (verbosityLevel>=4)
471 std::cout << mytag << "Cleanup..." << std::endl;
472 igos.postStep();
473 if (verbosityLevel>=4)
474 std::cout << mytag << "Cleanup... done." << std::endl;
475
476 step++;
477 return dt;
478 }
479
488 template <typename F>
489 T apply (T time, T dt, TrlV& xold, F& f, TrlV& xnew)
490 {
491 DefaultLimiter limiter;
492 return apply(time,dt,xold,f,xnew,limiter);
493 }
494
504 template<typename F, typename Limiter>
505 T apply (T time, T dt, TrlV& xold, F&f, TrlV& xnew, Limiter& limiter)
506 {
507 // save formatting attributes
508 ios_base_all_saver format_attribute_saver(std::cout);
509 LocalTag mytag;
510 mytag << "ExplicitOneStepMethod::apply(): ";
511
512 std::vector<TrlV*> x(1); // vector of pointers to all steps
513 x[0] = &xold; // initially we have only one
514 if(verbosityLevel>=4)
515 std::cout << mytag << "Creating residual vectors alpha and beta..."
516 << std::endl;
517 TstV alpha(igos.testGridFunctionSpace()), beta(igos.testGridFunctionSpace()); // split residual vectors
518 if(verbosityLevel>=4)
519 std::cout << mytag
520 << "Creating residual vectors alpha and beta... done."
521 << std::endl;
522
523 // In order to apply boundary constraints correctly we need to
524 // solve the linear system from each stage in residual form
525 // with appropriate constraints. 'residual' and 'update' are
526 // two vectors needed for the residual formulation.
527 if(verbosityLevel>=4)
528 std::cout << mytag << "Creating residual vector and update for residual formulation of linear problem per stage"
529 << std::endl;
530 TrlV residual(igos.testGridFunctionSpace());
531 TrlV update(igos.testGridFunctionSpace());
532 if(verbosityLevel>=4)
533 std::cout << mytag << "Creating residual vector and update for residual... done."
534 << std::endl;
535
536
537 if (verbosityLevel>=1){
538 std::ios_base::fmtflags oldflags = std::cout.flags();
539 std::cout << "TIME STEP [" << method->name() << "] "
540 << std::setw(6) << step
541 << " time (from): "
542 << std::setw(12) << std::setprecision(4) << std::scientific
543 << time
544 << " dt: "
545 << std::setw(12) << std::setprecision(4) << std::scientific
546 << dt
547 << " time (to): "
548 << std::setw(12) << std::setprecision(4) << std::scientific
549 << time+dt
550 << std::endl;
551 std::cout.flags(oldflags);
552 }
553
554 // prepare assembler
555 if(verbosityLevel>=4)
556 std::cout << mytag << "Preparing assembler..." << std::endl;
557 igos.preStep(*method,time,dt);
558 if(verbosityLevel>=4)
559 std::cout << mytag << "Preparing assembler... done." << std::endl;
560
561 // loop over all stages
562 for(unsigned r=1; r<=method->s(); ++r)
563 {
564 LocalTag stagetag(mytag);
565 stagetag << "stage " << r << ": ";
566 if (verbosityLevel>=4)
567 std::cout << stagetag << "Start." << std::endl;
568
569 if (verbosityLevel>=2){
570 std::ios_base::fmtflags oldflags = std::cout.flags();
571 std::cout << "STAGE "
572 << r
573 << " time (to): "
574 << std::setw(12) << std::setprecision(4) << std::scientific
575 << time+method->d(r)*dt
576 << "." << std::endl;
577 std::cout.flags(oldflags);
578 }
579
580 // get vector for current stage
581 if (r==method->s())
582 {
583 // last stage
584 x.push_back(&xnew);
585 }
586 else
587 {
588 // intermediate step
589 x.push_back(new TrlV(igos.trialGridFunctionSpace()));
590 }
591
592 // compute residuals and jacobian
593 if (verbosityLevel>=4) std::cout << "assembling D, alpha, beta ..." << std::endl;
594 D = Real(0.0);
595 alpha = 0.0;
596 beta = 0.0;
597
598 // apply slope limiter to old solution (e.g for finite volume reconstruction scheme)
599 limiter.prestage(*x[r-1]);
600
601 // set initial value from xnew or last stage
602 TrlV const * const init_guess = (r==1) ? &xnew : x[r-1];
603 // set boundary conditions
604 igos.interpolate(r,*init_guess,f,*x[r]);
605
606 if(verbosityLevel>=4)
607 std::cout << stagetag << "Assembling residual..." << std::endl;
608 igos.explicit_jacobian_residual(r,x,D,alpha,beta);
609 if(verbosityLevel>=4)
610 std::cout << stagetag << "Assembling residual... done."
611 << std::endl;
612
613 // let time controller compute the optimal dt in first stage
614 if (r==1)
615 {
616 T newdt = tc->suggestTimestep(time,dt);
617 newdt = std::min(newdt, dt);
618
619 if (verbosityLevel>=4){
620 std::ios_base::fmtflags oldflags = std::cout.flags();
621 std::cout << stagetag
622 << "current dt: "
623 << std::setw(12) << std::setprecision(4) << std::scientific
624 << dt
625 << " suggested dt: "
626 << std::setw(12) << std::setprecision(4) << std::scientific
627 << newdt
628 << std::endl;
629 std::cout.flags(oldflags);
630 }
631
632 if (verbosityLevel>=2 && newdt!=dt)
633 {
634 std::ios_base::fmtflags oldflags = std::cout.flags();
635 std::cout << "changed dt to "
636 << std::setw(12) << std::setprecision(4) << std::scientific
637 << newdt
638 << std::endl;
639 std::cout.flags(oldflags);
640 }
641 dt = newdt;
642 }
643
644 // combine residual with selected dt
645 if (verbosityLevel>=4)
646 std::cout << stagetag
647 << "Combining residuals with selected dt..."
648 << std::endl;
649 alpha.axpy(dt,beta);
650 if (verbosityLevel>=4)
651 std::cout << stagetag
652 << "Combining residuals with selected dt... done."
653 << std::endl;
654
655
656
657 // Set up residual formulation (for Dx[r]=alpha) and
658 // compute update by solving linear system
659 using Backend::native;
660 native(D).mv(native(*x[r]), native(residual));
661 residual -= alpha;
662 auto cc = igos.trialConstraints();
663 Dune::PDELab::set_constrained_dofs(cc, 0.0, residual);
664 if (verbosityLevel>=4)
665 std::cout << stagetag << "Solving linear system..."
666 << std::endl;
667 ls.apply(D, update, residual, ls_reduct);
668 if (verbosityLevel>=4)
669 std::cout << stagetag << "Solving linear system... done."
670 << std::endl;
671 *x[r] -= update;
672
673 // apply slope limiter to new solution (e.g DG scheme)
674 limiter.poststage(*x[r]);
675
676 // stage cleanup
677 if (verbosityLevel>=4)
678 std::cout << stagetag << "Cleanup..." << std::endl;
679 igos.postStage();
680 if (verbosityLevel>=4)
681 std::cout << stagetag << "Cleanup... done" << std::endl;
682
683 if (verbosityLevel>=4)
684 std::cout << stagetag << "Finished." << std::endl;
685 }
686
687 // delete intermediate steps
688 for(unsigned i=1; i<method->s(); ++i) delete x[i];
689
690 // step cleanup
691 if (verbosityLevel>=4)
692 std::cout << mytag << "Cleanup..." << std::endl;
693 igos.postStep();
694 if (verbosityLevel>=4)
695 std::cout << mytag << "Cleanup... done." << std::endl;
696
697 step++;
698 return dt;
699 }
700
701 private:
702
704 class DefaultLimiter
705 {
706 public:
707 template<typename V>
708 void prestage(V& v)
709 {}
710
711 template<typename V>
712 void poststage(V& v)
713 {}
714 };
715
716 const TimeSteppingParameterInterface<T> *method;
717 IGOS& igos;
718 LS& ls;
719 int verbosityLevel;
720 int step;
721 M D;
722 TimeControllerInterface<T> *tc;
723 bool allocated;
724 double ls_reduct;
725 };
726
728 } // end namespace PDELab
729} // end namespace Dune
730#endif // DUNE_PDELAB_INSTATIONARY_EXPLICITONESTEP_HH
limit time step to maximum dt * CFL number
Definition: explicitonestep.hh:65
virtual RealType suggestTimestep(RealType time, RealType givendt) override
Return name of the scheme.
Definition: explicitonestep.hh:82
Base class for all PDELab exceptions.
Definition: exceptions.hh:19
Do one step of an explicit time-stepping scheme.
Definition: explicitonestep.hh:110
T apply(T time, T dt, TrlV &xold, F &f, TrlV &xnew, Limiter &limiter)
do one step;
Definition: explicitonestep.hh:505
ExplicitOneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, LS &ls_, TC &tc_)
construct a new one step scheme
Definition: explicitonestep.hh:200
T apply(T time, T dt, TrlV &xold, F &f, TrlV &xnew)
do one step;
Definition: explicitonestep.hh:489
T apply(T time, T dt, TrlV &xold, TrlV &xnew)
do one step;
Definition: explicitonestep.hh:292
ExplicitOneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, LS &ls_, double ls_reduct_)
Construct a new Explicit One Step Method object.
Definition: explicitonestep.hh:158
T apply(T time, T dt, TrlV &xold, TrlV &xnew, Limiter &limiter)
do one step;
Definition: explicitonestep.hh:307
void setMethod(const TimeSteppingParameterInterface< T > &method_)
redefine the method to be used; can be done before every step
Definition: explicitonestep.hh:278
void setVerbosityLevel(int level)
change verbosity level; 0 means completely quiet
Definition: explicitonestep.hh:246
void setStepNumber(int newstep)
change number of current step
Definition: explicitonestep.hh:255
ExplicitOneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, LS &ls_, TC &tc_, double ls_reduct_)
Construct a new Explicit One Step Method object.
Definition: explicitonestep.hh:220
ExplicitOneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, LS &ls_)
construct a new one step scheme
Definition: explicitonestep.hh:136
void setReduction(const double &ls_reduction_)
Set the Reduction for linear system soltion.
Definition: explicitonestep.hh:268
Insert standard boilerplate into log messages.
Definition: logtag.hh:172
Default time controller; just returns given dt.
Definition: explicitonestep.hh:45
virtual RealType suggestTimestep(RealType time, RealType givendt) override
Return name of the scheme.
Definition: explicitonestep.hh:51
Controller interface for adaptive time stepping.
Definition: explicitonestep.hh:27
virtual ~TimeControllerInterface()
every abstract base class has a virtual destructor
Definition: explicitonestep.hh:36
virtual RealType suggestTimestep(RealType time, RealType givendt)=0
Return name of the scheme.
Utility class for storing and resetting stream attributes.
Definition: ios_state.hh:34
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:796
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
virtual bool implicit() const =0
Return true if method is implicit.
virtual std::string name() const =0
Return name of the scheme.
virtual R d(int r) const =0
Return entries of the d Vector.
virtual unsigned s() const =0
Return number of stages of the method.
Utility class for storing and resetting stream attributes.
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)