Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
theta_time_step.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 * Author: Andreas Vogel
4 *
5 * This file is part of UG4.
6 *
7 * UG4 is free software: you can redistribute it and/or modify it under the
8 * terms of the GNU Lesser General Public License version 3 (as published by the
9 * Free Software Foundation) with the following additional attribution
10 * requirements (according to LGPL/GPL v3 §7):
11 *
12 * (1) The following notice must be displayed in the Appropriate Legal Notices
13 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 *
15 * (2) The following notice must be displayed at a prominent place in the
16 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 *
18 * (3) The following bibliography is recommended for citation and must be
19 * preserved in all covered files:
20 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 * parallel geometric multigrid solver on hierarchically distributed grids.
22 * Computing and visualization in science 16, 4 (2013), 151-164"
23 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 * flexible software system for simulating pde based models on high performance
25 * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 *
27 * This program is distributed in the hope that it will be useful,
28 * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 * GNU Lesser General Public License for more details.
31 */
32
33#ifndef __H__UG__LIB_DISC__TIME_DISC__THETA_TIME_STEP__
34#define __H__UG__LIB_DISC__TIME_DISC__THETA_TIME_STEP__
35
36// extern libraries
37#include <deque>
38#include <cmath>
39
40// other ug libraries
42#include "common/common.h"
43
44// module-intern libraries
47
48namespace ug{
49
52
54template <typename TAlgebra>
56 : public ITimeDiscretization<TAlgebra>
57{
58 public:
60 typedef TAlgebra algebra_type;
61
63 typedef typename algebra_type::matrix_type matrix_type;
64
66 typedef typename algebra_type::vector_type vector_type;
67
70
73
74 public:
80
82
84 virtual size_t num_prev_steps() const {return m_prevSteps;}
85
88 number dt);
89
92 number dt, const GridLevel& gl);
93
96
99 const GridLevel& gl);
100
101 virtual number future_time() const {return m_futureTime;}
102
103 public:
104 void assemble_jacobian(matrix_type& J, const vector_type& u, const GridLevel& gl);
105
106 void assemble_defect(vector_type& d, const vector_type& u, const GridLevel& gl);
107
108 void assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl);
109
110 void assemble_rhs(vector_type& b, const vector_type& u, const GridLevel& gl);
111
112 void assemble_rhs(vector_type& b, const GridLevel& gl);
113
114 void adjust_solution(vector_type& u, const GridLevel& gl);
115
118
120 void calc_error(const vector_type& u, error_vector_type* u_vtk);
121 void calc_error(const vector_type& u) {calc_error(u, NULL);};
122 void calc_error(const vector_type& u, error_vector_type& u_vtk){calc_error(u, &u_vtk);};
123
126 void invalidate_error() {this->m_spDomDisc->invalidate_error();};
127
129 bool is_error_valid() {return this->m_spDomDisc->is_error_valid();};
130
133
134 protected:
136 virtual number update_scaling(std::vector<number>& vSM,
137 std::vector<number>& vSA,
138 number dt, number currentTime,
140
141 size_t m_prevSteps;
142 std::vector<number> m_vScaleMass;
143 std::vector<number> m_vScaleStiff;
144
148};
149
151
164template <typename TAlgebra>
166 : public MultiStepTimeDiscretization<TAlgebra>
167{
168 public:
171
173 typedef TAlgebra algebra_type;
174
176 typedef typename algebra_type::matrix_type matrix_type;
177
179 typedef typename algebra_type::vector_type vector_type;
180
181 public:
184 : MultiStepTimeDiscretization<TAlgebra>(spDD),
185 m_stage(1), m_scheme("Theta")
186 {
187 set_theta(1.0);
188 this->m_prevSteps = 1;
189 }
190
193 : MultiStepTimeDiscretization<TAlgebra>(spDD),
194 m_stage(1), m_scheme("Theta")
195 {
196 set_theta(theta);
197 this->m_prevSteps = 1;
198 }
199
202 : MultiStepTimeDiscretization<TAlgebra>(spDD),
203 m_stage(1), m_scheme(scheme)
204 {
205 set_theta(1.0);
206 this->m_prevSteps = 1;
207 }
208
209 virtual ~ThetaTimeStep() {};
210
212 void set_scheme(const char* scheme) {m_scheme = scheme;}
213
215 virtual size_t num_stages() const
216 {
217 if (m_scheme == "Theta") return 1;
218 else if (m_scheme == "Alexander") return 2;
219 else if (m_scheme == "FracStep") return 3;
220 else UG_THROW("Step Scheme not recognized.");
221 }
222
224 virtual void set_stage(size_t stage) {m_stage = stage;}
225
227 void set_theta(number theta) {m_theta = theta;}
228
229 protected:
230 virtual number update_scaling(std::vector<number>& vSM,
231 std::vector<number>& vSA,
232 number dt, number currentTime,
234 {
235 // resize scaling factors
236 vSM.resize(2);
237 vSM[0] = 1.;
238 vSM[1] = -1.;
239
240 if(m_scheme == "Theta")
241 {
242 vSA.resize(2);
243 vSA[0] = (m_theta) * dt;
244 vSA[1] = (1.- m_theta) * dt;
245 return currentTime + dt;
246 }
247 else if(m_scheme == "Alexander")
248 {
249 vSA.resize(2);
250 const number gamma = 1 - 1. / sqrt(2.);
251 switch(m_stage)
252 {
253 case 1:
254 vSA[0] = gamma * dt;
255 vSA[1] = 0;
256 return currentTime + gamma * dt;
257 case 2:
258 vSA[0] = gamma * dt;
259 vSA[1] = (1.- 2*gamma) * dt;
260 return currentTime + (1 - gamma) * dt;
261 default:
262 UG_THROW("Alexander scheme has only 2 stages")
263 }
264 }
265 else if(m_scheme == "FracStep")
266 {
267 vSA.resize(2);
268 switch(m_stage)
269 {
270 case 1:
271 vSA[0] = (2.-sqrt(2.)) * (1-1./sqrt(2.)) * dt;
272 vSA[1] = (1.- (2.-sqrt(2.)))* (1-1./sqrt(2.)) * dt;
273 return currentTime + (1-1./sqrt(2.)) * dt;
274 case 2:
275 vSA[0] = (sqrt(2.)-1) * (sqrt(2.)-1) * dt;
276 vSA[1] = (1.- (sqrt(2.)-1)) * (sqrt(2.)-1) * dt;
277 return currentTime + (sqrt(2.)-1) * dt;
278 case 3:
279 vSA[0] = (2.-sqrt(2.)) * (1-1./sqrt(2.)) * dt;
280 vSA[1] = (1.- (2.-sqrt(2.)))* (1-1./sqrt(2.)) * dt;
281 return currentTime + (1-1./sqrt(2.)) * dt;
282 default:
283 UG_THROW("FracStep scheme has only 3 stages")
284 }
285 }
286 else
287 UG_THROW("Unknown Multi-Stage Theta Scheme: "<< m_scheme<<".");
288
289 }
290
292
293 size_t m_stage;
294
295 std::string m_scheme;
296};
297
298
299template <typename TAlgebra>
300class BDF
301 : public MultiStepTimeDiscretization<TAlgebra>
302{
303 public:
306
308 typedef TAlgebra algebra_type;
309
311 typedef typename algebra_type::matrix_type matrix_type;
312
314 typedef typename algebra_type::vector_type vector_type;
315
316 public:
319 : MultiStepTimeDiscretization<TAlgebra>(spDD)
320 {
321 set_order(2);
322 }
323
326 : MultiStepTimeDiscretization<TAlgebra>(spDD)
327 {
328 set_order(order);
329 }
330
331 virtual ~BDF() {};
332
334 void set_order(size_t order) {m_order = order; this->m_prevSteps = order;}
335
337 virtual size_t num_stages() const {return 1;}
338
340 virtual void set_stage(size_t stage)
341 {
342 if(stage!=1) UG_THROW("BDF has only one stage.");
343 }
344
345 protected:
346 virtual number update_scaling(std::vector<number>& vSM,
347 std::vector<number>& vSA,
348 number dt, number currentTime,
350 {
351 // resize scaling factors
352 vSM.resize(m_order+1);
353 vSA.resize(m_order+1);
354
355 // future time
356 const number futureTime = currentTime + dt;
357
358 // get time points
359 if(prevSol->size() < m_order)
360 UG_THROW("BDF("<<m_order<<") needs at least "<< m_order <<
361 " previous solutions, but only "<<prevSol->size()<<"passed.");
362
363 std::vector<number> vTimePoint(m_order+1);
364 vTimePoint[0] = futureTime;
365 for(size_t i = 1; i <= m_order; ++i)
366 vTimePoint[i] = prevSol->time(i-1);
367
368 // evaluate derivative of Lagrange Polynoms at future time
369 for(size_t i = 0; i <= m_order; ++i)
370 {
371 vSM[i] = 0;
372
373 for(size_t j = 0; j < vTimePoint.size(); ++j)
374 {
375 if(j == i) continue;
376
377 number prod = 1;
378 for(size_t k = 0; k < vTimePoint.size(); ++k)
379 {
380 if(k == i) continue;
381 if(k == j) continue;
382
383 prod *= (vTimePoint[0]-vTimePoint[k])/
384 (vTimePoint[i]-vTimePoint[k]);
385 }
386 prod *= 1./(vTimePoint[i]-vTimePoint[j]);
387
388 vSM[i] += prod;
389 }
390 }
391
392 // only first value of vSA != 0
393 const number scale = 1.0 / vSM[0];
394 vSA[0] = scale;
395 for(size_t i = 1; i <= m_order; ++i) vSA[i] = 0;
396
397 // scale first s_m to 1.0
398 for(int i = m_order; i >= 0; --i) vSM[i] *= scale;
399
400 return futureTime;
401 }
402
403 size_t m_order;
404};
405
406
408template <typename TAlgebra>
409class SDIRK
410 : public MultiStepTimeDiscretization<TAlgebra>
411{
412 public:
415
417 typedef TAlgebra algebra_type;
418
420 typedef typename algebra_type::matrix_type matrix_type;
421
423 typedef typename algebra_type::vector_type vector_type;
424
425 public:
428 : MultiStepTimeDiscretization<TAlgebra>(spDD),
429 m_stage(1), m_order(1)
430 {
431 this->m_prevSteps = 1;
432 }
433
436 : MultiStepTimeDiscretization<TAlgebra>(spDD),
437 m_order(order)
438 {
440 this->m_prevSteps = 1;
441 }
442
443 virtual ~SDIRK() {};
444
446 void set_order(int order) {
447 if(order > 4) UG_THROW("DIRK: Only up to order 4.");
448 m_order = order;
449 }
450
452 virtual size_t num_stages() const {
453 switch(m_order){
454 case 1: return 2; // Midpoint
455 case 2: return 2; // Alexander(2)
456 case 3: return 3; // Alexander(3)
457 case 4: return 5; // HairerWanner(4)
458 default: UG_THROW("DIRK: Only up to order 4.")
459 }
460 }
461
463 virtual void set_stage(size_t stage);
464
465 public:
467 number dt);
468 /* Please overwrite any of the following methods, if applicable:
469 virtual void prepare_step_elem(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
470 number dt, const GridLevel& gl);
471 virtual void finish_step(SmartPtr<VectorTimeSeries<vector_type> > currSol);
472 virtual void finish_step_elem(SmartPtr<VectorTimeSeries<vector_type> > currSol,
473 const GridLevel& gl);
474 */
475
476 public:
477 void assemble_jacobian(matrix_type& J, const vector_type& u, const GridLevel& gl);
478
479 void assemble_defect(vector_type& d, const vector_type& u, const GridLevel& gl);
480
481 void assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl);
482
483 void assemble_rhs(vector_type& b, const vector_type& u, const GridLevel& gl);
484
485 void assemble_rhs(vector_type& b, const GridLevel& gl);
486
487 void adjust_solution(vector_type& u, const GridLevel& gl);
488
489 protected:
490 virtual number update_scaling(std::vector<number>& vSM,
491 std::vector<number>& vSA,
492 number dt);
493
494 virtual number update_scaling(std::vector<number>& vSM,
495 std::vector<number>& vSA,
496 number dt, number currentTime,
498 UG_THROW("Not used.");
499 }
500
501 size_t m_stage;
505};
506
507
508} // end namespace ug
509
510
512
513// include implementation
514#include "theta_time_step_impl.h"
515
516#endif /* __H__UG__LIB_DISC__TIME_DISC__THETA_TIME_STEP__ */
Definition smart_pointer.h:296
Definition smart_pointer.h:108
Definition theta_time_step.h:302
BDF(SmartPtr< IDomainDiscretization< TAlgebra > > spDD, size_t order)
theta = 0 -> Backward Euler
Definition theta_time_step.h:325
virtual ~BDF()
Definition theta_time_step.h:331
algebra_type::vector_type vector_type
Type of algebra vector.
Definition theta_time_step.h:314
algebra_type::matrix_type matrix_type
Type of algebra matrix.
Definition theta_time_step.h:311
virtual void set_stage(size_t stage)
sets the stage
Definition theta_time_step.h:340
BDF(SmartPtr< IDomainDiscretization< TAlgebra > > spDD)
constructor
Definition theta_time_step.h:318
TAlgebra algebra_type
Type of algebra.
Definition theta_time_step.h:308
size_t m_order
Definition theta_time_step.h:403
virtual number update_scaling(std::vector< number > &vSM, std::vector< number > &vSA, number dt, number currentTime, ConstSmartPtr< VectorTimeSeries< vector_type > > prevSol)
Error estimator ///.
Definition theta_time_step.h:346
void set_order(size_t order)
sets the theta value
Definition theta_time_step.h:334
virtual size_t num_stages() const
returns the number of stages
Definition theta_time_step.h:337
IDomainDiscretization< TAlgebra > domain_discretization_type
Domain Discretization type.
Definition theta_time_step.h:305
Definition grid_level.h:42
Interface for domain discretization.
Definition domain_disc_interface.h:136
Time Discretization Interface.
Definition time_disc_interface.h:65
SmartPtr< IDomainDiscretization< TAlgebra > > m_spDomDisc
Domain Discretization.
Definition time_disc_interface.h:173
multi step time stepping scheme
Definition theta_time_step.h:57
std::vector< number > m_vScaleStiff
Scaling for stiffness part.
Definition theta_time_step.h:143
algebra_type::vector_type vector_type
Type of algebra vector.
Definition theta_time_step.h:66
virtual size_t num_prev_steps() const
returns number of previous time steps needed
Definition theta_time_step.h:84
virtual number update_scaling(std::vector< number > &vSM, std::vector< number > &vSA, number dt, number currentTime, ConstSmartPtr< VectorTimeSeries< vector_type > > prevSol)=0
Error estimator ///.
void assemble_jacobian(matrix_type &J, const vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:113
bool is_error_valid()
returns whether error indicators are valid
Definition theta_time_step.h:129
TAlgebra algebra_type
Type of algebra.
Definition theta_time_step.h:60
virtual void finish_step(SmartPtr< VectorTimeSeries< vector_type > > currSol)
Definition theta_time_step_impl.h:295
number m_dt
Time Step size.
Definition theta_time_step.h:146
virtual void finish_step_elem(SmartPtr< VectorTimeSeries< vector_type > > currSol, const GridLevel &gl)
Definition theta_time_step_impl.h:313
virtual void prepare_step_elem(SmartPtr< VectorTimeSeries< vector_type > > prevSol, number dt, const GridLevel &gl)
prepares the assembling of Defect/Jacobian for a time step
Definition theta_time_step_impl.h:77
virtual void prepare_step(SmartPtr< VectorTimeSeries< vector_type > > prevSol, number dt)
prepares the assembling of Defect/Jacobian for a time step
Definition theta_time_step_impl.h:46
void adjust_solution(vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:173
SmartPtr< VectorTimeSeries< vector_type > > m_pPrevSol
Previous solutions.
Definition theta_time_step.h:145
CPUAlgebra::vector_type error_vector_type
Type of algebra vector.
Definition theta_time_step.h:69
IDomainDiscretization< algebra_type > domain_discretization_type
Domain Discretization type.
Definition theta_time_step.h:72
void assemble_defect(vector_type &d, const vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:143
void calc_error(const vector_type &u, error_vector_type &u_vtk)
Definition theta_time_step.h:122
void invalidate_error()
Definition theta_time_step.h:126
virtual ~MultiStepTimeDiscretization()
Definition theta_time_step.h:81
void calc_error(const vector_type &u, error_vector_type *u_vtk)
Error estimator ///.
Definition theta_time_step_impl.h:261
void calc_error(const vector_type &u)
Definition theta_time_step.h:121
algebra_type::matrix_type matrix_type
Type of algebra matrix.
Definition theta_time_step.h:63
number m_futureTime
Future Time.
Definition theta_time_step.h:147
virtual number future_time() const
returns the future time point (i.e. the one that will be computed)
Definition theta_time_step.h:101
MultiStepTimeDiscretization(SmartPtr< IDomainDiscretization< algebra_type > > spDD)
constructor
Definition theta_time_step.h:76
std::vector< number > m_vScaleMass
Scaling for mass part.
Definition theta_time_step.h:142
void assemble_linear(matrix_type &A, vector_type &b, const GridLevel &gl)
Definition theta_time_step_impl.h:184
size_t m_prevSteps
number of previous steps needed.
Definition theta_time_step.h:141
void assemble_rhs(vector_type &b, const vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:231
Singly Diagonal Implicit Runge Kutta Method.
Definition theta_time_step.h:411
number m_lastTime
Definition theta_time_step.h:504
TAlgebra algebra_type
Type of algebra.
Definition theta_time_step.h:417
number m_Time0
Definition theta_time_step.h:503
void assemble_jacobian(matrix_type &J, const vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:643
algebra_type::vector_type vector_type
Type of algebra vector.
Definition theta_time_step.h:423
SDIRK(SmartPtr< IDomainDiscretization< TAlgebra > > spDD)
default constructor (implicit Euler)
Definition theta_time_step.h:427
void assemble_linear(matrix_type &A, vector_type &b, const GridLevel &gl)
Definition theta_time_step_impl.h:708
void assemble_rhs(vector_type &b, const vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:722
virtual void prepare_step(SmartPtr< VectorTimeSeries< vector_type > > prevSol, number dt)
prepares the assembling of Defect/Jacobian for a time step
Definition theta_time_step_impl.h:612
virtual void set_stage(size_t stage)
sets the stage
Definition theta_time_step_impl.h:339
SDIRK(SmartPtr< IDomainDiscretization< TAlgebra > > spDD, int order)
theta = 1.0 -> Implicit Euler, 0.0 -> Explicit Euler
Definition theta_time_step.h:435
int m_order
Definition theta_time_step.h:502
IDomainDiscretization< TAlgebra > domain_discretization_type
Domain Discretization type.
Definition theta_time_step.h:414
void assemble_defect(vector_type &d, const vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:670
virtual ~SDIRK()
Definition theta_time_step.h:443
void adjust_solution(vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:697
size_t m_stage
Definition theta_time_step.h:501
virtual number update_scaling(std::vector< number > &vSM, std::vector< number > &vSA, number dt)
Definition theta_time_step_impl.h:345
virtual number update_scaling(std::vector< number > &vSM, std::vector< number > &vSA, number dt, number currentTime, ConstSmartPtr< VectorTimeSeries< vector_type > > prevSol)
Error estimator ///.
Definition theta_time_step.h:494
void set_order(int order)
sets the scheme
Definition theta_time_step.h:446
virtual size_t num_stages() const
returns number of stages
Definition theta_time_step.h:452
algebra_type::matrix_type matrix_type
Type of algebra matrix.
Definition theta_time_step.h:420
theta time stepping scheme
Definition theta_time_step.h:167
TAlgebra algebra_type
Type of algebra.
Definition theta_time_step.h:173
virtual void set_stage(size_t stage)
sets the stage
Definition theta_time_step.h:224
IDomainDiscretization< TAlgebra > domain_discretization_type
Domain Discretization type.
Definition theta_time_step.h:170
virtual size_t num_stages() const
returns number of stages
Definition theta_time_step.h:215
algebra_type::matrix_type matrix_type
Type of algebra matrix.
Definition theta_time_step.h:176
algebra_type::vector_type vector_type
Type of algebra vector.
Definition theta_time_step.h:179
number m_theta
Definition theta_time_step.h:291
void set_scheme(const char *scheme)
sets the scheme
Definition theta_time_step.h:212
virtual ~ThetaTimeStep()
Definition theta_time_step.h:209
ThetaTimeStep(SmartPtr< IDomainDiscretization< TAlgebra > > spDD, number theta)
theta = 1.0 -> Implicit Euler, 0.0 -> Explicit Euler
Definition theta_time_step.h:192
virtual number update_scaling(std::vector< number > &vSM, std::vector< number > &vSA, number dt, number currentTime, ConstSmartPtr< VectorTimeSeries< vector_type > > prevSol)
Error estimator ///.
Definition theta_time_step.h:230
ThetaTimeStep(SmartPtr< IDomainDiscretization< TAlgebra > > spDD)
default constructor (implicit Euler)
Definition theta_time_step.h:183
ThetaTimeStep(SmartPtr< IDomainDiscretization< TAlgebra > > spDD, const char *scheme)
theta = 1.0 -> Implicit Euler, 0.0 -> Explicit Euler
Definition theta_time_step.h:201
void set_theta(number theta)
sets the theta value
Definition theta_time_step.h:227
std::string m_scheme
Definition theta_time_step.h:295
size_t m_stage
Definition theta_time_step.h:293
time series of solutions and corresponding time point
Definition solution_time_series.h:59
#define UG_THROW(msg)
Definition error.h:57
double number
Definition types.h:124
the ug namespace