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