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 PetscReal *X, *Y, H[4], beta[2]; 211 PetscInt i, j, k; 212 PetscErrorCode ierr; 213 214 PetscFunctionBegin; 215 ierr = PetscMalloc2(n*2, &X, n*2, &Y);CHKERRQ(ierr); 216 for (k = 0; k < n; ++k) { 217 /* X[n,2] = [1, x] */ 218 X[k*2+0] = 1.0; 219 X[k*2+1] = x[k]; 220 } 221 /* H = X^T X */ 222 for (i = 0; i < 2; ++i) { 223 for (j = 0; j < 2; ++j) { 224 H[i*2+j] = 0.0; 225 for (k = 0; k < n; ++k) { 226 H[i*2+j] += X[k*2+i] * X[k*2+j]; 227 } 228 } 229 } 230 /* H = (X^T X)^{-1} */ 231 { 232 PetscBLASInt two = 2, ipiv[2], info; 233 PetscReal work[2]; 234 235 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 236 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&two, &two, H, &two, ipiv, &info)); 237 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&two, H, &two, ipiv, work, &two, &info)); 238 ierr = PetscFPTrapPop();CHKERRQ(ierr); 239 } 240 /* Y = H X^T */ 241 for (i = 0; i < 2; ++i) { 242 for (k = 0; k < n; ++k) { 243 Y[i*n+k] = 0.0; 244 for (j = 0; j < 2; ++j) { 245 Y[i*n+k] += H[i*2+j] * X[k*2+j]; 246 } 247 } 248 } 249 /* beta = Y error = [y-intercept, slope] */ 250 for (i = 0; i < 2; ++i) { 251 beta[i] = 0.0; 252 for (k = 0; k < n; ++k) { 253 beta[i] += Y[i*n+k] * y[k]; 254 } 255 } 256 ierr = PetscFree2(X, Y);CHKERRQ(ierr); 257 *intercept = beta[0]; 258 *slope = beta[1]; 259 PetscFunctionReturn(0); 260 } 261 262 /*@ 263 PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization 264 265 Not collective 266 267 Input Parameter: 268 . ce - The PetscConvEst object 269 270 Output Parameter: 271 . alpha - The convergence rate 272 273 Note: The convergence rate alpha is defined by 274 $ || u_h - u_exact || < C h^alpha 275 where u_h is the discrete solution, and h is a measure of the discretization size. 276 277 We solve a series of problems on refined meshes, calculate an error based upon the exact solution in the DS, 278 and then fit the result to our model above using linear regression. 279 280 Options database keys: 281 . -snes_convergence_estimate : Execute convergence estimation and print out the rate 282 283 Level: intermediate 284 285 .keywords: PetscConvEst, convergence 286 .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate() 287 @*/ 288 PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal *alpha) 289 { 290 DM *dm; 291 PetscDS prob; 292 PetscObject disc; 293 MPI_Comm comm; 294 const char *uname, *dmname; 295 void *ctx; 296 Vec u; 297 PetscReal t = 0.0, *x, *y, slope, intercept; 298 PetscInt *dof, dim, Nr = ce->Nr, r; 299 PetscErrorCode ierr; 300 301 PetscFunctionBegin; 302 ierr = PetscObjectGetComm((PetscObject) ce, &comm);CHKERRQ(ierr); 303 ierr = DMGetDimension(ce->idm, &dim);CHKERRQ(ierr); 304 ierr = DMGetApplicationContext(ce->idm, &ctx);CHKERRQ(ierr); 305 ierr = DMGetDS(ce->idm, &prob);CHKERRQ(ierr); 306 ierr = DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);CHKERRQ(ierr); 307 ierr = PetscDSGetDiscretization(prob, 0, &disc);CHKERRQ(ierr); 308 ierr = PetscObjectGetName(disc, &uname);CHKERRQ(ierr); 309 ierr = PetscMalloc2((Nr+1), &dm, (Nr+1), &dof);CHKERRQ(ierr); 310 dm[0] = ce->idm; 311 *alpha = 0.0; 312 /* Loop over meshes */ 313 for (r = 0; r <= Nr; ++r) { 314 if (r > 0) { 315 ierr = DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);CHKERRQ(ierr); 316 ierr = DMSetCoarseDM(dm[r], dm[r-1]);CHKERRQ(ierr); 317 ierr = DMSetDS(dm[r], prob);CHKERRQ(ierr); 318 ierr = PetscObjectGetName((PetscObject) dm[r-1], &dmname);CHKERRQ(ierr); 319 ierr = PetscObjectSetName((PetscObject) dm[r], dmname);CHKERRQ(ierr); 320 } 321 ierr = DMViewFromOptions(dm[r], NULL, "-conv_dm_view");CHKERRQ(ierr); 322 ierr = DMPlexGetHeightStratum(dm[r], 0, NULL, &dof[r]);CHKERRQ(ierr); 323 /* Create solution */ 324 ierr = DMCreateGlobalVector(dm[r], &u);CHKERRQ(ierr); 325 ierr = PetscObjectSetName((PetscObject) u, uname);CHKERRQ(ierr); 326 /* Setup solver */ 327 ierr = SNESReset(ce->snes);CHKERRQ(ierr); 328 ierr = SNESSetDM(ce->snes, dm[r]);CHKERRQ(ierr); 329 ierr = DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);CHKERRQ(ierr); 330 ierr = SNESSetFromOptions(ce->snes);CHKERRQ(ierr); 331 /* Create initial guess */ 332 ierr = DMProjectFunction(dm[r], t, ce->initGuess, NULL, INSERT_VALUES, u);CHKERRQ(ierr); 333 ierr = SNESSolve(ce->snes, NULL, u);CHKERRQ(ierr); 334 ierr = DMComputeL2FieldDiff(dm[r], t, ce->exactSol, NULL, u, &ce->errors[r*ce->Nf]);CHKERRQ(ierr); 335 /* Monitor */ 336 if (ce->monitor) { 337 PetscReal *errors = &ce->errors[r*ce->Nf]; 338 PetscInt f; 339 340 ierr = PetscPrintf(comm, "L_2 Error: [");CHKERRQ(ierr); 341 for (f = 0; f < ce->Nf; ++f) { 342 if (f > 0) {ierr = PetscPrintf(comm, ", ");CHKERRQ(ierr);} 343 if (errors[f] < 1.0e-11) {ierr = PetscPrintf(comm, "< 1e-11");CHKERRQ(ierr);} 344 else {ierr = PetscPrintf(comm, "%g", errors[f]);CHKERRQ(ierr);} 345 } 346 ierr = PetscPrintf(comm, "]\n");CHKERRQ(ierr); 347 } 348 /* Cleanup */ 349 ierr = VecDestroy(&u);CHKERRQ(ierr); 350 } 351 for (r = 1; r <= Nr; ++r) { 352 ierr = DMDestroy(&dm[r]);CHKERRQ(ierr); 353 } 354 /* Fit convergence rate */ 355 ierr = PetscMalloc2(Nr+1, &x, Nr+1, &y);CHKERRQ(ierr); 356 for (r = 0; r <= Nr; ++r) { 357 x[r] = PetscLog10Real(dof[r]); 358 y[r] = PetscLog10Real(ce->errors[r*ce->Nf+0]); 359 } 360 ierr = PetscConvEstLinearRegression_Private(ce, Nr+1, x, y, &slope, &intercept);CHKERRQ(ierr); 361 ierr = PetscFree2(x, y);CHKERRQ(ierr); 362 /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */ 363 *alpha = -slope * dim; 364 ierr = PetscFree2(dm, dof);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 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, 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 ierr = PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);CHKERRQ(ierr); 394 ierr = PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: %0.2g\n", (double) alpha);CHKERRQ(ierr); 395 ierr = PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);CHKERRQ(ierr); 396 } 397 PetscFunctionReturn(0); 398 } 399