v0.16.0
Loading...
Searching...
No Matches
Macros | Functions
ArcLengthTools.cpp File Reference
User modules » Arc-Length control
#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)
 

Detailed Description

Implementation of arc-length control method

FIXME: Some variables not comply with naming convention, need to be fixed.

Definition in file ArcLengthTools.cpp.

Macro Definition Documentation

◆ ArcFunctionBegin

#define ArcFunctionBegin
Value:
MoFEMFunctionBegin; \
MOFEM_LOG_CHANNEL("ARC_LENGTH"); \
MOFEM_LOG_FUNCTION(); \
MOFEM_LOG_TAG("ARC_LENGTH", "ArcLength");
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 11 of file ArcLengthTools.cpp.

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

Function Documentation

◆ ArcLengthMatMultShellOp()

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

mult operator for Arc Length Shell Mat

Examples
mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 174 of file ArcLengthTools.cpp.

174 {
175 MoFEMFunctionBegin;
176 void *void_ctx;
177 CHKERR MatShellGetContext(A, &void_ctx);
178 ArcLengthMatShell *ctx = static_cast<ArcLengthMatShell *>(void_ctx);
179 CHKERR MatMult(ctx->Aij, x, f);
180 double lambda;
181 CHKERR ctx->setLambda(x, &lambda, SCATTER_FORWARD);
182 double db_dot_x;
183 CHKERR VecDot(ctx->arcPtrRaw->db, x, &db_dot_x);
184 double f_lambda;
185 f_lambda = ctx->arcPtrRaw->dIag * lambda + db_dot_x;
186 CHKERR ctx->setLambda(f, &f_lambda, SCATTER_REVERSE);
187 CHKERR VecAXPY(f, lambda, ctx->arcPtrRaw->F_lambda);
188 MoFEMFunctionReturn(0);
189}

◆ 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

Examples
mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 229 of file ArcLengthTools.cpp.

229 {
230 ArcFunctionBegin;
231 void *void_ctx;
232 CHKERR PCShellGetContext(pc, &void_ctx);
233 PCArcLengthCtx *ctx = static_cast<PCArcLengthCtx *>(void_ctx);
234 void *void_MatCtx;
235 CHKERR MatShellGetContext(ctx->shellAij, &void_MatCtx);
236 ArcLengthMatShell *mat_ctx = static_cast<ArcLengthMatShell *>(void_MatCtx);
237 PetscBool same;
238 CHKERR PetscObjectTypeCompare((PetscObject)ctx->kSP, KSPPREONLY, &same);
239
240 double res_lambda;
241 CHKERR mat_ctx->setLambda(pc_f, &res_lambda, SCATTER_FORWARD);
242
243 // Solve residual
244 CHKERR KSPSetInitialGuessNonzero(ctx->kSP, PETSC_FALSE);
245 CHKERR KSPSetInitialGuessKnoll(ctx->kSP, PETSC_FALSE);
246 CHKERR KSPSolve(ctx->kSP, pc_f, pc_x);
247 double db_dot_pc_x;
248 CHKERR VecDot(ctx->arcPtrRaw->db, pc_x, &db_dot_pc_x);
249
250 // Solve for x_lambda
251 if (same != PETSC_TRUE) {
252 CHKERR KSPSetInitialGuessNonzero(ctx->kSP, PETSC_TRUE);
253 } else {
254 CHKERR KSPSetInitialGuessNonzero(ctx->kSP, PETSC_FALSE);
255 }
256 CHKERR KSPSolve(ctx->kSP, ctx->arcPtrRaw->F_lambda, ctx->arcPtrRaw->xLambda);
257 double db_dot_x_lambda;
258 CHKERR VecDot(ctx->arcPtrRaw->db, ctx->arcPtrRaw->xLambda, &db_dot_x_lambda);
259
260 // Calculate d_lambda
261 double denominator = ctx->arcPtrRaw->dIag + db_dot_x_lambda;
262 double ddlambda = -(res_lambda - db_dot_pc_x) / denominator;
263
264 // Update solution vector
265 CHKERR VecAXPY(pc_x, -ddlambda, ctx->arcPtrRaw->xLambda);
266 CHKERR mat_ctx->setLambda(pc_x, &ddlambda, SCATTER_REVERSE);
267
268 if (ddlambda != ddlambda || denominator == 0) {
269
270 double nrm2_pc_f, nrm2_db, nrm2_pc_x, nrm2_xLambda;
271 CHKERR VecNorm(pc_f, NORM_2, &nrm2_pc_f);
272 CHKERR VecNorm(ctx->arcPtrRaw->db, NORM_2, &nrm2_db);
273 CHKERR VecNorm(pc_x, NORM_2, &nrm2_pc_x);
274 CHKERR VecNorm(ctx->arcPtrRaw->xLambda, NORM_2, &nrm2_xLambda);
275
276 MOFEM_LOG("ARC_LENGTH", Sev::error)
277 << "problem with ddlambda=" << ddlambda;
278 MOFEM_LOG("ARC_LENGTH", Sev::error) << "res_lambda=" << res_lambda;
279 MOFEM_LOG("ARC_LENGTH", Sev::error) << "denominator=" << denominator;
280 MOFEM_LOG("ARC_LENGTH", Sev::error) << "db_dot_pc_x=" << db_dot_pc_x;
281 MOFEM_LOG("ARC_LENGTH", Sev::error)
282 << "db_dot_x_lambda=" << db_dot_x_lambda;
283 MOFEM_LOG("ARC_LENGTH", Sev::error)
284 << "diag=" << ctx->arcPtrRaw->dIag;
285 MOFEM_LOG("ARC_LENGTH", Sev::error) << "nrm2_db=" << nrm2_db;
286 MOFEM_LOG("ARC_LENGTH", Sev::error) << "nrm2_pc_f=" << nrm2_pc_f;
287 MOFEM_LOG("ARC_LENGTH", Sev::error) << "nrm2_pc_x=" << nrm2_pc_x;
288 MOFEM_LOG("ARC_LENGTH", Sev::error)
289 << "nrm2_xLambda=" << nrm2_xLambda;
290
291 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
292 "Increment of lambda is not number");
293 }
294
295 // Debugging PC
296 if (0) {
297 Vec y;
298 CHKERR VecDuplicate(pc_x, &y);
299 CHKERR MatMult(ctx->shellAij, pc_x, y);
300 CHKERR VecAXPY(y, -1, pc_f);
301 double res_lambda_y;
302 CHKERR mat_ctx->setLambda(y, &res_lambda_y, SCATTER_FORWARD);
303 double zero;
304 CHKERR mat_ctx->setLambda(y, &zero, SCATTER_REVERSE);
305 double norm_y;
306 CHKERR VecNorm(y, NORM_2, &norm_y);
307 MOFEM_LOG_C("ARC_LENGTH", Sev::noisy,
308 "Debug res y = %3.4e res_lambda_y = %3.4e", norm_y,
309 res_lambda_y);
310 CHKERR VecDestroy(&y);
311 }
312
313 MoFEMFunctionReturn(0);
314}

◆ PCSetupArcLength()

MoFEMErrorCode PCSetupArcLength ( PC  pc)

set up structure for Arc Length pre-conditioner

it sets pre-conditioner for matrix K

Examples
mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 316 of file ArcLengthTools.cpp.

316 {
317 MoFEMFunctionBegin;
318 void *void_ctx;
319 CHKERR PCShellGetContext(pc, &void_ctx);
320 PCArcLengthCtx *ctx = static_cast<PCArcLengthCtx *>(void_ctx);
321 auto get_pc_ops = [&](auto pc) {
322 MoFEMFunctionBegin;
323 Mat shell_aij_raw, aij_raw;
324 CHKERR PCGetOperators(pc, &shell_aij_raw, &aij_raw);
325 ctx->shellAij = SmartPetscObj<Mat>(shell_aij_raw, true);
326 ctx->Aij = SmartPetscObj<Mat>(aij_raw, true);
327 MoFEMFunctionReturn(0);
328 };
329 CHKERR get_pc_ops(pc);
330 CHKERR PCSetUseAmat(pc, PETSC_TRUE);
331 CHKERR PCSetOperators(ctx->pC, ctx->Aij, ctx->Aij);
332 CHKERR PCSetFromOptions(ctx->pC);
333 CHKERR PCSetUp(ctx->pC);
334#if PETSC_VERSION_LT(3, 12, 0)
335 CHKERR KSPSetTabLevel(ctx->kSP, 3);
336#else
337 CHKERR PetscObjectSetTabLevel((PetscObject)ctx->kSP, 3);
338#endif
339 CHKERR KSPSetFromOptions(ctx->kSP);
340 CHKERR KSPSetOperators(ctx->kSP, ctx->Aij, ctx->Aij);
341 CHKERR KSPSetPC(ctx->kSP, ctx->pC);
342 CHKERR KSPSetUp(ctx->kSP);
343 MoFEMFunctionReturn(0);
344}
Generated by Doxygen 1.9.8 and hosted at University of Glasgow