dune-pdelab  2.0.0
method.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_MULTISTEP_METHOD_HH
5 #define DUNE_PDELAB_MULTISTEP_METHOD_HH
6 
7 #include <iomanip>
8 #include <iostream>
9 
10 #include <dune/common/ios_state.hh>
11 #include <dune/common/shared_ptr.hh>
12 #include <dune/common/timer.hh>
13 
15 
16 namespace Dune {
17  namespace PDELab {
18 
21 
23 
30  template<class T, class MGOS, class PDESolver, class TrialV,
31  class TestV = TrialV>
34 
35  const Parameters* parameters;
36  MGOS& mgos;
37  PDESolver& pdeSolver;
38  unsigned verbosity;
39  unsigned step;
40 
41  public:
43 
52  MultiStepMethod(const Parameters& parameters_,
53  MGOS& mgos_, PDESolver& pdeSolver_)
54  : parameters(&parameters_), mgos(mgos_), pdeSolver(pdeSolver_),
55  verbosity(1), step(1)
56  {
57  if(mgos.trialGridFunctionSpace().gridView().comm().rank() > 0)
58  verbosity = 0;
59  }
60 
62  void setVerbosityLevel (int level)
63  {
64  if(mgos.trialGridFunctionSpace().gridView().comm().rank() == 0)
65  verbosity = level;
66  }
67 
69  void setMethod(const Parameters &parameters_) {
70  parameters = &parameters_;
71  }
72 
74  /*
75  * \param time Start of time step
76  * \param dt Suggested time step size
77  * \param oldValues Vector of pointers to the old values. Must support
78  * the expression *oldvalues[i], which should yield a
79  * reference to the old value at time time-i*dt.
80  * \param xnew Where to store the new value.
81  *
82  * \return Actual time step size
83  */
84  template<class OldValues>
85  T apply(T time, T dt, const OldValues& oldValues, TrialV& xnew)
86  {
87  Timer allTimer;
88  Timer subTimer;
89 
90  // save formatting attributes
91  ios_base_all_saver format_attribute_saver(std::cout);
92 
93  if(verbosity >= 1)
94  std::cout << "TIME STEP [" << parameters->name() << "] "
95  << std::setw(6) << step
96  << " time (from): "
97  << std::setw(12) << std::setprecision(4) << std::scientific
98  << time
99  << " dt: "
100  << std::setw(12) << std::setprecision(4) << std::scientific
101  << dt
102  << " time (to): "
103  << std::setw(12) << std::setprecision(4) << std::scientific
104  << time+dt
105  << std::endl;
106 
107  // prepare assembler
108  if(verbosity >= 2) {
109  std::cout << "== prepare assembler" << std::endl;
110  subTimer.reset();
111  }
112  mgos.preStep(time,dt, oldValues);
113  if(verbosity >= 2)
114  std::cout << "== prepare assembler (" << subTimer.elapsed() << "s)"
115  << std::endl;
116 
117  // solve
118  if(verbosity >= 2) {
119  std::cout << "== apply solver" << std::endl;
120  subTimer.reset();
121  }
122  pdeSolver.apply(xnew);
123  if(verbosity >= 2)
124  std::cout << "== apply solver (" << subTimer.elapsed() << "s)"
125  << std::endl;
126 
127  // postprocessing in the assembler
128  if(verbosity >= 2) {
129  std::cout << "== cleanup assembler" << std::endl;
130  subTimer.reset();
131  }
132  mgos.postStep();
133  if(verbosity >= 2)
134  std::cout << "== cleanup assembler (" << subTimer.elapsed() << "s)"
135  << std::endl;
136 
137  ++step;
138 
139  if(verbosity >= 2)
140  std::cout << "== time step done (" << allTimer.elapsed() << "s)"
141  << std::endl;
142 
143  return dt;
144  }
145 
147  /* This is a version which interpolates constraints at the start of each
148  * stage
149  *
150  * \param time Start of time step
151  * \param dt Suggested time step size
152  * \param oldValues Vector of pointers to the old values. Must support
153  * the expression *oldValues[i], which should yield a
154  * reference to the old value at time time-i*dt.
155  * \param f Function to interpolate boundary conditions from.
156  * Should support the method setTime().
157  * \param xnew Where to store the new value.
158  *
159  * \return Actual time step size
160  */
161  template<typename OldValues, typename F>
162  T apply(T time, T dt, const OldValues& oldValues, F& f, TrialV& xnew)
163  {
164  Timer allTimer;
165  Timer subTimer;
166 
167  // save formatting attributes
168  ios_base_all_saver format_attribute_saver(std::cout);
169 
170  if(verbosity >= 1)
171  std::cout << "TIME STEP [" << parameters->name() << "] "
172  << std::setw(6) << step
173  << " time (from): "
174  << std::setw(12) << std::setprecision(4) << std::scientific
175  << time
176  << " dt: "
177  << std::setw(12) << std::setprecision(4) << std::scientific
178  << dt
179  << " time (to): "
180  << std::setw(12) << std::setprecision(4) << std::scientific
181  << time+dt
182  << std::endl;
183 
184  // prepare assembler
185  if(verbosity >= 2) {
186  std::cout << "== prepare assembler" << std::endl;
187  subTimer.reset();
188  }
189  mgos.preStep(time, dt, oldValues);
190  if(verbosity >= 2)
191  std::cout << "== prepare assembler (" << subTimer.elapsed() << "s)"
192  << std::endl;
193 
194  // set boundary conditions and initial value
195  if(verbosity >= 2) {
196  std::cout << "== setup result vector" << std::endl;
197  subTimer.reset();
198  }
199  f.setTime(time+dt);
200  mgos.interpolate(*oldValues[0],f,xnew);
201  if(verbosity >= 2)
202  std::cout << "== setup result vector (" << subTimer.elapsed() << "s)"
203  << std::endl;
204 
205  // solve stage
206  if(verbosity >= 2) {
207  std::cout << "== apply solver" << std::endl;
208  subTimer.reset();
209  }
210  pdeSolver.apply(xnew);
211  if(verbosity >= 2)
212  std::cout << "== apply solver (" << subTimer.elapsed() << "s)"
213  << std::endl;
214 
215  // postprocessing in the assembler
216  if(verbosity >= 2) {
217  std::cout << "== cleanup assembler" << std::endl;
218  subTimer.reset();
219  }
220  mgos.postStep();
221  if(verbosity >= 2)
222  std::cout << "== cleanup assembler (" << subTimer.elapsed() << "s)"
223  << std::endl;
224 
225  step++;
226 
227  if(verbosity >= 2)
228  std::cout << "== time step done (" << allTimer.elapsed() << "s)"
229  << std::endl;
230 
231  return dt;
232  }
233 
235 
244  shared_ptr<const TrialV> apply(T time, T dt)
245  {
246  Timer allTimer;
247  Timer subTimer;
248 
249  // save formatting attributes
250  ios_base_all_saver format_attribute_saver(std::cout);
251 
252  if(verbosity >= 1)
253  std::cout << "TIME STEP [" << parameters->name() << "] "
254  << std::setw(6) << step
255  << " time (from): "
256  << std::setw(12) << std::setprecision(4) << std::scientific
257  << time
258  << " dt: "
259  << std::setw(12) << std::setprecision(4) << std::scientific
260  << dt
261  << " time (to): "
262  << std::setw(12) << std::setprecision(4) << std::scientific
263  << time+dt
264  << std::endl;
265 
266  // prepare assembler
267  if(verbosity >= 2) {
268  std::cout << "== prepare assembler" << std::endl;
269  subTimer.reset();
270  }
271  mgos.preStep(step, time, dt);
272  if(verbosity >= 2)
273  std::cout << "== prepare assembler (" << subTimer.elapsed() << "s)"
274  << std::endl;
275 
276  // create vector, using last time step as start
277  if(verbosity >= 2) {
278  std::cout << "== setup result vector" << std::endl;
279  subTimer.reset();
280  }
281  shared_ptr<TrialV>
282  xnew(new TrialV(*mgos.getCache()->getUnknowns(step-1)));
283  if(verbosity >= 2)
284  std::cout << "== setup result vector (" << subTimer.elapsed() << "s)"
285  << std::endl;
286 
287  // solve
288  if(verbosity >= 2) {
289  std::cout << "== apply solver" << std::endl;
290  subTimer.reset();
291  }
292  pdeSolver.apply(*xnew);
293  if(verbosity >= 2)
294  std::cout << "== apply solver (" << subTimer.elapsed() << "s)"
295  << std::endl;
296 
297  // postprocessing in the assembler
298  if(verbosity >= 2) {
299  std::cout << "== cleanup assembler" << std::endl;
300  subTimer.reset();
301  }
302  mgos.postStep();
303  if(verbosity >= 2)
304  std::cout << "== cleanup assembler (" << subTimer.elapsed() << "s)"
305  << std::endl;
306 
307  // store result for next step
308  mgos.getCache()->setUnknowns(step, xnew);
309 
310  ++step;
311 
312  if(verbosity >= 2)
313  std::cout << "== time step done (" << allTimer.elapsed() << "s)"
314  << std::endl;
315 
316  return xnew;
317  }
318 
320 
334  template<typename F>
335  shared_ptr<const TrialV> apply(T time, T dt, F& f)
336  {
337  Timer allTimer;
338  Timer subTimer;
339 
340  // save formatting attributes
341  ios_base_all_saver format_attribute_saver(std::cout);
342 
343  if(verbosity >= 1)
344  std::cout << "TIME STEP [" << parameters->name() << "] "
345  << std::setw(6) << step
346  << " time (from): "
347  << std::setw(12) << std::setprecision(4) << std::scientific
348  << time
349  << " dt: "
350  << std::setw(12) << std::setprecision(4) << std::scientific
351  << dt
352  << " time (to): "
353  << std::setw(12) << std::setprecision(4) << std::scientific
354  << time+dt
355  << std::endl;
356 
357  // prepare assembler
358  if(verbosity >= 2) {
359  std::cout << "== prepare assembler" << std::endl;
360  subTimer.reset();
361  }
362  mgos.preStep(step, time, dt);
363  if(verbosity >= 2)
364  std::cout << "== prepare assembler (" << subTimer.elapsed() << "s)"
365  << std::endl;
366 
367  // setup vector
368  if(verbosity >= 2) {
369  std::cout << "== setup result vector" << std::endl;
370  subTimer.reset();
371  }
372  shared_ptr<TrialV> xnew(new TrialV(mgos.trialGridFunctionSpace()));
373  // set boundary conditions and initial value
374  f.setTime(time+dt);
375  mgos.interpolate(*mgos.getCache()->getUnknowns(step-1),f,*xnew);
376  if(verbosity >= 2)
377  std::cout << "== setup result vector (" << subTimer.elapsed() << "s)"
378  << std::endl;
379 
380  // solve stage
381  if(verbosity >= 2) {
382  std::cout << "== apply solver" << std::endl;
383  subTimer.reset();
384  }
385  pdeSolver.apply(*xnew);
386  if(verbosity >= 2)
387  std::cout << "== apply solver (" << subTimer.elapsed() << "s)"
388  << std::endl;
389 
390  // postprocessing in the assembler
391  if(verbosity >= 2) {
392  std::cout << "== cleanup assembler" << std::endl;
393  subTimer.reset();
394  }
395  mgos.postStep();
396  if(verbosity >= 2)
397  std::cout << "== cleanup assembler (" << subTimer.elapsed() << "s)"
398  << std::endl;
399 
400  // store result for next step
401  mgos.getCache()->setUnknowns(step, xnew);
402 
403  ++step;
404 
405  if(verbosity >= 2)
406  std::cout << "== time step done (" << allTimer.elapsed() << "s)"
407  << std::endl;
408 
409  return xnew;
410  }
411  };
412 
414  } // namespace PDELab
415 } // namespace Dune
416 
417 #endif // DUNE_PDELAB_MULTISTEP_METHOD_HH
Do one step of a multi-step time-stepping scheme.
Definition: method.hh:32
shared_ptr< const TrialV > apply(T time, T dt, F &f)
do one step (with caching)
Definition: method.hh:335
T apply(T time, T dt, const OldValues &oldValues, F &f, TrialV &xnew)
do one step;
Definition: method.hh:162
T apply(T time, T dt, const OldValues &oldValues, TrialV &xnew)
do one step;
Definition: method.hh:85
void setVerbosityLevel(int level)
change verbosity level; 0 means completely quiet
Definition: method.hh:62
virtual std::string name() const =0
Return name of the scheme.
shared_ptr< const TrialV > apply(T time, T dt)
do one step (with caching)
Definition: method.hh:244
const F & f
Definition: common/constraints.hh:145
void setMethod(const Parameters &parameters_)
redefine the method to be used; can be done before every step
Definition: method.hh:69
MultiStepMethod(const Parameters &parameters_, MGOS &mgos_, PDESolver &pdeSolver_)
construct a new multi-step scheme
Definition: method.hh:52
Base parameter class for multi step time schemes.
Definition: parameter.hh:26