Line data Source code
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
41 : #include "lib_algebra/cpu_algebra_types.h"
42 : #include "common/common.h"
43 :
44 : // module-intern libraries
45 : #include "lib_disc/time_disc/time_disc_interface.h"
46 : #include "lib_disc/local_finite_element/common/lagrange1d.h"
47 :
48 : namespace ug{
49 :
50 : /// \ingroup lib_disc_time_assemble
51 : /// @{
52 :
53 : /// multi step time stepping scheme
54 : template <typename TAlgebra>
55 : class MultiStepTimeDiscretization
56 : : public ITimeDiscretization<TAlgebra>
57 : {
58 : public:
59 : /// Type of algebra
60 : typedef TAlgebra algebra_type;
61 :
62 : /// Type of algebra matrix
63 : typedef typename algebra_type::matrix_type matrix_type;
64 :
65 : /// Type of algebra vector
66 : typedef typename algebra_type::vector_type vector_type;
67 :
68 : /// Type of algebra vector
69 : typedef typename CPUAlgebra::vector_type error_vector_type;
70 :
71 : /// Domain Discretization type
72 : typedef IDomainDiscretization<algebra_type> domain_discretization_type;
73 :
74 : public:
75 : /// constructor
76 0 : MultiStepTimeDiscretization(SmartPtr<IDomainDiscretization<algebra_type> > spDD)
77 : : ITimeDiscretization<TAlgebra>(spDD),
78 0 : m_pPrevSol(NULL)
79 0 : {}
80 :
81 0 : virtual ~MultiStepTimeDiscretization(){};
82 :
83 : /// \copydoc ITimeDiscretization::num_prev_steps()
84 0 : virtual size_t num_prev_steps() const {return m_prevSteps;}
85 :
86 : /// \copydoc ITimeDiscretization::prepare_step()
87 : virtual void prepare_step(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
88 : number dt);
89 :
90 : /// \copydoc ITimeDiscretization::prepare_step_elem()
91 : virtual void prepare_step_elem(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
92 : number dt, const GridLevel& gl);
93 :
94 : /// \copydoc ITimeDiscretization::finish_step()
95 : virtual void finish_step(SmartPtr<VectorTimeSeries<vector_type> > currSol);
96 :
97 : /// \copydoc ITimeDiscretization::finish_step_elem()
98 : virtual void finish_step_elem(SmartPtr<VectorTimeSeries<vector_type> > currSol,
99 : const GridLevel& gl);
100 :
101 0 : 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 :
116 : ///////////////////////////////////////////////////////////////////
117 : /// Error estimator ///
118 :
119 : /// calculates error indicators for elements from error estimators
120 : void calc_error(const vector_type& u, error_vector_type* u_vtk);
121 0 : void calc_error(const vector_type& u) {calc_error(u, NULL);};
122 0 : void calc_error(const vector_type& u, error_vector_type& u_vtk){calc_error(u, &u_vtk);};
123 :
124 : /// marks error indicators as invalid; in order to revalidate them,
125 : /// they will have to be newly calculated by a call to calc_error
126 0 : void invalidate_error() {this->m_spDomDisc->invalidate_error();};
127 :
128 : /// returns whether error indicators are valid
129 0 : bool is_error_valid() {return this->m_spDomDisc->is_error_valid();};
130 :
131 : /// Error estimator ///
132 : ///////////////////////////////////////////////////////////////////
133 :
134 : protected:
135 : /// updates the scaling factors, returns the future time
136 : virtual number update_scaling(std::vector<number>& vSM,
137 : std::vector<number>& vSA,
138 : number dt, number currentTime,
139 : ConstSmartPtr<VectorTimeSeries<vector_type> > prevSol) = 0;
140 :
141 : size_t m_prevSteps; ///< number of previous steps needed.
142 : std::vector<number> m_vScaleMass; ///< Scaling for mass part
143 : std::vector<number> m_vScaleStiff; ///< Scaling for stiffness part
144 :
145 : SmartPtr<VectorTimeSeries<vector_type> > m_pPrevSol; ///< Previous solutions
146 : number m_dt; ///< Time Step size
147 : number m_futureTime; ///< Future Time
148 : };
149 :
150 : /// theta time stepping scheme
151 : /**
152 : * This time stepping scheme discretizes equations of the form
153 : * \f[
154 : * \partial_t u(t) = f(t)
155 : * \f]
156 : * as
157 : * \f[
158 : * \frac{u(t^{k+1}) - u(t^k)}{\Delta t} = (1-\theta) \cdot f(t^{k+1})
159 : * + \theta \cdot f(t^k)
160 : * \f]
161 : *
162 : * Thus, for \f$\theta = 1 \f$ this is the Backward-Euler time stepping.
163 : */
164 : template <typename TAlgebra>
165 : class ThetaTimeStep
166 : : public MultiStepTimeDiscretization<TAlgebra>
167 : {
168 : public:
169 : /// Domain Discretization type
170 : typedef IDomainDiscretization<TAlgebra> domain_discretization_type;
171 :
172 : /// Type of algebra
173 : typedef TAlgebra algebra_type;
174 :
175 : /// Type of algebra matrix
176 : typedef typename algebra_type::matrix_type matrix_type;
177 :
178 : /// Type of algebra vector
179 : typedef typename algebra_type::vector_type vector_type;
180 :
181 : public:
182 : /// default constructor (implicit Euler)
183 0 : ThetaTimeStep(SmartPtr<IDomainDiscretization<TAlgebra> > spDD)
184 : : MultiStepTimeDiscretization<TAlgebra>(spDD),
185 0 : m_stage(1), m_scheme("Theta")
186 : {
187 : set_theta(1.0);
188 0 : this->m_prevSteps = 1;
189 0 : }
190 :
191 : /// theta = 1.0 -> Implicit Euler, 0.0 -> Explicit Euler
192 0 : ThetaTimeStep(SmartPtr<IDomainDiscretization<TAlgebra> > spDD, number theta)
193 : : MultiStepTimeDiscretization<TAlgebra>(spDD),
194 0 : m_stage(1), m_scheme("Theta")
195 : {
196 : set_theta(theta);
197 0 : this->m_prevSteps = 1;
198 0 : }
199 :
200 : /// theta = 1.0 -> Implicit Euler, 0.0 -> Explicit Euler
201 0 : ThetaTimeStep(SmartPtr<IDomainDiscretization<TAlgebra> > spDD, const char* scheme)
202 : : MultiStepTimeDiscretization<TAlgebra>(spDD),
203 0 : m_stage(1), m_scheme(scheme)
204 : {
205 : set_theta(1.0);
206 0 : this->m_prevSteps = 1;
207 0 : }
208 :
209 0 : virtual ~ThetaTimeStep() {};
210 :
211 : /// sets the scheme
212 0 : void set_scheme(const char* scheme) {m_scheme = scheme;}
213 :
214 : /// returns number of stages
215 0 : virtual size_t num_stages() const
216 : {
217 0 : if (m_scheme == "Theta") return 1;
218 0 : else if (m_scheme == "Alexander") return 2;
219 0 : else if (m_scheme == "FracStep") return 3;
220 0 : else UG_THROW("Step Scheme not recognized.");
221 : }
222 :
223 : /// sets the stage
224 0 : virtual void set_stage(size_t stage) {m_stage = stage;}
225 :
226 : /// sets the theta value
227 0 : void set_theta(number theta) {m_theta = theta;}
228 :
229 : protected:
230 0 : virtual number update_scaling(std::vector<number>& vSM,
231 : std::vector<number>& vSA,
232 : number dt, number currentTime,
233 : ConstSmartPtr<VectorTimeSeries<vector_type> > prevSol)
234 : {
235 : // resize scaling factors
236 0 : vSM.resize(2);
237 0 : vSM[0] = 1.;
238 0 : vSM[1] = -1.;
239 :
240 0 : if(m_scheme == "Theta")
241 : {
242 0 : vSA.resize(2);
243 0 : vSA[0] = (m_theta) * dt;
244 0 : vSA[1] = (1.- m_theta) * dt;
245 0 : return currentTime + dt;
246 : }
247 0 : else if(m_scheme == "Alexander")
248 : {
249 0 : vSA.resize(2);
250 : const number gamma = 1 - 1. / sqrt(2.);
251 0 : switch(m_stage)
252 : {
253 : case 1:
254 0 : vSA[0] = gamma * dt;
255 0 : vSA[1] = 0;
256 0 : return currentTime + gamma * dt;
257 : case 2:
258 0 : vSA[0] = gamma * dt;
259 0 : vSA[1] = (1.- 2*gamma) * dt;
260 0 : return currentTime + (1 - gamma) * dt;
261 0 : default:
262 0 : UG_THROW("Alexander scheme has only 2 stages")
263 : }
264 : }
265 0 : else if(m_scheme == "FracStep")
266 : {
267 0 : vSA.resize(2);
268 0 : switch(m_stage)
269 : {
270 : case 1:
271 0 : vSA[0] = (2.-sqrt(2.)) * (1-1./sqrt(2.)) * dt;
272 0 : vSA[1] = (1.- (2.-sqrt(2.)))* (1-1./sqrt(2.)) * dt;
273 0 : return currentTime + (1-1./sqrt(2.)) * dt;
274 : case 2:
275 0 : vSA[0] = (sqrt(2.)-1) * (sqrt(2.)-1) * dt;
276 0 : vSA[1] = (1.- (sqrt(2.)-1)) * (sqrt(2.)-1) * dt;
277 0 : return currentTime + (sqrt(2.)-1) * dt;
278 : case 3:
279 0 : vSA[0] = (2.-sqrt(2.)) * (1-1./sqrt(2.)) * dt;
280 0 : vSA[1] = (1.- (2.-sqrt(2.)))* (1-1./sqrt(2.)) * dt;
281 0 : return currentTime + (1-1./sqrt(2.)) * dt;
282 0 : default:
283 0 : UG_THROW("FracStep scheme has only 3 stages")
284 : }
285 : }
286 : else
287 0 : UG_THROW("Unknown Multi-Stage Theta Scheme: "<< m_scheme<<".");
288 :
289 : }
290 :
291 : number m_theta;
292 :
293 : size_t m_stage;
294 :
295 : std::string m_scheme;
296 : };
297 :
298 :
299 : template <typename TAlgebra>
300 : class BDF
301 : : public MultiStepTimeDiscretization<TAlgebra>
302 : {
303 : public:
304 : /// Domain Discretization type
305 : typedef IDomainDiscretization<TAlgebra> domain_discretization_type;
306 :
307 : /// Type of algebra
308 : typedef TAlgebra algebra_type;
309 :
310 : /// Type of algebra matrix
311 : typedef typename algebra_type::matrix_type matrix_type;
312 :
313 : /// Type of algebra vector
314 : typedef typename algebra_type::vector_type vector_type;
315 :
316 : public:
317 : /// constructor
318 0 : BDF(SmartPtr<IDomainDiscretization<TAlgebra> > spDD)
319 0 : : MultiStepTimeDiscretization<TAlgebra>(spDD)
320 : {
321 : set_order(2);
322 0 : }
323 :
324 : /// theta = 0 -> Backward Euler
325 0 : BDF(SmartPtr<IDomainDiscretization<TAlgebra> > spDD, size_t order)
326 0 : : MultiStepTimeDiscretization<TAlgebra>(spDD)
327 : {
328 : set_order(order);
329 0 : }
330 :
331 0 : virtual ~BDF() {};
332 :
333 : /// sets the theta value
334 0 : void set_order(size_t order) {m_order = order; this->m_prevSteps = order;}
335 :
336 : /// returns the number of stages
337 0 : virtual size_t num_stages() const {return 1;}
338 :
339 : /// sets the stage
340 0 : virtual void set_stage(size_t stage)
341 : {
342 0 : if(stage!=1) UG_THROW("BDF has only one stage.");
343 0 : }
344 :
345 : protected:
346 0 : virtual number update_scaling(std::vector<number>& vSM,
347 : std::vector<number>& vSA,
348 : number dt, number currentTime,
349 : ConstSmartPtr<VectorTimeSeries<vector_type> > prevSol)
350 : {
351 : // resize scaling factors
352 0 : vSM.resize(m_order+1);
353 0 : vSA.resize(m_order+1);
354 :
355 : // future time
356 0 : const number futureTime = currentTime + dt;
357 :
358 : // get time points
359 0 : if(prevSol->size() < m_order)
360 0 : UG_THROW("BDF("<<m_order<<") needs at least "<< m_order <<
361 : " previous solutions, but only "<<prevSol->size()<<"passed.");
362 :
363 0 : std::vector<number> vTimePoint(m_order+1);
364 0 : vTimePoint[0] = futureTime;
365 0 : for(size_t i = 1; i <= m_order; ++i)
366 0 : vTimePoint[i] = prevSol->time(i-1);
367 :
368 : // evaluate derivative of Lagrange Polynoms at future time
369 0 : for(size_t i = 0; i <= m_order; ++i)
370 : {
371 0 : vSM[i] = 0;
372 :
373 0 : for(size_t j = 0; j < vTimePoint.size(); ++j)
374 : {
375 0 : if(j == i) continue;
376 :
377 : number prod = 1;
378 0 : for(size_t k = 0; k < vTimePoint.size(); ++k)
379 : {
380 0 : if(k == i) continue;
381 0 : if(k == j) continue;
382 :
383 0 : prod *= (vTimePoint[0]-vTimePoint[k])/
384 0 : (vTimePoint[i]-vTimePoint[k]);
385 : }
386 0 : prod *= 1./(vTimePoint[i]-vTimePoint[j]);
387 :
388 0 : vSM[i] += prod;
389 : }
390 : }
391 :
392 : // only first value of vSA != 0
393 0 : const number scale = 1.0 / vSM[0];
394 0 : vSA[0] = scale;
395 0 : for(size_t i = 1; i <= m_order; ++i) vSA[i] = 0;
396 :
397 : // scale first s_m to 1.0
398 0 : for(int i = m_order; i >= 0; --i) vSM[i] *= scale;
399 :
400 0 : return futureTime;
401 0 : }
402 :
403 : size_t m_order;
404 : };
405 :
406 :
407 : /// Singly Diagonal Implicit Runge Kutta Method
408 : template <typename TAlgebra>
409 : class SDIRK
410 : : public MultiStepTimeDiscretization<TAlgebra>
411 : {
412 : public:
413 : /// Domain Discretization type
414 : typedef IDomainDiscretization<TAlgebra> domain_discretization_type;
415 :
416 : /// Type of algebra
417 : typedef TAlgebra algebra_type;
418 :
419 : /// Type of algebra matrix
420 : typedef typename algebra_type::matrix_type matrix_type;
421 :
422 : /// Type of algebra vector
423 : typedef typename algebra_type::vector_type vector_type;
424 :
425 : public:
426 : /// default constructor (implicit Euler)
427 0 : SDIRK(SmartPtr<IDomainDiscretization<TAlgebra> > spDD)
428 : : MultiStepTimeDiscretization<TAlgebra>(spDD),
429 0 : m_stage(1), m_order(1)
430 : {
431 0 : this->m_prevSteps = 1;
432 0 : }
433 :
434 : /// theta = 1.0 -> Implicit Euler, 0.0 -> Explicit Euler
435 0 : SDIRK(SmartPtr<IDomainDiscretization<TAlgebra> > spDD, int order)
436 : : MultiStepTimeDiscretization<TAlgebra>(spDD),
437 0 : m_order(order)
438 : {
439 0 : set_order(m_order);
440 0 : this->m_prevSteps = 1;
441 0 : }
442 :
443 0 : virtual ~SDIRK() {};
444 :
445 : /// sets the scheme
446 0 : void set_order(int order) {
447 0 : if(order > 4) UG_THROW("DIRK: Only up to order 4.");
448 0 : m_order = order;
449 0 : }
450 :
451 : /// returns number of stages
452 0 : virtual size_t num_stages() const {
453 0 : 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 0 : default: UG_THROW("DIRK: Only up to order 4.")
459 : }
460 : }
461 :
462 : /// sets the stage
463 : virtual void set_stage(size_t stage);
464 :
465 : public:
466 : virtual void prepare_step(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
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 0 : virtual number update_scaling(std::vector<number>& vSM,
495 : std::vector<number>& vSA,
496 : number dt, number currentTime,
497 : ConstSmartPtr<VectorTimeSeries<vector_type> > prevSol) {
498 0 : UG_THROW("Not used.");
499 : }
500 :
501 : size_t m_stage;
502 : int m_order;
503 : number m_Time0;
504 : number m_lastTime;
505 : };
506 :
507 :
508 : } // end namespace ug
509 :
510 :
511 : /// }
512 :
513 : // include implementation
514 : #include "theta_time_step_impl.h"
515 :
516 : #endif /* __H__UG__LIB_DISC__TIME_DISC__THETA_TIME_STEP__ */
|
