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