v0.16.0
Loading...
Searching...
No Matches
Macros | Functions
ArcLengthTools.cpp File Reference
#include <MoFEM.hpp>
#include <ArcLengthTools.hpp>

Go to the source code of this file.

Macros

#define ArcFunctionBegin
 

Functions

MoFEMErrorCode ArcLengthMatMultShellOp (Mat A, Vec x, Vec f)
 
MoFEMErrorCode PCApplyArcLength (PC pc, Vec pc_f, Vec pc_x)
 
MoFEMErrorCode PCSetupArcLength (PC pc)
 

Macro Definition Documentation

◆ ArcFunctionBegin

#define ArcFunctionBegin
Value:
MoFEMFunctionBegin; \
MOFEM_LOG_CHANNEL("WORLD"); \
MOFEM_LOG_FUNCTION(); \
MOFEM_LOG_TAG("WORLD", "Arc");
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359

Definition at line 13 of file ArcLengthTools.cpp.

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

Function Documentation

◆ ArcLengthMatMultShellOp()

MoFEMErrorCode ArcLengthMatMultShellOp ( Mat  A,
Vec  x,
Vec  f 
)

mult operator for Arc Length Shell Mat

Definition at line 184 of file ArcLengthTools.cpp.

184 {
185 MoFEMFunctionBegin;
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);
198 MoFEMFunctionReturn(0);
199}

◆ PCApplyArcLength()

MoFEMErrorCode PCApplyArcLength ( PC  pc,
Vec  pc_f,
Vec  pc_x 
)

apply operator for Arc Length pre-conditioner solves K*pc_x = pc_f solves K*xLambda = -dF_lambda solves ddlambda = ( res_lambda - db*xLambda )/( diag + db*pc_x ) calculate pc_x = pc_x + ddlambda*xLambda

Definition at line 243 of file ArcLengthTools.cpp.

243 {
244 ArcFunctionBegin;
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
323 MoFEMFunctionReturn(0);
324}

◆ PCSetupArcLength()

MoFEMErrorCode PCSetupArcLength ( PC  pc)

set up structure for Arc Length pre-conditioner

it sets pre-conditioner for matrix K

Definition at line 326 of file ArcLengthTools.cpp.

326 {
327 MoFEMFunctionBegin;
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) {
332 MoFEMFunctionBegin;
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);
337 MoFEMFunctionReturn(0);
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);
353 MoFEMFunctionReturn(0);
354}
Generated by Doxygen 1.9.8 and hosted at University of Glasgow