v0.14.0
Loading...
Searching...
No Matches
TsCtx.cpp
Go to the documentation of this file.
1
2
3namespace MoFEM {
4
5TsCtx::TsCtx(MoFEM::Interface &m_field, const std::string &problem_name)
6 : mField(m_field), moab(m_field.get_moab()), problemName(problem_name),
7 bH(MF_EXIST), zeroMatrix(true) {
8 PetscLogEventRegister("LoopTsIFunction", 0, &MOFEM_EVENT_TsCtxIFunction);
9 PetscLogEventRegister("LoopTsIJacobian", 0, &MOFEM_EVENT_TsCtxIJacobian);
10 PetscLogEventRegister("LoopTsRHSFunction", 0, &MOFEM_EVENT_TsCtxRHSFunction);
11 PetscLogEventRegister("LoopTsRHSJacobian", 0, &MOFEM_EVENT_TsCtxRHSJacobian);
12 PetscLogEventRegister("LoopTsMonitor", 0, &MOFEM_EVENT_TsCtxMonitor);
13 PetscLogEventRegister("LoopTsI2Function", 0, &MOFEM_EVENT_TsCtxI2Function);
14 PetscLogEventRegister("LoopTsI2Jacobian", 0, &MOFEM_EVENT_TsCtxI2Jacobian);
15
16 if (!LogManager::checkIfChannelExist("TSWORLD")) {
17 auto core_log = logging::core::get();
18
19 core_log->add_sink(
21 core_log->add_sink(
23 core_log->add_sink(
25
26 LogManager::setLog("TSWORLD");
27 LogManager::setLog("TSSYNC");
28 LogManager::setLog("TSSELF");
29
30 MOFEM_LOG_TAG("TSWORLD", "TS");
31 MOFEM_LOG_TAG("TSSYNC", "TS");
32 MOFEM_LOG_TAG("TSSELF", "TS");
33 }
34}
35
55
56PetscErrorCode TsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec F,
57 void *ctx) {
59 TsCtx *ts_ctx = static_cast<TsCtx *>(ctx);
60 PetscLogEventBegin(ts_ctx->MOFEM_EVENT_TsCtxIFunction, 0, 0, 0, 0);
61 CHKERR VecGhostUpdateBegin(u, INSERT_VALUES, SCATTER_FORWARD);
62 CHKERR VecGhostUpdateEnd(u, INSERT_VALUES, SCATTER_FORWARD);
63 CHKERR VecGhostUpdateBegin(u_t, INSERT_VALUES, SCATTER_FORWARD);
64 CHKERR VecGhostUpdateEnd(u_t, INSERT_VALUES, SCATTER_FORWARD);
65 CHKERR ts_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
66 ts_ctx->problemName, COL, u, INSERT_VALUES, SCATTER_REVERSE);
67
68 auto zero_ghost_vec = [](Vec g) {
70 Vec l;
71 CHKERR VecGhostGetLocalForm(g, &l);
72 double *a;
73 CHKERR VecGetArray(l, &a);
74 int s;
75 CHKERR VecGetLocalSize(l, &s);
76 for (int i = 0; i != s; ++i)
77 a[i] = 0;
78 CHKERR VecRestoreArray(l, &a);
79 CHKERR VecGhostRestoreLocalForm(g, &l);
81 };
82 CHKERR zero_ghost_vec(F);
83
84 ts_ctx->vecAssembleSwitch = boost::movelib::make_unique<bool>(true);
85
86 int step;
87#if PETSC_VERSION_GE(3, 8, 0)
88 CHKERR TSGetStepNumber(ts, &step);
89#else
90 CHKERR TSGetTimeStepNumber(ts, &step);
91#endif
92
93 auto cache_ptr = boost::make_shared<CacheTuple>();
95
96 auto set = [&](auto &fe) {
97 fe.ts = ts;
98 fe.ts_u = u;
99 fe.ts_u_t = u_t;
100 fe.ts_F = F;
101 fe.ts_t = t;
102 fe.ts_step = step;
105 fe.ksp_ctx = KspMethod::CTX_SETFUNCTION;
106 fe.data_ctx = PetscData::CtxSetF | PetscData::CtxSetX |
108
109 CHKERR TSGetSNES(ts, &fe.snes);
110 CHKERR SNESGetKSP(fe.snes, &fe.ksp);
111
112 fe.cacheWeakPtr = cache_ptr;
113 CHK_THROW_MESSAGE(TSGetTimeStep(ts, &fe.ts_dt), "get time step failed");
114 };
115
116 auto unset = [&](auto &fe) {
117 fe.ts_ctx = TSMethod::CTX_TSNONE;
118 fe.snes_ctx = SnesMethod::CTX_SNESNONE;
119 fe.ksp_ctx = KspMethod::CTX_KSPNONE;
120 fe.data_ctx = PetscData::CtxSetNone;
121 };
122
123 // preprocess
124 for (auto &bit : ts_ctx->preProcessIFunction) {
125 bit->vecAssembleSwitch = boost::move(ts_ctx->vecAssembleSwitch);
126 set(*bit);
128 *bit);
129 unset(*bit);
130 ts_ctx->vecAssembleSwitch = boost::move(bit->vecAssembleSwitch);
131 }
132
133 // fe loops
134 for (auto &lit : ts_ctx->loopsIFunction) {
135 lit.second->vecAssembleSwitch = boost::move(ts_ctx->vecAssembleSwitch);
136 set(*lit.second);
138 *(lit.second), nullptr,
139 ts_ctx->bH, cache_ptr);
140 unset(*lit.second);
141 ts_ctx->vecAssembleSwitch = boost::move(lit.second->vecAssembleSwitch);
142 }
143
144 // post process
145 for (auto &bit : ts_ctx->postProcessIFunction) {
146 bit->vecAssembleSwitch = boost::move(ts_ctx->vecAssembleSwitch);
147 set(*bit);
149 *bit);
150 unset(*bit);
151 ts_ctx->vecAssembleSwitch = boost::move(bit->vecAssembleSwitch);
152 }
153
155 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
156 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
157 CHKERR VecAssemblyBegin(F);
158 CHKERR VecAssemblyEnd(F);
159 }
160
161 PetscLogEventEnd(ts_ctx->MOFEM_EVENT_TsCtxIFunction, 0, 0, 0, 0);
163}
164
165PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
166 Mat A, Mat B, void *ctx) {
168
169 TsCtx *ts_ctx = static_cast<TsCtx *>(ctx);
170 PetscLogEventBegin(ts_ctx->MOFEM_EVENT_TsCtxIFunction, 0, 0, 0, 0);
171 CHKERR VecGhostUpdateBegin(u, INSERT_VALUES, SCATTER_FORWARD);
172 CHKERR VecGhostUpdateEnd(u, INSERT_VALUES, SCATTER_FORWARD);
173 CHKERR VecGhostUpdateBegin(u_t, INSERT_VALUES, SCATTER_FORWARD);
174 CHKERR VecGhostUpdateEnd(u_t, INSERT_VALUES, SCATTER_FORWARD);
175 CHKERR ts_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
176 ts_ctx->problemName, COL, u, INSERT_VALUES, SCATTER_REVERSE);
177 if (ts_ctx->zeroMatrix) {
178 CHKERR MatZeroEntries(B);
179 }
180 int step;
181#if PETSC_VERSION_GE(3, 8, 0)
182 CHKERR TSGetStepNumber(ts, &step);
183#else
184 CHKERR TSGetTimeStepNumber(ts, &step);
185#endif
186
188 boost::movelib::make_unique<bool>(ts_ctx->zeroMatrix);
189 auto cache_ptr = boost::make_shared<CacheTuple>();
191
192 auto set = [&](auto &fe) {
193 fe.ts = ts;
194 fe.ts_u = u;
195 fe.ts_u_t = u_t;
196 fe.ts_A = A;
197 fe.ts_B = B;
198 fe.ts_t = t;
199 fe.ts_a = a;
200 fe.ts_step = step;
203 fe.ksp_ctx = KspMethod::CTX_OPERATORS;
206
207 CHKERR TSGetSNES(ts, &fe.snes);
208 CHKERR SNESGetKSP(fe.snes, &fe.ksp);
209
210 fe.cacheWeakPtr = cache_ptr;
211 CHK_THROW_MESSAGE(TSGetTimeStep(ts, &fe.ts_dt), "get time step failed");
212 };
213
214 auto unset = [&](auto &fe) {
215 fe.ts_ctx = TSMethod::CTX_TSNONE;
216 fe.snes_ctx = SnesMethod::CTX_SNESNONE;
217 fe.ksp_ctx = KspMethod::CTX_KSPNONE;
218 fe.data_ctx = PetscData::CtxSetNone;
219 };
220
221 // preproces
222 for (auto &bit : ts_ctx->preProcessIJacobian) {
223 bit->matAssembleSwitch = boost::move(ts_ctx->matAssembleSwitch);
224 set(*bit);
226 *bit);
227 unset(*bit);
228 ts_ctx->matAssembleSwitch = boost::move(bit->matAssembleSwitch);
229 }
230
231 for (auto &lit : ts_ctx->loopsIJacobian) {
232 lit.second->matAssembleSwitch = boost::move(ts_ctx->matAssembleSwitch);
233 set(*lit.second);
235 *(lit.second), nullptr,
236 ts_ctx->bH, cache_ptr);
237 unset(*lit.second);
238 ts_ctx->matAssembleSwitch = boost::move(lit.second->matAssembleSwitch);
239 }
240
241 // post process
242 for (auto &bit : ts_ctx->postProcessIJacobian) {
243 bit->matAssembleSwitch = boost::move(ts_ctx->matAssembleSwitch);
244 set(*bit);
246 *bit);
247 unset(*bit);
248 ts_ctx->matAssembleSwitch = boost::move(bit->matAssembleSwitch);
249 }
250
251 if (*(ts_ctx->matAssembleSwitch)) {
252 CHKERR MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
253 CHKERR MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
254 }
255 PetscLogEventEnd(ts_ctx->MOFEM_EVENT_TsCtxIFunction, 0, 0, 0, 0);
257}
258
259PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u,
260 void *ctx) {
262 TsCtx *ts_ctx = static_cast<TsCtx *>(ctx);
263 PetscLogEventBegin(ts_ctx->MOFEM_EVENT_TsCtxMonitor, 0, 0, 0, 0);
264 CHKERR VecGhostUpdateBegin(u, INSERT_VALUES, SCATTER_FORWARD);
265 CHKERR VecGhostUpdateEnd(u, INSERT_VALUES, SCATTER_FORWARD);
266 CHKERR ts_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
267 ts_ctx->problemName, COL, u, INSERT_VALUES, SCATTER_REVERSE);
268
269 auto cache_ptr = boost::make_shared<CacheTuple>();
271
272 auto set = [&](auto &fe) {
273 fe.ts = ts;
274 fe.ts_u = u;
275 fe.ts_t = t;
276 fe.ts_step = step;
277 fe.ts_F = PETSC_NULL;
279 fe.snes_ctx = SnesMethod::CTX_SNESNONE;
280 fe.ksp_ctx = KspMethod::CTX_KSPNONE;
282
283 CHKERR TSGetSNES(ts, &fe.snes);
284 CHKERR SNESGetKSP(fe.snes, &fe.ksp);
285
286 fe.cacheWeakPtr = cache_ptr;
287 CHK_THROW_MESSAGE(TSGetTimeStep(ts, &fe.ts_dt), "get time step failed");
288 };
289
290 auto unset = [&](auto &fe) {
291 fe.ts_ctx = TSMethod::CTX_TSNONE;
292 fe.data_ctx = PetscData::CtxSetNone;
293 };
294
295 // preprocess
296 for (auto &bit : ts_ctx->preProcessMonitor) {
297 set(*bit);
299 *bit);
300 unset(*bit);
301 }
302
303 for (auto &lit : ts_ctx->loopsMonitor) {
304 set(*lit.second);
306 *(lit.second), nullptr,
307 ts_ctx->bH, cache_ptr);
308 unset(*lit.second);
309 }
310
311 // post process
312 for (auto &bit : ts_ctx->postProcessMonitor) {
313 set(*bit);
315 *bit);
316 unset(*bit);
317 }
318
319 PetscLogEventEnd(ts_ctx->MOFEM_EVENT_TsCtxMonitor, 0, 0, 0, 0);
321}
322
323PetscErrorCode TsSetRHSFunction(TS ts, PetscReal t, Vec u, Vec F, void *ctx) {
325 TsCtx *ts_ctx = static_cast<TsCtx *>(ctx);
326 PetscLogEventBegin(ts_ctx->MOFEM_EVENT_TsCtxRHSFunction, 0, 0, 0, 0);
327 CHKERR VecGhostUpdateBegin(u, INSERT_VALUES, SCATTER_FORWARD);
328 CHKERR VecGhostUpdateEnd(u, INSERT_VALUES, SCATTER_FORWARD);
329 CHKERR ts_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
330 ts_ctx->problemName, COL, u, INSERT_VALUES, SCATTER_REVERSE);
331
332 auto zero_ghost_vec = [](Vec g) {
334 Vec l;
335 CHKERR VecGhostGetLocalForm(g, &l);
336 double *a;
337 CHKERR VecGetArray(l, &a);
338 int s;
339 CHKERR VecGetLocalSize(l, &s);
340 for (int i = 0; i != s; ++i)
341 a[i] = 0;
342 CHKERR VecRestoreArray(l, &a);
343 CHKERR VecGhostRestoreLocalForm(g, &l);
345 };
346 CHKERR zero_ghost_vec(F);
347
348 ts_ctx->vecAssembleSwitch = boost::movelib::make_unique<bool>(true);
349 auto cache_ptr = boost::make_shared<CacheTuple>();
351
352 int step;
353#if PETSC_VERSION_GE(3, 8, 0)
354 CHKERR TSGetStepNumber(ts, &step);
355#else
356 CHKERR TSGetTimeStepNumber(ts, &step);
357#endif
358
359 auto set = [&](auto &fe) {
360 fe.ts_u = u;
361 fe.ts_F = F;
362 fe.ts_t = t;
363 fe.ts = ts;
364 fe.ts_step = step;
367 fe.ksp_ctx = KspMethod::CTX_SETFUNCTION;
368 fe.data_ctx =
370 fe.cacheWeakPtr = cache_ptr;
371 CHK_THROW_MESSAGE(TSGetTimeStep(ts, &fe.ts_dt), "get time step failed");
372 };
373
374 auto unset = [&](auto &fe) {
375 fe.ts_ctx = TSMethod::CTX_TSNONE;
376 fe.snes_ctx = SnesMethod::CTX_SNESNONE;
377 fe.ksp_ctx = KspMethod::CTX_KSPNONE;
378 fe.data_ctx = PetscData::CtxSetNone;
379 };
380
381 for (auto &bit : ts_ctx->preProcessRHSFunction) {
382 bit->vecAssembleSwitch = boost::move(ts_ctx->vecAssembleSwitch);
383 set(*bit);
385 *bit);
386 unset(*bit);
387 ts_ctx->vecAssembleSwitch = boost::move(bit->vecAssembleSwitch);
388 }
389
390 // fe loops
391 for (auto &lit : ts_ctx->loopsRHSFunction) {
392 lit.second->vecAssembleSwitch = boost::move(ts_ctx->vecAssembleSwitch);
393 set(*lit.second);
395 *(lit.second), nullptr,
396 ts_ctx->bH, cache_ptr);
397 unset(*lit.second);
398 ts_ctx->vecAssembleSwitch = boost::move(lit.second->vecAssembleSwitch);
399 }
400
401 // post process
402 for (auto &bit : ts_ctx->postProcessRHSFunction) {
403 bit->vecAssembleSwitch = boost::move(ts_ctx->vecAssembleSwitch);
404 set(*bit);
406 *bit);
407 unset(*bit);
408 ts_ctx->vecAssembleSwitch = boost::move(bit->vecAssembleSwitch);
409 }
410
412 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
413 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
414 CHKERR VecAssemblyBegin(F);
415 CHKERR VecAssemblyEnd(F);
416 }
417
418 PetscLogEventEnd(ts_ctx->MOFEM_EVENT_TsCtxRHSFunction, 0, 0, 0, 0);
420}
421
422PetscErrorCode TsSetRHSJacobian(TS ts, PetscReal t, Vec u, Mat A, Mat B,
423 void *ctx) {
425 TsCtx *ts_ctx = static_cast<TsCtx *>(ctx);
426 PetscLogEventBegin(ts_ctx->MOFEM_EVENT_TsCtxRHSJacobian, 0, 0, 0, 0);
427 CHKERR VecGhostUpdateBegin(u, INSERT_VALUES, SCATTER_FORWARD);
428 CHKERR VecGhostUpdateEnd(u, INSERT_VALUES, SCATTER_FORWARD);
429 CHKERR ts_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
430 ts_ctx->problemName, COL, u, INSERT_VALUES, SCATTER_REVERSE);
431
432 if (ts_ctx->zeroMatrix) {
433 CHKERR MatZeroEntries(B);
434 }
435
437 boost::movelib::make_unique<bool>(ts_ctx->zeroMatrix);
438 auto cache_ptr = boost::make_shared<CacheTuple>();
440
441 int step;
442#if PETSC_VERSION_GE(3, 8, 0)
443 CHKERR TSGetStepNumber(ts, &step);
444#else
445 CHKERR TSGetTimeStepNumber(ts, &step);
446#endif
447
448 auto set = [&](auto &fe) {
449 fe.ts_u = u;
450 fe.ts_A = A;
451 fe.ts_B = B;
452 fe.ts_t = t;
453 fe.ts_step = step;
456 fe.ksp_ctx = KspMethod::CTX_OPERATORS;
459 fe.ts = ts;
460 fe.cacheWeakPtr = cache_ptr;
461 CHK_THROW_MESSAGE(TSGetTimeStep(ts, &fe.ts_dt), "get time step failed");
462 };
463
464 auto unset = [&](auto &fe) {
465 fe.ts_ctx = TSMethod::CTX_TSNONE;
466 fe.snes_ctx = SnesMethod::CTX_SNESNONE;
467 fe.ksp_ctx = KspMethod::CTX_KSPNONE;
468 fe.data_ctx = PetscData::CtxSetNone;
469 };
470
471 // preprocess
472 for (auto &bit : ts_ctx->preProcessRHSJacobian) {
473 bit->matAssembleSwitch = boost::move(ts_ctx->matAssembleSwitch);
474 set(*bit);
476 *bit);
477 unset(*bit);
478 ts_ctx->matAssembleSwitch = boost::move(bit->matAssembleSwitch);
479 }
480
481 // fe loops
482 for (auto &lit : ts_ctx->loopsRHSJacobian) {
483 lit.second->matAssembleSwitch = boost::move(ts_ctx->matAssembleSwitch);
484 set(*lit.second);
486 *(lit.second), nullptr,
487 ts_ctx->bH, cache_ptr);
488 unset(*lit.second);
489 ts_ctx->matAssembleSwitch = boost::move(lit.second->matAssembleSwitch);
490 }
491
492 // post process
493 for (auto &bit : ts_ctx->postProcessRHSJacobian) {
494 bit->matAssembleSwitch = boost::move(ts_ctx->matAssembleSwitch);
495 set(*bit);
497 *bit);
498 unset(*bit);
499 ts_ctx->matAssembleSwitch = boost::move(bit->matAssembleSwitch);
500 }
501
502 if (*(ts_ctx->matAssembleSwitch)) {
503 CHKERR MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
504 CHKERR MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
505 }
506
507 PetscLogEventEnd(ts_ctx->MOFEM_EVENT_TsCtxRHSJacobian, 0, 0, 0, 0);
509}
510
511PetscErrorCode TsSetI2Jacobian(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt,
512 PetscReal a, PetscReal aa, Mat A, Mat B,
513 void *ctx) {
515
516 TsCtx *ts_ctx = static_cast<TsCtx *>(ctx);
517 PetscLogEventBegin(ts_ctx->MOFEM_EVENT_TsCtxI2Function, 0, 0, 0, 0);
518 CHKERR VecGhostUpdateBegin(u, INSERT_VALUES, SCATTER_FORWARD);
519 CHKERR VecGhostUpdateEnd(u, INSERT_VALUES, SCATTER_FORWARD);
520 CHKERR VecGhostUpdateBegin(u_t, INSERT_VALUES, SCATTER_FORWARD);
521 CHKERR VecGhostUpdateEnd(u_t, INSERT_VALUES, SCATTER_FORWARD);
522 CHKERR VecGhostUpdateBegin(u_tt, INSERT_VALUES, SCATTER_FORWARD);
523 CHKERR VecGhostUpdateEnd(u_tt, INSERT_VALUES, SCATTER_FORWARD);
524 CHKERR ts_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
525 ts_ctx->problemName, COL, u, INSERT_VALUES, SCATTER_REVERSE);
526 if (ts_ctx->zeroMatrix) {
527 CHKERR MatZeroEntries(B);
528 }
529 int step;
530#if PETSC_VERSION_GE(3, 8, 0)
531 CHKERR TSGetStepNumber(ts, &step);
532#else
533 CHKERR TSGetTimeStepNumber(ts, &step);
534#endif
535
537 boost::movelib::make_unique<bool>(ts_ctx->zeroMatrix);
538 auto cache_ptr = boost::make_shared<CacheTuple>();
540
541 auto set = [&](auto &fe) {
542 fe.ts_u = u;
543 fe.ts_u_t = u_t;
544 fe.ts_u_tt = u_tt;
545 fe.ts_A = A;
546 fe.ts_B = B;
547 fe.ts_t = t;
548 fe.ts_a = a;
549 fe.ts_aa = aa;
550 fe.ts_step = step;
551
554 fe.ksp_ctx = KspMethod::CTX_OPERATORS;
558
559 CHKERR TSGetSNES(ts, &fe.snes);
560 CHKERR SNESGetKSP(fe.snes, &fe.ksp);
561
562 fe.ts = ts;
563 fe.cacheWeakPtr = cache_ptr;
564 CHK_THROW_MESSAGE(TSGetTimeStep(ts, &fe.ts_dt), "get time step failed");
565 };
566
567 auto unset = [&](auto &fe) {
568 fe.ts_ctx = TSMethod::CTX_TSNONE;
569 fe.snes_ctx = SnesMethod::CTX_SNESNONE;
570 fe.ksp_ctx = KspMethod::CTX_KSPNONE;
571 fe.data_ctx = PetscData::CtxSetNone;
572 };
573
574 // preprocess
575 for (auto &bit : ts_ctx->preProcessIJacobian) {
576 bit->matAssembleSwitch = boost::move(ts_ctx->matAssembleSwitch);
577 set(*bit);
579 *bit);
580 unset(*bit);
581 ts_ctx->matAssembleSwitch = boost::move(bit->matAssembleSwitch);
582 }
583
584 for (auto &lit : ts_ctx->loopsIJacobian) {
585 lit.second->matAssembleSwitch = boost::move(ts_ctx->matAssembleSwitch);
586 set(*lit.second);
588 *(lit.second), nullptr,
589 ts_ctx->bH, cache_ptr);
590 unset(*lit.second);
591 ts_ctx->matAssembleSwitch = boost::move(lit.second->matAssembleSwitch);
592 }
593
594 // post process
595 for (auto &bit : ts_ctx->postProcessIJacobian) {
596 bit->matAssembleSwitch = boost::move(ts_ctx->matAssembleSwitch);
597 set(*bit);
599 *bit);
600 unset(*bit);
601 ts_ctx->matAssembleSwitch = boost::move(bit->matAssembleSwitch);
602 }
603
604 if (*(ts_ctx->matAssembleSwitch)) {
605 CHKERR MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
606 CHKERR MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
607 }
608 PetscLogEventEnd(ts_ctx->MOFEM_EVENT_TsCtxI2Function, 0, 0, 0, 0);
610}
611
612PetscErrorCode TsSetI2Function(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt,
613 Vec F, void *ctx) {
615 TsCtx *ts_ctx = static_cast<TsCtx *>(ctx);
616 PetscLogEventBegin(ts_ctx->MOFEM_EVENT_TsCtxIFunction, 0, 0, 0, 0);
617 CHKERR VecGhostUpdateBegin(u, INSERT_VALUES, SCATTER_FORWARD);
618 CHKERR VecGhostUpdateEnd(u, INSERT_VALUES, SCATTER_FORWARD);
619 CHKERR VecGhostUpdateBegin(u_t, INSERT_VALUES, SCATTER_FORWARD);
620 CHKERR VecGhostUpdateEnd(u_t, INSERT_VALUES, SCATTER_FORWARD);
621 CHKERR VecGhostUpdateBegin(u_tt, INSERT_VALUES, SCATTER_FORWARD);
622 CHKERR VecGhostUpdateEnd(u_tt, INSERT_VALUES, SCATTER_FORWARD);
623 CHKERR ts_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
624 ts_ctx->problemName, COL, u, INSERT_VALUES, SCATTER_REVERSE);
625
626 auto zero_ghost_vec = [](Vec g) {
628 Vec l;
629 CHKERR VecGhostGetLocalForm(g, &l);
630 double *a;
631 CHKERR VecGetArray(l, &a);
632 int s;
633 CHKERR VecGetLocalSize(l, &s);
634 for (int i = 0; i != s; ++i)
635 a[i] = 0;
636 CHKERR VecRestoreArray(l, &a);
637 CHKERR VecGhostRestoreLocalForm(g, &l);
639 };
640 CHKERR zero_ghost_vec(F);
641
642 ts_ctx->vecAssembleSwitch = boost::movelib::make_unique<bool>(true);
643 auto cache_ptr = boost::make_shared<CacheTuple>();
645
646 int step;
647#if PETSC_VERSION_GE(3, 8, 0)
648 CHKERR TSGetStepNumber(ts, &step);
649#else
650 CHKERR TSGetTimeStepNumber(ts, &step);
651#endif
652
653 auto set = [&](auto &fe) {
654 fe.ts_u = u;
655 fe.ts_u_t = u_t;
656 fe.ts_u_tt = u_tt;
657 fe.ts_F = F;
658 fe.ts_t = t;
659 fe.ts_step = step;
662 fe.ksp_ctx = KspMethod::CTX_SETFUNCTION;
663 fe.data_ctx = PetscData::CtxSetF | PetscData::CtxSetX |
666
667 CHKERR TSGetSNES(ts, &fe.snes);
668 CHKERR SNESGetKSP(fe.snes, &fe.ksp);
669
670 fe.ts = ts;
671 fe.cacheWeakPtr = cache_ptr;
672 CHK_THROW_MESSAGE(TSGetTimeStep(ts, &fe.ts_dt), "get time step failed");
673 };
674
675 auto unset = [&](auto &fe) {
676 fe.ts_ctx = TSMethod::CTX_TSNONE;
677 fe.snes_ctx = SnesMethod::CTX_SNESNONE;
678 fe.ksp_ctx = KspMethod::CTX_KSPNONE;
679 fe.data_ctx = PetscData::CtxSetNone;
680 };
681
682 // preprocess
683 for (auto &bit : ts_ctx->preProcessIFunction) {
684 bit->vecAssembleSwitch = boost::move(ts_ctx->vecAssembleSwitch);
685 set(*bit);
687 *bit);
688 unset(*bit);
689 ts_ctx->vecAssembleSwitch = boost::move(bit->vecAssembleSwitch);
690 }
691
692 // fe loops
693 for (auto &lit : ts_ctx->loopsIFunction) {
694 lit.second->vecAssembleSwitch = boost::move(ts_ctx->vecAssembleSwitch);
695 set(*lit.second);
697 *(lit.second), nullptr,
698 ts_ctx->bH, cache_ptr);
699 unset(*lit.second);
700 ts_ctx->vecAssembleSwitch = boost::move(lit.second->vecAssembleSwitch);
701 }
702
703 // post process
704 for (auto &bit : ts_ctx->postProcessIFunction) {
705 bit->vecAssembleSwitch = boost::move(ts_ctx->vecAssembleSwitch);
706 set(*bit);
708 *bit);
709 unset(*bit);
710 ts_ctx->vecAssembleSwitch = boost::move(bit->vecAssembleSwitch);
711 }
712
714 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
715 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
716 CHKERR VecAssemblyBegin(F);
717 CHKERR VecAssemblyEnd(F);
718 }
719
720 PetscLogEventEnd(ts_ctx->MOFEM_EVENT_TsCtxIFunction, 0, 0, 0, 0);
722}
723
725 : alpha(0.75), gamma(0.5), desiredIt(6), offApat(PETSC_FALSE) {
726 CHKERR PetscOptionsGetScalar("", "-ts_mofem_adapt_alpha", &alpha, PETSC_NULL);
727 CHKERR PetscOptionsGetScalar("", "-ts_mofem_adapt_gamma", &gamma, PETSC_NULL);
728 CHKERR PetscOptionsGetInt("", "-ts_mofem_adapt_desired_it", &desiredIt,
729 PETSC_NULL);
730 CHKERR PetscOptionsGetBool("", "-ts_mofem_adapt_off", &offApat, PETSC_NULL);
731
732 MOFEM_LOG("TSWORLD", Sev::inform)
733 << "TS adaptivity: alpha = " << alpha << ", gamma = " << gamma
734 << ", desiredIt = " << desiredIt << ", offAdapt = " << offApat;
735}
736
737PetscErrorCode TSAdaptChooseMoFEM(TSAdapt adapt, TS ts, PetscReal h,
738 PetscInt *next_sc, PetscReal *next_h,
739 PetscBool *accept, PetscReal *wlte,
740 PetscReal *wltea, PetscReal *wlter) {
742
743 auto ts_adapt_mofem = boost::make_shared<TSAdaptMoFEM>();
744
745 *next_sc = 0; /* Reuse the same order scheme */
746 *wlte = -1; /* Weighted local truncation error was not evaluated */
747 *wltea = -1; /* Weighted absolute local truncation error is not used */
748 *wlter = -1; /* Weighted relative local truncation error is not used */
749
750 *accept = PETSC_TRUE;
751 *next_h = h; /* Reuse the old step */
752
753 if (!ts_adapt_mofem->offApat) {
754
755 SNES snes;
756 CHKERR TSGetSNES(ts, &snes);
757
758 SNESConvergedReason reason;
759 CHKERR SNESGetConvergedReason(snes, &reason);
760
761 int it;
762 CHKERR SNESGetIterationNumber(snes, &it);
763
764 if (reason < 0) {
765 h *= ts_adapt_mofem->alpha;
766 *next_h = h;
768 "TSWORLD", Sev::warning,
769 "\tDiverged set step length: it = %d, h = %3.4g set h = %3.4g \n", it,
770 h, *next_h);
771 } else if (reason > 0) {
772 h *= pow((static_cast<double>(ts_adapt_mofem->desiredIt) /
773 static_cast<double>(it + 1)),
774 static_cast<double>(ts_adapt_mofem->gamma));
775 *next_h = PetscClipInterval(h, adapt->dt_min, adapt->dt_max);
777 "TSWORLD", Sev::inform,
778 "\tConverged set step length: it = %d, h = %3.4g set h = %3.4g \n",
779 it, h, *next_h);
780 }
781 }
782
784}
785
786PetscErrorCode TSAdaptResetMoFEM(TSAdapt adapt) {
787 PetscFunctionBegin;
788 PetscFunctionReturn(0);
789}
790
791PetscErrorCode TSAdaptDestroyMoFEM(TSAdapt adapt) {
795}
796
797PetscErrorCode TSAdaptCreateMoFEM(TSAdapt adapt) {
798 PetscFunctionBegin;
799 adapt->ops->choose = TSAdaptChooseMoFEM;
800 adapt->ops->reset = TSAdaptResetMoFEM;
801 adapt->ops->destroy = TSAdaptDestroyMoFEM;
802 PetscFunctionReturn(0);
803}
804
805} // namespace MoFEM
#define MOFEM_LOG_C(channel, severity, format,...)
constexpr double a
@ COL
@ MF_EXIST
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ F
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
MoFEM::TsCtx * ts_ctx
FTensor::Index< 'l', 3 > l
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition TsCtx.cpp:165
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition TsCtx.cpp:259
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode TsSetI2Function(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, Vec F, void *ctx)
Calculation the right hand side for second order PDE in time.
Definition TsCtx.cpp:612
PetscErrorCode TsSetRHSFunction(TS ts, PetscReal t, Vec u, Vec F, void *ctx)
TS solver function.
Definition TsCtx.cpp:323
PetscErrorCode TSAdaptChooseMoFEM(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter)
Definition TsCtx.cpp:737
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
PetscErrorCode TsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
Set IFunction for TS solver.
Definition TsCtx.cpp:56
PetscErrorCode TSAdaptDestroyMoFEM(TSAdapt adapt)
Definition TsCtx.cpp:791
PetscErrorCode TsSetI2Jacobian(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, PetscReal a, PetscReal aa, Mat A, Mat B, void *ctx)
Calculation Jacobian for second order PDE in time.
Definition TsCtx.cpp:511
PetscErrorCode TSAdaptResetMoFEM(TSAdapt adapt)
Definition TsCtx.cpp:786
PetscErrorCode TsSetRHSJacobian(TS ts, PetscReal t, Vec u, Mat A, Mat B, void *ctx)
TS solver function.
Definition TsCtx.cpp:422
PetscErrorCode TSAdaptCreateMoFEM(TSAdapt adapt)
Craete MOFEM adapt.
Definition TsCtx.cpp:797
constexpr AssemblyType A
double h
constexpr double t
plate stiffness
Definition plate.cpp:59
constexpr double g
virtual MoFEMErrorCode cache_problem_entities(const std::string prb_name, CacheTupleWeakPtr cache_ptr)=0
Cache variables.
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
Deprecated interface functions.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
static bool checkIfChannelExist(const std::string channel)
Check if channel exist.
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
static constexpr Switches CtxSetA
static constexpr Switches CtxSetX
static constexpr Switches CtxSetX_TT
static constexpr Switches CtxSetNone
static constexpr Switches CtxSetF
static constexpr Switches CtxSetX_T
static constexpr Switches CtxSetB
static constexpr Switches CtxSetTime
PetscBool offApat
Definition TsCtx.hpp:369
Interface for Time Stepping (TS) solver.
Definition TsCtx.hpp:17
PetscLogEvent MOFEM_EVENT_TsCtxI2Jacobian
Definition TsCtx.hpp:214
FEMethodsSequence loopsIJacobian
Definition TsCtx.hpp:29
bool zeroMatrix
Definition TsCtx.hpp:45
MoFEMTypes bH
If set to MF_EXIST check if element exist.
Definition TsCtx.hpp:23
PetscLogEvent MOFEM_EVENT_TsCtxRHSJacobian
Definition TsCtx.hpp:209
FEMethodsSequence loopsRHSFunction
Definition TsCtx.hpp:33
PetscLogEvent MOFEM_EVENT_TsCtxMonitor
Definition TsCtx.hpp:212
BasicMethodsSequence postProcessRHSFunction
Definition TsCtx.hpp:43
BasicMethodsSequence preProcessMonitor
Definition TsCtx.hpp:38
PetscLogEvent MOFEM_EVENT_TsCtxRHSFunction
Definition TsCtx.hpp:208
BasicMethodsSequence preProcessIJacobian
Definition TsCtx.hpp:34
BasicMethodsSequence postProcessIFunction
Definition TsCtx.hpp:37
FEMethodsSequence loopsIFunction
Definition TsCtx.hpp:30
FEMethodsSequence loopsRHSJacobian
Definition TsCtx.hpp:32
BasicMethodsSequence preProcessIFunction
Definition TsCtx.hpp:36
FEMethodsSequence loopsMonitor
Definition TsCtx.hpp:31
BasicMethodsSequence postProcessMonitor
Definition TsCtx.hpp:39
BasicMethodsSequence preProcessRHSJacobian
Definition TsCtx.hpp:40
boost::movelib::unique_ptr< bool > vecAssembleSwitch
Definition TsCtx.hpp:216
PetscLogEvent MOFEM_EVENT_TsCtxIJacobian
Definition TsCtx.hpp:211
boost::movelib::unique_ptr< bool > matAssembleSwitch
Definition TsCtx.hpp:217
MoFEMErrorCode clearLoops()
Clear loops.
Definition TsCtx.cpp:36
std::string problemName
Definition TsCtx.hpp:22
TsCtx(MoFEM::Interface &m_field, const std::string &problem_name)
Definition TsCtx.cpp:5
MoFEM::Interface & mField
Definition TsCtx.hpp:19
BasicMethodsSequence preProcessRHSFunction
Definition TsCtx.hpp:41
BasicMethodsSequence postProcessIJacobian
Definition TsCtx.hpp:35
PetscLogEvent MOFEM_EVENT_TsCtxIFunction
Definition TsCtx.hpp:210
BasicMethodsSequence postProcessRHSJacobian
Definition TsCtx.hpp:42
PetscLogEvent MOFEM_EVENT_TsCtxI2Function
Definition TsCtx.hpp:213
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.