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(); 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 ierr = PetscLogStageRegister(stageName, &stage);CHKERRQ(ierr); 348 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 349 if (r > 0) { 350 ierr = DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);CHKERRQ(ierr); 351 ierr = DMSetCoarseDM(dm[r], dm[r-1]);CHKERRQ(ierr); 352 ierr = DMCopyTransform(ce->idm, dm[r]);CHKERRQ(ierr); 353 ierr = PetscObjectGetName((PetscObject) dm[r-1], &dmname);CHKERRQ(ierr); 354 ierr = PetscObjectSetName((PetscObject) dm[r], dmname);CHKERRQ(ierr); 355 for (f = 0; f < ce->Nf; ++f) { 356 PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *); 357 358 ierr = DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);CHKERRQ(ierr); 359 ierr = DMSetNullSpaceConstructor(dm[r], f, nspconstr);CHKERRQ(ierr); 360 } 361 } 362 ierr = DMViewFromOptions(dm[r], NULL, "-conv_dm_view");CHKERRQ(ierr); 363 /* Create solution */ 364 ierr = DMCreateGlobalVector(dm[r], &u);CHKERRQ(ierr); 365 ierr = DMGetField(dm[r], 0, NULL, &disc);CHKERRQ(ierr); 366 ierr = PetscObjectGetName(disc, &uname);CHKERRQ(ierr); 367 ierr = PetscObjectSetName((PetscObject) u, uname);CHKERRQ(ierr); 368 /* Setup solver */ 369 ierr = SNESReset(snes);CHKERRQ(ierr); 370 ierr = SNESSetDM(snes, dm[r]);CHKERRQ(ierr); 371 ierr = DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);CHKERRQ(ierr); 372 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 373 /* Set nullspace for Jacobian */ 374 ierr = PetscConvEstSetJacobianNullspace_Private(ce, snes);CHKERRQ(ierr); 375 /* Create initial guess */ 376 ierr = PetscConvEstComputeInitialGuess(ce, r, dm[r], u);CHKERRQ(ierr); 377 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 378 ierr = PetscLogEventBegin(ce->event, ce, 0, 0, 0);CHKERRQ(ierr); 379 ierr = PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r*ce->Nf]);CHKERRQ(ierr); 380 ierr = PetscLogEventEnd(ce->event, ce, 0, 0, 0);CHKERRQ(ierr); 381 for (f = 0; f < ce->Nf; ++f) { 382 PetscSection s, fs; 383 PetscInt lsize; 384 385 /* Could use DMGetOutputDM() to add in Dirichlet dofs */ 386 ierr = DMGetLocalSection(dm[r], &s);CHKERRQ(ierr); 387 ierr = PetscSectionGetField(s, f, &fs);CHKERRQ(ierr); 388 ierr = PetscSectionGetConstrainedStorageSize(fs, &lsize);CHKERRQ(ierr); 389 ierr = MPI_Allreduce(&lsize, &ce->dofs[r*ce->Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) snes));CHKERRMPI(ierr); 390 ierr = PetscLogEventSetDof(ce->event, f, ce->dofs[r*ce->Nf+f]);CHKERRQ(ierr); 391 ierr = PetscLogEventSetError(ce->event, f, ce->errors[r*ce->Nf+f]);CHKERRQ(ierr); 392 } 393 /* Monitor */ 394 ierr = PetscConvEstMonitorDefault(ce, r);CHKERRQ(ierr); 395 if (!r) { 396 /* PCReset() does not wipe out the level structure */ 397 KSP ksp; 398 PC pc; 399 400 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr); 401 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 402 ierr = PCMGGetLevels(pc, &oldnlev);CHKERRQ(ierr); 403 } 404 /* Cleanup */ 405 ierr = VecDestroy(&u);CHKERRQ(ierr); 406 ierr = PetscLogStagePop();CHKERRQ(ierr); 407 } 408 for (r = 1; r <= Nr; ++r) { 409 ierr = DMDestroy(&dm[r]);CHKERRQ(ierr); 410 } 411 /* Fit convergence rate */ 412 ierr = PetscMalloc2(Nr+1, &x, Nr+1, &y);CHKERRQ(ierr); 413 for (f = 0; f < ce->Nf; ++f) { 414 for (r = 0; r <= Nr; ++r) { 415 x[r] = PetscLog10Real(ce->dofs[r*ce->Nf+f]); 416 y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]); 417 } 418 ierr = PetscLinearRegression(Nr+1, x, y, &slope, &intercept);CHKERRQ(ierr); 419 /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */ 420 alpha[f] = -slope * dim; 421 } 422 ierr = PetscFree2(x, y);CHKERRQ(ierr); 423 ierr = PetscFree(dm);CHKERRQ(ierr); 424 /* Restore solver */ 425 ierr = SNESReset(snes);CHKERRQ(ierr); 426 { 427 /* PCReset() does not wipe out the level structure */ 428 KSP ksp; 429 PC pc; 430 431 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr); 432 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 433 ierr = PCMGSetLevels(pc, oldnlev, NULL);CHKERRQ(ierr); 434 ierr = DMSetRefineLevel(ce->idm, oldlevel);CHKERRQ(ierr); /* The damn DMCoarsen() calls in PCMG can reset this */ 435 } 436 ierr = SNESSetDM(snes, ce->idm);CHKERRQ(ierr); 437 ierr = DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);CHKERRQ(ierr); 438 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 439 ierr = PetscConvEstSetJacobianNullspace_Private(ce, snes);CHKERRQ(ierr); 440 PetscFunctionReturn(0); 441 } 442 443 /*@ 444 PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization 445 446 Not collective 447 448 Input Parameter: 449 . ce - The PetscConvEst object 450 451 Output Parameter: 452 . alpha - The convergence rate for each field 453 454 Note: The convergence rate alpha is defined by 455 $ || u_\Delta - u_exact || < C \Delta^alpha 456 where u_\Delta is the discrete solution, and Delta is a measure of the discretization size. We usually use h for the 457 spatial resolution and \Delta t for the temporal resolution. 458 459 We solve a series of problems using increasing resolution (refined meshes or decreased timesteps), calculate an error 460 based upon the exact solution in the DS, and then fit the result to our model above using linear regression. 461 462 Options database keys: 463 + -snes_convergence_estimate : Execute convergence estimation inside SNESSolve() and print out the rate 464 - -ts_convergence_estimate : Execute convergence estimation inside TSSolve() and print out the rate 465 466 Level: intermediate 467 468 .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate(), SNESSolve(), TSSolve() 469 @*/ 470 PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[]) 471 { 472 PetscInt f; 473 PetscErrorCode ierr; 474 475 PetscFunctionBegin; 476 if (ce->event < 0) {ierr = PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &ce->event);CHKERRQ(ierr);} 477 for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0; 478 ierr = (*ce->ops->getconvrate)(ce, alpha);CHKERRQ(ierr); 479 PetscFunctionReturn(0); 480 } 481 482 /*@ 483 PetscConvEstRateView - Displays the convergence rate to a viewer 484 485 Collective on SNES 486 487 Parameter: 488 + snes - iterative context obtained from SNESCreate() 489 . alpha - the convergence rate for each field 490 - viewer - the viewer to display the reason 491 492 Options Database Keys: 493 . -snes_convergence_estimate - print the convergence rate 494 495 Level: developer 496 497 .seealso: PetscConvEstGetRate() 498 @*/ 499 PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer) 500 { 501 PetscBool isAscii; 502 PetscErrorCode ierr; 503 504 PetscFunctionBegin; 505 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);CHKERRQ(ierr); 506 if (isAscii) { 507 PetscInt Nf = ce->Nf, f; 508 509 ierr = PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);CHKERRQ(ierr); 510 ierr = PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: ");CHKERRQ(ierr); 511 if (Nf > 1) {ierr = PetscViewerASCIIPrintf(viewer, "[");CHKERRQ(ierr);} 512 for (f = 0; f < Nf; ++f) { 513 if (f > 0) {ierr = PetscViewerASCIIPrintf(viewer, ", ");CHKERRQ(ierr);} 514 ierr = PetscViewerASCIIPrintf(viewer, "%#.2g", (double) alpha[f]);CHKERRQ(ierr); 515 } 516 if (Nf > 1) {ierr = PetscViewerASCIIPrintf(viewer, "]");CHKERRQ(ierr);} 517 ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr); 518 ierr = PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);CHKERRQ(ierr); 519 } 520 PetscFunctionReturn(0); 521 } 522 523 /*@ 524 PetscConvEstCreate - Create a PetscConvEst object 525 526 Collective 527 528 Input Parameter: 529 . comm - The communicator for the PetscConvEst object 530 531 Output Parameter: 532 . ce - The PetscConvEst object 533 534 Level: beginner 535 536 .seealso: PetscConvEstDestroy(), PetscConvEstGetConvRate() 537 @*/ 538 PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce) 539 { 540 PetscErrorCode ierr; 541 542 PetscFunctionBegin; 543 PetscValidPointer(ce, 2); 544 ierr = PetscSysInitializePackage();CHKERRQ(ierr); 545 ierr = PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView);CHKERRQ(ierr); 546 (*ce)->monitor = PETSC_FALSE; 547 (*ce)->r = 2.0; 548 (*ce)->Nr = 4; 549 (*ce)->event = -1; 550 (*ce)->ops->setsolver = PetscConvEstSetSNES_Private; 551 (*ce)->ops->initguess = PetscConvEstInitGuessSNES_Private; 552 (*ce)->ops->computeerror = PetscConvEstComputeErrorSNES_Private; 553 (*ce)->ops->getconvrate = PetscConvEstGetConvRateSNES_Private; 554 PetscFunctionReturn(0); 555 } 556