dune-pdelab  2.4.1
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),
34  linear_solver_iterations(0),
35  nonlinear_solver_iterations(0)
36  {}
37  };
38 
40  {
43  OneStepMethodResult() : total(), successful()
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 
121  {
122  method = &method_;
123  }
124 
132  T apply (T time, T dt, TrlV& xold, TrlV& xnew)
133  {
134  // save formatting attributes
135  ios_base_all_saver format_attribute_saver(std::cout);
136 
137  // do statistics
138  OneStepMethodPartialResult step_result;
139 
140  std::vector<TrlV*> x(1); // vector of pointers to all steps
141  x[0] = &xold; // initially we have only one
142 
143  if (verbosityLevel>=1){
144  std::ios_base::fmtflags oldflags = std::cout.flags();
145  std::cout << "TIME STEP [" << method->name() << "] "
146  << std::setw(6) << step
147  << " time (from): "
148  << std::setw(12) << std::setprecision(4) << std::scientific
149  << time
150  << " dt: "
151  << std::setw(12) << std::setprecision(4) << std::scientific
152  << dt
153  << " time (to): "
154  << std::setw(12) << std::setprecision(4) << std::scientific
155  << time+dt
156  << std::endl;
157  std::cout.flags(oldflags);
158  }
159 
160  // prepare assembler
161  igos.preStep(*method,time,dt);
162 
163  // loop over all stages
164  for (unsigned r=1; r<=method->s(); ++r)
165  {
166  if (verbosityLevel>=2){
167  std::ios_base::fmtflags oldflags = std::cout.flags();
168  std::cout << "STAGE "
169  << r
170  << " time (to): "
171  << std::setw(12) << std::setprecision(4) << std::scientific
172  << time+method->d(r)*dt
173  << "." << std::endl;
174  std::cout.flags(oldflags);
175  }
176 
177  // prepare stage
178  igos.preStage(r,x);
179 
180  // get vector for current stage
181  if (r==method->s())
182  {
183  // last stage
184  x.push_back(&xnew);
185  if (r>1) xnew = *(x[r-1]); // if r=1 then xnew has already initial guess
186  }
187  else
188  {
189  // intermediate step
190  x.push_back(new TrlV(igos.trialGridFunctionSpace()));
191  if (r>1)
192  *(x[r]) = *(x[r-1]); // use result of last stage as initial guess
193  else
194  *(x[r]) = xnew;
195  }
196 
197  // solve stage
198  try {
199  pdesolver.apply(*x[r]);
200  }
201  catch (...)
202  {
203  // time step failed -> accumulate to total only
204  PDESolverResult pderes = pdesolver.result();
205  step_result.assembler_time += pderes.assembler_time;
206  step_result.linear_solver_time += pderes.linear_solver_time;
207  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
208  step_result.nonlinear_solver_iterations += pderes.iterations;
209  res.total.assembler_time += step_result.assembler_time;
210  res.total.linear_solver_time += step_result.linear_solver_time;
211  res.total.linear_solver_iterations += step_result.linear_solver_iterations;
212  res.total.nonlinear_solver_iterations += step_result.nonlinear_solver_iterations;
213  res.total.timesteps += 1;
214  throw;
215  }
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 
222  // stage cleanup
223  igos.postStage();
224  }
225 
226  // delete intermediate steps
227  for (unsigned i=1; i<method->s(); ++i) delete x[i];
228 
229  // step cleanup
230  igos.postStep();
231 
232  // update statistics
233  res.total.assembler_time += step_result.assembler_time;
234  res.total.linear_solver_time += step_result.linear_solver_time;
235  res.total.linear_solver_iterations += step_result.linear_solver_iterations;
236  res.total.nonlinear_solver_iterations += step_result.nonlinear_solver_iterations;
237  res.total.timesteps += 1;
238  res.successful.assembler_time += step_result.assembler_time;
239  res.successful.linear_solver_time += step_result.linear_solver_time;
240  res.successful.linear_solver_iterations += step_result.linear_solver_iterations;
241  res.successful.nonlinear_solver_iterations += step_result.nonlinear_solver_iterations;
242  res.successful.timesteps += 1;
243  if (verbosityLevel>=1){
244  std::ios_base::fmtflags oldflags = std::cout.flags();
245  std::cout << "::: timesteps " << std::setw(6) << res.successful.timesteps
246  << " (" << res.total.timesteps << ")" << std::endl;
247  std::cout << "::: nl iterations " << std::setw(6) << res.successful.nonlinear_solver_iterations
248  << " (" << res.total.nonlinear_solver_iterations << ")" << std::endl;
249  std::cout << "::: lin iterations " << std::setw(6) << res.successful.linear_solver_iterations
250  << " (" << res.total.linear_solver_iterations << ")" << std::endl;
251  std::cout << "::: assemble time " << std::setw(12) << std::setprecision(4) << std::scientific
252  << res.successful.assembler_time << " (" << res.total.assembler_time << ")" << std::endl;
253  std::cout << "::: lin solve time " << std::setw(12) << std::setprecision(4) << std::scientific
254  << res.successful.linear_solver_time << " (" << res.total.linear_solver_time << ")" << std::endl;
255  std::cout.flags(oldflags);
256  }
257 
258  step++;
259  return dt;
260  }
261 
272  template<typename F>
273  T apply (T time, T dt, TrlV& xold, F& f, TrlV& xnew)
274  {
275  // do statistics
276  OneStepMethodPartialResult step_result;
277 
278  // save formatting attributes
279  ios_base_all_saver format_attribute_saver(std::cout);
280 
281  std::vector<TrlV*> x(1); // vector of pointers to all steps
282  x[0] = &xold; // initially we have only one
283 
284  if (verbosityLevel>=1){
285  std::ios_base::fmtflags oldflags = std::cout.flags();
286  std::cout << "TIME STEP [" << method->name() << "] "
287  << std::setw(6) << step
288  << " time (from): "
289  << std::setw(12) << std::setprecision(4) << std::scientific
290  << time
291  << " dt: "
292  << std::setw(12) << std::setprecision(4) << std::scientific
293  << dt
294  << " time (to): "
295  << std::setw(12) << std::setprecision(4) << std::scientific
296  << time+dt
297  << std::endl;
298  std::cout.flags(oldflags);
299  }
300 
301  // prepare assembler
302  igos.preStep(*method,time,dt);
303 
304  // loop over all stages
305  for (unsigned r=1; r<=method->s(); ++r)
306  {
307  if (verbosityLevel>=2){
308  std::ios_base::fmtflags oldflags = std::cout.flags();
309  std::cout << "STAGE "
310  << r
311  << " time (to): "
312  << std::setw(12) << std::setprecision(4) << std::scientific
313  << time+method->d(r)*dt
314  << "." << std::endl;
315  std::cout.flags(oldflags);
316  }
317 
318  // prepare stage
319  igos.preStage(r,x);
320 
321  // get vector for current stage
322  if (r==method->s())
323  {
324  // last stage
325  x.push_back(&xnew);
326  }
327  else
328  {
329  // intermediate step
330  x.push_back(new TrlV(igos.trialGridFunctionSpace()));
331  }
332 
333  // set boundary conditions and initial value
334  igos.interpolate(r,*x[r-1],f,*x[r]);
335 
336  // solve stage
337  try {
338  pdesolver.apply(*x[r]);
339  }
340  catch (...)
341  {
342  // time step failed -> accumulate to total only
343  PDESolverResult pderes = pdesolver.result();
344  step_result.assembler_time += pderes.assembler_time;
345  step_result.linear_solver_time += pderes.linear_solver_time;
346  step_result.linear_solver_iterations += pderes.linear_solver_iterations;
347  step_result.nonlinear_solver_iterations += pderes.iterations;
348  res.total.assembler_time += step_result.assembler_time;
349  res.total.linear_solver_time += step_result.linear_solver_time;
350  res.total.linear_solver_iterations += step_result.linear_solver_iterations;
351  res.total.nonlinear_solver_iterations += step_result.nonlinear_solver_iterations;
352  res.total.timesteps += 1;
353  throw;
354  }
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 
361  // stage cleanup
362  igos.postStage();
363  }
364 
365  // delete intermediate steps
366  for (unsigned i=1; i<method->s(); ++i) delete x[i];
367 
368  // step cleanup
369  igos.postStep();
370 
371  // update statistics
372  res.total.assembler_time += step_result.assembler_time;
373  res.total.linear_solver_time += step_result.linear_solver_time;
374  res.total.linear_solver_iterations += step_result.linear_solver_iterations;
375  res.total.nonlinear_solver_iterations += step_result.nonlinear_solver_iterations;
376  res.total.timesteps += 1;
377  res.successful.assembler_time += step_result.assembler_time;
378  res.successful.linear_solver_time += step_result.linear_solver_time;
379  res.successful.linear_solver_iterations += step_result.linear_solver_iterations;
380  res.successful.nonlinear_solver_iterations += step_result.nonlinear_solver_iterations;
381  res.successful.timesteps += 1;
382  if (verbosityLevel>=1){
383  std::ios_base::fmtflags oldflags = std::cout.flags();
384  std::cout << "::: timesteps " << std::setw(6) << res.successful.timesteps
385  << " (" << res.total.timesteps << ")" << std::endl;
386  std::cout << "::: nl iterations " << std::setw(6) << res.successful.nonlinear_solver_iterations
387  << " (" << res.total.nonlinear_solver_iterations << ")" << std::endl;
388  std::cout << "::: lin iterations " << std::setw(6) << res.successful.linear_solver_iterations
389  << " (" << res.total.linear_solver_iterations << ")" << std::endl;
390  std::cout << "::: assemble time " << std::setw(12) << std::setprecision(4) << std::scientific
391  << res.successful.assembler_time << " (" << res.total.assembler_time << ")" << std::endl;
392  std::cout << "::: lin solve time " << std::setw(12) << std::setprecision(4) << std::scientific
393  << res.successful.linear_solver_time << " (" << res.total.linear_solver_time << ")" << std::endl;
394  std::cout.flags(oldflags);
395  }
396 
397  step++;
398  return dt;
399  }
400 
401  private:
402  const TimeSteppingParameterInterface<T> *method;
403  IGOS& igos;
404  PDESOLVER& pdesolver;
405  int verbosityLevel;
406  int step;
407  Result res;
408  };
409 
411  } // end namespace PDELab
412 } // end namespace Dune
413 #endif
void setMethod(const TimeSteppingParameterInterface< T > &method_)
redefine the method to be used; can be done before every step
Definition: implicitonestep.hh:120
OneStepMethodResult Result
Definition: implicitonestep.hh:61
int linear_solver_iterations
Definition: implicitonestep.hh:27
OneStepMethodResult()
Definition: implicitonestep.hh:43
OneStepMethod(const TimeSteppingParameterInterface< T > &method_, IGOS &igos_, PDESOLVER &pdesolver_)
construct a new one step scheme
Definition: implicitonestep.hh:75
Do one step of a time-stepping scheme.
Definition: implicitonestep.hh:56
Definition: adaptivity.hh:27
void setVerbosityLevel(int level)
change verbosity level; 0 means completely quiet
Definition: implicitonestep.hh:84
T apply(T time, T dt, TrlV &xold, TrlV &xnew)
do one step;
Definition: implicitonestep.hh:132
double assembler_time
Definition: implicitonestep.hh:25
const Result & result() const
Definition: implicitonestep.hh:107
PDESOLVER & getPDESolver()
Access to the (non) linear solver.
Definition: implicitonestep.hh:102
double linear_solver_time
Definition: implicitonestep.hh:26
OneStepMethodPartialResult successful
Definition: implicitonestep.hh:42
OneStepMethodPartialResult()
Definition: implicitonestep.hh:30
Definition: implicitonestep.hh:39
OneStepMethodPartialResult total
Definition: implicitonestep.hh:41
void setStepNumber(int newstep)
change number of current step
Definition: implicitonestep.hh:93
unsigned int timesteps
Definition: implicitonestep.hh:24
const PDESOLVER & getPDESolver() const
Access to the (non) linear solver.
Definition: implicitonestep.hh:96
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:273
int nonlinear_solver_iterations
Definition: implicitonestep.hh:28
Definition: implicitonestep.hh:22