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