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