v0.14.0
Public Member Functions | Public Attributes | Private Attributes | Friends | List of all members
PCArcLengthCtx Struct Reference

structure for Arc Length pre-conditioner More...

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

Collaboration diagram for PCArcLengthCtx:
[legend]

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.

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

◆ 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.

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

◆ 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.

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

◆ 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.

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

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.

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)
friend

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 }

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:
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
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:257
PCArcLengthCtx::arcPtr
boost::shared_ptr< ArcLengthCtx > arcPtr
Definition: ArcLengthTools.hpp:260
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
ArcLengthMatShell::setLambda
MoFEMErrorCode setLambda(Vec ksp_x, double *lambda, ScatterMode scattermode)
Definition: ArcLengthTools.cpp:143
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
MoFEM::createPC
auto createPC(MPI_Comm comm)
Definition: PetscSmartObj.hpp:263
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
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
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:416
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:346
Generated by Doxygen 1.8.17 and hosted at University of Glasgow