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