dune-pdelab  2.0.0
instationary/onestep.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 
4 #ifndef DUNE_PDELAB_ONESTEP_HH
5 #define DUNE_PDELAB_ONESTEP_HH
6 
7 #include <iomanip>
8 #include <iostream>
9 #include <ostream>
10 #include <vector>
11 
12 #include <stdio.h>
13 
14 #include <dune/common/exceptions.hh>
15 #include <dune/common/fmatrix.hh>
16 #include <dune/common/fvector.hh>
17 #include <dune/common/ios_state.hh>
18 
21 
22 namespace Dune {
23  namespace PDELab {
24 
35  template<class R>
37  {
38  public:
39 
41  {
42  D[0] = 0.0; D[1] = 1.0;
43  A[0][0] = -1.0; A[0][1] = 1.0;
44  B[0][0] = 1.0; B[0][1] = 0.0;
45  }
46 
49  virtual bool implicit () const
50  {
51  return false;
52  }
53 
56  virtual unsigned s () const
57  {
58  return 1;
59  }
60 
64  virtual R a (int r, int i) const
65  {
66  return A[r-1][i];
67  }
68 
72  virtual R b (int r, int i) const
73  {
74  return B[r-1][i];
75  }
76 
80  virtual R d (int i) const
81  {
82  return D[i];
83  }
84 
87  virtual std::string name () const
88  {
89  return std::string("explicit Euler");
90  }
91 
92  private:
93  Dune::FieldVector<R,2> D;
94  Dune::FieldMatrix<R,1,2> A;
95  Dune::FieldMatrix<R,1,2> B;
96  };
97 
106  template<class R>
108  {
111 
112  public:
115  : theta(theta_)
116  {
117  D[0] = 0.0; D[1] = 1.0;
118  A[0][0] = -1.0; A[0][1] = 1.0;
119  B[0][0] = 1.0-theta; B[0][1] = theta;
120  }
121 
124  virtual bool implicit () const
125  {
126  if (theta>0.0)
127  return true;
128  else
129  return false;
130  }
131 
134  virtual unsigned s () const
135  {
136  return 1;
137  }
138 
142  virtual R a (int r, int i) const
143  {
144  return A[r-1][i];
145  }
146 
150  virtual R b (int r, int i) const
151  {
152  return B[r-1][i];
153  }
154 
158  virtual R d (int i) const
159  {
160  return D[i];
161  }
162 
165  virtual std::string name () const
166  {
167  return std::string("one step theta");
168  }
169 
170  private:
171  R theta;
172  Dune::FieldVector<R,2> D;
173  Dune::FieldMatrix<R,1,2> A;
174  Dune::FieldMatrix<R,1,2> B;
175  };
176 
183  template<class R>
185  {
186  public:
187 
189  {
190  D[0] = 0.0; D[1] = 1.0; D[2] = 1.0;
191 
192  A[0][0] = -1.0; A[0][1] = 1.0; A[0][2] = 0.0;
193  A[1][0] = -0.5; A[1][1] = -0.5; A[1][2] = 1.0;
194 
195  B[0][0] = 1.0; B[0][1] = 0.0; B[0][2] = 0.0;
196  B[1][0] = 0.0; B[1][1] = 0.5; B[1][2] = 0.0;
197  }
198 
201  virtual bool implicit () const
202  {
203  return false;
204  }
205 
208  virtual unsigned s () const
209  {
210  return 2;
211  }
212 
216  virtual R a (int r, int i) const
217  {
218  return A[r-1][i];
219  }
220 
224  virtual R b (int r, int i) const
225  {
226  return B[r-1][i];
227  }
228 
232  virtual R d (int i) const
233  {
234  return D[i];
235  }
236 
239  virtual std::string name () const
240  {
241  return std::string("Heun");
242  }
243 
244  private:
245  Dune::FieldVector<R,3> D;
246  Dune::FieldMatrix<R,2,3> A;
247  Dune::FieldMatrix<R,2,3> B;
248  };
249 
256  template<class R>
258  {
259  public:
260 
262  {
263  D[0] = 0.0; D[1] = 1.0; D[2] = 0.5; D[3] = 1.0;
264 
265  A[0][0] = -1.0; A[0][1] = 1.0; A[0][2] = 0.0; A[0][3] = 0.0;
266  A[1][0] = -0.75; A[1][1] = -0.25; A[1][2] = 1.0; A[1][3] = 0.0;
267  A[2][0] = -1.0/3.0; A[2][1] = 0.0; A[2][2] = -2.0/3.0; A[2][3] = 1.0;
268 
269  B[0][0] = 1.0; B[0][1] = 0.0; B[0][2] = 0.0; B[0][3] = 0.0;
270  B[1][0] = 0.0; B[1][1] = 0.25; B[1][2] = 0.0; B[1][3] = 0.0;
271  B[2][0] = 0.0; B[2][1] = 0.0; B[2][2] = 2.0/3.0; B[2][3] = 0.0;
272 
273  }
274 
277  virtual bool implicit () const
278  {
279  return false;
280  }
281 
284  virtual unsigned s () const
285  {
286  return 3;
287  }
288 
292  virtual R a (int r, int i) const
293  {
294  return A[r-1][i];
295  }
296 
300  virtual R b (int r, int i) const
301  {
302  return B[r-1][i];
303  }
304 
308  virtual R d (int i) const
309  {
310  return D[i];
311  }
312 
315  virtual std::string name () const
316  {
317  return std::string("Shu's third order method");
318  }
319 
320  private:
321  Dune::FieldVector<R,4> D;
322  Dune::FieldMatrix<R,3,4> A;
323  Dune::FieldMatrix<R,3,4> B;
324  };
325 
326 
333  template<class R>
335  {
336  public:
337 
339  {
340  D[0] = 0.0; D[1] = 0.5; D[2] = 0.5; D[3] = 1.0; D[4] = 1.0;
341 
342  A[0][0] = -1.0; A[0][1] = 1.0; A[0][2] = 0.0; A[0][3] = 0.0; A[0][4] = 0.0;
343  A[1][0] = -1.0; A[1][1] = 0.0; A[1][2] = 1.0; A[1][3] = 0.0; A[1][4] = 0.0;
344  A[2][0] = -1.0; A[2][1] = 0.0; A[2][2] = 0.0; A[2][3] = 1.0; A[2][4] = 0.0;
345  A[3][0] = -1.0; A[3][1] = 0.0; A[3][2] = 0.0; A[3][3] = 0.0; A[3][4] = 1.0;
346 
347  B[0][0] = 0.5; B[0][1] = 0.0; B[0][2] = 0.0; B[0][3] = 0.0; B[0][4] = 0.0;
348  B[1][0] = 0.0; B[1][1] = 0.5; B[1][2] = 0.0; B[1][3] = 0.0; B[1][4] = 0.0;
349  B[2][0] = 0.0; B[2][1] = 0.0; B[2][2] = 1.0; B[2][3] = 0.0; B[2][4] = 0.0;
350  B[3][0] = 1.0/6.0; B[3][1] = 1.0/3.0; B[3][2] = 1.0/3.0; B[3][3] = 1.0/6.0; B[3][4] = 0.0;
351 
352  }
353 
356  virtual bool implicit () const
357  {
358  return false;
359  }
360 
363  virtual unsigned s () const
364  {
365  return 4;
366  }
367 
371  virtual R a (int r, int i) const
372  {
373  return A[r-1][i];
374  }
375 
379  virtual R b (int r, int i) const
380  {
381  return B[r-1][i];
382  }
383 
387  virtual R d (int i) const
388  {
389  return D[i];
390  }
391 
394  virtual std::string name () const
395  {
396  return std::string("RK4");
397  }
398 
399  private:
400  Dune::FieldVector<R,5> D;
401  Dune::FieldMatrix<R,4,5> A;
402  Dune::FieldMatrix<R,4,5> B;
403  };
404 
405 
406 
407 
414  template<class R>
416  {
417  public:
418 
420  {
421  alpha = 1.0 - 0.5*sqrt(2.0);
422 
423  D[0] = 0.0; D[1] = alpha; D[2] = 1.0;
424 
425  A[0][0] = -1.0; A[0][1] = 1.0; A[0][2] = 0.0;
426  A[1][0] = -1.0; A[1][1] = 0.0; A[1][2] = 1.0;
427 
428  B[0][0] = 0.0; B[0][1] = alpha; B[0][2] = 0.0;
429  B[1][0] = 0.0; B[1][1] = 1.0-alpha; B[1][2] = alpha;
430  }
431 
434  virtual bool implicit () const
435  {
436  return true;
437  }
438 
441  virtual unsigned s () const
442  {
443  return 2;
444  }
445 
449  virtual R a (int r, int i) const
450  {
451  return A[r-1][i];
452  }
453 
457  virtual R b (int r, int i) const
458  {
459  return B[r-1][i];
460  }
461 
465  virtual R d (int i) const
466  {
467  return D[i];
468  }
469 
472  virtual std::string name () const
473  {
474  return std::string("Alexander (order 2)");
475  }
476 
477  private:
478  R alpha;
479  Dune::FieldVector<R,3> D;
480  Dune::FieldMatrix<R,2,3> A;
481  Dune::FieldMatrix<R,2,3> B;
482  };
483 
490  template<class R>
492  {
493  public:
494 
496  {
497  R alpha, theta, thetap, beta;
498  theta = 1.0 - 0.5*sqrt(2.0);
499  thetap = 1.0-2.0*theta;
500  alpha = 2.0-sqrt(2.0);
501  beta = 1.0-alpha;
502 
503  D[0] = 0.0; D[1] = theta; D[2] = 1.0-theta; D[3] = 1.0;
504 
505  A[0][0] = -1.0; A[0][1] = 1.0; A[0][2] = 0.0; A[0][3] = 0.0;
506  A[1][0] = 0.0; A[1][1] = -1.0; A[1][2] = 1.0; A[1][3] = 0.0;
507  A[2][0] = 0.0; A[2][1] = 0.0; A[2][2] = -1.0; A[2][3] = 1.0;
508 
509  B[0][0] = beta*theta; B[0][1] = alpha*theta; B[0][2] = 0.0; B[0][3] = 0.0;
510  B[1][0] = 0.0; B[1][1] = alpha*thetap; B[1][2] = alpha*theta; B[1][3] = 0.0;
511  B[2][0] = 0.0; B[2][1] = 0.0; B[2][2] = beta*theta; B[2][3] = alpha*theta;
512  }
513 
516  virtual bool implicit () const
517  {
518  return true;
519  }
520 
523  virtual unsigned s () const
524  {
525  return 3;
526  }
527 
531  virtual R a (int r, int i) const
532  {
533  return A[r-1][i];
534  }
535 
539  virtual R b (int r, int i) const
540  {
541  return B[r-1][i];
542  }
543 
547  virtual R d (int i) const
548  {
549  return D[i];
550  }
551 
554  virtual std::string name () const
555  {
556  return std::string("Fractional step theta");
557  }
558 
559  private:
560  Dune::FieldVector<R,4> D;
561  Dune::FieldMatrix<R,3,4> A;
562  Dune::FieldMatrix<R,3,4> B;
563  };
564 
572  template<class R>
574  {
575  public:
576 
578  {
579  R alpha = 0.4358665215;
580 
581  // Newton iteration for alpha
582  for (int i=1; i<=10; i++)
583  {
584  alpha = alpha - (alpha*(alpha*alpha-3.0*(alpha-0.5))-1.0/6.0)/(3.0*alpha*(alpha-2.0)+1.5);
585  // std::cout.precision(16);
586  // std::cout << "alpha "
587  // << std::setw(8) << i << " "
588  // << std::scientific << alpha << std::endl;
589  }
590 
591  R tau2 = (1.0+alpha)*0.5;
592  R b1 = -(6.0*alpha*alpha -16.0*alpha + 1.0)*0.25;
593  R b2 = (6*alpha*alpha - 20.0*alpha + 5.0)*0.25;
594 
595  // std::cout.precision(16);
596  // std::cout << "alpha " << std::scientific << alpha << std::endl;
597  // std::cout << "tau2 " << std::scientific << tau2 << std::endl;
598  // std::cout << "b1 " << std::scientific << b1 << std::endl;
599  // std::cout << "b2 " << std::scientific << b2 << std::endl;
600 
601  D[0] = 0.0; D[1] = alpha; D[2] = tau2; D[3] = 1.0;
602 
603  A[0][0] = -1.0; A[0][1] = 1.0; A[0][2] = 0.0; A[0][3] = 0.0;
604  A[1][0] = -1.0; A[1][1] = 0.0; A[1][2] = 1.0; A[1][3] = 0.0;
605  A[2][0] = -1.0; A[2][1] = 0.0; A[2][2] = 0.0; A[2][3] = 1.0;
606 
607  B[0][0] = 0.0; B[0][1] = alpha; B[0][2] = 0.0; B[0][3] = 0.0;
608  B[1][0] = 0.0; B[1][1] = tau2-alpha; B[1][2] = alpha; B[1][3] = 0.0;
609  B[2][0] = 0.0; B[2][1] = b1; B[2][2] = b2; B[2][3] = alpha;
610  }
611 
614  virtual bool implicit () const
615  {
616  return true;
617  }
618 
621  virtual unsigned s () const
622  {
623  return 3;
624  }
625 
629  virtual R a (int r, int i) const
630  {
631  return A[r-1][i];
632  }
633 
637  virtual R b (int r, int i) const
638  {
639  return B[r-1][i];
640  }
641 
645  virtual R d (int i) const
646  {
647  return D[i];
648  }
649 
652  virtual std::string name () const
653  {
654  return std::string("Alexander (claims order 3)");
655  }
656 
657  private:
658  R alpha, theta, thetap, beta;
659  Dune::FieldVector<R,4> D;
660  Dune::FieldMatrix<R,3,4> A;
661  Dune::FieldMatrix<R,3,4> B;
662  };
663 
664 
669  template<class R>
671  {
672  public:
673  typedef R RealType;
674 
677  virtual RealType suggestTimestep (RealType time, RealType givendt) = 0;
678 
681  };
682 
684 
687  template<class R>
689  {
690  public:
691  typedef R RealType;
692 
695  virtual RealType suggestTimestep (RealType time, RealType givendt)
696  {
697  return givendt;
698  }
699  };
700 
701 
703 
707  template<class R, class IGOS>
709  {
710  public:
711  typedef R RealType;
712 
713  CFLTimeController (R cfl_, const IGOS& igos_) : cfl(cfl_), target(1e100), igos(igos_)
714  {}
715 
716  CFLTimeController (R cfl_, R target_, const IGOS& igos_) : cfl(cfl_), target(target_), igos(igos_)
717  {}
718 
719  void setTarget (R target_)
720  {
721  target = target_;
722  }
723 
726  virtual RealType suggestTimestep (RealType time, RealType givendt)
727  {
728  RealType suggested = cfl*igos.suggestTimestep(givendt);
729  if (time+2.0*suggested<target)
730  return suggested;
731  if (time+suggested<target)
732  return 0.5*(target-time);
733  return target-time;
734  }
735 
736  private:
737  R cfl;
738  R target;
739  const IGOS& igos;
740  };
741 
742 
743  // Status information of Newton's method
745  {
746  unsigned int timesteps;
751 
753  timesteps(0),
754  assembler_time(0.0),
755  linear_solver_time(0.0),
758  {}
759  };
760 
762  {
766  {}
767  };
768 
770 
777  template<class T, class IGOS, class PDESOLVER, class TrlV, class TstV = TrlV>
779  {
780  typedef typename PDESOLVER::Result PDESolverResult;
781 
782  public:
784 
786 
798  IGOS& igos_, PDESOLVER& pdesolver_)
799  : method(&method_), igos(igos_), pdesolver(pdesolver_), verbosityLevel(1), step(1), res()
800  {
801  if (igos.trialGridFunctionSpace().gridView().comm().rank()>0)
802  verbosityLevel = 0;
803  }
804 
806  void setVerbosityLevel (int level)
807  {
808  if (igos.trialGridFunctionSpace().gridView().comm().rank()>0)
809  verbosityLevel = 0;
810  else
811  verbosityLevel = level;
812  }
813 
815  void setStepNumber(int newstep) { step = newstep; }
816 
818  const PDESOLVER & getPDESolver() const
819  {
820  return pdesolver;
821  }
822 
824  PDESOLVER & getPDESolver()
825  {
826  return pdesolver;
827  }
828 
829  const Result& result() const
830  {
831  return res;
832  }
833 
835 
843  {
844  method = &method_;
845  }
846 
854  T apply (T time, T dt, TrlV& xold, TrlV& xnew)
855  {
856  // save formatting attributes
857  ios_base_all_saver format_attribute_saver(std::cout);
858 
859  // do statistics
860  OneStepMethodPartialResult step_result;
861 
862  std::vector<TrlV*> x(1); // vector of pointers to all steps
863  x[0] = &xold; // initially we have only one
864 
865  if (verbosityLevel>=1){
866  std::ios_base::fmtflags oldflags = std::cout.flags();
867  std::cout << "TIME STEP [" << method->name() << "] "
868  << std::setw(6) << step
869  << " time (from): "
870  << std::setw(12) << std::setprecision(4) << std::scientific
871  << time
872  << " dt: "
873  << std::setw(12) << std::setprecision(4) << std::scientific
874  << dt
875  << " time (to): "
876  << std::setw(12) << std::setprecision(4) << std::scientific
877  << time+dt
878  << std::endl;
879  std::cout.flags(oldflags);
880  }
881 
882  // prepare assembler
883  igos.preStep(*method,time,dt);
884 
885  // loop over all stages
886  for (unsigned r=1; r<=method->s(); ++r)
887  {
888  if (verbosityLevel>=2){
889  std::ios_base::fmtflags oldflags = std::cout.flags();
890  std::cout << "STAGE "
891  << r
892  << " time (to): "
893  << std::setw(12) << std::setprecision(4) << std::scientific
894  << time+method->d(r)*dt
895  << "." << std::endl;
896  std::cout.flags(oldflags);
897  }
898 
899  // prepare stage
900  igos.preStage(r,x);
901 
902  // get vector for current stage
903  if (r==method->s())
904  {
905  // last stage
906  x.push_back(&xnew);
907  if (r>1) xnew = *(x[r-1]); // if r=1 then xnew has already initial guess
908  }
909  else
910  {
911  // intermediate step
912  x.push_back(new TrlV(igos.trialGridFunctionSpace()));
913  if (r>1)
914  *(x[r]) = *(x[r-1]); // use result of last stage as initial guess
915  else
916  *(x[r]) = xnew;
917  }
918 
919  // solve stage
920  try {
921  pdesolver.apply(*x[r]);
922  }
923  catch (...)
924  {
925  // time step failed -> accumulate to total only
926  PDESolverResult pderes = pdesolver.result();
927  step_result.assembler_time += pderes.assembler_time;
928  step_result.linear_solver_time += pderes.linear_solver_time;
929  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
930  step_result.nonlinear_solver_iterations += pderes.iterations;
931  res.total.assembler_time += step_result.assembler_time;
932  res.total.linear_solver_time += step_result.linear_solver_time;
935  res.total.timesteps += 1;
936  throw;
937  }
938  PDESolverResult pderes = pdesolver.result();
939  step_result.assembler_time += pderes.assembler_time;
940  step_result.linear_solver_time += pderes.linear_solver_time;
941  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
942  step_result.nonlinear_solver_iterations += pderes.iterations;
943 
944  // stage cleanup
945  igos.postStage();
946  }
947 
948  // delete intermediate steps
949  for (unsigned i=1; i<method->s(); ++i) delete x[i];
950 
951  // step cleanup
952  igos.postStep();
953 
954  // update statistics
955  res.total.assembler_time += step_result.assembler_time;
956  res.total.linear_solver_time += step_result.linear_solver_time;
959  res.total.timesteps += 1;
960  res.successful.assembler_time += step_result.assembler_time;
964  res.successful.timesteps += 1;
965  if (verbosityLevel>=1){
966  std::ios_base::fmtflags oldflags = std::cout.flags();
967  std::cout << "::: timesteps " << std::setw(6) << res.successful.timesteps
968  << " (" << res.total.timesteps << ")" << std::endl;
969  std::cout << "::: nl iterations " << std::setw(6) << res.successful.nonlinear_solver_iterations
970  << " (" << res.total.nonlinear_solver_iterations << ")" << std::endl;
971  std::cout << "::: lin iterations " << std::setw(6) << res.successful.linear_solver_iterations
972  << " (" << res.total.linear_solver_iterations << ")" << std::endl;
973  std::cout << "::: assemble time " << std::setw(12) << std::setprecision(4) << std::scientific
974  << res.successful.assembler_time << " (" << res.total.assembler_time << ")" << std::endl;
975  std::cout << "::: lin solve time " << std::setw(12) << std::setprecision(4) << std::scientific
976  << res.successful.linear_solver_time << " (" << res.total.linear_solver_time << ")" << std::endl;
977  std::cout.flags(oldflags);
978  }
979 
980  step++;
981  return dt;
982  }
983 
994  template<typename F>
995  T apply (T time, T dt, TrlV& xold, F& f, TrlV& xnew)
996  {
997  // do statistics
998  OneStepMethodPartialResult step_result;
999 
1000  // save formatting attributes
1001  ios_base_all_saver format_attribute_saver(std::cout);
1002 
1003  std::vector<TrlV*> x(1); // vector of pointers to all steps
1004  x[0] = &xold; // initially we have only one
1005 
1006  if (verbosityLevel>=1){
1007  std::ios_base::fmtflags oldflags = std::cout.flags();
1008  std::cout << "TIME STEP [" << method->name() << "] "
1009  << std::setw(6) << step
1010  << " time (from): "
1011  << std::setw(12) << std::setprecision(4) << std::scientific
1012  << time
1013  << " dt: "
1014  << std::setw(12) << std::setprecision(4) << std::scientific
1015  << dt
1016  << " time (to): "
1017  << std::setw(12) << std::setprecision(4) << std::scientific
1018  << time+dt
1019  << std::endl;
1020  std::cout.flags(oldflags);
1021  }
1022 
1023  // prepare assembler
1024  igos.preStep(*method,time,dt);
1025 
1026  // loop over all stages
1027  for (unsigned r=1; r<=method->s(); ++r)
1028  {
1029  if (verbosityLevel>=2){
1030  std::ios_base::fmtflags oldflags = std::cout.flags();
1031  std::cout << "STAGE "
1032  << r
1033  << " time (to): "
1034  << std::setw(12) << std::setprecision(4) << std::scientific
1035  << time+method->d(r)*dt
1036  << "." << std::endl;
1037  std::cout.flags(oldflags);
1038  }
1039 
1040  // prepare stage
1041  igos.preStage(r,x);
1042 
1043  // get vector for current stage
1044  if (r==method->s())
1045  {
1046  // last stage
1047  x.push_back(&xnew);
1048  }
1049  else
1050  {
1051  // intermediate step
1052  x.push_back(new TrlV(igos.trialGridFunctionSpace()));
1053  }
1054 
1055  // set boundary conditions and initial value
1056  igos.interpolate(r,*x[r-1],f,*x[r]);
1057 
1058  // solve stage
1059  try {
1060  pdesolver.apply(*x[r]);
1061  }
1062  catch (...)
1063  {
1064  // time step failed -> accumulate to total only
1065  PDESolverResult pderes = pdesolver.result();
1066  step_result.assembler_time += pderes.assembler_time;
1067  step_result.linear_solver_time += pderes.linear_solver_time;
1068  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
1069  step_result.nonlinear_solver_iterations += pderes.iterations;
1070  res.total.assembler_time += step_result.assembler_time;
1071  res.total.linear_solver_time += step_result.linear_solver_time;
1074  res.total.timesteps += 1;
1075  throw;
1076  }
1077  PDESolverResult pderes = pdesolver.result();
1078  step_result.assembler_time += pderes.assembler_time;
1079  step_result.linear_solver_time += pderes.linear_solver_time;
1080  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
1081  step_result.nonlinear_solver_iterations += pderes.iterations;
1082 
1083  // stage cleanup
1084  igos.postStage();
1085  }
1086 
1087  // delete intermediate steps
1088  for (unsigned i=1; i<method->s(); ++i) delete x[i];
1089 
1090  // step cleanup
1091  igos.postStep();
1092 
1093  // update statistics
1094  res.total.assembler_time += step_result.assembler_time;
1095  res.total.linear_solver_time += step_result.linear_solver_time;
1098  res.total.timesteps += 1;
1099  res.successful.assembler_time += step_result.assembler_time;
1103  res.successful.timesteps += 1;
1104  if (verbosityLevel>=1){
1105  std::ios_base::fmtflags oldflags = std::cout.flags();
1106  std::cout << "::: timesteps " << std::setw(6) << res.successful.timesteps
1107  << " (" << res.total.timesteps << ")" << std::endl;
1108  std::cout << "::: nl iterations " << std::setw(6) << res.successful.nonlinear_solver_iterations
1109  << " (" << res.total.nonlinear_solver_iterations << ")" << std::endl;
1110  std::cout << "::: lin iterations " << std::setw(6) << res.successful.linear_solver_iterations
1111  << " (" << res.total.linear_solver_iterations << ")" << std::endl;
1112  std::cout << "::: assemble time " << std::setw(12) << std::setprecision(4) << std::scientific
1113  << res.successful.assembler_time << " (" << res.total.assembler_time << ")" << std::endl;
1114  std::cout << "::: lin solve time " << std::setw(12) << std::setprecision(4) << std::scientific
1115  << res.successful.linear_solver_time << " (" << res.total.linear_solver_time << ")" << std::endl;
1116  std::cout.flags(oldflags);
1117  }
1118 
1119  step++;
1120  return dt;
1121  }
1122 
1123  private:
1124  const TimeSteppingParameterInterface<T> *method;
1125  IGOS& igos;
1126  PDESOLVER& pdesolver;
1127  int verbosityLevel;
1128  int step;
1129  Result res;
1130  };
1131 
1133 
1141  template<class T, class IGOS, class LS, class TrlV, class TstV = TrlV, class TC = SimpleTimeController<T> >
1143  {
1144  typedef typename TrlV::ElementType Real;
1145  typedef typename IGOS::template MatrixContainer<Real>::Type M;
1146 
1147  public:
1149 
1160  ExplicitOneStepMethod(const TimeSteppingParameterInterface<T>& method_, IGOS& igos_, LS& ls_)
1161  : method(&method_), igos(igos_), ls(ls_), verbosityLevel(1), step(1), D(igos),
1162  tc(new SimpleTimeController<T>()), allocated(true)
1163  {
1164  if (method->implicit())
1165  DUNE_THROW(Exception,"explicit one step method called with implicit scheme");
1166  if (igos.trialGridFunctionSpace().gridView().comm().rank()>0)
1167  verbosityLevel = 0;
1168  }
1169 
1171 
1182  ExplicitOneStepMethod(const TimeSteppingParameterInterface<T>& method_, IGOS& igos_, LS& ls_, TC& tc_)
1183  : method(&method_), igos(igos_), ls(ls_), verbosityLevel(1), step(1), D(igos),
1184  tc(&tc_), allocated(false)
1185  {
1186  if (method->implicit())
1187  DUNE_THROW(Exception,"explicit one step method called with implicit scheme");
1188  }
1189 
1191  {
1192  if (allocated) delete tc;
1193  }
1194 
1196  void setVerbosityLevel (int level)
1197  {
1198  if (igos.trialGridFunctionSpace().gridView().comm().rank()>0)
1199  verbosityLevel = 0;
1200  else
1201  verbosityLevel = level;
1202  }
1203 
1205  void setStepNumber(int newstep) { step = newstep; }
1206 
1208 
1216  {
1217  method = &method_;
1218  if (method->implicit())
1219  DUNE_THROW(Exception,"explicit one step method called with implicit scheme");
1220  }
1221 
1230  T apply (T time, T dt, TrlV& xold, TrlV& xnew)
1231  {
1232  DefaultLimiter limiter;
1233  return apply(time,dt,xold,xnew,limiter);
1234  }
1235 
1236  template<typename Limiter>
1237  T apply (T time, T dt, TrlV& xold, TrlV& xnew, Limiter& limiter)
1238  {
1239  // save formatting attributes
1240  ios_base_all_saver format_attribute_saver(std::cout);
1241  LocalTag mytag;
1242  mytag << "ExplicitOneStepMethod::apply(): ";
1243 
1244  std::vector<TrlV*> x(1); // vector of pointers to all steps
1245  x[0] = &xold; // initially we have only one
1246  if(verbosityLevel>=4)
1247  std::cout << mytag << "Creating residual vectors alpha and beta..."
1248  << std::endl;
1249  TstV alpha(igos.testGridFunctionSpace()), beta(igos.testGridFunctionSpace()); // split residual vectors
1250  if(verbosityLevel>=4)
1251  std::cout << mytag
1252  << "Creating residual vectors alpha and beta... done."
1253  << std::endl;
1254 
1255  if (verbosityLevel>=1){
1256  std::ios_base::fmtflags oldflags = std::cout.flags();
1257  std::cout << "TIME STEP [" << method->name() << "] "
1258  << std::setw(6) << step
1259  << " time (from): "
1260  << std::setw(12) << std::setprecision(4) << std::scientific
1261  << time
1262  << " dt: "
1263  << std::setw(12) << std::setprecision(4) << std::scientific
1264  << dt
1265  << " time (to): "
1266  << std::setw(12) << std::setprecision(4) << std::scientific
1267  << time+dt
1268  << std::endl;
1269  std::cout.flags(oldflags);
1270  }
1271 
1272  // prepare assembler
1273  if(verbosityLevel>=4)
1274  std::cout << mytag << "Preparing assembler..." << std::endl;
1275  igos.preStep(*method,time,dt);
1276  if(verbosityLevel>=4)
1277  std::cout << mytag << "Preparing assembler... done." << std::endl;
1278 
1279  // loop over all stages
1280  for(unsigned r=1; r<=method->s(); ++r)
1281  {
1282  LocalTag stagetag(mytag);
1283  stagetag << "stage " << r << ": ";
1284  if (verbosityLevel>=4)
1285  std::cout << stagetag << "Start." << std::endl;
1286 
1287  if (verbosityLevel>=2){
1288  std::ios_base::fmtflags oldflags = std::cout.flags();
1289  std::cout << "STAGE "
1290  << r
1291  << " time (to): "
1292  << std::setw(12) << std::setprecision(4) << std::scientific
1293  << time+method->d(r)*dt
1294  << "." << std::endl;
1295  std::cout.flags(oldflags);
1296  }
1297 
1298  // get vector for current stage
1299  if (r==method->s())
1300  {
1301  // last stage
1302  x.push_back(&xnew);
1303  if (r>1) xnew = *(x[r-1]); // if r=1 then xnew has already initial guess
1304  }
1305  else
1306  {
1307  // intermediate step
1308  x.push_back(new TrlV(igos.trialGridFunctionSpace()));
1309  if (r>1)
1310  *(x[r]) = *(x[r-1]); // use result of last stage as initial guess
1311  else
1312  *(x[r]) = xnew;
1313  }
1314 
1315  // compute residuals and jacobian
1316  if (verbosityLevel>=4) std::cout << "assembling D, alpha, beta ..." << std::endl;
1317  D = Real(0.0);
1318  alpha = 0.0;
1319  beta = 0.0;
1320 
1321  //apply slope limiter to old solution (e.g for finite volume reconstruction scheme)
1322  limiter.prestage(*x[r-1]);
1323 
1324  if(verbosityLevel>=4)
1325  std::cout << stagetag << "Assembling residual..." << std::endl;
1326  igos.explicit_jacobian_residual(r,x,D,alpha,beta);
1327  if(verbosityLevel>=4)
1328  std::cout << stagetag << "Assembling residual... done."
1329  << std::endl;
1330 
1331  // let time controller compute the optimal dt in first stage
1332  if (r==1)
1333  {
1334  T newdt = tc->suggestTimestep(time,dt);
1335  newdt = std::min(newdt, dt);
1336 
1337  if (verbosityLevel>=4){
1338  std::ios_base::fmtflags oldflags = std::cout.flags();
1339  std::cout << stagetag
1340  << "current dt: "
1341  << std::setw(12) << std::setprecision(4) << std::scientific
1342  << dt
1343  << " suggested dt: "
1344  << std::setw(12) << std::setprecision(4) << std::scientific
1345  << newdt
1346  << std::endl;
1347  std::cout.flags(oldflags);
1348  }
1349 
1350  if (verbosityLevel>=2 && newdt!=dt)
1351  {
1352  std::ios_base::fmtflags oldflags = std::cout.flags();
1353  std::cout << "changed dt to "
1354  << std::setw(12) << std::setprecision(4) << std::scientific
1355  << newdt
1356  << std::endl;
1357  std::cout.flags(oldflags);
1358  }
1359  dt = newdt;
1360  }
1361 
1362  // combine residual with selected dt
1363  if (verbosityLevel>=4)
1364  std::cout << stagetag
1365  << "Combining residuals with selected dt..."
1366  << std::endl;
1367  alpha.axpy(dt,beta);
1368  if (verbosityLevel>=4)
1369  std::cout << stagetag
1370  << "Combining residuals with selected dt... done."
1371  << std::endl;
1372 
1373  // solve diagonal system
1374  if (verbosityLevel>=4)
1375  std::cout << stagetag << "Solving diagonal system..."
1376  << std::endl;
1377  ls.apply(D,*x[r],alpha,0.99); // dummy reduction
1378  if (verbosityLevel>=4)
1379  std::cout << stagetag << "Solving diagonal system... done."
1380  << std::endl;
1381 
1382  // apply slope limiter to new solution (e.g DG scheme)
1383  limiter.poststage(*x[r]);
1384 
1385  // stage cleanup
1386  if (verbosityLevel>=4)
1387  std::cout << stagetag << "Cleanup..." << std::endl;
1388  igos.postStage();
1389  if (verbosityLevel>=4)
1390  std::cout << stagetag << "Cleanup... done" << std::endl;
1391 
1392  if (verbosityLevel>=4)
1393  std::cout << stagetag << "Finished." << std::endl;
1394  }
1395 
1396  // delete intermediate steps
1397  for(unsigned i=1; i<method->s(); ++i) delete x[i];
1398 
1399  // step cleanup
1400  if (verbosityLevel>=4)
1401  std::cout << mytag << "Cleanup..." << std::endl;
1402  igos.postStep();
1403  if (verbosityLevel>=4)
1404  std::cout << mytag << "Cleanup... done." << std::endl;
1405 
1406  step++;
1407  return dt;
1408  }
1409 
1410  private:
1411 
1413  class DefaultLimiter
1414  {
1415  public:
1416  template<typename V>
1417  void prestage(V& v)
1418  {}
1419 
1420  template<typename V>
1421  void poststage(V& v)
1422  {}
1423  };
1424 
1425  const TimeSteppingParameterInterface<T> *method;
1426  IGOS& igos;
1427  LS& ls;
1428  int verbosityLevel;
1429  int step;
1430  M D;
1431  TimeControllerInterface<T> *tc;
1432  bool allocated;
1433  };
1434 
1436  {
1437  public:
1438  FilenameHelper(const char *basename_, int i_=0)
1439  : i(i_)
1440  {
1441  sprintf(basename,"%s",basename_);
1442  }
1443 
1444  FilenameHelper(const std::string & basename_, int i_=0)
1445  : i(i_)
1446  {
1447  sprintf(basename,"%s",basename_.c_str());
1448  }
1449 
1450  const char *getName (int i_)
1451  {
1452  sprintf(fname,"%s-%05d",basename,i_);
1453  return fname;
1454  }
1455 
1456  const char *getName ()
1457  {
1458  sprintf(fname,"%s-%05d",basename,i);
1459  return fname;
1460  }
1461 
1462  void increment ()
1463  {
1464  i++;
1465  }
1466 
1467  private:
1468  char fname[255];
1469  char basename[255];
1470  int i;
1471  };
1472 
1474  } // namespace PDELab
1475 } // namespace Dune
1476 
1477 #endif
CFLTimeController(R cfl_, R target_, const IGOS &igos_)
Definition: instationary/onestep.hh:716
virtual R a(int r, int i) const
Return entries of the A matrix.
Definition: instationary/onestep.hh:371
Parameters to turn the OneStepMethod into an Alexander3 scheme.
Definition: instationary/onestep.hh:573
virtual std::string name() const
Return name of the scheme.
Definition: instationary/onestep.hh:239
virtual R a(int r, int i) const
Return entries of the A matrix.
Definition: instationary/onestep.hh:64
virtual std::string name() const
Return name of the scheme.
Definition: instationary/onestep.hh:394
virtual std::string name() const =0
Return name of the scheme.
void setMethod(const TimeSteppingParameterInterface< T > &method_)
redefine the method to be used; can be done before every step
Definition: instationary/onestep.hh:842
Alexander2Parameter()
Definition: instationary/onestep.hh:419
OneStepThetaParameter(R theta_)
construct OneStepThetaParameter class
Definition: instationary/onestep.hh:114
virtual R b(int r, int i) const
Return entries of the B matrix.
Definition: instationary/onestep.hh:150
virtual R b(int r, int i) const
Return entries of the B matrix.
Definition: instationary/onestep.hh:72
virtual unsigned s() const
Return number of stages s of the method.
Definition: instationary/onestep.hh:208
OneStepMethodResult()
Definition: instationary/onestep.hh:765
OneStepMethodResult Result
Definition: instationary/onestep.hh:783
virtual std::string name() const
Return name of the scheme.
Definition: instationary/onestep.hh:472
void increment()
Definition: instationary/onestep.hh:1462
virtual R b(int r, int i) const
Return entries of the B matrix.
Definition: instationary/onestep.hh:224
virtual unsigned s() const
Return number of stages s of the method.
Definition: instationary/onestep.hh:284
Base parameter class for time stepping scheme parameters.
Definition: timesteppingparameterinterface.hh:19
OneStepMethodPartialResult()
Definition: instationary/onestep.hh:752
Controller interface for adaptive time stepping.
Definition: instationary/onestep.hh:670
virtual R a(int r, int i) const
Return entries of the A matrix.
Definition: instationary/onestep.hh:531
virtual R a(int r, int i) const
Return entries of the A matrix.
Definition: instationary/onestep.hh:142
virtual std::string name() const
Return name of the scheme.
Definition: instationary/onestep.hh:165
virtual std::string name() const
Return name of the scheme.
Definition: instationary/onestep.hh:554
virtual bool implicit() const
Return true if method is implicit.
Definition: instationary/onestep.hh:49
Do one step of an explicit time-stepping scheme.
Definition: instationary/onestep.hh:1142
virtual RealType suggestTimestep(RealType time, RealType givendt)
Return name of the scheme.
Definition: instationary/onestep.hh:726
unsigned int timesteps
Definition: instationary/onestep.hh:746
virtual bool implicit() const
Return true if method is implicit.
Definition: instationary/onestep.hh:201
Default time controller; just returns given dt.
Definition: instationary/onestep.hh:688
~ExplicitOneStepMethod()
Definition: instationary/onestep.hh:1190
virtual unsigned s() const
Return number of stages s of the method.
Definition: instationary/onestep.hh:621
void setStepNumber(int newstep)
change number of current step
Definition: instationary/onestep.hh:1205
Parameters to turn the OneStepMethod into an Alexander scheme.
Definition: instationary/onestep.hh:415
ExplicitEulerParameter()
Definition: instationary/onestep.hh:40
void setVerbosityLevel(int level)
change verbosity level; 0 means completely quiet
Definition: instationary/onestep.hh:1196
virtual bool implicit() const
Return true if method is implicit.
Definition: instationary/onestep.hh:277
Shu3Parameter()
Definition: instationary/onestep.hh:261
R RealType
Definition: instationary/onestep.hh:711
HeunParameter()
Definition: instationary/onestep.hh:188
virtual unsigned s() const =0
Return number of stages of the method.
Insert standard boilerplate into log messages.
Definition: logtag.hh:181
virtual R b(int r, int i) const
Return entries of the B matrix.
Definition: instationary/onestep.hh:379
double linear_solver_time
Definition: instationary/onestep.hh:748
virtual bool implicit() const
Return true if method is implicit.
Definition: instationary/onestep.hh:516
Alexander3Parameter()
Definition: instationary/onestep.hh:577
T apply(T time, T dt, TrlV &xold, TrlV &xnew, Limiter &limiter)
Definition: instationary/onestep.hh:1237
void setStepNumber(int newstep)
change number of current step
Definition: instationary/onestep.hh:815
OneStepMethodPartialResult total
Definition: instationary/onestep.hh:763
virtual unsigned s() const
Return number of stages s of the method.
Definition: instationary/onestep.hh:523
T apply(T time, T dt, TrlV &xold, TrlV &xnew)
do one step;
Definition: instationary/onestep.hh:854
const PDESOLVER & getPDESolver() const
Access to the (non) linear solver.
Definition: instationary/onestep.hh:818
virtual unsigned s() const
Return number of stages s of the method.
Definition: instationary/onestep.hh:363
virtual R d(int i) const
Return entries of the d Vector.
Definition: instationary/onestep.hh:158
T apply(T time, T dt, TrlV &xold, TrlV &xnew)
do one step;
Definition: instationary/onestep.hh:1230
T apply(T time, T dt, TrlV &xold, F &f, TrlV &xnew)
do one step; This is a version which interpolates constraints at the start of each stage ...
Definition: instationary/onestep.hh:995
virtual std::string name() const
Return name of the scheme.
Definition: instationary/onestep.hh:315
virtual R b(int r, int i) const
Return entries of the B matrix.
Definition: instationary/onestep.hh:457
Parameters to turn the OneStepMethod into an one step theta method.
Definition: instationary/onestep.hh:107
FractionalStepParameter()
Definition: instationary/onestep.hh:495
void setMethod(const TimeSteppingParameterInterface< T > &method_)
redefine the method to be used; can be done before every step
Definition: instationary/onestep.hh:1215
Parameters to turn the ExplicitOneStepMethod into a Heun scheme.
Definition: instationary/onestep.hh:184
OneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, PDESOLVER &pdesolver_)
construct a new one step scheme
Definition: instationary/onestep.hh:797
virtual R d(int i) const
Return entries of the d Vector.
Definition: instationary/onestep.hh:387
R RealType
Definition: instationary/onestep.hh:673
ExplicitOneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, LS &ls_)
construct a new one step scheme
Definition: instationary/onestep.hh:1160
virtual R b(int r, int i) const
Return entries of the B matrix.
Definition: instationary/onestep.hh:637
virtual R d(int i) const
Return entries of the d Vector.
Definition: instationary/onestep.hh:645
Parameters to turn the ExplicitOneStepMethod into an explicite Euler method.
Definition: instationary/onestep.hh:36
PDESOLVER & getPDESolver()
Access to the (non) linear solver.
Definition: instationary/onestep.hh:824
ExplicitOneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, LS &ls_, TC &tc_)
construct a new one step scheme
Definition: instationary/onestep.hh:1182
virtual R d(int i) const
Return entries of the d Vector.
Definition: instationary/onestep.hh:80
virtual R b(int r, int i) const
Return entries of the B matrix.
Definition: instationary/onestep.hh:300
int linear_solver_iterations
Definition: instationary/onestep.hh:749
virtual ~TimeControllerInterface()
every abstract base class has a virtual destructor
Definition: instationary/onestep.hh:680
virtual std::string name() const
Return name of the scheme.
Definition: instationary/onestep.hh:652
virtual bool implicit() const =0
Return true if method is implicit.
OneStepMethodPartialResult successful
Definition: instationary/onestep.hh:764
virtual bool implicit() const
Return true if method is implicit.
Definition: instationary/onestep.hh:124
virtual RealType suggestTimestep(RealType time, RealType givendt)=0
Return name of the scheme.
virtual R d(int i) const
Return entries of the d Vector.
Definition: instationary/onestep.hh:465
const char * getName()
Definition: instationary/onestep.hh:1456
virtual R d(int r) const =0
Return entries of the d Vector.
virtual unsigned s() const
Return number of stages s of the method.
Definition: instationary/onestep.hh:441
Parameters to turn the ExplicitOneStepMethod into a third order strong stability preserving (SSP) sch...
Definition: instationary/onestep.hh:257
limit time step to maximum dt * CFL number
Definition: instationary/onestep.hh:708
R RealType
Definition: instationary/onestep.hh:691
Parameters to turn the OneStepMethod into a fractional step theta scheme.
Definition: instationary/onestep.hh:491
virtual R a(int r, int i) const
Return entries of the A matrix.
Definition: instationary/onestep.hh:216
int nonlinear_solver_iterations
Definition: instationary/onestep.hh:750
const F & f
Definition: common/constraints.hh:145
virtual R a(int r, int i) const
Return entries of the A matrix.
Definition: instationary/onestep.hh:292
virtual R d(int i) const
Return entries of the d Vector.
Definition: instationary/onestep.hh:547
CFLTimeController(R cfl_, const IGOS &igos_)
Definition: instationary/onestep.hh:713
Definition: instationary/onestep.hh:761
virtual R d(int i) const
Return entries of the d Vector.
Definition: instationary/onestep.hh:308
Do one step of a time-stepping scheme.
Definition: instationary/onestep.hh:778
virtual bool implicit() const
Return true if method is implicit.
Definition: instationary/onestep.hh:614
Parameters to turn the ExplicitOneStepMethod into a classical fourth order Runge-Kutta method...
Definition: instationary/onestep.hh:334
void setVerbosityLevel(int level)
change verbosity level; 0 means completely quiet
Definition: instationary/onestep.hh:806
void setTarget(R target_)
Definition: instationary/onestep.hh:719
virtual unsigned s() const
Return number of stages s of the method.
Definition: instationary/onestep.hh:56
virtual R d(int i) const
Return entries of the d Vector.
Definition: instationary/onestep.hh:232
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
virtual bool implicit() const
Return true if method is implicit.
Definition: instationary/onestep.hh:356
FilenameHelper(const char *basename_, int i_=0)
Definition: instationary/onestep.hh:1438
virtual R a(int r, int i) const
Return entries of the A matrix.
Definition: instationary/onestep.hh:449
virtual R a(int r, int i) const
Return entries of the A matrix.
Definition: instationary/onestep.hh:629
virtual bool implicit() const
Return true if method is implicit.
Definition: instationary/onestep.hh:434
virtual RealType suggestTimestep(RealType time, RealType givendt)
Return name of the scheme.
Definition: instationary/onestep.hh:695
const char * getName(int i_)
Definition: instationary/onestep.hh:1450
RK4Parameter()
Definition: instationary/onestep.hh:338
virtual R b(int r, int i) const
Return entries of the B matrix.
Definition: instationary/onestep.hh:539
Definition: instationary/onestep.hh:1435
double assembler_time
Definition: instationary/onestep.hh:747
FilenameHelper(const std::string &basename_, int i_=0)
Definition: instationary/onestep.hh:1444
Definition: instationary/onestep.hh:744
const Result & result() const
Definition: instationary/onestep.hh:829
virtual unsigned s() const
Return number of stages s of the method.
Definition: instationary/onestep.hh:134
virtual std::string name() const
Return name of the scheme.
Definition: instationary/onestep.hh:87