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, Vec xLocal, Vec v) 524 { 525 PetscReal v0, J, invJ, detJ; 526 const PetscInt dof = ctx->dof; 527 const PetscScalar *coords; 528 PetscScalar *a; 529 PetscInt p; 530 531 PetscFunctionBegin; 532 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 533 PetscCall(VecGetArray(v, &a)); 534 for (p = 0; p < ctx->n; ++p) { 535 PetscInt c = ctx->cells[p]; 536 PetscScalar *x = NULL; 537 PetscReal xir[1]; 538 PetscInt xSize, comp; 539 540 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ)); 541 PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c); 542 xir[0] = invJ * PetscRealPart(coords[p] - v0); 543 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x)); 544 if (2 * dof == xSize) { 545 for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) + x[1 * dof + comp] * xir[0]; 546 } else if (dof == xSize) { 547 for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp]; 548 } 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); 549 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x)); 550 } 551 PetscCall(VecRestoreArray(v, &a)); 552 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 553 PetscFunctionReturn(PETSC_SUCCESS); 554 } 555 556 static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v) 557 { 558 PetscReal *v0, *J, *invJ, detJ; 559 const PetscScalar *coords; 560 PetscScalar *a; 561 PetscInt p; 562 563 PetscFunctionBegin; 564 PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ)); 565 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 566 PetscCall(VecGetArray(v, &a)); 567 for (p = 0; p < ctx->n; ++p) { 568 PetscInt c = ctx->cells[p]; 569 PetscScalar *x = NULL; 570 PetscReal xi[4]; 571 PetscInt d, f, comp; 572 573 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 574 PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c); 575 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x)); 576 for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp]; 577 578 for (d = 0; d < ctx->dim; ++d) { 579 xi[d] = 0.0; 580 for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]); 581 for (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]; 582 } 583 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x)); 584 } 585 PetscCall(VecRestoreArray(v, &a)); 586 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 587 PetscCall(PetscFree3(v0, J, invJ)); 588 PetscFunctionReturn(PETSC_SUCCESS); 589 } 590 591 static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v) 592 { 593 PetscReal *v0, *J, *invJ, detJ; 594 const PetscScalar *coords; 595 PetscScalar *a; 596 PetscInt p; 597 598 PetscFunctionBegin; 599 PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ)); 600 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 601 PetscCall(VecGetArray(v, &a)); 602 for (p = 0; p < ctx->n; ++p) { 603 PetscInt c = ctx->cells[p]; 604 const PetscInt order[3] = {2, 1, 3}; 605 PetscScalar *x = NULL; 606 PetscReal xi[4]; 607 PetscInt d, f, comp; 608 609 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 610 PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c); 611 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x)); 612 for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp]; 613 614 for (d = 0; d < ctx->dim; ++d) { 615 xi[d] = 0.0; 616 for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]); 617 for (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]; 618 } 619 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x)); 620 } 621 PetscCall(VecRestoreArray(v, &a)); 622 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 623 PetscCall(PetscFree3(v0, J, invJ)); 624 PetscFunctionReturn(PETSC_SUCCESS); 625 } 626 627 static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx) 628 { 629 const PetscScalar *vertices = (const PetscScalar *)ctx; 630 const PetscScalar x0 = vertices[0]; 631 const PetscScalar y0 = vertices[1]; 632 const PetscScalar x1 = vertices[2]; 633 const PetscScalar y1 = vertices[3]; 634 const PetscScalar x2 = vertices[4]; 635 const PetscScalar y2 = vertices[5]; 636 const PetscScalar x3 = vertices[6]; 637 const PetscScalar y3 = vertices[7]; 638 const PetscScalar f_1 = x1 - x0; 639 const PetscScalar g_1 = y1 - y0; 640 const PetscScalar f_3 = x3 - x0; 641 const PetscScalar g_3 = y3 - y0; 642 const PetscScalar f_01 = x2 - x1 - x3 + x0; 643 const PetscScalar g_01 = y2 - y1 - y3 + y0; 644 const PetscScalar *ref; 645 PetscScalar *real; 646 647 PetscFunctionBegin; 648 PetscCall(VecGetArrayRead(Xref, &ref)); 649 PetscCall(VecGetArray(Xreal, &real)); 650 { 651 const PetscScalar p0 = ref[0]; 652 const PetscScalar p1 = ref[1]; 653 654 real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1; 655 real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1; 656 } 657 PetscCall(PetscLogFlops(28)); 658 PetscCall(VecRestoreArrayRead(Xref, &ref)); 659 PetscCall(VecRestoreArray(Xreal, &real)); 660 PetscFunctionReturn(PETSC_SUCCESS); 661 } 662 663 #include <petsc/private/dmimpl.h> 664 static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx) 665 { 666 const PetscScalar *vertices = (const PetscScalar *)ctx; 667 const PetscScalar x0 = vertices[0]; 668 const PetscScalar y0 = vertices[1]; 669 const PetscScalar x1 = vertices[2]; 670 const PetscScalar y1 = vertices[3]; 671 const PetscScalar x2 = vertices[4]; 672 const PetscScalar y2 = vertices[5]; 673 const PetscScalar x3 = vertices[6]; 674 const PetscScalar y3 = vertices[7]; 675 const PetscScalar f_01 = x2 - x1 - x3 + x0; 676 const PetscScalar g_01 = y2 - y1 - y3 + y0; 677 const PetscScalar *ref; 678 679 PetscFunctionBegin; 680 PetscCall(VecGetArrayRead(Xref, &ref)); 681 { 682 const PetscScalar x = ref[0]; 683 const PetscScalar y = ref[1]; 684 const PetscInt rows[2] = {0, 1}; 685 PetscScalar values[4]; 686 687 values[0] = (x1 - x0 + f_01 * y) * 0.5; 688 values[1] = (x3 - x0 + f_01 * x) * 0.5; 689 values[2] = (y1 - y0 + g_01 * y) * 0.5; 690 values[3] = (y3 - y0 + g_01 * x) * 0.5; 691 PetscCall(MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES)); 692 } 693 PetscCall(PetscLogFlops(30)); 694 PetscCall(VecRestoreArrayRead(Xref, &ref)); 695 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 696 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 697 PetscFunctionReturn(PETSC_SUCCESS); 698 } 699 700 static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v) 701 { 702 DM dmCoord; 703 PetscFE fem = NULL; 704 SNES snes; 705 KSP ksp; 706 PC pc; 707 Vec coordsLocal, r, ref, real; 708 Mat J; 709 PetscTabulation T = NULL; 710 const PetscScalar *coords; 711 PetscScalar *a; 712 PetscReal xir[2] = {0., 0.}; 713 PetscInt Nf, p; 714 const PetscInt dof = ctx->dof; 715 716 PetscFunctionBegin; 717 PetscCall(DMGetNumFields(dm, &Nf)); 718 if (Nf) { 719 PetscObject obj; 720 PetscClassId id; 721 722 PetscCall(DMGetField(dm, 0, NULL, &obj)); 723 PetscCall(PetscObjectGetClassId(obj, &id)); 724 if (id == PETSCFE_CLASSID) { 725 fem = (PetscFE)obj; 726 PetscCall(PetscFECreateTabulation(fem, 1, 1, xir, 0, &T)); 727 } 728 } 729 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 730 PetscCall(DMGetCoordinateDM(dm, &dmCoord)); 731 PetscCall(SNESCreate(PETSC_COMM_SELF, &snes)); 732 PetscCall(SNESSetOptionsPrefix(snes, "quad_interp_")); 733 PetscCall(VecCreate(PETSC_COMM_SELF, &r)); 734 PetscCall(VecSetSizes(r, 2, 2)); 735 PetscCall(VecSetType(r, dm->vectype)); 736 PetscCall(VecDuplicate(r, &ref)); 737 PetscCall(VecDuplicate(r, &real)); 738 PetscCall(MatCreate(PETSC_COMM_SELF, &J)); 739 PetscCall(MatSetSizes(J, 2, 2, 2, 2)); 740 PetscCall(MatSetType(J, MATSEQDENSE)); 741 PetscCall(MatSetUp(J)); 742 PetscCall(SNESSetFunction(snes, r, QuadMap_Private, NULL)); 743 PetscCall(SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL)); 744 PetscCall(SNESGetKSP(snes, &ksp)); 745 PetscCall(KSPGetPC(ksp, &pc)); 746 PetscCall(PCSetType(pc, PCLU)); 747 PetscCall(SNESSetFromOptions(snes)); 748 749 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 750 PetscCall(VecGetArray(v, &a)); 751 for (p = 0; p < ctx->n; ++p) { 752 PetscScalar *x = NULL, *vertices = NULL; 753 PetscScalar *xi; 754 PetscInt c = ctx->cells[p], comp, coordSize, xSize; 755 756 /* Can make this do all points at once */ 757 PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 758 PetscCheck(4 * 2 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %" PetscInt_FMT " should be %d", coordSize, 4 * 2); 759 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x)); 760 PetscCall(SNESSetFunction(snes, NULL, NULL, vertices)); 761 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices)); 762 PetscCall(VecGetArray(real, &xi)); 763 xi[0] = coords[p * ctx->dim + 0]; 764 xi[1] = coords[p * ctx->dim + 1]; 765 PetscCall(VecRestoreArray(real, &xi)); 766 PetscCall(SNESSolve(snes, real, ref)); 767 PetscCall(VecGetArray(ref, &xi)); 768 xir[0] = PetscRealPart(xi[0]); 769 xir[1] = PetscRealPart(xi[1]); 770 if (4 * dof == xSize) { 771 for (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]; 772 } else if (dof == xSize) { 773 for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp]; 774 } else { 775 PetscInt d; 776 777 PetscCheck(fem, ctx->comm, PETSC_ERR_ARG_WRONG, "Cannot have a higher order interpolant if the discretization is not PetscFE"); 778 xir[0] = 2.0 * xir[0] - 1.0; 779 xir[1] = 2.0 * xir[1] - 1.0; 780 PetscCall(PetscFEComputeTabulation(fem, 1, xir, 0, T)); 781 for (comp = 0; comp < dof; ++comp) { 782 a[p * dof + comp] = 0.0; 783 for (d = 0; d < xSize / dof; ++d) a[p * dof + comp] += x[d * dof + comp] * T->T[0][d * dof + comp]; 784 } 785 } 786 PetscCall(VecRestoreArray(ref, &xi)); 787 PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 788 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x)); 789 } 790 PetscCall(PetscTabulationDestroy(&T)); 791 PetscCall(VecRestoreArray(v, &a)); 792 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 793 794 PetscCall(SNESDestroy(&snes)); 795 PetscCall(VecDestroy(&r)); 796 PetscCall(VecDestroy(&ref)); 797 PetscCall(VecDestroy(&real)); 798 PetscCall(MatDestroy(&J)); 799 PetscFunctionReturn(PETSC_SUCCESS); 800 } 801 802 static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx) 803 { 804 const PetscScalar *vertices = (const PetscScalar *)ctx; 805 const PetscScalar x0 = vertices[0]; 806 const PetscScalar y0 = vertices[1]; 807 const PetscScalar z0 = vertices[2]; 808 const PetscScalar x1 = vertices[9]; 809 const PetscScalar y1 = vertices[10]; 810 const PetscScalar z1 = vertices[11]; 811 const PetscScalar x2 = vertices[6]; 812 const PetscScalar y2 = vertices[7]; 813 const PetscScalar z2 = vertices[8]; 814 const PetscScalar x3 = vertices[3]; 815 const PetscScalar y3 = vertices[4]; 816 const PetscScalar z3 = vertices[5]; 817 const PetscScalar x4 = vertices[12]; 818 const PetscScalar y4 = vertices[13]; 819 const PetscScalar z4 = vertices[14]; 820 const PetscScalar x5 = vertices[15]; 821 const PetscScalar y5 = vertices[16]; 822 const PetscScalar z5 = vertices[17]; 823 const PetscScalar x6 = vertices[18]; 824 const PetscScalar y6 = vertices[19]; 825 const PetscScalar z6 = vertices[20]; 826 const PetscScalar x7 = vertices[21]; 827 const PetscScalar y7 = vertices[22]; 828 const PetscScalar z7 = vertices[23]; 829 const PetscScalar f_1 = x1 - x0; 830 const PetscScalar g_1 = y1 - y0; 831 const PetscScalar h_1 = z1 - z0; 832 const PetscScalar f_3 = x3 - x0; 833 const PetscScalar g_3 = y3 - y0; 834 const PetscScalar h_3 = z3 - z0; 835 const PetscScalar f_4 = x4 - x0; 836 const PetscScalar g_4 = y4 - y0; 837 const PetscScalar h_4 = z4 - z0; 838 const PetscScalar f_01 = x2 - x1 - x3 + x0; 839 const PetscScalar g_01 = y2 - y1 - y3 + y0; 840 const PetscScalar h_01 = z2 - z1 - z3 + z0; 841 const PetscScalar f_12 = x7 - x3 - x4 + x0; 842 const PetscScalar g_12 = y7 - y3 - y4 + y0; 843 const PetscScalar h_12 = z7 - z3 - z4 + z0; 844 const PetscScalar f_02 = x5 - x1 - x4 + x0; 845 const PetscScalar g_02 = y5 - y1 - y4 + y0; 846 const PetscScalar h_02 = z5 - z1 - z4 + z0; 847 const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7; 848 const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7; 849 const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7; 850 const PetscScalar *ref; 851 PetscScalar *real; 852 853 PetscFunctionBegin; 854 PetscCall(VecGetArrayRead(Xref, &ref)); 855 PetscCall(VecGetArray(Xreal, &real)); 856 { 857 const PetscScalar p0 = ref[0]; 858 const PetscScalar p1 = ref[1]; 859 const PetscScalar p2 = ref[2]; 860 861 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; 862 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; 863 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; 864 } 865 PetscCall(PetscLogFlops(114)); 866 PetscCall(VecRestoreArrayRead(Xref, &ref)); 867 PetscCall(VecRestoreArray(Xreal, &real)); 868 PetscFunctionReturn(PETSC_SUCCESS); 869 } 870 871 static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx) 872 { 873 const PetscScalar *vertices = (const PetscScalar *)ctx; 874 const PetscScalar x0 = vertices[0]; 875 const PetscScalar y0 = vertices[1]; 876 const PetscScalar z0 = vertices[2]; 877 const PetscScalar x1 = vertices[9]; 878 const PetscScalar y1 = vertices[10]; 879 const PetscScalar z1 = vertices[11]; 880 const PetscScalar x2 = vertices[6]; 881 const PetscScalar y2 = vertices[7]; 882 const PetscScalar z2 = vertices[8]; 883 const PetscScalar x3 = vertices[3]; 884 const PetscScalar y3 = vertices[4]; 885 const PetscScalar z3 = vertices[5]; 886 const PetscScalar x4 = vertices[12]; 887 const PetscScalar y4 = vertices[13]; 888 const PetscScalar z4 = vertices[14]; 889 const PetscScalar x5 = vertices[15]; 890 const PetscScalar y5 = vertices[16]; 891 const PetscScalar z5 = vertices[17]; 892 const PetscScalar x6 = vertices[18]; 893 const PetscScalar y6 = vertices[19]; 894 const PetscScalar z6 = vertices[20]; 895 const PetscScalar x7 = vertices[21]; 896 const PetscScalar y7 = vertices[22]; 897 const PetscScalar z7 = vertices[23]; 898 const PetscScalar f_xy = x2 - x1 - x3 + x0; 899 const PetscScalar g_xy = y2 - y1 - y3 + y0; 900 const PetscScalar h_xy = z2 - z1 - z3 + z0; 901 const PetscScalar f_yz = x7 - x3 - x4 + x0; 902 const PetscScalar g_yz = y7 - y3 - y4 + y0; 903 const PetscScalar h_yz = z7 - z3 - z4 + z0; 904 const PetscScalar f_xz = x5 - x1 - x4 + x0; 905 const PetscScalar g_xz = y5 - y1 - y4 + y0; 906 const PetscScalar h_xz = z5 - z1 - z4 + z0; 907 const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7; 908 const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7; 909 const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7; 910 const PetscScalar *ref; 911 912 PetscFunctionBegin; 913 PetscCall(VecGetArrayRead(Xref, &ref)); 914 { 915 const PetscScalar x = ref[0]; 916 const PetscScalar y = ref[1]; 917 const PetscScalar z = ref[2]; 918 const PetscInt rows[3] = {0, 1, 2}; 919 PetscScalar values[9]; 920 921 values[0] = (x1 - x0 + f_xy * y + f_xz * z + f_xyz * y * z) / 2.0; 922 values[1] = (x3 - x0 + f_xy * x + f_yz * z + f_xyz * x * z) / 2.0; 923 values[2] = (x4 - x0 + f_yz * y + f_xz * x + f_xyz * x * y) / 2.0; 924 values[3] = (y1 - y0 + g_xy * y + g_xz * z + g_xyz * y * z) / 2.0; 925 values[4] = (y3 - y0 + g_xy * x + g_yz * z + g_xyz * x * z) / 2.0; 926 values[5] = (y4 - y0 + g_yz * y + g_xz * x + g_xyz * x * y) / 2.0; 927 values[6] = (z1 - z0 + h_xy * y + h_xz * z + h_xyz * y * z) / 2.0; 928 values[7] = (z3 - z0 + h_xy * x + h_yz * z + h_xyz * x * z) / 2.0; 929 values[8] = (z4 - z0 + h_yz * y + h_xz * x + h_xyz * x * y) / 2.0; 930 931 PetscCall(MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES)); 932 } 933 PetscCall(PetscLogFlops(152)); 934 PetscCall(VecRestoreArrayRead(Xref, &ref)); 935 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 936 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 937 PetscFunctionReturn(PETSC_SUCCESS); 938 } 939 940 static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v) 941 { 942 DM dmCoord; 943 SNES snes; 944 KSP ksp; 945 PC pc; 946 Vec coordsLocal, r, ref, real; 947 Mat J; 948 const PetscScalar *coords; 949 PetscScalar *a; 950 PetscInt p; 951 952 PetscFunctionBegin; 953 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 954 PetscCall(DMGetCoordinateDM(dm, &dmCoord)); 955 PetscCall(SNESCreate(PETSC_COMM_SELF, &snes)); 956 PetscCall(SNESSetOptionsPrefix(snes, "hex_interp_")); 957 PetscCall(VecCreate(PETSC_COMM_SELF, &r)); 958 PetscCall(VecSetSizes(r, 3, 3)); 959 PetscCall(VecSetType(r, dm->vectype)); 960 PetscCall(VecDuplicate(r, &ref)); 961 PetscCall(VecDuplicate(r, &real)); 962 PetscCall(MatCreate(PETSC_COMM_SELF, &J)); 963 PetscCall(MatSetSizes(J, 3, 3, 3, 3)); 964 PetscCall(MatSetType(J, MATSEQDENSE)); 965 PetscCall(MatSetUp(J)); 966 PetscCall(SNESSetFunction(snes, r, HexMap_Private, NULL)); 967 PetscCall(SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL)); 968 PetscCall(SNESGetKSP(snes, &ksp)); 969 PetscCall(KSPGetPC(ksp, &pc)); 970 PetscCall(PCSetType(pc, PCLU)); 971 PetscCall(SNESSetFromOptions(snes)); 972 973 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 974 PetscCall(VecGetArray(v, &a)); 975 for (p = 0; p < ctx->n; ++p) { 976 PetscScalar *x = NULL, *vertices = NULL; 977 PetscScalar *xi; 978 PetscReal xir[3]; 979 PetscInt c = ctx->cells[p], comp, coordSize, xSize; 980 981 /* Can make this do all points at once */ 982 PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 983 PetscCheck(8 * 3 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid coordinate closure size %" PetscInt_FMT " should be %d", coordSize, 8 * 3); 984 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x)); 985 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); 986 PetscCall(SNESSetFunction(snes, NULL, NULL, vertices)); 987 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices)); 988 PetscCall(VecGetArray(real, &xi)); 989 xi[0] = coords[p * ctx->dim + 0]; 990 xi[1] = coords[p * ctx->dim + 1]; 991 xi[2] = coords[p * ctx->dim + 2]; 992 PetscCall(VecRestoreArray(real, &xi)); 993 PetscCall(SNESSolve(snes, real, ref)); 994 PetscCall(VecGetArray(ref, &xi)); 995 xir[0] = PetscRealPart(xi[0]); 996 xir[1] = PetscRealPart(xi[1]); 997 xir[2] = PetscRealPart(xi[2]); 998 if (8 * ctx->dof == xSize) { 999 for (comp = 0; comp < ctx->dof; ++comp) { 1000 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]) + 1001 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]; 1002 } 1003 } else { 1004 for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp]; 1005 } 1006 PetscCall(VecRestoreArray(ref, &xi)); 1007 PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 1008 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x)); 1009 } 1010 PetscCall(VecRestoreArray(v, &a)); 1011 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 1012 1013 PetscCall(SNESDestroy(&snes)); 1014 PetscCall(VecDestroy(&r)); 1015 PetscCall(VecDestroy(&ref)); 1016 PetscCall(VecDestroy(&real)); 1017 PetscCall(MatDestroy(&J)); 1018 PetscFunctionReturn(PETSC_SUCCESS); 1019 } 1020 1021 /*@C 1022 DMInterpolationEvaluate - Using the input from dm and x, calculates interpolated field values at the interpolation points. 1023 1024 Input Parameters: 1025 + ctx - The `DMInterpolationInfo` context 1026 . dm - The `DM` 1027 - x - The local vector containing the field to be interpolated 1028 1029 Output Parameter: 1030 . v - The vector containing the interpolated values 1031 1032 Level: beginner 1033 1034 Note: 1035 A suitable `v` can be obtained using `DMInterpolationGetVector()`. 1036 1037 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 1038 @*/ 1039 PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v) 1040 { 1041 PetscDS ds; 1042 PetscInt n, p, Nf, field; 1043 PetscBool useDS = PETSC_FALSE; 1044 1045 PetscFunctionBegin; 1046 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1047 PetscValidHeaderSpecific(x, VEC_CLASSID, 3); 1048 PetscValidHeaderSpecific(v, VEC_CLASSID, 4); 1049 PetscCall(VecGetLocalSize(v, &n)); 1050 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); 1051 if (!n) PetscFunctionReturn(PETSC_SUCCESS); 1052 PetscCall(DMGetDS(dm, &ds)); 1053 if (ds) { 1054 useDS = PETSC_TRUE; 1055 PetscCall(PetscDSGetNumFields(ds, &Nf)); 1056 for (field = 0; field < Nf; ++field) { 1057 PetscObject obj; 1058 PetscClassId id; 1059 1060 PetscCall(PetscDSGetDiscretization(ds, field, &obj)); 1061 PetscCall(PetscObjectGetClassId(obj, &id)); 1062 if (id != PETSCFE_CLASSID) { 1063 useDS = PETSC_FALSE; 1064 break; 1065 } 1066 } 1067 } 1068 if (useDS) { 1069 const PetscScalar *coords; 1070 PetscScalar *interpolant; 1071 PetscInt cdim, d; 1072 1073 PetscCall(DMGetCoordinateDim(dm, &cdim)); 1074 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 1075 PetscCall(VecGetArrayWrite(v, &interpolant)); 1076 for (p = 0; p < ctx->n; ++p) { 1077 PetscReal pcoords[3], xi[3]; 1078 PetscScalar *xa = NULL; 1079 PetscInt coff = 0, foff = 0, clSize; 1080 1081 if (ctx->cells[p] < 0) continue; 1082 for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]); 1083 PetscCall(DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi)); 1084 PetscCall(DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa)); 1085 for (field = 0; field < Nf; ++field) { 1086 PetscTabulation T; 1087 PetscFE fe; 1088 1089 PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe)); 1090 PetscCall(PetscFECreateTabulation(fe, 1, 1, xi, 0, &T)); 1091 { 1092 const PetscReal *basis = T->T[0]; 1093 const PetscInt Nb = T->Nb; 1094 const PetscInt Nc = T->Nc; 1095 PetscInt f, fc; 1096 1097 for (fc = 0; fc < Nc; ++fc) { 1098 interpolant[p * ctx->dof + coff + fc] = 0.0; 1099 for (f = 0; f < Nb; ++f) interpolant[p * ctx->dof + coff + fc] += xa[foff + f] * basis[(0 * Nb + f) * Nc + fc]; 1100 } 1101 coff += Nc; 1102 foff += Nb; 1103 } 1104 PetscCall(PetscTabulationDestroy(&T)); 1105 } 1106 PetscCall(DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa)); 1107 PetscCheck(coff == ctx->dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total components %" PetscInt_FMT " != %" PetscInt_FMT " dof specified for interpolation", coff, ctx->dof); 1108 PetscCheck(foff == clSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total FE space size %" PetscInt_FMT " != %" PetscInt_FMT " closure size", foff, clSize); 1109 } 1110 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 1111 PetscCall(VecRestoreArrayWrite(v, &interpolant)); 1112 } else { 1113 DMPolytopeType ct; 1114 1115 /* TODO Check each cell individually */ 1116 PetscCall(DMPlexGetCellType(dm, ctx->cells[0], &ct)); 1117 switch (ct) { 1118 case DM_POLYTOPE_SEGMENT: 1119 PetscCall(DMInterpolate_Segment_Private(ctx, dm, x, v)); 1120 break; 1121 case DM_POLYTOPE_TRIANGLE: 1122 PetscCall(DMInterpolate_Triangle_Private(ctx, dm, x, v)); 1123 break; 1124 case DM_POLYTOPE_QUADRILATERAL: 1125 PetscCall(DMInterpolate_Quad_Private(ctx, dm, x, v)); 1126 break; 1127 case DM_POLYTOPE_TETRAHEDRON: 1128 PetscCall(DMInterpolate_Tetrahedron_Private(ctx, dm, x, v)); 1129 break; 1130 case DM_POLYTOPE_HEXAHEDRON: 1131 PetscCall(DMInterpolate_Hex_Private(ctx, dm, x, v)); 1132 break; 1133 default: 1134 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for cell type %s", DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]); 1135 } 1136 } 1137 PetscFunctionReturn(PETSC_SUCCESS); 1138 } 1139 1140 /*@C 1141 DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context 1142 1143 Collective 1144 1145 Input Parameter: 1146 . ctx - the context 1147 1148 Level: beginner 1149 1150 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 1151 @*/ 1152 PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx) 1153 { 1154 PetscFunctionBegin; 1155 PetscAssertPointer(ctx, 1); 1156 PetscCall(VecDestroy(&(*ctx)->coords)); 1157 PetscCall(PetscFree((*ctx)->points)); 1158 PetscCall(PetscFree((*ctx)->cells)); 1159 PetscCall(PetscFree(*ctx)); 1160 *ctx = NULL; 1161 PetscFunctionReturn(PETSC_SUCCESS); 1162 } 1163 1164 /*@C 1165 SNESMonitorFields - Monitors the residual for each field separately 1166 1167 Collective 1168 1169 Input Parameters: 1170 + snes - the `SNES` context, must have an attached `DM` 1171 . its - iteration number 1172 . fgnorm - 2-norm of residual 1173 - vf - `PetscViewerAndFormat` of `PetscViewerType` `PETSCVIEWERASCII` 1174 1175 Level: intermediate 1176 1177 Note: 1178 This routine prints the residual norm at each iteration. 1179 1180 .seealso: [](ch_snes), `SNES`, `SNESMonitorSet()`, `SNESMonitorDefault()` 1181 @*/ 1182 PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf) 1183 { 1184 PetscViewer viewer = vf->viewer; 1185 Vec res; 1186 DM dm; 1187 PetscSection s; 1188 const PetscScalar *r; 1189 PetscReal *lnorms, *norms; 1190 PetscInt numFields, f, pStart, pEnd, p; 1191 1192 PetscFunctionBegin; 1193 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 4); 1194 PetscCall(SNESGetFunction(snes, &res, NULL, NULL)); 1195 PetscCall(SNESGetDM(snes, &dm)); 1196 PetscCall(DMGetLocalSection(dm, &s)); 1197 PetscCall(PetscSectionGetNumFields(s, &numFields)); 1198 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 1199 PetscCall(PetscCalloc2(numFields, &lnorms, numFields, &norms)); 1200 PetscCall(VecGetArrayRead(res, &r)); 1201 for (p = pStart; p < pEnd; ++p) { 1202 for (f = 0; f < numFields; ++f) { 1203 PetscInt fdof, foff, d; 1204 1205 PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof)); 1206 PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff)); 1207 for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff + d])); 1208 } 1209 } 1210 PetscCall(VecRestoreArrayRead(res, &r)); 1211 PetscCall(MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm))); 1212 PetscCall(PetscViewerPushFormat(viewer, vf->format)); 1213 PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel)); 1214 PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " SNES Function norm %14.12e [", its, (double)fgnorm)); 1215 for (f = 0; f < numFields; ++f) { 1216 if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", ")); 1217 PetscCall(PetscViewerASCIIPrintf(viewer, "%14.12e", (double)PetscSqrtReal(norms[f]))); 1218 } 1219 PetscCall(PetscViewerASCIIPrintf(viewer, "]\n")); 1220 PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel)); 1221 PetscCall(PetscViewerPopFormat(viewer)); 1222 PetscCall(PetscFree2(lnorms, norms)); 1223 PetscFunctionReturn(PETSC_SUCCESS); 1224 } 1225 1226 /********************* Residual Computation **************************/ 1227 1228 /*@ 1229 DMPlexSNESComputeResidualFEM - Sums the local residual into vector `F` from the local input `X` using pointwise functions specified by the user 1230 1231 Input Parameters: 1232 + dm - The mesh 1233 . X - Local solution 1234 - user - The user context 1235 1236 Output Parameter: 1237 . F - Local output vector 1238 1239 Level: developer 1240 1241 Note: 1242 The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional. 1243 1244 .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()` 1245 @*/ 1246 PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user) 1247 { 1248 DM plex; 1249 IS allcellIS; 1250 PetscInt Nds, s; 1251 1252 PetscFunctionBegin; 1253 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1254 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS)); 1255 PetscCall(DMGetNumDS(dm, &Nds)); 1256 for (s = 0; s < Nds; ++s) { 1257 PetscDS ds; 1258 IS cellIS; 1259 PetscFormKey key; 1260 1261 PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL)); 1262 key.value = 0; 1263 key.field = 0; 1264 key.part = 0; 1265 if (!key.label) { 1266 PetscCall(PetscObjectReference((PetscObject)allcellIS)); 1267 cellIS = allcellIS; 1268 } else { 1269 IS pointIS; 1270 1271 key.value = 1; 1272 PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS)); 1273 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS)); 1274 PetscCall(ISDestroy(&pointIS)); 1275 } 1276 PetscCall(DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user)); 1277 PetscCall(ISDestroy(&cellIS)); 1278 } 1279 PetscCall(ISDestroy(&allcellIS)); 1280 PetscCall(DMDestroy(&plex)); 1281 PetscFunctionReturn(PETSC_SUCCESS); 1282 } 1283 1284 /*@ 1285 DMPlexSNESComputeResidualDS - Sums the local residual into vector `F` from the local input `X` using all pointwise functions with unique keys in the `PetscDS` 1286 1287 Input Parameters: 1288 + dm - The mesh 1289 . X - Local solution 1290 - user - The user context 1291 1292 Output Parameter: 1293 . F - Local output vector 1294 1295 Level: developer 1296 1297 Note: 1298 The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional. 1299 1300 .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()` 1301 @*/ 1302 PetscErrorCode DMPlexSNESComputeResidualDS(DM dm, Vec X, Vec F, void *user) 1303 { 1304 DM plex; 1305 IS allcellIS; 1306 PetscInt Nds, s; 1307 1308 PetscFunctionBegin; 1309 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1310 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS)); 1311 PetscCall(DMGetNumDS(dm, &Nds)); 1312 for (s = 0; s < Nds; ++s) { 1313 PetscDS ds; 1314 DMLabel label; 1315 IS cellIS; 1316 1317 PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL)); 1318 { 1319 PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1}; 1320 PetscWeakForm wf; 1321 PetscInt Nm = 2, m, Nk = 0, k, kp, off = 0; 1322 PetscFormKey *reskeys; 1323 1324 /* Get unique residual keys */ 1325 for (m = 0; m < Nm; ++m) { 1326 PetscInt Nkm; 1327 PetscCall(PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm)); 1328 Nk += Nkm; 1329 } 1330 PetscCall(PetscMalloc1(Nk, &reskeys)); 1331 for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys)); 1332 PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk); 1333 PetscCall(PetscFormKeySort(Nk, reskeys)); 1334 for (k = 0, kp = 1; kp < Nk; ++kp) { 1335 if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) { 1336 ++k; 1337 if (kp != k) reskeys[k] = reskeys[kp]; 1338 } 1339 } 1340 Nk = k; 1341 1342 PetscCall(PetscDSGetWeakForm(ds, &wf)); 1343 for (k = 0; k < Nk; ++k) { 1344 DMLabel label = reskeys[k].label; 1345 PetscInt val = reskeys[k].value; 1346 1347 if (!label) { 1348 PetscCall(PetscObjectReference((PetscObject)allcellIS)); 1349 cellIS = allcellIS; 1350 } else { 1351 IS pointIS; 1352 1353 PetscCall(DMLabelGetStratumIS(label, val, &pointIS)); 1354 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS)); 1355 PetscCall(ISDestroy(&pointIS)); 1356 } 1357 PetscCall(DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user)); 1358 PetscCall(ISDestroy(&cellIS)); 1359 } 1360 PetscCall(PetscFree(reskeys)); 1361 } 1362 } 1363 PetscCall(ISDestroy(&allcellIS)); 1364 PetscCall(DMDestroy(&plex)); 1365 PetscFunctionReturn(PETSC_SUCCESS); 1366 } 1367 1368 #ifdef PETSC_HAVE_LIBCEED 1369 PetscErrorCode DMPlexSNESComputeResidualCEED(DM dm, Vec locX, Vec locF, void *user) 1370 { 1371 Ceed ceed; 1372 DMCeed sd = dm->dmceed; 1373 CeedVector clocX, clocF; 1374 1375 PetscFunctionBegin; 1376 PetscCall(DMGetCeed(dm, &ceed)); 1377 PetscCheck(sd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This DM has no CEED data. Call DMCeedCreate() before computing the residual."); 1378 PetscCall(DMCeedComputeGeometry(dm, sd)); 1379 1380 PetscCall(VecGetCeedVectorRead(locX, ceed, &clocX)); 1381 PetscCall(VecGetCeedVector(locF, ceed, &clocF)); 1382 PetscCallCEED(CeedOperatorApplyAdd(sd->op, clocX, clocF, CEED_REQUEST_IMMEDIATE)); 1383 PetscCall(VecRestoreCeedVectorRead(locX, &clocX)); 1384 PetscCall(VecRestoreCeedVector(locF, &clocF)); 1385 1386 { 1387 DM_Plex *mesh = (DM_Plex *)dm->data; 1388 1389 if (mesh->printFEM) { 1390 PetscSection section; 1391 Vec locFbc; 1392 PetscInt pStart, pEnd, p, maxDof; 1393 PetscScalar *zeroes; 1394 1395 PetscCall(DMGetLocalSection(dm, §ion)); 1396 PetscCall(VecDuplicate(locF, &locFbc)); 1397 PetscCall(VecCopy(locF, locFbc)); 1398 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 1399 PetscCall(PetscSectionGetMaxDof(section, &maxDof)); 1400 PetscCall(PetscCalloc1(maxDof, &zeroes)); 1401 for (p = pStart; p < pEnd; ++p) PetscCall(VecSetValuesSection(locFbc, section, p, zeroes, INSERT_BC_VALUES)); 1402 PetscCall(PetscFree(zeroes)); 1403 PetscCall(DMPrintLocalVec(dm, "Residual", mesh->printTol, locFbc)); 1404 PetscCall(VecDestroy(&locFbc)); 1405 } 1406 } 1407 PetscFunctionReturn(PETSC_SUCCESS); 1408 } 1409 #endif 1410 1411 /*@ 1412 DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input `X` 1413 1414 Input Parameters: 1415 + dm - The mesh 1416 - user - The user context 1417 1418 Output Parameter: 1419 . X - Local solution 1420 1421 Level: developer 1422 1423 .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()` 1424 @*/ 1425 PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user) 1426 { 1427 DM plex; 1428 1429 PetscFunctionBegin; 1430 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1431 PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL)); 1432 PetscCall(DMDestroy(&plex)); 1433 PetscFunctionReturn(PETSC_SUCCESS); 1434 } 1435 1436 /*@ 1437 DMSNESComputeJacobianAction - Compute the action of the Jacobian J(X) on Y 1438 1439 Input Parameters: 1440 + dm - The `DM` 1441 . X - Local solution vector 1442 . Y - Local input vector 1443 - user - The user context 1444 1445 Output Parameter: 1446 . F - local output vector 1447 1448 Level: developer 1449 1450 Notes: 1451 Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly. 1452 1453 .seealso: [](ch_snes), `DM`, ``DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()` 1454 @*/ 1455 PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, void *user) 1456 { 1457 DM plex; 1458 IS allcellIS; 1459 PetscInt Nds, s; 1460 1461 PetscFunctionBegin; 1462 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1463 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS)); 1464 PetscCall(DMGetNumDS(dm, &Nds)); 1465 for (s = 0; s < Nds; ++s) { 1466 PetscDS ds; 1467 DMLabel label; 1468 IS cellIS; 1469 1470 PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL)); 1471 { 1472 PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3}; 1473 PetscWeakForm wf; 1474 PetscInt Nm = 4, m, Nk = 0, k, kp, off = 0; 1475 PetscFormKey *jackeys; 1476 1477 /* Get unique Jacobian keys */ 1478 for (m = 0; m < Nm; ++m) { 1479 PetscInt Nkm; 1480 PetscCall(PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm)); 1481 Nk += Nkm; 1482 } 1483 PetscCall(PetscMalloc1(Nk, &jackeys)); 1484 for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys)); 1485 PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk); 1486 PetscCall(PetscFormKeySort(Nk, jackeys)); 1487 for (k = 0, kp = 1; kp < Nk; ++kp) { 1488 if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) { 1489 ++k; 1490 if (kp != k) jackeys[k] = jackeys[kp]; 1491 } 1492 } 1493 Nk = k; 1494 1495 PetscCall(PetscDSGetWeakForm(ds, &wf)); 1496 for (k = 0; k < Nk; ++k) { 1497 DMLabel label = jackeys[k].label; 1498 PetscInt val = jackeys[k].value; 1499 1500 if (!label) { 1501 PetscCall(PetscObjectReference((PetscObject)allcellIS)); 1502 cellIS = allcellIS; 1503 } else { 1504 IS pointIS; 1505 1506 PetscCall(DMLabelGetStratumIS(label, val, &pointIS)); 1507 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS)); 1508 PetscCall(ISDestroy(&pointIS)); 1509 } 1510 PetscCall(DMPlexComputeJacobian_Action_Internal(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, user)); 1511 PetscCall(ISDestroy(&cellIS)); 1512 } 1513 PetscCall(PetscFree(jackeys)); 1514 } 1515 } 1516 PetscCall(ISDestroy(&allcellIS)); 1517 PetscCall(DMDestroy(&plex)); 1518 PetscFunctionReturn(PETSC_SUCCESS); 1519 } 1520 1521 /*@ 1522 DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix `Jac` at the local solution `X` using pointwise functions specified by the user. 1523 1524 Input Parameters: 1525 + dm - The `DM` 1526 . X - Local input vector 1527 - user - The user context 1528 1529 Output Parameters: 1530 + Jac - Jacobian matrix 1531 - JacP - approximate Jacobian from which the preconditioner will be built, often `Jac` 1532 1533 Level: developer 1534 1535 Note: 1536 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1537 like a GPU, or vectorize on a multicore machine. 1538 1539 .seealso: [](ch_snes), `DMPLEX`, `Mat` 1540 @*/ 1541 PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, void *user) 1542 { 1543 DM plex; 1544 IS allcellIS; 1545 PetscBool hasJac, hasPrec; 1546 PetscInt Nds, s; 1547 1548 PetscFunctionBegin; 1549 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1550 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS)); 1551 PetscCall(DMGetNumDS(dm, &Nds)); 1552 for (s = 0; s < Nds; ++s) { 1553 PetscDS ds; 1554 IS cellIS; 1555 PetscFormKey key; 1556 1557 PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL)); 1558 key.value = 0; 1559 key.field = 0; 1560 key.part = 0; 1561 if (!key.label) { 1562 PetscCall(PetscObjectReference((PetscObject)allcellIS)); 1563 cellIS = allcellIS; 1564 } else { 1565 IS pointIS; 1566 1567 key.value = 1; 1568 PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS)); 1569 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS)); 1570 PetscCall(ISDestroy(&pointIS)); 1571 } 1572 if (!s) { 1573 PetscCall(PetscDSHasJacobian(ds, &hasJac)); 1574 PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec)); 1575 if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac)); 1576 PetscCall(MatZeroEntries(JacP)); 1577 } 1578 PetscCall(DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user)); 1579 PetscCall(ISDestroy(&cellIS)); 1580 } 1581 PetscCall(ISDestroy(&allcellIS)); 1582 PetscCall(DMDestroy(&plex)); 1583 PetscFunctionReturn(PETSC_SUCCESS); 1584 } 1585 1586 struct _DMSNESJacobianMFCtx { 1587 DM dm; 1588 Vec X; 1589 void *ctx; 1590 }; 1591 1592 static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A) 1593 { 1594 struct _DMSNESJacobianMFCtx *ctx; 1595 1596 PetscFunctionBegin; 1597 PetscCall(MatShellGetContext(A, &ctx)); 1598 PetscCall(MatShellSetContext(A, NULL)); 1599 PetscCall(DMDestroy(&ctx->dm)); 1600 PetscCall(VecDestroy(&ctx->X)); 1601 PetscCall(PetscFree(ctx)); 1602 PetscFunctionReturn(PETSC_SUCCESS); 1603 } 1604 1605 static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z) 1606 { 1607 struct _DMSNESJacobianMFCtx *ctx; 1608 1609 PetscFunctionBegin; 1610 PetscCall(MatShellGetContext(A, &ctx)); 1611 PetscCall(DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx)); 1612 PetscFunctionReturn(PETSC_SUCCESS); 1613 } 1614 1615 /*@ 1616 DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free 1617 1618 Collective 1619 1620 Input Parameters: 1621 + dm - The `DM` 1622 . X - The evaluation point for the Jacobian 1623 - user - A user context, or `NULL` 1624 1625 Output Parameter: 1626 . J - The `Mat` 1627 1628 Level: advanced 1629 1630 Note: 1631 Vec `X` is kept in `J`, so updating `X` then updates the evaluation point. 1632 1633 .seealso: [](ch_snes), `DM`, `SNES`, `DMSNESComputeJacobianAction()` 1634 @*/ 1635 PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, void *user, Mat *J) 1636 { 1637 struct _DMSNESJacobianMFCtx *ctx; 1638 PetscInt n, N; 1639 1640 PetscFunctionBegin; 1641 PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J)); 1642 PetscCall(MatSetType(*J, MATSHELL)); 1643 PetscCall(VecGetLocalSize(X, &n)); 1644 PetscCall(VecGetSize(X, &N)); 1645 PetscCall(MatSetSizes(*J, n, n, N, N)); 1646 PetscCall(PetscObjectReference((PetscObject)dm)); 1647 PetscCall(PetscObjectReference((PetscObject)X)); 1648 PetscCall(PetscMalloc1(1, &ctx)); 1649 ctx->dm = dm; 1650 ctx->X = X; 1651 ctx->ctx = user; 1652 PetscCall(MatShellSetContext(*J, ctx)); 1653 PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))DMSNESJacobianMF_Destroy_Private)); 1654 PetscCall(MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))DMSNESJacobianMF_Mult_Private)); 1655 PetscFunctionReturn(PETSC_SUCCESS); 1656 } 1657 1658 /* 1659 MatComputeNeumannOverlap - Computes an unassembled (Neumann) local overlapping Mat in nonlinear context. 1660 1661 Input Parameters: 1662 + X - `SNES` linearization point 1663 . ovl - index set of overlapping subdomains 1664 1665 Output Parameter: 1666 . J - unassembled (Neumann) local matrix 1667 1668 Level: intermediate 1669 1670 .seealso: [](ch_snes), `DMCreateNeumannOverlap()`, `MATIS`, `PCHPDDMSetAuxiliaryMat()` 1671 */ 1672 static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx) 1673 { 1674 SNES snes; 1675 Mat pJ; 1676 DM ovldm, origdm; 1677 DMSNES sdm; 1678 PetscErrorCode (*bfun)(DM, Vec, void *); 1679 PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *); 1680 void *bctx, *jctx; 1681 1682 PetscFunctionBegin; 1683 PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ)); 1684 PetscCheck(pJ, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing overlapping Mat"); 1685 PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm)); 1686 PetscCheck(origdm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing original DM"); 1687 PetscCall(MatGetDM(pJ, &ovldm)); 1688 PetscCall(DMSNESGetBoundaryLocal(origdm, &bfun, &bctx)); 1689 PetscCall(DMSNESSetBoundaryLocal(ovldm, bfun, bctx)); 1690 PetscCall(DMSNESGetJacobianLocal(origdm, &jfun, &jctx)); 1691 PetscCall(DMSNESSetJacobianLocal(ovldm, jfun, jctx)); 1692 PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes)); 1693 if (!snes) { 1694 PetscCall(SNESCreate(PetscObjectComm((PetscObject)ovl), &snes)); 1695 PetscCall(SNESSetDM(snes, ovldm)); 1696 PetscCall(PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes)); 1697 PetscCall(PetscObjectDereference((PetscObject)snes)); 1698 } 1699 PetscCall(DMGetDMSNES(ovldm, &sdm)); 1700 PetscCall(VecLockReadPush(X)); 1701 { 1702 void *ctx; 1703 PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *); 1704 PetscCall(DMSNESGetJacobian(ovldm, &J, &ctx)); 1705 PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx)); 1706 } 1707 PetscCall(VecLockReadPop(X)); 1708 /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */ 1709 { 1710 Mat locpJ; 1711 1712 PetscCall(MatISGetLocalMat(pJ, &locpJ)); 1713 PetscCall(MatCopy(locpJ, J, SAME_NONZERO_PATTERN)); 1714 } 1715 PetscFunctionReturn(PETSC_SUCCESS); 1716 } 1717 1718 /*@ 1719 DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, residual, and Jacobian. 1720 1721 Input Parameters: 1722 + dm - The `DM` object 1723 . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see `PetscDSAddBoundary()`) 1724 . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see `PetscDSSetResidual()`) 1725 - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see `PetscDSSetJacobian()`) 1726 1727 Level: developer 1728 1729 .seealso: [](ch_snes), `DMPLEX`, `SNES` 1730 @*/ 1731 PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx) 1732 { 1733 PetscBool useCeed; 1734 1735 PetscFunctionBegin; 1736 PetscCall(DMPlexGetUseCeed(dm, &useCeed)); 1737 PetscCall(DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, boundaryctx)); 1738 if (useCeed) { 1739 #ifdef PETSC_HAVE_LIBCEED 1740 PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualCEED, residualctx)); 1741 #else 1742 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot use CEED traversals without LibCEED. Rerun configure with --download-ceed"); 1743 #endif 1744 } else PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, residualctx)); 1745 PetscCall(DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, jacobianctx)); 1746 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex)); 1747 PetscFunctionReturn(PETSC_SUCCESS); 1748 } 1749 1750 /*@C 1751 DMSNESCheckDiscretization - Check the discretization error of the exact solution 1752 1753 Input Parameters: 1754 + snes - the `SNES` object 1755 . dm - the `DM` 1756 . t - the time 1757 . u - a `DM` vector 1758 - tol - A tolerance for the check, or -1 to print the results instead 1759 1760 Output Parameter: 1761 . error - An array which holds the discretization error in each field, or `NULL` 1762 1763 Level: developer 1764 1765 Note: 1766 The user must call `PetscDSSetExactSolution()` beforehand 1767 1768 Developer Note: 1769 How is this related to `PetscConvEst`? 1770 1771 .seealso: [](ch_snes), `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()` 1772 @*/ 1773 PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[]) 1774 { 1775 PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 1776 void **ectxs; 1777 PetscReal *err; 1778 MPI_Comm comm; 1779 PetscInt Nf, f; 1780 1781 PetscFunctionBegin; 1782 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1783 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1784 PetscValidHeaderSpecific(u, VEC_CLASSID, 4); 1785 if (error) PetscAssertPointer(error, 6); 1786 1787 PetscCall(DMComputeExactSolution(dm, t, u, NULL)); 1788 PetscCall(VecViewFromOptions(u, NULL, "-vec_view")); 1789 1790 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 1791 PetscCall(DMGetNumFields(dm, &Nf)); 1792 PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err)); 1793 { 1794 PetscInt Nds, s; 1795 1796 PetscCall(DMGetNumDS(dm, &Nds)); 1797 for (s = 0; s < Nds; ++s) { 1798 PetscDS ds; 1799 DMLabel label; 1800 IS fieldIS; 1801 const PetscInt *fields; 1802 PetscInt dsNf, f; 1803 1804 PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL)); 1805 PetscCall(PetscDSGetNumFields(ds, &dsNf)); 1806 PetscCall(ISGetIndices(fieldIS, &fields)); 1807 for (f = 0; f < dsNf; ++f) { 1808 const PetscInt field = fields[f]; 1809 PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field])); 1810 } 1811 PetscCall(ISRestoreIndices(fieldIS, &fields)); 1812 } 1813 } 1814 if (Nf > 1) { 1815 PetscCall(DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err)); 1816 if (tol >= 0.0) { 1817 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); 1818 } else if (error) { 1819 for (f = 0; f < Nf; ++f) error[f] = err[f]; 1820 } else { 1821 PetscCall(PetscPrintf(comm, "L_2 Error: [")); 1822 for (f = 0; f < Nf; ++f) { 1823 if (f) PetscCall(PetscPrintf(comm, ", ")); 1824 PetscCall(PetscPrintf(comm, "%g", (double)err[f])); 1825 } 1826 PetscCall(PetscPrintf(comm, "]\n")); 1827 } 1828 } else { 1829 PetscCall(DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0])); 1830 if (tol >= 0.0) { 1831 PetscCheck(err[0] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double)err[0], (double)tol); 1832 } else if (error) { 1833 error[0] = err[0]; 1834 } else { 1835 PetscCall(PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0])); 1836 } 1837 } 1838 PetscCall(PetscFree3(exacts, ectxs, err)); 1839 PetscFunctionReturn(PETSC_SUCCESS); 1840 } 1841 1842 /*@C 1843 DMSNESCheckResidual - Check the residual of the exact solution 1844 1845 Input Parameters: 1846 + snes - the `SNES` object 1847 . dm - the `DM` 1848 . u - a `DM` vector 1849 - tol - A tolerance for the check, or -1 to print the results instead 1850 1851 Output Parameter: 1852 . residual - The residual norm of the exact solution, or `NULL` 1853 1854 Level: developer 1855 1856 .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()` 1857 @*/ 1858 PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual) 1859 { 1860 MPI_Comm comm; 1861 Vec r; 1862 PetscReal res; 1863 1864 PetscFunctionBegin; 1865 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1866 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1867 PetscValidHeaderSpecific(u, VEC_CLASSID, 3); 1868 if (residual) PetscAssertPointer(residual, 5); 1869 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 1870 PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL)); 1871 PetscCall(VecDuplicate(u, &r)); 1872 PetscCall(SNESComputeFunction(snes, u, r)); 1873 PetscCall(VecNorm(r, NORM_2, &res)); 1874 if (tol >= 0.0) { 1875 PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol); 1876 } else if (residual) { 1877 *residual = res; 1878 } else { 1879 PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res)); 1880 PetscCall(VecFilter(r, 1.0e-10)); 1881 PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual")); 1882 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_")); 1883 PetscCall(VecViewFromOptions(r, NULL, "-vec_view")); 1884 } 1885 PetscCall(VecDestroy(&r)); 1886 PetscFunctionReturn(PETSC_SUCCESS); 1887 } 1888 1889 /*@C 1890 DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test 1891 1892 Input Parameters: 1893 + snes - the `SNES` object 1894 . dm - the `DM` 1895 . u - a `DM` vector 1896 - tol - A tolerance for the check, or -1 to print the results instead 1897 1898 Output Parameters: 1899 + isLinear - Flag indicaing that the function looks linear, or `NULL` 1900 - convRate - The rate of convergence of the linear model, or `NULL` 1901 1902 Level: developer 1903 1904 .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()` 1905 @*/ 1906 PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate) 1907 { 1908 MPI_Comm comm; 1909 PetscDS ds; 1910 Mat J, M; 1911 MatNullSpace nullspace; 1912 PetscReal slope, intercept; 1913 PetscBool hasJac, hasPrec, isLin = PETSC_FALSE; 1914 1915 PetscFunctionBegin; 1916 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1917 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1918 PetscValidHeaderSpecific(u, VEC_CLASSID, 3); 1919 if (isLinear) PetscAssertPointer(isLinear, 5); 1920 if (convRate) PetscAssertPointer(convRate, 6); 1921 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 1922 PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL)); 1923 /* Create and view matrices */ 1924 PetscCall(DMCreateMatrix(dm, &J)); 1925 PetscCall(DMGetDS(dm, &ds)); 1926 PetscCall(PetscDSHasJacobian(ds, &hasJac)); 1927 PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec)); 1928 if (hasJac && hasPrec) { 1929 PetscCall(DMCreateMatrix(dm, &M)); 1930 PetscCall(SNESComputeJacobian(snes, u, J, M)); 1931 PetscCall(PetscObjectSetName((PetscObject)M, "Preconditioning Matrix")); 1932 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_")); 1933 PetscCall(MatViewFromOptions(M, NULL, "-mat_view")); 1934 PetscCall(MatDestroy(&M)); 1935 } else { 1936 PetscCall(SNESComputeJacobian(snes, u, J, J)); 1937 } 1938 PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian")); 1939 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_")); 1940 PetscCall(MatViewFromOptions(J, NULL, "-mat_view")); 1941 /* Check nullspace */ 1942 PetscCall(MatGetNullSpace(J, &nullspace)); 1943 if (nullspace) { 1944 PetscBool isNull; 1945 PetscCall(MatNullSpaceTest(nullspace, J, &isNull)); 1946 PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid."); 1947 } 1948 /* Taylor test */ 1949 { 1950 PetscRandom rand; 1951 Vec du, uhat, r, rhat, df; 1952 PetscReal h; 1953 PetscReal *es, *hs, *errors; 1954 PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1; 1955 PetscInt Nv, v; 1956 1957 /* Choose a perturbation direction */ 1958 PetscCall(PetscRandomCreate(comm, &rand)); 1959 PetscCall(VecDuplicate(u, &du)); 1960 PetscCall(VecSetRandom(du, rand)); 1961 PetscCall(PetscRandomDestroy(&rand)); 1962 PetscCall(VecDuplicate(u, &df)); 1963 PetscCall(MatMult(J, du, df)); 1964 /* Evaluate residual at u, F(u), save in vector r */ 1965 PetscCall(VecDuplicate(u, &r)); 1966 PetscCall(SNESComputeFunction(snes, u, r)); 1967 /* Look at the convergence of our Taylor approximation as we approach u */ 1968 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) 1969 ; 1970 PetscCall(PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors)); 1971 PetscCall(VecDuplicate(u, &uhat)); 1972 PetscCall(VecDuplicate(u, &rhat)); 1973 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) { 1974 PetscCall(VecWAXPY(uhat, h, du, u)); 1975 /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */ 1976 PetscCall(SNESComputeFunction(snes, uhat, rhat)); 1977 PetscCall(VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df)); 1978 PetscCall(VecNorm(rhat, NORM_2, &errors[Nv])); 1979 1980 es[Nv] = errors[Nv] == 0 ? -16.0 : PetscLog10Real(errors[Nv]); 1981 hs[Nv] = PetscLog10Real(h); 1982 } 1983 PetscCall(VecDestroy(&uhat)); 1984 PetscCall(VecDestroy(&rhat)); 1985 PetscCall(VecDestroy(&df)); 1986 PetscCall(VecDestroy(&r)); 1987 PetscCall(VecDestroy(&du)); 1988 for (v = 0; v < Nv; ++v) { 1989 if ((tol >= 0) && (errors[v] > tol)) break; 1990 else if (errors[v] > PETSC_SMALL) break; 1991 } 1992 if (v == Nv) isLin = PETSC_TRUE; 1993 PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept)); 1994 PetscCall(PetscFree3(es, hs, errors)); 1995 /* Slope should be about 2 */ 1996 if (tol >= 0) { 1997 PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope); 1998 } else if (isLinear || convRate) { 1999 if (isLinear) *isLinear = isLin; 2000 if (convRate) *convRate = slope; 2001 } else { 2002 if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope)); 2003 else PetscCall(PetscPrintf(comm, "Function appears to be linear\n")); 2004 } 2005 } 2006 PetscCall(MatDestroy(&J)); 2007 PetscFunctionReturn(PETSC_SUCCESS); 2008 } 2009 2010 PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u) 2011 { 2012 PetscFunctionBegin; 2013 PetscCall(DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL)); 2014 PetscCall(DMSNESCheckResidual(snes, dm, u, -1.0, NULL)); 2015 PetscCall(DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL)); 2016 PetscFunctionReturn(PETSC_SUCCESS); 2017 } 2018 2019 /*@C 2020 DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information 2021 2022 Input Parameters: 2023 + snes - the `SNES` object 2024 - u - representative `SNES` vector 2025 2026 Level: developer 2027 2028 Note: 2029 The user must call `PetscDSSetExactSolution()` before this call 2030 2031 .seealso: [](ch_snes), `SNES`, `DM` 2032 @*/ 2033 PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u) 2034 { 2035 DM dm; 2036 Vec sol; 2037 PetscBool check; 2038 2039 PetscFunctionBegin; 2040 PetscCall(PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check)); 2041 if (!check) PetscFunctionReturn(PETSC_SUCCESS); 2042 PetscCall(SNESGetDM(snes, &dm)); 2043 PetscCall(VecDuplicate(u, &sol)); 2044 PetscCall(SNESSetSolution(snes, sol)); 2045 PetscCall(DMSNESCheck_Internal(snes, dm, sol)); 2046 PetscCall(VecDestroy(&sol)); 2047 PetscFunctionReturn(PETSC_SUCCESS); 2048 } 2049