GCC Code Coverage Report

Directory: ./
File: Solver/bdf.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 142 397 35.8%
Branches: 36 138 26.1%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors:
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/bdf.hpp"
21
22 namespace felisce {
23 102 Bdf::Bdf():
24 102 m_order(0),
25 102 m_numberOfComp(1),
26 102 m_size(0),
27 102 m_build_order_2(false),
28 102 m_build_order_3(false),
29 102 m_build_RHS(false) {
30
31 102 }
32
33 102 Bdf::~Bdf() {
34 102 m_alpha.clear();
35 102 m_beta.clear();
36
2/2
✓ Branch 0 taken 102 times.
✓ Branch 1 taken 102 times.
 204 for (felInt i=0 ; i< m_numberOfComp; i++) {
37
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 90 times.
 102 if (m_build_order_2)
38 12 m_sol_n_1[i].destroy();
39
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 102 times.
 102 if (m_build_order_3)
40 ✗ m_sol_n_2[i].destroy();
41
2/2
✓ Branch 0 taken 98 times.
✓ Branch 1 taken 4 times.
 102 if (m_build_RHS) {
42 98 m_sol_n[i].destroy();
43 98 m_rhs[i].destroy();
44 }
45 }
46 102 }
47
48 116 void Bdf::defineOrder(int n, int nComp) {
49 116 m_order = n;
50 116 m_numberOfComp=nComp;
51 116 m_alpha.resize(n+1);
52 116 m_beta.resize(n);
53 116 m_sol_n.resize(nComp);
54 116 m_sol_n_1.resize(nComp);
55 116 m_sol_n_2.resize(nComp);
56 116 m_rhs.resize(nComp);
57 116 m_solExtrapol.resize(nComp);
58
59
2/4
✓ Branch 0 taken 116 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 116 times.
 116 if ( n <= 0 || n > BDF_MAX_ORDER ) {
60 ✗ FEL_ERROR(" We support BDF order from 1 to 3 \n");
61 }
62
63
2/4
✓ Branch 0 taken 108 times.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 116 switch (n) {
64 108 case 1:
65 108 m_alpha[0] = 1.0; // Backward Euler
66 108 m_alpha[1] = 1.0;
67 108 m_beta[0] = 1.0;// u^{n+1} \approx u^n
68 108 break;
69 8 case 2:
70 8 m_alpha[0] = 3.0 / 2.0;
71 8 m_alpha[1] = 2.0;
72 8 m_alpha[2] = -1.0 / 2.0;
73 8 m_beta[0] = 2.0;
74 8 m_beta[1] = -1.0;
75 8 break;
76 ✗ case 3:
77 ✗ m_alpha[0] = 11.0 / 6.0;
78 ✗ m_alpha[1] = 3.0;
79 ✗ m_alpha[2] = -3.0 / 2.0;
80 ✗ m_alpha[3] = 1.0 / 3.0;
81 ✗ m_beta[0] = 3.0;
82 ✗ m_beta[1] = -3.0;
83 ✗ m_beta[2] = 1.0;
84 ✗ break;
85 }
86 116 }
87
88
89 94 void Bdf::initialize(const PetscVector& sol_0) {
90
2/4
✓ Branch 0 taken 86 times.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 94 switch (m_order) {
91 86 case 1:
92 //!Get size.
93 86 sol_0.getSize(&m_size);
94 86 m_sol_n[0].duplicateFrom(sol_0);
95 //!Copy because you want to keep sol_0.
96 86 m_sol_n[0].copyFrom(sol_0);
97
98 86 m_rhs[0].duplicateFrom(sol_0);
99 86 m_build_RHS = true;
100 86 break;
101 8 case 2:
102 8 sol_0.getSize(&m_size);
103 8 m_sol_n[0].duplicateFrom(sol_0);
104 8 m_sol_n[0].copyFrom(sol_0);
105
106 8 m_sol_n_1[0].duplicateFrom(sol_0);
107 8 m_sol_n_1[0].copyFrom(sol_0);
108 8 m_build_order_2 = true;
109
110 8 m_rhs[0].duplicateFrom(sol_0);
111 8 m_build_RHS = true;
112 8 break;
113 ✗ case 3:
114 ✗ sol_0.getSize(&m_size);
115 ✗ m_sol_n[0].duplicateFrom(sol_0);
116 ✗ m_sol_n[0].copyFrom(sol_0);
117
118 ✗ m_sol_n_1[0].duplicateFrom(sol_0);
119 ✗ m_sol_n_1[0].copyFrom(sol_0);
120 ✗ m_build_order_2 = true;
121
122 ✗ m_sol_n_2[0].duplicateFrom(sol_0);
123 ✗ m_sol_n_2[0].copyFrom(sol_0);
124 ✗ m_build_order_3 = true;
125
126 ✗ m_rhs[0].duplicateFrom(sol_0);
127 ✗ m_build_RHS = true;
128 ✗ break;
129 ✗ default:
130 ✗ FEL_ERROR("unknown BDF order");
131 ✗ break;
132 }
133 94 }
134
135 ✗ void Bdf::initialize(std::vector<PetscVector>& sol_0) {
136 ✗ switch (m_order) {
137 ✗ case 1:
138 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
139 //!Get size.
140 ✗ sol_0[i].getSize(&m_size);
141 ✗ m_sol_n[i].duplicateFrom(sol_0[i]);
142 //!Copy because you want to keep sol_0.
143 ✗ m_sol_n[i].copyFrom(sol_0[i]);
144
145 ✗ m_rhs[i].duplicateFrom(sol_0[i]);
146 }
147 ✗ m_build_RHS = true;
148 ✗ break;
149 ✗ case 2:
150 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
151 ✗ sol_0[i].getSize(&m_size);
152 ✗ m_sol_n[i].duplicateFrom(sol_0[i]);
153 ✗ m_sol_n[i].copyFrom(sol_0[i]);
154
155 ✗ m_sol_n_1[i].duplicateFrom(sol_0[i]);
156 ✗ m_sol_n_1[i].copyFrom(sol_0[i]);
157 }
158 ✗ m_build_order_2 = true;
159 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
160 ✗ m_rhs[i].duplicateFrom(sol_0[i]);
161 }
162 ✗ m_build_RHS = true;
163 ✗ break;
164 ✗ case 3:
165 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
166 ✗ sol_0[i].getSize(&m_size);
167 ✗ m_sol_n[i].duplicateFrom(sol_0[i]);
168 ✗ m_sol_n[i].copyFrom(sol_0[i]);
169
170 ✗ m_sol_n_1[i].duplicateFrom(sol_0[i]);
171 ✗ m_sol_n_1[i].copyFrom(sol_0[i]);
172 }
173 ✗ m_build_order_2 = true;
174 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
175 ✗ m_sol_n_2[i].duplicateFrom(sol_0[i]);
176 ✗ m_sol_n_2[i].copyFrom(sol_0[i]);
177 }
178 ✗ m_build_order_3 = true;
179 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
180 ✗ m_rhs[i].duplicateFrom(sol_0[i]);
181 }
182 ✗ m_build_RHS = true;
183 ✗ break;
184 ✗ default:
185 ✗ FEL_ERROR("unknown BDF order");
186 ✗ break;
187 }
188 }
189
190
191 20 void Bdf::initialize(const PetscVector& sol_0, const PetscVector& sol_1) {
192 20 sol_0.getSize(&m_size);
193
194 20 m_sol_n[0].duplicateFrom(sol_1);
195 //VecSwap(m_sol_n,sol_1);
196 20 m_sol_n[0].copyFrom(sol_1);
197
198 20 m_sol_n_1[0].duplicateFrom(sol_0);
199 //VecSwap(m_sol_n_1,sol_0);
200 20 m_sol_n_1[0].copyFrom(sol_0);
201
202 20 m_build_order_2 = true;
203
204 20 m_rhs[0].duplicateFrom(sol_0);
205 20 m_build_RHS = true;
206 20 }
207
208
209 ✗ void Bdf::initialize(std::vector<PetscVector>& sol_0, std::vector<PetscVector>& sol_1) {
210
211 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
212 ✗ sol_0[i].getSize(&m_size);
213
214 ✗ m_sol_n[i].duplicateFrom(sol_1[i]);
215 ✗ m_sol_n[i].copyFrom(sol_1[i]);
216
217 ✗ m_sol_n_1[i].duplicateFrom(sol_0[i]);
218 ✗ m_sol_n_1[i].copyFrom(sol_0[i]);
219 }
220 ✗ m_build_order_2 = true;
221 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
222 ✗ m_rhs[i].duplicateFrom(sol_0[i]);
223 }
224 ✗ m_build_RHS = true;
225 }
226
227 ✗ void Bdf::initialize(const PetscVector& sol_0, const PetscVector& sol_1, const PetscVector& sol_2) {
228
229 ✗ sol_0.getSize(&m_size);
230
231 ✗ m_sol_n[0].duplicateFrom(sol_2);
232 //VecSwap(m_sol_n,sol_2);
233 ✗ m_sol_n[0].copyFrom(sol_2);
234
235 ✗ m_sol_n_1[0].duplicateFrom(sol_1);
236 //VecSwap(m_sol_n_1,sol_1);
237 ✗ m_sol_n_1[0].copyFrom(sol_1);
238
239 ✗ m_sol_n_2[0].duplicateFrom(sol_0);
240 //VecSwap(m_sol_n_2,sol_0);
241 ✗ m_sol_n_2[0].copyFrom(sol_0);
242
243 ✗ m_build_order_3 = true;
244
245 ✗ m_rhs[0].duplicateFrom(sol_0);
246 ✗ m_build_RHS = true;
247 }
248
249 ✗ void Bdf::initialize(std::vector<PetscVector>& sol_0, std::vector<PetscVector>& sol_1, std::vector<PetscVector>& sol_2) {
250
251 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
252 ✗ sol_0[i].getSize(&m_size);
253
254 ✗ m_sol_n[i].duplicateFrom(sol_2[i]);
255 //VecSwap(m_sol_n,sol_2);
256 ✗ m_sol_n[i].copyFrom(sol_2[i]);
257
258 ✗ m_sol_n_1[i].duplicateFrom(sol_1[i]);
259 //VecSwap(m_sol_n_1,sol_1);
260 ✗ m_sol_n_1[i].copyFrom(sol_1[i]);
261
262 ✗ m_sol_n_2[i].duplicateFrom(sol_0[i]);
263 //VecSwap(m_sol_n_2,sol_0);
264 ✗ m_sol_n_2[i].copyFrom(sol_0[i]);
265 }
266 ✗ m_build_order_3 = true;
267 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
268 ✗ m_rhs[i].duplicateFrom(sol_0[i]);
269 }
270 ✗ m_build_RHS = true;
271 }
272
273 18 void Bdf::reinitialize(int order, const PetscVector& sol0, const PetscVector& sol1, const PetscVector& sol2) {
274 // destroy everything
275 18 m_alpha.clear();
276 18 m_beta.clear();
277
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 2 times.
 18 if (m_build_order_2)
278 16 m_sol_n_1[0].destroy();
279
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
 18 if (m_build_order_3)
280 ✗ m_sol_n_2[0].destroy();
281
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 2 times.
 18 if (m_build_RHS) {
282 16 m_sol_n[0].destroy();
283 16 m_rhs[0].destroy();
284 }
285
286 18 m_build_order_2 = false;
287 18 m_build_order_3 = false;
288 18 m_build_RHS = false;
289
290 // reallocate everything
291 18 defineOrder(order);
292
293
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
 18 if(sol2.isNotNull())
294 ✗ initialize(sol0, sol1, sol2);
295
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 else if(sol1.isNotNull())
296 18 initialize(sol0, sol1);
297 else
298 ✗ initialize(sol0);
299 18 }
300
301 2129 void Bdf::update(PetscVector& new_sol_n) {
302 /*switch (m_order) {
303 case 1:
304 sol_n[0].copyFrom(new_sol_n);
305 break;
306 case 2:
307 sol_n_1[0].copyFrom(m_sol_n[0]);
308 sol_n[0].copyFrom(new_sol_n);
309 break;
310 case 3:
311 sol_n_2[0].copyFrom(m_sol_n_1[0]);
312 sol_n_1[0].copyFrom(m_sol_n[0]);
313 sol_n[0].copyFrom(new_sol_n);
314 break;
315 default:
316 FEL_ERROR("BDF not yet implemented for order grater than 3.");
317 break;
318 }*/
319
320 // use m_build_order instead of m_order in case someone wants to initialize the bdf with more time step
321
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2129 times.
 2129 if(m_build_order_3) {
322 ✗ m_sol_n_2[0].copyFrom(m_sol_n_1[0]);
323 ✗ m_sol_n_1[0].copyFrom(m_sol_n[0]);
324 ✗ m_sol_n[0].copyFrom(new_sol_n);
325
2/2
✓ Branch 0 taken 127 times.
✓ Branch 1 taken 2002 times.
 2129 } else if(m_build_order_2) {
326 127 m_sol_n_1[0].copyFrom(m_sol_n[0]);
327 127 m_sol_n[0].copyFrom(new_sol_n);
328
1/2
✓ Branch 0 taken 2002 times.
✗ Branch 1 not taken.
 2002 } else if(m_build_RHS) {
329 2002 m_sol_n[0].copyFrom(new_sol_n);
330 } else {
331 ✗ FEL_ERROR("BDF not yet implemented for order grater than 3.");
332 }
333 2129 }
334
335 ✗ void Bdf::update(std::vector<PetscVector>& new_sol_n) {
336 /*switch (m_order) {
337 case 1:
338 for (felInt i=0 ; i< m_numberOfComp; i++) {
339 m_sol_n[i].copyFrom(new_sol_n[i]);
340 }
341 break;
342 case 2:
343 for (felInt i=0 ; i< m_numberOfComp; i++) {
344 m_sol_n_1[i].copyFrom(m_sol_n[i]);
345 m_sol_n[i].copyFrom(new_sol_n[i]);
346 }
347 break;
348 case 3:
349 for (felInt i=0 ; i< m_numberOfComp; i++) {
350 m_sol_n_2[i].copyFrom(m_sol_n_1[i]);
351 m_sol_n_1[i].copyFrom(m_sol_n[i]);
352 m_sol_n[i].copyFrom(new_sol_n[i]);
353 }
354 break;
355 default:
356 FEL_ERROR("BDF not yet implemented for order grater than 3.");
357 break;
358 }*/
359
360 // use m_build_order instead of m_order in case someone wants to initialize the bdf with more time step
361 ✗ if(m_build_order_3) {
362 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
363 ✗ m_sol_n_2[i].copyFrom(m_sol_n_1[i]);
364 ✗ m_sol_n_1[i].copyFrom(m_sol_n[i]);
365 ✗ m_sol_n[i].copyFrom(new_sol_n[i]);
366 }
367 ✗ } else if(m_build_order_2) {
368 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
369 ✗ m_sol_n_1[i].copyFrom(m_sol_n[i]);
370 ✗ m_sol_n[i].copyFrom(new_sol_n[i]);
371 }
372 ✗ } else if(m_build_RHS) {
373 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
374 ✗ m_sol_n[i].copyFrom(new_sol_n[i]);
375 }
376 } else {
377 ✗ FEL_ERROR("BDF not yet implemented for order grater than 3.");
378 }
379 }
380
381
382 2248 void Bdf::computeRHSTime(double dt, PetscVector& RHSTime) {
383 2248 RHSTime.zeroEntries();
384
2/4
✓ Branch 0 taken 2128 times.
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2248 switch (m_order) {
385 2128 case 1:
386 //RHSTime = m_alpha[0](1)*m_sol_n
387 2128 RHSTime.axpy(m_alpha[1]/dt,m_sol_n[0]);
388 2128 break;
389 120 case 2:
390 //RHSTime = m_alpha[0](1)*m_sol_n
391 120 RHSTime.axpy(m_alpha[1]/dt,m_sol_n[0]);
392 //RHSTime = m_alpha[0](1)*m_sol_n + m_alpha[0](2)*m_sol_n_1
393 120 RHSTime.axpy(m_alpha[2]/dt,m_sol_n_1[0]);
394 120 break;
395 ✗ case 3:
396 ✗ RHSTime.axpy(m_alpha[1]/dt,m_sol_n[0]);
397 ✗ RHSTime.axpy(m_alpha[2]/dt,m_sol_n_1[0]);
398 ✗ RHSTime.axpy(m_alpha[3]/dt,m_sol_n_2[0]);
399 ✗ break;
400 ✗ default:
401 ✗ FEL_ERROR("BDF not yet implemented for order grater than 3.");
402 ✗ break;
403 }
404 2248 }
405
406 ✗ void Bdf::computeRHSTime(double dt, std::vector<PetscVector>& RHSTime) {
407 ✗ for (felInt i=0 ; i< m_numberOfComp; i++)
408 ✗ RHSTime[i].zeroEntries();
409 ✗ switch (m_order) {
410 ✗ case 1:
411 //RHSTime = m_alpha[0](1)*m_sol_n
412 ✗ for (felInt i=0 ; i< m_numberOfComp; i++)
413 ✗ RHSTime[i].axpy(m_alpha[1]/dt,m_sol_n[i]);
414 ✗ break;
415 ✗ case 2:
416 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
417 //RHSTime = m_alpha[0](1)*m_sol_n
418 ✗ RHSTime[i].axpy(m_alpha[1]/dt,m_sol_n[i]);
419 //RHSTime = m_alpha[0](1)*m_sol_n + m_alpha[0](2)*m_sol_n_1
420 ✗ RHSTime[i].axpy(m_alpha[2]/dt,m_sol_n_1[i]);
421 }
422 ✗ break;
423 ✗ case 3:
424 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
425 ✗ RHSTime[i].axpy(m_alpha[1]/dt,m_sol_n[i]);
426 ✗ RHSTime[i].axpy(m_alpha[2]/dt,m_sol_n_1[i]);
427 ✗ RHSTime[i].axpy(m_alpha[3]/dt,m_sol_n_2[i]);
428 }
429 ✗ break;
430 ✗ default:
431 ✗ FEL_ERROR("BDF not yet implemented for order grater than 3.");
432 ✗ break;
433 }
434 }
435
436 20 void Bdf::computeRHSTimeByPart(double dt, std::vector<PetscVector>& RHSTimeByPart) {
437
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 20 times.
 20 if(m_order != static_cast<felInt>(RHSTimeByPart.size()))
438 ✗ FEL_ERROR("The bdf order and the size of the rhs part std::vector are not the same");
439
440
2/2
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 20 times.
 40 for(std::size_t i=0; i<RHSTimeByPart.size(); ++i)
441 20 RHSTimeByPart[i].zeroEntries();
442
443
1/4
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 20 switch (m_order) {
444 20 case 1:
445 // RHSTime[0] = m_alpha[0](1)*m_sol_n
446 20 RHSTimeByPart[0].axpy(m_alpha[1]/dt, m_sol_n[0]);
447 20 break;
448
449 ✗ case 2:
450 // RHSTime[0] = m_alpha[0](1)*m_sol_n
451 ✗ RHSTimeByPart[0].axpy(m_alpha[1]/dt, m_sol_n[0]);
452
453 //RHSTime[1] = m_alpha[0](2)*m_sol_n_1
454 ✗ RHSTimeByPart[1].axpy(m_alpha[2]/dt, m_sol_n_1[0]);
455 ✗ break;
456
457 ✗ case 3:
458 ✗ RHSTimeByPart[0].axpy(m_alpha[1]/dt, m_sol_n[0]);
459 ✗ RHSTimeByPart[1].axpy(m_alpha[2]/dt, m_sol_n_1[0]);
460 ✗ RHSTimeByPart[2].axpy(m_alpha[3]/dt, m_sol_n_2[0]);
461 ✗ break;
462
463 ✗ default:
464 ✗ FEL_ERROR("BDF not yet implemented for order grater than 3.");
465 ✗ break;
466 }
467 20 }
468
469
470 2248 void Bdf::computeRHSTime(double dt) {
471 2248 computeRHSTime(dt,m_rhs[0]);
472 2248 }
473
474 // Evaluates extrapolation
475 2068 void Bdf::extrapolate( PetscVector& extrap) {
476 2068 extrap.zeroEntries();
477
2/4
✓ Branch 0 taken 1948 times.
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2068 switch (m_order) {
478 1948 case 1:
479 1948 extrap.axpy(m_beta[0],m_sol_n[0]);
480 1948 break;
481 120 case 2:
482 120 extrap.axpy(m_beta[0],m_sol_n[0]);
483 120 extrap.axpy(m_beta[1],m_sol_n_1[0]);
484 120 break;
485 ✗ case 3:
486 ✗ extrap.axpy(m_beta[0],m_sol_n[0]);
487 ✗ extrap.axpy(m_beta[1],m_sol_n_1[0]);
488 ✗ extrap.axpy(m_beta[2],m_sol_n_2[0]);
489 ✗ break;
490 ✗ default:
491 ✗ FEL_ERROR("unknown BDF order");
492 ✗ break;
493 }
494 2068 }
495
496 ✗ void Bdf::extrapolate( std::vector<PetscVector>& extrap) {
497 ✗ for (felInt i=0 ; i< m_numberOfComp; i++)
498 ✗ extrap[i].zeroEntries();
499 ✗ switch (m_order) {
500 ✗ case 1:
501 ✗ for (felInt i=0 ; i< m_numberOfComp; i++)
502 ✗ extrap[i].axpy(m_beta[0],m_sol_n[i]);
503 ✗ break;
504 ✗ case 2:
505 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
506 ✗ extrap[i].axpy(m_beta[0],m_sol_n[i]);
507 ✗ extrap[i].axpy(m_beta[1],m_sol_n_1[i]);
508 }
509 ✗ break;
510 ✗ case 3:
511 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
512 ✗ extrap[i].axpy(m_beta[0],m_sol_n[i]);
513 ✗ extrap[i].axpy(m_beta[1],m_sol_n_1[i]);
514 ✗ extrap[i].axpy(m_beta[2],m_sol_n_2[i]);
515 }
516 ✗ break;
517 ✗ default:
518 ✗ FEL_ERROR("unknown BDF order");
519 ✗ break;
520 }
521 }
522
523 20 void Bdf::extrapolateByPart(std::vector<PetscVector>& partExtrap) {
524
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 20 times.
 20 if(m_order != static_cast<felInt>(partExtrap.size()))
525 ✗ FEL_ERROR("The bdf order and the size of the extrapolation part std::vector are not the same");
526
527
2/2
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 20 times.
 40 for(std::size_t i=0; i<partExtrap.size(); i++)
528 20 partExtrap[i].zeroEntries();
529
530
1/4
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 20 switch (m_order) {
531 20 case 1:
532 20 partExtrap[0].axpy(m_beta[0], m_sol_n[0]);
533 20 break;
534 ✗ case 2:
535 ✗ partExtrap[0].axpy(m_beta[0], m_sol_n[0]);
536 ✗ partExtrap[1].axpy(m_beta[1], m_sol_n_1[0]);
537 ✗ break;
538 ✗ case 3:
539 ✗ partExtrap[0].axpy(m_beta[0], m_sol_n[0]);
540 ✗ partExtrap[1].axpy(m_beta[1], m_sol_n_1[0]);
541 ✗ partExtrap[2].axpy(m_beta[2], m_sol_n_2[0]);
542 ✗ break;
543 ✗ default:
544 ✗ FEL_ERROR("unknown BDF order");
545 ✗ break;
546 }
547 20 }
548
549 // evaluates time derivative std::vector a posteriori; to call before update
550 ✗ void Bdf::time_der( double dt, const PetscVector& m_vecSol, PetscVector& deriv ) const {
551 ✗ deriv.zeroEntries();
552 ✗ deriv.axpy(m_alpha[0],m_vecSol);
553
554 ✗ switch (m_order) {
555 ✗ case 1:
556 ✗ deriv.axpy(-1.0*m_alpha[1],m_sol_n[0]);
557 ✗ break;
558 ✗ case 2:
559 ✗ deriv.axpy(-1.0*m_alpha[1],m_sol_n[0]);
560 ✗ deriv.axpy(-1.0*m_alpha[2],m_sol_n_1[0]);
561 ✗ break;
562 ✗ case 3:
563 ✗ deriv.axpy(-1.0*m_alpha[1],m_sol_n[0]);
564 ✗ deriv.axpy(-1.0*m_alpha[2],m_sol_n_1[0]);
565 ✗ deriv.axpy(-1.0*m_alpha[3],m_sol_n_2[0]);
566 ✗ break;
567 ✗ default:
568 ✗ FEL_ERROR("unknown BDF order");
569 ✗ break;
570 }
571 ✗ deriv.scale(1./dt);
572 }
573
574 ✗ void Bdf::time_der( double dt, const std::vector<PetscVector>& m_vecSol, std::vector<PetscVector>& deriv ) const {
575 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
576 ✗ deriv[i].zeroEntries();
577 ✗ deriv[i].axpy(m_alpha[0],m_vecSol[i]);
578 }
579 ✗ switch (m_order) {
580 ✗ case 1:
581 ✗ for (felInt i=0 ; i< m_numberOfComp; i++)
582 ✗ deriv[i].axpy(-1.0*m_alpha[1],m_sol_n[i]);
583 ✗ break;
584 ✗ case 2:
585 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
586 ✗ deriv[i].axpy(-1.0*m_alpha[1],m_sol_n[i]);
587 ✗ deriv[i].axpy(-1.0*m_alpha[2],m_sol_n_1[i]);
588 }
589 ✗ break;
590 ✗ case 3:
591 ✗ for (felInt i=0 ; i< m_numberOfComp; i++) {
592 ✗ deriv[i].axpy(-1.0*m_alpha[1],m_sol_n[i]);
593 ✗ deriv[i].axpy(-1.0*m_alpha[2],m_sol_n_1[i]);
594 ✗ deriv[i].axpy(-1.0*m_alpha[3],m_sol_n_2[i]);
595 }
596 ✗ break;
597 ✗ default:
598 ✗ FEL_ERROR("unknown BDF order");
599 ✗ break;
600 }
601 ✗ for (felInt i=0 ; i< m_numberOfComp; i++)
602 ✗ deriv[i].scale(1./dt);
603 }
604 }
605
606