structure for Arc Length pre-conditioner More...

#include <users_modules/basic_finite_elements/src/ArcLengthTools.hpp>

Collaboration diagram for PCArcLengthCtx:

Public Member Functions
	PCArcLengthCtx (Mat shell_Aij, Mat aij, boost::shared_ptr< ArcLengthCtx > arc_ptr)

	PCArcLengthCtx (PC pc, Mat shell_Aij, Mat aij, boost::shared_ptr< ArcLengthCtx > arc_ptr)

DEPRECATED	PCArcLengthCtx (Mat shell_Aij, Mat aij, ArcLengthCtx *arc_ptr_raw)

DEPRECATED	PCArcLengthCtx (PC pc, Mat shell_Aij, Mat aij, ArcLengthCtx *arc_ptr_raw)

Public Attributes
SmartPetscObj< KSP >	kSP

SmartPetscObj< PC >	pC

SmartPetscObj< Mat >	shellAij

SmartPetscObj< Mat >	Aij

ArcLengthCtx *	arcPtrRaw

Private Attributes
boost::shared_ptr< ArcLengthCtx >	arcPtr

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

MoFEMErrorCode	PCSetupArcLength (PC pc)

Detailed Description

structure for Arc Length pre-conditioner

Definition at line 235 of file ArcLengthTools.hpp.

Constructor & Destructor Documentation

◆ PCArcLengthCtx() [1/4]

PCArcLengthCtx::PCArcLengthCtx	(	Mat	shell_Aij,
		Mat	aij,
		boost::shared_ptr< ArcLengthCtx >	arc_ptr
	)

Definition at line 222 of file ArcLengthTools.cpp.

    : shellAij(shell_Aij, true), Aij(aij, true), arcPtrRaw(arc_ptr.get()),
      arcPtr(arc_ptr) {
  auto comm = PetscObjectComm((PetscObject)aij);
  pC = createPC(comm);
  kSP = createKSP(comm);
  ierr = KSPAppendOptionsPrefix(kSP, "arc_length_");
  CHKERRABORT(PETSC_COMM_WORLD, ierr);
}

◆ PCArcLengthCtx() [2/4]

PCArcLengthCtx::PCArcLengthCtx	(	PC	pc,
		Mat	shell_Aij,
		Mat	aij,
		boost::shared_ptr< ArcLengthCtx >	arc_ptr
	)

Definition at line 233 of file ArcLengthTools.cpp.

    : pC(pc, true), shellAij(shell_Aij, true), Aij(aij, true),
      arcPtrRaw(arc_ptr.get()), arcPtr(arc_ptr) {
  auto comm = PetscObjectComm((PetscObject)aij);
  kSP = createKSP(comm);
  ierr = KSPAppendOptionsPrefix(kSP, "arc_length_");
  CHKERRABORT(PETSC_COMM_WORLD, ierr);
}

◆ PCArcLengthCtx() [3/4]

PCArcLengthCtx::PCArcLengthCtx	(	Mat	shell_Aij,
		Mat	aij,
		ArcLengthCtx *	arc_ptr_raw
	)

Deprecated:: use with shared_ptr

Definition at line 203 of file ArcLengthTools.cpp.

    : shellAij(shell_Aij, true), Aij(aij, true), arcPtrRaw(arc_ptr) {
  auto comm = PetscObjectComm((PetscObject)aij);
  pC = createPC(comm);
  kSP = createKSP(comm);
  ierr = KSPAppendOptionsPrefix(kSP, "arc_length_");
  CHKERRABORT(PETSC_COMM_WORLD, ierr);
}

◆ PCArcLengthCtx() [4/4]

PCArcLengthCtx::PCArcLengthCtx	(	PC	pc,
		Mat	shell_Aij,
		Mat	aij,
		ArcLengthCtx *	arc_ptr_raw
	)

Deprecated:: use with shared_ptr

Definition at line 212 of file ArcLengthTools.cpp.

    : pC(pc, true), shellAij(shell_Aij, true), Aij(aij, true),
      arcPtrRaw(arc_ptr) {
  auto comm = PetscObjectComm((PetscObject)aij);
  kSP = createKSP(comm);
  ierr = KSPAppendOptionsPrefix(kSP, "arc_length_");
  CHKERRABORT(PETSC_COMM_WORLD, ierr);
}

Friends And Related Function Documentation

◆ PCApplyArcLength

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

friend

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 )

friend

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);
}

Member Data Documentation

◆ Aij

SmartPetscObj<Mat> PCArcLengthCtx::Aij

Definition at line 240 of file ArcLengthTools.hpp.

◆ arcPtr

boost::shared_ptr<ArcLengthCtx> PCArcLengthCtx::arcPtr

private

Definition at line 260 of file ArcLengthTools.hpp.

◆ arcPtrRaw

ArcLengthCtx* PCArcLengthCtx::arcPtrRaw

Definition at line 248 of file ArcLengthTools.hpp.

◆ kSP

SmartPetscObj<KSP> PCArcLengthCtx::kSP

Definition at line 237 of file ArcLengthTools.hpp.

◆ pC

SmartPetscObj<PC> PCArcLengthCtx::pC

Definition at line 238 of file ArcLengthTools.hpp.

◆ shellAij

SmartPetscObj<Mat> PCArcLengthCtx::shellAij

Definition at line 239 of file ArcLengthTools.hpp.

The documentation for this struct was generated from the following files:

users_modules/basic_finite_elements/src/ArcLengthTools.hpp
users_modules/basic_finite_elements/src/impl/ArcLengthTools.cpp

Public Member Functions

Public Attributes

Private Attributes

Friends