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
229 {
231 void *void_ctx;
232 CHKERR PCShellGetContext(pc, &void_ctx);
234 void *void_MatCtx;
237 PetscBool same;
238 CHKERR PetscObjectTypeCompare((PetscObject)ctx->
kSP, KSPPREONLY, &same);
239
240 double res_lambda;
242
243
244 CHKERR KSPSetInitialGuessNonzero(ctx->
kSP, PETSC_FALSE);
245 CHKERR KSPSetInitialGuessKnoll(ctx->
kSP, PETSC_FALSE);
247 double db_dot_pc_x;
249
250
251 if (same != PETSC_TRUE) {
252 CHKERR KSPSetInitialGuessNonzero(ctx->
kSP, PETSC_TRUE);
253 } else {
254 CHKERR KSPSetInitialGuessNonzero(ctx->
kSP, PETSC_FALSE);
255 }
257 double db_dot_x_lambda;
259
260
262 double ddlambda = -(res_lambda - db_dot_pc_x) / denominator;
263
264
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);
273 CHKERR VecNorm(pc_x, NORM_2, &nrm2_pc_x);
275
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;
282 << "db_dot_x_lambda=" << db_dot_x_lambda;
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;
289 << "nrm2_xLambda=" << nrm2_xLambda;
290
292 "Increment of lambda is not number");
293 }
294
295
296 if (0) {
298 CHKERR VecDuplicate(pc_x, &y);
300 CHKERR VecAXPY(y, -1, pc_f);
301 double res_lambda_y;
303 double zero;
305 double norm_y;
306 CHKERR VecNorm(y, NORM_2, &norm_y);
308 "Debug res y = %3.4e res_lambda_y = %3.4e", norm_y,
309 res_lambda_y);
311 }
312
314}
#define MOFEM_LOG_C(channel, severity, format,...)
@ MOFEM_DATA_INCONSISTENCY
#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