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