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