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_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 :
42 : namespace ug{
43 :
44 : template <typename TAlgebra>
45 0 : void MultiStepTimeDiscretization<TAlgebra>::
46 : prepare_step(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
47 : number dt)
48 : {
49 : PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_prepare_step, "discretization MultiStepTimeDiscretization");
50 : // perform checks
51 0 : if(prevSol->size() < m_prevSteps)
52 0 : 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 0 : m_pPrevSol = prevSol;
58 :
59 : // remember time step size
60 0 : m_dt = dt;
61 :
62 : // update scalings
63 0 : 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 0 : this->m_spDomDisc->prepare_timestep(m_pPrevSol, m_futureTime);
71 : }
72 0 : UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot prepare time step.");
73 0 : }
74 :
75 : template <typename TAlgebra>
76 0 : void MultiStepTimeDiscretization<TAlgebra>::
77 : prepare_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 0 : if(prevSol->size() < m_prevSteps)
83 0 : 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 0 : m_pPrevSol = prevSol;
89 :
90 : // remember time step size
91 0 : m_dt = dt;
92 :
93 : // update scalings
94 0 : 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 0 : this->m_spDomDisc->prepare_timestep(m_pPrevSol, m_futureTime, gl);
102 : }
103 0 : UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot prepare time step.");
104 :
105 : // prepare timestep element-wise
106 : try{
107 0 : this->m_spDomDisc->prepare_timestep_elem(m_pPrevSol, gl);
108 0 : }UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot prepare timestep element-wise.");
109 0 : }
110 :
111 : template <typename TAlgebra>
112 0 : void MultiStepTimeDiscretization<TAlgebra>::
113 : assemble_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 0 : if(m_pPrevSol->size() < m_prevSteps)
118 0 : 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 0 : m_pPrevSol->push(pU, m_futureTime);
131 :
132 : // assemble jacobian using current iterate
133 : try{
134 0 : this->m_spDomDisc->assemble_jacobian(J, m_pPrevSol, m_vScaleStiff[0], gl);
135 0 : }UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot assemble jacobian.");
136 :
137 : // pop unknown solution to solution time series
138 : m_pPrevSol->remove_latest();
139 0 : }
140 :
141 : template <typename TAlgebra>
142 0 : void MultiStepTimeDiscretization<TAlgebra>::
143 : assemble_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 0 : if(m_pPrevSol->size() < m_prevSteps)
148 0 : 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 0 : m_pPrevSol->push(pU, m_futureTime);
161 :
162 : // future solution part
163 : try{
164 0 : this->m_spDomDisc->assemble_defect(d, m_pPrevSol, m_vScaleMass, m_vScaleStiff, gl);
165 0 : }UG_CATCH_THROW("ThetaTimeStep: Cannot assemble defect.");
166 :
167 : // pop unknown solution to solution time series
168 : m_pPrevSol->remove_latest();
169 0 : }
170 :
171 : template <typename TAlgebra>
172 0 : void MultiStepTimeDiscretization<TAlgebra>::
173 : adjust_solution(vector_type& u, const GridLevel& gl)
174 : {
175 : PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_adjust_solution, "discretization MultiStepTimeDiscretization");
176 : // adjust solution
177 : try{
178 0 : this->m_spDomDisc->adjust_solution(u, m_futureTime, gl);
179 0 : }UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot adjust solution.");
180 0 : }
181 :
182 : template <typename TAlgebra>
183 0 : void MultiStepTimeDiscretization<TAlgebra>::
184 : assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl)
185 : {
186 : PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_assemble_linear, "discretization MultiStepTimeDiscretization");
187 : // perform checks
188 0 : if(m_pPrevSol->size() < m_prevSteps)
189 0 : 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 0 : m_pPrevSol->push(m_pPrevSol->latest(), m_futureTime);
196 :
197 : // assemble jacobian using current iterate
198 : try{
199 0 : this->m_spDomDisc->assemble_linear(A, b, m_pPrevSol, m_vScaleMass, m_vScaleStiff, gl);
200 0 : }UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot assemble jacobian.");
201 :
202 : // pop unknown solution from solution time series
203 : m_pPrevSol->remove_latest();
204 0 : }
205 :
206 : template <typename TAlgebra>
207 0 : void MultiStepTimeDiscretization<TAlgebra>::
208 : assemble_rhs(vector_type& b, const GridLevel& gl)
209 : {
210 : PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_assemble_rhs, "discretization MultiStepTimeDiscretization");
211 : // perform checks
212 0 : if(m_pPrevSol->size() < m_prevSteps)
213 0 : 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 0 : m_pPrevSol->push(m_pPrevSol->latest(), m_futureTime);
219 :
220 : // assemble jacobian using current iterate
221 : try{
222 0 : this->m_spDomDisc->assemble_rhs(b, m_pPrevSol, m_vScaleMass, m_vScaleStiff, gl);
223 0 : }UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot assemble jacobian.");
224 :
225 : // pop unknown solution from solution time series
226 : m_pPrevSol->remove_latest();
227 0 : }
228 :
229 : template <typename TAlgebra>
230 0 : void MultiStepTimeDiscretization<TAlgebra>::
231 : assemble_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 0 : if(m_pPrevSol->size() < m_prevSteps)
236 0 : 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 0 : m_pPrevSol->push(pU, m_futureTime);
249 :
250 : // assemble jacobian using current iterate
251 : try{
252 0 : this->m_spDomDisc->assemble_rhs(b, m_pPrevSol, m_vScaleMass, m_vScaleStiff, gl);
253 0 : }UG_CATCH_THROW("ThetaTimeStep: Cannot assemble jacobian.");
254 :
255 : // pop unknown solution to solution time series
256 : m_pPrevSol->remove_latest();
257 0 : }
258 :
259 : template<typename TAlgebra>
260 0 : void MultiStepTimeDiscretization<TAlgebra>::
261 : calc_error(const vector_type& u, error_vector_type* u_vtk)
262 : {
263 : PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_calc_error, "discretization MultiStepTimeDiscretization");
264 : // perform checks
265 0 : if(m_pPrevSol->size() < m_prevSteps)
266 0 : 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 0 : 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 0 : this->m_spDomDisc->calc_error(m_pPrevSol, m_vScaleMass, m_vScaleStiff, gl, u_vtk);
285 : }
286 0 : UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot assemble error estimators.");
287 :
288 : // pop unknown solution to solution time series
289 : m_pPrevSol->remove_latest();
290 0 : }
291 :
292 :
293 : template <typename TAlgebra>
294 0 : void MultiStepTimeDiscretization<TAlgebra>::
295 : finish_step(SmartPtr<VectorTimeSeries<vector_type> > currSol)
296 : {
297 : // perform checks whether 'currSol' is a solutionTimeSeries only with the new values
298 0 : if (currSol->time(0) != m_futureTime)
299 0 : 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 0 : this->m_spDomDisc->finish_timestep(currSol);
307 : }
308 0 : UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot finish timestep.");
309 0 : }
310 :
311 : template <typename TAlgebra>
312 0 : void MultiStepTimeDiscretization<TAlgebra>::
313 : finish_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 0 : if(currSol->time(0) != m_futureTime)
318 0 : 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 0 : this->m_spDomDisc->finish_timestep(currSol, gl);
326 : }
327 0 : UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot finish timestep.");
328 :
329 : // finish timestep element-wise using the current solution
330 : try{
331 0 : this->m_spDomDisc->finish_timestep_elem(currSol, gl);
332 0 : }UG_CATCH_THROW("MultiStepTimeDiscretization: Cannot finish timestep element-wise.");
333 0 : }
334 :
335 : //////////////////////////////////////////////////////////////////////////////
336 : //////////////////////////////////////////////////////////////////////////////
337 :
338 : template <typename TAlgebra>
339 0 : void SDIRK<TAlgebra>::set_stage(size_t stage)
340 : {
341 0 : m_stage = stage;
342 0 : }
343 :
344 : template <typename TAlgebra>
345 0 : number SDIRK<TAlgebra>::update_scaling(std::vector<number>& vSM,
346 : std::vector<number>& vSA,
347 : number dt)
348 : {
349 0 : if(m_order == 1) // Mittelpunkt
350 : {
351 0 : switch(m_stage)
352 : {
353 0 : case 1:
354 0 : vSM.resize(2);
355 0 : vSA.resize(2);
356 0 : vSM[0] = 1.;
357 0 : vSM[1] = -1.;
358 0 : vSA[0] = dt * 1./2.;
359 0 : vSA[1] = 0;
360 0 : return m_Time0 + 1./2. * dt;
361 0 : case 2:
362 0 : vSM.resize(3);
363 0 : vSA.resize(3);
364 0 : vSM[0] = 1.;
365 0 : vSM[1] = 0;
366 0 : vSM[2] = -1.;
367 0 : vSA[0] = 0;
368 0 : vSA[1] = dt;
369 0 : vSA[2] = 0;
370 0 : return m_Time0 + dt;
371 0 : default:
372 0 : UG_THROW("Midpoint scheme has only 2 stages")
373 : }
374 : }
375 0 : else if(m_order == 2) // Alexander, order 2
376 : {
377 : const number gamma = 1 - 1. / sqrt(2.);
378 0 : switch(m_stage)
379 : {
380 0 : case 1:
381 0 : vSM.resize(2);
382 0 : vSA.resize(2);
383 0 : vSM[0] = 1.;
384 0 : vSM[1] = -1.;
385 0 : vSA[0] = dt * gamma;
386 0 : vSA[1] = 0;
387 0 : return m_Time0 + gamma * dt;
388 0 : case 2:
389 0 : vSM.resize(3);
390 0 : vSA.resize(3);
391 0 : vSM[0] = 1.;
392 0 : vSM[1] = 0;
393 0 : vSM[2] = -1;
394 0 : vSA[0] = dt * gamma;
395 0 : vSA[1] = dt * (1. - gamma);
396 0 : vSA[2] = 0;
397 0 : return m_Time0 + dt;
398 0 : default:
399 0 : UG_THROW("Alexander(2) scheme has only 2 stages")
400 : }
401 : }
402 0 : 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 0 : switch(m_stage)
409 : {
410 0 : case 1:
411 0 : vSM.resize(2);
412 0 : vSA.resize(2);
413 0 : vSM[0] = 1.;
414 0 : vSM[1] = -1.;
415 0 : vSA[0] = dt * alpha;
416 0 : vSA[1] = 0;
417 0 : return m_Time0 + alpha * dt;
418 0 : case 2:
419 0 : vSM.resize(3);
420 0 : vSA.resize(3);
421 0 : vSM[0] = 1.;
422 0 : vSM[1] = 0;
423 0 : vSM[2] = -1;
424 0 : vSA[0] = dt * alpha;
425 0 : vSA[1] = dt * (tau-alpha);
426 0 : vSA[2] = 0;
427 0 : return m_Time0 + tau * dt;
428 0 : case 3:
429 0 : vSM.resize(4);
430 0 : vSA.resize(4);
431 0 : vSM[0] = 1.;
432 0 : vSM[1] = 0;
433 0 : vSM[2] = 0;
434 0 : vSM[3] = -1;
435 0 : vSA[0] = dt * alpha;
436 0 : vSA[1] = dt * b1;
437 0 : vSA[2] = dt * b2;
438 0 : vSA[3] = 0;
439 0 : return m_Time0 + dt;
440 0 : default:
441 0 : UG_THROW("Alexander(3) scheme has only 3 stages")
442 : }
443 : }
444 0 : else if(m_order == 4) // Hairer,Wanner, order 4, L-stable DIRK
445 : {
446 0 : switch(m_stage)
447 : {
448 0 : case 1:
449 0 : vSM.resize(2);
450 0 : vSA.resize(2);
451 0 : vSM[0] = 1.;
452 0 : vSM[1] = -1.;
453 0 : vSA[0] = dt * 1./4.;
454 0 : vSA[1] = 0;
455 0 : return m_Time0 + 1./4. * dt;
456 0 : case 2:
457 0 : vSM.resize(3);
458 0 : vSA.resize(3);
459 0 : vSM[0] = 1.;
460 0 : vSM[1] = 0;
461 0 : vSM[2] = -1;
462 0 : vSA[0] = dt * 1./4.;
463 0 : vSA[1] = dt * 1./2.;
464 0 : vSA[2] = 0;
465 0 : return m_Time0 + 3./4. * dt;
466 0 : case 3:
467 0 : vSM.resize(4);
468 0 : vSA.resize(4);
469 0 : vSM[0] = 1.;
470 0 : vSM[1] = 0;
471 0 : vSM[2] = 0;
472 0 : vSM[3] = -1;
473 0 : vSA[0] = dt * 1./4.;
474 0 : vSA[1] = dt * (-1./25.);
475 0 : vSA[2] = dt * 17./50.;
476 0 : vSA[3] = 0;
477 0 : return m_Time0 + 11./20. * dt;
478 0 : case 4:
479 0 : vSM.resize(5);
480 0 : vSA.resize(5);
481 0 : vSM[0] = 1.;
482 0 : vSM[1] = 0;
483 0 : vSM[2] = 0;
484 0 : vSM[3] = 0;
485 0 : vSM[4] = -1;
486 0 : vSA[0] = dt * 1./4.;
487 0 : vSA[1] = dt * (15./544.);
488 0 : vSA[2] = dt * (-137./2720.);
489 0 : vSA[3] = dt * (371./1360.);
490 0 : vSA[4] = 0;
491 0 : return m_Time0 + 1./2. * dt;
492 0 : case 5:
493 0 : vSM.resize(6);
494 0 : vSA.resize(6);
495 0 : vSM[0] = 1.;
496 0 : vSM[1] = 0;
497 0 : vSM[2] = 0;
498 0 : vSM[3] = 0;
499 0 : vSM[4] = 0;
500 0 : vSM[5] = -1;
501 0 : vSA[0] = dt * 1./4.;
502 0 : vSA[1] = dt * (-85./12.);
503 0 : vSA[2] = dt * (125./16.);
504 0 : vSA[3] = dt * (-49./48.);
505 0 : vSA[4] = dt * (25./24.);
506 0 : vSA[5] = 0;
507 0 : return m_Time0 + 1. * dt;
508 0 : default:
509 0 : 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 0 : UG_THROW("SDIRK: "<< m_order <<" missing.");
607 :
608 : }
609 :
610 : template <typename TAlgebra>
611 0 : void SDIRK<TAlgebra>::
612 : prepare_step(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
613 : number dt)
614 : {
615 : // remember old values
616 0 : if(m_stage == 1){
617 0 : this->m_pPrevSol = SmartPtr<VectorTimeSeries<vector_type> >(
618 0 : new VectorTimeSeries<vector_type>);
619 0 : m_Time0 = prevSol->time(0);
620 0 : this->m_futureTime = m_Time0;
621 : }
622 0 : this->m_pPrevSol->push(prevSol->solution(0)->clone(), prevSol->time(0));
623 :
624 : // remember time step size
625 0 : this->m_dt = dt;
626 :
627 : // update scalings
628 0 : m_lastTime = this->m_futureTime;
629 :
630 0 : 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 0 : this->m_spDomDisc->prepare_timestep(this->m_pPrevSol, this->m_futureTime);
637 : }
638 0 : UG_CATCH_THROW("ThetaTimeStep: Cannot prepare time step.");
639 0 : }
640 :
641 : template <typename TAlgebra>
642 0 : void SDIRK<TAlgebra>::
643 : assemble_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 0 : this->m_pPrevSol->push(pU, this->m_futureTime);
658 :
659 : // assemble jacobian using current iterate
660 : try{
661 0 : this->m_spDomDisc->assemble_jacobian(J, this->m_pPrevSol, this->m_vScaleStiff[0], gl);
662 0 : }UG_CATCH_THROW("SDIRK: Cannot assemble jacobian.");
663 :
664 : // pop unknown solution to solution time series
665 : this->m_pPrevSol->remove_latest();
666 0 : }
667 :
668 : template <typename TAlgebra>
669 0 : void SDIRK<TAlgebra>::
670 : assemble_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 0 : this->m_pPrevSol->push(pU, this->m_futureTime);
685 :
686 : // future solution part
687 : try{
688 0 : this->m_spDomDisc->assemble_defect(d, this->m_pPrevSol, this->m_vScaleMass, this->m_vScaleStiff, gl);
689 0 : }UG_CATCH_THROW("SDIRK: Cannot assemble defect.");
690 :
691 : // pop unknown solution to solution time series
692 : this->m_pPrevSol->remove_latest();
693 0 : }
694 :
695 : template <typename TAlgebra>
696 0 : void SDIRK<TAlgebra>::
697 : adjust_solution(vector_type& u, const GridLevel& gl)
698 : {
699 : PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_adjust_solution, "discretization MultiStepTimeDiscretization");
700 : // adjust solution
701 : try{
702 0 : this->m_spDomDisc->adjust_solution(u, this->m_futureTime, gl);
703 0 : }UG_CATCH_THROW("ThetaTimeStep: Cannot adjust solution.");
704 0 : }
705 :
706 : template <typename TAlgebra>
707 0 : void SDIRK<TAlgebra>::
708 : assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl)
709 : {
710 0 : UG_THROW("Not implemented")
711 : }
712 :
713 : template <typename TAlgebra>
714 0 : void SDIRK<TAlgebra>::
715 : assemble_rhs(vector_type& b, const GridLevel& gl)
716 : {
717 0 : UG_THROW("Not implemented")
718 : }
719 :
720 : template <typename TAlgebra>
721 0 : void SDIRK<TAlgebra>::
722 : assemble_rhs(vector_type& b, const vector_type& u, const GridLevel& gl)
723 : {
724 0 : UG_THROW("Not implemented")
725 : }
726 :
727 : /* Please overwrite any of the following methods, if applicable:
728 : template <typename TAlgebra>
729 : void SDIRK<TAlgebra>::
730 : prepare_step_elem(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
731 : number dt, const GridLevel& gl)
732 : {
733 : UG_THROW("Not implemented")
734 : }
735 :
736 : template <typename TAlgebra>
737 : void SDIRK<TAlgebra>::
738 : finish_step(SmartPtr<VectorTimeSeries<vector_type> > currSol)
739 : {
740 : UG_THROW("Not implemented")
741 : }
742 :
743 : template <typename TAlgebra>
744 : void SDIRK<TAlgebra>::
745 : finish_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__ */
|
