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 /*@ 16 PetscConvEstCreate - Create a PetscConvEst object 17 18 Collective on MPI_Comm 19 20 Input Parameter: 21 . comm - The communicator for the PetscConvEst object 22 23 Output Parameter: 24 . ce - The PetscConvEst object 25 26 Level: beginner 27 28 .keywords: PetscConvEst, convergence, create 29 .seealso: PetscConvEstDestroy(), PetscConvEstGetConvRate() 30 @*/ 31 PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce) 32 { 33 PetscErrorCode ierr; 34 35 PetscFunctionBegin; 36 PetscValidPointer(ce, 2); 37 ierr = PetscSysInitializePackage();CHKERRQ(ierr); 38 ierr = PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView);CHKERRQ(ierr); 39 (*ce)->monitor = PETSC_FALSE; 40 (*ce)->Nr = 4; 41 PetscFunctionReturn(0); 42 } 43 44 /*@ 45 PetscConvEstDestroy - Destroys a PetscConvEst object 46 47 Collective on PetscConvEst 48 49 Input Parameter: 50 . ce - The PetscConvEst object 51 52 Level: beginner 53 54 .keywords: PetscConvEst, convergence, destroy 55 .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate() 56 @*/ 57 PetscErrorCode PetscConvEstDestroy(PetscConvEst *ce) 58 { 59 PetscErrorCode ierr; 60 61 PetscFunctionBegin; 62 if (!*ce) PetscFunctionReturn(0); 63 PetscValidHeaderSpecific((*ce),PETSC_OBJECT_CLASSID,1); 64 if (--((PetscObject)(*ce))->refct > 0) { 65 *ce = NULL; 66 PetscFunctionReturn(0); 67 } 68 ierr = PetscFree3((*ce)->initGuess, (*ce)->exactSol, (*ce)->ctxs);CHKERRQ(ierr); 69 ierr = PetscFree((*ce)->errors);CHKERRQ(ierr); 70 ierr = PetscHeaderDestroy(ce);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 /*@ 75 PetscConvEstSetFromOptions - Sets a PetscConvEst object from options 76 77 Collective on PetscConvEst 78 79 Input Parameters: 80 . ce - The PetscConvEst object 81 82 Level: beginner 83 84 .keywords: PetscConvEst, convergence, options 85 .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate() 86 @*/ 87 PetscErrorCode PetscConvEstSetFromOptions(PetscConvEst ce) 88 { 89 PetscErrorCode ierr; 90 91 PetscFunctionBegin; 92 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject) ce), "", "Convergence Estimator Options", "PetscConvEst");CHKERRQ(ierr); 93 ierr = PetscOptionsInt("-convest_num_refine", "The number of refinements for the convergence check", "PetscConvEst", ce->Nr, &ce->Nr, NULL);CHKERRQ(ierr); 94 ierr = PetscOptionsBool("-convest_monitor", "Monitor the error for each convergence check", "PetscConvEst", ce->monitor, &ce->monitor, NULL);CHKERRQ(ierr); 95 ierr = PetscOptionsEnd(); 96 PetscFunctionReturn(0); 97 } 98 99 /*@ 100 PetscConvEstView - Views a PetscConvEst object 101 102 Collective on PetscConvEst 103 104 Input Parameters: 105 + ce - The PetscConvEst object 106 - viewer - The PetscViewer object 107 108 Level: beginner 109 110 .keywords: PetscConvEst, convergence, view 111 .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate() 112 @*/ 113 PetscErrorCode PetscConvEstView(PetscConvEst ce, PetscViewer viewer) 114 { 115 PetscErrorCode ierr; 116 117 PetscFunctionBegin; 118 ierr = PetscObjectPrintClassNamePrefixType((PetscObject) ce, viewer);CHKERRQ(ierr); 119 ierr = PetscViewerASCIIPrintf(viewer, "ConvEst with %D levels\n", ce->Nr+1);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 123 /*@ 124 PetscConvEstGetSolver - Gets the solver used to produce discrete solutions 125 126 Not collective 127 128 Input Parameter: 129 . ce - The PetscConvEst object 130 131 Output Parameter: 132 . snes - The solver 133 134 Level: intermediate 135 136 .keywords: PetscConvEst, convergence 137 .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate() 138 @*/ 139 PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, SNES *snes) 140 { 141 PetscFunctionBegin; 142 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1); 143 PetscValidPointer(snes, 2); 144 *snes = ce->snes; 145 PetscFunctionReturn(0); 146 } 147 148 /*@ 149 PetscConvEstSetSolver - Sets the solver used to produce discrete solutions 150 151 Not collective 152 153 Input Parameters: 154 + ce - The PetscConvEst object 155 - snes - The solver 156 157 Level: intermediate 158 159 Note: The solver MUST have an attached DM/DS, so that we know the exact solution 160 161 .keywords: PetscConvEst, convergence 162 .seealso: PetscConvEstGetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate() 163 @*/ 164 PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, SNES snes) 165 { 166 PetscErrorCode ierr; 167 168 PetscFunctionBegin; 169 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1); 170 PetscValidHeaderSpecific(snes, SNES_CLASSID, 2); 171 ce->snes = snes; 172 ierr = SNESGetDM(ce->snes, &ce->idm);CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175 176 /*@ 177 PetscConvEstSetUp - After the solver is specified, we create structures for estimating convergence 178 179 Collective on PetscConvEst 180 181 Input Parameters: 182 . ce - The PetscConvEst object 183 184 Level: beginner 185 186 .keywords: PetscConvEst, convergence, setup 187 .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate() 188 @*/ 189 PetscErrorCode PetscConvEstSetUp(PetscConvEst ce) 190 { 191 PetscDS prob; 192 PetscInt f; 193 PetscErrorCode ierr; 194 195 PetscFunctionBegin; 196 ierr = DMGetDS(ce->idm, &prob);CHKERRQ(ierr); 197 ierr = PetscDSGetNumFields(prob, &ce->Nf);CHKERRQ(ierr); 198 ierr = PetscMalloc1((ce->Nr+1)*ce->Nf, &ce->errors);CHKERRQ(ierr); 199 ierr = PetscMalloc3(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol, ce->Nf, &ce->ctxs);CHKERRQ(ierr); 200 for (f = 0; f < ce->Nf; ++f) ce->initGuess[f] = zero_private; 201 for (f = 0; f < ce->Nf; ++f) { 202 ierr = PetscDSGetExactSolution(prob, f, &ce->exactSol[f], &ce->ctxs[f]);CHKERRQ(ierr); 203 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); 204 } 205 PetscFunctionReturn(0); 206 } 207 208 /*@ 209 PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization 210 211 Not collective 212 213 Input Parameter: 214 . ce - The PetscConvEst object 215 216 Output Parameter: 217 . alpha - The convergence rate for each field 218 219 Note: The convergence rate alpha is defined by 220 $ || u_h - u_exact || < C h^alpha 221 where u_h is the discrete solution, and h is a measure of the discretization size. 222 223 We solve a series of problems on refined meshes, calculate an error based upon the exact solution in the DS, 224 and then fit the result to our model above using linear regression. 225 226 Options database keys: 227 . -snes_convergence_estimate : Execute convergence estimation and print out the rate 228 229 Level: intermediate 230 231 .keywords: PetscConvEst, convergence 232 .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate() 233 @*/ 234 PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[]) 235 { 236 DM *dm; 237 PetscObject disc; 238 MPI_Comm comm; 239 const char *uname, *dmname; 240 void *ctx; 241 Vec u; 242 PetscReal t = 0.0, *x, *y, slope, intercept; 243 PetscInt *dof, dim, Nr = ce->Nr, r, f, oldlevel, oldnlev; 244 PetscLogEvent event; 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 ierr = PetscObjectGetComm((PetscObject) ce, &comm);CHKERRQ(ierr); 249 ierr = DMGetDimension(ce->idm, &dim);CHKERRQ(ierr); 250 ierr = DMGetApplicationContext(ce->idm, &ctx);CHKERRQ(ierr); 251 ierr = DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);CHKERRQ(ierr); 252 ierr = DMGetRefineLevel(ce->idm, &oldlevel);CHKERRQ(ierr); 253 ierr = PetscMalloc2((Nr+1), &dm, (Nr+1)*ce->Nf, &dof);CHKERRQ(ierr); 254 dm[0] = ce->idm; 255 for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0; 256 /* Loop over meshes */ 257 ierr = PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &event);CHKERRQ(ierr); 258 for (r = 0; r <= Nr; ++r) { 259 PetscLogStage stage; 260 char stageName[PETSC_MAX_PATH_LEN]; 261 262 ierr = PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);CHKERRQ(ierr); 263 ierr = PetscLogStageRegister(stageName, &stage);CHKERRQ(ierr); 264 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 265 if (r > 0) { 266 ierr = DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);CHKERRQ(ierr); 267 ierr = DMSetCoarseDM(dm[r], dm[r-1]);CHKERRQ(ierr); 268 ierr = DMCopyDisc(ce->idm, dm[r]);CHKERRQ(ierr); 269 ierr = DMCopyTransform(ce->idm, dm[r]);CHKERRQ(ierr); 270 ierr = PetscObjectGetName((PetscObject) dm[r-1], &dmname);CHKERRQ(ierr); 271 ierr = PetscObjectSetName((PetscObject) dm[r], dmname);CHKERRQ(ierr); 272 for (f = 0; f <= ce->Nf; ++f) { 273 PetscErrorCode (*nspconstr)(DM, PetscInt, MatNullSpace *); 274 ierr = DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);CHKERRQ(ierr); 275 ierr = DMSetNullSpaceConstructor(dm[r], f, nspconstr);CHKERRQ(ierr); 276 } 277 } 278 ierr = DMViewFromOptions(dm[r], NULL, "-conv_dm_view");CHKERRQ(ierr); 279 /* Create solution */ 280 ierr = DMCreateGlobalVector(dm[r], &u);CHKERRQ(ierr); 281 ierr = DMGetField(dm[r], 0, NULL, &disc);CHKERRQ(ierr); 282 ierr = PetscObjectGetName(disc, &uname);CHKERRQ(ierr); 283 ierr = PetscObjectSetName((PetscObject) u, uname);CHKERRQ(ierr); 284 /* Setup solver */ 285 ierr = SNESReset(ce->snes);CHKERRQ(ierr); 286 ierr = SNESSetDM(ce->snes, dm[r]);CHKERRQ(ierr); 287 ierr = DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);CHKERRQ(ierr); 288 ierr = SNESSetFromOptions(ce->snes);CHKERRQ(ierr); 289 /* Create initial guess */ 290 ierr = DMProjectFunction(dm[r], t, ce->initGuess, ce->ctxs, INSERT_VALUES, u);CHKERRQ(ierr); 291 ierr = SNESSolve(ce->snes, NULL, u);CHKERRQ(ierr); 292 ierr = PetscLogEventBegin(event, ce, 0, 0, 0);CHKERRQ(ierr); 293 ierr = DMComputeL2FieldDiff(dm[r], t, ce->exactSol, ce->ctxs, u, &ce->errors[r*ce->Nf]);CHKERRQ(ierr); 294 ierr = PetscLogEventEnd(event, ce, 0, 0, 0);CHKERRQ(ierr); 295 for (f = 0; f < ce->Nf; ++f) { 296 PetscSection s, fs; 297 PetscInt lsize; 298 299 /* Could use DMGetOutputDM() to add in Dirichlet dofs */ 300 ierr = DMGetSection(dm[r], &s);CHKERRQ(ierr); 301 ierr = PetscSectionGetField(s, f, &fs);CHKERRQ(ierr); 302 ierr = PetscSectionGetConstrainedStorageSize(fs, &lsize);CHKERRQ(ierr); 303 ierr = MPI_Allreduce(&lsize, &dof[r*ce->Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) ce->snes));CHKERRQ(ierr); 304 ierr = PetscLogEventSetDof(event, f, dof[r*ce->Nf+f]);CHKERRQ(ierr); 305 ierr = PetscLogEventSetError(event, f, ce->errors[r*ce->Nf+f]);CHKERRQ(ierr); 306 } 307 /* Monitor */ 308 if (ce->monitor) { 309 PetscReal *errors = &ce->errors[r*ce->Nf]; 310 311 ierr = PetscPrintf(comm, "L_2 Error: ");CHKERRQ(ierr); 312 if (ce->Nf > 1) {ierr = PetscPrintf(comm, "[");CHKERRQ(ierr);} 313 for (f = 0; f < ce->Nf; ++f) { 314 if (f > 0) {ierr = PetscPrintf(comm, ", ");CHKERRQ(ierr);} 315 if (errors[f] < 1.0e-11) {ierr = PetscPrintf(comm, "< 1e-11");CHKERRQ(ierr);} 316 else {ierr = PetscPrintf(comm, "%g", (double)errors[f]);CHKERRQ(ierr);} 317 } 318 if (ce->Nf > 1) {ierr = PetscPrintf(comm, "]");CHKERRQ(ierr);} 319 ierr = PetscPrintf(comm, "\n");CHKERRQ(ierr); 320 } 321 if (!r) { 322 /* PCReset() does not wipe out the level structure */ 323 KSP ksp; 324 PC pc; 325 326 ierr = SNESGetKSP(ce->snes, &ksp);CHKERRQ(ierr); 327 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 328 ierr = PCMGGetLevels(pc, &oldnlev);CHKERRQ(ierr); 329 } 330 /* Cleanup */ 331 ierr = VecDestroy(&u);CHKERRQ(ierr); 332 ierr = PetscLogStagePop();CHKERRQ(ierr); 333 } 334 for (r = 1; r <= Nr; ++r) { 335 ierr = DMDestroy(&dm[r]);CHKERRQ(ierr); 336 } 337 /* Fit convergence rate */ 338 ierr = PetscMalloc2(Nr+1, &x, Nr+1, &y);CHKERRQ(ierr); 339 for (f = 0; f < ce->Nf; ++f) { 340 for (r = 0; r <= Nr; ++r) { 341 x[r] = PetscLog10Real(dof[r*ce->Nf+f]); 342 y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]); 343 } 344 ierr = PetscLinearRegression(Nr+1, x, y, &slope, &intercept);CHKERRQ(ierr); 345 /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */ 346 alpha[f] = -slope * dim; 347 } 348 ierr = PetscFree2(x, y);CHKERRQ(ierr); 349 ierr = PetscFree2(dm, dof);CHKERRQ(ierr); 350 /* Restore solver */ 351 ierr = SNESReset(ce->snes);CHKERRQ(ierr); 352 { 353 /* PCReset() does not wipe out the level structure */ 354 KSP ksp; 355 PC pc; 356 357 ierr = SNESGetKSP(ce->snes, &ksp);CHKERRQ(ierr); 358 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 359 ierr = PCMGSetLevels(pc, oldnlev, NULL);CHKERRQ(ierr); 360 ierr = DMSetRefineLevel(ce->idm, oldlevel);CHKERRQ(ierr); /* The damn DMCoarsen() calls in PCMG can reset this */ 361 } 362 ierr = SNESSetDM(ce->snes, ce->idm);CHKERRQ(ierr); 363 ierr = DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);CHKERRQ(ierr); 364 ierr = SNESSetFromOptions(ce->snes);CHKERRQ(ierr); 365 PetscFunctionReturn(0); 366 } 367 368 /*@ 369 PetscConvEstRateView - Displays the convergence rate to a viewer 370 371 Collective on SNES 372 373 Parameter: 374 + snes - iterative context obtained from SNESCreate() 375 . alpha - the convergence rate for each field 376 - viewer - the viewer to display the reason 377 378 Options Database Keys: 379 . -snes_convergence_estimate - print the convergence rate 380 381 Level: developer 382 383 .seealso: PetscConvEstGetRate() 384 @*/ 385 PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer) 386 { 387 PetscBool isAscii; 388 PetscErrorCode ierr; 389 390 PetscFunctionBegin; 391 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);CHKERRQ(ierr); 392 if (isAscii) { 393 DM dm; 394 PetscInt Nf, f; 395 396 ierr = SNESGetDM(ce->snes, &dm);CHKERRQ(ierr); 397 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 398 ierr = PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);CHKERRQ(ierr); 399 ierr = PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: ");CHKERRQ(ierr); 400 if (Nf > 1) {ierr = PetscViewerASCIIPrintf(viewer, "[");CHKERRQ(ierr);} 401 for (f = 0; f < Nf; ++f) { 402 if (f > 0) {ierr = PetscViewerASCIIPrintf(viewer, ", ");CHKERRQ(ierr);} 403 ierr = PetscViewerASCIIPrintf(viewer, "%#.2g", (double) alpha[f]);CHKERRQ(ierr); 404 } 405 if (Nf > 1) {ierr = PetscViewerASCIIPrintf(viewer, "]");CHKERRQ(ierr);} 406 ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr); 407 ierr = PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);CHKERRQ(ierr); 408 } 409 PetscFunctionReturn(0); 410 } 411