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 - Form the local residual 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 Level: developer 1204 1205 .seealso: DMPlexComputeJacobianAction() 1206 @*/ 1207 PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user) 1208 { 1209 DM plex; 1210 IS allcellIS; 1211 PetscInt Nds, s; 1212 PetscErrorCode ierr; 1213 1214 PetscFunctionBegin; 1215 ierr = DMSNESConvertPlex(dm, &plex, PETSC_TRUE);CHKERRQ(ierr); 1216 ierr = DMPlexGetAllCells_Internal(plex, &allcellIS);CHKERRQ(ierr); 1217 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr); 1218 for (s = 0; s < Nds; ++s) { 1219 PetscDS ds; 1220 IS cellIS; 1221 PetscHashFormKey key; 1222 1223 ierr = DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);CHKERRQ(ierr); 1224 key.value = 0; 1225 key.field = 0; 1226 if (!key.label) { 1227 ierr = PetscObjectReference((PetscObject) allcellIS);CHKERRQ(ierr); 1228 cellIS = allcellIS; 1229 } else { 1230 IS pointIS; 1231 1232 key.value = 1; 1233 ierr = DMLabelGetStratumIS(key.label, key.value, &pointIS);CHKERRQ(ierr); 1234 ierr = ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);CHKERRQ(ierr); 1235 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1236 } 1237 ierr = DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user);CHKERRQ(ierr); 1238 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 1239 } 1240 ierr = ISDestroy(&allcellIS);CHKERRQ(ierr); 1241 ierr = DMDestroy(&plex);CHKERRQ(ierr); 1242 PetscFunctionReturn(0); 1243 } 1244 1245 PetscErrorCode DMSNESComputeResidual(DM dm, Vec X, Vec F, void *user) 1246 { 1247 DM plex; 1248 IS allcellIS; 1249 PetscInt Nds, s; 1250 PetscErrorCode ierr; 1251 1252 PetscFunctionBegin; 1253 ierr = DMSNESConvertPlex(dm, &plex, PETSC_TRUE);CHKERRQ(ierr); 1254 ierr = DMPlexGetAllCells_Internal(plex, &allcellIS);CHKERRQ(ierr); 1255 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr); 1256 for (s = 0; s < Nds; ++s) { 1257 PetscDS ds; 1258 DMLabel label; 1259 IS cellIS; 1260 1261 ierr = DMGetRegionNumDS(dm, s, &label, NULL, &ds);CHKERRQ(ierr); 1262 { 1263 PetscHMapForm resmap[2] = {ds->wf->f0, ds->wf->f1}; 1264 PetscWeakForm wf; 1265 PetscInt Nm = 2, m, Nk = 0, k, kp, off = 0; 1266 PetscHashFormKey *reskeys; 1267 1268 /* Get unique residual keys */ 1269 for (m = 0; m < Nm; ++m) { 1270 PetscInt Nkm; 1271 ierr = PetscHMapFormGetSize(resmap[m], &Nkm);CHKERRQ(ierr); 1272 Nk += Nkm; 1273 } 1274 ierr = PetscMalloc1(Nk, &reskeys);CHKERRQ(ierr); 1275 for (m = 0; m < Nm; ++m) { 1276 ierr = PetscHMapFormGetKeys(resmap[m], &off, reskeys);CHKERRQ(ierr); 1277 } 1278 if (off != Nk) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %D should be %D", off, Nk); 1279 ierr = PetscHashFormKeySort(Nk, reskeys);CHKERRQ(ierr); 1280 for (k = 0, kp = 1; kp < Nk; ++kp) { 1281 if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) { 1282 ++k; 1283 if (kp != k) reskeys[k] = reskeys[kp]; 1284 } 1285 } 1286 Nk = k; 1287 1288 ierr = PetscDSGetWeakForm(ds, &wf);CHKERRQ(ierr); 1289 for (k = 0; k < Nk; ++k) { 1290 DMLabel label = reskeys[k].label; 1291 PetscInt val = reskeys[k].value; 1292 1293 if (!label) { 1294 ierr = PetscObjectReference((PetscObject) allcellIS);CHKERRQ(ierr); 1295 cellIS = allcellIS; 1296 } else { 1297 IS pointIS; 1298 1299 ierr = DMLabelGetStratumIS(label, val, &pointIS);CHKERRQ(ierr); 1300 ierr = ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);CHKERRQ(ierr); 1301 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1302 } 1303 ierr = DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user);CHKERRQ(ierr); 1304 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 1305 } 1306 ierr = PetscFree(reskeys);CHKERRQ(ierr); 1307 } 1308 } 1309 ierr = ISDestroy(&allcellIS);CHKERRQ(ierr); 1310 ierr = DMDestroy(&plex);CHKERRQ(ierr); 1311 PetscFunctionReturn(0); 1312 } 1313 1314 /*@ 1315 DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input X 1316 1317 Input Parameters: 1318 + dm - The mesh 1319 - user - The user context 1320 1321 Output Parameter: 1322 . X - Local solution 1323 1324 Level: developer 1325 1326 .seealso: DMPlexComputeJacobianAction() 1327 @*/ 1328 PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user) 1329 { 1330 DM plex; 1331 PetscErrorCode ierr; 1332 1333 PetscFunctionBegin; 1334 ierr = DMSNESConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr); 1335 ierr = DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL);CHKERRQ(ierr); 1336 ierr = DMDestroy(&plex);CHKERRQ(ierr); 1337 PetscFunctionReturn(0); 1338 } 1339 1340 /*@ 1341 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. 1342 1343 Input Parameters: 1344 + dm - The mesh 1345 . cellIS - 1346 . t - The time 1347 . X_tShift - The multiplier for the Jacobian with repsect to X_t 1348 . X - Local solution vector 1349 . X_t - Time-derivative of the local solution vector 1350 . Y - Local input vector 1351 - user - The user context 1352 1353 Output Parameter: 1354 . Z - Local output vector 1355 1356 Note: 1357 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1358 like a GPU, or vectorize on a multicore machine. 1359 1360 Level: developer 1361 1362 .seealso: FormFunctionLocal() 1363 @*/ 1364 PetscErrorCode DMPlexComputeJacobianAction(DM dm, IS cellIS, PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Vec Y, Vec Z, void *user) 1365 { 1366 DM_Plex *mesh = (DM_Plex *) dm->data; 1367 const char *name = "Jacobian"; 1368 DM dmAux, plex, plexAux = NULL; 1369 DMEnclosureType encAux; 1370 Vec A; 1371 PetscDS prob, probAux = NULL; 1372 PetscQuadrature quad; 1373 PetscSection section, globalSection, sectionAux; 1374 PetscScalar *elemMat, *elemMatD, *u, *u_t, *a = NULL, *y, *z; 1375 PetscInt Nf, fieldI, fieldJ; 1376 PetscInt totDim, totDimAux = 0; 1377 const PetscInt *cells; 1378 PetscInt cStart, cEnd, numCells, c; 1379 PetscBool hasDyn; 1380 DMField coordField; 1381 PetscHashFormKey key; 1382 PetscErrorCode ierr; 1383 1384 PetscFunctionBegin; 1385 ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 1386 ierr = DMSNESConvertPlex(dm, &plex, PETSC_TRUE);CHKERRQ(ierr); 1387 if (!cellIS) { 1388 PetscInt depth; 1389 1390 ierr = DMPlexGetDepth(plex, &depth);CHKERRQ(ierr); 1391 ierr = DMGetStratumIS(plex, "dim", depth, &cellIS);CHKERRQ(ierr); 1392 if (!cellIS) {ierr = DMGetStratumIS(plex, "depth", depth, &cellIS);CHKERRQ(ierr);} 1393 } else { 1394 ierr = PetscObjectReference((PetscObject) cellIS);CHKERRQ(ierr); 1395 } 1396 key.label = NULL; 1397 key.value = 0; 1398 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 1399 ierr = DMGetGlobalSection(dm, &globalSection);CHKERRQ(ierr); 1400 ierr = ISGetLocalSize(cellIS, &numCells);CHKERRQ(ierr); 1401 ierr = ISGetPointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 1402 ierr = DMGetCellDS(dm, cells ? cells[cStart] : cStart, &prob);CHKERRQ(ierr); 1403 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1404 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1405 ierr = PetscDSHasDynamicJacobian(prob, &hasDyn);CHKERRQ(ierr); 1406 hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE; 1407 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1408 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1409 if (dmAux) { 1410 ierr = DMGetEnclosureRelation(dmAux, dm, &encAux);CHKERRQ(ierr); 1411 ierr = DMConvert(dmAux, DMPLEX, &plexAux);CHKERRQ(ierr); 1412 ierr = DMGetLocalSection(plexAux, §ionAux);CHKERRQ(ierr); 1413 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1414 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1415 } 1416 ierr = VecSet(Z, 0.0);CHKERRQ(ierr); 1417 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); 1418 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1419 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 1420 for (c = cStart; c < cEnd; ++c) { 1421 const PetscInt cell = cells ? cells[c] : c; 1422 const PetscInt cind = c - cStart; 1423 PetscScalar *x = NULL, *x_t = NULL; 1424 PetscInt i; 1425 1426 ierr = DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);CHKERRQ(ierr); 1427 for (i = 0; i < totDim; ++i) u[cind*totDim+i] = x[i]; 1428 ierr = DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);CHKERRQ(ierr); 1429 if (X_t) { 1430 ierr = DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);CHKERRQ(ierr); 1431 for (i = 0; i < totDim; ++i) u_t[cind*totDim+i] = x_t[i]; 1432 ierr = DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);CHKERRQ(ierr); 1433 } 1434 if (dmAux) { 1435 PetscInt subcell; 1436 ierr = DMGetEnclosurePoint(dmAux, dm, encAux, cell, &subcell);CHKERRQ(ierr); 1437 ierr = DMPlexVecGetClosure(plexAux, sectionAux, A, subcell, NULL, &x);CHKERRQ(ierr); 1438 for (i = 0; i < totDimAux; ++i) a[cind*totDimAux+i] = x[i]; 1439 ierr = DMPlexVecRestoreClosure(plexAux, sectionAux, A, subcell, NULL, &x);CHKERRQ(ierr); 1440 } 1441 ierr = DMPlexVecGetClosure(dm, section, Y, cell, NULL, &x);CHKERRQ(ierr); 1442 for (i = 0; i < totDim; ++i) y[cind*totDim+i] = x[i]; 1443 ierr = DMPlexVecRestoreClosure(dm, section, Y, cell, NULL, &x);CHKERRQ(ierr); 1444 } 1445 ierr = PetscArrayzero(elemMat, numCells*totDim*totDim);CHKERRQ(ierr); 1446 if (hasDyn) {ierr = PetscArrayzero(elemMatD, numCells*totDim*totDim);CHKERRQ(ierr);} 1447 for (fieldI = 0; fieldI < Nf; ++fieldI) { 1448 PetscFE fe; 1449 PetscInt Nb; 1450 /* Conforming batches */ 1451 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1452 /* Remainder */ 1453 PetscInt Nr, offset, Nq; 1454 PetscQuadrature qGeom = NULL; 1455 PetscInt maxDegree; 1456 PetscFEGeom *cgeomFEM, *chunkGeom = NULL, *remGeom = NULL; 1457 1458 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); 1459 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1460 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1461 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1462 ierr = DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);CHKERRQ(ierr); 1463 if (maxDegree <= 1) {ierr = DMFieldCreateDefaultQuadrature(coordField,cellIS,&qGeom);CHKERRQ(ierr);} 1464 if (!qGeom) { 1465 ierr = PetscFEGetQuadrature(fe,&qGeom);CHKERRQ(ierr); 1466 ierr = PetscObjectReference((PetscObject)qGeom);CHKERRQ(ierr); 1467 } 1468 ierr = PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1469 ierr = DMSNESGetFEGeom(coordField,cellIS,qGeom,PETSC_FALSE,&cgeomFEM);CHKERRQ(ierr); 1470 blockSize = Nb; 1471 batchSize = numBlocks * blockSize; 1472 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1473 numChunks = numCells / (numBatches*batchSize); 1474 Ne = numChunks*numBatches*batchSize; 1475 Nr = numCells % (numBatches*batchSize); 1476 offset = numCells - Nr; 1477 ierr = PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);CHKERRQ(ierr); 1478 ierr = PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&remGeom);CHKERRQ(ierr); 1479 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) { 1480 key.field = fieldI*Nf + fieldJ; 1481 ierr = PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, key, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMat);CHKERRQ(ierr); 1482 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); 1483 if (hasDyn) { 1484 ierr = PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, key, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMatD);CHKERRQ(ierr); 1485 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); 1486 } 1487 } 1488 ierr = PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&remGeom);CHKERRQ(ierr); 1489 ierr = PetscFEGeomRestoreChunk(cgeomFEM,0,offset,&chunkGeom);CHKERRQ(ierr); 1490 ierr = DMSNESRestoreFEGeom(coordField,cellIS,qGeom,PETSC_FALSE,&cgeomFEM);CHKERRQ(ierr); 1491 ierr = PetscQuadratureDestroy(&qGeom);CHKERRQ(ierr); 1492 } 1493 if (hasDyn) { 1494 for (c = 0; c < numCells*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c]; 1495 } 1496 for (c = cStart; c < cEnd; ++c) { 1497 const PetscInt cell = cells ? cells[c] : c; 1498 const PetscInt cind = c - cStart; 1499 const PetscBLASInt M = totDim, one = 1; 1500 const PetscScalar a = 1.0, b = 0.0; 1501 1502 PetscStackCallBLAS("BLASgemv", BLASgemv_("N", &M, &M, &a, &elemMat[cind*totDim*totDim], &M, &y[cind*totDim], &one, &b, z, &one)); 1503 if (mesh->printFEM > 1) { 1504 ierr = DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[cind*totDim*totDim]);CHKERRQ(ierr); 1505 ierr = DMPrintCellVector(c, "Y", totDim, &y[cind*totDim]);CHKERRQ(ierr); 1506 ierr = DMPrintCellVector(c, "Z", totDim, z);CHKERRQ(ierr); 1507 } 1508 ierr = DMPlexVecSetClosure(dm, section, Z, cell, z, ADD_VALUES);CHKERRQ(ierr); 1509 } 1510 ierr = PetscFree6(u,u_t,elemMat,elemMatD,y,z);CHKERRQ(ierr); 1511 if (mesh->printFEM) { 1512 ierr = PetscPrintf(PetscObjectComm((PetscObject)Z), "Z:\n");CHKERRQ(ierr); 1513 ierr = VecView(Z, NULL);CHKERRQ(ierr); 1514 } 1515 ierr = PetscFree(a);CHKERRQ(ierr); 1516 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 1517 ierr = DMDestroy(&plexAux);CHKERRQ(ierr); 1518 ierr = DMDestroy(&plex);CHKERRQ(ierr); 1519 ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 1520 PetscFunctionReturn(0); 1521 } 1522 1523 /*@ 1524 DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user. 1525 1526 Input Parameters: 1527 + dm - The mesh 1528 . X - Local input vector 1529 - user - The user context 1530 1531 Output Parameter: 1532 . Jac - Jacobian matrix 1533 1534 Note: 1535 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1536 like a GPU, or vectorize on a multicore machine. 1537 1538 Level: developer 1539 1540 .seealso: FormFunctionLocal() 1541 @*/ 1542 PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user) 1543 { 1544 DM plex; 1545 IS allcellIS; 1546 PetscBool hasJac, hasPrec; 1547 PetscInt Nds, s; 1548 PetscErrorCode ierr; 1549 1550 PetscFunctionBegin; 1551 ierr = DMSNESConvertPlex(dm, &plex, PETSC_TRUE);CHKERRQ(ierr); 1552 ierr = DMPlexGetAllCells_Internal(plex, &allcellIS);CHKERRQ(ierr); 1553 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr); 1554 for (s = 0; s < Nds; ++s) { 1555 PetscDS ds; 1556 IS cellIS; 1557 PetscHashFormKey key; 1558 1559 ierr = DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);CHKERRQ(ierr); 1560 key.value = 0; 1561 key.field = 0; 1562 if (!key.label) { 1563 ierr = PetscObjectReference((PetscObject) allcellIS);CHKERRQ(ierr); 1564 cellIS = allcellIS; 1565 } else { 1566 IS pointIS; 1567 1568 key.value = 1; 1569 ierr = DMLabelGetStratumIS(key.label, key.value, &pointIS);CHKERRQ(ierr); 1570 ierr = ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);CHKERRQ(ierr); 1571 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1572 } 1573 if (!s) { 1574 ierr = PetscDSHasJacobian(ds, &hasJac);CHKERRQ(ierr); 1575 ierr = PetscDSHasJacobianPreconditioner(ds, &hasPrec);CHKERRQ(ierr); 1576 if (hasJac && hasPrec) {ierr = MatZeroEntries(Jac);CHKERRQ(ierr);} 1577 ierr = MatZeroEntries(JacP);CHKERRQ(ierr); 1578 } 1579 ierr = DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user);CHKERRQ(ierr); 1580 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 1581 } 1582 ierr = ISDestroy(&allcellIS);CHKERRQ(ierr); 1583 ierr = DMDestroy(&plex);CHKERRQ(ierr); 1584 PetscFunctionReturn(0); 1585 } 1586 1587 /* 1588 MatComputeNeumannOverlap - Computes an unassembled (Neumann) local overlapping Mat in nonlinear context. 1589 1590 Input Parameters: 1591 + X - SNES linearization point 1592 . ovl - index set of overlapping subdomains 1593 1594 Output Parameter: 1595 . J - unassembled (Neumann) local matrix 1596 1597 Level: intermediate 1598 1599 .seealso: DMCreateNeumannOverlap(), MATIS, PCHPDDMSetAuxiliaryMat() 1600 */ 1601 static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx) 1602 { 1603 SNES snes; 1604 Mat pJ; 1605 DM ovldm,origdm; 1606 DMSNES sdm; 1607 PetscErrorCode (*bfun)(DM,Vec,void*); 1608 PetscErrorCode (*jfun)(DM,Vec,Mat,Mat,void*); 1609 void *bctx,*jctx; 1610 PetscErrorCode ierr; 1611 1612 PetscFunctionBegin; 1613 ierr = PetscObjectQuery((PetscObject)ovl,"_DM_Overlap_HPDDM_MATIS",(PetscObject*)&pJ);CHKERRQ(ierr); 1614 if (!pJ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing overlapping Mat");CHKERRQ(ierr); 1615 ierr = PetscObjectQuery((PetscObject)ovl,"_DM_Original_HPDDM",(PetscObject*)&origdm);CHKERRQ(ierr); 1616 if (!origdm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing original DM");CHKERRQ(ierr); 1617 ierr = MatGetDM(pJ,&ovldm);CHKERRQ(ierr); 1618 ierr = DMSNESGetBoundaryLocal(origdm,&bfun,&bctx);CHKERRQ(ierr); 1619 ierr = DMSNESSetBoundaryLocal(ovldm,bfun,bctx);CHKERRQ(ierr); 1620 ierr = DMSNESGetJacobianLocal(origdm,&jfun,&jctx);CHKERRQ(ierr); 1621 ierr = DMSNESSetJacobianLocal(ovldm,jfun,jctx);CHKERRQ(ierr); 1622 ierr = PetscObjectQuery((PetscObject)ovl,"_DM_Overlap_HPDDM_SNES",(PetscObject*)&snes);CHKERRQ(ierr); 1623 if (!snes) { 1624 ierr = SNESCreate(PetscObjectComm((PetscObject)ovl),&snes);CHKERRQ(ierr); 1625 ierr = SNESSetDM(snes,ovldm);CHKERRQ(ierr); 1626 ierr = PetscObjectCompose((PetscObject)ovl,"_DM_Overlap_HPDDM_SNES",(PetscObject)snes);CHKERRQ(ierr); 1627 ierr = PetscObjectDereference((PetscObject)snes);CHKERRQ(ierr); 1628 } 1629 ierr = DMGetDMSNES(ovldm,&sdm);CHKERRQ(ierr); 1630 ierr = VecLockReadPush(X);CHKERRQ(ierr); 1631 PetscStackPush("SNES user Jacobian function"); 1632 ierr = (*sdm->ops->computejacobian)(snes,X,pJ,pJ,sdm->jacobianctx);CHKERRQ(ierr); 1633 PetscStackPop; 1634 ierr = VecLockReadPop(X);CHKERRQ(ierr); 1635 /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */ 1636 { 1637 Mat locpJ; 1638 1639 ierr = MatISGetLocalMat(pJ,&locpJ);CHKERRQ(ierr); 1640 ierr = MatCopy(locpJ,J,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 1641 } 1642 PetscFunctionReturn(0); 1643 } 1644 1645 /*@ 1646 DMPlexSetSNESLocalFEM - Use DMPlex's internal FEM routines to compute SNES boundary values, residual, and Jacobian. 1647 1648 Input Parameters: 1649 + dm - The DM object 1650 . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see PetscDSAddBoundary()) 1651 . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see PetscDSSetResidual()) 1652 - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see PetscDSSetJacobian()) 1653 1654 Level: developer 1655 @*/ 1656 PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx) 1657 { 1658 PetscErrorCode ierr; 1659 1660 PetscFunctionBegin; 1661 ierr = DMSNESSetBoundaryLocal(dm,DMPlexSNESComputeBoundaryFEM,boundaryctx);CHKERRQ(ierr); 1662 ierr = DMSNESSetFunctionLocal(dm,DMPlexSNESComputeResidualFEM,residualctx);CHKERRQ(ierr); 1663 ierr = DMSNESSetJacobianLocal(dm,DMPlexSNESComputeJacobianFEM,jacobianctx);CHKERRQ(ierr); 1664 ierr = PetscObjectComposeFunction((PetscObject)dm,"MatComputeNeumannOverlap_C",MatComputeNeumannOverlap_Plex);CHKERRQ(ierr); 1665 PetscFunctionReturn(0); 1666 } 1667 1668 /*@C 1669 DMSNESCheckDiscretization - Check the discretization error of the exact solution 1670 1671 Input Parameters: 1672 + snes - the SNES object 1673 . dm - the DM 1674 . t - the time 1675 . u - a DM vector 1676 - tol - A tolerance for the check, or -1 to print the results instead 1677 1678 Output Parameters: 1679 . error - An array which holds the discretization error in each field, or NULL 1680 1681 Note: The user must call PetscDSSetExactSolution() beforehand 1682 1683 Level: developer 1684 1685 .seealso: DNSNESCheckFromOptions(), DMSNESCheckResidual(), DMSNESCheckJacobian(), PetscDSSetExactSolution() 1686 @*/ 1687 PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[]) 1688 { 1689 PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 1690 void **ectxs; 1691 PetscReal *err; 1692 MPI_Comm comm; 1693 PetscInt Nf, f; 1694 PetscErrorCode ierr; 1695 1696 PetscFunctionBegin; 1697 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1698 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1699 PetscValidHeaderSpecific(u, VEC_CLASSID, 3); 1700 if (error) PetscValidRealPointer(error, 6); 1701 1702 ierr = DMComputeExactSolution(dm, t, u, NULL);CHKERRQ(ierr); 1703 ierr = VecViewFromOptions(u, NULL, "-vec_view");CHKERRQ(ierr); 1704 1705 ierr = PetscObjectGetComm((PetscObject) snes, &comm);CHKERRQ(ierr); 1706 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1707 ierr = PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err);CHKERRQ(ierr); 1708 { 1709 PetscInt Nds, s; 1710 1711 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr); 1712 for (s = 0; s < Nds; ++s) { 1713 PetscDS ds; 1714 DMLabel label; 1715 IS fieldIS; 1716 const PetscInt *fields; 1717 PetscInt dsNf, f; 1718 1719 ierr = DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds);CHKERRQ(ierr); 1720 ierr = PetscDSGetNumFields(ds, &dsNf);CHKERRQ(ierr); 1721 ierr = ISGetIndices(fieldIS, &fields);CHKERRQ(ierr); 1722 for (f = 0; f < dsNf; ++f) { 1723 const PetscInt field = fields[f]; 1724 ierr = PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]);CHKERRQ(ierr); 1725 } 1726 ierr = ISRestoreIndices(fieldIS, &fields);CHKERRQ(ierr); 1727 } 1728 } 1729 if (Nf > 1) { 1730 ierr = DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err);CHKERRQ(ierr); 1731 if (tol >= 0.0) { 1732 for (f = 0; f < Nf; ++f) { 1733 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); 1734 } 1735 } else if (error) { 1736 for (f = 0; f < Nf; ++f) error[f] = err[f]; 1737 } else { 1738 ierr = PetscPrintf(comm, "L_2 Error: [");CHKERRQ(ierr); 1739 for (f = 0; f < Nf; ++f) { 1740 if (f) {ierr = PetscPrintf(comm, ", ");CHKERRQ(ierr);} 1741 ierr = PetscPrintf(comm, "%g", (double)err[f]);CHKERRQ(ierr); 1742 } 1743 ierr = PetscPrintf(comm, "]\n");CHKERRQ(ierr); 1744 } 1745 } else { 1746 ierr = DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]);CHKERRQ(ierr); 1747 if (tol >= 0.0) { 1748 if (err[0] > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double) err[0], (double) tol); 1749 } else if (error) { 1750 error[0] = err[0]; 1751 } else { 1752 ierr = PetscPrintf(comm, "L_2 Error: %g\n", (double) err[0]);CHKERRQ(ierr); 1753 } 1754 } 1755 ierr = PetscFree3(exacts, ectxs, err);CHKERRQ(ierr); 1756 PetscFunctionReturn(0); 1757 } 1758 1759 /*@C 1760 DMSNESCheckResidual - Check the residual of the exact solution 1761 1762 Input Parameters: 1763 + snes - the SNES object 1764 . dm - the DM 1765 . u - a DM vector 1766 - tol - A tolerance for the check, or -1 to print the results instead 1767 1768 Output Parameters: 1769 . residual - The residual norm of the exact solution, or NULL 1770 1771 Level: developer 1772 1773 .seealso: DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckJacobian() 1774 @*/ 1775 PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual) 1776 { 1777 MPI_Comm comm; 1778 Vec r; 1779 PetscReal res; 1780 PetscErrorCode ierr; 1781 1782 PetscFunctionBegin; 1783 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1784 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1785 PetscValidHeaderSpecific(u, VEC_CLASSID, 3); 1786 if (residual) PetscValidRealPointer(residual, 5); 1787 ierr = PetscObjectGetComm((PetscObject) snes, &comm);CHKERRQ(ierr); 1788 ierr = DMComputeExactSolution(dm, 0.0, u, NULL);CHKERRQ(ierr); 1789 ierr = VecDuplicate(u, &r);CHKERRQ(ierr); 1790 ierr = SNESComputeFunction(snes, u, r);CHKERRQ(ierr); 1791 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr); 1792 if (tol >= 0.0) { 1793 if (res > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double) res, (double) tol); 1794 } else if (residual) { 1795 *residual = res; 1796 } else { 1797 ierr = PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);CHKERRQ(ierr); 1798 ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr); 1799 ierr = PetscObjectSetName((PetscObject) r, "Initial Residual");CHKERRQ(ierr); 1800 ierr = PetscObjectSetOptionsPrefix((PetscObject)r,"res_");CHKERRQ(ierr); 1801 ierr = VecViewFromOptions(r, NULL, "-vec_view");CHKERRQ(ierr); 1802 } 1803 ierr = VecDestroy(&r);CHKERRQ(ierr); 1804 PetscFunctionReturn(0); 1805 } 1806 1807 /*@C 1808 DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test 1809 1810 Input Parameters: 1811 + snes - the SNES object 1812 . dm - the DM 1813 . u - a DM vector 1814 - tol - A tolerance for the check, or -1 to print the results instead 1815 1816 Output Parameters: 1817 + isLinear - Flag indicaing that the function looks linear, or NULL 1818 - convRate - The rate of convergence of the linear model, or NULL 1819 1820 Level: developer 1821 1822 .seealso: DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckResidual() 1823 @*/ 1824 PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate) 1825 { 1826 MPI_Comm comm; 1827 PetscDS ds; 1828 Mat J, M; 1829 MatNullSpace nullspace; 1830 PetscReal slope, intercept; 1831 PetscBool hasJac, hasPrec, isLin = PETSC_FALSE; 1832 PetscErrorCode ierr; 1833 1834 PetscFunctionBegin; 1835 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1836 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1837 PetscValidHeaderSpecific(u, VEC_CLASSID, 3); 1838 if (isLinear) PetscValidBoolPointer(isLinear, 5); 1839 if (convRate) PetscValidRealPointer(convRate, 5); 1840 ierr = PetscObjectGetComm((PetscObject) snes, &comm);CHKERRQ(ierr); 1841 ierr = DMComputeExactSolution(dm, 0.0, u, NULL);CHKERRQ(ierr); 1842 /* Create and view matrices */ 1843 ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr); 1844 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 1845 ierr = PetscDSHasJacobian(ds, &hasJac);CHKERRQ(ierr); 1846 ierr = PetscDSHasJacobianPreconditioner(ds, &hasPrec);CHKERRQ(ierr); 1847 if (hasJac && hasPrec) { 1848 ierr = DMCreateMatrix(dm, &M);CHKERRQ(ierr); 1849 ierr = SNESComputeJacobian(snes, u, J, M);CHKERRQ(ierr); 1850 ierr = PetscObjectSetName((PetscObject) M, "Preconditioning Matrix");CHKERRQ(ierr); 1851 ierr = PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");CHKERRQ(ierr); 1852 ierr = MatViewFromOptions(M, NULL, "-mat_view");CHKERRQ(ierr); 1853 ierr = MatDestroy(&M);CHKERRQ(ierr); 1854 } else { 1855 ierr = SNESComputeJacobian(snes, u, J, J);CHKERRQ(ierr); 1856 } 1857 ierr = PetscObjectSetName((PetscObject) J, "Jacobian");CHKERRQ(ierr); 1858 ierr = PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");CHKERRQ(ierr); 1859 ierr = MatViewFromOptions(J, NULL, "-mat_view");CHKERRQ(ierr); 1860 /* Check nullspace */ 1861 ierr = MatGetNullSpace(J, &nullspace);CHKERRQ(ierr); 1862 if (nullspace) { 1863 PetscBool isNull; 1864 ierr = MatNullSpaceTest(nullspace, J, &isNull);CHKERRQ(ierr); 1865 if (!isNull) SETERRQ(comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid."); 1866 } 1867 /* Taylor test */ 1868 { 1869 PetscRandom rand; 1870 Vec du, uhat, r, rhat, df; 1871 PetscReal h; 1872 PetscReal *es, *hs, *errors; 1873 PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1; 1874 PetscInt Nv, v; 1875 1876 /* Choose a perturbation direction */ 1877 ierr = PetscRandomCreate(comm, &rand);CHKERRQ(ierr); 1878 ierr = VecDuplicate(u, &du);CHKERRQ(ierr); 1879 ierr = VecSetRandom(du, rand);CHKERRQ(ierr); 1880 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr); 1881 ierr = VecDuplicate(u, &df);CHKERRQ(ierr); 1882 ierr = MatMult(J, du, df);CHKERRQ(ierr); 1883 /* Evaluate residual at u, F(u), save in vector r */ 1884 ierr = VecDuplicate(u, &r);CHKERRQ(ierr); 1885 ierr = SNESComputeFunction(snes, u, r);CHKERRQ(ierr); 1886 /* Look at the convergence of our Taylor approximation as we approach u */ 1887 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv); 1888 ierr = PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);CHKERRQ(ierr); 1889 ierr = VecDuplicate(u, &uhat);CHKERRQ(ierr); 1890 ierr = VecDuplicate(u, &rhat);CHKERRQ(ierr); 1891 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) { 1892 ierr = VecWAXPY(uhat, h, du, u);CHKERRQ(ierr); 1893 /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */ 1894 ierr = SNESComputeFunction(snes, uhat, rhat);CHKERRQ(ierr); 1895 ierr = VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);CHKERRQ(ierr); 1896 ierr = VecNorm(rhat, NORM_2, &errors[Nv]);CHKERRQ(ierr); 1897 1898 es[Nv] = PetscLog10Real(errors[Nv]); 1899 hs[Nv] = PetscLog10Real(h); 1900 } 1901 ierr = VecDestroy(&uhat);CHKERRQ(ierr); 1902 ierr = VecDestroy(&rhat);CHKERRQ(ierr); 1903 ierr = VecDestroy(&df);CHKERRQ(ierr); 1904 ierr = VecDestroy(&r);CHKERRQ(ierr); 1905 ierr = VecDestroy(&du);CHKERRQ(ierr); 1906 for (v = 0; v < Nv; ++v) { 1907 if ((tol >= 0) && (errors[v] > tol)) break; 1908 else if (errors[v] > PETSC_SMALL) break; 1909 } 1910 if (v == Nv) isLin = PETSC_TRUE; 1911 ierr = PetscLinearRegression(Nv, hs, es, &slope, &intercept);CHKERRQ(ierr); 1912 ierr = PetscFree3(es, hs, errors);CHKERRQ(ierr); 1913 /* Slope should be about 2 */ 1914 if (tol >= 0) { 1915 if (!isLin && PetscAbsReal(2 - slope) > tol) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double) slope); 1916 } else if (isLinear || convRate) { 1917 if (isLinear) *isLinear = isLin; 1918 if (convRate) *convRate = slope; 1919 } else { 1920 if (!isLin) {ierr = PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double) slope);CHKERRQ(ierr);} 1921 else {ierr = PetscPrintf(comm, "Function appears to be linear\n");CHKERRQ(ierr);} 1922 } 1923 } 1924 ierr = MatDestroy(&J);CHKERRQ(ierr); 1925 PetscFunctionReturn(0); 1926 } 1927 1928 PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u) 1929 { 1930 PetscErrorCode ierr; 1931 1932 PetscFunctionBegin; 1933 ierr = DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL);CHKERRQ(ierr); 1934 ierr = DMSNESCheckResidual(snes, dm, u, -1.0, NULL);CHKERRQ(ierr); 1935 ierr = DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL);CHKERRQ(ierr); 1936 PetscFunctionReturn(0); 1937 } 1938 1939 /*@C 1940 DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information 1941 1942 Input Parameters: 1943 + snes - the SNES object 1944 - u - representative SNES vector 1945 1946 Note: The user must call PetscDSSetExactSolution() beforehand 1947 1948 Level: developer 1949 @*/ 1950 PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u) 1951 { 1952 DM dm; 1953 Vec sol; 1954 PetscBool check; 1955 PetscErrorCode ierr; 1956 1957 PetscFunctionBegin; 1958 ierr = PetscOptionsHasName(((PetscObject)snes)->options,((PetscObject)snes)->prefix, "-dmsnes_check", &check);CHKERRQ(ierr); 1959 if (!check) PetscFunctionReturn(0); 1960 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 1961 ierr = VecDuplicate(u, &sol);CHKERRQ(ierr); 1962 ierr = SNESSetSolution(snes, sol);CHKERRQ(ierr); 1963 ierr = DMSNESCheck_Internal(snes, dm, sol);CHKERRQ(ierr); 1964 ierr = VecDestroy(&sol);CHKERRQ(ierr); 1965 PetscFunctionReturn(0); 1966 } 1967