DUNE PDELab (git)

implicitonestep.hh
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
11#include <dune/pdelab/instationary/onestepparameter.hh>
12
13namespace Dune {
14 namespace PDELab {
15
21 // Status information of Newton's method
22 struct OneStepMethodPartialResult
23 {
24 unsigned int timesteps;
25 double assembler_time;
26 double linear_solver_time;
27 int linear_solver_iterations;
28 int nonlinear_solver_iterations;
29
30 OneStepMethodPartialResult() :
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
39 struct OneStepMethodResult
40 {
41 OneStepMethodPartialResult total;
42 OneStepMethodPartialResult successful;
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:
61 typedef OneStepMethodResult Result;
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_;
121 setStepNumber(res.successful.timesteps+1);
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;
223 res.total.linear_solver_iterations += step_result.linear_solver_iterations;
224 res.total.nonlinear_solver_iterations += step_result.nonlinear_solver_iterations;
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;
253 res.total.linear_solver_iterations += step_result.linear_solver_iterations;
254 res.total.nonlinear_solver_iterations += step_result.nonlinear_solver_iterations;
255 res.total.timesteps += 1;
256 res.successful.assembler_time += step_result.assembler_time;
257 res.successful.linear_solver_time += step_result.linear_solver_time;
258 res.successful.linear_solver_iterations += step_result.linear_solver_iterations;
259 res.successful.nonlinear_solver_iterations += step_result.nonlinear_solver_iterations;
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 initial value from xnew or last stage
352 TrlV const * const init_guess = (r==1) ? &xnew : x[r-1];
353 // set boundary conditions
354 igos.interpolate(r,*init_guess,f,*x[r]);
355
356 // solve stage
357 try {
358 pdesolver.apply(*x[r]);
359 }
360 catch (...)
361 {
362 // time step failed -> accumulate to total only
363 PDESolverResult pderes = pdesolver.result();
364 step_result.assembler_time += pderes.assembler_time;
365 step_result.linear_solver_time += pderes.linear_solver_time;
366 step_result.linear_solver_iterations += pderes.linear_solver_iterations;
367 step_result.nonlinear_solver_iterations += pderes.iterations;
368 res.total.assembler_time += step_result.assembler_time;
369 res.total.linear_solver_time += step_result.linear_solver_time;
370 res.total.linear_solver_iterations += step_result.linear_solver_iterations;
371 res.total.nonlinear_solver_iterations += step_result.nonlinear_solver_iterations;
372 res.total.timesteps += 1;
373
374 // delete intermediate steps
375 for (unsigned i=1; i<r; ++i) delete x[i];
376 if (r < method->s())
377 delete x[r];
378
379 throw;
380 }
381 PDESolverResult pderes = pdesolver.result();
382 step_result.assembler_time += pderes.assembler_time;
383 step_result.linear_solver_time += pderes.linear_solver_time;
384 step_result.linear_solver_iterations += pderes.linear_solver_iterations;
385 step_result.nonlinear_solver_iterations += pderes.iterations;
386
387 // stage cleanup
388 igos.postStage();
389 }
390
391 // delete intermediate steps
392 for (unsigned i=1; i<method->s(); ++i) delete x[i];
393
394 // step cleanup
395 igos.postStep();
396
397 // update statistics
398 res.total.assembler_time += step_result.assembler_time;
399 res.total.linear_solver_time += step_result.linear_solver_time;
400 res.total.linear_solver_iterations += step_result.linear_solver_iterations;
401 res.total.nonlinear_solver_iterations += step_result.nonlinear_solver_iterations;
402 res.total.timesteps += 1;
403 res.successful.assembler_time += step_result.assembler_time;
404 res.successful.linear_solver_time += step_result.linear_solver_time;
405 res.successful.linear_solver_iterations += step_result.linear_solver_iterations;
406 res.successful.nonlinear_solver_iterations += step_result.nonlinear_solver_iterations;
407 res.successful.timesteps += 1;
408 if (verbosityLevel>=1){
409 std::ios_base::fmtflags oldflags = std::cout.flags();
410 std::cout << "::: timesteps " << std::setw(6) << res.successful.timesteps
411 << " (" << res.total.timesteps << ")" << std::endl;
412 std::cout << "::: nl iterations " << std::setw(6) << res.successful.nonlinear_solver_iterations
413 << " (" << res.total.nonlinear_solver_iterations << ")" << std::endl;
414 std::cout << "::: lin iterations " << std::setw(6) << res.successful.linear_solver_iterations
415 << " (" << res.total.linear_solver_iterations << ")" << std::endl;
416 std::cout << "::: assemble time " << std::setw(12) << std::setprecision(4) << std::scientific
417 << res.successful.assembler_time << " (" << res.total.assembler_time << ")" << std::endl;
418 std::cout << "::: lin solve time " << std::setw(12) << std::setprecision(4) << std::scientific
419 << res.successful.linear_solver_time << " (" << res.total.linear_solver_time << ")" << std::endl;
420 std::cout.flags(oldflags);
421 }
422
423 step++;
424 return dt;
425 }
426
427 private:
429 IGOS& igos;
430 PDESOLVER& pdesolver;
431 int verbosityLevel;
432 int step;
433 Result res;
434 };
435
437 } // end namespace PDELab
438} // end namespace Dune
439#endif // DUNE_PDELAB_INSTATIONARY_IMPLICITONESTEP_HH
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
void setStepNumber(int newstep)
change number of current step
Definition: implicitonestep.hh:93
void setMethod(const TimeSteppingParameterInterface< T > &method_)
redefine the method to be used; can be done before every step
Definition: implicitonestep.hh:132
PDESOLVER & getPDESolver()
Access to the (non) linear solver.
Definition: implicitonestep.hh:102
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 PDESOLVER & getPDESolver() const
Access to the (non) linear solver.
Definition: implicitonestep.hh:96
void setVerbosityLevel(int level)
change verbosity level; 0 means completely quiet
Definition: implicitonestep.hh:84
Utility class for storing and resetting stream attributes.
Definition: ios_state.hh:34
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.
Utility class for storing and resetting stream attributes.
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)