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 /********************* Residual Computation **************************/ 1217 1218 /*@ 1219 DMPlexSNESComputeResidualFEM - Sums the local residual into vector `F` 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 . F - Local output vector 1228 1229 Level: developer 1230 1231 Note: 1232 The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional. 1233 1234 .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()` 1235 @*/ 1236 PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user) 1237 { 1238 DM plex; 1239 IS allcellIS; 1240 PetscInt Nds, s; 1241 1242 PetscFunctionBegin; 1243 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1244 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS)); 1245 PetscCall(DMGetNumDS(dm, &Nds)); 1246 for (s = 0; s < Nds; ++s) { 1247 PetscDS ds; 1248 IS cellIS; 1249 PetscFormKey key; 1250 1251 PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL)); 1252 key.value = 0; 1253 key.field = 0; 1254 key.part = 0; 1255 if (!key.label) { 1256 PetscCall(PetscObjectReference((PetscObject)allcellIS)); 1257 cellIS = allcellIS; 1258 } else { 1259 IS pointIS; 1260 1261 key.value = 1; 1262 PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS)); 1263 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS)); 1264 PetscCall(ISDestroy(&pointIS)); 1265 } 1266 PetscCall(DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user)); 1267 PetscCall(ISDestroy(&cellIS)); 1268 } 1269 PetscCall(ISDestroy(&allcellIS)); 1270 PetscCall(DMDestroy(&plex)); 1271 PetscFunctionReturn(PETSC_SUCCESS); 1272 } 1273 1274 /*@ 1275 DMPlexSNESComputeResidualDS - Sums the local residual into vector `F` from the local input `X` using all pointwise functions with unique keys in the `PetscDS` 1276 1277 Input Parameters: 1278 + dm - The mesh 1279 . X - Local solution 1280 - user - The user context 1281 1282 Output Parameter: 1283 . F - Local output vector 1284 1285 Level: developer 1286 1287 Note: 1288 The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional. 1289 1290 .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()` 1291 @*/ 1292 PetscErrorCode DMPlexSNESComputeResidualDS(DM dm, Vec X, Vec F, void *user) 1293 { 1294 DM plex; 1295 IS allcellIS; 1296 PetscInt Nds, s; 1297 1298 PetscFunctionBegin; 1299 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1300 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS)); 1301 PetscCall(DMGetNumDS(dm, &Nds)); 1302 for (s = 0; s < Nds; ++s) { 1303 PetscDS ds; 1304 DMLabel label; 1305 IS cellIS; 1306 1307 PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL)); 1308 { 1309 PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1}; 1310 PetscWeakForm wf; 1311 PetscInt Nm = 2, m, Nk = 0, k, kp, off = 0; 1312 PetscFormKey *reskeys; 1313 1314 /* Get unique residual keys */ 1315 for (m = 0; m < Nm; ++m) { 1316 PetscInt Nkm; 1317 PetscCall(PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm)); 1318 Nk += Nkm; 1319 } 1320 PetscCall(PetscMalloc1(Nk, &reskeys)); 1321 for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys)); 1322 PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk); 1323 PetscCall(PetscFormKeySort(Nk, reskeys)); 1324 for (k = 0, kp = 1; kp < Nk; ++kp) { 1325 if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) { 1326 ++k; 1327 if (kp != k) reskeys[k] = reskeys[kp]; 1328 } 1329 } 1330 Nk = k; 1331 1332 PetscCall(PetscDSGetWeakForm(ds, &wf)); 1333 for (k = 0; k < Nk; ++k) { 1334 DMLabel label = reskeys[k].label; 1335 PetscInt val = reskeys[k].value; 1336 1337 if (!label) { 1338 PetscCall(PetscObjectReference((PetscObject)allcellIS)); 1339 cellIS = allcellIS; 1340 } else { 1341 IS pointIS; 1342 1343 PetscCall(DMLabelGetStratumIS(label, val, &pointIS)); 1344 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS)); 1345 PetscCall(ISDestroy(&pointIS)); 1346 } 1347 PetscCall(DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user)); 1348 PetscCall(ISDestroy(&cellIS)); 1349 } 1350 PetscCall(PetscFree(reskeys)); 1351 } 1352 } 1353 PetscCall(ISDestroy(&allcellIS)); 1354 PetscCall(DMDestroy(&plex)); 1355 PetscFunctionReturn(PETSC_SUCCESS); 1356 } 1357 1358 #ifdef PETSC_HAVE_LIBCEED 1359 PetscErrorCode DMPlexSNESComputeResidualCEED(DM dm, Vec locX, Vec locF, void *user) 1360 { 1361 Ceed ceed; 1362 DMCeed sd = dm->dmceed; 1363 CeedVector clocX, clocF; 1364 1365 PetscFunctionBegin; 1366 PetscCall(DMGetCeed(dm, &ceed)); 1367 PetscCheck(sd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This DM has no CEED data. Call DMCeedCreate() before computing the residual."); 1368 PetscCall(DMCeedComputeGeometry(dm, sd)); 1369 1370 PetscCall(VecGetCeedVectorRead(locX, ceed, &clocX)); 1371 PetscCall(VecGetCeedVector(locF, ceed, &clocF)); 1372 PetscCallCEED(CeedOperatorApplyAdd(sd->op, clocX, clocF, CEED_REQUEST_IMMEDIATE)); 1373 PetscCall(VecRestoreCeedVectorRead(locX, &clocX)); 1374 PetscCall(VecRestoreCeedVector(locF, &clocF)); 1375 1376 { 1377 DM_Plex *mesh = (DM_Plex *)dm->data; 1378 1379 if (mesh->printFEM) { 1380 PetscSection section; 1381 Vec locFbc; 1382 PetscInt pStart, pEnd, p, maxDof; 1383 PetscScalar *zeroes; 1384 1385 PetscCall(DMGetLocalSection(dm, §ion)); 1386 PetscCall(VecDuplicate(locF, &locFbc)); 1387 PetscCall(VecCopy(locF, locFbc)); 1388 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 1389 PetscCall(PetscSectionGetMaxDof(section, &maxDof)); 1390 PetscCall(PetscCalloc1(maxDof, &zeroes)); 1391 for (p = pStart; p < pEnd; ++p) PetscCall(VecSetValuesSection(locFbc, section, p, zeroes, INSERT_BC_VALUES)); 1392 PetscCall(PetscFree(zeroes)); 1393 PetscCall(DMPrintLocalVec(dm, "Residual", mesh->printTol, locFbc)); 1394 PetscCall(VecDestroy(&locFbc)); 1395 } 1396 } 1397 PetscFunctionReturn(PETSC_SUCCESS); 1398 } 1399 #endif 1400 1401 /*@ 1402 DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input `X` 1403 1404 Input Parameters: 1405 + dm - The mesh 1406 - user - The user context 1407 1408 Output Parameter: 1409 . X - Local solution 1410 1411 Level: developer 1412 1413 .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()` 1414 @*/ 1415 PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user) 1416 { 1417 DM plex; 1418 1419 PetscFunctionBegin; 1420 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1421 PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL)); 1422 PetscCall(DMDestroy(&plex)); 1423 PetscFunctionReturn(PETSC_SUCCESS); 1424 } 1425 1426 /*@ 1427 DMSNESComputeJacobianAction - Compute the action of the Jacobian J(X) on Y 1428 1429 Input Parameters: 1430 + dm - The `DM` 1431 . X - Local solution vector 1432 . Y - Local input vector 1433 - user - The user context 1434 1435 Output Parameter: 1436 . F - local output vector 1437 1438 Level: developer 1439 1440 Notes: 1441 Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly. 1442 1443 .seealso: [](ch_snes), `DM`, ``DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()` 1444 @*/ 1445 PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, void *user) 1446 { 1447 DM plex; 1448 IS allcellIS; 1449 PetscInt Nds, s; 1450 1451 PetscFunctionBegin; 1452 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1453 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS)); 1454 PetscCall(DMGetNumDS(dm, &Nds)); 1455 for (s = 0; s < Nds; ++s) { 1456 PetscDS ds; 1457 DMLabel label; 1458 IS cellIS; 1459 1460 PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL)); 1461 { 1462 PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3}; 1463 PetscWeakForm wf; 1464 PetscInt Nm = 4, m, Nk = 0, k, kp, off = 0; 1465 PetscFormKey *jackeys; 1466 1467 /* Get unique Jacobian keys */ 1468 for (m = 0; m < Nm; ++m) { 1469 PetscInt Nkm; 1470 PetscCall(PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm)); 1471 Nk += Nkm; 1472 } 1473 PetscCall(PetscMalloc1(Nk, &jackeys)); 1474 for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys)); 1475 PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk); 1476 PetscCall(PetscFormKeySort(Nk, jackeys)); 1477 for (k = 0, kp = 1; kp < Nk; ++kp) { 1478 if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) { 1479 ++k; 1480 if (kp != k) jackeys[k] = jackeys[kp]; 1481 } 1482 } 1483 Nk = k; 1484 1485 PetscCall(PetscDSGetWeakForm(ds, &wf)); 1486 for (k = 0; k < Nk; ++k) { 1487 DMLabel label = jackeys[k].label; 1488 PetscInt val = jackeys[k].value; 1489 1490 if (!label) { 1491 PetscCall(PetscObjectReference((PetscObject)allcellIS)); 1492 cellIS = allcellIS; 1493 } else { 1494 IS pointIS; 1495 1496 PetscCall(DMLabelGetStratumIS(label, val, &pointIS)); 1497 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS)); 1498 PetscCall(ISDestroy(&pointIS)); 1499 } 1500 PetscCall(DMPlexComputeJacobian_Action_Internal(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, user)); 1501 PetscCall(ISDestroy(&cellIS)); 1502 } 1503 PetscCall(PetscFree(jackeys)); 1504 } 1505 } 1506 PetscCall(ISDestroy(&allcellIS)); 1507 PetscCall(DMDestroy(&plex)); 1508 PetscFunctionReturn(PETSC_SUCCESS); 1509 } 1510 1511 /*@ 1512 DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix `Jac` at the local solution `X` using pointwise functions specified by the user. 1513 1514 Input Parameters: 1515 + dm - The `DM` 1516 . X - Local input vector 1517 - user - The user context 1518 1519 Output Parameters: 1520 + Jac - Jacobian matrix 1521 - JacP - approximate Jacobian from which the preconditioner will be built, often `Jac` 1522 1523 Level: developer 1524 1525 Note: 1526 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1527 like a GPU, or vectorize on a multicore machine. 1528 1529 .seealso: [](ch_snes), `DMPLEX`, `Mat` 1530 @*/ 1531 PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, void *user) 1532 { 1533 DM plex; 1534 IS allcellIS; 1535 PetscBool hasJac, hasPrec; 1536 PetscInt Nds, s; 1537 1538 PetscFunctionBegin; 1539 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE)); 1540 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS)); 1541 PetscCall(DMGetNumDS(dm, &Nds)); 1542 for (s = 0; s < Nds; ++s) { 1543 PetscDS ds; 1544 IS cellIS; 1545 PetscFormKey key; 1546 1547 PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL)); 1548 key.value = 0; 1549 key.field = 0; 1550 key.part = 0; 1551 if (!key.label) { 1552 PetscCall(PetscObjectReference((PetscObject)allcellIS)); 1553 cellIS = allcellIS; 1554 } else { 1555 IS pointIS; 1556 1557 key.value = 1; 1558 PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS)); 1559 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS)); 1560 PetscCall(ISDestroy(&pointIS)); 1561 } 1562 if (!s) { 1563 PetscCall(PetscDSHasJacobian(ds, &hasJac)); 1564 PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec)); 1565 if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac)); 1566 PetscCall(MatZeroEntries(JacP)); 1567 } 1568 PetscCall(DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user)); 1569 PetscCall(ISDestroy(&cellIS)); 1570 } 1571 PetscCall(ISDestroy(&allcellIS)); 1572 PetscCall(DMDestroy(&plex)); 1573 PetscFunctionReturn(PETSC_SUCCESS); 1574 } 1575 1576 struct _DMSNESJacobianMFCtx { 1577 DM dm; 1578 Vec X; 1579 void *ctx; 1580 }; 1581 1582 static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A) 1583 { 1584 struct _DMSNESJacobianMFCtx *ctx; 1585 1586 PetscFunctionBegin; 1587 PetscCall(MatShellGetContext(A, &ctx)); 1588 PetscCall(MatShellSetContext(A, NULL)); 1589 PetscCall(DMDestroy(&ctx->dm)); 1590 PetscCall(VecDestroy(&ctx->X)); 1591 PetscCall(PetscFree(ctx)); 1592 PetscFunctionReturn(PETSC_SUCCESS); 1593 } 1594 1595 static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z) 1596 { 1597 struct _DMSNESJacobianMFCtx *ctx; 1598 1599 PetscFunctionBegin; 1600 PetscCall(MatShellGetContext(A, &ctx)); 1601 PetscCall(DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx)); 1602 PetscFunctionReturn(PETSC_SUCCESS); 1603 } 1604 1605 /*@ 1606 DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free 1607 1608 Collective 1609 1610 Input Parameters: 1611 + dm - The `DM` 1612 . X - The evaluation point for the Jacobian 1613 - user - A user context, or `NULL` 1614 1615 Output Parameter: 1616 . J - The `Mat` 1617 1618 Level: advanced 1619 1620 Note: 1621 Vec `X` is kept in `J`, so updating `X` then updates the evaluation point. 1622 1623 .seealso: [](ch_snes), `DM`, `SNES`, `DMSNESComputeJacobianAction()` 1624 @*/ 1625 PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, void *user, Mat *J) 1626 { 1627 struct _DMSNESJacobianMFCtx *ctx; 1628 PetscInt n, N; 1629 1630 PetscFunctionBegin; 1631 PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J)); 1632 PetscCall(MatSetType(*J, MATSHELL)); 1633 PetscCall(VecGetLocalSize(X, &n)); 1634 PetscCall(VecGetSize(X, &N)); 1635 PetscCall(MatSetSizes(*J, n, n, N, N)); 1636 PetscCall(PetscObjectReference((PetscObject)dm)); 1637 PetscCall(PetscObjectReference((PetscObject)X)); 1638 PetscCall(PetscMalloc1(1, &ctx)); 1639 ctx->dm = dm; 1640 ctx->X = X; 1641 ctx->ctx = user; 1642 PetscCall(MatShellSetContext(*J, ctx)); 1643 PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))DMSNESJacobianMF_Destroy_Private)); 1644 PetscCall(MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))DMSNESJacobianMF_Mult_Private)); 1645 PetscFunctionReturn(PETSC_SUCCESS); 1646 } 1647 1648 /* 1649 MatComputeNeumannOverlap - Computes an unassembled (Neumann) local overlapping Mat in nonlinear context. 1650 1651 Input Parameters: 1652 + X - `SNES` linearization point 1653 . ovl - index set of overlapping subdomains 1654 1655 Output Parameter: 1656 . J - unassembled (Neumann) local matrix 1657 1658 Level: intermediate 1659 1660 .seealso: [](ch_snes), `DMCreateNeumannOverlap()`, `MATIS`, `PCHPDDMSetAuxiliaryMat()` 1661 */ 1662 static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx) 1663 { 1664 SNES snes; 1665 Mat pJ; 1666 DM ovldm, origdm; 1667 DMSNES sdm; 1668 PetscErrorCode (*bfun)(DM, Vec, void *); 1669 PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *); 1670 void *bctx, *jctx; 1671 1672 PetscFunctionBegin; 1673 PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ)); 1674 PetscCheck(pJ, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing overlapping Mat"); 1675 PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm)); 1676 PetscCheck(origdm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing original DM"); 1677 PetscCall(MatGetDM(pJ, &ovldm)); 1678 PetscCall(DMSNESGetBoundaryLocal(origdm, &bfun, &bctx)); 1679 PetscCall(DMSNESSetBoundaryLocal(ovldm, bfun, bctx)); 1680 PetscCall(DMSNESGetJacobianLocal(origdm, &jfun, &jctx)); 1681 PetscCall(DMSNESSetJacobianLocal(ovldm, jfun, jctx)); 1682 PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes)); 1683 if (!snes) { 1684 PetscCall(SNESCreate(PetscObjectComm((PetscObject)ovl), &snes)); 1685 PetscCall(SNESSetDM(snes, ovldm)); 1686 PetscCall(PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes)); 1687 PetscCall(PetscObjectDereference((PetscObject)snes)); 1688 } 1689 PetscCall(DMGetDMSNES(ovldm, &sdm)); 1690 PetscCall(VecLockReadPush(X)); 1691 { 1692 void *ctx; 1693 PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *); 1694 PetscCall(DMSNESGetJacobian(ovldm, &J, &ctx)); 1695 PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx)); 1696 } 1697 PetscCall(VecLockReadPop(X)); 1698 /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */ 1699 { 1700 Mat locpJ; 1701 1702 PetscCall(MatISGetLocalMat(pJ, &locpJ)); 1703 PetscCall(MatCopy(locpJ, J, SAME_NONZERO_PATTERN)); 1704 } 1705 PetscFunctionReturn(PETSC_SUCCESS); 1706 } 1707 1708 /*@ 1709 DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, residual, and Jacobian. 1710 1711 Input Parameters: 1712 + dm - The `DM` object 1713 . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see `PetscDSAddBoundary()`) 1714 . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see `PetscDSSetResidual()`) 1715 - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see `PetscDSSetJacobian()`) 1716 1717 Level: developer 1718 1719 .seealso: [](ch_snes), `DMPLEX`, `SNES` 1720 @*/ 1721 PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx) 1722 { 1723 PetscBool useCeed; 1724 1725 PetscFunctionBegin; 1726 PetscCall(DMPlexGetUseCeed(dm, &useCeed)); 1727 PetscCall(DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, boundaryctx)); 1728 if (useCeed) { 1729 #ifdef PETSC_HAVE_LIBCEED 1730 PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualCEED, residualctx)); 1731 #else 1732 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot use CEED traversals without LibCEED. Rerun configure with --download-ceed"); 1733 #endif 1734 } else PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, residualctx)); 1735 PetscCall(DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, jacobianctx)); 1736 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex)); 1737 PetscFunctionReturn(PETSC_SUCCESS); 1738 } 1739 1740 /*@C 1741 DMSNESCheckDiscretization - Check the discretization error of the exact solution 1742 1743 Input Parameters: 1744 + snes - the `SNES` object 1745 . dm - the `DM` 1746 . t - the time 1747 . u - a `DM` vector 1748 - tol - A tolerance for the check, or -1 to print the results instead 1749 1750 Output Parameter: 1751 . error - An array which holds the discretization error in each field, or `NULL` 1752 1753 Level: developer 1754 1755 Note: 1756 The user must call `PetscDSSetExactSolution()` beforehand 1757 1758 Developer Note: 1759 How is this related to `PetscConvEst`? 1760 1761 .seealso: [](ch_snes), `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()` 1762 @*/ 1763 PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[]) 1764 { 1765 PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 1766 void **ectxs; 1767 PetscReal *err; 1768 MPI_Comm comm; 1769 PetscInt Nf, f; 1770 1771 PetscFunctionBegin; 1772 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1773 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1774 PetscValidHeaderSpecific(u, VEC_CLASSID, 4); 1775 if (error) PetscAssertPointer(error, 6); 1776 1777 PetscCall(DMComputeExactSolution(dm, t, u, NULL)); 1778 PetscCall(VecViewFromOptions(u, NULL, "-vec_view")); 1779 1780 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 1781 PetscCall(DMGetNumFields(dm, &Nf)); 1782 PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err)); 1783 { 1784 PetscInt Nds, s; 1785 1786 PetscCall(DMGetNumDS(dm, &Nds)); 1787 for (s = 0; s < Nds; ++s) { 1788 PetscDS ds; 1789 DMLabel label; 1790 IS fieldIS; 1791 const PetscInt *fields; 1792 PetscInt dsNf, f; 1793 1794 PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL)); 1795 PetscCall(PetscDSGetNumFields(ds, &dsNf)); 1796 PetscCall(ISGetIndices(fieldIS, &fields)); 1797 for (f = 0; f < dsNf; ++f) { 1798 const PetscInt field = fields[f]; 1799 PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field])); 1800 } 1801 PetscCall(ISRestoreIndices(fieldIS, &fields)); 1802 } 1803 } 1804 if (Nf > 1) { 1805 PetscCall(DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err)); 1806 if (tol >= 0.0) { 1807 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); 1808 } else if (error) { 1809 for (f = 0; f < Nf; ++f) error[f] = err[f]; 1810 } else { 1811 PetscCall(PetscPrintf(comm, "L_2 Error: [")); 1812 for (f = 0; f < Nf; ++f) { 1813 if (f) PetscCall(PetscPrintf(comm, ", ")); 1814 PetscCall(PetscPrintf(comm, "%g", (double)err[f])); 1815 } 1816 PetscCall(PetscPrintf(comm, "]\n")); 1817 } 1818 } else { 1819 PetscCall(DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0])); 1820 if (tol >= 0.0) { 1821 PetscCheck(err[0] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double)err[0], (double)tol); 1822 } else if (error) { 1823 error[0] = err[0]; 1824 } else { 1825 PetscCall(PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0])); 1826 } 1827 } 1828 PetscCall(PetscFree3(exacts, ectxs, err)); 1829 PetscFunctionReturn(PETSC_SUCCESS); 1830 } 1831 1832 /*@C 1833 DMSNESCheckResidual - Check the residual of the exact solution 1834 1835 Input Parameters: 1836 + snes - the `SNES` object 1837 . dm - the `DM` 1838 . u - a `DM` vector 1839 - tol - A tolerance for the check, or -1 to print the results instead 1840 1841 Output Parameter: 1842 . residual - The residual norm of the exact solution, or `NULL` 1843 1844 Level: developer 1845 1846 .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()` 1847 @*/ 1848 PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual) 1849 { 1850 MPI_Comm comm; 1851 Vec r; 1852 PetscReal res; 1853 1854 PetscFunctionBegin; 1855 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1856 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1857 PetscValidHeaderSpecific(u, VEC_CLASSID, 3); 1858 if (residual) PetscAssertPointer(residual, 5); 1859 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 1860 PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL)); 1861 PetscCall(VecDuplicate(u, &r)); 1862 PetscCall(SNESComputeFunction(snes, u, r)); 1863 PetscCall(VecNorm(r, NORM_2, &res)); 1864 if (tol >= 0.0) { 1865 PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol); 1866 } else if (residual) { 1867 *residual = res; 1868 } else { 1869 PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res)); 1870 PetscCall(VecFilter(r, 1.0e-10)); 1871 PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual")); 1872 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_")); 1873 PetscCall(VecViewFromOptions(r, NULL, "-vec_view")); 1874 } 1875 PetscCall(VecDestroy(&r)); 1876 PetscFunctionReturn(PETSC_SUCCESS); 1877 } 1878 1879 /*@C 1880 DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test 1881 1882 Input Parameters: 1883 + snes - the `SNES` object 1884 . dm - the `DM` 1885 . u - a `DM` vector 1886 - tol - A tolerance for the check, or -1 to print the results instead 1887 1888 Output Parameters: 1889 + isLinear - Flag indicaing that the function looks linear, or `NULL` 1890 - convRate - The rate of convergence of the linear model, or `NULL` 1891 1892 Level: developer 1893 1894 .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()` 1895 @*/ 1896 PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate) 1897 { 1898 MPI_Comm comm; 1899 PetscDS ds; 1900 Mat J, M; 1901 MatNullSpace nullspace; 1902 PetscReal slope, intercept; 1903 PetscBool hasJac, hasPrec, isLin = PETSC_FALSE; 1904 1905 PetscFunctionBegin; 1906 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1907 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1908 PetscValidHeaderSpecific(u, VEC_CLASSID, 3); 1909 if (isLinear) PetscAssertPointer(isLinear, 5); 1910 if (convRate) PetscAssertPointer(convRate, 6); 1911 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 1912 PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL)); 1913 /* Create and view matrices */ 1914 PetscCall(DMCreateMatrix(dm, &J)); 1915 PetscCall(DMGetDS(dm, &ds)); 1916 PetscCall(PetscDSHasJacobian(ds, &hasJac)); 1917 PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec)); 1918 if (hasJac && hasPrec) { 1919 PetscCall(DMCreateMatrix(dm, &M)); 1920 PetscCall(SNESComputeJacobian(snes, u, J, M)); 1921 PetscCall(PetscObjectSetName((PetscObject)M, "Preconditioning Matrix")); 1922 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_")); 1923 PetscCall(MatViewFromOptions(M, NULL, "-mat_view")); 1924 PetscCall(MatDestroy(&M)); 1925 } else { 1926 PetscCall(SNESComputeJacobian(snes, u, J, J)); 1927 } 1928 PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian")); 1929 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_")); 1930 PetscCall(MatViewFromOptions(J, NULL, "-mat_view")); 1931 /* Check nullspace */ 1932 PetscCall(MatGetNullSpace(J, &nullspace)); 1933 if (nullspace) { 1934 PetscBool isNull; 1935 PetscCall(MatNullSpaceTest(nullspace, J, &isNull)); 1936 PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid."); 1937 } 1938 /* Taylor test */ 1939 { 1940 PetscRandom rand; 1941 Vec du, uhat, r, rhat, df; 1942 PetscReal h; 1943 PetscReal *es, *hs, *errors; 1944 PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1; 1945 PetscInt Nv, v; 1946 1947 /* Choose a perturbation direction */ 1948 PetscCall(PetscRandomCreate(comm, &rand)); 1949 PetscCall(VecDuplicate(u, &du)); 1950 PetscCall(VecSetRandom(du, rand)); 1951 PetscCall(PetscRandomDestroy(&rand)); 1952 PetscCall(VecDuplicate(u, &df)); 1953 PetscCall(MatMult(J, du, df)); 1954 /* Evaluate residual at u, F(u), save in vector r */ 1955 PetscCall(VecDuplicate(u, &r)); 1956 PetscCall(SNESComputeFunction(snes, u, r)); 1957 /* Look at the convergence of our Taylor approximation as we approach u */ 1958 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) 1959 ; 1960 PetscCall(PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors)); 1961 PetscCall(VecDuplicate(u, &uhat)); 1962 PetscCall(VecDuplicate(u, &rhat)); 1963 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) { 1964 PetscCall(VecWAXPY(uhat, h, du, u)); 1965 /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */ 1966 PetscCall(SNESComputeFunction(snes, uhat, rhat)); 1967 PetscCall(VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df)); 1968 PetscCall(VecNorm(rhat, NORM_2, &errors[Nv])); 1969 1970 es[Nv] = errors[Nv] == 0 ? -16.0 : PetscLog10Real(errors[Nv]); 1971 hs[Nv] = PetscLog10Real(h); 1972 } 1973 PetscCall(VecDestroy(&uhat)); 1974 PetscCall(VecDestroy(&rhat)); 1975 PetscCall(VecDestroy(&df)); 1976 PetscCall(VecDestroy(&r)); 1977 PetscCall(VecDestroy(&du)); 1978 for (v = 0; v < Nv; ++v) { 1979 if ((tol >= 0) && (errors[v] > tol)) break; 1980 else if (errors[v] > PETSC_SMALL) break; 1981 } 1982 if (v == Nv) isLin = PETSC_TRUE; 1983 PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept)); 1984 PetscCall(PetscFree3(es, hs, errors)); 1985 /* Slope should be about 2 */ 1986 if (tol >= 0) { 1987 PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope); 1988 } else if (isLinear || convRate) { 1989 if (isLinear) *isLinear = isLin; 1990 if (convRate) *convRate = slope; 1991 } else { 1992 if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope)); 1993 else PetscCall(PetscPrintf(comm, "Function appears to be linear\n")); 1994 } 1995 } 1996 PetscCall(MatDestroy(&J)); 1997 PetscFunctionReturn(PETSC_SUCCESS); 1998 } 1999 2000 PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u) 2001 { 2002 PetscFunctionBegin; 2003 PetscCall(DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL)); 2004 PetscCall(DMSNESCheckResidual(snes, dm, u, -1.0, NULL)); 2005 PetscCall(DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL)); 2006 PetscFunctionReturn(PETSC_SUCCESS); 2007 } 2008 2009 /*@C 2010 DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information 2011 2012 Input Parameters: 2013 + snes - the `SNES` object 2014 - u - representative `SNES` vector 2015 2016 Level: developer 2017 2018 Note: 2019 The user must call `PetscDSSetExactSolution()` before this call 2020 2021 .seealso: [](ch_snes), `SNES`, `DM` 2022 @*/ 2023 PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u) 2024 { 2025 DM dm; 2026 Vec sol; 2027 PetscBool check; 2028 2029 PetscFunctionBegin; 2030 PetscCall(PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check)); 2031 if (!check) PetscFunctionReturn(PETSC_SUCCESS); 2032 PetscCall(SNESGetDM(snes, &dm)); 2033 PetscCall(VecDuplicate(u, &sol)); 2034 PetscCall(SNESSetSolution(snes, sol)); 2035 PetscCall(DMSNESCheck_Internal(snes, dm, sol)); 2036 PetscCall(VecDestroy(&sol)); 2037 PetscFunctionReturn(PETSC_SUCCESS); 2038 } 2039