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