v0.14.0
Searching...
No Matches
ArcLengthTools.cpp File Reference
#include <MoFEM.hpp>
#include <ArcLengthTools.hpp>

Go to the source code of this file.

## Macros

#define ArcFunctionBegin

## 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 control method

FIXME: Some variables not comply with naming convention, need to be fixed.

Definition in file ArcLengthTools.cpp.

## ◆ ArcFunctionBegin

 #define ArcFunctionBegin
Value:
MOFEM_LOG_CHANNEL("WORLD"); \
MOFEM_LOG_FUNCTION(); \
MOFEM_LOG_TAG("WORLD", "Arc");
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346

Definition at line 13 of file ArcLengthTools.cpp.

## ◆ ArcLengthMatMultShellOp()

 MoFEMErrorCode ArcLengthMatMultShellOp ( Mat A, Vec x, Vec f )

mult operator for Arc Length Shell Mat

Definition at line 184 of file ArcLengthTools.cpp.

184 {
186 void *void_ctx;
187 CHKERR MatShellGetContext(A, &void_ctx);
188 ArcLengthMatShell *ctx = static_cast<ArcLengthMatShell *>(void_ctx);
189 CHKERR MatMult(ctx->Aij, x, f);
190 double lambda;
191 CHKERR ctx->setLambda(x, &lambda, SCATTER_FORWARD);
192 double db_dot_x;
193 CHKERR VecDot(ctx->arcPtrRaw->db, x, &db_dot_x);
194 double f_lambda;
195 f_lambda = ctx->arcPtrRaw->dIag * lambda + db_dot_x;
196 CHKERR ctx->setLambda(f, &f_lambda, SCATTER_REVERSE);
197 CHKERR VecAXPY(f, lambda, ctx->arcPtrRaw->F_lambda);
199}
#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
static double lambda
SmartPetscObj< Vec > db
db derivative of f(dx*dx), i.e. db = d[ f(dx*dx) ]/dx
SmartPetscObj< Vec > F_lambda
double dIag
diagonal value
shell matrix for arc-length method
ArcLengthCtx * arcPtrRaw
SmartPetscObj< Mat > Aij
MoFEMErrorCode setLambda(Vec ksp_x, double *lambda, ScatterMode scattermode)

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

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:311
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
const FTensor::Tensor2< T, Dim, Dim > Vec
SmartPetscObj< Vec > xLambda
solution of eq. K*xLambda = F_lambda
structure for Arc Length pre-conditioner
ArcLengthCtx * arcPtrRaw
SmartPetscObj< Mat > shellAij
SmartPetscObj< KSP > kSP

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

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}
intrusive_ptr for managing petsc objects
SmartPetscObj< Mat > Aij
SmartPetscObj< PC > pC