v0.14.0
Classes | Functions
ArcLengthTools.hpp File Reference

Go to the source code of this file.

Classes

struct  ArcLengthCtx
 Store variables for ArcLength analysis. More...
 
struct  ArcLengthMatShell
 shell matrix for arc-length method More...
 
struct  PCArcLengthCtx
 structure for Arc Length pre-conditioner More...
 
struct  ZeroFLmabda
 Zero F_lambda. More...
 
struct  SimpleArcLengthControl
 
struct  SphericalArcLengthControl
 Implementation of spherical arc-length method. More...
 

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 element

Definition in file ArcLengthTools.hpp.

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  {
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 }

◆ 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  {
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 }

◆ 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  {
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 }
ArcLengthMatShell
shell matrix for arc-length method
Definition: ArcLengthTools.hpp:205
PCArcLengthCtx
structure for Arc Length pre-conditioner
Definition: ArcLengthTools.hpp:235
ArcFunctionBegin
#define ArcFunctionBegin
Definition: ArcLengthTools.cpp:13
ArcLengthCtx::db
SmartPetscObj< Vec > db
db derivative of f(dx*dx), i.e. db = d[ f(dx*dx) ]/dx
Definition: ArcLengthTools.hpp:85
PCArcLengthCtx::pC
SmartPetscObj< PC > pC
Definition: ArcLengthTools.hpp:238
PCArcLengthCtx::Aij
SmartPetscObj< Mat > Aij
Definition: ArcLengthTools.hpp:240
PCArcLengthCtx::kSP
SmartPetscObj< KSP > kSP
Definition: ArcLengthTools.hpp:237
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
ArcLengthMatShell::setLambda
MoFEMErrorCode setLambda(Vec ksp_x, double *lambda, ScatterMode scattermode)
Definition: ArcLengthTools.cpp:143
ArcLengthMatShell::arcPtrRaw
ArcLengthCtx * arcPtrRaw
Definition: ArcLengthTools.hpp:209
ArcLengthCtx::xLambda
SmartPetscObj< Vec > xLambda
solution of eq. K*xLambda = F_lambda
Definition: ArcLengthTools.hpp:86
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
ArcLengthMatShell::Aij
SmartPetscObj< Mat > Aij
Definition: ArcLengthTools.hpp:207
ArcLengthCtx::dIag
double dIag
diagonal value
Definition: ArcLengthTools.hpp:78
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
PCArcLengthCtx::arcPtrRaw
ArcLengthCtx * arcPtrRaw
Definition: ArcLengthTools.hpp:248
PCArcLengthCtx::shellAij
SmartPetscObj< Mat > shellAij
Definition: ArcLengthTools.hpp:239
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::SmartPetscObj< Mat >
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
ArcLengthCtx::F_lambda
SmartPetscObj< Vec > F_lambda
F_lambda reference load vector.
Definition: ArcLengthTools.hpp:83
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359