1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 3 #include <petscds.h> 4 #include <petsc/private/petscimpl.h> 5 #include <petsc/private/petscfeimpl.h> 6 7 #ifdef PETSC_HAVE_LIBCEED 8 #include <petscdmceed.h> 9 #include <petscdmplexceed.h> 10 #endif 11 12 static void pressure_Private(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar p[]) 13 { 14 p[0] = u[uOff[1]]; 15 } 16 17 /* 18 SNESCorrectDiscretePressure_Private - Add a vector in the nullspace to make the continuum integral of the pressure field equal to zero. 19 This is normally used only to evaluate convergence rates for the pressure accurately. 20 21 Collective 22 23 Input Parameters: 24 + snes - The `SNES` 25 . pfield - The field number for pressure 26 . nullspace - The pressure nullspace 27 . u - The solution vector 28 - ctx - An optional user context 29 30 Output Parameter: 31 . u - The solution with a continuum pressure integral of zero 32 33 Level: developer 34 35 Note: 36 If int(u) = a and int(n) = b, then int(u - a/b n) = a - a/b b = 0. We assume that the nullspace is a single vector given explicitly. 37 38 .seealso: [](ch_snes), `SNESConvergedCorrectPressure()` 39 */ 40 static PetscErrorCode SNESCorrectDiscretePressure_Private(SNES snes, PetscInt pfield, MatNullSpace nullspace, Vec u, void *ctx) 41 { 42 DM dm; 43 PetscDS ds; 44 const Vec *nullvecs; 45 PetscScalar pintd, *intc, *intn; 46 MPI_Comm comm; 47 PetscInt Nf, Nv; 48 49 PetscFunctionBegin; 50 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 51 PetscCall(SNESGetDM(snes, &dm)); 52 PetscCheck(dm, comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a SNES DM"); 53 PetscCheck(nullspace, comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a Jacobian nullspace"); 54 PetscCall(DMGetDS(dm, &ds)); 55 PetscCall(PetscDSSetObjective(ds, pfield, pressure_Private)); 56 PetscCall(MatNullSpaceGetVecs(nullspace, NULL, &Nv, &nullvecs)); 57 PetscCheck(Nv == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Can only handle a single null vector for pressure, not %" PetscInt_FMT, Nv); 58 PetscCall(VecDot(nullvecs[0], u, &pintd)); 59 PetscCheck(PetscAbsScalar(pintd) <= PETSC_SMALL, comm, PETSC_ERR_ARG_WRONG, "Discrete integral of pressure: %g", (double)PetscRealPart(pintd)); 60 PetscCall(PetscDSGetNumFields(ds, &Nf)); 61 PetscCall(PetscMalloc2(Nf, &intc, Nf, &intn)); 62 PetscCall(DMPlexComputeIntegralFEM(dm, nullvecs[0], intn, ctx)); 63 PetscCall(DMPlexComputeIntegralFEM(dm, u, intc, ctx)); 64 PetscCall(VecAXPY(u, -intc[pfield] / intn[pfield], nullvecs[0])); 65 #if defined(PETSC_USE_DEBUG) 66 PetscCall(DMPlexComputeIntegralFEM(dm, u, intc, ctx)); 67 PetscCheck(PetscAbsScalar(intc[pfield]) <= PETSC_SMALL, comm, PETSC_ERR_ARG_WRONG, "Continuum integral of pressure after correction: %g", (double)PetscRealPart(intc[pfield])); 68 #endif 69 PetscCall(PetscFree2(intc, intn)); 70 PetscFunctionReturn(PETSC_SUCCESS); 71 } 72 73 /*@C 74 SNESConvergedCorrectPressure - The regular `SNES` convergence test that, up on convergence, adds a vector in the nullspace 75 to make the continuum integral of the pressure field equal to zero. 76 77 Logically Collective 78 79 Input Parameters: 80 + snes - the `SNES` context 81 . it - the iteration (0 indicates before any Newton steps) 82 . xnorm - 2-norm of current iterate 83 . gnorm - 2-norm of current step 84 . f - 2-norm of function at current iterate 85 - ctx - Optional user context 86 87 Output Parameter: 88 . reason - `SNES_CONVERGED_ITERATING`, `SNES_CONVERGED_ITS`, or `SNES_DIVERGED_FNORM_NAN` 89 90 Options Database Key: 91 . -snes_convergence_test correct_pressure - see `SNESSetFromOptions()` 92 93 Level: advanced 94 95 Notes: 96 In order to use this convergence test, you must set up several PETSc structures. First fields must be added to the `DM`, and a `PetscDS` 97 must be created with discretizations of those fields. We currently assume that the pressure field has index 1. 98 The pressure field must have a nullspace, likely created using the `DMSetNullSpaceConstructor()` interface. 99 Last we must be able to integrate the pressure over the domain, so the `DM` attached to the SNES `must` be a `DMPLEX` at this time. 100 101 Developer Note: 102 This is a total misuse of the `SNES` convergence test handling system. It should be removed. Perhaps a `SNESSetPostSolve()` could 103 be constructed to handle this process. 104 105 .seealso: [](ch_snes), `SNES`, `DM`, `SNESConvergedDefault()`, `SNESSetConvergenceTest()`, `DMSetNullSpaceConstructor()` 106 @*/ 107 PetscErrorCode SNESConvergedCorrectPressure(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal f, SNESConvergedReason *reason, void *ctx) 108 { 109 PetscBool monitorIntegral = PETSC_FALSE; 110 111 PetscFunctionBegin; 112 PetscCall(SNESConvergedDefault(snes, it, xnorm, gnorm, f, reason, ctx)); 113 if (monitorIntegral) { 114 Mat J; 115 Vec u; 116 MatNullSpace nullspace; 117 const Vec *nullvecs; 118 PetscScalar pintd; 119 120 PetscCall(SNESGetSolution(snes, &u)); 121 PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL)); 122 PetscCall(MatGetNullSpace(J, &nullspace)); 123 PetscCall(MatNullSpaceGetVecs(nullspace, NULL, NULL, &nullvecs)); 124 PetscCall(VecDot(nullvecs[0], u, &pintd)); 125 PetscCall(PetscInfo(snes, "SNES: Discrete integral of pressure: %g\n", (double)PetscRealPart(pintd))); 126 } 127 if (*reason > 0) { 128 Mat J; 129 Vec u; 130 MatNullSpace nullspace; 131 PetscInt pfield = 1; 132 133 PetscCall(SNESGetSolution(snes, &u)); 134 PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL)); 135 PetscCall(MatGetNullSpace(J, &nullspace)); 136 PetscCall(SNESCorrectDiscretePressure_Private(snes, pfield, nullspace, u, ctx)); 137 } 138 PetscFunctionReturn(PETSC_SUCCESS); 139 } 140 141 /************************** Interpolation *******************************/ 142 143 static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy) 144 { 145 PetscBool isPlex; 146 147 PetscFunctionBegin; 148 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex)); 149 if (isPlex) { 150 *plex = dm; 151 PetscCall(PetscObjectReference((PetscObject)dm)); 152 } else { 153 PetscCall(PetscObjectQuery((PetscObject)dm, "dm_plex", (PetscObject *)plex)); 154 if (!*plex) { 155 PetscCall(DMConvert(dm, DMPLEX, plex)); 156 PetscCall(PetscObjectCompose((PetscObject)dm, "dm_plex", (PetscObject)*plex)); 157 if (copy) { 158 PetscCall(DMCopyDMSNES(dm, *plex)); 159 PetscCall(DMCopyAuxiliaryVec(dm, *plex)); 160 } 161 } else { 162 PetscCall(PetscObjectReference((PetscObject)*plex)); 163 } 164 } 165 PetscFunctionReturn(PETSC_SUCCESS); 166 } 167 168 /*@C 169 DMInterpolationCreate - Creates a `DMInterpolationInfo` context 170 171 Collective 172 173 Input Parameter: 174 . comm - the communicator 175 176 Output Parameter: 177 . ctx - the context 178 179 Level: beginner 180 181 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationDestroy()` 182 @*/ 183 PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx) 184 { 185 PetscFunctionBegin; 186 PetscAssertPointer(ctx, 2); 187 PetscCall(PetscNew(ctx)); 188 189 (*ctx)->comm = comm; 190 (*ctx)->dim = -1; 191 (*ctx)->nInput = 0; 192 (*ctx)->points = NULL; 193 (*ctx)->cells = NULL; 194 (*ctx)->n = -1; 195 (*ctx)->coords = NULL; 196 PetscFunctionReturn(PETSC_SUCCESS); 197 } 198 199 /*@C 200 DMInterpolationSetDim - Sets the spatial dimension for the interpolation context 201 202 Not Collective 203 204 Input Parameters: 205 + ctx - the context 206 - dim - the spatial dimension 207 208 Level: intermediate 209 210 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()` 211 @*/ 212 PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim) 213 { 214 PetscFunctionBegin; 215 PetscCheck(!(dim < 1) && !(dim > 3), ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %" PetscInt_FMT, dim); 216 ctx->dim = dim; 217 PetscFunctionReturn(PETSC_SUCCESS); 218 } 219 220 /*@C 221 DMInterpolationGetDim - Gets the spatial dimension for the interpolation context 222 223 Not Collective 224 225 Input Parameter: 226 . ctx - the context 227 228 Output Parameter: 229 . dim - the spatial dimension 230 231 Level: intermediate 232 233 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()` 234 @*/ 235 PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim) 236 { 237 PetscFunctionBegin; 238 PetscAssertPointer(dim, 2); 239 *dim = ctx->dim; 240 PetscFunctionReturn(PETSC_SUCCESS); 241 } 242 243 /*@C 244 DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context 245 246 Not Collective 247 248 Input Parameters: 249 + ctx - the context 250 - dof - the number of fields 251 252 Level: intermediate 253 254 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()` 255 @*/ 256 PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof) 257 { 258 PetscFunctionBegin; 259 PetscCheck(dof >= 1, ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %" PetscInt_FMT, dof); 260 ctx->dof = dof; 261 PetscFunctionReturn(PETSC_SUCCESS); 262 } 263 264 /*@C 265 DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context 266 267 Not Collective 268 269 Input Parameter: 270 . ctx - the context 271 272 Output Parameter: 273 . dof - the number of fields 274 275 Level: intermediate 276 277 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()` 278 @*/ 279 PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof) 280 { 281 PetscFunctionBegin; 282 PetscAssertPointer(dof, 2); 283 *dof = ctx->dof; 284 PetscFunctionReturn(PETSC_SUCCESS); 285 } 286 287 /*@C 288 DMInterpolationAddPoints - Add points at which we will interpolate the fields 289 290 Not Collective 291 292 Input Parameters: 293 + ctx - the context 294 . n - the number of points 295 - points - the coordinates for each point, an array of size n * dim 296 297 Level: intermediate 298 299 Note: 300 The coordinate information is copied. 301 302 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationCreate()` 303 @*/ 304 PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[]) 305 { 306 PetscFunctionBegin; 307 PetscCheck(ctx->dim >= 0, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set"); 308 PetscCheck(!ctx->points, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet"); 309 ctx->nInput = n; 310 311 PetscCall(PetscMalloc1(n * ctx->dim, &ctx->points)); 312 PetscCall(PetscArraycpy(ctx->points, points, n * ctx->dim)); 313 PetscFunctionReturn(PETSC_SUCCESS); 314 } 315 316 /*@C 317 DMInterpolationSetUp - Compute spatial indices for point location during interpolation 318 319 Collective 320 321 Input Parameters: 322 + ctx - the context 323 . dm - the `DM` for the function space used for interpolation 324 . redundantPoints - If `PETSC_TRUE`, all processes are passing in the same array of points. Otherwise, points need to be communicated among processes. 325 - ignoreOutsideDomain - If `PETSC_TRUE`, ignore points outside the domain, otherwise return an error 326 327 Level: intermediate 328 329 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 330 @*/ 331 PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain) 332 { 333 MPI_Comm comm = ctx->comm; 334 PetscScalar *a; 335 PetscInt p, q, i; 336 PetscMPIInt rank, size; 337 Vec pointVec; 338 PetscSF cellSF; 339 PetscLayout layout; 340 PetscReal *globalPoints; 341 PetscScalar *globalPointsScalar; 342 const PetscInt *ranges; 343 PetscMPIInt *counts, *displs; 344 const PetscSFNode *foundCells; 345 const PetscInt *foundPoints; 346 PetscMPIInt *foundProcs, *globalProcs; 347 PetscInt n, N, numFound; 348 349 PetscFunctionBegin; 350 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 351 PetscCallMPI(MPI_Comm_size(comm, &size)); 352 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 353 PetscCheck(ctx->dim >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set"); 354 /* Locate points */ 355 n = ctx->nInput; 356 if (!redundantPoints) { 357 PetscCall(PetscLayoutCreate(comm, &layout)); 358 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 359 PetscCall(PetscLayoutSetLocalSize(layout, n)); 360 PetscCall(PetscLayoutSetUp(layout)); 361 PetscCall(PetscLayoutGetSize(layout, &N)); 362 /* Communicate all points to all processes */ 363 PetscCall(PetscMalloc3(N * ctx->dim, &globalPoints, size, &counts, size, &displs)); 364 PetscCall(PetscLayoutGetRanges(layout, &ranges)); 365 for (p = 0; p < size; ++p) { 366 counts[p] = (ranges[p + 1] - ranges[p]) * ctx->dim; 367 displs[p] = ranges[p] * ctx->dim; 368 } 369 PetscCallMPI(MPI_Allgatherv(ctx->points, n * ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm)); 370 } else { 371 N = n; 372 globalPoints = ctx->points; 373 counts = displs = NULL; 374 layout = NULL; 375 } 376 #if 0 377 PetscCall(PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs)); 378 /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */ 379 #else 380 #if defined(PETSC_USE_COMPLEX) 381 PetscCall(PetscMalloc1(N * ctx->dim, &globalPointsScalar)); 382 for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i]; 383 #else 384 globalPointsScalar = globalPoints; 385 #endif 386 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec)); 387 PetscCall(PetscMalloc2(N, &foundProcs, N, &globalProcs)); 388 for (p = 0; p < N; ++p) foundProcs[p] = size; 389 cellSF = NULL; 390 PetscCall(DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF)); 391 PetscCall(PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells)); 392 #endif 393 for (p = 0; p < numFound; ++p) { 394 if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank; 395 } 396 /* Let the lowest rank process own each point */ 397 PetscCall(MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm)); 398 ctx->n = 0; 399 for (p = 0; p < N; ++p) { 400 if (globalProcs[p] == size) { 401 PetscCheck(ignoreOutsideDomain, comm, PETSC_ERR_PLIB, "Point %" PetscInt_FMT ": %g %g %g not located in mesh", p, (double)globalPoints[p * ctx->dim + 0], (double)(ctx->dim > 1 ? globalPoints[p * ctx->dim + 1] : 0.0), 402 (double)(ctx->dim > 2 ? globalPoints[p * ctx->dim + 2] : 0.0)); 403 if (rank == 0) ++ctx->n; 404 } else if (globalProcs[p] == rank) ++ctx->n; 405 } 406 /* Create coordinates vector and array of owned cells */ 407 PetscCall(PetscMalloc1(ctx->n, &ctx->cells)); 408 PetscCall(VecCreate(comm, &ctx->coords)); 409 PetscCall(VecSetSizes(ctx->coords, ctx->n * ctx->dim, PETSC_DECIDE)); 410 PetscCall(VecSetBlockSize(ctx->coords, ctx->dim)); 411 PetscCall(VecSetType(ctx->coords, VECSTANDARD)); 412 PetscCall(VecGetArray(ctx->coords, &a)); 413 for (p = 0, q = 0, i = 0; p < N; ++p) { 414 if (globalProcs[p] == rank) { 415 PetscInt d; 416 417 for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p * ctx->dim + d]; 418 ctx->cells[q] = foundCells[q].index; 419 ++q; 420 } 421 if (globalProcs[p] == size && rank == 0) { 422 PetscInt d; 423 424 for (d = 0; d < ctx->dim; ++d, ++i) a[i] = 0.; 425 ctx->cells[q] = -1; 426 ++q; 427 } 428 } 429 PetscCall(VecRestoreArray(ctx->coords, &a)); 430 #if 0 431 PetscCall(PetscFree3(foundCells,foundProcs,globalProcs)); 432 #else 433 PetscCall(PetscFree2(foundProcs, globalProcs)); 434 PetscCall(PetscSFDestroy(&cellSF)); 435 PetscCall(VecDestroy(&pointVec)); 436 #endif 437 if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar)); 438 if (!redundantPoints) PetscCall(PetscFree3(globalPoints, counts, displs)); 439 PetscCall(PetscLayoutDestroy(&layout)); 440 PetscFunctionReturn(PETSC_SUCCESS); 441 } 442 443 /*@C 444 DMInterpolationGetCoordinates - Gets a `Vec` with the coordinates of each interpolation point 445 446 Collective 447 448 Input Parameter: 449 . ctx - the context 450 451 Output Parameter: 452 . coordinates - the coordinates of interpolation points 453 454 Level: intermediate 455 456 Note: 457 The local vector entries correspond to interpolation points lying on this process, according to the associated `DM`. 458 This is a borrowed vector that the user should not destroy. 459 460 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 461 @*/ 462 PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates) 463 { 464 PetscFunctionBegin; 465 PetscAssertPointer(coordinates, 2); 466 PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup."); 467 *coordinates = ctx->coords; 468 PetscFunctionReturn(PETSC_SUCCESS); 469 } 470 471 /*@C 472 DMInterpolationGetVector - Gets a `Vec` which can hold all the interpolated field values 473 474 Collective 475 476 Input Parameter: 477 . ctx - the context 478 479 Output Parameter: 480 . v - a vector capable of holding the interpolated field values 481 482 Level: intermediate 483 484 Note: 485 This vector should be returned using `DMInterpolationRestoreVector()`. 486 487 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationRestoreVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 488 @*/ 489 PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v) 490 { 491 PetscFunctionBegin; 492 PetscAssertPointer(v, 2); 493 PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup."); 494 PetscCall(VecCreate(ctx->comm, v)); 495 PetscCall(VecSetSizes(*v, ctx->n * ctx->dof, PETSC_DECIDE)); 496 PetscCall(VecSetBlockSize(*v, ctx->dof)); 497 PetscCall(VecSetType(*v, VECSTANDARD)); 498 PetscFunctionReturn(PETSC_SUCCESS); 499 } 500 501 /*@C 502 DMInterpolationRestoreVector - Returns a `Vec` which can hold all the interpolated field values 503 504 Collective 505 506 Input Parameters: 507 + ctx - the context 508 - v - a vector capable of holding the interpolated field values 509 510 Level: intermediate 511 512 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 513 @*/ 514 PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v) 515 { 516 PetscFunctionBegin; 517 PetscAssertPointer(v, 2); 518 PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup."); 519 PetscCall(VecDestroy(v)); 520 PetscFunctionReturn(PETSC_SUCCESS); 521 } 522 523 static inline PetscErrorCode DMInterpolate_Segment_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 524 { 525 const PetscInt c = ctx->cells[p]; 526 const PetscInt dof = ctx->dof; 527 PetscScalar *x = NULL; 528 const PetscScalar *coords; 529 PetscScalar *a; 530 PetscReal v0, J, invJ, detJ, xir[1]; 531 PetscInt xSize, comp; 532 533 PetscFunctionBegin; 534 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 535 PetscCall(VecGetArray(v, &a)); 536 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ)); 537 PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c); 538 xir[0] = invJ * PetscRealPart(coords[p] - v0); 539 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x)); 540 if (2 * dof == xSize) { 541 for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) + x[1 * dof + comp] * xir[0]; 542 } else if (dof == xSize) { 543 for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp]; 544 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Input closure size %" PetscInt_FMT " must be either %" PetscInt_FMT " or %" PetscInt_FMT, xSize, 2 * dof, dof); 545 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x)); 546 PetscCall(VecRestoreArray(v, &a)); 547 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 548 PetscFunctionReturn(PETSC_SUCCESS); 549 } 550 551 static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 552 { 553 const PetscInt c = ctx->cells[p]; 554 PetscScalar *x = NULL; 555 const PetscScalar *coords; 556 PetscScalar *a; 557 PetscReal *v0, *J, *invJ, detJ; 558 PetscReal xi[4]; 559 560 PetscFunctionBegin; 561 PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ)); 562 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 563 PetscCall(VecGetArray(v, &a)); 564 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 565 PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c); 566 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x)); 567 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp]; 568 for (PetscInt d = 0; d < ctx->dim; ++d) { 569 xi[d] = 0.0; 570 for (PetscInt f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]); 571 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] += PetscRealPart(x[(d + 1) * ctx->dof + comp] - x[0 * ctx->dof + comp]) * xi[d]; 572 } 573 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x)); 574 PetscCall(VecRestoreArray(v, &a)); 575 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 576 PetscCall(PetscFree3(v0, J, invJ)); 577 PetscFunctionReturn(PETSC_SUCCESS); 578 } 579 580 static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 581 { 582 const PetscInt c = ctx->cells[p]; 583 const PetscInt order[3] = {2, 1, 3}; 584 PetscScalar *x = NULL; 585 PetscReal *v0, *J, *invJ, detJ; 586 const PetscScalar *coords; 587 PetscScalar *a; 588 PetscReal xi[4]; 589 590 PetscFunctionBegin; 591 PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ)); 592 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 593 PetscCall(VecGetArray(v, &a)); 594 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 595 PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c); 596 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x)); 597 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp]; 598 for (PetscInt d = 0; d < ctx->dim; ++d) { 599 xi[d] = 0.0; 600 for (PetscInt f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]); 601 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] += PetscRealPart(x[order[d] * ctx->dof + comp] - x[0 * ctx->dof + comp]) * xi[d]; 602 } 603 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x)); 604 PetscCall(VecRestoreArray(v, &a)); 605 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 606 PetscCall(PetscFree3(v0, J, invJ)); 607 PetscFunctionReturn(PETSC_SUCCESS); 608 } 609 610 static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx) 611 { 612 const PetscScalar *vertices = (const PetscScalar *)ctx; 613 const PetscScalar x0 = vertices[0]; 614 const PetscScalar y0 = vertices[1]; 615 const PetscScalar x1 = vertices[2]; 616 const PetscScalar y1 = vertices[3]; 617 const PetscScalar x2 = vertices[4]; 618 const PetscScalar y2 = vertices[5]; 619 const PetscScalar x3 = vertices[6]; 620 const PetscScalar y3 = vertices[7]; 621 const PetscScalar f_1 = x1 - x0; 622 const PetscScalar g_1 = y1 - y0; 623 const PetscScalar f_3 = x3 - x0; 624 const PetscScalar g_3 = y3 - y0; 625 const PetscScalar f_01 = x2 - x1 - x3 + x0; 626 const PetscScalar g_01 = y2 - y1 - y3 + y0; 627 const PetscScalar *ref; 628 PetscScalar *real; 629 630 PetscFunctionBegin; 631 PetscCall(VecGetArrayRead(Xref, &ref)); 632 PetscCall(VecGetArray(Xreal, &real)); 633 { 634 const PetscScalar p0 = ref[0]; 635 const PetscScalar p1 = ref[1]; 636 637 real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1; 638 real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1; 639 } 640 PetscCall(PetscLogFlops(28)); 641 PetscCall(VecRestoreArrayRead(Xref, &ref)); 642 PetscCall(VecRestoreArray(Xreal, &real)); 643 PetscFunctionReturn(PETSC_SUCCESS); 644 } 645 646 #include <petsc/private/dmimpl.h> 647 static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx) 648 { 649 const PetscScalar *vertices = (const PetscScalar *)ctx; 650 const PetscScalar x0 = vertices[0]; 651 const PetscScalar y0 = vertices[1]; 652 const PetscScalar x1 = vertices[2]; 653 const PetscScalar y1 = vertices[3]; 654 const PetscScalar x2 = vertices[4]; 655 const PetscScalar y2 = vertices[5]; 656 const PetscScalar x3 = vertices[6]; 657 const PetscScalar y3 = vertices[7]; 658 const PetscScalar f_01 = x2 - x1 - x3 + x0; 659 const PetscScalar g_01 = y2 - y1 - y3 + y0; 660 const PetscScalar *ref; 661 662 PetscFunctionBegin; 663 PetscCall(VecGetArrayRead(Xref, &ref)); 664 { 665 const PetscScalar x = ref[0]; 666 const PetscScalar y = ref[1]; 667 const PetscInt rows[2] = {0, 1}; 668 PetscScalar values[4]; 669 670 values[0] = (x1 - x0 + f_01 * y) * 0.5; 671 values[1] = (x3 - x0 + f_01 * x) * 0.5; 672 values[2] = (y1 - y0 + g_01 * y) * 0.5; 673 values[3] = (y3 - y0 + g_01 * x) * 0.5; 674 PetscCall(MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES)); 675 } 676 PetscCall(PetscLogFlops(30)); 677 PetscCall(VecRestoreArrayRead(Xref, &ref)); 678 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 679 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 680 PetscFunctionReturn(PETSC_SUCCESS); 681 } 682 683 static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 684 { 685 const PetscInt c = ctx->cells[p]; 686 PetscFE fem = NULL; 687 DM dmCoord; 688 SNES snes; 689 KSP ksp; 690 PC pc; 691 Vec coordsLocal, r, ref, real; 692 Mat J; 693 PetscTabulation T = NULL; 694 const PetscScalar *coords; 695 PetscScalar *a; 696 PetscReal xir[2] = {0., 0.}; 697 PetscInt Nf; 698 const PetscInt dof = ctx->dof; 699 PetscScalar *x = NULL, *vertices = NULL; 700 PetscScalar *xi; 701 PetscInt coordSize, xSize; 702 703 PetscFunctionBegin; 704 PetscCall(DMGetNumFields(dm, &Nf)); 705 if (Nf) { 706 PetscObject obj; 707 PetscClassId id; 708 709 PetscCall(DMGetField(dm, 0, NULL, &obj)); 710 PetscCall(PetscObjectGetClassId(obj, &id)); 711 if (id == PETSCFE_CLASSID) { 712 fem = (PetscFE)obj; 713 PetscCall(PetscFECreateTabulation(fem, 1, 1, xir, 0, &T)); 714 } 715 } 716 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 717 PetscCall(DMGetCoordinateDM(dm, &dmCoord)); 718 PetscCall(SNESCreate(PETSC_COMM_SELF, &snes)); 719 PetscCall(SNESSetOptionsPrefix(snes, "quad_interp_")); 720 PetscCall(VecCreate(PETSC_COMM_SELF, &r)); 721 PetscCall(VecSetSizes(r, 2, 2)); 722 PetscCall(VecSetType(r, dm->vectype)); 723 PetscCall(VecDuplicate(r, &ref)); 724 PetscCall(VecDuplicate(r, &real)); 725 PetscCall(MatCreate(PETSC_COMM_SELF, &J)); 726 PetscCall(MatSetSizes(J, 2, 2, 2, 2)); 727 PetscCall(MatSetType(J, MATSEQDENSE)); 728 PetscCall(MatSetUp(J)); 729 PetscCall(SNESSetFunction(snes, r, QuadMap_Private, NULL)); 730 PetscCall(SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL)); 731 PetscCall(SNESGetKSP(snes, &ksp)); 732 PetscCall(KSPGetPC(ksp, &pc)); 733 PetscCall(PCSetType(pc, PCLU)); 734 PetscCall(SNESSetFromOptions(snes)); 735 736 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 737 PetscCall(VecGetArray(v, &a)); 738 PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 739 PetscCheck(4 * 2 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %" PetscInt_FMT " should be %d", coordSize, 4 * 2); 740 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x)); 741 PetscCall(SNESSetFunction(snes, NULL, NULL, vertices)); 742 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices)); 743 PetscCall(VecGetArray(real, &xi)); 744 xi[0] = coords[p * ctx->dim + 0]; 745 xi[1] = coords[p * ctx->dim + 1]; 746 PetscCall(VecRestoreArray(real, &xi)); 747 PetscCall(SNESSolve(snes, real, ref)); 748 PetscCall(VecGetArray(ref, &xi)); 749 xir[0] = PetscRealPart(xi[0]); 750 xir[1] = PetscRealPart(xi[1]); 751 if (4 * dof == xSize) { 752 for (PetscInt comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) * (1 - xir[1]) + x[1 * dof + comp] * xir[0] * (1 - xir[1]) + x[2 * dof + comp] * xir[0] * xir[1] + x[3 * dof + comp] * (1 - xir[0]) * xir[1]; 753 } else if (dof == xSize) { 754 for (PetscInt comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp]; 755 } else { 756 PetscCheck(fem, ctx->comm, PETSC_ERR_ARG_WRONG, "Cannot have a higher order interpolant if the discretization is not PetscFE"); 757 xir[0] = 2.0 * xir[0] - 1.0; 758 xir[1] = 2.0 * xir[1] - 1.0; 759 PetscCall(PetscFEComputeTabulation(fem, 1, xir, 0, T)); 760 for (PetscInt comp = 0; comp < dof; ++comp) { 761 a[p * dof + comp] = 0.0; 762 for (PetscInt d = 0; d < xSize / dof; ++d) a[p * dof + comp] += x[d * dof + comp] * T->T[0][d * dof + comp]; 763 } 764 } 765 PetscCall(VecRestoreArray(ref, &xi)); 766 PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 767 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x)); 768 PetscCall(PetscTabulationDestroy(&T)); 769 PetscCall(VecRestoreArray(v, &a)); 770 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 771 772 PetscCall(SNESDestroy(&snes)); 773 PetscCall(VecDestroy(&r)); 774 PetscCall(VecDestroy(&ref)); 775 PetscCall(VecDestroy(&real)); 776 PetscCall(MatDestroy(&J)); 777 PetscFunctionReturn(PETSC_SUCCESS); 778 } 779 780 static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx) 781 { 782 const PetscScalar *vertices = (const PetscScalar *)ctx; 783 const PetscScalar x0 = vertices[0]; 784 const PetscScalar y0 = vertices[1]; 785 const PetscScalar z0 = vertices[2]; 786 const PetscScalar x1 = vertices[9]; 787 const PetscScalar y1 = vertices[10]; 788 const PetscScalar z1 = vertices[11]; 789 const PetscScalar x2 = vertices[6]; 790 const PetscScalar y2 = vertices[7]; 791 const PetscScalar z2 = vertices[8]; 792 const PetscScalar x3 = vertices[3]; 793 const PetscScalar y3 = vertices[4]; 794 const PetscScalar z3 = vertices[5]; 795 const PetscScalar x4 = vertices[12]; 796 const PetscScalar y4 = vertices[13]; 797 const PetscScalar z4 = vertices[14]; 798 const PetscScalar x5 = vertices[15]; 799 const PetscScalar y5 = vertices[16]; 800 const PetscScalar z5 = vertices[17]; 801 const PetscScalar x6 = vertices[18]; 802 const PetscScalar y6 = vertices[19]; 803 const PetscScalar z6 = vertices[20]; 804 const PetscScalar x7 = vertices[21]; 805 const PetscScalar y7 = vertices[22]; 806 const PetscScalar z7 = vertices[23]; 807 const PetscScalar f_1 = x1 - x0; 808 const PetscScalar g_1 = y1 - y0; 809 const PetscScalar h_1 = z1 - z0; 810 const PetscScalar f_3 = x3 - x0; 811 const PetscScalar g_3 = y3 - y0; 812 const PetscScalar h_3 = z3 - z0; 813 const PetscScalar f_4 = x4 - x0; 814 const PetscScalar g_4 = y4 - y0; 815 const PetscScalar h_4 = z4 - z0; 816 const PetscScalar f_01 = x2 - x1 - x3 + x0; 817 const PetscScalar g_01 = y2 - y1 - y3 + y0; 818 const PetscScalar h_01 = z2 - z1 - z3 + z0; 819 const PetscScalar f_12 = x7 - x3 - x4 + x0; 820 const PetscScalar g_12 = y7 - y3 - y4 + y0; 821 const PetscScalar h_12 = z7 - z3 - z4 + z0; 822 const PetscScalar f_02 = x5 - x1 - x4 + x0; 823 const PetscScalar g_02 = y5 - y1 - y4 + y0; 824 const PetscScalar h_02 = z5 - z1 - z4 + z0; 825 const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7; 826 const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7; 827 const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7; 828 const PetscScalar *ref; 829 PetscScalar *real; 830 831 PetscFunctionBegin; 832 PetscCall(VecGetArrayRead(Xref, &ref)); 833 PetscCall(VecGetArray(Xreal, &real)); 834 { 835 const PetscScalar p0 = ref[0]; 836 const PetscScalar p1 = ref[1]; 837 const PetscScalar p2 = ref[2]; 838 839 real[0] = x0 + f_1 * p0 + f_3 * p1 + f_4 * p2 + f_01 * p0 * p1 + f_12 * p1 * p2 + f_02 * p0 * p2 + f_012 * p0 * p1 * p2; 840 real[1] = y0 + g_1 * p0 + g_3 * p1 + g_4 * p2 + g_01 * p0 * p1 + g_01 * p0 * p1 + g_12 * p1 * p2 + g_02 * p0 * p2 + g_012 * p0 * p1 * p2; 841 real[2] = z0 + h_1 * p0 + h_3 * p1 + h_4 * p2 + h_01 * p0 * p1 + h_01 * p0 * p1 + h_12 * p1 * p2 + h_02 * p0 * p2 + h_012 * p0 * p1 * p2; 842 } 843 PetscCall(PetscLogFlops(114)); 844 PetscCall(VecRestoreArrayRead(Xref, &ref)); 845 PetscCall(VecRestoreArray(Xreal, &real)); 846 PetscFunctionReturn(PETSC_SUCCESS); 847 } 848 849 static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx) 850 { 851 const PetscScalar *vertices = (const PetscScalar *)ctx; 852 const PetscScalar x0 = vertices[0]; 853 const PetscScalar y0 = vertices[1]; 854 const PetscScalar z0 = vertices[2]; 855 const PetscScalar x1 = vertices[9]; 856 const PetscScalar y1 = vertices[10]; 857 const PetscScalar z1 = vertices[11]; 858 const PetscScalar x2 = vertices[6]; 859 const PetscScalar y2 = vertices[7]; 860 const PetscScalar z2 = vertices[8]; 861 const PetscScalar x3 = vertices[3]; 862 const PetscScalar y3 = vertices[4]; 863 const PetscScalar z3 = vertices[5]; 864 const PetscScalar x4 = vertices[12]; 865 const PetscScalar y4 = vertices[13]; 866 const PetscScalar z4 = vertices[14]; 867 const PetscScalar x5 = vertices[15]; 868 const PetscScalar y5 = vertices[16]; 869 const PetscScalar z5 = vertices[17]; 870 const PetscScalar x6 = vertices[18]; 871 const PetscScalar y6 = vertices[19]; 872 const PetscScalar z6 = vertices[20]; 873 const PetscScalar x7 = vertices[21]; 874 const PetscScalar y7 = vertices[22]; 875 const PetscScalar z7 = vertices[23]; 876 const PetscScalar f_xy = x2 - x1 - x3 + x0; 877 const PetscScalar g_xy = y2 - y1 - y3 + y0; 878 const PetscScalar h_xy = z2 - z1 - z3 + z0; 879 const PetscScalar f_yz = x7 - x3 - x4 + x0; 880 const PetscScalar g_yz = y7 - y3 - y4 + y0; 881 const PetscScalar h_yz = z7 - z3 - z4 + z0; 882 const PetscScalar f_xz = x5 - x1 - x4 + x0; 883 const PetscScalar g_xz = y5 - y1 - y4 + y0; 884 const PetscScalar h_xz = z5 - z1 - z4 + z0; 885 const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7; 886 const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7; 887 const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7; 888 const PetscScalar *ref; 889 890 PetscFunctionBegin; 891 PetscCall(VecGetArrayRead(Xref, &ref)); 892 { 893 const PetscScalar x = ref[0]; 894 const PetscScalar y = ref[1]; 895 const PetscScalar z = ref[2]; 896 const PetscInt rows[3] = {0, 1, 2}; 897 PetscScalar values[9]; 898 899 values[0] = (x1 - x0 + f_xy * y + f_xz * z + f_xyz * y * z) / 2.0; 900 values[1] = (x3 - x0 + f_xy * x + f_yz * z + f_xyz * x * z) / 2.0; 901 values[2] = (x4 - x0 + f_yz * y + f_xz * x + f_xyz * x * y) / 2.0; 902 values[3] = (y1 - y0 + g_xy * y + g_xz * z + g_xyz * y * z) / 2.0; 903 values[4] = (y3 - y0 + g_xy * x + g_yz * z + g_xyz * x * z) / 2.0; 904 values[5] = (y4 - y0 + g_yz * y + g_xz * x + g_xyz * x * y) / 2.0; 905 values[6] = (z1 - z0 + h_xy * y + h_xz * z + h_xyz * y * z) / 2.0; 906 values[7] = (z3 - z0 + h_xy * x + h_yz * z + h_xyz * x * z) / 2.0; 907 values[8] = (z4 - z0 + h_yz * y + h_xz * x + h_xyz * x * y) / 2.0; 908 909 PetscCall(MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES)); 910 } 911 PetscCall(PetscLogFlops(152)); 912 PetscCall(VecRestoreArrayRead(Xref, &ref)); 913 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 914 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 915 PetscFunctionReturn(PETSC_SUCCESS); 916 } 917 918 static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 919 { 920 const PetscInt c = ctx->cells[p]; 921 DM dmCoord; 922 SNES snes; 923 KSP ksp; 924 PC pc; 925 Vec coordsLocal, r, ref, real; 926 Mat J; 927 const PetscScalar *coords; 928 PetscScalar *a; 929 PetscScalar *x = NULL, *vertices = NULL; 930 PetscScalar *xi; 931 PetscReal xir[3]; 932 PetscInt coordSize, xSize; 933 934 PetscFunctionBegin; 935 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 936 PetscCall(DMGetCoordinateDM(dm, &dmCoord)); 937 PetscCall(SNESCreate(PETSC_COMM_SELF, &snes)); 938 PetscCall(SNESSetOptionsPrefix(snes, "hex_interp_")); 939 PetscCall(VecCreate(PETSC_COMM_SELF, &r)); 940 PetscCall(VecSetSizes(r, 3, 3)); 941 PetscCall(VecSetType(r, dm->vectype)); 942 PetscCall(VecDuplicate(r, &ref)); 943 PetscCall(VecDuplicate(r, &real)); 944 PetscCall(MatCreate(PETSC_COMM_SELF, &J)); 945 PetscCall(MatSetSizes(J, 3, 3, 3, 3)); 946 PetscCall(MatSetType(J, MATSEQDENSE)); 947 PetscCall(MatSetUp(J)); 948 PetscCall(SNESSetFunction(snes, r, HexMap_Private, NULL)); 949 PetscCall(SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL)); 950 PetscCall(SNESGetKSP(snes, &ksp)); 951 PetscCall(KSPGetPC(ksp, &pc)); 952 PetscCall(PCSetType(pc, PCLU)); 953 PetscCall(SNESSetFromOptions(snes)); 954 955 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 956 PetscCall(VecGetArray(v, &a)); 957 PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 958 PetscCheck(8 * 3 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid coordinate closure size %" PetscInt_FMT " should be %d", coordSize, 8 * 3); 959 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x)); 960 PetscCheck((8 * ctx->dof == xSize) || (ctx->dof == xSize), ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input closure size %" PetscInt_FMT " should be %" PetscInt_FMT " or %" PetscInt_FMT, xSize, 8 * ctx->dof, ctx->dof); 961 PetscCall(SNESSetFunction(snes, NULL, NULL, vertices)); 962 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices)); 963 PetscCall(VecGetArray(real, &xi)); 964 xi[0] = coords[p * ctx->dim + 0]; 965 xi[1] = coords[p * ctx->dim + 1]; 966 xi[2] = coords[p * ctx->dim + 2]; 967 PetscCall(VecRestoreArray(real, &xi)); 968 PetscCall(SNESSolve(snes, real, ref)); 969 PetscCall(VecGetArray(ref, &xi)); 970 xir[0] = PetscRealPart(xi[0]); 971 xir[1] = PetscRealPart(xi[1]); 972 xir[2] = PetscRealPart(xi[2]); 973 if (8 * ctx->dof == xSize) { 974 for (PetscInt comp = 0; comp < ctx->dof; ++comp) { 975 a[p * ctx->dof + comp] = x[0 * ctx->dof + comp] * (1 - xir[0]) * (1 - xir[1]) * (1 - xir[2]) + x[3 * ctx->dof + comp] * xir[0] * (1 - xir[1]) * (1 - xir[2]) + x[2 * ctx->dof + comp] * xir[0] * xir[1] * (1 - xir[2]) + x[1 * ctx->dof + comp] * (1 - xir[0]) * xir[1] * (1 - xir[2]) + 976 x[4 * ctx->dof + comp] * (1 - xir[0]) * (1 - xir[1]) * xir[2] + x[5 * ctx->dof + comp] * xir[0] * (1 - xir[1]) * xir[2] + x[6 * ctx->dof + comp] * xir[0] * xir[1] * xir[2] + x[7 * ctx->dof + comp] * (1 - xir[0]) * xir[1] * xir[2]; 977 } 978 } else { 979 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp]; 980 } 981 PetscCall(VecRestoreArray(ref, &xi)); 982 PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 983 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x)); 984 PetscCall(VecRestoreArray(v, &a)); 985 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 986 987 PetscCall(SNESDestroy(&snes)); 988 PetscCall(VecDestroy(&r)); 989 PetscCall(VecDestroy(&ref)); 990 PetscCall(VecDestroy(&real)); 991 PetscCall(MatDestroy(&J)); 992 PetscFunctionReturn(PETSC_SUCCESS); 993 } 994 995 /*@C 996 DMInterpolationEvaluate - Using the input from dm and x, calculates interpolated field values at the interpolation points. 997 998 Input Parameters: 999 + ctx - The `DMInterpolationInfo` context 1000 . dm - The `DM` 1001 - x - The local vector containing the field to be interpolated 1002 1003 Output Parameter: 1004 . v - The vector containing the interpolated values 1005 1006 Level: beginner 1007 1008 Note: 1009 A suitable `v` can be obtained using `DMInterpolationGetVector()`. 1010 1011 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 1012 @*/ 1013 PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v) 1014 { 1015 PetscDS ds; 1016 PetscInt n, p, Nf, field; 1017 PetscBool useDS = PETSC_FALSE; 1018 1019 PetscFunctionBegin; 1020 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1021 PetscValidHeaderSpecific(x, VEC_CLASSID, 3); 1022 PetscValidHeaderSpecific(v, VEC_CLASSID, 4); 1023 PetscCall(VecGetLocalSize(v, &n)); 1024 PetscCheck(n == ctx->n * ctx->dof, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input vector size %" PetscInt_FMT " should be %" PetscInt_FMT, n, ctx->n * ctx->dof); 1025 if (!n) PetscFunctionReturn(PETSC_SUCCESS); 1026 PetscCall(DMGetDS(dm, &ds)); 1027 if (ds) { 1028 useDS = PETSC_TRUE; 1029 PetscCall(PetscDSGetNumFields(ds, &Nf)); 1030 for (field = 0; field < Nf; ++field) { 1031 PetscObject obj; 1032 PetscClassId id; 1033 1034 PetscCall(PetscDSGetDiscretization(ds, field, &obj)); 1035 PetscCall(PetscObjectGetClassId(obj, &id)); 1036 if (id != PETSCFE_CLASSID && id != PETSCFV_CLASSID) { 1037 useDS = PETSC_FALSE; 1038 break; 1039 } 1040 } 1041 } 1042 if (useDS) { 1043 const PetscScalar *coords; 1044 PetscScalar *interpolant; 1045 PetscInt cdim, d; 1046 1047 PetscCall(DMGetCoordinateDim(dm, &cdim)); 1048 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 1049 PetscCall(VecGetArrayWrite(v, &interpolant)); 1050 for (p = 0; p < ctx->n; ++p) { 1051 PetscReal pcoords[3], xi[3]; 1052 PetscScalar *xa = NULL; 1053 PetscInt coff = 0, foff = 0, clSize; 1054 1055 if (ctx->cells[p] < 0) continue; 1056 for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]); 1057 PetscCall(DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi)); 1058 PetscCall(DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa)); 1059 for (field = 0; field < Nf; ++field) { 1060 PetscTabulation T; 1061 PetscObject obj; 1062 PetscClassId id; 1063 1064 PetscCall(PetscDSGetDiscretization(ds, field, &obj)); 1065 PetscCall(PetscObjectGetClassId(obj, &id)); 1066 if (id == PETSCFE_CLASSID) { 1067 PetscFE fe = (PetscFE)obj; 1068 1069 PetscCall(PetscFECreateTabulation(fe, 1, 1, xi, 0, &T)); 1070 { 1071 const PetscReal *basis = T->T[0]; 1072 const PetscInt Nb = T->Nb; 1073 const PetscInt Nc = T->Nc; 1074 1075 for (PetscInt fc = 0; fc < Nc; ++fc) { 1076 interpolant[p * ctx->dof + coff + fc] = 0.0; 1077 for (PetscInt f = 0; f < Nb; ++f) interpolant[p * ctx->dof + coff + fc] += xa[foff + f] * basis[(0 * Nb + f) * Nc + fc]; 1078 } 1079 coff += Nc; 1080 foff += Nb; 1081 } 1082 PetscCall(PetscTabulationDestroy(&T)); 1083 } else if (id == PETSCFV_CLASSID) { 1084 PetscFV fv = (PetscFV)obj; 1085 PetscInt Nc; 1086 1087 // TODO Could use reconstruction if available 1088 PetscCall(PetscFVGetNumComponents(fv, &Nc)); 1089 for (PetscInt fc = 0; fc < Nc; ++fc) interpolant[p * ctx->dof + coff + fc] = xa[foff + fc]; 1090 coff += Nc; 1091 foff += Nc; 1092 } 1093 } 1094 PetscCall(DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa)); 1095 PetscCheck(coff == ctx->dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total components %" PetscInt_FMT " != %" PetscInt_FMT " dof specified for interpolation", coff, ctx->dof); 1096 PetscCheck(foff == clSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total FE/FV space size %" PetscInt_FMT " != %" PetscInt_FMT " closure size", foff, clSize); 1097 } 1098 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 1099 PetscCall(VecRestoreArrayWrite(v, &interpolant)); 1100 } else { 1101 for (PetscInt p = 0; p < ctx->n; ++p) { 1102 const PetscInt cell = ctx->cells[p]; 1103 DMPolytopeType ct; 1104 1105 PetscCall(DMPlexGetCellType(dm, cell, &ct)); 1106 switch (ct) { 1107 case DM_POLYTOPE_SEGMENT: 1108 PetscCall(DMInterpolate_Segment_Private(ctx, dm, p, x, v)); 1109 break; 1110 case DM_POLYTOPE_TRIANGLE: 1111 PetscCall(DMInterpolate_Triangle_Private(ctx, dm, p, x, v)); 1112 break; 1113 case DM_POLYTOPE_QUADRILATERAL: 1114 PetscCall(DMInterpolate_Quad_Private(ctx, dm, p, x, v)); 1115 break; 1116 case DM_POLYTOPE_TETRAHEDRON: 1117 PetscCall(DMInterpolate_Tetrahedron_Private(ctx, dm, p, x, v)); 1118 break; 1119 case DM_POLYTOPE_HEXAHEDRON: 1120 PetscCall(DMInterpolate_Hex_Private(ctx, dm, cell, x, v)); 1121 break; 1122 default: 1123 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for cell type %s", DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]); 1124 } 1125 } 1126 } 1127 PetscFunctionReturn(PETSC_SUCCESS); 1128 } 1129 1130 /*@C 1131 DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context 1132 1133 Collective 1134 1135 Input Parameter: 1136 . ctx - the context 1137 1138 Level: beginner 1139 1140 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 1141 @*/ 1142 PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx) 1143 { 1144 PetscFunctionBegin; 1145 PetscAssertPointer(ctx, 1); 1146 PetscCall(VecDestroy(&(*ctx)->coords)); 1147 PetscCall(PetscFree((*ctx)->points)); 1148 PetscCall(PetscFree((*ctx)->cells)); 1149 PetscCall(PetscFree(*ctx)); 1150 *ctx = NULL; 1151 PetscFunctionReturn(PETSC_SUCCESS); 1152 } 1153 1154 /*@C 1155 SNESMonitorFields - Monitors the residual for each field separately 1156 1157 Collective 1158 1159 Input Parameters: 1160 + snes - the `SNES` context, must have an attached `DM` 1161 . its - iteration number 1162 . fgnorm - 2-norm of residual 1163 - vf - `PetscViewerAndFormat` of `PetscViewerType` `PETSCVIEWERASCII` 1164 1165 Level: intermediate 1166 1167 Note: 1168 This routine prints the residual norm at each iteration. 1169 1170 .seealso: [](ch_snes), `SNES`, `SNESMonitorSet()`, `SNESMonitorDefault()` 1171 @*/ 1172 PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf) 1173 { 1174 PetscViewer viewer = vf->viewer; 1175 Vec res; 1176 DM dm; 1177 PetscSection s; 1178 const PetscScalar *r; 1179 PetscReal *lnorms, *norms; 1180 PetscInt numFields, f, pStart, pEnd, p; 1181 1182 PetscFunctionBegin; 1183 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 4); 1184 PetscCall(SNESGetFunction(snes, &res, NULL, NULL)); 1185 PetscCall(SNESGetDM(snes, &dm)); 1186 PetscCall(DMGetLocalSection(dm, &s)); 1187 PetscCall(PetscSectionGetNumFields(s, &numFields)); 1188 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 1189 PetscCall(PetscCalloc2(numFields, &lnorms, numFields, &norms)); 1190 PetscCall(VecGetArrayRead(res, &r)); 1191 for (p = pStart; p < pEnd; ++p) { 1192 for (f = 0; f < numFields; ++f) { 1193 PetscInt fdof, foff, d; 1194 1195 PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof)); 1196 PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff)); 1197 for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff + d])); 1198 } 1199 } 1200 PetscCall(VecRestoreArrayRead(res, &r)); 1201 PetscCall(MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm))); 1202 PetscCall(PetscViewerPushFormat(viewer, vf->format)); 1203 PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel)); 1204 PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " SNES Function norm %14.12e [", its, (double)fgnorm)); 1205 for (f = 0; f < numFields; ++f) { 1206 if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", ")); 1207 PetscCall(PetscViewerASCIIPrintf(viewer, "%14.12e", (double)PetscSqrtReal(norms[f]))); 1208 } 1209 PetscCall(PetscViewerASCIIPrintf(viewer, "]\n")); 1210 PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel)); 1211 PetscCall(PetscViewerPopFormat(viewer)); 1212 PetscCall(PetscFree2(lnorms, norms)); 1213 PetscFunctionReturn(PETSC_SUCCESS); 1214 } 1215 1216 /********************* SNES callbacks **************************/ 1217 1218 /*@ 1219 DMPlexSNESComputeObjectiveFEM - Sums the local objectives from the local input X using pointwise functions specified by the user 1220 1221 Input Parameters: 1222 + dm - The mesh 1223 . X - Local solution 1224 - user - The user context 1225 1226 Output Parameter: 1227 . obj - Local objective value 1228 1229 Level: developer 1230 1231 .seealso: `DM`, `DMPlexSNESComputeResidualFEM()` 1232 @*/ 1233 PetscErrorCode DMPlexSNESComputeObjectiveFEM(DM dm, Vec X, PetscReal *obj, void *user) 1234 { 1235 PetscInt Nf, cellHeight, cStart, cEnd; 1236 PetscScalar *cintegral; 1237 1238 PetscFunctionBegin; 1239 PetscCall(DMGetNumFields(dm, &Nf)); 1240 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 1241 PetscCall(DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd)); 1242 PetscCall(PetscCalloc1((cEnd - cStart) * Nf, &cintegral)); 1243 PetscCall(PetscLogEventBegin(DMPLEX_IntegralFEM, dm, 0, 0, 0)); 1244 PetscCall(DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user)); 1245 /* Sum up values */ 1246 *obj = 0; 1247 for (PetscInt c = cStart; c < cEnd; ++c) 1248 for (PetscInt f = 0; f < Nf; ++f) *obj += PetscRealPart(cintegral[(c - cStart) * Nf + f]); 1249 PetscCall(PetscLogEventBegin(DMPLEX_IntegralFEM, dm, 0, 0, 0)); 1250 PetscCall(PetscFree(cintegral)); 1251 PetscFunctionReturn(PETSC_SUCCESS); 1252 } 1253 1254 /*@ 1255 DMPlexSNESComputeResidualFEM - Sums the local residual into vector `F` from the local input `X` using pointwise functions specified by the user 1256 1257 Input Parameters: 1258 + dm - The mesh 1259 . X - Local solution 1260 - user - The user context 1261 1262 Output Parameter: 1263 . F - Local output vector 1264 1265 Level: developer 1266 1267 Note: 1268 The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional. 1269 1270 .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMSNESComputeJacobianAction()` 1271 @*/ 1272 PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user) 1273 { 1274 DM plex; 1275 IS allcellIS; 1276 PetscInt Nds, s; 1277 1278 PetscFunctionBegin; 1279 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1280 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS)); 1281 PetscCall(DMGetNumDS(dm, &Nds)); 1282 for (s = 0; s < Nds; ++s) { 1283 PetscDS ds; 1284 IS cellIS; 1285 PetscFormKey key; 1286 1287 PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL)); 1288 key.value = 0; 1289 key.field = 0; 1290 key.part = 0; 1291 if (!key.label) { 1292 PetscCall(PetscObjectReference((PetscObject)allcellIS)); 1293 cellIS = allcellIS; 1294 } else { 1295 IS pointIS; 1296 1297 key.value = 1; 1298 PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS)); 1299 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS)); 1300 PetscCall(ISDestroy(&pointIS)); 1301 } 1302 PetscCall(DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user)); 1303 PetscCall(ISDestroy(&cellIS)); 1304 } 1305 PetscCall(ISDestroy(&allcellIS)); 1306 PetscCall(DMDestroy(&plex)); 1307 PetscFunctionReturn(PETSC_SUCCESS); 1308 } 1309 1310 /*@ 1311 DMPlexSNESComputeResidualDS - Sums the local residual into vector `F` from the local input `X` using all pointwise functions with unique keys in the `PetscDS` 1312 1313 Input Parameters: 1314 + dm - The mesh 1315 . X - Local solution 1316 - user - The user context 1317 1318 Output Parameter: 1319 . F - Local output vector 1320 1321 Level: developer 1322 1323 Note: 1324 The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional. 1325 1326 .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()` 1327 @*/ 1328 PetscErrorCode DMPlexSNESComputeResidualDS(DM dm, Vec X, Vec F, void *user) 1329 { 1330 DM plex; 1331 IS allcellIS; 1332 PetscInt Nds, s; 1333 1334 PetscFunctionBegin; 1335 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1336 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS)); 1337 PetscCall(DMGetNumDS(dm, &Nds)); 1338 for (s = 0; s < Nds; ++s) { 1339 PetscDS ds; 1340 DMLabel label; 1341 IS cellIS; 1342 1343 PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL)); 1344 { 1345 PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1}; 1346 PetscWeakForm wf; 1347 PetscInt Nm = 2, m, Nk = 0, k, kp, off = 0; 1348 PetscFormKey *reskeys; 1349 1350 /* Get unique residual keys */ 1351 for (m = 0; m < Nm; ++m) { 1352 PetscInt Nkm; 1353 PetscCall(PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm)); 1354 Nk += Nkm; 1355 } 1356 PetscCall(PetscMalloc1(Nk, &reskeys)); 1357 for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys)); 1358 PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk); 1359 PetscCall(PetscFormKeySort(Nk, reskeys)); 1360 for (k = 0, kp = 1; kp < Nk; ++kp) { 1361 if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) { 1362 ++k; 1363 if (kp != k) reskeys[k] = reskeys[kp]; 1364 } 1365 } 1366 Nk = k; 1367 1368 PetscCall(PetscDSGetWeakForm(ds, &wf)); 1369 for (k = 0; k < Nk; ++k) { 1370 DMLabel label = reskeys[k].label; 1371 PetscInt val = reskeys[k].value; 1372 1373 if (!label) { 1374 PetscCall(PetscObjectReference((PetscObject)allcellIS)); 1375 cellIS = allcellIS; 1376 } else { 1377 IS pointIS; 1378 1379 PetscCall(DMLabelGetStratumIS(label, val, &pointIS)); 1380 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS)); 1381 PetscCall(ISDestroy(&pointIS)); 1382 } 1383 PetscCall(DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user)); 1384 PetscCall(ISDestroy(&cellIS)); 1385 } 1386 PetscCall(PetscFree(reskeys)); 1387 } 1388 } 1389 PetscCall(ISDestroy(&allcellIS)); 1390 PetscCall(DMDestroy(&plex)); 1391 PetscFunctionReturn(PETSC_SUCCESS); 1392 } 1393 1394 #ifdef PETSC_HAVE_LIBCEED 1395 PetscErrorCode DMPlexSNESComputeResidualCEED(DM dm, Vec locX, Vec locF, void *user) 1396 { 1397 Ceed ceed; 1398 DMCeed sd = dm->dmceed; 1399 CeedVector clocX, clocF; 1400 1401 PetscFunctionBegin; 1402 PetscCall(DMGetCeed(dm, &ceed)); 1403 PetscCheck(sd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This DM has no CEED data. Call DMCeedCreate() before computing the residual."); 1404 PetscCall(DMCeedComputeGeometry(dm, sd)); 1405 1406 PetscCall(VecGetCeedVectorRead(locX, ceed, &clocX)); 1407 PetscCall(VecGetCeedVector(locF, ceed, &clocF)); 1408 PetscCallCEED(CeedOperatorApplyAdd(sd->op, clocX, clocF, CEED_REQUEST_IMMEDIATE)); 1409 PetscCall(VecRestoreCeedVectorRead(locX, &clocX)); 1410 PetscCall(VecRestoreCeedVector(locF, &clocF)); 1411 1412 { 1413 DM_Plex *mesh = (DM_Plex *)dm->data; 1414 1415 if (mesh->printFEM) { 1416 PetscSection section; 1417 Vec locFbc; 1418 PetscInt pStart, pEnd, p, maxDof; 1419 PetscScalar *zeroes; 1420 1421 PetscCall(DMGetLocalSection(dm, §ion)); 1422 PetscCall(VecDuplicate(locF, &locFbc)); 1423 PetscCall(VecCopy(locF, locFbc)); 1424 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 1425 PetscCall(PetscSectionGetMaxDof(section, &maxDof)); 1426 PetscCall(PetscCalloc1(maxDof, &zeroes)); 1427 for (p = pStart; p < pEnd; ++p) PetscCall(VecSetValuesSection(locFbc, section, p, zeroes, INSERT_BC_VALUES)); 1428 PetscCall(PetscFree(zeroes)); 1429 PetscCall(DMPrintLocalVec(dm, "Residual", mesh->printTol, locFbc)); 1430 PetscCall(VecDestroy(&locFbc)); 1431 } 1432 } 1433 PetscFunctionReturn(PETSC_SUCCESS); 1434 } 1435 #endif 1436 1437 /*@ 1438 DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input `X` 1439 1440 Input Parameters: 1441 + dm - The mesh 1442 - user - The user context 1443 1444 Output Parameter: 1445 . X - Local solution 1446 1447 Level: developer 1448 1449 .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()` 1450 @*/ 1451 PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user) 1452 { 1453 DM plex; 1454 1455 PetscFunctionBegin; 1456 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1457 PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL)); 1458 PetscCall(DMDestroy(&plex)); 1459 PetscFunctionReturn(PETSC_SUCCESS); 1460 } 1461 1462 /*@ 1463 DMSNESComputeJacobianAction - Compute the action of the Jacobian J(X) on Y 1464 1465 Input Parameters: 1466 + dm - The `DM` 1467 . X - Local solution vector 1468 . Y - Local input vector 1469 - user - The user context 1470 1471 Output Parameter: 1472 . F - local output vector 1473 1474 Level: developer 1475 1476 Notes: 1477 Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly. 1478 1479 .seealso: [](ch_snes), `DM`, ``DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()` 1480 @*/ 1481 PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, void *user) 1482 { 1483 DM plex; 1484 IS allcellIS; 1485 PetscInt Nds, s; 1486 1487 PetscFunctionBegin; 1488 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1489 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS)); 1490 PetscCall(DMGetNumDS(dm, &Nds)); 1491 for (s = 0; s < Nds; ++s) { 1492 PetscDS ds; 1493 DMLabel label; 1494 IS cellIS; 1495 1496 PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL)); 1497 { 1498 PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3}; 1499 PetscWeakForm wf; 1500 PetscInt Nm = 4, m, Nk = 0, k, kp, off = 0; 1501 PetscFormKey *jackeys; 1502 1503 /* Get unique Jacobian keys */ 1504 for (m = 0; m < Nm; ++m) { 1505 PetscInt Nkm; 1506 PetscCall(PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm)); 1507 Nk += Nkm; 1508 } 1509 PetscCall(PetscMalloc1(Nk, &jackeys)); 1510 for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys)); 1511 PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk); 1512 PetscCall(PetscFormKeySort(Nk, jackeys)); 1513 for (k = 0, kp = 1; kp < Nk; ++kp) { 1514 if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) { 1515 ++k; 1516 if (kp != k) jackeys[k] = jackeys[kp]; 1517 } 1518 } 1519 Nk = k; 1520 1521 PetscCall(PetscDSGetWeakForm(ds, &wf)); 1522 for (k = 0; k < Nk; ++k) { 1523 DMLabel label = jackeys[k].label; 1524 PetscInt val = jackeys[k].value; 1525 1526 if (!label) { 1527 PetscCall(PetscObjectReference((PetscObject)allcellIS)); 1528 cellIS = allcellIS; 1529 } else { 1530 IS pointIS; 1531 1532 PetscCall(DMLabelGetStratumIS(label, val, &pointIS)); 1533 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS)); 1534 PetscCall(ISDestroy(&pointIS)); 1535 } 1536 PetscCall(DMPlexComputeJacobian_Action_Internal(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, user)); 1537 PetscCall(ISDestroy(&cellIS)); 1538 } 1539 PetscCall(PetscFree(jackeys)); 1540 } 1541 } 1542 PetscCall(ISDestroy(&allcellIS)); 1543 PetscCall(DMDestroy(&plex)); 1544 PetscFunctionReturn(PETSC_SUCCESS); 1545 } 1546 1547 /*@ 1548 DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix `Jac` at the local solution `X` using pointwise functions specified by the user. 1549 1550 Input Parameters: 1551 + dm - The `DM` 1552 . X - Local input vector 1553 - user - The user context 1554 1555 Output Parameters: 1556 + Jac - Jacobian matrix 1557 - JacP - approximate Jacobian from which the preconditioner will be built, often `Jac` 1558 1559 Level: developer 1560 1561 Note: 1562 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1563 like a GPU, or vectorize on a multicore machine. 1564 1565 .seealso: [](ch_snes), `DMPLEX`, `Mat` 1566 @*/ 1567 PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, void *user) 1568 { 1569 DM plex; 1570 IS allcellIS; 1571 PetscBool hasJac, hasPrec; 1572 PetscInt Nds, s; 1573 1574 PetscFunctionBegin; 1575 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1576 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS)); 1577 PetscCall(DMGetNumDS(dm, &Nds)); 1578 for (s = 0; s < Nds; ++s) { 1579 PetscDS ds; 1580 IS cellIS; 1581 PetscFormKey key; 1582 1583 PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL)); 1584 key.value = 0; 1585 key.field = 0; 1586 key.part = 0; 1587 if (!key.label) { 1588 PetscCall(PetscObjectReference((PetscObject)allcellIS)); 1589 cellIS = allcellIS; 1590 } else { 1591 IS pointIS; 1592 1593 key.value = 1; 1594 PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS)); 1595 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS)); 1596 PetscCall(ISDestroy(&pointIS)); 1597 } 1598 if (!s) { 1599 PetscCall(PetscDSHasJacobian(ds, &hasJac)); 1600 PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec)); 1601 if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac)); 1602 PetscCall(MatZeroEntries(JacP)); 1603 } 1604 PetscCall(DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user)); 1605 PetscCall(ISDestroy(&cellIS)); 1606 } 1607 PetscCall(ISDestroy(&allcellIS)); 1608 PetscCall(DMDestroy(&plex)); 1609 PetscFunctionReturn(PETSC_SUCCESS); 1610 } 1611 1612 struct _DMSNESJacobianMFCtx { 1613 DM dm; 1614 Vec X; 1615 void *ctx; 1616 }; 1617 1618 static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A) 1619 { 1620 struct _DMSNESJacobianMFCtx *ctx; 1621 1622 PetscFunctionBegin; 1623 PetscCall(MatShellGetContext(A, &ctx)); 1624 PetscCall(MatShellSetContext(A, NULL)); 1625 PetscCall(DMDestroy(&ctx->dm)); 1626 PetscCall(VecDestroy(&ctx->X)); 1627 PetscCall(PetscFree(ctx)); 1628 PetscFunctionReturn(PETSC_SUCCESS); 1629 } 1630 1631 static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z) 1632 { 1633 struct _DMSNESJacobianMFCtx *ctx; 1634 1635 PetscFunctionBegin; 1636 PetscCall(MatShellGetContext(A, &ctx)); 1637 PetscCall(DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx)); 1638 PetscFunctionReturn(PETSC_SUCCESS); 1639 } 1640 1641 /*@ 1642 DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free 1643 1644 Collective 1645 1646 Input Parameters: 1647 + dm - The `DM` 1648 . X - The evaluation point for the Jacobian 1649 - user - A user context, or `NULL` 1650 1651 Output Parameter: 1652 . J - The `Mat` 1653 1654 Level: advanced 1655 1656 Note: 1657 Vec `X` is kept in `J`, so updating `X` then updates the evaluation point. 1658 1659 .seealso: [](ch_snes), `DM`, `SNES`, `DMSNESComputeJacobianAction()` 1660 @*/ 1661 PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, void *user, Mat *J) 1662 { 1663 struct _DMSNESJacobianMFCtx *ctx; 1664 PetscInt n, N; 1665 1666 PetscFunctionBegin; 1667 PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J)); 1668 PetscCall(MatSetType(*J, MATSHELL)); 1669 PetscCall(VecGetLocalSize(X, &n)); 1670 PetscCall(VecGetSize(X, &N)); 1671 PetscCall(MatSetSizes(*J, n, n, N, N)); 1672 PetscCall(PetscObjectReference((PetscObject)dm)); 1673 PetscCall(PetscObjectReference((PetscObject)X)); 1674 PetscCall(PetscMalloc1(1, &ctx)); 1675 ctx->dm = dm; 1676 ctx->X = X; 1677 ctx->ctx = user; 1678 PetscCall(MatShellSetContext(*J, ctx)); 1679 PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))DMSNESJacobianMF_Destroy_Private)); 1680 PetscCall(MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))DMSNESJacobianMF_Mult_Private)); 1681 PetscFunctionReturn(PETSC_SUCCESS); 1682 } 1683 1684 static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx) 1685 { 1686 SNES snes; 1687 Mat pJ; 1688 DM ovldm, origdm; 1689 DMSNES sdm; 1690 PetscErrorCode (*bfun)(DM, Vec, void *); 1691 PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *); 1692 void *bctx, *jctx; 1693 1694 PetscFunctionBegin; 1695 PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ)); 1696 PetscCheck(pJ, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing overlapping Mat"); 1697 PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm)); 1698 PetscCheck(origdm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing original DM"); 1699 PetscCall(MatGetDM(pJ, &ovldm)); 1700 PetscCall(DMSNESGetBoundaryLocal(origdm, &bfun, &bctx)); 1701 PetscCall(DMSNESSetBoundaryLocal(ovldm, bfun, bctx)); 1702 PetscCall(DMSNESGetJacobianLocal(origdm, &jfun, &jctx)); 1703 PetscCall(DMSNESSetJacobianLocal(ovldm, jfun, jctx)); 1704 PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes)); 1705 if (!snes) { 1706 PetscCall(SNESCreate(PetscObjectComm((PetscObject)ovl), &snes)); 1707 PetscCall(SNESSetDM(snes, ovldm)); 1708 PetscCall(PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes)); 1709 PetscCall(PetscObjectDereference((PetscObject)snes)); 1710 } 1711 PetscCall(DMGetDMSNES(ovldm, &sdm)); 1712 PetscCall(VecLockReadPush(X)); 1713 { 1714 void *ctx; 1715 PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *); 1716 PetscCall(DMSNESGetJacobian(ovldm, &J, &ctx)); 1717 PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx)); 1718 } 1719 PetscCall(VecLockReadPop(X)); 1720 /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */ 1721 { 1722 Mat locpJ; 1723 1724 PetscCall(MatISGetLocalMat(pJ, &locpJ)); 1725 PetscCall(MatCopy(locpJ, J, SAME_NONZERO_PATTERN)); 1726 } 1727 PetscFunctionReturn(PETSC_SUCCESS); 1728 } 1729 1730 /*@ 1731 DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, objective, residual, and Jacobian. 1732 1733 Input Parameters: 1734 + dm - The `DM` object 1735 . use_obj - Use the objective function callback 1736 - ctx - The user context that will be passed to pointwise evaluation routines 1737 1738 Level: developer 1739 1740 .seealso: [](ch_snes),`DMPLEX`, `SNES`, `PetscDSAddBoundary()`, `PetscDSSetObjective()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()` 1741 @*/ 1742 PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, PetscBool use_obj, void *ctx) 1743 { 1744 PetscBool useCeed; 1745 1746 PetscFunctionBegin; 1747 PetscCall(DMPlexGetUseCeed(dm, &useCeed)); 1748 PetscCall(DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, ctx)); 1749 if (use_obj) PetscCall(DMSNESSetObjectiveLocal(dm, DMPlexSNESComputeObjectiveFEM, ctx)); 1750 if (useCeed) { 1751 #ifdef PETSC_HAVE_LIBCEED 1752 PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualCEED, ctx)); 1753 #else 1754 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot use CEED traversals without LibCEED. Rerun configure with --download-ceed"); 1755 #endif 1756 } else PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, ctx)); 1757 PetscCall(DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, ctx)); 1758 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex)); 1759 PetscFunctionReturn(PETSC_SUCCESS); 1760 } 1761 1762 /*@C 1763 DMSNESCheckDiscretization - Check the discretization error of the exact solution 1764 1765 Input Parameters: 1766 + snes - the `SNES` object 1767 . dm - the `DM` 1768 . t - the time 1769 . u - a `DM` vector 1770 - tol - A tolerance for the check, or -1 to print the results instead 1771 1772 Output Parameter: 1773 . error - An array which holds the discretization error in each field, or `NULL` 1774 1775 Level: developer 1776 1777 Note: 1778 The user must call `PetscDSSetExactSolution()` beforehand 1779 1780 Developer Note: 1781 How is this related to `PetscConvEst`? 1782 1783 .seealso: [](ch_snes), `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()` 1784 @*/ 1785 PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[]) 1786 { 1787 PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 1788 void **ectxs; 1789 PetscReal *err; 1790 MPI_Comm comm; 1791 PetscInt Nf, f; 1792 1793 PetscFunctionBegin; 1794 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1795 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1796 PetscValidHeaderSpecific(u, VEC_CLASSID, 4); 1797 if (error) PetscAssertPointer(error, 6); 1798 1799 PetscCall(DMComputeExactSolution(dm, t, u, NULL)); 1800 PetscCall(VecViewFromOptions(u, NULL, "-vec_view")); 1801 1802 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 1803 PetscCall(DMGetNumFields(dm, &Nf)); 1804 PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err)); 1805 { 1806 PetscInt Nds, s; 1807 1808 PetscCall(DMGetNumDS(dm, &Nds)); 1809 for (s = 0; s < Nds; ++s) { 1810 PetscDS ds; 1811 DMLabel label; 1812 IS fieldIS; 1813 const PetscInt *fields; 1814 PetscInt dsNf, f; 1815 1816 PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL)); 1817 PetscCall(PetscDSGetNumFields(ds, &dsNf)); 1818 PetscCall(ISGetIndices(fieldIS, &fields)); 1819 for (f = 0; f < dsNf; ++f) { 1820 const PetscInt field = fields[f]; 1821 PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field])); 1822 } 1823 PetscCall(ISRestoreIndices(fieldIS, &fields)); 1824 } 1825 } 1826 if (Nf > 1) { 1827 PetscCall(DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err)); 1828 if (tol >= 0.0) { 1829 for (f = 0; f < Nf; ++f) PetscCheck(err[f] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g for field %" PetscInt_FMT " exceeds tolerance %g", (double)err[f], f, (double)tol); 1830 } else if (error) { 1831 for (f = 0; f < Nf; ++f) error[f] = err[f]; 1832 } else { 1833 PetscCall(PetscPrintf(comm, "L_2 Error: [")); 1834 for (f = 0; f < Nf; ++f) { 1835 if (f) PetscCall(PetscPrintf(comm, ", ")); 1836 PetscCall(PetscPrintf(comm, "%g", (double)err[f])); 1837 } 1838 PetscCall(PetscPrintf(comm, "]\n")); 1839 } 1840 } else { 1841 PetscCall(DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0])); 1842 if (tol >= 0.0) { 1843 PetscCheck(err[0] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double)err[0], (double)tol); 1844 } else if (error) { 1845 error[0] = err[0]; 1846 } else { 1847 PetscCall(PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0])); 1848 } 1849 } 1850 PetscCall(PetscFree3(exacts, ectxs, err)); 1851 PetscFunctionReturn(PETSC_SUCCESS); 1852 } 1853 1854 /*@C 1855 DMSNESCheckResidual - Check the residual of the exact solution 1856 1857 Input Parameters: 1858 + snes - the `SNES` object 1859 . dm - the `DM` 1860 . u - a `DM` vector 1861 - tol - A tolerance for the check, or -1 to print the results instead 1862 1863 Output Parameter: 1864 . residual - The residual norm of the exact solution, or `NULL` 1865 1866 Level: developer 1867 1868 .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()` 1869 @*/ 1870 PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual) 1871 { 1872 MPI_Comm comm; 1873 Vec r; 1874 PetscReal res; 1875 1876 PetscFunctionBegin; 1877 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1878 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1879 PetscValidHeaderSpecific(u, VEC_CLASSID, 3); 1880 if (residual) PetscAssertPointer(residual, 5); 1881 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 1882 PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL)); 1883 PetscCall(VecDuplicate(u, &r)); 1884 PetscCall(SNESComputeFunction(snes, u, r)); 1885 PetscCall(VecNorm(r, NORM_2, &res)); 1886 if (tol >= 0.0) { 1887 PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol); 1888 } else if (residual) { 1889 *residual = res; 1890 } else { 1891 PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res)); 1892 PetscCall(VecFilter(r, 1.0e-10)); 1893 PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual")); 1894 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_")); 1895 PetscCall(VecViewFromOptions(r, NULL, "-vec_view")); 1896 } 1897 PetscCall(VecDestroy(&r)); 1898 PetscFunctionReturn(PETSC_SUCCESS); 1899 } 1900 1901 /*@C 1902 DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test 1903 1904 Input Parameters: 1905 + snes - the `SNES` object 1906 . dm - the `DM` 1907 . u - a `DM` vector 1908 - tol - A tolerance for the check, or -1 to print the results instead 1909 1910 Output Parameters: 1911 + isLinear - Flag indicaing that the function looks linear, or `NULL` 1912 - convRate - The rate of convergence of the linear model, or `NULL` 1913 1914 Level: developer 1915 1916 .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()` 1917 @*/ 1918 PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate) 1919 { 1920 MPI_Comm comm; 1921 PetscDS ds; 1922 Mat J, M; 1923 MatNullSpace nullspace; 1924 PetscReal slope, intercept; 1925 PetscBool hasJac, hasPrec, isLin = PETSC_FALSE; 1926 1927 PetscFunctionBegin; 1928 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1929 if (dm) PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1930 if (u) PetscValidHeaderSpecific(u, VEC_CLASSID, 3); 1931 if (isLinear) PetscAssertPointer(isLinear, 5); 1932 if (convRate) PetscAssertPointer(convRate, 6); 1933 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 1934 if (!dm) PetscCall(SNESGetDM(snes, &dm)); 1935 if (u) PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL)); 1936 else PetscCall(SNESGetSolution(snes, &u)); 1937 /* Create and view matrices */ 1938 PetscCall(DMCreateMatrix(dm, &J)); 1939 PetscCall(DMGetDS(dm, &ds)); 1940 PetscCall(PetscDSHasJacobian(ds, &hasJac)); 1941 PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec)); 1942 if (hasJac && hasPrec) { 1943 PetscCall(DMCreateMatrix(dm, &M)); 1944 PetscCall(SNESComputeJacobian(snes, u, J, M)); 1945 PetscCall(PetscObjectSetName((PetscObject)M, "Preconditioning Matrix")); 1946 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_")); 1947 PetscCall(MatViewFromOptions(M, NULL, "-mat_view")); 1948 PetscCall(MatDestroy(&M)); 1949 } else { 1950 PetscCall(SNESComputeJacobian(snes, u, J, J)); 1951 } 1952 PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian")); 1953 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_")); 1954 PetscCall(MatViewFromOptions(J, NULL, "-mat_view")); 1955 /* Check nullspace */ 1956 PetscCall(MatGetNullSpace(J, &nullspace)); 1957 if (nullspace) { 1958 PetscBool isNull; 1959 PetscCall(MatNullSpaceTest(nullspace, J, &isNull)); 1960 PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid."); 1961 } 1962 /* Taylor test */ 1963 { 1964 PetscRandom rand; 1965 Vec du, uhat, r, rhat, df; 1966 PetscReal h; 1967 PetscReal *es, *hs, *errors; 1968 PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1; 1969 PetscInt Nv, v; 1970 1971 /* Choose a perturbation direction */ 1972 PetscCall(PetscRandomCreate(comm, &rand)); 1973 PetscCall(VecDuplicate(u, &du)); 1974 PetscCall(VecSetRandom(du, rand)); 1975 PetscCall(PetscRandomDestroy(&rand)); 1976 PetscCall(VecDuplicate(u, &df)); 1977 PetscCall(MatMult(J, du, df)); 1978 /* Evaluate residual at u, F(u), save in vector r */ 1979 PetscCall(VecDuplicate(u, &r)); 1980 PetscCall(SNESComputeFunction(snes, u, r)); 1981 /* Look at the convergence of our Taylor approximation as we approach u */ 1982 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) 1983 ; 1984 PetscCall(PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors)); 1985 PetscCall(VecDuplicate(u, &uhat)); 1986 PetscCall(VecDuplicate(u, &rhat)); 1987 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) { 1988 PetscCall(VecWAXPY(uhat, h, du, u)); 1989 /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */ 1990 PetscCall(SNESComputeFunction(snes, uhat, rhat)); 1991 PetscCall(VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df)); 1992 PetscCall(VecNorm(rhat, NORM_2, &errors[Nv])); 1993 1994 es[Nv] = errors[Nv] == 0 ? -16.0 : PetscLog10Real(errors[Nv]); 1995 hs[Nv] = PetscLog10Real(h); 1996 } 1997 PetscCall(VecDestroy(&uhat)); 1998 PetscCall(VecDestroy(&rhat)); 1999 PetscCall(VecDestroy(&df)); 2000 PetscCall(VecDestroy(&r)); 2001 PetscCall(VecDestroy(&du)); 2002 for (v = 0; v < Nv; ++v) { 2003 if ((tol >= 0) && (errors[v] > tol)) break; 2004 else if (errors[v] > PETSC_SMALL) break; 2005 } 2006 if (v == Nv) isLin = PETSC_TRUE; 2007 PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept)); 2008 PetscCall(PetscFree3(es, hs, errors)); 2009 /* Slope should be about 2 */ 2010 if (tol >= 0) { 2011 PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope); 2012 } else if (isLinear || convRate) { 2013 if (isLinear) *isLinear = isLin; 2014 if (convRate) *convRate = slope; 2015 } else { 2016 if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope)); 2017 else PetscCall(PetscPrintf(comm, "Function appears to be linear\n")); 2018 } 2019 } 2020 PetscCall(MatDestroy(&J)); 2021 PetscFunctionReturn(PETSC_SUCCESS); 2022 } 2023 2024 PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u) 2025 { 2026 PetscFunctionBegin; 2027 PetscCall(DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL)); 2028 PetscCall(DMSNESCheckResidual(snes, dm, u, -1.0, NULL)); 2029 PetscCall(DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL)); 2030 PetscFunctionReturn(PETSC_SUCCESS); 2031 } 2032 2033 /*@C 2034 DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information 2035 2036 Input Parameters: 2037 + snes - the `SNES` object 2038 - u - representative `SNES` vector 2039 2040 Level: developer 2041 2042 Note: 2043 The user must call `PetscDSSetExactSolution()` before this call 2044 2045 .seealso: [](ch_snes), `SNES`, `DM` 2046 @*/ 2047 PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u) 2048 { 2049 DM dm; 2050 Vec sol; 2051 PetscBool check; 2052 2053 PetscFunctionBegin; 2054 PetscCall(PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check)); 2055 if (!check) PetscFunctionReturn(PETSC_SUCCESS); 2056 PetscCall(SNESGetDM(snes, &dm)); 2057 PetscCall(VecDuplicate(u, &sol)); 2058 PetscCall(SNESSetSolution(snes, sol)); 2059 PetscCall(DMSNESCheck_Internal(snes, dm, sol)); 2060 PetscCall(VecDestroy(&sol)); 2061 PetscFunctionReturn(PETSC_SUCCESS); 2062 } 2063