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