v0.14.0
Loading...
Searching...
No Matches
ArcLengthTools.cpp
Go to the documentation of this file.
1/** \file ArcLengthTools.hpp
2
3 Implementation of Arc Length element
4
5 */
6
7
8
9#include <MoFEM.hpp>
10using namespace MoFEM;
11#include <ArcLengthTools.hpp>
12
13#define ArcFunctionBegin \
14 MoFEMFunctionBegin; \
15 MOFEM_LOG_CHANNEL("WORLD"); \
16 MOFEM_LOG_FUNCTION(); \
17 MOFEM_LOG_TAG("WORLD", "Arc");
18
19// ********************
20// Arc-length ctx class
21
24 this->s = s;
25 MOFEM_LOG_C("WORLD", Sev::inform, "\tSet s = %6.4e", this->s);
27}
28
29MoFEMErrorCode ArcLengthCtx::setAlphaBeta(double alpha, double beta) {
31 this->alpha = alpha;
32 this->beta = beta;
33 MOFEM_LOG_C("WORLD", Sev::inform, "\tSet alpha = %6.4e beta = %6.4e",
34 this->alpha, this->beta);
36}
37
39 const std::string &problem_name,
40 const std::string &field_name)
41 : mField(m_field), dx2(0), F_lambda2(0), res_lambda(0) {
42
43 auto create_f_lambda = [&]() {
45 CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(problem_name, ROW,
46 F_lambda);
47 CHKERR VecSetOption(F_lambda, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
49 };
50
51 auto vec_duplicate = [&]() {
58 };
59
60 auto zero_vectors = [&]() {
62 CHKERR VecZeroEntries(F_lambda);
63 CHKERR VecZeroEntries(db);
64 CHKERR VecZeroEntries(xLambda);
65 CHKERR VecZeroEntries(x0);
66 CHKERR VecZeroEntries(dx);
68 };
69
70 auto find_lambda_dof = [&]() {
72
73 const Problem *problem_ptr;
74 CHKERR m_field.get_problem(problem_name, &problem_ptr);
75 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr_no_const =
76 problem_ptr->getNumeredRowDofsPtr();
77 auto bit_number = m_field.get_field_bit_number(field_name);
78 auto dIt = dofs_ptr_no_const->get<Unique_mi_tag>().lower_bound(
80 auto hi_dit = dofs_ptr_no_const->get<Unique_mi_tag>().upper_bound(
82 if (std::distance(dIt, hi_dit) != 1)
83 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
84 ("can not find unique LAMBDA (load factor) but found " +
85 boost::lexical_cast<std::string>(std::distance(dIt, hi_dit)))
86 .c_str());
87 arcDofRawPtr = (*dIt).get();
89 };
90
91 auto create_ghost_vecs = [&]() {
93 Vec ghost_d_lambda, ghost_diag;
94 if ((unsigned int)mField.get_comm_rank() == arcDofRawPtr->getPart()) {
95 CHKERR VecCreateGhostWithArray(mField.get_comm(), 1, 1, 0, PETSC_NULL,
96 &dLambda, &ghost_d_lambda);
97
98 CHKERR VecCreateGhostWithArray(mField.get_comm(), 1, 1, 0, PETSC_NULL,
99 &dIag, &ghost_diag);
100 } else {
101 int one[] = {0};
102 CHKERR VecCreateGhostWithArray(mField.get_comm(), 0, 1, 1, one, &dLambda,
103 &ghost_d_lambda);
104 CHKERR VecCreateGhostWithArray(mField.get_comm(), 0, 1, 1, one, &dIag,
105 &ghost_diag);
106 }
107 dLambda = 0;
108 dIag = 0;
109 ghosTdLambda = SmartPetscObj<Vec>(ghost_d_lambda);
110 ghostDiag = SmartPetscObj<Vec>(ghost_diag);
112 };
113
114 ierr = create_f_lambda();
115 CHKERRABORT(PETSC_COMM_SELF, ierr);
116
117 ierr = vec_duplicate();
118 CHKERRABORT(PETSC_COMM_SELF, ierr);
119
120 ierr = zero_vectors();
121 CHKERRABORT(PETSC_COMM_SELF, ierr);
122
123 ierr = find_lambda_dof();
124 CHKERRABORT(PETSC_COMM_SELF, ierr);
125
126 ierr = create_ghost_vecs();
127 CHKERRABORT(PETSC_COMM_SELF, ierr);
128}
129
130// ***********************
131// Arc-length shell matrix
132
134 string problem_name)
135 : Aij(aij, true), problemName(problem_name), arcPtrRaw(arc_ptr_raw) {}
136
138 boost::shared_ptr<ArcLengthCtx> arc_ptr,
139 string problem_name)
140 : Aij(aij, true), problemName(problem_name), arcPtrRaw(arc_ptr.get()),
141 arcPtr(arc_ptr) {}
142
144 ScatterMode scattermode) {
146
147 int part = arcPtrRaw->getPart();
148 int rank = arcPtrRaw->mField.get_comm_rank();
149
150 switch (scattermode) {
151 case SCATTER_FORWARD: {
152 Vec lambda_ghost;
153 if (rank == part) {
154 CHKERR VecCreateGhostWithArray(arcPtrRaw->mField.get_comm(), 1, 1, 0,
155 PETSC_NULL, lambda, &lambda_ghost);
156 } else {
157 int one[] = {0};
158 CHKERR VecCreateGhostWithArray(arcPtrRaw->mField.get_comm(), 0, 1, 1, one,
159 lambda, &lambda_ghost);
160 }
161 int idx = arcPtrRaw->getPetscGlobalDofIdx();
162 if (part == rank) {
163 CHKERR VecGetValues(ksp_x, 1, &idx, lambda);
164 }
165 CHKERR VecGhostUpdateBegin(lambda_ghost, INSERT_VALUES, SCATTER_FORWARD);
166 CHKERR VecGhostUpdateEnd(lambda_ghost, INSERT_VALUES, SCATTER_FORWARD);
167 CHKERR VecDestroy(&lambda_ghost);
168 } break;
169 case SCATTER_REVERSE: {
170 if (arcPtrRaw->getPetscLocalDofIdx() != -1) {
171 PetscScalar *array;
172 CHKERR VecGetArray(ksp_x, &array);
174 CHKERR VecRestoreArray(ksp_x, &array);
175 }
176 } break;
177 default:
178 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
179 }
180
182}
183
186 void *void_ctx;
187 CHKERR MatShellGetContext(A, &void_ctx);
188 ArcLengthMatShell *ctx = static_cast<ArcLengthMatShell *>(void_ctx);
189 CHKERR MatMult(ctx->Aij, x, f);
190 double lambda;
191 CHKERR ctx->setLambda(x, &lambda, SCATTER_FORWARD);
192 double db_dot_x;
193 CHKERR VecDot(ctx->arcPtrRaw->db, x, &db_dot_x);
194 double f_lambda;
195 f_lambda = ctx->arcPtrRaw->dIag * lambda + db_dot_x;
196 CHKERR ctx->setLambda(f, &f_lambda, SCATTER_REVERSE);
197 CHKERR VecAXPY(f, lambda, ctx->arcPtrRaw->F_lambda);
199}
200
201// arc-length preconditioner
202
203PCArcLengthCtx::PCArcLengthCtx(Mat shell_Aij, Mat aij, ArcLengthCtx *arc_ptr)
204 : shellAij(shell_Aij, true), Aij(aij, true), arcPtrRaw(arc_ptr) {
205 auto comm = PetscObjectComm((PetscObject)aij);
206 pC = createPC(comm);
207 kSP = createKSP(comm);
208 ierr = KSPAppendOptionsPrefix(kSP, "arc_length_");
209 CHKERRABORT(PETSC_COMM_WORLD, ierr);
210}
211
212PCArcLengthCtx::PCArcLengthCtx(PC pc, Mat shell_Aij, Mat aij,
213 ArcLengthCtx *arc_ptr)
214 : pC(pc, true), shellAij(shell_Aij, true), Aij(aij, true),
215 arcPtrRaw(arc_ptr) {
216 auto comm = PetscObjectComm((PetscObject)aij);
217 kSP = createKSP(comm);
218 ierr = KSPAppendOptionsPrefix(kSP, "arc_length_");
219 CHKERRABORT(PETSC_COMM_WORLD, ierr);
220}
221
222PCArcLengthCtx::PCArcLengthCtx(Mat shell_Aij, Mat aij,
223 boost::shared_ptr<ArcLengthCtx> arc_ptr)
224 : shellAij(shell_Aij, true), Aij(aij, true), arcPtrRaw(arc_ptr.get()),
225 arcPtr(arc_ptr) {
226 auto comm = PetscObjectComm((PetscObject)aij);
227 pC = createPC(comm);
228 kSP = createKSP(comm);
229 ierr = KSPAppendOptionsPrefix(kSP, "arc_length_");
230 CHKERRABORT(PETSC_COMM_WORLD, ierr);
231}
232
233PCArcLengthCtx::PCArcLengthCtx(PC pc, Mat shell_Aij, Mat aij,
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);
238 kSP = createKSP(comm);
239 ierr = KSPAppendOptionsPrefix(kSP, "arc_length_");
240 CHKERRABORT(PETSC_COMM_WORLD, ierr);
241}
242
243MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x) {
245 void *void_ctx;
246 CHKERR PCShellGetContext(pc, &void_ctx);
247 PCArcLengthCtx *ctx = static_cast<PCArcLengthCtx *>(void_ctx);
248 void *void_MatCtx;
249 MatShellGetContext(ctx->shellAij, &void_MatCtx);
250 ArcLengthMatShell *mat_ctx = static_cast<ArcLengthMatShell *>(void_MatCtx);
251 PetscBool same;
252 PetscObjectTypeCompare((PetscObject)ctx->kSP, KSPPREONLY, &same);
253
254 double res_lambda;
255 CHKERR mat_ctx->setLambda(pc_f, &res_lambda, SCATTER_FORWARD);
256
257 // Solve residual
258 CHKERR KSPSetInitialGuessNonzero(ctx->kSP, PETSC_FALSE);
259 CHKERR KSPSetInitialGuessKnoll(ctx->kSP, PETSC_FALSE);
260 CHKERR KSPSolve(ctx->kSP, pc_f, pc_x);
261 double db_dot_pc_x;
262 CHKERR VecDot(ctx->arcPtrRaw->db, pc_x, &db_dot_pc_x);
263
264 // Solve for x_lambda
265 if (same != PETSC_TRUE) {
266 CHKERR KSPSetInitialGuessNonzero(ctx->kSP, PETSC_TRUE);
267 } else {
268 CHKERR KSPSetInitialGuessNonzero(ctx->kSP, PETSC_FALSE);
269 }
270 CHKERR KSPSolve(ctx->kSP, ctx->arcPtrRaw->F_lambda, ctx->arcPtrRaw->xLambda);
271 double db_dot_x_lambda;
272 CHKERR VecDot(ctx->arcPtrRaw->db, ctx->arcPtrRaw->xLambda, &db_dot_x_lambda);
273
274 // Calculate d_lambda
275 double denominator = ctx->arcPtrRaw->dIag + db_dot_x_lambda;
276 double ddlambda = -(res_lambda - db_dot_pc_x) / denominator;
277
278 // Update solution vector
279 CHKERR VecAXPY(pc_x, -ddlambda, ctx->arcPtrRaw->xLambda);
280 CHKERR mat_ctx->setLambda(pc_x, &ddlambda, SCATTER_REVERSE);
281
282 if (ddlambda != ddlambda || denominator == 0) {
283
284 double nrm2_pc_f, nrm2_db, nrm2_pc_x, nrm2_xLambda;
285 CHKERR VecNorm(pc_f, NORM_2, &nrm2_pc_f);
286 CHKERR VecNorm(ctx->arcPtrRaw->db, NORM_2, &nrm2_db);
287 CHKERR VecNorm(pc_x, NORM_2, &nrm2_pc_x);
288 CHKERR VecNorm(ctx->arcPtrRaw->xLambda, NORM_2, &nrm2_xLambda);
289
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;
295 MOFEM_LOG("WORLD", Sev::error) << "diag=" << ctx->arcPtrRaw->dIag;
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;
300
301 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
302 "Increment of lambda is not number");
303
304 }
305
306 // Debugging PC
307 if (0) {
308 Vec y;
309 CHKERR VecDuplicate(pc_x, &y);
310 CHKERR MatMult(ctx->shellAij, pc_x, y);
311 CHKERR VecAXPY(y, -1, pc_f);
312 double res_lambda_y;
313 CHKERR mat_ctx->setLambda(y, &res_lambda_y, SCATTER_FORWARD);
314 double zero;
315 CHKERR mat_ctx->setLambda(y, &zero, SCATTER_REVERSE);
316 double norm_y;
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);
320 CHKERR VecDestroy(&y);
321 }
322
324}
325
328 void *void_ctx;
329 CHKERR PCShellGetContext(pc, &void_ctx);
330 PCArcLengthCtx *ctx = static_cast<PCArcLengthCtx *>(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);
335 ctx->shellAij = SmartPetscObj<Mat>(shell_aij_raw, true);
336 ctx->Aij = SmartPetscObj<Mat>(aij_raw, true);
338 };
339 CHKERR get_pc_ops(pc);
340 CHKERR PCSetUseAmat(pc, PETSC_TRUE);
341 CHKERR PCSetOperators(ctx->pC, ctx->Aij, ctx->Aij);
342 CHKERR PCSetFromOptions(ctx->pC);
343 CHKERR PCSetUp(ctx->pC);
344#if PETSC_VERSION_LT(3, 12, 0)
345 CHKERR KSPSetTabLevel(ctx->kSP, 3);
346#else
347 CHKERR PetscObjectSetTabLevel((PetscObject)ctx->kSP, 3);
348#endif
349 CHKERR KSPSetFromOptions(ctx->kSP);
350 CHKERR KSPSetOperators(ctx->kSP, ctx->Aij, ctx->Aij);
351 CHKERR KSPSetPC(ctx->kSP, ctx->pC);
352 CHKERR KSPSetUp(ctx->kSP);
354}
355
356// ***********************
357// Zero F_lambda vector
358
359ZeroFLmabda::ZeroFLmabda(boost::shared_ptr<ArcLengthCtx> arc_ptr)
360 : arcPtr(arc_ptr) {}
361
364 switch (snes_ctx) {
365 case CTX_SNESSETFUNCTION: {
366
367 auto zero_vals = [&](auto v) {
369 int size = problemPtr->getNbLocalDofsRow();
370 int ghosts = problemPtr->getNbGhostDofsRow();
371 double *array;
372 CHKERR VecGetArray(v, &array);
373 for (int i = 0; i != size + ghosts; ++i)
374 array[i] = 0;
375 CHKERR VecRestoreArray(v, &array);
377 };
378
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);
386
387 } break;
388 default:
389 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
390 "Lambda can be zeroed ONLY when the right hand side is evaluated.");
391 }
393}
394
395AssembleFlambda::AssembleFlambda(boost::shared_ptr<ArcLengthCtx> arc_ptr,
396 boost::shared_ptr<DirichletDisplacementBc> bc)
397 : arcPtr(arc_ptr) {
398 bCs.push_back(bc);
399}
400
401MoFEMErrorCode AssembleFlambda::preProcess() {
404}
405MoFEMErrorCode AssembleFlambda::operator()() {
408}
409
410MoFEMErrorCode AssembleFlambda::postProcess() {
412 switch (snes_ctx) {
413 case CTX_SNESSETFUNCTION: {
414
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);
423
424 auto set_bc = [&](auto l_snes_f, auto l_f_lambda) {
426 if (!bCs.empty()) {
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;
436 }
437 }
438 }
439 CHKERR VecRestoreArray(l_snes_f, &f_array);
440 CHKERR VecRestoreArray(l_f_lambda, &f_lambda_array);
441 }
443 };
444
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];
456 }
457 CHKERR VecRestoreArray(l_snes_f, &f_array);
458 CHKERR VecRestoreArray(l_f_lambda, &f_lambda_array);
460 };
461
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) {
473 f_array[i] = 0;
474 f_lambda_array[i] = 0;
475 }
476 CHKERR VecRestoreArray(l_snes_f, &f_array);
477 CHKERR VecRestoreArray(l_f_lambda, &f_lambda_array);
479 };
480
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);
487
488 CHKERR VecGhostRestoreLocalForm(snes_f, &l_snes_f);
489 CHKERR VecGhostRestoreLocalForm(arcPtr->F_lambda, &l_f_lambda);
490
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);
495
496 MOFEM_LOG_C("WORLD", Sev::inform, "\tF_lambda2 = %6.4g lambda = %6.4g",
497 arcPtr->F_lambda2, arcPtr->getFieldData());
498 MOFEM_LOG_C("WORLD", Sev::verbose,
499 "\tsnes_f norm = %6.4e snes_x norm = %6.4g", snes_fnorm,
500 snes_xnorm);
501
502 if (!boost::math::isfinite(snes_fnorm)) {
503 CHKERR arcPtr->mField.getInterface<Tools>()->checkVectorForNotANumber(
504 problemPtr, ROW, snes_f);
505 }
506
507 } break;
508 default:
509 SETERRQ(
510 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
511 "Lambda can be assembled only when the right hand side is evaluated.");
512
513 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Impossible case");
514 }
516}
517
518// ************************
519// Simple arc-length method
520
522 boost::shared_ptr<ArcLengthCtx> &arc_ptr, const bool assemble)
523 : FEMethod(), arcPtr(arc_ptr), aSsemble(assemble) {}
524
526
529 switch (snes_ctx) {
530 case CTX_SNESSETFUNCTION: {
531 if (aSsemble) {
532 CHKERR VecAssemblyBegin(snes_f);
533 CHKERR VecAssemblyEnd(snes_f);
534 }
537 } break;
538 case CTX_SNESSETJACOBIAN: {
539 if (aSsemble) {
540 CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
541 CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
542 }
543 } break;
544 default:
545 break;
546 }
548}
549
552 switch (snes_ctx) {
553 case CTX_SNESSETFUNCTION: {
554 arcPtr->res_lambda = calculateLambdaInt() - arcPtr->s;
555 CHKERR VecSetValue(snes_f, arcPtr->getPetscGlobalDofIdx(),
556 arcPtr->res_lambda, ADD_VALUES);
557 } break;
558 case CTX_SNESSETJACOBIAN: {
559 arcPtr->dIag = arcPtr->beta;
560 CHKERR MatSetValue(snes_B, arcPtr->getPetscGlobalDofIdx(),
561 arcPtr->getPetscGlobalDofIdx(), 1, ADD_VALUES);
562 } break;
563 default:
564 break;
565 }
567}
568
571 switch (snes_ctx) {
572 case CTX_SNESSETFUNCTION: {
573 if (aSsemble) {
574 CHKERR VecAssemblyBegin(snes_f);
575 CHKERR VecAssemblyEnd(snes_f);
576 }
577 } break;
578 case CTX_SNESSETJACOBIAN: {
579 if (aSsemble) {
580 CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
581 CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
582 }
583 CHKERR VecGhostUpdateBegin(arcPtr->ghostDiag, INSERT_VALUES,
584 SCATTER_FORWARD);
585 CHKERR VecGhostUpdateEnd(arcPtr->ghostDiag, INSERT_VALUES, SCATTER_FORWARD);
586 } break;
587 default:
588 break;
589 }
591}
592
594 return arcPtr->beta * arcPtr->dLambda;
595}
596
599 CHKERR VecZeroEntries(arcPtr->db);
600 CHKERR VecGhostUpdateBegin(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
601 CHKERR VecGhostUpdateEnd(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
603}
604
607 // Calculate dx
608 CHKERR VecGhostUpdateBegin(arcPtr->x0, INSERT_VALUES, SCATTER_FORWARD);
609 CHKERR VecGhostUpdateEnd(arcPtr->x0, INSERT_VALUES, SCATTER_FORWARD);
610
611 Vec l_x, l_x0, l_dx;
612 CHKERR VecGhostGetLocalForm(x, &l_x);
613 CHKERR VecGhostGetLocalForm(arcPtr->x0, &l_x0);
614 CHKERR VecGhostGetLocalForm(arcPtr->dx, &l_dx);
615 {
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);
620 int size =
621 problemPtr->getNbLocalDofsRow() + problemPtr->getNbGhostDofsRow();
622 for (int i = 0; i != size; ++i) {
623 dx_array[i] = x_array[i] - x0_array[i];
624 }
625 CHKERR VecRestoreArray(l_x, &x_array);
626 CHKERR VecRestoreArray(l_x0, &x0_array);
627 CHKERR VecRestoreArray(l_dx, &dx_array);
628 }
629 CHKERR VecGhostRestoreLocalForm(x, &l_x);
630 CHKERR VecGhostRestoreLocalForm(arcPtr->x0, &l_x0);
631 CHKERR VecGhostRestoreLocalForm(arcPtr->dx, &l_dx);
632
633 // Calculate dlambda
634 if (arcPtr->getPetscLocalDofIdx() != -1) {
635 double *array;
636 CHKERR VecGetArray(arcPtr->dx, &array);
637 arcPtr->dLambda = array[arcPtr->getPetscLocalDofIdx()];
638 array[arcPtr->getPetscLocalDofIdx()] = 0;
639 CHKERR VecRestoreArray(arcPtr->dx, &array);
640 }
641 CHKERR VecGhostUpdateBegin(arcPtr->ghosTdLambda, INSERT_VALUES,
642 SCATTER_FORWARD);
643 CHKERR VecGhostUpdateEnd(arcPtr->ghosTdLambda, INSERT_VALUES,
644 SCATTER_FORWARD);
645
646 // Calculate dx2
647 double x_nrm, x0_nrm;
648 CHKERR VecNorm(x, NORM_2, &x_nrm);
649 CHKERR VecNorm(arcPtr->x0, NORM_2, &x0_nrm);
650 CHKERR VecDot(arcPtr->dx, arcPtr->dx, &arcPtr->dx2);
651
652 MOFEM_LOG_C("WORLD", Sev::verbose,
653 "\tx norm = %6.4e x0 norm = %6.4e dx2 = %6.4e", x_nrm, x0_nrm,
654 arcPtr->dx2);
656}
657
658// ***************************
659// Spherical arc-length control
660
662 : FEMethod(), arcPtrRaw(arc_ptr_raw) {}
663
665 boost::shared_ptr<ArcLengthCtx> &arc_ptr)
666 : FEMethod(), arcPtrRaw(arc_ptr.get()), arcPtr(arc_ptr) {}
667
669
672 switch (snes_ctx) {
673 case CTX_SNESSETFUNCTION: {
676 } break;
677 case CTX_SNESSETJACOBIAN: {
678 } break;
679 default:
680 break;
681 }
683}
684
686 return arcPtrRaw->alpha * arcPtrRaw->dx2 + pow(arcPtrRaw->dLambda, 2) *
687 pow(arcPtrRaw->beta, 2) *
689}
690
693 CHKERR VecCopy(arcPtrRaw->dx, arcPtrRaw->db);
694 CHKERR VecScale(arcPtrRaw->db, 2 * arcPtrRaw->alpha);
696}
697
700 switch (snes_ctx) {
701 case CTX_SNESSETFUNCTION: {
703 CHKERR VecSetValue(snes_f, arcPtrRaw->getPetscGlobalDofIdx(),
704 arcPtrRaw->res_lambda, ADD_VALUES);
705 MOFEM_LOG_C("WORLD", Sev::verbose, "\tres_lambda = %6.4e\n",
707 } break;
708 case CTX_SNESSETJACOBIAN: {
709 arcPtrRaw->dIag =
711 CHKERR MatSetValue(snes_B, arcPtrRaw->getPetscGlobalDofIdx(),
712 arcPtrRaw->getPetscGlobalDofIdx(), 1, ADD_VALUES);
713 } break;
714 default:
715 break;
716 }
718}
719
722 switch (snes_ctx) {
723 case CTX_SNESSETFUNCTION: {
724 MOFEM_LOG_C("WORLD", Sev::verbose, "\tlambda = %6.4e\n",
726 } break;
727 case CTX_SNESSETJACOBIAN: {
728 CHKERR VecGhostUpdateBegin(arcPtrRaw->ghostDiag, INSERT_VALUES,
729 SCATTER_FORWARD);
730 CHKERR VecGhostUpdateEnd(arcPtrRaw->ghostDiag, INSERT_VALUES,
731 SCATTER_FORWARD);
732 CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
733 CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
734 MOFEM_LOG_C("WORLD", Sev::verbose, "\tdiag = %6.4e", arcPtrRaw->dIag);
735 } break;
736 default:
737 break;
738 }
740}
741
744 // dx
745 CHKERR VecCopy(x, arcPtrRaw->dx);
746 CHKERR VecAXPY(arcPtrRaw->dx, -1, arcPtrRaw->x0);
747 CHKERR VecGhostUpdateBegin(arcPtrRaw->dx, INSERT_VALUES, SCATTER_FORWARD);
748 CHKERR VecGhostUpdateEnd(arcPtrRaw->dx, INSERT_VALUES, SCATTER_FORWARD);
749 // dlambda
750 if (arcPtrRaw->getPetscLocalDofIdx() != -1) {
751 double *array;
752 CHKERR VecGetArray(arcPtrRaw->dx, &array);
754 array[arcPtrRaw->getPetscLocalDofIdx()] = 0;
755 CHKERR VecRestoreArray(arcPtrRaw->dx, &array);
756 }
757 CHKERR VecGhostUpdateBegin(arcPtrRaw->ghosTdLambda, INSERT_VALUES,
758 SCATTER_FORWARD);
759 CHKERR VecGhostUpdateEnd(arcPtrRaw->ghosTdLambda, INSERT_VALUES,
760 SCATTER_FORWARD);
761 // dx2
763 MOFEM_LOG_C("WORLD", Sev::verbose, "\tdlambda = %6.4e dx2 = %6.4e\n",
766}
767
771 *dlambda = std::sqrt(pow(arcPtrRaw->s, 2) /
772 (pow(arcPtrRaw->beta, 2) * arcPtrRaw->F_lambda2));
773 if (!(*dlambda == *dlambda)) {
774 MOFEM_LOG("WORLD", Sev::error)
775 << "s " << arcPtrRaw->s << " " << arcPtrRaw->beta << " "
777 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
778 "Increment of lambda is not a number");
779 }
781}
782
785 // check if local dof idx is non zero, i.e. that lambda is accessible from
786 // this processor
787 if (arcPtrRaw->getPetscLocalDofIdx() != -1) {
788 double *array;
789 CHKERR VecGetArray(x, &array);
790 double lambda_old = array[arcPtrRaw->getPetscLocalDofIdx()];
791 if (!(dlambda == dlambda)) {
792 MOFEM_LOG("WORLD", Sev::error)
793 << "s " << arcPtrRaw->s << " " << arcPtrRaw->beta << " "
795 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
796 "Increment of lambda is not a number");
797 }
798 array[arcPtrRaw->getPetscLocalDofIdx()] = lambda_old + dlambda;
799 MOFEM_LOG_C("WORLD", Sev::verbose, "\tlambda = %6.4e, %6.4e (%6.4e)\n",
800 lambda_old, array[arcPtrRaw->getPetscLocalDofIdx()], dlambda);
801 CHKERR VecRestoreArray(x, &array);
802 }
804}
MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x)
MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
MoFEMErrorCode PCSetupArcLength(PC pc)
#define ArcFunctionBegin
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
@ ROW
Definition: definitions.h:123
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
const double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto createKSP(MPI_Comm comm)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createPC(MPI_Comm comm)
constexpr AssemblyType A
constexpr auto field_name
Store variables for ArcLength analysis.
SmartPetscObj< Vec > db
db derivative of f(dx*dx), i.e. db = d[ f(dx*dx) ]/dx
double s
arc length radius
double res_lambda
f_lambda - s
DofIdx getPetscLocalDofIdx()
Get local index of load factor.
SmartPetscObj< Vec > dx
dx = x-x0
MoFEMErrorCode setAlphaBeta(double alpha, double beta)
set parameters controlling arc-length equations alpha controls off diagonal therms beta controls diag...
ArcLengthCtx(MoFEM::Interface &m_field, const std::string &problem_name, const std::string &field_name="LAMBDA")
double alpha
displacement scaling factor
SmartPetscObj< Vec > F_lambda
F_lambda reference load vector.
MoFEMErrorCode setS(double s)
set arc radius
MoFEM::Interface & mField
double F_lambda2
inner_prod(F_lambda,F_lambda);
double dIag
diagonal value
double dx2
inner_prod(dX,dX)
FieldData & getFieldData()
Get value of load factor.
SmartPetscObj< Vec > ghostDiag
SmartPetscObj< Vec > x0
displacement vector at beginning of step
double dLambda
increment of load factor
double beta
force scaling factor
SmartPetscObj< Vec > ghosTdLambda
int getPart()
Get proc owning lambda dof.
NumeredDofEntity * arcDofRawPtr
DofIdx getPetscGlobalDofIdx()
Get global index of load factor.
SmartPetscObj< Vec > xLambda
solution of eq. K*xLambda = F_lambda
shell matrix for arc-length method
ArcLengthCtx * arcPtrRaw
SmartPetscObj< Mat > Aij
ArcLengthMatShell(Mat aij, boost::shared_ptr< ArcLengthCtx > arc_ptr, string problem_name)
MoFEMErrorCode setLambda(Vec ksp_x, double *lambda, ScatterMode scattermode)
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
structure for User Loop Methods on finite elements
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
unsigned int getPart() const
keeps basic data about problem
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
intrusive_ptr for managing petsc objects
Auxiliary tools.
Definition: Tools.hpp:19
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
structure for Arc Length pre-conditioner
ArcLengthCtx * arcPtrRaw
SmartPetscObj< Mat > Aij
SmartPetscObj< Mat > shellAij
PCArcLengthCtx(Mat shell_Aij, Mat aij, boost::shared_ptr< ArcLengthCtx > arc_ptr)
SmartPetscObj< KSP > kSP
SmartPetscObj< PC > pC
MoFEMErrorCode preProcess()
MoFEMErrorCode calculateDb()
Calculate db.
MoFEMErrorCode operator()()
SimpleArcLengthControl(boost::shared_ptr< ArcLengthCtx > &arc_ptr, const bool assemble=false)
double calculateLambdaInt()
Calculate internal lambda.
MoFEMErrorCode calculateDxAndDlambda(Vec x)
boost::shared_ptr< ArcLengthCtx > arcPtr
MoFEMErrorCode postProcess()
MoFEMErrorCode postProcess()
virtual MoFEMErrorCode calculateDxAndDlambda(Vec x)
virtual MoFEMErrorCode calculateDb()
Calculate db.
DEPRECATED SphericalArcLengthControl(ArcLengthCtx *arc_ptr_raw)
virtual double calculateLambdaInt()
Calculate f_lambda(dx,lambda)
virtual MoFEMErrorCode calculateInitDlambda(double *dlambda)
virtual MoFEMErrorCode setDlambdaToX(Vec x, double dlambda)
ZeroFLmabda(boost::shared_ptr< ArcLengthCtx > arc_ptr)
boost::shared_ptr< ArcLengthCtx > arcPtr
MoFEMErrorCode preProcess()