10 using namespace MoFEM;
13 #define ArcFunctionBegin \
15 MOFEM_LOG_CHANNEL("WORLD"); \
16 MOFEM_LOG_FUNCTION(); \
17 MOFEM_LOG_TAG("WORLD", "Arc");
25 MOFEM_LOG_C(
"WORLD", Sev::inform,
"\tSet s = %6.4e", this->s);
33 MOFEM_LOG_C(
"WORLD", Sev::inform,
"\tSet alpha = %6.4e beta = %6.4e",
34 this->alpha, this->beta);
39 const std::string &problem_name,
41 : mField(m_field), dx2(0), F_lambda2(0), res_lambda(0) {
43 auto create_f_lambda = [&]() {
47 CHKERR VecSetOption(
F_lambda, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
51 auto vec_duplicate = [&]() {
60 auto zero_vectors = [&]() {
70 auto find_lambda_dof = [&]() {
75 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr_no_const =
78 auto dIt = dofs_ptr_no_const->get<
Unique_mi_tag>().lower_bound(
79 FieldEntity::getLoBitNumberUId(bit_number));
80 auto hi_dit = dofs_ptr_no_const->get<
Unique_mi_tag>().upper_bound(
81 FieldEntity::getHiBitNumberUId(bit_number));
82 if (std::distance(dIt, hi_dit) != 1)
84 (
"can not find unique LAMBDA (load factor) but found " +
85 boost::lexical_cast<std::string>(std::distance(dIt, hi_dit)))
91 auto create_ghost_vecs = [&]() {
93 Vec ghost_d_lambda, ghost_diag;
114 ierr = create_f_lambda();
115 CHKERRABORT(PETSC_COMM_SELF,
ierr);
117 ierr = vec_duplicate();
118 CHKERRABORT(PETSC_COMM_SELF,
ierr);
120 ierr = zero_vectors();
121 CHKERRABORT(PETSC_COMM_SELF,
ierr);
123 ierr = find_lambda_dof();
124 CHKERRABORT(PETSC_COMM_SELF,
ierr);
126 ierr = create_ghost_vecs();
127 CHKERRABORT(PETSC_COMM_SELF,
ierr);
135 : Aij(aij, true), problemName(problem_name), arcPtrRaw(arc_ptr_raw) {}
138 boost::shared_ptr<ArcLengthCtx> arc_ptr,
140 : Aij(aij, true), problemName(problem_name), arcPtrRaw(arc_ptr.get()),
144 ScatterMode scattermode) {
150 switch (scattermode) {
151 case SCATTER_FORWARD: {
155 PETSC_NULL,
lambda, &lambda_ghost);
165 CHKERR VecGhostUpdateBegin(lambda_ghost, INSERT_VALUES, SCATTER_FORWARD);
166 CHKERR VecGhostUpdateEnd(lambda_ghost, INSERT_VALUES, SCATTER_FORWARD);
167 CHKERR VecDestroy(&lambda_ghost);
169 case SCATTER_REVERSE: {
172 CHKERR VecGetArray(ksp_x, &array);
174 CHKERR VecRestoreArray(ksp_x, &array);
187 CHKERR MatShellGetContext(
A, &void_ctx);
204 : shellAij(shell_Aij, true), Aij(aij, true), arcPtrRaw(arc_ptr) {
205 auto comm = PetscObjectComm((PetscObject)aij);
208 ierr = KSPAppendOptionsPrefix(
kSP,
"arc_length_");
209 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
214 : pC(pc, true), shellAij(shell_Aij, true), Aij(aij, true),
216 auto comm = PetscObjectComm((PetscObject)aij);
218 ierr = KSPAppendOptionsPrefix(
kSP,
"arc_length_");
219 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
223 boost::shared_ptr<ArcLengthCtx> arc_ptr)
224 : shellAij(shell_Aij, true), Aij(aij, true), arcPtrRaw(arc_ptr.get()),
226 auto comm = PetscObjectComm((PetscObject)aij);
229 ierr = KSPAppendOptionsPrefix(
kSP,
"arc_length_");
230 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
234 boost::shared_ptr<ArcLengthCtx> arc_ptr)
235 : pC(pc, true), shellAij(shell_Aij, true), Aij(aij, true),
236 arcPtrRaw(arc_ptr.get()), arcPtr(arc_ptr) {
237 auto comm = PetscObjectComm((PetscObject)aij);
239 ierr = KSPAppendOptionsPrefix(
kSP,
"arc_length_");
240 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
246 CHKERR PCShellGetContext(pc, &void_ctx);
249 MatShellGetContext(ctx->
shellAij, &void_MatCtx);
252 PetscObjectTypeCompare((PetscObject)ctx->
kSP, KSPPREONLY, &same);
258 CHKERR KSPSetInitialGuessNonzero(ctx->
kSP, PETSC_FALSE);
259 CHKERR KSPSetInitialGuessKnoll(ctx->
kSP, PETSC_FALSE);
265 if (same != PETSC_TRUE) {
266 CHKERR KSPSetInitialGuessNonzero(ctx->
kSP, PETSC_TRUE);
268 CHKERR KSPSetInitialGuessNonzero(ctx->
kSP, PETSC_FALSE);
271 double db_dot_x_lambda;
276 double ddlambda = -(res_lambda - db_dot_pc_x) / denominator;
282 if (ddlambda != ddlambda || denominator == 0) {
284 double nrm2_pc_f, nrm2_db, nrm2_pc_x, nrm2_xLambda;
285 CHKERR VecNorm(pc_f, NORM_2, &nrm2_pc_f);
287 CHKERR VecNorm(pc_x, NORM_2, &nrm2_pc_x);
290 MOFEM_LOG(
"WORLD", Sev::error) <<
"problem with ddlambda=" << ddlambda;
291 MOFEM_LOG(
"WORLD", Sev::error) <<
"res_lambda=" << res_lambda;
292 MOFEM_LOG(
"WORLD", Sev::error) <<
"denominator=" << denominator;
293 MOFEM_LOG(
"WORLD", Sev::error) <<
"db_dot_pc_x=" << db_dot_pc_x;
294 MOFEM_LOG(
"WORLD", Sev::error) <<
"db_dot_x_lambda=" << db_dot_x_lambda;
296 MOFEM_LOG(
"WORLD", Sev::error) <<
"nrm2_db=" << nrm2_db;
297 MOFEM_LOG(
"WORLD", Sev::error) <<
"nrm2_pc_f=" << nrm2_pc_f;
298 MOFEM_LOG(
"WORLD", Sev::error) <<
"nrm2_pc_x=" << nrm2_pc_x;
299 MOFEM_LOG(
"WORLD", Sev::error) <<
"nrm2_xLambda=" << nrm2_xLambda;
302 "Increment of lambda is not number");
309 CHKERR VecDuplicate(pc_x, &y);
311 CHKERR VecAXPY(y, -1, pc_f);
317 CHKERR VecNorm(y, NORM_2, &norm_y);
318 MOFEM_LOG_C(
"WORLD", Sev::noisy,
"Debug res y = %3.4e res_lambda_y = %3.4e",
319 norm_y, res_lambda_y);
329 CHKERR PCShellGetContext(pc, &void_ctx);
331 auto get_pc_ops = [&](
auto pc) {
333 Mat shell_aij_raw, aij_raw;
334 CHKERR PCGetOperators(pc, &shell_aij_raw, &aij_raw);
340 CHKERR PCSetUseAmat(pc, PETSC_TRUE);
344 #if PETSC_VERSION_LT(3, 12, 0)
347 CHKERR PetscObjectSetTabLevel((PetscObject)ctx->
kSP, 3);
365 case CTX_SNESSETFUNCTION: {
367 auto zero_vals = [&](
auto v) {
369 int size = problemPtr->getNbLocalDofsRow();
370 int ghosts = problemPtr->getNbGhostDofsRow();
373 for (
int i = 0;
i != size + ghosts; ++
i)
375 CHKERR VecRestoreArray(
v, &array);
379 Vec l_x_lambda, l_f_lambda;
380 CHKERR VecGhostGetLocalForm(
arcPtr->xLambda, &l_x_lambda);
381 CHKERR VecGhostGetLocalForm(
arcPtr->F_lambda, &l_f_lambda);
382 CHKERR zero_vals(l_x_lambda);
383 CHKERR zero_vals(l_f_lambda);
384 CHKERR VecGhostRestoreLocalForm(
arcPtr->xLambda, &l_x_lambda);
385 CHKERR VecGhostRestoreLocalForm(
arcPtr->F_lambda, &l_f_lambda);
390 "Lambda can be zeroed ONLY when the right hand side is evaluated.");
395 AssembleFlambda::AssembleFlambda(boost::shared_ptr<ArcLengthCtx> arc_ptr,
396 boost::shared_ptr<DirichletDisplacementBc> bc)
413 case CTX_SNESSETFUNCTION: {
415 CHKERR VecAssemblyBegin(arcPtr->F_lambda);
416 CHKERR VecAssemblyEnd(arcPtr->F_lambda);
417 CHKERR VecAssemblyBegin(snes_f);
418 CHKERR VecAssemblyEnd(snes_f);
419 CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, ADD_VALUES, SCATTER_REVERSE);
420 CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, ADD_VALUES, SCATTER_REVERSE);
421 CHKERR VecGhostUpdateBegin(snes_f, ADD_VALUES, SCATTER_REVERSE);
422 CHKERR VecGhostUpdateEnd(snes_f, ADD_VALUES, SCATTER_REVERSE);
424 auto set_bc = [&](
auto l_snes_f,
auto l_f_lambda) {
427 double *f_array, *f_lambda_array;
428 CHKERR VecGetArray(l_snes_f, &f_array);
429 CHKERR VecGetArray(l_f_lambda, &f_lambda_array);
430 for (
auto &bc : bCs) {
431 for (
auto idx : bc->dofsIndices) {
432 auto weak_dof = problemPtr->getRowDofsByPetscGlobalDofIdx(idx);
433 if (
auto shared_dof = weak_dof.lock()) {
434 f_array[shared_dof->getPetscLocalDofIdx()] = 0;
435 f_lambda_array[shared_dof->getPetscLocalDofIdx()] = 0;
439 CHKERR VecRestoreArray(l_snes_f, &f_array);
440 CHKERR VecRestoreArray(l_f_lambda, &f_lambda_array);
445 auto add_f_lambda = [&](
auto l_snes_f,
auto l_f_lambda) {
447 int size = problemPtr->getNbLocalDofsRow();
448 int ghosts = problemPtr->getNbGhostDofsRow();
449 double lambda = arcPtr->getFieldData();
450 int local_lambda_idx = arcPtr->getPetscLocalDofIdx();
451 double *f_array, *f_lambda_array;
452 CHKERR VecGetArray(l_snes_f, &f_array);
453 CHKERR VecGetArray(l_f_lambda, &f_lambda_array);
454 for (
int i = 0;
i != size; ++
i) {
455 f_array[
i] +=
lambda * f_lambda_array[
i];
457 CHKERR VecRestoreArray(l_snes_f, &f_array);
458 CHKERR VecRestoreArray(l_f_lambda, &f_lambda_array);
462 auto zero_ghost = [&](
auto l_snes_f,
auto l_f_lambda) {
464 int size = problemPtr->getNbLocalDofsRow();
465 int ghosts = problemPtr->getNbGhostDofsRow();
466 double lambda = arcPtr->getFieldData();
467 int local_lambda_idx = arcPtr->getPetscLocalDofIdx();
468 double *f_array, *f_lambda_array;
469 CHKERR VecGetArray(l_snes_f, &f_array);
470 CHKERR VecGetArray(l_f_lambda, &f_lambda_array);
471 f_lambda_array[local_lambda_idx] = 0;
472 for (
int i = size;
i != size + ghosts; ++
i) {
474 f_lambda_array[
i] = 0;
476 CHKERR VecRestoreArray(l_snes_f, &f_array);
477 CHKERR VecRestoreArray(l_f_lambda, &f_lambda_array);
481 Vec l_snes_f, l_f_lambda;
482 CHKERR VecGhostGetLocalForm(snes_f, &l_snes_f);
483 CHKERR VecGhostGetLocalForm(arcPtr->F_lambda, &l_f_lambda);
484 CHKERR add_f_lambda(l_snes_f, l_f_lambda);
485 CHKERR set_bc(l_snes_f, l_f_lambda);
486 CHKERR zero_ghost(l_snes_f, l_f_lambda);
488 CHKERR VecGhostRestoreLocalForm(snes_f, &l_snes_f);
489 CHKERR VecGhostRestoreLocalForm(arcPtr->F_lambda, &l_f_lambda);
491 double snes_fnorm, snes_xnorm;
492 CHKERR VecNorm(snes_f, NORM_2, &snes_fnorm);
493 CHKERR VecNorm(snes_x, NORM_2, &snes_xnorm);
494 CHKERR VecDot(arcPtr->F_lambda, arcPtr->F_lambda, &arcPtr->F_lambda2);
496 MOFEM_LOG_C(
"WORLD", Sev::inform,
"\tF_lambda2 = %6.4g lambda = %6.4g",
497 arcPtr->F_lambda2, arcPtr->getFieldData());
499 "\tsnes_f norm = %6.4e snes_x norm = %6.4g", snes_fnorm,
502 if (!boost::math::isfinite(snes_fnorm)) {
503 CHKERR arcPtr->mField.getInterface<
Tools>()->checkVectorForNotANumber(
504 problemPtr,
ROW, snes_f);
511 "Lambda can be assembled only when the right hand side is evaluated.");
522 boost::shared_ptr<ArcLengthCtx> &arc_ptr,
const bool assemble)
523 :
FEMethod(), arcPtr(arc_ptr), aSsemble(assemble) {}
530 case CTX_SNESSETFUNCTION: {
532 CHKERR VecAssemblyBegin(snes_f);
533 CHKERR VecAssemblyEnd(snes_f);
538 case CTX_SNESSETJACOBIAN: {
540 CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
541 CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
553 case CTX_SNESSETFUNCTION: {
555 CHKERR VecSetValue(snes_f,
arcPtr->getPetscGlobalDofIdx(),
556 arcPtr->res_lambda, ADD_VALUES);
558 case CTX_SNESSETJACOBIAN: {
560 CHKERR MatSetValue(snes_B,
arcPtr->getPetscGlobalDofIdx(),
561 arcPtr->getPetscGlobalDofIdx(), 1, ADD_VALUES);
572 case CTX_SNESSETFUNCTION: {
574 CHKERR VecAssemblyBegin(snes_f);
575 CHKERR VecAssemblyEnd(snes_f);
578 case CTX_SNESSETJACOBIAN: {
580 CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
581 CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
583 CHKERR VecGhostUpdateBegin(
arcPtr->ghostDiag, INSERT_VALUES,
585 CHKERR VecGhostUpdateEnd(
arcPtr->ghostDiag, INSERT_VALUES, SCATTER_FORWARD);
600 CHKERR VecGhostUpdateBegin(
arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
601 CHKERR VecGhostUpdateEnd(
arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
608 CHKERR VecGhostUpdateBegin(
arcPtr->x0, INSERT_VALUES, SCATTER_FORWARD);
609 CHKERR VecGhostUpdateEnd(
arcPtr->x0, INSERT_VALUES, SCATTER_FORWARD);
612 CHKERR VecGhostGetLocalForm(x, &l_x);
616 double *x_array, *x0_array, *dx_array;
617 CHKERR VecGetArray(l_x, &x_array);
618 CHKERR VecGetArray(l_x0, &x0_array);
619 CHKERR VecGetArray(l_dx, &dx_array);
621 problemPtr->getNbLocalDofsRow() + problemPtr->getNbGhostDofsRow();
622 for (
int i = 0;
i != size; ++
i) {
623 dx_array[
i] = x_array[
i] - x0_array[
i];
625 CHKERR VecRestoreArray(l_x, &x_array);
626 CHKERR VecRestoreArray(l_x0, &x0_array);
627 CHKERR VecRestoreArray(l_dx, &dx_array);
629 CHKERR VecGhostRestoreLocalForm(x, &l_x);
634 if (
arcPtr->getPetscLocalDofIdx() != -1) {
638 array[
arcPtr->getPetscLocalDofIdx()] = 0;
641 CHKERR VecGhostUpdateBegin(
arcPtr->ghosTdLambda, INSERT_VALUES,
643 CHKERR VecGhostUpdateEnd(
arcPtr->ghosTdLambda, INSERT_VALUES,
647 double x_nrm, x0_nrm;
648 CHKERR VecNorm(x, NORM_2, &x_nrm);
653 "\tx norm = %6.4e x0 norm = %6.4e dx2 = %6.4e", x_nrm, x0_nrm,
662 :
FEMethod(), arcPtrRaw(arc_ptr_raw) {}
665 boost::shared_ptr<ArcLengthCtx> &arc_ptr)
666 :
FEMethod(), arcPtrRaw(arc_ptr.get()), arcPtr(arc_ptr) {}
673 case CTX_SNESSETFUNCTION: {
677 case CTX_SNESSETJACOBIAN: {
701 case CTX_SNESSETFUNCTION: {
705 MOFEM_LOG_C(
"WORLD", Sev::verbose,
"\tres_lambda = %6.4e\n",
708 case CTX_SNESSETJACOBIAN: {
723 case CTX_SNESSETFUNCTION: {
724 MOFEM_LOG_C(
"WORLD", Sev::verbose,
"\tlambda = %6.4e\n",
727 case CTX_SNESSETJACOBIAN: {
732 CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
733 CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
763 MOFEM_LOG_C(
"WORLD", Sev::verbose,
"\tdlambda = %6.4e dx2 = %6.4e\n",
773 if (!(*dlambda == *dlambda)) {
778 "Increment of lambda is not a number");
789 CHKERR VecGetArray(x, &array);
791 if (!(dlambda == dlambda)) {
796 "Increment of lambda is not a number");
799 MOFEM_LOG_C(
"WORLD", Sev::verbose,
"\tlambda = %6.4e, %6.4e (%6.4e)\n",
801 CHKERR VecRestoreArray(x, &array);