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