dune-pdelab  2.7-git
implicitonestep.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_IMPLICITONESTEP_HH
5 #define DUNE_PDELAB_INSTATIONARY_IMPLICITONESTEP_HH
6 
7 #include <iostream>
8 #include <iomanip>
9 
10 #include <dune/common/ios_state.hh>
12 
13 namespace Dune {
14  namespace PDELab {
15 
21  // Status information of Newton's method
23  {
24  unsigned int timesteps;
29 
31  timesteps(0),
32  assembler_time(0.0),
33  linear_solver_time(0.0),
36  {}
37  };
38 
40  {
44  {}
45  };
46 
48 
55  template<class T, class IGOS, class PDESOLVER, class TrlV, class TstV = TrlV>
57  {
58  typedef typename PDESOLVER::Result PDESolverResult;
59 
60  public:
62 
64 
76  IGOS& igos_, PDESOLVER& pdesolver_)
77  : method(&method_), igos(igos_), pdesolver(pdesolver_), verbosityLevel(1), step(1), res()
78  {
79  if (igos.trialGridFunctionSpace().gridView().comm().rank()>0)
80  verbosityLevel = 0;
81  }
82 
84  void setVerbosityLevel (int level)
85  {
86  if (igos.trialGridFunctionSpace().gridView().comm().rank()>0)
87  verbosityLevel = 0;
88  else
89  verbosityLevel = level;
90  }
91 
93  void setStepNumber(int newstep) { step = newstep; }
94 
96  const PDESOLVER & getPDESolver() const
97  {
98  return pdesolver;
99  }
100 
102  PDESOLVER & getPDESolver()
103  {
104  return pdesolver;
105  }
106 
107  const Result& result() const
108  {
109  return res;
110  }
111 
113 
118  void setResult (const OneStepMethodResult& result_)
119  {
120  res = result_;
122  }
123 
125 
133  {
134  method = &method_;
135  }
136 
144  T apply (T time, T dt, TrlV& xold, TrlV& xnew)
145  {
146  // save formatting attributes
147  ios_base_all_saver format_attribute_saver(std::cout);
148 
149  // do statistics
150  OneStepMethodPartialResult step_result;
151 
152  std::vector<TrlV*> x(1); // vector of pointers to all steps
153  x[0] = &xold; // initially we have only one
154 
155  if (verbosityLevel>=1){
156  std::ios_base::fmtflags oldflags = std::cout.flags();
157  std::cout << "TIME STEP [" << method->name() << "] "
158  << std::setw(6) << step
159  << " time (from): "
160  << std::setw(12) << std::setprecision(4) << std::scientific
161  << time
162  << " dt: "
163  << std::setw(12) << std::setprecision(4) << std::scientific
164  << dt
165  << " time (to): "
166  << std::setw(12) << std::setprecision(4) << std::scientific
167  << time+dt
168  << std::endl;
169  std::cout.flags(oldflags);
170  }
171 
172  // prepare assembler
173  igos.preStep(*method,time,dt);
174 
175  // loop over all stages
176  for (unsigned r=1; r<=method->s(); ++r)
177  {
178  if (verbosityLevel>=2){
179  std::ios_base::fmtflags oldflags = std::cout.flags();
180  std::cout << "STAGE "
181  << r
182  << " time (to): "
183  << std::setw(12) << std::setprecision(4) << std::scientific
184  << time+method->d(r)*dt
185  << "." << std::endl;
186  std::cout.flags(oldflags);
187  }
188 
189  // prepare stage
190  igos.preStage(r,x);
191 
192  // get vector for current stage
193  if (r==method->s())
194  {
195  // last stage
196  x.push_back(&xnew);
197  if (r>1) xnew = *(x[r-1]); // if r=1 then xnew has already initial guess
198  }
199  else
200  {
201  // intermediate step
202  x.push_back(new TrlV(igos.trialGridFunctionSpace()));
203  if (r>1)
204  *(x[r]) = *(x[r-1]); // use result of last stage as initial guess
205  else
206  *(x[r]) = xnew;
207  }
208 
209  // solve stage
210  try {
211  pdesolver.apply(*x[r]);
212  }
213  catch (...)
214  {
215  // time step failed -> accumulate to total only
216  PDESolverResult pderes = pdesolver.result();
217  step_result.assembler_time += pderes.assembler_time;
218  step_result.linear_solver_time += pderes.linear_solver_time;
219  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
220  step_result.nonlinear_solver_iterations += pderes.iterations;
221  res.total.assembler_time += step_result.assembler_time;
222  res.total.linear_solver_time += step_result.linear_solver_time;
225  res.total.timesteps += 1;
226  throw;
227  }
228  PDESolverResult pderes = pdesolver.result();
229  step_result.assembler_time += pderes.assembler_time;
230  step_result.linear_solver_time += pderes.linear_solver_time;
231  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
232  step_result.nonlinear_solver_iterations += pderes.iterations;
233 
234  // stage cleanup
235  igos.postStage();
236  }
237 
238  // delete intermediate steps
239  for (unsigned i=1; i<method->s(); ++i) delete x[i];
240 
241  // step cleanup
242  igos.postStep();
243 
244  // update statistics
245  res.total.assembler_time += step_result.assembler_time;
246  res.total.linear_solver_time += step_result.linear_solver_time;
249  res.total.timesteps += 1;
250  res.successful.assembler_time += step_result.assembler_time;
254  res.successful.timesteps += 1;
255  if (verbosityLevel>=1){
256  std::ios_base::fmtflags oldflags = std::cout.flags();
257  std::cout << "::: timesteps " << std::setw(6) << res.successful.timesteps
258  << " (" << res.total.timesteps << ")" << std::endl;
259  std::cout << "::: nl iterations " << std::setw(6) << res.successful.nonlinear_solver_iterations
260  << " (" << res.total.nonlinear_solver_iterations << ")" << std::endl;
261  std::cout << "::: lin iterations " << std::setw(6) << res.successful.linear_solver_iterations
262  << " (" << res.total.linear_solver_iterations << ")" << std::endl;
263  std::cout << "::: assemble time " << std::setw(12) << std::setprecision(4) << std::scientific
264  << res.successful.assembler_time << " (" << res.total.assembler_time << ")" << std::endl;
265  std::cout << "::: lin solve time " << std::setw(12) << std::setprecision(4) << std::scientific
266  << res.successful.linear_solver_time << " (" << res.total.linear_solver_time << ")" << std::endl;
267  std::cout.flags(oldflags);
268  }
269 
270  step++;
271  return dt;
272  }
273 
284  template<typename F>
285  T apply (T time, T dt, TrlV& xold, F& f, TrlV& xnew)
286  {
287  // do statistics
288  OneStepMethodPartialResult step_result;
289 
290  // save formatting attributes
291  ios_base_all_saver format_attribute_saver(std::cout);
292 
293  std::vector<TrlV*> x(1); // vector of pointers to all steps
294  x[0] = &xold; // initially we have only one
295 
296  if (verbosityLevel>=1){
297  std::ios_base::fmtflags oldflags = std::cout.flags();
298  std::cout << "TIME STEP [" << method->name() << "] "
299  << std::setw(6) << step
300  << " time (from): "
301  << std::setw(12) << std::setprecision(4) << std::scientific
302  << time
303  << " dt: "
304  << std::setw(12) << std::setprecision(4) << std::scientific
305  << dt
306  << " time (to): "
307  << std::setw(12) << std::setprecision(4) << std::scientific
308  << time+dt
309  << std::endl;
310  std::cout.flags(oldflags);
311  }
312 
313  // prepare assembler
314  igos.preStep(*method,time,dt);
315 
316  // loop over all stages
317  for (unsigned r=1; r<=method->s(); ++r)
318  {
319  if (verbosityLevel>=2){
320  std::ios_base::fmtflags oldflags = std::cout.flags();
321  std::cout << "STAGE "
322  << r
323  << " time (to): "
324  << std::setw(12) << std::setprecision(4) << std::scientific
325  << time+method->d(r)*dt
326  << "." << std::endl;
327  std::cout.flags(oldflags);
328  }
329 
330  // prepare stage
331  igos.preStage(r,x);
332 
333  // get vector for current stage
334  if (r==method->s())
335  {
336  // last stage
337  x.push_back(&xnew);
338  }
339  else
340  {
341  // intermediate step
342  x.push_back(new TrlV(igos.trialGridFunctionSpace()));
343  }
344 
345  // set boundary conditions and initial value
346  igos.interpolate(r,*x[r-1],f,*x[r]);
347 
348  // solve stage
349  try {
350  pdesolver.apply(*x[r]);
351  }
352  catch (...)
353  {
354  // time step failed -> accumulate to total only
355  PDESolverResult pderes = pdesolver.result();
356  step_result.assembler_time += pderes.assembler_time;
357  step_result.linear_solver_time += pderes.linear_solver_time;
358  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
359  step_result.nonlinear_solver_iterations += pderes.iterations;
360  res.total.assembler_time += step_result.assembler_time;
361  res.total.linear_solver_time += step_result.linear_solver_time;
364  res.total.timesteps += 1;
365  throw;
366  }
367  PDESolverResult pderes = pdesolver.result();
368  step_result.assembler_time += pderes.assembler_time;
369  step_result.linear_solver_time += pderes.linear_solver_time;
370  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
371  step_result.nonlinear_solver_iterations += pderes.iterations;
372 
373  // stage cleanup
374  igos.postStage();
375  }
376 
377  // delete intermediate steps
378  for (unsigned i=1; i<method->s(); ++i) delete x[i];
379 
380  // step cleanup
381  igos.postStep();
382 
383  // update statistics
384  res.total.assembler_time += step_result.assembler_time;
385  res.total.linear_solver_time += step_result.linear_solver_time;
388  res.total.timesteps += 1;
389  res.successful.assembler_time += step_result.assembler_time;
393  res.successful.timesteps += 1;
394  if (verbosityLevel>=1){
395  std::ios_base::fmtflags oldflags = std::cout.flags();
396  std::cout << "::: timesteps " << std::setw(6) << res.successful.timesteps
397  << " (" << res.total.timesteps << ")" << std::endl;
398  std::cout << "::: nl iterations " << std::setw(6) << res.successful.nonlinear_solver_iterations
399  << " (" << res.total.nonlinear_solver_iterations << ")" << std::endl;
400  std::cout << "::: lin iterations " << std::setw(6) << res.successful.linear_solver_iterations
401  << " (" << res.total.linear_solver_iterations << ")" << std::endl;
402  std::cout << "::: assemble time " << std::setw(12) << std::setprecision(4) << std::scientific
403  << res.successful.assembler_time << " (" << res.total.assembler_time << ")" << std::endl;
404  std::cout << "::: lin solve time " << std::setw(12) << std::setprecision(4) << std::scientific
405  << res.successful.linear_solver_time << " (" << res.total.linear_solver_time << ")" << std::endl;
406  std::cout.flags(oldflags);
407  }
408 
409  step++;
410  return dt;
411  }
412 
413  private:
414  const TimeSteppingParameterInterface<T> *method;
415  IGOS& igos;
416  PDESOLVER& pdesolver;
417  int verbosityLevel;
418  int step;
419  Result res;
420  };
421 
423  } // end namespace PDELab
424 } // end namespace Dune
425 #endif // DUNE_PDELAB_INSTATIONARY_IMPLICITONESTEP_HH
Dune::PDELab::OneStepMethod::apply
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: implicitonestep.hh:285
Dune::PDELab::OneStepMethodResult::total
OneStepMethodPartialResult total
Definition: implicitonestep.hh:41
Dune::PDELab::OneStepMethodPartialResult::nonlinear_solver_iterations
int nonlinear_solver_iterations
Definition: implicitonestep.hh:28
Dune::PDELab::OneStepMethod
Do one step of a time-stepping scheme.
Definition: implicitonestep.hh:56
Dune::PDELab::OneStepMethodPartialResult::linear_solver_time
double linear_solver_time
Definition: implicitonestep.hh:26
Dune::PDELab::OneStepMethodPartialResult
Definition: implicitonestep.hh:22
Dune
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Dune::PDELab::OneStepMethod::apply
T apply(T time, T dt, TrlV &xold, TrlV &xnew)
do one step;
Definition: implicitonestep.hh:144
Dune::PDELab::OneStepMethod::setStepNumber
void setStepNumber(int newstep)
change number of current step
Definition: implicitonestep.hh:93
Dune::PDELab::OneStepMethodResult::successful
OneStepMethodPartialResult successful
Definition: implicitonestep.hh:42
Dune::PDELab::OneStepMethod::getPDESolver
const PDESOLVER & getPDESolver() const
Access to the (non) linear solver.
Definition: implicitonestep.hh:96
Dune::PDELab::OneStepMethod::getPDESolver
PDESOLVER & getPDESolver()
Access to the (non) linear solver.
Definition: implicitonestep.hh:102
Dune::PDELab::OneStepMethodResult::OneStepMethodResult
OneStepMethodResult()
Definition: implicitonestep.hh:43
Dune::PDELab::TimeSteppingParameterInterface::d
virtual R d(int r) const =0
Return entries of the d Vector.
Dune::PDELab::OneStepMethod::OneStepMethod
OneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, PDESOLVER &pdesolver_)
construct a new one step scheme
Definition: implicitonestep.hh:75
Dune::PDELab::OneStepMethodPartialResult::timesteps
unsigned int timesteps
Definition: implicitonestep.hh:24
Dune::PDELab::OneStepMethod::setVerbosityLevel
void setVerbosityLevel(int level)
change verbosity level; 0 means completely quiet
Definition: implicitonestep.hh:84
Dune::PDELab::OneStepMethodPartialResult::assembler_time
double assembler_time
Definition: implicitonestep.hh:25
Dune::PDELab::OneStepMethod::Result
OneStepMethodResult Result
Definition: implicitonestep.hh:61
Dune::PDELab::TimeSteppingParameterInterface::s
virtual unsigned s() const =0
Return number of stages of the method.
Dune::PDELab::OneStepMethodResult
Definition: implicitonestep.hh:39
Dune::PDELab::TimeSteppingParameterInterface::name
virtual std::string name() const =0
Return name of the scheme.
Dune::PDELab::OneStepMethod::result
const Result & result() const
Definition: implicitonestep.hh:107
Dune::PDELab::TimeSteppingParameterInterface< T >
Dune::PDELab::OneStepMethod::setMethod
void setMethod(const TimeSteppingParameterInterface< T > &method_)
redefine the method to be used; can be done before every step
Definition: implicitonestep.hh:132
Dune::PDELab::OneStepMethodPartialResult::OneStepMethodPartialResult
OneStepMethodPartialResult()
Definition: implicitonestep.hh:30
onestepparameter.hh
Dune::PDELab::OneStepMethodPartialResult::linear_solver_iterations
int linear_solver_iterations
Definition: implicitonestep.hh:27
Dune::PDELab::OneStepMethod::setResult
void setResult(const OneStepMethodResult &result_)
Set a new result.
Definition: implicitonestep.hh:118