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("-num_refine", "The number of refinements for the convergence check", "PetscConvEst", ce->Nr, &ce->Nr, 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 = PetscMalloc2(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol);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]);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 static PetscErrorCode PetscConvEstLinearRegression_Private(PetscConvEst ce, PetscInt n, const PetscReal x[], const PetscReal y[], PetscReal *slope, PetscReal *intercept) 209 { 210 PetscScalar H[4]; 211 PetscReal *X, *Y, beta[2]; 212 PetscInt i, j, k; 213 PetscErrorCode ierr; 214 215 PetscFunctionBegin; 216 *slope = *intercept = 0.0; 217 ierr = PetscMalloc2(n*2, &X, n*2, &Y);CHKERRQ(ierr); 218 for (k = 0; k < n; ++k) { 219 /* X[n,2] = [1, x] */ 220 X[k*2+0] = 1.0; 221 X[k*2+1] = x[k]; 222 } 223 /* H = X^T X */ 224 for (i = 0; i < 2; ++i) { 225 for (j = 0; j < 2; ++j) { 226 H[i*2+j] = 0.0; 227 for (k = 0; k < n; ++k) { 228 H[i*2+j] += X[k*2+i] * X[k*2+j]; 229 } 230 } 231 } 232 /* H = (X^T X)^{-1} */ 233 { 234 PetscBLASInt two = 2, ipiv[2], info; 235 PetscScalar work[2]; 236 237 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 238 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&two, &two, H, &two, ipiv, &info)); 239 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&two, H, &two, ipiv, work, &two, &info)); 240 ierr = PetscFPTrapPop();CHKERRQ(ierr); 241 } 242 /* Y = H X^T */ 243 for (i = 0; i < 2; ++i) { 244 for (k = 0; k < n; ++k) { 245 Y[i*n+k] = 0.0; 246 for (j = 0; j < 2; ++j) { 247 Y[i*n+k] += PetscRealPart(H[i*2+j]) * X[k*2+j]; 248 } 249 } 250 } 251 /* beta = Y error = [y-intercept, slope] */ 252 for (i = 0; i < 2; ++i) { 253 beta[i] = 0.0; 254 for (k = 0; k < n; ++k) { 255 beta[i] += Y[i*n+k] * y[k]; 256 } 257 } 258 ierr = PetscFree2(X, Y);CHKERRQ(ierr); 259 *intercept = beta[0]; 260 *slope = beta[1]; 261 PetscFunctionReturn(0); 262 } 263 264 /*@ 265 PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization 266 267 Not collective 268 269 Input Parameter: 270 . ce - The PetscConvEst object 271 272 Output Parameter: 273 . alpha - The convergence rate for each field 274 275 Note: The convergence rate alpha is defined by 276 $ || u_h - u_exact || < C h^alpha 277 where u_h is the discrete solution, and h is a measure of the discretization size. 278 279 We solve a series of problems on refined meshes, calculate an error based upon the exact solution in the DS, 280 and then fit the result to our model above using linear regression. 281 282 Options database keys: 283 . -snes_convergence_estimate : Execute convergence estimation and print out the rate 284 285 Level: intermediate 286 287 .keywords: PetscConvEst, convergence 288 .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate() 289 @*/ 290 PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal *alpha) 291 { 292 DM *dm; 293 PetscDS prob; 294 PetscObject disc; 295 MPI_Comm comm; 296 const char *uname, *dmname; 297 void *ctx; 298 Vec u; 299 PetscReal t = 0.0, *x, *y, slope, intercept; 300 PetscInt *dof, dim, Nf, f, Nr = ce->Nr, r; 301 PetscLogEvent event; 302 PetscErrorCode ierr; 303 304 PetscFunctionBegin; 305 ierr = PetscObjectGetComm((PetscObject) ce, &comm);CHKERRQ(ierr); 306 ierr = DMGetDimension(ce->idm, &dim);CHKERRQ(ierr); 307 ierr = DMGetApplicationContext(ce->idm, &ctx);CHKERRQ(ierr); 308 ierr = DMGetDS(ce->idm, &prob);CHKERRQ(ierr); 309 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 310 ierr = DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);CHKERRQ(ierr); 311 ierr = PetscMalloc2((Nr+1), &dm, (Nr+1)*Nf, &dof);CHKERRQ(ierr); 312 dm[0] = ce->idm; 313 *alpha = 0.0; 314 /* Loop over meshes */ 315 ierr = PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &event);CHKERRQ(ierr); 316 for (r = 0; r <= Nr; ++r) { 317 PetscLogStage stage; 318 char stageName[PETSC_MAX_PATH_LEN]; 319 PetscInt f; 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 } 331 ierr = DMViewFromOptions(dm[r], NULL, "-conv_dm_view");CHKERRQ(ierr); 332 /* Create solution */ 333 ierr = DMCreateGlobalVector(dm[r], &u);CHKERRQ(ierr); 334 ierr = PetscDSGetDiscretization(prob, 0, &disc);CHKERRQ(ierr); 335 ierr = PetscObjectGetName(disc, &uname);CHKERRQ(ierr); 336 ierr = PetscObjectSetName((PetscObject) u, uname);CHKERRQ(ierr); 337 /* Setup solver */ 338 ierr = SNESReset(ce->snes);CHKERRQ(ierr); 339 ierr = SNESSetDM(ce->snes, dm[r]);CHKERRQ(ierr); 340 ierr = DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);CHKERRQ(ierr); 341 ierr = SNESSetFromOptions(ce->snes);CHKERRQ(ierr); 342 /* Create initial guess */ 343 ierr = DMProjectFunction(dm[r], t, ce->initGuess, NULL, INSERT_VALUES, u);CHKERRQ(ierr); 344 ierr = SNESSolve(ce->snes, NULL, u);CHKERRQ(ierr); 345 ierr = PetscLogEventBegin(event, ce, 0, 0, 0);CHKERRQ(ierr); 346 ierr = DMComputeL2FieldDiff(dm[r], t, ce->exactSol, NULL, u, &ce->errors[r*ce->Nf]);CHKERRQ(ierr); 347 ierr = PetscLogEventEnd(event, ce, 0, 0, 0);CHKERRQ(ierr); 348 for (f = 0; f < ce->Nf; ++f) { 349 PetscSection s, fs; 350 PetscInt lsize; 351 352 /* Could use DMGetOutputDM() to add in Dirichlet dofs */ 353 ierr = DMGetSection(dm[r], &s);CHKERRQ(ierr); 354 ierr = PetscSectionGetField(s, f, &fs);CHKERRQ(ierr); 355 ierr = PetscSectionGetConstrainedStorageSize(fs, &lsize);CHKERRQ(ierr); 356 ierr = MPI_Allreduce(&lsize, &dof[r*ce->Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) ce->snes));CHKERRQ(ierr); 357 ierr = PetscLogEventSetDof(event, f, dof[r*ce->Nf+f]);CHKERRQ(ierr); 358 ierr = PetscLogEventSetError(event, f, ce->errors[r*ce->Nf+f]);CHKERRQ(ierr); 359 } 360 /* Monitor */ 361 if (ce->monitor) { 362 PetscReal *errors = &ce->errors[r*ce->Nf]; 363 364 ierr = PetscPrintf(comm, "L_2 Error: [");CHKERRQ(ierr); 365 for (f = 0; f < ce->Nf; ++f) { 366 if (f > 0) {ierr = PetscPrintf(comm, ", ");CHKERRQ(ierr);} 367 if (errors[f] < 1.0e-11) {ierr = PetscPrintf(comm, "< 1e-11");CHKERRQ(ierr);} 368 else {ierr = PetscPrintf(comm, "%g", (double)errors[f]);CHKERRQ(ierr);} 369 } 370 ierr = PetscPrintf(comm, "]\n");CHKERRQ(ierr); 371 } 372 /* Cleanup */ 373 ierr = VecDestroy(&u);CHKERRQ(ierr); 374 ierr = PetscLogStagePop();CHKERRQ(ierr); 375 } 376 for (r = 1; r <= Nr; ++r) { 377 ierr = DMDestroy(&dm[r]);CHKERRQ(ierr); 378 } 379 /* Fit convergence rate */ 380 ierr = PetscMalloc2(Nr+1, &x, Nr+1, &y);CHKERRQ(ierr); 381 for (f = 0; f < ce->Nf; ++f) { 382 for (r = 0; r <= Nr; ++r) { 383 x[r] = PetscLog10Real(dof[r*ce->Nf+f]); 384 y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]); 385 } 386 ierr = PetscConvEstLinearRegression_Private(ce, Nr+1, x, y, &slope, &intercept);CHKERRQ(ierr); 387 /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */ 388 alpha[f] = -slope * dim; 389 } 390 ierr = PetscFree2(x, y);CHKERRQ(ierr); 391 ierr = PetscFree2(dm, dof);CHKERRQ(ierr); 392 /* Restore solver */ 393 ierr = SNESReset(ce->snes);CHKERRQ(ierr); 394 ierr = SNESSetDM(ce->snes, ce->idm);CHKERRQ(ierr); 395 ierr = DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);CHKERRQ(ierr); 396 ierr = SNESSetFromOptions(ce->snes);CHKERRQ(ierr); 397 PetscFunctionReturn(0); 398 } 399 400 /*@ 401 PetscConvEstRateView - Displays the convergence rate to a viewer 402 403 Collective on SNES 404 405 Parameter: 406 + snes - iterative context obtained from SNESCreate() 407 . alpha - the convergence rate for each field 408 - viewer - the viewer to display the reason 409 410 Options Database Keys: 411 . -snes_convergence_estimate - print the convergence rate 412 413 Level: developer 414 415 .seealso: PetscConvEstGetRate() 416 @*/ 417 PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer) 418 { 419 PetscBool isAscii; 420 PetscErrorCode ierr; 421 422 PetscFunctionBegin; 423 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);CHKERRQ(ierr); 424 if (isAscii) { 425 DM dm; 426 PetscInt Nf, f; 427 428 ierr = SNESGetDM(ce->snes, &dm);CHKERRQ(ierr); 429 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 430 ierr = PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);CHKERRQ(ierr); 431 ierr = PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: ");CHKERRQ(ierr); 432 if (Nf > 1) {ierr = PetscViewerASCIIPrintf(viewer, "[");CHKERRQ(ierr);} 433 for (f = 0; f < Nf; ++f) { 434 if (f > 0) {ierr = PetscViewerASCIIPrintf(viewer, ", ");CHKERRQ(ierr);} 435 ierr = PetscViewerASCIIPrintf(viewer, "%.2g", (double) alpha[f]);CHKERRQ(ierr); 436 } 437 if (Nf > 1) {ierr = PetscViewerASCIIPrintf(viewer, "]");CHKERRQ(ierr);} 438 ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr); 439 ierr = PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);CHKERRQ(ierr); 440 } 441 PetscFunctionReturn(0); 442 } 443