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