v0.16.0
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | Private Attributes | Friends | List of all members
PCArcLengthCtx Struct Reference

structure for Arc Length pre-conditioner More...

#include "tutorials/cor-12_cohesive_interface/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)
 
 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
 
ArcLengthCtxarcPtrRaw
 

Private Attributes

boost::shared_ptr< ArcLengthCtxarcPtr
 

Friends

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

Detailed Description

structure for Arc Length pre-conditioner

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

Definition at line 233 of file ArcLengthTools.hpp.

Constructor & Destructor Documentation

◆ PCArcLengthCtx() [1/8]

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

Definition at line 210 of file ArcLengthTools.cpp.

212 : shellAij(shell_Aij, true), Aij(aij, true), arcPtrRaw(arc_ptr.get()),
213 arcPtr(arc_ptr) {
214 auto comm = PetscObjectComm((PetscObject)aij);
215 pC = createPC(comm);
216 kSP = createKSP(comm);
217 CHKERRABORT(PETSC_COMM_WORLD, KSPAppendOptionsPrefix(kSP, "arc_length_"));
218}
auto createKSP(MPI_Comm comm)
auto createPC(MPI_Comm comm)
SmartPetscObj< Mat > Aij
ArcLengthCtx * arcPtrRaw
SmartPetscObj< KSP > kSP
boost::shared_ptr< ArcLengthCtx > arcPtr
SmartPetscObj< Mat > shellAij
SmartPetscObj< PC > pC

◆ PCArcLengthCtx() [2/8]

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

Definition at line 220 of file ArcLengthTools.cpp.

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

◆ PCArcLengthCtx() [3/8]

PCArcLengthCtx::PCArcLengthCtx ( Mat  shell_Aij,
Mat  aij,
ArcLengthCtx arc_ptr_raw 
)
Deprecated:
use with shared_ptr

Definition at line 193 of file ArcLengthTools.cpp.

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

◆ PCArcLengthCtx() [4/8]

PCArcLengthCtx::PCArcLengthCtx ( PC  pc,
Mat  shell_Aij,
Mat  aij,
ArcLengthCtx arc_ptr_raw 
)
Deprecated:
use with shared_ptr

Definition at line 201 of file ArcLengthTools.cpp.

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

◆ PCArcLengthCtx() [5/8]

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

◆ PCArcLengthCtx() [6/8]

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

◆ PCArcLengthCtx() [7/8]

DEPRECATED PCArcLengthCtx::PCArcLengthCtx ( Mat  shell_Aij,
Mat  aij,
ArcLengthCtx arc_ptr_raw 
)
Deprecated:
use with shared_ptr

◆ PCArcLengthCtx() [8/8]

DEPRECATED PCArcLengthCtx::PCArcLengthCtx ( PC  pc,
Mat  shell_Aij,
Mat  aij,
ArcLengthCtx arc_ptr_raw 
)
Deprecated:
use with shared_ptr

Friends And Related Symbol Documentation

◆ PCApplyArcLength [1/2]

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 229 of file ArcLengthTools.cpp.

229 {
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
314}
#define MOFEM_LOG_C(channel, severity, format,...)
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
const FTensor::Tensor2< T, Dim, Dim > Vec
SmartPetscObj< Vec > xLambda
solution of eq. K*xLambda = F_lambda
double dIag
diagonal value
SmartPetscObj< Vec > F_lambda
F_lambda reference load vector.
SmartPetscObj< Vec > db
db derivative of f(dx*dx), i.e. db = d[ f(dx*dx) ]/dx
shell matrix for arc-length method
MoFEMErrorCode setLambda(Vec ksp_x, double *lambda, ScatterMode scattermode)
structure for Arc Length pre-conditioner
#define ArcFunctionBegin

◆ PCApplyArcLength [2/2]

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 229 of file ArcLengthTools.cpp.

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

◆ PCSetupArcLength [1/2]

MoFEMErrorCode PCSetupArcLength ( PC  pc)
friend

set up structure for Arc Length pre-conditioner

it sets pre-conditioner for matrix K

Definition at line 316 of file ArcLengthTools.cpp.

316 {
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) {
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);
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);
344}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
intrusive_ptr for managing petsc objects

◆ PCSetupArcLength [2/2]

MoFEMErrorCode PCSetupArcLength ( PC  pc)
friend

set up structure for Arc Length pre-conditioner

it sets pre-conditioner for matrix K

Definition at line 316 of file ArcLengthTools.cpp.

316 {
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) {
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);
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);
344}

Member Data Documentation

◆ Aij

SmartPetscObj< Mat > PCArcLengthCtx::Aij

Definition at line 238 of file ArcLengthTools.hpp.

◆ arcPtr

boost::shared_ptr< ArcLengthCtx > PCArcLengthCtx::arcPtr
private

Definition at line 258 of file ArcLengthTools.hpp.

◆ arcPtrRaw

ArcLengthCtx * PCArcLengthCtx::arcPtrRaw

Definition at line 246 of file ArcLengthTools.hpp.

◆ kSP

SmartPetscObj< KSP > PCArcLengthCtx::kSP

Definition at line 235 of file ArcLengthTools.hpp.

◆ pC

SmartPetscObj< PC > PCArcLengthCtx::pC

Definition at line 236 of file ArcLengthTools.hpp.

◆ shellAij

SmartPetscObj< Mat > PCArcLengthCtx::shellAij

Definition at line 237 of file ArcLengthTools.hpp.


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