52 friend PetscErrorCode
SnesRhs(SNES snes, Vec x, Vec f,
void *ctx);
54 friend PetscErrorCode
SnesMat(SNES snes, Vec x, Mat A, Mat
B,
void *ctx);
57 MatAssemblyType type);
61 MatAssemblyType type);
127 m_field, problem_name)) {}
182 : mField(m_field), moab(m_field.get_moab()), problemName(problem_name) {
187 auto core_log = logging::core::get();
211 loopsOperator.clear();
213 preProcessOperator.clear();
214 postProcessOperator.clear();
215 preProcessRhs.clear();
216 postProcessRhs.clear();
217 loopsLoadTangent.clear();
218 preProcessLoadTangent.clear();
219 postProcessLoadTangent.clear();
223PetscErrorCode
SnesRhs(SNES snes, Vec x, Vec f,
void *ctx) {
226 PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesRhs, 0, 0, 0, 0);
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) {
233 CHKERR VecAssemblyBegin(f);
235 MPI_Comm comm = PetscObjectComm((PetscObject)f);
236 PetscSynchronizedPrintf(comm,
"SNES Verify x\n");
238 CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
239 CHKERR snes_ctx->mField.getInterface<
Tools>()->checkVectorForNotANumber(
243 snes_ctx->problemName,
COL, x, INSERT_VALUES, SCATTER_REVERSE);
245 auto zero_ghost_vec = [](Vec
g) {
253 for (
int i = 0;
i != s; ++
i)
256 CHKERR VecGhostRestoreLocalForm(
g, &
l);
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,
266 auto set = [&](
auto &fe) {
276 CHKERR SNESGetKSP(snes, &fe.ksp);
281 CHKERR PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONAL, &same);
283 CHKERR SNESNewtonALGetLoadParameter(snes, &fe.ts_t);
287 fe.cacheWeakPtr = cache_ptr;
291 auto unset = [&](
auto &fe) {
297 for (
auto &
bit : snes_ctx->preProcessRhs) {
298 bit->vecAssembleSwitch = boost::move(vec_assemble_switch);
300 CHKERR snes_ctx->mField.problem_basic_method_preProcess(
301 snes_ctx->problemName, *
bit);
303 vec_assemble_switch = boost::move(
bit->vecAssembleSwitch);
306 for (
auto &lit : snes_ctx->loopsRhs) {
307 lit.second->vecAssembleSwitch = boost::move(vec_assemble_switch);
309 CHKERR snes_ctx->mField.loop_finite_elements(
310 snes_ctx->problemName, lit.first, *lit.second,
nullptr, snes_ctx->bH,
313 if (snes_ctx->vErify) {
315 CHKERR VecAssemblyBegin(f);
317 MPI_Comm comm = PetscObjectComm((PetscObject)f);
318 PetscSynchronizedPrintf(comm,
"SNES Verify f FE < %s >\n",
321 CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
322 CHKERR snes_ctx->mField.getInterface<
Tools>()->checkVectorForNotANumber(
326 vec_assemble_switch = boost::move(lit.second->vecAssembleSwitch);
329 for (
auto &
bit : snes_ctx->postProcessRhs) {
330 bit->vecAssembleSwitch = boost::move(vec_assemble_switch);
332 CHKERR snes_ctx->mField.problem_basic_method_postProcess(
333 snes_ctx->problemName, *
bit);
335 vec_assemble_switch = boost::move(
bit->vecAssembleSwitch);
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);
345 if (snes_ctx->rhsHook) {
346 CHKERR snes_ctx->rhsHook(snes, x, f, cache_ptr, ctx);
349 PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesRhs, 0, 0, 0, 0);
356 PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesLoadTangent, 0, 0, 0, 0);
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) {
363 CHKERR VecAssemblyBegin(f);
365 MPI_Comm comm = PetscObjectComm((PetscObject)f);
366 PetscSynchronizedPrintf(comm,
"SNES Verify x\n");
368 CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
369 CHKERR snes_ctx->mField.getInterface<
Tools>()->checkVectorForNotANumber(
373 snes_ctx->problemName,
COL, x, INSERT_VALUES, SCATTER_REVERSE);
375 auto zero_ghost_vec = [](Vec
g) {
383 for (
int i = 0;
i != s; ++
i)
386 CHKERR VecGhostRestoreLocalForm(
g, &
l);
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,
396 auto set = [&](
auto &fe) {
406 CHKERR SNESGetKSP(snes, &fe.ksp);
411 CHKERR PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONAL, &same);
413 CHKERR SNESNewtonALGetLoadParameter(snes, &fe.ts_t);
417 fe.cacheWeakPtr = cache_ptr;
421 auto unset = [&](
auto &fe) {
427 for (
auto &
bit : snes_ctx->preProcessLoadTangent) {
428 bit->vecAssembleSwitch = boost::move(vec_assemble_switch);
430 CHKERR snes_ctx->mField.problem_basic_method_preProcess(
431 snes_ctx->problemName, *
bit);
433 vec_assemble_switch = boost::move(
bit->vecAssembleSwitch);
436 for (
auto &lit : snes_ctx->loopsLoadTangent) {
437 lit.second->vecAssembleSwitch = boost::move(vec_assemble_switch);
439 CHKERR snes_ctx->mField.loop_finite_elements(
440 snes_ctx->problemName, lit.first, *lit.second,
nullptr, snes_ctx->bH,
443 if (snes_ctx->vErify) {
445 CHKERR VecAssemblyBegin(f);
447 MPI_Comm comm = PetscObjectComm((PetscObject)f);
448 PetscSynchronizedPrintf(comm,
"SNES Verify f FE < %s >\n",
451 CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
452 CHKERR snes_ctx->mField.getInterface<
Tools>()->checkVectorForNotANumber(
456 vec_assemble_switch = boost::move(lit.second->vecAssembleSwitch);
459 for (
auto &
bit : snes_ctx->postProcessLoadTangent) {
460 bit->vecAssembleSwitch = boost::move(vec_assemble_switch);
462 CHKERR snes_ctx->mField.problem_basic_method_postProcess(
463 snes_ctx->problemName, *
bit);
465 vec_assemble_switch = boost::move(
bit->vecAssembleSwitch);
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);
475 if (snes_ctx->loadTangentHook) {
476 CHKERR snes_ctx->loadTangentHook(snes, x, f, cache_ptr, ctx);
479 PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesLoadTangent, 0, 0, 0, 0);
483PetscErrorCode
SnesMat(SNES snes, Vec x, Mat A, Mat
B,
void *ctx) {
486 PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesMat, 0, 0, 0, 0);
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,
493 CHKERR SNESGetSolutionUpdate(snes, &dx);
495 auto set = [&](
auto &fe) {
507 CHKERR SNESGetKSP(snes, &fe.ksp);
509 fe.cacheWeakPtr = cache_ptr;
513 auto unset = [&](
auto &fe) {
519 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
520 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
522 snes_ctx->problemName,
COL, x, INSERT_VALUES, SCATTER_REVERSE);
524 if (*mat_assemble_switch) {
525 if (snes_ctx->zeroPreCondMatrixB)
528 for (
auto &
bit : snes_ctx->preProcessOperator) {
529 bit->matAssembleSwitch = boost::move(mat_assemble_switch);
531 CHKERR snes_ctx->mField.problem_basic_method_preProcess(
532 snes_ctx->problemName, *
bit);
534 mat_assemble_switch = boost::move(
bit->matAssembleSwitch);
537 for (
auto &lit : snes_ctx->loopsOperator) {
538 lit.second->matAssembleSwitch = boost::move(mat_assemble_switch);
540 CHKERR snes_ctx->mField.loop_finite_elements(
541 snes_ctx->problemName, lit.first, *(lit.second),
nullptr, snes_ctx->bH,
544 mat_assemble_switch = boost::move(lit.second->matAssembleSwitch);
547 for (
auto &
bit : snes_ctx->postProcessOperator) {
548 bit->matAssembleSwitch = boost::move(mat_assemble_switch);
550 CHKERR snes_ctx->mField.problem_basic_method_postProcess(
551 snes_ctx->problemName, *
bit);
553 mat_assemble_switch = boost::move(
bit->matAssembleSwitch);
556 if (*mat_assemble_switch) {
557 CHKERR MatAssemblyBegin(
B, snes_ctx->typeOfAssembly);
558 CHKERR MatAssemblyEnd(
B, snes_ctx->typeOfAssembly);
560 PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesMat, 0, 0, 0, 0);
566 auto get_snes_ctx = [](SNES snes) {
569 "Cannot get SNES context");
574 get_snes_ctx(snes)->typeOfAssembly = type;
580 auto get_snes_ctx = [](SNES snes) {
583 "Cannot get SNES context");
588 get_snes_ctx(snes)->bH = bh;
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;
601 std::vector<double> lnorms(fields_ptr->size(), 0),
602 norms(fields_ptr->size(), 0);
605 CHKERR SNESGetFunction(snes, &res, NULL, NULL);
608 CHKERR VecGetArrayRead(res, &r);
611 for (
auto fi : *fields_ptr) {
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;
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]));
625 CHKERR VecRestoreArrayRead(res, &r);
627 MPIU_Allreduce(&*lnorms.begin(), &*norms.begin(), lnorms.size(), MPIU_REAL,
628 MPIU_SUM, PetscObjectComm((PetscObject)snes));
632 << its <<
" Function norm " << boost::format(
"%14.12e") % (
double)fgnorm
635 for (
auto fi : *fields_ptr) {
637 << its <<
"\t Field norm "
638 << boost::format(
"%14.12e") % (
double)PetscSqrtReal(norms[f]) <<
" "
642 MOFEM_LOG(
"SNES_WORLD", Sev::inform) << its <<
"]";
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;
657 std::vector<double> lnorms(3 * fields_ptr->size(), 0),
658 norms(3 * fields_ptr->size(), 0);
662 CHKERR SNESGetSolution(snes, &x_update);
664 CHKERR SNESGetSolutionUpdate(snes, &x_update);
665 CHKERR SNESGetFunction(snes, &res, PETSC_NULLPTR, PETSC_NULLPTR);
669 CHKERR VecGetArrayRead(res, &r);
670 CHKERR VecGetArrayRead(x_update, &x);
674 for (
auto fi : *fields_ptr) {
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;
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];
690 CHKERR VecRestoreArrayRead(res, &r);
691 CHKERR VecRestoreArrayRead(x_update, &x);
693 MPIU_Allreduce(&*lnorms.begin(), &*norms.begin(), lnorms.size(), MPIU_REAL,
694 MPIU_SUM, PetscObjectComm((PetscObject)snes));
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];
702 energy_norm = std::sqrt(energy_norm);
703 residual_norm = std::sqrt(residual_norm);
704 increment_norm = std::sqrt(increment_norm);
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 <<
"[";
713 for (
auto fi : *fields_ptr) {
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 "
722 MOFEM_LOG(
"SNES_WORLD", Sev::inform) << its <<
" ]";
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
#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 ...
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.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
std::deque< BasicMethodPtr > BasicMethodsSequence
std::deque< PairNameFEMethodPtr > FEMethodsSequence
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
FEMethodsSequence & getComputeRhs()
MatAssemblyType typeOfAssembly
type of assembly at the end
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...
FEMethodsSequence & getLoadTangent()
BasicMethodsSequence & getPreProcLoadTangent()
friend MoFEMErrorCode SnesMoFEMSetBehavior(SNES snes, MoFEMTypes bh)
Set behavior if finite element in sequence does not exist.
PetscLogEvent MOFEM_EVENT_SnesMat
Log events to assemble tangent matrix.
BasicMethodsSequence preProcessOperator
friend MoFEMErrorCode MoFEMSNESMonitorEnergy(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
FEMethodsSequence loopsRhs
BasicMethodsSequence & getPostProcComputeRhs()
moab::Interface & moab
moab Interface
BasicMethodsSequence postProcessLoadTangent
MoFEM::BasicMethodsSequence BasicMethodsSequence
std::string problemName
problem name
friend MoFEMErrorCode SNESMoFEMSetAssmblyType(SNES snes, MatAssemblyType type)
virtual ~SnesCtxImpl()=default
MoFEMErrorCode copyLoops(const SnesCtxImpl &snes_ctx)
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.
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.
HookFunction & getRhsHook()
FEMethodsSequence loopsOperator
HookFunction loadTangentHook
SnesCtx::HookFunction HookFunction
BasicMethodsSequence & getPreProcComputeRhs()
FEMethodsSequence loopsLoadTangent
FEMethodsSequence & getSetOperators()
MoFEM::Interface & mField
database Interface
BasicMethodsSequence postProcessRhs
Sequence of methods run after residual is assembled.
BasicMethodsSequence & getPostProcSetOperators()
HookFunction & getLoadTangentHook()
MoFEMTypes bH
If set to MF_EXIST check if element exist, default MF_EXIST.
PetscLogEvent MOFEM_EVENT_SnesLoadTangent
SnesCtxImpl(Interface &m_field, std::string problem_name)
BasicMethodsSequence preProcessRhs
Sequence of methods run before residual is assembled.
MoFEMErrorCode clearLoops()
BasicMethodsSequence & getPostProcLoadTangent()
BasicMethodsSequence postProcessOperator
MoFEM::FEMethodsSequence FEMethodsSequence
PetscLogEvent MOFEM_EVENT_SnesRhs
Log events to assemble residual.
friend MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
bool vErify
If true verify vector.
BasicMethodsSequence & getPreProcSetOperators()
friend MoFEMErrorCode SnesMoFEMSetAssemblyType(SNES snes, MatAssemblyType type)
Set assembly type at the end of SnesMat.
Interface for nonlinear (SNES) solver.
MoFEM::BasicMethodsSequence BasicMethodsSequence
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...
FEMethodsSequence & getSetOperators()
friend MoFEMErrorCode SnesMoFEMSetBehavior(SNES snes, MoFEMTypes bh)
Set behavior if finite element in sequence does not exist.
friend MoFEMErrorCode MoFEMSNESMonitorEnergy(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
SnesCtx(Interface &m_field, std::string problem_name)
BasicMethodsSequence & getPreProcLoadTangent()
Get the BasicMethod sequence for preprocessing of FunctionFn.
BasicMethodsSequence & getPreProcSetOperators()
MoFEM::FEMethodsSequence FEMethodsSequence
BasicMethodsSequence & getPostProcLoadTangent()
Get the BasicMethod sequence for postprocessing of FunctionFn.
MoFEMErrorCode copyLoops(const SnesCtx &snes_ctx)
Copy sequences from other SNES context.
boost::function< MoFEMErrorCode( SNES snes, Vec x, Vec F, boost::shared_ptr< CacheTuple > cache_ptr, void *ctx)> HookFunction
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.
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.
BasicMethodsSequence & getPostProcComputeRhs()
HookFunction & getLoadTangentHook()
Get the Load Tangent Hook function.
BasicMethodsSequence & getPostProcSetOperators()
boost::movelib::unique_ptr< SnesCtxImpl > snesCtxImpl
BasicMethodsSequence & getPreProcComputeRhs()
FEMethodsSequence & getLoadTangent()
Get the finite element pipeline for FunctionFn, that calculate tangent of function load for arc lengt...
FEMethodsSequence & getComputeRhs()
HookFunction & getRhsHook()
Get the Right Hand Side Hook function.
friend MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
MoFEMErrorCode clearLoops()
Clear loops.
friend MoFEMErrorCode SnesMoFEMSetAssemblyType(SNES snes, MatAssemblyType type)
Set assembly type at the end of SnesMat.
Vector manager is used to create vectors \mofem_vectors.