1 #include <petscconvest.h> /*I "petscconvest.h" I*/ 2 #include <petscdmplex.h> 3 #include <petscds.h> 4 5 #include <petsc/private/petscconvestimpl.h> 6 7 static PetscErrorCode zero_private(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 8 { 9 PetscInt c; 10 for (c = 0; c < Nc; ++c) u[c] = 0.0; 11 return 0; 12 } 13 14 /*@ 15 PetscConvEstDestroy - Destroys a PetscConvEst object 16 17 Collective on PetscConvEst 18 19 Input Parameter: 20 . ce - The PetscConvEst object 21 22 Level: beginner 23 24 .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate() 25 @*/ 26 PetscErrorCode PetscConvEstDestroy(PetscConvEst *ce) 27 { 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 if (!*ce) PetscFunctionReturn(0); 32 PetscValidHeaderSpecific((*ce),PETSC_OBJECT_CLASSID,1); 33 if (--((PetscObject)(*ce))->refct > 0) { 34 *ce = NULL; 35 PetscFunctionReturn(0); 36 } 37 ierr = PetscFree3((*ce)->initGuess, (*ce)->exactSol, (*ce)->ctxs);CHKERRQ(ierr); 38 ierr = PetscFree2((*ce)->dofs, (*ce)->errors);CHKERRQ(ierr); 39 ierr = PetscHeaderDestroy(ce);CHKERRQ(ierr); 40 PetscFunctionReturn(0); 41 } 42 43 /*@ 44 PetscConvEstSetFromOptions - Sets a PetscConvEst object from options 45 46 Collective on PetscConvEst 47 48 Input Parameters: 49 . ce - The PetscConvEst object 50 51 Level: beginner 52 53 .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate() 54 @*/ 55 PetscErrorCode PetscConvEstSetFromOptions(PetscConvEst ce) 56 { 57 PetscErrorCode ierr; 58 59 PetscFunctionBegin; 60 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject) ce), "", "Convergence Estimator Options", "PetscConvEst");CHKERRQ(ierr); 61 ierr = PetscOptionsInt("-convest_num_refine", "The number of refinements for the convergence check", "PetscConvEst", ce->Nr, &ce->Nr, NULL);CHKERRQ(ierr); 62 ierr = PetscOptionsReal("-convest_refine_factor", "The increase in resolution in each dimension", "PetscConvEst", ce->r, &ce->r, NULL);CHKERRQ(ierr); 63 ierr = PetscOptionsBool("-convest_monitor", "Monitor the error for each convergence check", "PetscConvEst", ce->monitor, &ce->monitor, NULL);CHKERRQ(ierr); 64 ierr = PetscOptionsEnd();CHKERRQ(ierr); 65 PetscFunctionReturn(0); 66 } 67 68 /*@ 69 PetscConvEstView - Views a PetscConvEst object 70 71 Collective on PetscConvEst 72 73 Input Parameters: 74 + ce - The PetscConvEst object 75 - viewer - The PetscViewer object 76 77 Level: beginner 78 79 .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate() 80 @*/ 81 PetscErrorCode PetscConvEstView(PetscConvEst ce, PetscViewer viewer) 82 { 83 PetscErrorCode ierr; 84 85 PetscFunctionBegin; 86 ierr = PetscObjectPrintClassNamePrefixType((PetscObject) ce, viewer);CHKERRQ(ierr); 87 ierr = PetscViewerASCIIPrintf(viewer, "ConvEst with %D levels\n", ce->Nr+1);CHKERRQ(ierr); 88 PetscFunctionReturn(0); 89 } 90 91 /*@ 92 PetscConvEstGetSolver - Gets the solver used to produce discrete solutions 93 94 Not collective 95 96 Input Parameter: 97 . ce - The PetscConvEst object 98 99 Output Parameter: 100 . solver - The solver 101 102 Level: intermediate 103 104 .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate() 105 @*/ 106 PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, PetscObject *solver) 107 { 108 PetscFunctionBegin; 109 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1); 110 PetscValidPointer(solver, 2); 111 *solver = ce->solver; 112 PetscFunctionReturn(0); 113 } 114 115 /*@ 116 PetscConvEstSetSolver - Sets the solver used to produce discrete solutions 117 118 Not collective 119 120 Input Parameters: 121 + ce - The PetscConvEst object 122 - solver - The solver 123 124 Level: intermediate 125 126 Note: The solver MUST have an attached DM/DS, so that we know the exact solution 127 128 .seealso: PetscConvEstGetSNES(), PetscConvEstCreate(), PetscConvEstGetConvRate() 129 @*/ 130 PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, PetscObject solver) 131 { 132 PetscErrorCode ierr; 133 134 PetscFunctionBegin; 135 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1); 136 PetscValidHeader(solver, 2); 137 ce->solver = solver; 138 ierr = (*ce->ops->setsolver)(ce, solver);CHKERRQ(ierr); 139 PetscFunctionReturn(0); 140 } 141 142 /*@ 143 PetscConvEstSetUp - After the solver is specified, we create structures for estimating convergence 144 145 Collective on PetscConvEst 146 147 Input Parameters: 148 . ce - The PetscConvEst object 149 150 Level: beginner 151 152 .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate() 153 @*/ 154 PetscErrorCode PetscConvEstSetUp(PetscConvEst ce) 155 { 156 PetscInt Nf, f, Nds, s; 157 PetscErrorCode ierr; 158 159 PetscFunctionBegin; 160 ierr = DMGetNumFields(ce->idm, &Nf);CHKERRQ(ierr); 161 ce->Nf = PetscMax(Nf, 1); 162 ierr = PetscMalloc2((ce->Nr+1)*ce->Nf, &ce->dofs, (ce->Nr+1)*ce->Nf, &ce->errors);CHKERRQ(ierr); 163 ierr = PetscCalloc3(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol, ce->Nf, &ce->ctxs);CHKERRQ(ierr); 164 for (f = 0; f < Nf; ++f) ce->initGuess[f] = zero_private; 165 ierr = DMGetNumDS(ce->idm, &Nds);CHKERRQ(ierr); 166 for (s = 0; s < Nds; ++s) { 167 PetscDS ds; 168 DMLabel label; 169 IS fieldIS; 170 const PetscInt *fields; 171 PetscInt dsNf; 172 173 ierr = DMGetRegionNumDS(ce->idm, s, &label, &fieldIS, &ds);CHKERRQ(ierr); 174 ierr = PetscDSGetNumFields(ds, &dsNf);CHKERRQ(ierr); 175 if (fieldIS) {ierr = ISGetIndices(fieldIS, &fields);CHKERRQ(ierr);} 176 for (f = 0; f < dsNf; ++f) { 177 const PetscInt field = fields[f]; 178 ierr = PetscDSGetExactSolution(ds, field, &ce->exactSol[field], &ce->ctxs[field]);CHKERRQ(ierr); 179 } 180 if (fieldIS) {ierr = ISRestoreIndices(fieldIS, &fields);CHKERRQ(ierr);} 181 } 182 for (f = 0; f < Nf; ++f) { 183 if (!ce->exactSol[f]) SETERRQ1(PetscObjectComm((PetscObject) ce), PETSC_ERR_ARG_WRONG, "DS must contain exact solution functions in order to estimate convergence, missing for field %D", f); 184 } 185 PetscFunctionReturn(0); 186 } 187 188 PetscErrorCode PetscConvEstComputeInitialGuess(PetscConvEst ce, PetscInt r, DM dm, Vec u) 189 { 190 PetscErrorCode ierr; 191 192 PetscFunctionBegin; 193 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1); 194 if (dm) PetscValidHeaderSpecific(dm, DM_CLASSID, 3); 195 PetscValidHeaderSpecific(u, VEC_CLASSID, 4); 196 ierr = (*ce->ops->initguess)(ce, r, dm, u);CHKERRQ(ierr); 197 PetscFunctionReturn(0); 198 } 199 200 PetscErrorCode PetscConvEstComputeError(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[]) 201 { 202 PetscErrorCode ierr; 203 204 PetscFunctionBegin; 205 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1); 206 if (dm) PetscValidHeaderSpecific(dm, DM_CLASSID, 3); 207 PetscValidHeaderSpecific(u, VEC_CLASSID, 4); 208 PetscValidRealPointer(errors, 5); 209 ierr = (*ce->ops->computeerror)(ce, r, dm, u, errors);CHKERRQ(ierr); 210 PetscFunctionReturn(0); 211 } 212 213 /*@ 214 PetscConvEstMonitorDefault - Monitors the convergence estimation loop 215 216 Collective on PetscConvEst 217 218 Input Parameter: 219 + ce - The PetscConvEst object 220 - r - The refinement level 221 222 Options database keys: 223 . -convest_monitor - Activate the monitor 224 225 Level: intermediate 226 227 .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate(), SNESSolve(), TSSolve() 228 @*/ 229 PetscErrorCode PetscConvEstMonitorDefault(PetscConvEst ce, PetscInt r) 230 { 231 MPI_Comm comm; 232 PetscInt f; 233 PetscErrorCode ierr; 234 235 PetscFunctionBegin; 236 if (ce->monitor) { 237 PetscInt *dofs = &ce->dofs[r*ce->Nf]; 238 PetscReal *errors = &ce->errors[r*ce->Nf]; 239 240 ierr = PetscObjectGetComm((PetscObject) ce, &comm);CHKERRQ(ierr); 241 ierr = PetscPrintf(comm, "N: ");CHKERRQ(ierr); 242 if (ce->Nf > 1) {ierr = PetscPrintf(comm, "[");CHKERRQ(ierr);} 243 for (f = 0; f < ce->Nf; ++f) { 244 if (f > 0) {ierr = PetscPrintf(comm, ", ");CHKERRQ(ierr);} 245 ierr = PetscPrintf(comm, "%7D", dofs[f]);CHKERRQ(ierr); 246 } 247 if (ce->Nf > 1) {ierr = PetscPrintf(comm, "]");CHKERRQ(ierr);} 248 ierr = PetscPrintf(comm, " ");CHKERRQ(ierr); 249 ierr = PetscPrintf(comm, "L_2 Error: ");CHKERRQ(ierr); 250 if (ce->Nf > 1) {ierr = PetscPrintf(comm, "[");CHKERRQ(ierr);} 251 for (f = 0; f < ce->Nf; ++f) { 252 if (f > 0) {ierr = PetscPrintf(comm, ", ");CHKERRQ(ierr);} 253 if (errors[f] < 1.0e-11) {ierr = PetscPrintf(comm, "< 1e-11");CHKERRQ(ierr);} 254 else {ierr = PetscPrintf(comm, "%g", (double) errors[f]);CHKERRQ(ierr);} 255 } 256 if (ce->Nf > 1) {ierr = PetscPrintf(comm, "]");CHKERRQ(ierr);} 257 ierr = PetscPrintf(comm, "\n");CHKERRQ(ierr); 258 } 259 PetscFunctionReturn(0); 260 } 261 262 static PetscErrorCode PetscConvEstSetSNES_Private(PetscConvEst ce, PetscObject solver) 263 { 264 PetscClassId id; 265 PetscErrorCode ierr; 266 267 PetscFunctionBegin; 268 ierr = PetscObjectGetClassId(ce->solver, &id);CHKERRQ(ierr); 269 if (id != SNES_CLASSID) SETERRQ(PetscObjectComm((PetscObject) ce), PETSC_ERR_ARG_WRONG, "Solver was not a SNES"); 270 ierr = SNESGetDM((SNES) ce->solver, &ce->idm);CHKERRQ(ierr); 271 PetscFunctionReturn(0); 272 } 273 274 static PetscErrorCode PetscConvEstInitGuessSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u) 275 { 276 PetscErrorCode ierr; 277 278 PetscFunctionBegin; 279 ierr = DMProjectFunction(dm, 0.0, ce->initGuess, ce->ctxs, INSERT_VALUES, u);CHKERRQ(ierr); 280 PetscFunctionReturn(0); 281 } 282 283 static PetscErrorCode PetscConvEstComputeErrorSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[]) 284 { 285 PetscErrorCode ierr; 286 287 PetscFunctionBegin; 288 ierr = DMComputeL2FieldDiff(dm, 0.0, ce->exactSol, ce->ctxs, u, errors);CHKERRQ(ierr); 289 PetscFunctionReturn(0); 290 } 291 292 static PetscErrorCode PetscConvEstSetJacobianNullspace_Private(PetscConvEst ce, SNES snes) 293 { 294 DM dm; 295 PetscInt f; 296 PetscErrorCode ierr; 297 298 PetscFunctionBegin; 299 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 300 for (f = 0; f < ce->Nf; ++f) { 301 PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *); 302 303 ierr = DMGetNullSpaceConstructor(dm, f, &nspconstr);CHKERRQ(ierr); 304 if (nspconstr) { 305 MatNullSpace nullsp; 306 Mat J; 307 308 ierr = (*nspconstr)(dm, f, f,&nullsp);CHKERRQ(ierr); 309 ierr = SNESSetUp(snes);CHKERRQ(ierr); 310 ierr = SNESGetJacobian(snes, &J, NULL, NULL, NULL);CHKERRQ(ierr); 311 ierr = MatSetNullSpace(J, nullsp);CHKERRQ(ierr); 312 ierr = MatNullSpaceDestroy(&nullsp);CHKERRQ(ierr); 313 break; 314 } 315 } 316 PetscFunctionReturn(0); 317 } 318 319 static PetscErrorCode PetscConvEstGetConvRateSNES_Private(PetscConvEst ce, PetscReal alpha[]) 320 { 321 SNES snes = (SNES) ce->solver; 322 DM *dm; 323 PetscObject disc; 324 PetscReal *x, *y, slope, intercept; 325 PetscInt Nr = ce->Nr, r, f, dim, oldlevel, oldnlev; 326 void *ctx; 327 PetscErrorCode ierr; 328 329 PetscFunctionBegin; 330 if (ce->r != 2.0) SETERRQ1(PetscObjectComm((PetscObject) ce), PETSC_ERR_SUP, "Only refinement factor 2 is currently supported (not %g)", (double) ce->r); 331 ierr = DMGetDimension(ce->idm, &dim);CHKERRQ(ierr); 332 ierr = DMGetApplicationContext(ce->idm, &ctx);CHKERRQ(ierr); 333 ierr = DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);CHKERRQ(ierr); 334 ierr = DMGetRefineLevel(ce->idm, &oldlevel);CHKERRQ(ierr); 335 ierr = PetscMalloc1((Nr+1), &dm);CHKERRQ(ierr); 336 /* Loop over meshes */ 337 dm[0] = ce->idm; 338 for (r = 0; r <= Nr; ++r) { 339 Vec u; 340 #if defined(PETSC_USE_LOG) 341 PetscLogStage stage; 342 #endif 343 char stageName[PETSC_MAX_PATH_LEN]; 344 const char *dmname, *uname; 345 346 ierr = PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);CHKERRQ(ierr); 347 #if defined(PETSC_USE_LOG) 348 ierr = PetscLogStageGetId(stageName, &stage);CHKERRQ(ierr); 349 if (stage < 0) {ierr = PetscLogStageRegister(stageName, &stage);CHKERRQ(ierr);} 350 #endif 351 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 352 if (r > 0) { 353 ierr = DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);CHKERRQ(ierr); 354 ierr = DMSetCoarseDM(dm[r], dm[r-1]);CHKERRQ(ierr); 355 ierr = DMCopyTransform(ce->idm, dm[r]);CHKERRQ(ierr); 356 ierr = PetscObjectGetName((PetscObject) dm[r-1], &dmname);CHKERRQ(ierr); 357 ierr = PetscObjectSetName((PetscObject) dm[r], dmname);CHKERRQ(ierr); 358 for (f = 0; f < ce->Nf; ++f) { 359 PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *); 360 361 ierr = DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);CHKERRQ(ierr); 362 ierr = DMSetNullSpaceConstructor(dm[r], f, nspconstr);CHKERRQ(ierr); 363 } 364 } 365 ierr = DMViewFromOptions(dm[r], NULL, "-conv_dm_view");CHKERRQ(ierr); 366 /* Create solution */ 367 ierr = DMCreateGlobalVector(dm[r], &u);CHKERRQ(ierr); 368 ierr = DMGetField(dm[r], 0, NULL, &disc);CHKERRQ(ierr); 369 ierr = PetscObjectGetName(disc, &uname);CHKERRQ(ierr); 370 ierr = PetscObjectSetName((PetscObject) u, uname);CHKERRQ(ierr); 371 /* Setup solver */ 372 ierr = SNESReset(snes);CHKERRQ(ierr); 373 ierr = SNESSetDM(snes, dm[r]);CHKERRQ(ierr); 374 ierr = DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);CHKERRQ(ierr); 375 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 376 /* Set nullspace for Jacobian */ 377 ierr = PetscConvEstSetJacobianNullspace_Private(ce, snes);CHKERRQ(ierr); 378 /* Create initial guess */ 379 ierr = PetscConvEstComputeInitialGuess(ce, r, dm[r], u);CHKERRQ(ierr); 380 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 381 ierr = PetscLogEventBegin(ce->event, ce, 0, 0, 0);CHKERRQ(ierr); 382 ierr = PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r*ce->Nf]);CHKERRQ(ierr); 383 ierr = PetscLogEventEnd(ce->event, ce, 0, 0, 0);CHKERRQ(ierr); 384 for (f = 0; f < ce->Nf; ++f) { 385 PetscSection s, fs; 386 PetscInt lsize; 387 388 /* Could use DMGetOutputDM() to add in Dirichlet dofs */ 389 ierr = DMGetLocalSection(dm[r], &s);CHKERRQ(ierr); 390 ierr = PetscSectionGetField(s, f, &fs);CHKERRQ(ierr); 391 ierr = PetscSectionGetConstrainedStorageSize(fs, &lsize);CHKERRQ(ierr); 392 ierr = MPI_Allreduce(&lsize, &ce->dofs[r*ce->Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) snes));CHKERRMPI(ierr); 393 ierr = PetscLogEventSetDof(ce->event, f, ce->dofs[r*ce->Nf+f]);CHKERRQ(ierr); 394 ierr = PetscLogEventSetError(ce->event, f, ce->errors[r*ce->Nf+f]);CHKERRQ(ierr); 395 } 396 /* Monitor */ 397 ierr = PetscConvEstMonitorDefault(ce, r);CHKERRQ(ierr); 398 if (!r) { 399 /* PCReset() does not wipe out the level structure */ 400 KSP ksp; 401 PC pc; 402 403 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr); 404 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 405 ierr = PCMGGetLevels(pc, &oldnlev);CHKERRQ(ierr); 406 } 407 /* Cleanup */ 408 ierr = VecDestroy(&u);CHKERRQ(ierr); 409 ierr = PetscLogStagePop();CHKERRQ(ierr); 410 } 411 for (r = 1; r <= Nr; ++r) { 412 ierr = DMDestroy(&dm[r]);CHKERRQ(ierr); 413 } 414 /* Fit convergence rate */ 415 ierr = PetscMalloc2(Nr+1, &x, Nr+1, &y);CHKERRQ(ierr); 416 for (f = 0; f < ce->Nf; ++f) { 417 for (r = 0; r <= Nr; ++r) { 418 x[r] = PetscLog10Real(ce->dofs[r*ce->Nf+f]); 419 y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]); 420 } 421 ierr = PetscLinearRegression(Nr+1, x, y, &slope, &intercept);CHKERRQ(ierr); 422 /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */ 423 alpha[f] = -slope * dim; 424 } 425 ierr = PetscFree2(x, y);CHKERRQ(ierr); 426 ierr = PetscFree(dm);CHKERRQ(ierr); 427 /* Restore solver */ 428 ierr = SNESReset(snes);CHKERRQ(ierr); 429 { 430 /* PCReset() does not wipe out the level structure */ 431 KSP ksp; 432 PC pc; 433 434 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr); 435 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 436 ierr = PCMGSetLevels(pc, oldnlev, NULL);CHKERRQ(ierr); 437 ierr = DMSetRefineLevel(ce->idm, oldlevel);CHKERRQ(ierr); /* The damn DMCoarsen() calls in PCMG can reset this */ 438 } 439 ierr = SNESSetDM(snes, ce->idm);CHKERRQ(ierr); 440 ierr = DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);CHKERRQ(ierr); 441 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 442 ierr = PetscConvEstSetJacobianNullspace_Private(ce, snes);CHKERRQ(ierr); 443 PetscFunctionReturn(0); 444 } 445 446 /*@ 447 PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization 448 449 Not collective 450 451 Input Parameter: 452 . ce - The PetscConvEst object 453 454 Output Parameter: 455 . alpha - The convergence rate for each field 456 457 Note: The convergence rate alpha is defined by 458 $ || u_\Delta - u_exact || < C \Delta^alpha 459 where u_\Delta is the discrete solution, and Delta is a measure of the discretization size. We usually use h for the 460 spatial resolution and \Delta t for the temporal resolution. 461 462 We solve a series of problems using increasing resolution (refined meshes or decreased timesteps), calculate an error 463 based upon the exact solution in the DS, and then fit the result to our model above using linear regression. 464 465 Options database keys: 466 + -snes_convergence_estimate : Execute convergence estimation inside SNESSolve() and print out the rate 467 - -ts_convergence_estimate : Execute convergence estimation inside TSSolve() and print out the rate 468 469 Level: intermediate 470 471 .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate(), SNESSolve(), TSSolve() 472 @*/ 473 PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[]) 474 { 475 PetscInt f; 476 PetscErrorCode ierr; 477 478 PetscFunctionBegin; 479 if (ce->event < 0) {ierr = PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &ce->event);CHKERRQ(ierr);} 480 for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0; 481 ierr = (*ce->ops->getconvrate)(ce, alpha);CHKERRQ(ierr); 482 PetscFunctionReturn(0); 483 } 484 485 /*@ 486 PetscConvEstRateView - Displays the convergence rate to a viewer 487 488 Collective on SNES 489 490 Parameter: 491 + snes - iterative context obtained from SNESCreate() 492 . alpha - the convergence rate for each field 493 - viewer - the viewer to display the reason 494 495 Options Database Keys: 496 . -snes_convergence_estimate - print the convergence rate 497 498 Level: developer 499 500 .seealso: PetscConvEstGetRate() 501 @*/ 502 PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer) 503 { 504 PetscBool isAscii; 505 PetscErrorCode ierr; 506 507 PetscFunctionBegin; 508 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);CHKERRQ(ierr); 509 if (isAscii) { 510 PetscInt Nf = ce->Nf, f; 511 512 ierr = PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);CHKERRQ(ierr); 513 ierr = PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: ");CHKERRQ(ierr); 514 if (Nf > 1) {ierr = PetscViewerASCIIPrintf(viewer, "[");CHKERRQ(ierr);} 515 for (f = 0; f < Nf; ++f) { 516 if (f > 0) {ierr = PetscViewerASCIIPrintf(viewer, ", ");CHKERRQ(ierr);} 517 ierr = PetscViewerASCIIPrintf(viewer, "%#.2g", (double) alpha[f]);CHKERRQ(ierr); 518 } 519 if (Nf > 1) {ierr = PetscViewerASCIIPrintf(viewer, "]");CHKERRQ(ierr);} 520 ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr); 521 ierr = PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);CHKERRQ(ierr); 522 } 523 PetscFunctionReturn(0); 524 } 525 526 /*@ 527 PetscConvEstCreate - Create a PetscConvEst object 528 529 Collective 530 531 Input Parameter: 532 . comm - The communicator for the PetscConvEst object 533 534 Output Parameter: 535 . ce - The PetscConvEst object 536 537 Level: beginner 538 539 .seealso: PetscConvEstDestroy(), PetscConvEstGetConvRate() 540 @*/ 541 PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce) 542 { 543 PetscErrorCode ierr; 544 545 PetscFunctionBegin; 546 PetscValidPointer(ce, 2); 547 ierr = PetscSysInitializePackage();CHKERRQ(ierr); 548 ierr = PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView);CHKERRQ(ierr); 549 (*ce)->monitor = PETSC_FALSE; 550 (*ce)->r = 2.0; 551 (*ce)->Nr = 4; 552 (*ce)->event = -1; 553 (*ce)->ops->setsolver = PetscConvEstSetSNES_Private; 554 (*ce)->ops->initguess = PetscConvEstInitGuessSNES_Private; 555 (*ce)->ops->computeerror = PetscConvEstComputeErrorSNES_Private; 556 (*ce)->ops->getconvrate = PetscConvEstGetConvRateSNES_Private; 557 PetscFunctionReturn(0); 558 } 559