Loading [MathJax]/extensions/tex2jax.js
ug4
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
theta_time_step_impl.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_IMPL__
34#define __H__UG__LIB_DISC__TIME_DISC__THETA_TIME_STEP_IMPL__
35
36#include "theta_time_step.h"
37
38#ifndef M_PI
39#define M_PI 3.14159265358979323846264338327950288 /* pi */
40#endif
41
42namespace ug{
43
44template <typename TAlgebra>
45void MultiStepTimeDiscretization<TAlgebra>::
46prepare_step(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
47 number dt)
48{
49 PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_prepare_step, "discretization MultiStepTimeDiscretization");
50// perform checks
51 if(prevSol->size() < m_prevSteps)
52 UG_THROW("MultiStepTimeDiscretization::prepare_step:"
53 " Number of previous solutions must be at least "<<
54 m_prevSteps <<", but only "<< prevSol->size() << " passed.\n");
55
56// remember old values
57 m_pPrevSol = prevSol;
58
59// remember time step size
60 m_dt = dt;
61
62// update scalings
63 m_futureTime = update_scaling(m_vScaleMass, m_vScaleStiff,
64 m_dt, m_pPrevSol->time(0),
65 m_pPrevSol);
66
67// prepare time step (elemDisc-wise)
68 try
69 {
70 this->m_spDomDisc->prepare_timestep(m_pPrevSol, m_futureTime);
71 }
72 UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot prepare time step.");
73}
74
75template <typename TAlgebra>
76void MultiStepTimeDiscretization<TAlgebra>::
77prepare_step_elem(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
78 number dt, const GridLevel& gl)
79{
80 PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_step_elem, "discretization MultiStepTimeDiscretization");
81// perform checks
82 if(prevSol->size() < m_prevSteps)
83 UG_THROW("MultiStepTimeDiscretization::prepare_step_elem:"
84 " Number of previous solutions must be at least "<<
85 m_prevSteps <<", but only "<< prevSol->size() << " passed.\n");
86
87// remember old values
88 m_pPrevSol = prevSol;
89
90// remember time step size
91 m_dt = dt;
92
93// update scalings
94 m_futureTime = update_scaling(m_vScaleMass, m_vScaleStiff,
95 m_dt, m_pPrevSol->time(0),
96 m_pPrevSol);
97
98// prepare time step (elemDisc-wise)
99 try
100 {
101 this->m_spDomDisc->prepare_timestep(m_pPrevSol, m_futureTime, gl);
102 }
103 UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot prepare time step.");
104
105// prepare timestep element-wise
106 try{
107 this->m_spDomDisc->prepare_timestep_elem(m_pPrevSol, gl);
108 }UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot prepare timestep element-wise.");
109}
110
111template <typename TAlgebra>
112void MultiStepTimeDiscretization<TAlgebra>::
113assemble_jacobian(matrix_type& J, const vector_type& u, const GridLevel& gl)
114{
115 PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_assemble_jacobian, "discretization MultiStepTimeDiscretization");
116// perform checks
117 if(m_pPrevSol->size() < m_prevSteps)
118 UG_THROW("MultiStepTimeDiscretization::assemble_jacobian:"
119 " Number of previous solutions must be at least "<<
120 m_prevSteps <<", but only "<< m_pPrevSol->size() << " passed.");
121
122// push unknown solution to solution time series
123// ATTENTION: Here, we must cast away the constness of the solution, but note,
124// that we pass pPrevSol as a const object in assemble_... Thus,
125// the solution will not be changed there and we pop it from the
126// Solution list afterwards, such that nothing happens to u
127 // \todo: avoid this hack, use smart ptr properly
128 int DummyRefCount = 2;
129 SmartPtr<vector_type> pU(const_cast<vector_type*>(&u), &DummyRefCount);
130 m_pPrevSol->push(pU, m_futureTime);
131
132// assemble jacobian using current iterate
133 try{
134 this->m_spDomDisc->assemble_jacobian(J, m_pPrevSol, m_vScaleStiff[0], gl);
135 }UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot assemble jacobian.");
136
137// pop unknown solution to solution time series
138 m_pPrevSol->remove_latest();
139}
140
141template <typename TAlgebra>
142void MultiStepTimeDiscretization<TAlgebra>::
143assemble_defect(vector_type& d, const vector_type& u, const GridLevel& gl)
144{
145 PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_assemble_defect, "discretization MultiStepTimeDiscretization");
146// perform checks
147 if(m_pPrevSol->size() < m_prevSteps)
148 UG_THROW("MultiStepTimeDiscretization::assemble_defect:"
149 " Number of previous solutions must be at least "<<
150 m_prevSteps <<", but only "<< m_pPrevSol->size() << " passed.");
151
152// push unknown solution to solution time series
153// ATTENTION: Here, we must cast away the constness of the solution, but note,
154// that we pass pPrevSol as a const object in assemble_... Thus,
155// the solution will not be changed there and we pop it from the
156// Solution list afterwards, such that nothing happens to u
157 // \todo: avoid this hack, use smart ptr properly
158 int DummyRefCount = 2;
159 SmartPtr<vector_type> pU(const_cast<vector_type*>(&u), &DummyRefCount);
160 m_pPrevSol->push(pU, m_futureTime);
161
162// future solution part
163 try{
164 this->m_spDomDisc->assemble_defect(d, m_pPrevSol, m_vScaleMass, m_vScaleStiff, gl);
165 }UG_CATCH_THROW("ThetaTimeStep: Cannot assemble defect.");
166
167// pop unknown solution to solution time series
168 m_pPrevSol->remove_latest();
169}
170
171template <typename TAlgebra>
172void MultiStepTimeDiscretization<TAlgebra>::
173adjust_solution(vector_type& u, const GridLevel& gl)
174{
175 PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_adjust_solution, "discretization MultiStepTimeDiscretization");
176// adjust solution
177 try{
178 this->m_spDomDisc->adjust_solution(u, m_futureTime, gl);
179 }UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot adjust solution.");
180}
181
182template <typename TAlgebra>
183void MultiStepTimeDiscretization<TAlgebra>::
184assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl)
185{
186 PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_assemble_linear, "discretization MultiStepTimeDiscretization");
187// perform checks
188 if(m_pPrevSol->size() < m_prevSteps)
189 UG_THROW("MultiStepTimeDiscretization::assemble_linear:"
190 " Number of previous solutions must be at least "<<
191 m_prevSteps <<", but only "<< m_pPrevSol->size() << " passed.");
192
193
194// push unknown solution to solution time series (not used, but formally needed)
195 m_pPrevSol->push(m_pPrevSol->latest(), m_futureTime);
196
197// assemble jacobian using current iterate
198 try{
199 this->m_spDomDisc->assemble_linear(A, b, m_pPrevSol, m_vScaleMass, m_vScaleStiff, gl);
200 }UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot assemble jacobian.");
201
202// pop unknown solution from solution time series
203 m_pPrevSol->remove_latest();
204}
205
206template <typename TAlgebra>
207void MultiStepTimeDiscretization<TAlgebra>::
208assemble_rhs(vector_type& b, const GridLevel& gl)
209{
210 PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_assemble_rhs, "discretization MultiStepTimeDiscretization");
211// perform checks
212 if(m_pPrevSol->size() < m_prevSteps)
213 UG_THROW("MultiStepTimeDiscretization::assemble_linear:"
214 " Number of previous solutions must be at least "<<
215 m_prevSteps <<", but only "<< m_pPrevSol->size() << " passed.");
216
217// push unknown solution to solution time series (not used, but formally needed)
218 m_pPrevSol->push(m_pPrevSol->latest(), m_futureTime);
219
220// assemble jacobian using current iterate
221 try{
222 this->m_spDomDisc->assemble_rhs(b, m_pPrevSol, m_vScaleMass, m_vScaleStiff, gl);
223 }UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot assemble jacobian.");
224
225// pop unknown solution from solution time series
226 m_pPrevSol->remove_latest();
227}
228
229template <typename TAlgebra>
230void MultiStepTimeDiscretization<TAlgebra>::
231assemble_rhs(vector_type& b, const vector_type& u, const GridLevel& gl)
232{
233 PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_assemble_rhs, "discretization MultiStepTimeDiscretization");
234// perform checks
235 if(m_pPrevSol->size() < m_prevSteps)
236 UG_THROW("MultiStepTimeDiscretization::assemble_linear:"
237 " Number of previous solutions must be at least "<<
238 m_prevSteps <<", but only "<< m_pPrevSol->size() << " passed.");
239
240// push unknown solution to solution time series
241// ATTENTION: Here, we must cast away the constness of the solution, but note,
242// that we pass pPrevSol as a const object in assemble_... Thus,
243// the solution will not be changed there and we pop it from the
244// Solution list afterwards, such that nothing happens to u
245 // \todo: avoid this hack, use smart ptr properly
246 int DummyRefCount = 2;
247 SmartPtr<vector_type> pU(const_cast<vector_type*>(&u), &DummyRefCount);
248 m_pPrevSol->push(pU, m_futureTime);
249
250// assemble jacobian using current iterate
251 try{
252 this->m_spDomDisc->assemble_rhs(b, m_pPrevSol, m_vScaleMass, m_vScaleStiff, gl);
253 }UG_CATCH_THROW("ThetaTimeStep: Cannot assemble jacobian.");
254
255// pop unknown solution to solution time series
256 m_pPrevSol->remove_latest();
257}
258
259template<typename TAlgebra>
260void MultiStepTimeDiscretization<TAlgebra>::
261calc_error(const vector_type& u, error_vector_type* u_vtk)
262{
263 PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_calc_error, "discretization MultiStepTimeDiscretization");
264// perform checks
265 if(m_pPrevSol->size() < m_prevSteps)
266 UG_THROW("MultiStepTimeDiscretization::calc_error:"
267 " Number of previous solutions must be at least "<<
268 m_prevSteps <<", but only "<< m_pPrevSol->size() << " passed.");
269
270// push unknown solution to solution time series
271// ATTENTION: Here, we must cast away the constness of the solution, but note,
272// that we pass pPrevSol as a const object in assemble_... Thus,
273// the solution will not be changed there and we pop it from the
274// Solution list afterwards, such that nothing happens to u
275 // \todo: avoid this hack, use smart ptr properly
276 int DummyRefCount = 2;
277 SmartPtr<vector_type> pU(const_cast<vector_type*>(&u), &DummyRefCount);
278 m_pPrevSol->push(pU, m_futureTime);
279
280// assemble error estimators using current iterate
281 try
282 {
283 GridLevel gl = GridLevel(); // new TOP-SURFACE grid level
284 this->m_spDomDisc->calc_error(m_pPrevSol, m_vScaleMass, m_vScaleStiff, gl, u_vtk);
285 }
286 UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot assemble error estimators.");
287
288// pop unknown solution to solution time series
289 m_pPrevSol->remove_latest();
290}
291
292
293template <typename TAlgebra>
294void MultiStepTimeDiscretization<TAlgebra>::
295finish_step(SmartPtr<VectorTimeSeries<vector_type> > currSol)
296{
297// perform checks whether 'currSol' is a solutionTimeSeries only with the new values
298 if (currSol->time(0) != m_futureTime)
299 UG_THROW("MultiStepTimeDiscretization::finish_step:"
300 " The solution of the SolutionTimeSeries used in this function"
301 " does not coincide with the current solution! ");
302
303// finish timestep using the current solution
304 try
305 {
306 this->m_spDomDisc->finish_timestep(currSol);
307 }
308 UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot finish timestep.");
309}
310
311template <typename TAlgebra>
312void MultiStepTimeDiscretization<TAlgebra>::
313finish_step_elem(SmartPtr<VectorTimeSeries<vector_type> > currSol,
314 const GridLevel& gl)
315{
316// perform checks whether 'currSol' is a solutionTimeSeries only with the new values
317 if(currSol->time(0) != m_futureTime)
318 UG_THROW("MultiStepTimeDiscretization::finish_step_elem:"
319 " The solution of the SolutionTimeSeries used in this function"
320 " does not coincide with the current solution! ");
321
322// finish timestep using the current solution
323 try
324 {
325 this->m_spDomDisc->finish_timestep(currSol, gl);
326 }
327 UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot finish timestep.");
328
329// finish timestep element-wise using the current solution
330 try{
331 this->m_spDomDisc->finish_timestep_elem(currSol, gl);
332 }UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot finish timestep element-wise.");
333}
334
337
338template <typename TAlgebra>
339void SDIRK<TAlgebra>::set_stage(size_t stage)
340{
341 m_stage = stage;
342}
343
344template <typename TAlgebra>
345number SDIRK<TAlgebra>::update_scaling(std::vector<number>& vSM,
346 std::vector<number>& vSA,
347 number dt)
348{
349 if(m_order == 1) // Mittelpunkt
350 {
351 switch(m_stage)
352 {
353 case 1:
354 vSM.resize(2);
355 vSA.resize(2);
356 vSM[0] = 1.;
357 vSM[1] = -1.;
358 vSA[0] = dt * 1./2.;
359 vSA[1] = 0;
360 return m_Time0 + 1./2. * dt;
361 case 2:
362 vSM.resize(3);
363 vSA.resize(3);
364 vSM[0] = 1.;
365 vSM[1] = 0;
366 vSM[2] = -1.;
367 vSA[0] = 0;
368 vSA[1] = dt;
369 vSA[2] = 0;
370 return m_Time0 + dt;
371 default:
372 UG_THROW("Midpoint scheme has only 2 stages")
373 }
374 }
375 else if(m_order == 2) // Alexander, order 2
376 {
377 const number gamma = 1 - 1. / sqrt(2.);
378 switch(m_stage)
379 {
380 case 1:
381 vSM.resize(2);
382 vSA.resize(2);
383 vSM[0] = 1.;
384 vSM[1] = -1.;
385 vSA[0] = dt * gamma;
386 vSA[1] = 0;
387 return m_Time0 + gamma * dt;
388 case 2:
389 vSM.resize(3);
390 vSA.resize(3);
391 vSM[0] = 1.;
392 vSM[1] = 0;
393 vSM[2] = -1;
394 vSA[0] = dt * gamma;
395 vSA[1] = dt * (1. - gamma);
396 vSA[2] = 0;
397 return m_Time0 + dt;
398 default:
399 UG_THROW("Alexander(2) scheme has only 2 stages")
400 }
401 }
402 else if(m_order == 3) // Alexander, order 3
403 {
404 const number alpha = 0.4358665215;// root of x^3-3x^2+3/2x-1/6 = 0
405 const number tau = (1. + alpha)/2.;
406 const number b1 = -(6*alpha*alpha-16*alpha+1.)/4.;
407 const number b2 = (6*alpha*alpha-20*alpha+5.)/4.;
408 switch(m_stage)
409 {
410 case 1:
411 vSM.resize(2);
412 vSA.resize(2);
413 vSM[0] = 1.;
414 vSM[1] = -1.;
415 vSA[0] = dt * alpha;
416 vSA[1] = 0;
417 return m_Time0 + alpha * dt;
418 case 2:
419 vSM.resize(3);
420 vSA.resize(3);
421 vSM[0] = 1.;
422 vSM[1] = 0;
423 vSM[2] = -1;
424 vSA[0] = dt * alpha;
425 vSA[1] = dt * (tau-alpha);
426 vSA[2] = 0;
427 return m_Time0 + tau * dt;
428 case 3:
429 vSM.resize(4);
430 vSA.resize(4);
431 vSM[0] = 1.;
432 vSM[1] = 0;
433 vSM[2] = 0;
434 vSM[3] = -1;
435 vSA[0] = dt * alpha;
436 vSA[1] = dt * b1;
437 vSA[2] = dt * b2;
438 vSA[3] = 0;
439 return m_Time0 + dt;
440 default:
441 UG_THROW("Alexander(3) scheme has only 3 stages")
442 }
443 }
444 else if(m_order == 4) // Hairer,Wanner, order 4, L-stable DIRK
445 {
446 switch(m_stage)
447 {
448 case 1:
449 vSM.resize(2);
450 vSA.resize(2);
451 vSM[0] = 1.;
452 vSM[1] = -1.;
453 vSA[0] = dt * 1./4.;
454 vSA[1] = 0;
455 return m_Time0 + 1./4. * dt;
456 case 2:
457 vSM.resize(3);
458 vSA.resize(3);
459 vSM[0] = 1.;
460 vSM[1] = 0;
461 vSM[2] = -1;
462 vSA[0] = dt * 1./4.;
463 vSA[1] = dt * 1./2.;
464 vSA[2] = 0;
465 return m_Time0 + 3./4. * dt;
466 case 3:
467 vSM.resize(4);
468 vSA.resize(4);
469 vSM[0] = 1.;
470 vSM[1] = 0;
471 vSM[2] = 0;
472 vSM[3] = -1;
473 vSA[0] = dt * 1./4.;
474 vSA[1] = dt * (-1./25.);
475 vSA[2] = dt * 17./50.;
476 vSA[3] = 0;
477 return m_Time0 + 11./20. * dt;
478 case 4:
479 vSM.resize(5);
480 vSA.resize(5);
481 vSM[0] = 1.;
482 vSM[1] = 0;
483 vSM[2] = 0;
484 vSM[3] = 0;
485 vSM[4] = -1;
486 vSA[0] = dt * 1./4.;
487 vSA[1] = dt * (15./544.);
488 vSA[2] = dt * (-137./2720.);
489 vSA[3] = dt * (371./1360.);
490 vSA[4] = 0;
491 return m_Time0 + 1./2. * dt;
492 case 5:
493 vSM.resize(6);
494 vSA.resize(6);
495 vSM[0] = 1.;
496 vSM[1] = 0;
497 vSM[2] = 0;
498 vSM[3] = 0;
499 vSM[4] = 0;
500 vSM[5] = -1;
501 vSA[0] = dt * 1./4.;
502 vSA[1] = dt * (-85./12.);
503 vSA[2] = dt * (125./16.);
504 vSA[3] = dt * (-49./48.);
505 vSA[4] = dt * (25./24.);
506 vSA[5] = 0;
507 return m_Time0 + 1. * dt;
508 default:
509 UG_THROW("HairerWanner(4) scheme has only 5 stages")
510 }
511 }
512 else if(m_order == 3) // Crouzeix, order 3
513 {
514 const number gamma = (3. + sqrt(3.))/6.;
515 switch(m_stage)
516 {
517 case 1:
518 vSM.resize(2);
519 vSA.resize(2);
520 vSM[0] = 1.;
521 vSM[1] = -1.;
522 vSA[0] = dt * gamma;
523 vSA[1] = 0;
524 return m_Time0 + gamma * dt;
525 case 2:
526 vSM.resize(3);
527 vSA.resize(3);
528 vSM[0] = 1.;
529 vSM[1] = 0;
530 vSM[2] = -1.;
531 vSA[0] = dt * gamma;
532 vSA[1] = dt * (1. - 2*gamma);
533 vSA[2] = 0;
534 return m_Time0 + (1-gamma)*dt;
535 case 3:
536 vSM.resize(4);
537 vSA.resize(4);
538 vSM[0] = 1.;
539 vSM[1] = 0;
540 vSM[2] = 0;
541 vSM[3] = -1.;
542 vSA[0] = 0;
543 vSA[1] = dt * 1./2.;
544 vSA[2] = dt * 1./2.;
545 vSA[3] = 0;
546 return m_Time0 + dt;
547 default:
548 UG_THROW("Crouzeix 3 scheme has only 3 stages")
549 }
550 }
551 else if(m_order == 4) // Crouzeix, order 4
552 {
553 const number gamma = 1./2. + cos(M_PI/18.) / sqrt(3.);
554 const number delta = 1./(6. * (2*gamma-1) * (2*gamma-1));
555 switch(m_stage)
556 {
557 case 1:
558 vSM.resize(2);
559 vSA.resize(2);
560 vSM[0] = 1.;
561 vSM[1] = -1.;
562 vSA[0] = dt * gamma;
563 vSA[1] = 0;
564 return m_Time0 + gamma * dt;
565 case 2:
566 vSM.resize(3);
567 vSA.resize(3);
568 vSM[0] = 1.;
569 vSM[1] = 0;
570 vSM[2] = -1.;
571 vSA[0] = dt * gamma;
572 vSA[1] = dt * (1./2. - gamma);
573 vSA[2] = 0;
574 return m_Time0 + 1./2.*dt;
575 case 3:
576 vSM.resize(4);
577 vSA.resize(4);
578 vSM[0] = 1.;
579 vSM[1] = 0;
580 vSM[2] = 0;
581 vSM[3] = -1.;
582 vSA[0] = dt * gamma;
583 vSA[1] = dt * (1-4*gamma);
584 vSA[2] = dt * 2*gamma;
585 vSA[3] = 0;
586 return m_Time0 + (1-gamma)* dt;
587 case 4:
588 vSM.resize(5);
589 vSA.resize(5);
590 vSM[0] = 1.;
591 vSM[1] = 0;
592 vSM[2] = 0;
593 vSM[3] = 0;
594 vSM[4] = -1.;
595 vSA[0] = 0;
596 vSA[1] = dt * delta;
597 vSA[2] = dt * (1-2*delta);
598 vSA[3] = dt * delta;
599 vSA[4] = 0;
600 return m_Time0 + dt;
601 default:
602 UG_THROW("Crouzeix 4 scheme has only 4 stages")
603 }
604 }
605 else
606 UG_THROW("SDIRK: "<< m_order <<" missing.");
607
608}
609
610template <typename TAlgebra>
611void SDIRK<TAlgebra>::
612prepare_step(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
613 number dt)
614{
615// remember old values
616 if(m_stage == 1){
617 this->m_pPrevSol = SmartPtr<VectorTimeSeries<vector_type> >(
618 new VectorTimeSeries<vector_type>);
619 m_Time0 = prevSol->time(0);
620 this->m_futureTime = m_Time0;
621 }
622 this->m_pPrevSol->push(prevSol->solution(0)->clone(), prevSol->time(0));
623
624// remember time step size
625 this->m_dt = dt;
626
627// update scalings
628 m_lastTime = this->m_futureTime;
629
630 this->m_futureTime = update_scaling(this->m_vScaleMass, this->m_vScaleStiff,
631 this->m_dt);
632
633// prepare time step (elemDisc-wise)
634 try
635 {
636 this->m_spDomDisc->prepare_timestep(this->m_pPrevSol, this->m_futureTime);
637 }
638 UG_CATCH_THROW("ThetaTimeStep: Cannot prepare time step.");
639}
640
641template <typename TAlgebra>
642void SDIRK<TAlgebra>::
643assemble_jacobian(matrix_type& J, const vector_type& u, const GridLevel& gl)
644{
645// if(this->m_pPrevSol->size() < m_stage /*&& m_stage != num_stages()*/){
646// this->m_pPrevSol->push(u.clone(), m_lastTime);
647// }
648
649// push unknown solution to solution time series
650// ATTENTION: Here, we must cast away the constness of the solution, but note,
651// that we pass pPrevSol as a const object in assemble_... Thus,
652// the solution will not be changed there and we pop it from the
653// Solution list afterwards, such that nothing happens to u
654 // \todo: avoid this hack, use smart ptr properly
655 int DummyRefCount = 2;
656 SmartPtr<vector_type> pU(const_cast<vector_type*>(&u), &DummyRefCount);
657 this->m_pPrevSol->push(pU, this->m_futureTime);
658
659// assemble jacobian using current iterate
660 try{
661 this->m_spDomDisc->assemble_jacobian(J, this->m_pPrevSol, this->m_vScaleStiff[0], gl);
662 }UG_CATCH_THROW("SDIRK: Cannot assemble jacobian.");
663
664// pop unknown solution to solution time series
665 this->m_pPrevSol->remove_latest();
666}
667
668template <typename TAlgebra>
669void SDIRK<TAlgebra>::
670assemble_defect(vector_type& d, const vector_type& u, const GridLevel& gl)
671{
672// if(this->m_pPrevSol->size() < m_stage /*&& m_stage != num_stages()*/){
673// this->m_pPrevSol->push(u.clone(), m_lastTime);
674// }
675
676// push unknown solution to solution time series
677// ATTENTION: Here, we must cast away the constness of the solution, but note,
678// that we pass pPrevSol as a const object in assemble_... Thus,
679// the solution will not be changed there and we pop it from the
680// Solution list afterwards, such that nothing happens to u
681 // \todo: avoid this hack, use smart ptr properly
682 int DummyRefCount = 2;
683 SmartPtr<vector_type> pU(const_cast<vector_type*>(&u), &DummyRefCount);
684 this->m_pPrevSol->push(pU, this->m_futureTime);
685
686// future solution part
687 try{
688 this->m_spDomDisc->assemble_defect(d, this->m_pPrevSol, this->m_vScaleMass, this->m_vScaleStiff, gl);
689 }UG_CATCH_THROW("SDIRK: Cannot assemble defect.");
690
691// pop unknown solution to solution time series
692 this->m_pPrevSol->remove_latest();
693}
694
695template <typename TAlgebra>
696void SDIRK<TAlgebra>::
697adjust_solution(vector_type& u, const GridLevel& gl)
698{
699 PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_adjust_solution, "discretization MultiStepTimeDiscretization");
700// adjust solution
701 try{
702 this->m_spDomDisc->adjust_solution(u, this->m_futureTime, gl);
703 }UG_CATCH_THROW("ThetaTimeStep: Cannot adjust solution.");
704}
705
706template <typename TAlgebra>
707void SDIRK<TAlgebra>::
708assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl)
709{
710 UG_THROW("Not implemented")
711}
712
713template <typename TAlgebra>
714void SDIRK<TAlgebra>::
715assemble_rhs(vector_type& b, const GridLevel& gl)
716{
717 UG_THROW("Not implemented")
718}
719
720template <typename TAlgebra>
721void SDIRK<TAlgebra>::
722assemble_rhs(vector_type& b, const vector_type& u, const GridLevel& gl)
723{
724 UG_THROW("Not implemented")
725}
726
727/* Please overwrite any of the following methods, if applicable:
728template <typename TAlgebra>
729void SDIRK<TAlgebra>::
730prepare_step_elem(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
731 number dt, const GridLevel& gl)
732{
733 UG_THROW("Not implemented")
734}
735
736template <typename TAlgebra>
737void SDIRK<TAlgebra>::
738finish_step(SmartPtr<VectorTimeSeries<vector_type> > currSol)
739{
740 UG_THROW("Not implemented")
741}
742
743template <typename TAlgebra>
744void SDIRK<TAlgebra>::
745finish_step_elem(SmartPtr<VectorTimeSeries<vector_type> > currSol,
746 const GridLevel& gl)
747{
748 UG_THROW("Not implemented")
749}
750*/
751
752
753} // end namespace ug
754
755
756#endif /* __H__UG__LIB_DISC__TIME_DISC__THETA_TIME_STEP_IMPL__ */
SmartPtr
Definition smart_pointer.h:108
ug::GridLevel
Definition grid_level.h:42
ug::MultiStepTimeDiscretization::vector_type
algebra_type::vector_type vector_type
Type of algebra vector.
Definition theta_time_step.h:66
ug::MultiStepTimeDiscretization::assemble_jacobian
void assemble_jacobian(matrix_type &J, const vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:113
ug::MultiStepTimeDiscretization::finish_step
virtual void finish_step(SmartPtr< VectorTimeSeries< vector_type > > currSol)
Definition theta_time_step_impl.h:295
ug::MultiStepTimeDiscretization::finish_step_elem
virtual void finish_step_elem(SmartPtr< VectorTimeSeries< vector_type > > currSol, const GridLevel &gl)
Definition theta_time_step_impl.h:313
ug::MultiStepTimeDiscretization::prepare_step_elem
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
ug::MultiStepTimeDiscretization::prepare_step
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
ug::MultiStepTimeDiscretization::adjust_solution
void adjust_solution(vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:173
ug::MultiStepTimeDiscretization::assemble_defect
void assemble_defect(vector_type &d, const vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:143
ug::MultiStepTimeDiscretization::calc_error
void calc_error(const vector_type &u, error_vector_type *u_vtk)
Error estimator ///.
Definition theta_time_step_impl.h:261
ug::MultiStepTimeDiscretization::matrix_type
algebra_type::matrix_type matrix_type
Type of algebra matrix.
Definition theta_time_step.h:63
ug::MultiStepTimeDiscretization::assemble_linear
void assemble_linear(matrix_type &A, vector_type &b, const GridLevel &gl)
Definition theta_time_step_impl.h:184
ug::MultiStepTimeDiscretization::assemble_rhs
void assemble_rhs(vector_type &b, const vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:231
ug::SDIRK::assemble_jacobian
void assemble_jacobian(matrix_type &J, const vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:643
ug::SDIRK::vector_type
algebra_type::vector_type vector_type
Type of algebra vector.
Definition theta_time_step.h:423
ug::SDIRK::assemble_linear
void assemble_linear(matrix_type &A, vector_type &b, const GridLevel &gl)
Definition theta_time_step_impl.h:708
ug::SDIRK::assemble_rhs
void assemble_rhs(vector_type &b, const vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:722
ug::SDIRK::prepare_step
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
ug::SDIRK::set_stage
virtual void set_stage(size_t stage)
sets the stage
Definition theta_time_step_impl.h:339
ug::SDIRK::assemble_defect
void assemble_defect(vector_type &d, const vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:670
ug::SDIRK::adjust_solution
void adjust_solution(vector_type &u, const GridLevel &gl)
Definition theta_time_step_impl.h:697
ug::SDIRK::update_scaling
virtual number update_scaling(std::vector< number > &vSM, std::vector< number > &vSA, number dt)
Definition theta_time_step_impl.h:345
ug::SDIRK::matrix_type
algebra_type::matrix_type matrix_type
Type of algebra matrix.
Definition theta_time_step.h:420
ug::VectorTimeSeries
time series of solutions and corresponding time point
Definition solution_time_series.h:59
UG_CATCH_THROW
#define UG_CATCH_THROW(msg)
Definition error.h:64
UG_THROW
#define UG_THROW(msg)
Definition error.h:57
number
double number
Definition types.h:124
ug
the ug namespace
PROFILE_BEGIN_GROUP
#define PROFILE_BEGIN_GROUP(name, groups)
Definition profiler.h:255
theta_time_step.h