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