v0.13.2
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 <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
 
ArcLengthCtxarcPtrRaw
 

Private Attributes

boost::shared_ptr< ArcLengthCtxarcPtr
 

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}
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
auto createKSP(MPI_Comm comm)
auto createPC(MPI_Comm comm)
ArcLengthCtx * arcPtrRaw
boost::shared_ptr< ArcLengthCtx > arcPtr
SmartPetscObj< Mat > Aij
SmartPetscObj< Mat > shellAij
SmartPetscObj< KSP > kSP
SmartPetscObj< PC > pC

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

◆ 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 {
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}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
intrusive_ptr for managing petsc objects

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: