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