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 
227  // delete intermediate steps
228  for (unsigned i=1; i<r; ++i) delete x[i];
229  if (r < method->s())
230  delete x[r];
231 
232  throw;
233  }
234  PDESolverResult pderes = pdesolver.result();
235  step_result.assembler_time += pderes.assembler_time;
236  step_result.linear_solver_time += pderes.linear_solver_time;
237  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
238  step_result.nonlinear_solver_iterations += pderes.iterations;
239 
240  // stage cleanup
241  igos.postStage();
242  }
243 
244  // delete intermediate steps
245  for (unsigned i=1; i<method->s(); ++i) delete x[i];
246 
247  // step cleanup
248  igos.postStep();
249 
250  // update statistics
251  res.total.assembler_time += step_result.assembler_time;
252  res.total.linear_solver_time += step_result.linear_solver_time;
255  res.total.timesteps += 1;
256  res.successful.assembler_time += step_result.assembler_time;
260  res.successful.timesteps += 1;
261  if (verbosityLevel>=1){
262  std::ios_base::fmtflags oldflags = std::cout.flags();
263  std::cout << "::: timesteps " << std::setw(6) << res.successful.timesteps
264  << " (" << res.total.timesteps << ")" << std::endl;
265  std::cout << "::: nl iterations " << std::setw(6) << res.successful.nonlinear_solver_iterations
266  << " (" << res.total.nonlinear_solver_iterations << ")" << std::endl;
267  std::cout << "::: lin iterations " << std::setw(6) << res.successful.linear_solver_iterations
268  << " (" << res.total.linear_solver_iterations << ")" << std::endl;
269  std::cout << "::: assemble time " << std::setw(12) << std::setprecision(4) << std::scientific
270  << res.successful.assembler_time << " (" << res.total.assembler_time << ")" << std::endl;
271  std::cout << "::: lin solve time " << std::setw(12) << std::setprecision(4) << std::scientific
272  << res.successful.linear_solver_time << " (" << res.total.linear_solver_time << ")" << std::endl;
273  std::cout.flags(oldflags);
274  }
275 
276  step++;
277  return dt;
278  }
279 
290  template<typename F>
291  T apply (T time, T dt, TrlV& xold, F& f, TrlV& xnew)
292  {
293  // do statistics
294  OneStepMethodPartialResult step_result;
295 
296  // save formatting attributes
297  ios_base_all_saver format_attribute_saver(std::cout);
298 
299  std::vector<TrlV*> x(1); // vector of pointers to all steps
300  x[0] = &xold; // initially we have only one
301 
302  if (verbosityLevel>=1){
303  std::ios_base::fmtflags oldflags = std::cout.flags();
304  std::cout << "TIME STEP [" << method->name() << "] "
305  << std::setw(6) << step
306  << " time (from): "
307  << std::setw(12) << std::setprecision(4) << std::scientific
308  << time
309  << " dt: "
310  << std::setw(12) << std::setprecision(4) << std::scientific
311  << dt
312  << " time (to): "
313  << std::setw(12) << std::setprecision(4) << std::scientific
314  << time+dt
315  << std::endl;
316  std::cout.flags(oldflags);
317  }
318 
319  // prepare assembler
320  igos.preStep(*method,time,dt);
321 
322  // loop over all stages
323  for (unsigned r=1; r<=method->s(); ++r)
324  {
325  if (verbosityLevel>=2){
326  std::ios_base::fmtflags oldflags = std::cout.flags();
327  std::cout << "STAGE "
328  << r
329  << " time (to): "
330  << std::setw(12) << std::setprecision(4) << std::scientific
331  << time+method->d(r)*dt
332  << "." << std::endl;
333  std::cout.flags(oldflags);
334  }
335 
336  // prepare stage
337  igos.preStage(r,x);
338 
339  // get vector for current stage
340  if (r==method->s())
341  {
342  // last stage
343  x.push_back(&xnew);
344  }
345  else
346  {
347  // intermediate step
348  x.push_back(new TrlV(igos.trialGridFunctionSpace()));
349  }
350 
351  // set boundary conditions and initial value
352  igos.interpolate(r,*x[r-1],f,*x[r]);
353 
354  // solve stage
355  try {
356  pdesolver.apply(*x[r]);
357  }
358  catch (...)
359  {
360  // time step failed -> accumulate to total only
361  PDESolverResult pderes = pdesolver.result();
362  step_result.assembler_time += pderes.assembler_time;
363  step_result.linear_solver_time += pderes.linear_solver_time;
364  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
365  step_result.nonlinear_solver_iterations += pderes.iterations;
366  res.total.assembler_time += step_result.assembler_time;
367  res.total.linear_solver_time += step_result.linear_solver_time;
370  res.total.timesteps += 1;
371 
372  // delete intermediate steps
373  for (unsigned i=1; i<r; ++i) delete x[i];
374  if (r < method->s())
375  delete x[r];
376 
377  throw;
378  }
379  PDESolverResult pderes = pdesolver.result();
380  step_result.assembler_time += pderes.assembler_time;
381  step_result.linear_solver_time += pderes.linear_solver_time;
382  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
383  step_result.nonlinear_solver_iterations += pderes.iterations;
384 
385  // stage cleanup
386  igos.postStage();
387  }
388 
389  // delete intermediate steps
390  for (unsigned i=1; i<method->s(); ++i) delete x[i];
391 
392  // step cleanup
393  igos.postStep();
394 
395  // update statistics
396  res.total.assembler_time += step_result.assembler_time;
397  res.total.linear_solver_time += step_result.linear_solver_time;
400  res.total.timesteps += 1;
401  res.successful.assembler_time += step_result.assembler_time;
405  res.successful.timesteps += 1;
406  if (verbosityLevel>=1){
407  std::ios_base::fmtflags oldflags = std::cout.flags();
408  std::cout << "::: timesteps " << std::setw(6) << res.successful.timesteps
409  << " (" << res.total.timesteps << ")" << std::endl;
410  std::cout << "::: nl iterations " << std::setw(6) << res.successful.nonlinear_solver_iterations
411  << " (" << res.total.nonlinear_solver_iterations << ")" << std::endl;
412  std::cout << "::: lin iterations " << std::setw(6) << res.successful.linear_solver_iterations
413  << " (" << res.total.linear_solver_iterations << ")" << std::endl;
414  std::cout << "::: assemble time " << std::setw(12) << std::setprecision(4) << std::scientific
415  << res.successful.assembler_time << " (" << res.total.assembler_time << ")" << std::endl;
416  std::cout << "::: lin solve time " << std::setw(12) << std::setprecision(4) << std::scientific
417  << res.successful.linear_solver_time << " (" << res.total.linear_solver_time << ")" << std::endl;
418  std::cout.flags(oldflags);
419  }
420 
421  step++;
422  return dt;
423  }
424 
425  private:
426  const TimeSteppingParameterInterface<T> *method;
427  IGOS& igos;
428  PDESOLVER& pdesolver;
429  int verbosityLevel;
430  int step;
431  Result res;
432  };
433 
435  } // end namespace PDELab
436 } // end namespace Dune
437 #endif // DUNE_PDELAB_INSTATIONARY_IMPLICITONESTEP_HH
const std::string s
Definition: function.hh:843
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.
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Definition: implicitonestep.hh:23
int linear_solver_iterations
Definition: implicitonestep.hh:27
double linear_solver_time
Definition: implicitonestep.hh:26
int nonlinear_solver_iterations
Definition: implicitonestep.hh:28
unsigned int timesteps
Definition: implicitonestep.hh:24
OneStepMethodPartialResult()
Definition: implicitonestep.hh:30
double assembler_time
Definition: implicitonestep.hh:25
Definition: implicitonestep.hh:40
OneStepMethodResult()
Definition: implicitonestep.hh:43
OneStepMethodPartialResult total
Definition: implicitonestep.hh:41
OneStepMethodPartialResult successful
Definition: implicitonestep.hh:42
Do one step of a time-stepping scheme.
Definition: implicitonestep.hh:57
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:291
const PDESOLVER & getPDESolver() const
Access to the (non) linear solver.
Definition: implicitonestep.hh:96
void setStepNumber(int newstep)
change number of current step
Definition: implicitonestep.hh:93
OneStepMethodResult Result
Definition: implicitonestep.hh:61
void setMethod(const TimeSteppingParameterInterface< T > &method_)
redefine the method to be used; can be done before every step
Definition: implicitonestep.hh:132
void setResult(const OneStepMethodResult &result_)
Set a new result.
Definition: implicitonestep.hh:118
T apply(T time, T dt, TrlV &xold, TrlV &xnew)
do one step;
Definition: implicitonestep.hh:144
OneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, PDESOLVER &pdesolver_)
construct a new one step scheme
Definition: implicitonestep.hh:75
const Result & result() const
Definition: implicitonestep.hh:107
PDESOLVER & getPDESolver()
Access to the (non) linear solver.
Definition: implicitonestep.hh:102
void setVerbosityLevel(int level)
change verbosity level; 0 means completely quiet
Definition: implicitonestep.hh:84