v0.15.0
Loading...
Searching...
No Matches
SnesCtx.cpp
Go to the documentation of this file.
1/**
2 * @brief Context for PETSc SNES, i.e. nonlinear solver
3 * @author Anonymous contributor under MIT license
4 *
5 */
6
7namespace MoFEM {
8
10
15
16 SnesCtxImpl(Interface &m_field, std::string problem_name);
17
18 virtual ~SnesCtxImpl() = default;
19
21
23
25
27
29
33
35
39
43
44 MoFEMErrorCode copyLoops(const SnesCtxImpl &snes_ctx);
45
47
49
51
52 friend PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx);
53
54 friend PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx);
55
57 MatAssemblyType type);
59
61 MatAssemblyType type);
62
64
65 friend MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its,
66 PetscReal fgnorm,
67 SnesCtx *snes_ctx);
68
69 friend MoFEMErrorCode MoFEMSNESMonitorEnergy(SNES snes, PetscInt its,
70 PetscReal fgnorm,
71 SnesCtx *snes_ctx);
72
73 friend PetscErrorCode SnesLoadTangent(SNES snes, Vec u, Vec F, void *ctx);
74
75protected:
76 PetscLogEvent MOFEM_EVENT_SnesRhs; ///< Log events to assemble residual
77 PetscLogEvent MOFEM_EVENT_SnesMat; ///< Log events to assemble tangent matrix
78 PetscLogEvent MOFEM_EVENT_SnesLoadTangent; ///< Log events to assemble tangent
79 ///< of load vector
80
81 MoFEM::Interface &mField; ///< database Interface
82 moab::Interface &moab; ///< moab Interface
83
84 std::string problemName; ///< problem name
85
86 FEMethodsSequence loopsOperator; ///< Sequence of finite elements instances
87 ///< assembling tangent matrix
88 FEMethodsSequence loopsRhs; ///< Sequence of finite elements instances
89 ///< assembling residual vector
91 loopsLoadTangent; ///< Sequence of finite elements instances assembling
92 ///< tangent of load vector in arc length control
93
94 BasicMethodsSequence preProcessOperator; ///< Sequence of methods run before
95 ///< tangent matrix is assembled
96 BasicMethodsSequence postProcessOperator; ///< Sequence of methods run after
97 ///< tangent matrix is assembled
99 preProcessLoadTangent; ///< Sequence of methods run before
100 ///< tangent of load vector in arc
101 ///< length control is assembled
103 preProcessRhs; ///< Sequence of methods run before residual is assembled
105 postProcessRhs; ///< Sequence of methods run after residual is assembled
107 postProcessLoadTangent; ///< Sequence of methods run after tangent of load
108 ///< vector in arc length control is assembled
109
111 MF_EXIST; ///< If set to MF_EXIST check if element exist, default MF_EXIST
112 bool zeroPreCondMatrixB = true; ///< If true zero matrix, otherwise user need
113 ///< to do it, default true
114 bool vErify = false; ///< If true verify vector
115 MatAssemblyType typeOfAssembly =
116 MAT_FINAL_ASSEMBLY; ///< type of assembly at the end
117
118 HookFunction loadTangentHook; ///< Hook function called at the end of load
119 ///< tangent evaluation
120
121 HookFunction rhsHook; ///< Hook function called at the end of residual
122 ///< evaluation
123};
124
125SnesCtx::SnesCtx(Interface &m_field, std::string problem_name)
126 : snesCtxImpl(boost::movelib::make_unique<SnesCtx::SnesCtxImpl>(
127 m_field, problem_name)) {}
128
129SnesCtx::~SnesCtx() = default;
130
132 return snesCtxImpl->getSetOperators();
133}
134
136 return snesCtxImpl->getComputeRhs();
137}
138
140 return snesCtxImpl->getPreProcComputeRhs();
141}
142
144 return snesCtxImpl->getPostProcComputeRhs();
145}
146
148 return snesCtxImpl->getPreProcSetOperators();
149}
150
152 return snesCtxImpl->getPostProcSetOperators();
153}
154
156 return snesCtxImpl->getLoadTangent();
157}
158
160 return snesCtxImpl->getPreProcLoadTangent();
161}
162
164 return snesCtxImpl->getPostProcLoadTangent();
165}
166
168 return snesCtxImpl->copyLoops(*snes_ctx.snesCtxImpl);
169}
170
172
174 return snesCtxImpl->getLoadTangentHook();
175}
176
178 return snesCtxImpl->getRhsHook();
179}
180
181SnesCtx::SnesCtxImpl::SnesCtxImpl(Interface &m_field, std::string problem_name)
182 : mField(m_field), moab(m_field.get_moab()), problemName(problem_name) {
183 PetscLogEventRegister("LoopSNESRhs", 0, &MOFEM_EVENT_SnesRhs);
184 PetscLogEventRegister("LoopSNESMat", 0, &MOFEM_EVENT_SnesMat);
185 PetscLogEventRegister("LoopSNESLoadTangent", 0, &MOFEM_EVENT_SnesLoadTangent);
186 if (!LogManager::checkIfChannelExist("SNES_WORLD")) {
187 auto core_log = logging::core::get();
188 core_log->add_sink(
190 LogManager::setLog("SNES_WORLD");
191 MOFEM_LOG_TAG("SNES_WORLD", "SNES");
192 }
193}
194
197 loopsOperator = snes_ctx.loopsOperator;
198 loopsRhs = snes_ctx.loopsRhs;
199 preProcessOperator = snes_ctx.preProcessOperator;
200 postProcessOperator = snes_ctx.postProcessOperator;
201 preProcessRhs = snes_ctx.preProcessRhs;
202 postProcessRhs = snes_ctx.postProcessRhs;
203 loopsLoadTangent = snes_ctx.loopsLoadTangent;
204 preProcessLoadTangent = snes_ctx.preProcessLoadTangent;
205 postProcessLoadTangent = snes_ctx.postProcessLoadTangent;
207}
208
211 loopsOperator.clear();
212 loopsRhs.clear();
213 preProcessOperator.clear();
214 postProcessOperator.clear();
215 preProcessRhs.clear();
216 postProcessRhs.clear();
217 loopsLoadTangent.clear();
218 preProcessLoadTangent.clear();
219 postProcessLoadTangent.clear();
221}
222
223PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx) {
224 auto *snes_ctx = static_cast<SnesCtx *>(ctx)->snesCtxImpl.get();
226 PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesRhs, 0, 0, 0, 0);
227 Vec dx;
228 CHKERR SNESGetSolutionUpdate(snes, &dx);
229 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
230 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
231 if (snes_ctx->vErify) {
232 // Verify finite elements, check for not a number
233 CHKERR VecAssemblyBegin(f);
234 CHKERR VecAssemblyEnd(f);
235 MPI_Comm comm = PetscObjectComm((PetscObject)f);
236 PetscSynchronizedPrintf(comm, "SNES Verify x\n");
237 const Problem *prb_ptr;
238 CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
239 CHKERR snes_ctx->mField.getInterface<Tools>()->checkVectorForNotANumber(
240 prb_ptr, COL, x);
241 }
242 CHKERR snes_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
243 snes_ctx->problemName, COL, x, INSERT_VALUES, SCATTER_REVERSE);
244
245 auto zero_ghost_vec = [](Vec g) {
247 Vec l;
248 CHKERR VecGhostGetLocalForm(g, &l);
249 double *a;
250 CHKERR VecGetArray(l, &a);
251 int s;
252 CHKERR VecGetLocalSize(l, &s);
253 for (int i = 0; i != s; ++i)
254 a[i] = 0;
255 CHKERR VecRestoreArray(l, &a);
256 CHKERR VecGhostRestoreLocalForm(g, &l);
258 };
259 CHKERR zero_ghost_vec(f);
260
261 auto vec_assemble_switch = boost::movelib::make_unique<bool>(true);
262 auto cache_ptr = boost::make_shared<CacheTuple>();
263 CHKERR snes_ctx->mField.cache_problem_entities(snes_ctx->problemName,
264 cache_ptr);
265
266 auto set = [&](auto &fe) {
268 fe.snes = snes;
269 fe.snes_x = x;
270 fe.snes_dx = dx;
271 fe.snes_f = f;
273 fe.ksp_ctx = KspMethod::CTX_SETFUNCTION;
275
276 CHKERR SNESGetKSP(snes, &fe.ksp);
277
278 // If SNES is of type SNESNEWTONAL then get load parameter, will act as a
279 // psued time
280 PetscBool same;
281 CHKERR PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONAL, &same);
282 if (same) {
283 CHKERR SNESNewtonALGetLoadParameter(snes, &fe.ts_t);
284 fe.data_ctx |= PetscData::CtxSetTime;
285 }
286
287 fe.cacheWeakPtr = cache_ptr;
289 };
290
291 auto unset = [&](auto &fe) {
292 fe.snes_ctx = SnesMethod::CTX_SNESNONE;
293 fe.ksp_ctx = KspMethod::CTX_KSPNONE;
294 fe.data_ctx = PetscData::CtxSetNone;
295 };
296
297 for (auto &bit : snes_ctx->preProcessRhs) {
298 bit->vecAssembleSwitch = boost::move(vec_assemble_switch);
299 CHKERR set(*bit);
300 CHKERR snes_ctx->mField.problem_basic_method_preProcess(
301 snes_ctx->problemName, *bit);
302 unset(*bit);
303 vec_assemble_switch = boost::move(bit->vecAssembleSwitch);
304 }
305
306 for (auto &lit : snes_ctx->loopsRhs) {
307 lit.second->vecAssembleSwitch = boost::move(vec_assemble_switch);
308 CHKERR set(*lit.second);
309 CHKERR snes_ctx->mField.loop_finite_elements(
310 snes_ctx->problemName, lit.first, *lit.second, nullptr, snes_ctx->bH,
311 cache_ptr);
312 unset(*lit.second);
313 if (snes_ctx->vErify) {
314 // Verify finite elements, check for not a number
315 CHKERR VecAssemblyBegin(f);
316 CHKERR VecAssemblyEnd(f);
317 MPI_Comm comm = PetscObjectComm((PetscObject)f);
318 PetscSynchronizedPrintf(comm, "SNES Verify f FE < %s >\n",
319 lit.first.c_str());
320 const Problem *prb_ptr;
321 CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
322 CHKERR snes_ctx->mField.getInterface<Tools>()->checkVectorForNotANumber(
323 prb_ptr, ROW, f);
324 }
325
326 vec_assemble_switch = boost::move(lit.second->vecAssembleSwitch);
327 }
328
329 for (auto &bit : snes_ctx->postProcessRhs) {
330 bit->vecAssembleSwitch = boost::move(vec_assemble_switch);
331 CHKERR set(*bit);
332 CHKERR snes_ctx->mField.problem_basic_method_postProcess(
333 snes_ctx->problemName, *bit);
334 unset(*bit);
335 vec_assemble_switch = boost::move(bit->vecAssembleSwitch);
336 }
337
338 if (*(vec_assemble_switch)) {
339 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
340 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
341 CHKERR VecAssemblyBegin(f);
342 CHKERR VecAssemblyEnd(f);
343 }
344
345 if (snes_ctx->rhsHook) {
346 CHKERR snes_ctx->rhsHook(snes, x, f, cache_ptr, ctx);
347 }
348
349 PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesRhs, 0, 0, 0, 0);
351}
352
353PetscErrorCode SnesLoadTangent(SNES snes, Vec x, Vec f, void *ctx) {
354 auto *snes_ctx = static_cast<SnesCtx *>(ctx)->snesCtxImpl.get();
356 PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesLoadTangent, 0, 0, 0, 0);
357 Vec dx;
358 CHKERR SNESGetSolutionUpdate(snes, &dx);
359 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
360 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
361 if (snes_ctx->vErify) {
362 // Verify finite elements, check for not a number
363 CHKERR VecAssemblyBegin(f);
364 CHKERR VecAssemblyEnd(f);
365 MPI_Comm comm = PetscObjectComm((PetscObject)f);
366 PetscSynchronizedPrintf(comm, "SNES Verify x\n");
367 const Problem *prb_ptr;
368 CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
369 CHKERR snes_ctx->mField.getInterface<Tools>()->checkVectorForNotANumber(
370 prb_ptr, COL, x);
371 }
372 CHKERR snes_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
373 snes_ctx->problemName, COL, x, INSERT_VALUES, SCATTER_REVERSE);
374
375 auto zero_ghost_vec = [](Vec g) {
377 Vec l;
378 CHKERR VecGhostGetLocalForm(g, &l);
379 double *a;
380 CHKERR VecGetArray(l, &a);
381 int s;
382 CHKERR VecGetLocalSize(l, &s);
383 for (int i = 0; i != s; ++i)
384 a[i] = 0;
385 CHKERR VecRestoreArray(l, &a);
386 CHKERR VecGhostRestoreLocalForm(g, &l);
388 };
389 CHKERR zero_ghost_vec(f);
390
391 auto vec_assemble_switch = boost::movelib::make_unique<bool>(true);
392 auto cache_ptr = boost::make_shared<CacheTuple>();
393 CHKERR snes_ctx->mField.cache_problem_entities(snes_ctx->problemName,
394 cache_ptr);
395
396 auto set = [&](auto &fe) {
398 fe.snes = snes;
399 fe.snes_x = x;
400 fe.snes_dx = dx;
401 fe.snes_f = f;
403 fe.ksp_ctx = KspMethod::CTX_SETFUNCTION;
405
406 CHKERR SNESGetKSP(snes, &fe.ksp);
407
408 // If SNES is of type SNESNEWTONAL then get load parameter, will act as a
409 // psued time
410 PetscBool same;
411 CHKERR PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONAL, &same);
412 if (same) {
413 CHKERR SNESNewtonALGetLoadParameter(snes, &fe.ts_t);
414 fe.data_ctx |= PetscData::CtxSetTime;
415 }
416
417 fe.cacheWeakPtr = cache_ptr;
419 };
420
421 auto unset = [&](auto &fe) {
422 fe.snes_ctx = SnesMethod::CTX_SNESNONE;
423 fe.ksp_ctx = KspMethod::CTX_KSPNONE;
424 fe.data_ctx = PetscData::CtxSetNone;
425 };
426
427 for (auto &bit : snes_ctx->preProcessLoadTangent) {
428 bit->vecAssembleSwitch = boost::move(vec_assemble_switch);
429 CHKERR set(*bit);
430 CHKERR snes_ctx->mField.problem_basic_method_preProcess(
431 snes_ctx->problemName, *bit);
432 unset(*bit);
433 vec_assemble_switch = boost::move(bit->vecAssembleSwitch);
434 }
435
436 for (auto &lit : snes_ctx->loopsLoadTangent) {
437 lit.second->vecAssembleSwitch = boost::move(vec_assemble_switch);
438 CHKERR set(*lit.second);
439 CHKERR snes_ctx->mField.loop_finite_elements(
440 snes_ctx->problemName, lit.first, *lit.second, nullptr, snes_ctx->bH,
441 cache_ptr);
442 unset(*lit.second);
443 if (snes_ctx->vErify) {
444 // Verify finite elements, check for not a number
445 CHKERR VecAssemblyBegin(f);
446 CHKERR VecAssemblyEnd(f);
447 MPI_Comm comm = PetscObjectComm((PetscObject)f);
448 PetscSynchronizedPrintf(comm, "SNES Verify f FE < %s >\n",
449 lit.first.c_str());
450 const Problem *prb_ptr;
451 CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
452 CHKERR snes_ctx->mField.getInterface<Tools>()->checkVectorForNotANumber(
453 prb_ptr, ROW, f);
454 }
455
456 vec_assemble_switch = boost::move(lit.second->vecAssembleSwitch);
457 }
458
459 for (auto &bit : snes_ctx->postProcessLoadTangent) {
460 bit->vecAssembleSwitch = boost::move(vec_assemble_switch);
461 CHKERR set(*bit);
462 CHKERR snes_ctx->mField.problem_basic_method_postProcess(
463 snes_ctx->problemName, *bit);
464 unset(*bit);
465 vec_assemble_switch = boost::move(bit->vecAssembleSwitch);
466 }
467
468 if (*(vec_assemble_switch)) {
469 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
470 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
471 CHKERR VecAssemblyBegin(f);
472 CHKERR VecAssemblyEnd(f);
473 }
474
475 if (snes_ctx->loadTangentHook) {
476 CHKERR snes_ctx->loadTangentHook(snes, x, f, cache_ptr, ctx);
477 }
478
479 PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesLoadTangent, 0, 0, 0, 0);
481}
482
483PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx) {
484 auto *snes_ctx = static_cast<SnesCtx *>(ctx)->snesCtxImpl.get();
486 PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesMat, 0, 0, 0, 0);
487
488 auto mat_assemble_switch = boost::movelib::make_unique<bool>(true);
489 auto cache_ptr = boost::make_shared<CacheTuple>();
490 CHKERR snes_ctx->mField.cache_problem_entities(snes_ctx->problemName,
491 cache_ptr);
492 Vec dx;
493 CHKERR SNESGetSolutionUpdate(snes, &dx);
494
495 auto set = [&](auto &fe) {
497 fe.snes = snes;
498 fe.snes_x = x;
499 fe.snes_dx = dx;
500 fe.snes_A = A;
501 fe.snes_B = B;
503 fe.ksp_ctx = KspMethod::CTX_OPERATORS;
506
507 CHKERR SNESGetKSP(snes, &fe.ksp);
508
509 fe.cacheWeakPtr = cache_ptr;
511 };
512
513 auto unset = [&](auto &fe) {
514 fe.snes_ctx = SnesMethod::CTX_SNESNONE;
515 fe.ksp_ctx = KspMethod::CTX_KSPNONE;
516 fe.data_ctx = PetscData::CtxSetNone;
517 };
518
519 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
520 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
521 CHKERR snes_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
522 snes_ctx->problemName, COL, x, INSERT_VALUES, SCATTER_REVERSE);
523
524 if (*mat_assemble_switch) {
525 if (snes_ctx->zeroPreCondMatrixB)
526 CHKERR MatZeroEntries(B);
527 }
528 for (auto &bit : snes_ctx->preProcessOperator) {
529 bit->matAssembleSwitch = boost::move(mat_assemble_switch);
530 CHKERR set(*bit);
531 CHKERR snes_ctx->mField.problem_basic_method_preProcess(
532 snes_ctx->problemName, *bit);
533 unset(*bit);
534 mat_assemble_switch = boost::move(bit->matAssembleSwitch);
535 }
536
537 for (auto &lit : snes_ctx->loopsOperator) {
538 lit.second->matAssembleSwitch = boost::move(mat_assemble_switch);
539 CHKERR set(*lit.second);
540 CHKERR snes_ctx->mField.loop_finite_elements(
541 snes_ctx->problemName, lit.first, *(lit.second), nullptr, snes_ctx->bH,
542 cache_ptr);
543 unset(*lit.second);
544 mat_assemble_switch = boost::move(lit.second->matAssembleSwitch);
545 }
546
547 for (auto &bit : snes_ctx->postProcessOperator) {
548 bit->matAssembleSwitch = boost::move(mat_assemble_switch);
549 CHKERR set(*bit);
550 CHKERR snes_ctx->mField.problem_basic_method_postProcess(
551 snes_ctx->problemName, *bit);
552 unset(*bit);
553 mat_assemble_switch = boost::move(bit->matAssembleSwitch);
554 }
555
556 if (*mat_assemble_switch) {
557 CHKERR MatAssemblyBegin(B, snes_ctx->typeOfAssembly);
558 CHKERR MatAssemblyEnd(B, snes_ctx->typeOfAssembly);
559 }
560 PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesMat, 0, 0, 0, 0);
562}
563
564MoFEMErrorCode SnesMoFEMSetAssemblyType(SNES snes, MatAssemblyType type) {
565
566 auto get_snes_ctx = [](SNES snes) {
567 SnesCtx *snes_ctx;
568 CHK_THROW_MESSAGE(SNESGetApplicationContext(snes, &snes_ctx),
569 "Cannot get SNES context");
570 return snes_ctx->snesCtxImpl.get();
571 };
572
574 get_snes_ctx(snes)->typeOfAssembly = type;
576}
577
579
580 auto get_snes_ctx = [](SNES snes) {
581 SnesCtx *snes_ctx;
582 CHK_THROW_MESSAGE(SNESGetApplicationContext(snes, &snes_ctx),
583 "Cannot get SNES context");
584 return snes_ctx->snesCtxImpl.get();
585 };
586
588 get_snes_ctx(snes)->bH = bh;
590}
591
592MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm,
593 SnesCtx *ctx) {
595 auto snes_ctx = ctx->snesCtxImpl.get();
596 auto &m_field = snes_ctx->mField;
597 auto problem_ptr = m_field.get_problem(snes_ctx->problemName);
598 auto fields_ptr = m_field.get_fields();
599 auto dofs = problem_ptr->numeredRowDofsPtr;
600
601 std::vector<double> lnorms(fields_ptr->size(), 0),
602 norms(fields_ptr->size(), 0);
603
604 Vec res;
605 CHKERR SNESGetFunction(snes, &res, NULL, NULL);
606
607 const double *r;
608 CHKERR VecGetArrayRead(res, &r);
609 {
610 int f = 0;
611 for (auto fi : *fields_ptr) {
612 const auto lo_uid = FieldEntity::getLoBitNumberUId(fi->bitNumber);
613 const auto hi_uid = FieldEntity::getHiBitNumberUId(fi->bitNumber);
614 const auto hi = dofs->get<Unique_mi_tag>().upper_bound(hi_uid);
615 for (auto lo = dofs->get<Unique_mi_tag>().lower_bound(lo_uid); lo != hi;
616 ++lo) {
617 const DofIdx loc_idx = (*lo)->getPetscLocalDofIdx();
618 if (loc_idx >= 0 && loc_idx < problem_ptr->nbLocDofsRow) {
619 lnorms[f] += PetscRealPart(PetscSqr(r[loc_idx]));
620 }
621 }
622 ++f;
623 }
624 }
625 CHKERR VecRestoreArrayRead(res, &r);
626
627 MPIU_Allreduce(&*lnorms.begin(), &*norms.begin(), lnorms.size(), MPIU_REAL,
628 MPIU_SUM, PetscObjectComm((PetscObject)snes));
629
630 {
631 MOFEM_LOG("SNES_WORLD", Sev::inform)
632 << its << " Function norm " << boost::format("%14.12e") % (double)fgnorm
633 << " [";
634 int f = 0;
635 for (auto fi : *fields_ptr) {
636 MOFEM_LOG("SNES_WORLD", Sev::inform)
637 << its << "\t Field norm "
638 << boost::format("%14.12e") % (double)PetscSqrtReal(norms[f]) << " "
639 << fi->getName();
640 ++f;
641 }
642 MOFEM_LOG("SNES_WORLD", Sev::inform) << its << "]";
643 }
644
646}
647
648MoFEMErrorCode MoFEMSNESMonitorEnergy(SNES snes, PetscInt its, PetscReal fgnorm,
649 SnesCtx *ctx) {
651 auto snes_ctx = ctx->snesCtxImpl.get();
652 auto &m_field = snes_ctx->mField;
653 auto problem_ptr = m_field.get_problem(snes_ctx->problemName);
654 auto fields_ptr = m_field.get_fields();
655 auto dofs = problem_ptr->numeredRowDofsPtr;
656
657 std::vector<double> lnorms(3 * fields_ptr->size(), 0),
658 norms(3 * fields_ptr->size(), 0);
659
660 Vec x_update, res;
661 if (its == 0)
662 CHKERR SNESGetSolution(snes, &x_update);
663 else
664 CHKERR SNESGetSolutionUpdate(snes, &x_update);
665 CHKERR SNESGetFunction(snes, &res, PETSC_NULLPTR, PETSC_NULLPTR);
666
667 const double *x, *r;
668 ;
669 CHKERR VecGetArrayRead(res, &r);
670 CHKERR VecGetArrayRead(x_update, &x);
671
672 {
673 int f = 0;
674 for (auto fi : *fields_ptr) {
675 const auto lo_uid = FieldEntity::getLoBitNumberUId(fi->bitNumber);
676 const auto hi_uid = FieldEntity::getHiBitNumberUId(fi->bitNumber);
677 const auto hi = dofs->get<Unique_mi_tag>().upper_bound(hi_uid);
678 for (auto lo = dofs->get<Unique_mi_tag>().lower_bound(lo_uid); lo != hi;
679 ++lo) {
680 const DofIdx loc_idx = (*lo)->getPetscLocalDofIdx();
681 if (loc_idx >= 0 && loc_idx < problem_ptr->nbLocDofsRow) {
682 lnorms[3 * f] += std::abs(x[loc_idx] * r[loc_idx]);
683 lnorms[3 * f + 1] += r[loc_idx] * r[loc_idx];
684 lnorms[3 * f + 2] += x[loc_idx] * x[loc_idx];
685 }
686 }
687 ++f;
688 }
689 }
690 CHKERR VecRestoreArrayRead(res, &r);
691 CHKERR VecRestoreArrayRead(x_update, &x);
692
693 MPIU_Allreduce(&*lnorms.begin(), &*norms.begin(), lnorms.size(), MPIU_REAL,
694 MPIU_SUM, PetscObjectComm((PetscObject)snes));
695
696 double energy_norm = 0, residual_norm = 0, increment_norm = 0;
697 for (int f = 0; f != fields_ptr->size(); ++f) {
698 energy_norm += norms[3 * f];
699 residual_norm += norms[3 * f + 1];
700 increment_norm += norms[3 * f + 2];
701 }
702 energy_norm = std::sqrt(energy_norm);
703 residual_norm = std::sqrt(residual_norm);
704 increment_norm = std::sqrt(increment_norm);
705
706 {
707 MOFEM_LOG("SNES_WORLD", Sev::inform)
708 << its << " Function norm Enrgy/Residual/Increment"
709 << boost::format(" %14.12e ") % energy_norm
710 << boost::format(" %14.12e ") % residual_norm
711 << boost::format(" %14.12e ") % increment_norm << "[";
712 int f = 0;
713 for (auto fi : *fields_ptr) {
714 MOFEM_LOG("SNES_WORLD", Sev::inform)
715 << its << "\t Energy/Residual/Increment norm "
716 << boost::format("%14.12e") % std::sqrt(norms[3 * f]) << "\t"
717 << boost::format("%14.12e") % std::sqrt(norms[3 * f + 1]) << " "
718 << boost::format("%14.12e") % std::sqrt(norms[3 * f + 2]) << "\t "
719 << fi->getName();
720 ++f;
721 }
722 MOFEM_LOG("SNES_WORLD", Sev::inform) << its << " ]";
723 }
724
726}
727
728} // namespace MoFEM
constexpr double a
@ COL
@ ROW
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
@ 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.
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
int DofIdx
Index of DOF.
Definition Types.hpp:18
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
std::deque< BasicMethodPtr > BasicMethodsSequence
Definition AuxPETSc.hpp:57
std::deque< PairNameFEMethodPtr > FEMethodsSequence
Definition AuxPETSc.hpp:56
constexpr AssemblyType A
constexpr double g
Deprecated interface functions.
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
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 bool checkIfChannelExist(const std::string channel)
Check if channel exist.
static constexpr Switches CtxSetA
static constexpr Switches CtxSetX
static constexpr Switches CtxSetNone
static constexpr Switches CtxSetF
static constexpr Switches CtxSetDX
static constexpr Switches CtxSetB
static constexpr Switches CtxSetTime
keeps basic data about problem
BasicMethodsSequence preProcessLoadTangent
Definition SnesCtx.cpp:99
FEMethodsSequence & getComputeRhs()
Definition SnesCtx.cpp:22
MatAssemblyType typeOfAssembly
type of assembly at the end
Definition SnesCtx.cpp:115
friend PetscErrorCode SnesLoadTangent(SNES snes, Vec u, Vec F, void *ctx)
This function calls finite element pipeline to compute tangent of of load vector in arc length contro...
Definition SnesCtx.cpp:353
FEMethodsSequence & getLoadTangent()
Definition SnesCtx.cpp:34
BasicMethodsSequence & getPreProcLoadTangent()
Definition SnesCtx.cpp:36
friend MoFEMErrorCode SnesMoFEMSetBehavior(SNES snes, MoFEMTypes bh)
Set behavior if finite element in sequence does not exist.
Definition SnesCtx.cpp:578
PetscLogEvent MOFEM_EVENT_SnesMat
Log events to assemble tangent matrix.
Definition SnesCtx.cpp:77
BasicMethodsSequence preProcessOperator
Definition SnesCtx.cpp:94
friend MoFEMErrorCode MoFEMSNESMonitorEnergy(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition SnesCtx.cpp:648
FEMethodsSequence loopsRhs
Definition SnesCtx.cpp:88
BasicMethodsSequence & getPostProcComputeRhs()
Definition SnesCtx.cpp:26
moab::Interface & moab
moab Interface
Definition SnesCtx.cpp:82
BasicMethodsSequence postProcessLoadTangent
Definition SnesCtx.cpp:107
MoFEM::BasicMethodsSequence BasicMethodsSequence
Definition SnesCtx.cpp:13
std::string problemName
problem name
Definition SnesCtx.cpp:84
friend MoFEMErrorCode SNESMoFEMSetAssmblyType(SNES snes, MatAssemblyType type)
virtual ~SnesCtxImpl()=default
MoFEMErrorCode copyLoops(const SnesCtxImpl &snes_ctx)
Definition SnesCtx.cpp:195
friend PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition SnesCtx.cpp:223
friend PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition SnesCtx.cpp:483
HookFunction & getRhsHook()
Definition SnesCtx.cpp:50
FEMethodsSequence loopsOperator
Definition SnesCtx.cpp:86
SnesCtx::HookFunction HookFunction
Definition SnesCtx.cpp:14
BasicMethodsSequence & getPreProcComputeRhs()
Definition SnesCtx.cpp:24
FEMethodsSequence loopsLoadTangent
Definition SnesCtx.cpp:91
FEMethodsSequence & getSetOperators()
Definition SnesCtx.cpp:20
MoFEM::Interface & mField
database Interface
Definition SnesCtx.cpp:81
BasicMethodsSequence postProcessRhs
Sequence of methods run after residual is assembled.
Definition SnesCtx.cpp:105
BasicMethodsSequence & getPostProcSetOperators()
Definition SnesCtx.cpp:30
HookFunction & getLoadTangentHook()
Definition SnesCtx.cpp:48
MoFEMTypes bH
If set to MF_EXIST check if element exist, default MF_EXIST.
Definition SnesCtx.cpp:110
PetscLogEvent MOFEM_EVENT_SnesLoadTangent
Definition SnesCtx.cpp:78
SnesCtxImpl(Interface &m_field, std::string problem_name)
Definition SnesCtx.cpp:181
BasicMethodsSequence preProcessRhs
Sequence of methods run before residual is assembled.
Definition SnesCtx.cpp:103
MoFEMErrorCode clearLoops()
Definition SnesCtx.cpp:209
BasicMethodsSequence & getPostProcLoadTangent()
Definition SnesCtx.cpp:40
BasicMethodsSequence postProcessOperator
Definition SnesCtx.cpp:96
MoFEM::FEMethodsSequence FEMethodsSequence
Definition SnesCtx.cpp:12
PetscLogEvent MOFEM_EVENT_SnesRhs
Log events to assemble residual.
Definition SnesCtx.cpp:76
friend MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition SnesCtx.cpp:592
bool vErify
If true verify vector.
Definition SnesCtx.cpp:114
BasicMethodsSequence & getPreProcSetOperators()
Definition SnesCtx.cpp:28
friend MoFEMErrorCode SnesMoFEMSetAssemblyType(SNES snes, MatAssemblyType type)
Set assembly type at the end of SnesMat.
Definition SnesCtx.cpp:564
Interface for nonlinear (SNES) solver.
Definition SnesCtx.hpp:15
MoFEM::BasicMethodsSequence BasicMethodsSequence
Definition SnesCtx.hpp:19
friend PetscErrorCode SnesLoadTangent(SNES snes, Vec u, Vec F, void *ctx)
This function calls finite element pipeline to compute tangent of of load vector in arc length contro...
Definition SnesCtx.cpp:353
FEMethodsSequence & getSetOperators()
Definition SnesCtx.cpp:131
friend MoFEMErrorCode SnesMoFEMSetBehavior(SNES snes, MoFEMTypes bh)
Set behavior if finite element in sequence does not exist.
Definition SnesCtx.cpp:578
friend MoFEMErrorCode MoFEMSNESMonitorEnergy(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition SnesCtx.cpp:648
SnesCtx(Interface &m_field, std::string problem_name)
Definition SnesCtx.cpp:125
BasicMethodsSequence & getPreProcLoadTangent()
Get the BasicMethod sequence for preprocessing of FunctionFn.
Definition SnesCtx.cpp:159
BasicMethodsSequence & getPreProcSetOperators()
Definition SnesCtx.cpp:147
MoFEM::FEMethodsSequence FEMethodsSequence
Definition SnesCtx.hpp:18
BasicMethodsSequence & getPostProcLoadTangent()
Get the BasicMethod sequence for postprocessing of FunctionFn.
Definition SnesCtx.cpp:163
MoFEMErrorCode copyLoops(const SnesCtx &snes_ctx)
Copy sequences from other SNES context.
Definition SnesCtx.cpp:167
boost::function< MoFEMErrorCode( SNES snes, Vec x, Vec F, boost::shared_ptr< CacheTuple > cache_ptr, void *ctx)> HookFunction
Definition SnesCtx.hpp:20
friend PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition SnesCtx.cpp:223
friend PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition SnesCtx.cpp:483
BasicMethodsSequence & getPostProcComputeRhs()
Definition SnesCtx.cpp:143
HookFunction & getLoadTangentHook()
Get the Load Tangent Hook function.
Definition SnesCtx.cpp:173
BasicMethodsSequence & getPostProcSetOperators()
Definition SnesCtx.cpp:151
boost::movelib::unique_ptr< SnesCtxImpl > snesCtxImpl
Definition SnesCtx.hpp:150
BasicMethodsSequence & getPreProcComputeRhs()
Definition SnesCtx.cpp:139
FEMethodsSequence & getLoadTangent()
Get the finite element pipeline for FunctionFn, that calculate tangent of function load for arc lengt...
Definition SnesCtx.cpp:155
FEMethodsSequence & getComputeRhs()
Definition SnesCtx.cpp:135
HookFunction & getRhsHook()
Get the Right Hand Side Hook function.
Definition SnesCtx.cpp:177
virtual ~SnesCtx()
friend MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition SnesCtx.cpp:592
MoFEMErrorCode clearLoops()
Clear loops.
Definition SnesCtx.cpp:171
friend MoFEMErrorCode SnesMoFEMSetAssemblyType(SNES snes, MatAssemblyType type)
Set assembly type at the end of SnesMat.
Definition SnesCtx.cpp:564
Auxiliary tools.
Definition Tools.hpp:19
Vector manager is used to create vectors \mofem_vectors.