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.

                                                            {
  MoFEMFunctionBegin;
  void *void_ctx;
  CHKERR MatShellGetContext(A, &void_ctx);
  ArcLengthMatShell *ctx = static_cast<ArcLengthMatShell *>(void_ctx);
  CHKERR MatMult(ctx->Aij, x, f);
  double lambda;
  CHKERR ctx->setLambda(x, &lambda, SCATTER_FORWARD);
  double db_dot_x;
  CHKERR VecDot(ctx->arcPtrRaw->db, x, &db_dot_x);
  double f_lambda;
  f_lambda = ctx->arcPtrRaw->dIag * lambda + db_dot_x;
  CHKERR ctx->setLambda(f, &f_lambda, SCATTER_REVERSE);
  CHKERR VecAXPY(f, lambda, ctx->arcPtrRaw->F_lambda);
  MoFEMFunctionReturn(0);
}

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.

                                                           {
  ArcFunctionBegin;
  void *void_ctx;
  CHKERR PCShellGetContext(pc, &void_ctx);
  PCArcLengthCtx *ctx = static_cast<PCArcLengthCtx *>(void_ctx);
  void *void_MatCtx;
  MatShellGetContext(ctx->shellAij, &void_MatCtx);
  ArcLengthMatShell *mat_ctx = static_cast<ArcLengthMatShell *>(void_MatCtx);
  PetscBool same;
  PetscObjectTypeCompare((PetscObject)ctx->kSP, KSPPREONLY, &same);
 
  double res_lambda;
  CHKERR mat_ctx->setLambda(pc_f, &res_lambda, SCATTER_FORWARD);
 
  // Solve residual
  CHKERR KSPSetInitialGuessNonzero(ctx->kSP, PETSC_FALSE);
  CHKERR KSPSetInitialGuessKnoll(ctx->kSP, PETSC_FALSE);
  CHKERR KSPSolve(ctx->kSP, pc_f, pc_x);
  double db_dot_pc_x;
  CHKERR VecDot(ctx->arcPtrRaw->db, pc_x, &db_dot_pc_x);
 
  // Solve for x_lambda
  if (same != PETSC_TRUE) {
    CHKERR KSPSetInitialGuessNonzero(ctx->kSP, PETSC_TRUE);
  } else {
    CHKERR KSPSetInitialGuessNonzero(ctx->kSP, PETSC_FALSE);
  }
  CHKERR KSPSolve(ctx->kSP, ctx->arcPtrRaw->F_lambda, ctx->arcPtrRaw->xLambda);
  double db_dot_x_lambda;
  CHKERR VecDot(ctx->arcPtrRaw->db, ctx->arcPtrRaw->xLambda, &db_dot_x_lambda);
 
  // Calculate d_lambda
  double denominator = ctx->arcPtrRaw->dIag + db_dot_x_lambda;
  double ddlambda = -(res_lambda - db_dot_pc_x) / denominator;
 
  // Update solution vector
  CHKERR VecAXPY(pc_x, -ddlambda, ctx->arcPtrRaw->xLambda);
  CHKERR mat_ctx->setLambda(pc_x, &ddlambda, SCATTER_REVERSE);
 
  if (ddlambda != ddlambda || denominator == 0) {
 
    double nrm2_pc_f, nrm2_db, nrm2_pc_x, nrm2_xLambda;
    CHKERR VecNorm(pc_f, NORM_2, &nrm2_pc_f);
    CHKERR VecNorm(ctx->arcPtrRaw->db, NORM_2, &nrm2_db);
    CHKERR VecNorm(pc_x, NORM_2, &nrm2_pc_x);
    CHKERR VecNorm(ctx->arcPtrRaw->xLambda, NORM_2, &nrm2_xLambda);
 
    MOFEM_LOG("WORLD", Sev::error) << "problem with ddlambda=" << ddlambda;
    MOFEM_LOG("WORLD", Sev::error) << "res_lambda=" << res_lambda;
    MOFEM_LOG("WORLD", Sev::error) << "denominator=" << denominator;
    MOFEM_LOG("WORLD", Sev::error) << "db_dot_pc_x=" << db_dot_pc_x;
    MOFEM_LOG("WORLD", Sev::error) << "db_dot_x_lambda=" << db_dot_x_lambda;
    MOFEM_LOG("WORLD", Sev::error) << "diag=" << ctx->arcPtrRaw->dIag;
    MOFEM_LOG("WORLD", Sev::error) << "nrm2_db=" << nrm2_db;
    MOFEM_LOG("WORLD", Sev::error) << "nrm2_pc_f=" << nrm2_pc_f;
    MOFEM_LOG("WORLD", Sev::error) << "nrm2_pc_x=" << nrm2_pc_x;
    MOFEM_LOG("WORLD", Sev::error) << "nrm2_xLambda=" << nrm2_xLambda;
 
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Increment of lambda is not number");
 
  }
 
  // Debugging PC
  if (0) {
    Vec y;
    CHKERR VecDuplicate(pc_x, &y);
    CHKERR MatMult(ctx->shellAij, pc_x, y);
    CHKERR VecAXPY(y, -1, pc_f);
    double res_lambda_y;
    CHKERR mat_ctx->setLambda(y, &res_lambda_y, SCATTER_FORWARD);
    double zero;
    CHKERR mat_ctx->setLambda(y, &zero, SCATTER_REVERSE);
    double norm_y;
    CHKERR VecNorm(y, NORM_2, &norm_y);
    MOFEM_LOG_C("WORLD", Sev::noisy, "Debug res y = %3.4e res_lambda_y = %3.4e",
                norm_y, res_lambda_y);
    CHKERR VecDestroy(&y);
  }
 
  MoFEMFunctionReturn(0);
}

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.

                                       {
  MoFEMFunctionBegin;
  void *void_ctx;
  CHKERR PCShellGetContext(pc, &void_ctx);
  PCArcLengthCtx *ctx = static_cast<PCArcLengthCtx *>(void_ctx);
  auto get_pc_ops = [&](auto pc) {
    MoFEMFunctionBegin;
    Mat shell_aij_raw, aij_raw;
    CHKERR PCGetOperators(pc, &shell_aij_raw, &aij_raw);
    ctx->shellAij = SmartPetscObj<Mat>(shell_aij_raw, true);
    ctx->Aij = SmartPetscObj<Mat>(aij_raw, true);
    MoFEMFunctionReturn(0);
  };
  CHKERR get_pc_ops(pc);
  CHKERR PCSetUseAmat(pc, PETSC_TRUE);
  CHKERR PCSetOperators(ctx->pC, ctx->Aij, ctx->Aij);
  CHKERR PCSetFromOptions(ctx->pC);
  CHKERR PCSetUp(ctx->pC);
#if PETSC_VERSION_LT(3, 12, 0)
  CHKERR KSPSetTabLevel(ctx->kSP, 3);
#else
  CHKERR PetscObjectSetTabLevel((PetscObject)ctx->kSP, 3);
#endif
  CHKERR KSPSetFromOptions(ctx->kSP);
  CHKERR KSPSetOperators(ctx->kSP, ctx->Aij, ctx->Aij);
  CHKERR KSPSetPC(ctx->kSP, ctx->pC);
  CHKERR KSPSetUp(ctx->kSP);
  MoFEMFunctionReturn(0);
}