Actual source code: dmplexsnes.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/snesimpl.h>
  3: #include <petscds.h>
  4: #include <petsc/private/petscimpl.h>
  5: #include <petsc/private/petscfeimpl.h>

  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: }

 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.

 18:   Collective on SNES

 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

 27:   Output Parameter:
 28: . u         - The solution with a continuum pressure integral of zero

 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.

 33:   Level: developer

 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;

 46:   PetscObjectGetComm((PetscObject) snes, &comm);
 47:   SNESGetDM(snes, &dm);
 50:   DMGetDS(dm, &ds);
 51:   PetscDSSetObjective(ds, pfield, pressure_Private);
 52:   MatNullSpaceGetVecs(nullspace, NULL, &Nv, &nullvecs);
 54:   VecDot(nullvecs[0], u, &pintd);
 56:   PetscDSGetNumFields(ds, &Nf);
 57:   PetscMalloc2(Nf, &intc, Nf, &intn);
 58:   DMPlexComputeIntegralFEM(dm, nullvecs[0], intn, ctx);
 59:   DMPlexComputeIntegralFEM(dm, u, intc, ctx);
 60:   VecAXPY(u, -intc[pfield]/intn[pfield], nullvecs[0]);
 61: #if defined (PETSC_USE_DEBUG)
 62:   DMPlexComputeIntegralFEM(dm, u, intc, ctx);
 64: #endif
 65:   PetscFree2(intc, intn);
 66:   return 0;
 67: }

 69: /*@C
 70:    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().

 72:    Logically Collective on SNES

 74:    Input Parameters:
 75: +  snes - the SNES context
 76: .  it - the iteration (0 indicates before any Newton steps)
 77: .  xnorm - 2-norm of current iterate
 78: .  snorm - 2-norm of current step
 79: .  fnorm - 2-norm of function at current iterate
 80: -  ctx   - Optional user context

 82:    Output Parameter:
 83: .  reason  - SNES_CONVERGED_ITERATING, SNES_CONVERGED_ITS, or SNES_DIVERGED_FNORM_NAN

 85:    Notes:
 86:    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.

 88:    Level: advanced

 90: .seealso: SNESConvergedDefault(), SNESSetConvergenceTest(), DMSetNullSpaceConstructor()
 91: @*/
 92: PetscErrorCode SNESConvergedCorrectPressure(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal f, SNESConvergedReason *reason, void *ctx)
 93: {
 94:   PetscBool      monitorIntegral = PETSC_FALSE;

 96:   SNESConvergedDefault(snes, it, xnorm, gnorm, f, reason, ctx);
 97:   if (monitorIntegral) {
 98:     Mat          J;
 99:     Vec          u;
100:     MatNullSpace nullspace;
101:     const Vec   *nullvecs;
102:     PetscScalar  pintd;

104:     SNESGetSolution(snes, &u);
105:     SNESGetJacobian(snes, &J, NULL, NULL, NULL);
106:     MatGetNullSpace(J, &nullspace);
107:     MatNullSpaceGetVecs(nullspace, NULL, NULL, &nullvecs);
108:     VecDot(nullvecs[0], u, &pintd);
109:     PetscInfo(snes, "SNES: Discrete integral of pressure: %g\n", (double) PetscRealPart(pintd));
110:   }
111:   if (*reason > 0) {
112:     Mat          J;
113:     Vec          u;
114:     MatNullSpace nullspace;
115:     PetscInt     pfield = 1;

117:     SNESGetSolution(snes, &u);
118:     SNESGetJacobian(snes, &J, NULL, NULL, NULL);
119:     MatGetNullSpace(J, &nullspace);
120:     SNESCorrectDiscretePressure_Private(snes, pfield, nullspace, u, ctx);
121:   }
122:   return 0;
123: }

125: /************************** Interpolation *******************************/

127: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
128: {
129:   PetscBool      isPlex;

131:   PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
132:   if (isPlex) {
133:     *plex = dm;
134:     PetscObjectReference((PetscObject) dm);
135:   } else {
136:     PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
137:     if (!*plex) {
138:       DMConvert(dm,DMPLEX,plex);
139:       PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
140:       if (copy) {
141:         DMCopyDMSNES(dm, *plex);
142:         DMCopyAuxiliaryVec(dm, *plex);
143:       }
144:     } else {
145:       PetscObjectReference((PetscObject) *plex);
146:     }
147:   }
148:   return 0;
149: }

151: /*@C
152:   DMInterpolationCreate - Creates a DMInterpolationInfo context

154:   Collective

156:   Input Parameter:
157: . comm - the communicator

159:   Output Parameter:
160: . ctx - the context

162:   Level: beginner

164: .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationDestroy()
165: @*/
166: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
167: {
169:   PetscNew(ctx);

171:   (*ctx)->comm   = comm;
172:   (*ctx)->dim    = -1;
173:   (*ctx)->nInput = 0;
174:   (*ctx)->points = NULL;
175:   (*ctx)->cells  = NULL;
176:   (*ctx)->n      = -1;
177:   (*ctx)->coords = NULL;
178:   return 0;
179: }

181: /*@C
182:   DMInterpolationSetDim - Sets the spatial dimension for the interpolation context

184:   Not collective

186:   Input Parameters:
187: + ctx - the context
188: - dim - the spatial dimension

190:   Level: intermediate

192: .seealso: DMInterpolationGetDim(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
193: @*/
194: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
195: {
197:   ctx->dim = dim;
198:   return 0;
199: }

201: /*@C
202:   DMInterpolationGetDim - Gets the spatial dimension for the interpolation context

204:   Not collective

206:   Input Parameter:
207: . ctx - the context

209:   Output Parameter:
210: . dim - the spatial dimension

212:   Level: intermediate

214: .seealso: DMInterpolationSetDim(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
215: @*/
216: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
217: {
219:   *dim = ctx->dim;
220:   return 0;
221: }

223: /*@C
224:   DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context

226:   Not collective

228:   Input Parameters:
229: + ctx - the context
230: - dof - the number of fields

232:   Level: intermediate

234: .seealso: DMInterpolationGetDof(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
235: @*/
236: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
237: {
239:   ctx->dof = dof;
240:   return 0;
241: }

243: /*@C
244:   DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context

246:   Not collective

248:   Input Parameter:
249: . ctx - the context

251:   Output Parameter:
252: . dof - the number of fields

254:   Level: intermediate

256: .seealso: DMInterpolationSetDof(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
257: @*/
258: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
259: {
261:   *dof = ctx->dof;
262:   return 0;
263: }

265: /*@C
266:   DMInterpolationAddPoints - Add points at which we will interpolate the fields

268:   Not collective

270:   Input Parameters:
271: + ctx    - the context
272: . n      - the number of points
273: - points - the coordinates for each point, an array of size n * dim

275:   Note: The coordinate information is copied.

277:   Level: intermediate

279: .seealso: DMInterpolationSetDim(), DMInterpolationEvaluate(), DMInterpolationCreate()
280: @*/
281: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
282: {
285:   ctx->nInput = n;

287:   PetscMalloc1(n*ctx->dim, &ctx->points);
288:   PetscArraycpy(ctx->points, points, n*ctx->dim);
289:   return 0;
290: }

292: /*@C
293:   DMInterpolationSetUp - Compute spatial indices for point location during interpolation

295:   Collective on ctx

297:   Input Parameters:
298: + ctx - the context
299: . dm  - the DM for the function space used for interpolation
300: . redundantPoints - If PETSC_TRUE, all processes are passing in the same array of points. Otherwise, points need to be communicated among processes.
301: - ignoreOutsideDomain - If PETSC_TRUE, ignore points outside the domain, otherwise return an error

303:   Level: intermediate

305: .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
306: @*/
307: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain)
308: {
309:   MPI_Comm          comm = ctx->comm;
310:   PetscScalar       *a;
311:   PetscInt          p, q, i;
312:   PetscMPIInt       rank, size;
313:   Vec               pointVec;
314:   PetscSF           cellSF;
315:   PetscLayout       layout;
316:   PetscReal         *globalPoints;
317:   PetscScalar       *globalPointsScalar;
318:   const PetscInt    *ranges;
319:   PetscMPIInt       *counts, *displs;
320:   const PetscSFNode *foundCells;
321:   const PetscInt    *foundPoints;
322:   PetscMPIInt       *foundProcs, *globalProcs;
323:   PetscInt          n, N, numFound;

326:   MPI_Comm_size(comm, &size);
327:   MPI_Comm_rank(comm, &rank);
329:   /* Locate points */
330:   n = ctx->nInput;
331:   if (!redundantPoints) {
332:     PetscLayoutCreate(comm, &layout);
333:     PetscLayoutSetBlockSize(layout, 1);
334:     PetscLayoutSetLocalSize(layout, n);
335:     PetscLayoutSetUp(layout);
336:     PetscLayoutGetSize(layout, &N);
337:     /* Communicate all points to all processes */
338:     PetscMalloc3(N*ctx->dim,&globalPoints,size,&counts,size,&displs);
339:     PetscLayoutGetRanges(layout, &ranges);
340:     for (p = 0; p < size; ++p) {
341:       counts[p] = (ranges[p+1] - ranges[p])*ctx->dim;
342:       displs[p] = ranges[p]*ctx->dim;
343:     }
344:     MPI_Allgatherv(ctx->points, n*ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);
345:   } else {
346:     N = n;
347:     globalPoints = ctx->points;
348:     counts = displs = NULL;
349:     layout = NULL;
350:   }
351: #if 0
352:   PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);
353:   /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
354: #else
355: #if defined(PETSC_USE_COMPLEX)
356:   PetscMalloc1(N*ctx->dim,&globalPointsScalar);
357:   for (i=0; i<N*ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
358: #else
359:   globalPointsScalar = globalPoints;
360: #endif
361:   VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N*ctx->dim, globalPointsScalar, &pointVec);
362:   PetscMalloc2(N,&foundProcs,N,&globalProcs);
363:   for (p = 0; p < N; ++p) {foundProcs[p] = size;}
364:   cellSF = NULL;
365:   DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF);
366:   PetscSFGetGraph(cellSF,NULL,&numFound,&foundPoints,&foundCells);
367: #endif
368:   for (p = 0; p < numFound; ++p) {
369:     if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
370:   }
371:   /* Let the lowest rank process own each point */
372:   MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);
373:   ctx->n = 0;
374:   for (p = 0; p < N; ++p) {
375:     if (globalProcs[p] == size) {
377:       else if (rank == 0) ++ctx->n;
378:     } else if (globalProcs[p] == rank) ++ctx->n;
379:   }
380:   /* Create coordinates vector and array of owned cells */
381:   PetscMalloc1(ctx->n, &ctx->cells);
382:   VecCreate(comm, &ctx->coords);
383:   VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);
384:   VecSetBlockSize(ctx->coords, ctx->dim);
385:   VecSetType(ctx->coords,VECSTANDARD);
386:   VecGetArray(ctx->coords, &a);
387:   for (p = 0, q = 0, i = 0; p < N; ++p) {
388:     if (globalProcs[p] == rank) {
389:       PetscInt d;

391:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p*ctx->dim+d];
392:       ctx->cells[q] = foundCells[q].index;
393:       ++q;
394:     }
395:     if (globalProcs[p] == size && rank == 0) {
396:       PetscInt d;

398:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = 0.;
399:       ctx->cells[q] = -1;
400:       ++q;
401:     }
402:   }
403:   VecRestoreArray(ctx->coords, &a);
404: #if 0
405:   PetscFree3(foundCells,foundProcs,globalProcs);
406: #else
407:   PetscFree2(foundProcs,globalProcs);
408:   PetscSFDestroy(&cellSF);
409:   VecDestroy(&pointVec);
410: #endif
411:   if ((void*)globalPointsScalar != (void*)globalPoints) PetscFree(globalPointsScalar);
412:   if (!redundantPoints) PetscFree3(globalPoints,counts,displs);
413:   PetscLayoutDestroy(&layout);
414:   return 0;
415: }

417: /*@C
418:   DMInterpolationGetCoordinates - Gets a Vec with the coordinates of each interpolation point

420:   Collective on ctx

422:   Input Parameter:
423: . ctx - the context

425:   Output Parameter:
426: . coordinates  - the coordinates of interpolation points

428:   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.

430:   Level: intermediate

432: .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
433: @*/
434: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
435: {
438:   *coordinates = ctx->coords;
439:   return 0;
440: }

442: /*@C
443:   DMInterpolationGetVector - Gets a Vec which can hold all the interpolated field values

445:   Collective on ctx

447:   Input Parameter:
448: . ctx - the context

450:   Output Parameter:
451: . v  - a vector capable of holding the interpolated field values

453:   Note: This vector should be returned using DMInterpolationRestoreVector().

455:   Level: intermediate

457: .seealso: DMInterpolationRestoreVector(), DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
458: @*/
459: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
460: {
463:   VecCreate(ctx->comm, v);
464:   VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);
465:   VecSetBlockSize(*v, ctx->dof);
466:   VecSetType(*v,VECSTANDARD);
467:   return 0;
468: }

470: /*@C
471:   DMInterpolationRestoreVector - Returns a Vec which can hold all the interpolated field values

473:   Collective on ctx

475:   Input Parameters:
476: + ctx - the context
477: - v  - a vector capable of holding the interpolated field values

479:   Level: intermediate

481: .seealso: DMInterpolationGetVector(), DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
482: @*/
483: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
484: {
487:   VecDestroy(v);
488:   return 0;
489: }

491: static inline PetscErrorCode DMInterpolate_Segment_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
492: {
493:   PetscReal          v0, J, invJ, detJ;
494:   const PetscInt     dof = ctx->dof;
495:   const PetscScalar *coords;
496:   PetscScalar       *a;
497:   PetscInt           p;

499:   VecGetArrayRead(ctx->coords, &coords);
500:   VecGetArray(v, &a);
501:   for (p = 0; p < ctx->n; ++p) {
502:     PetscInt     c = ctx->cells[p];
503:     PetscScalar *x = NULL;
504:     PetscReal    xir[1];
505:     PetscInt     xSize, comp;

507:     DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ);
509:     xir[0] = invJ*PetscRealPart(coords[p] - v0);
510:     DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
511:     if (2*dof == xSize) {
512:       for (comp = 0; comp < dof; ++comp) a[p*dof+comp] = x[0*dof+comp]*(1 - xir[0]) + x[1*dof+comp]*xir[0];
513:     } else if (dof == xSize) {
514:       for (comp = 0; comp < dof; ++comp) a[p*dof+comp] = x[0*dof+comp];
515:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Input closure size %" PetscInt_FMT " must be either %" PetscInt_FMT " or %" PetscInt_FMT, xSize, 2*dof, dof);
516:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
517:   }
518:   VecRestoreArray(v, &a);
519:   VecRestoreArrayRead(ctx->coords, &coords);
520:   return 0;
521: }

523: static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
524: {
525:   PetscReal      *v0, *J, *invJ, detJ;
526:   const PetscScalar *coords;
527:   PetscScalar    *a;
528:   PetscInt       p;

530:   PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
531:   VecGetArrayRead(ctx->coords, &coords);
532:   VecGetArray(v, &a);
533:   for (p = 0; p < ctx->n; ++p) {
534:     PetscInt     c = ctx->cells[p];
535:     PetscScalar *x = NULL;
536:     PetscReal    xi[4];
537:     PetscInt     d, f, comp;

539:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
541:     DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
542:     for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];

544:     for (d = 0; d < ctx->dim; ++d) {
545:       xi[d] = 0.0;
546:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
547:       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];
548:     }
549:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
550:   }
551:   VecRestoreArray(v, &a);
552:   VecRestoreArrayRead(ctx->coords, &coords);
553:   PetscFree3(v0, J, invJ);
554:   return 0;
555: }

557: static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
558: {
559:   PetscReal      *v0, *J, *invJ, detJ;
560:   const PetscScalar *coords;
561:   PetscScalar    *a;
562:   PetscInt       p;

564:   PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
565:   VecGetArrayRead(ctx->coords, &coords);
566:   VecGetArray(v, &a);
567:   for (p = 0; p < ctx->n; ++p) {
568:     PetscInt       c = ctx->cells[p];
569:     const PetscInt order[3] = {2, 1, 3};
570:     PetscScalar   *x = NULL;
571:     PetscReal      xi[4];
572:     PetscInt       d, f, comp;

574:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
576:     DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
577:     for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];

579:     for (d = 0; d < ctx->dim; ++d) {
580:       xi[d] = 0.0;
581:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
582:       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];
583:     }
584:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
585:   }
586:   VecRestoreArray(v, &a);
587:   VecRestoreArrayRead(ctx->coords, &coords);
588:   PetscFree3(v0, J, invJ);
589:   return 0;
590: }

592: static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
593: {
594:   const PetscScalar *vertices = (const PetscScalar*) ctx;
595:   const PetscScalar x0        = vertices[0];
596:   const PetscScalar y0        = vertices[1];
597:   const PetscScalar x1        = vertices[2];
598:   const PetscScalar y1        = vertices[3];
599:   const PetscScalar x2        = vertices[4];
600:   const PetscScalar y2        = vertices[5];
601:   const PetscScalar x3        = vertices[6];
602:   const PetscScalar y3        = vertices[7];
603:   const PetscScalar f_1       = x1 - x0;
604:   const PetscScalar g_1       = y1 - y0;
605:   const PetscScalar f_3       = x3 - x0;
606:   const PetscScalar g_3       = y3 - y0;
607:   const PetscScalar f_01      = x2 - x1 - x3 + x0;
608:   const PetscScalar g_01      = y2 - y1 - y3 + y0;
609:   const PetscScalar *ref;
610:   PetscScalar       *real;

612:   VecGetArrayRead(Xref,  &ref);
613:   VecGetArray(Xreal, &real);
614:   {
615:     const PetscScalar p0 = ref[0];
616:     const PetscScalar p1 = ref[1];

618:     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
619:     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
620:   }
621:   PetscLogFlops(28);
622:   VecRestoreArrayRead(Xref,  &ref);
623:   VecRestoreArray(Xreal, &real);
624:   return 0;
625: }

627: #include <petsc/private/dmimpl.h>
628: static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
629: {
630:   const PetscScalar *vertices = (const PetscScalar*) ctx;
631:   const PetscScalar x0        = vertices[0];
632:   const PetscScalar y0        = vertices[1];
633:   const PetscScalar x1        = vertices[2];
634:   const PetscScalar y1        = vertices[3];
635:   const PetscScalar x2        = vertices[4];
636:   const PetscScalar y2        = vertices[5];
637:   const PetscScalar x3        = vertices[6];
638:   const PetscScalar y3        = vertices[7];
639:   const PetscScalar f_01      = x2 - x1 - x3 + x0;
640:   const PetscScalar g_01      = y2 - y1 - y3 + y0;
641:   const PetscScalar *ref;

643:   VecGetArrayRead(Xref,  &ref);
644:   {
645:     const PetscScalar x       = ref[0];
646:     const PetscScalar y       = ref[1];
647:     const PetscInt    rows[2] = {0, 1};
648:     PetscScalar       values[4];

650:     values[0] = (x1 - x0 + f_01*y) * 0.5; values[1] = (x3 - x0 + f_01*x) * 0.5;
651:     values[2] = (y1 - y0 + g_01*y) * 0.5; values[3] = (y3 - y0 + g_01*x) * 0.5;
652:     MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);
653:   }
654:   PetscLogFlops(30);
655:   VecRestoreArrayRead(Xref,  &ref);
656:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
657:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
658:   return 0;
659: }

661: static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
662: {
663:   DM                 dmCoord;
664:   PetscFE            fem = NULL;
665:   SNES               snes;
666:   KSP                ksp;
667:   PC                 pc;
668:   Vec                coordsLocal, r, ref, real;
669:   Mat                J;
670:   PetscTabulation    T = NULL;
671:   const PetscScalar *coords;
672:   PetscScalar        *a;
673:   PetscReal          xir[2] = {0., 0.};
674:   PetscInt           Nf, p;
675:   const PetscInt     dof = ctx->dof;

677:   DMGetNumFields(dm, &Nf);
678:   if (Nf) {
679:     PetscObject  obj;
680:     PetscClassId id;

682:     DMGetField(dm, 0, NULL, &obj);
683:     PetscObjectGetClassId(obj, &id);
684:     if (id == PETSCFE_CLASSID) {fem = (PetscFE) obj; PetscFECreateTabulation(fem, 1, 1, xir, 0, &T);}
685:   }
686:   DMGetCoordinatesLocal(dm, &coordsLocal);
687:   DMGetCoordinateDM(dm, &dmCoord);
688:   SNESCreate(PETSC_COMM_SELF, &snes);
689:   SNESSetOptionsPrefix(snes, "quad_interp_");
690:   VecCreate(PETSC_COMM_SELF, &r);
691:   VecSetSizes(r, 2, 2);
692:   VecSetType(r,dm->vectype);
693:   VecDuplicate(r, &ref);
694:   VecDuplicate(r, &real);
695:   MatCreate(PETSC_COMM_SELF, &J);
696:   MatSetSizes(J, 2, 2, 2, 2);
697:   MatSetType(J, MATSEQDENSE);
698:   MatSetUp(J);
699:   SNESSetFunction(snes, r, QuadMap_Private, NULL);
700:   SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);
701:   SNESGetKSP(snes, &ksp);
702:   KSPGetPC(ksp, &pc);
703:   PCSetType(pc, PCLU);
704:   SNESSetFromOptions(snes);

706:   VecGetArrayRead(ctx->coords, &coords);
707:   VecGetArray(v, &a);
708:   for (p = 0; p < ctx->n; ++p) {
709:     PetscScalar *x = NULL, *vertices = NULL;
710:     PetscScalar *xi;
711:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

713:     /* Can make this do all points at once */
714:     DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
716:     DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
717:     SNESSetFunction(snes, NULL, NULL, vertices);
718:     SNESSetJacobian(snes, NULL, NULL, NULL, vertices);
719:     VecGetArray(real, &xi);
720:     xi[0]  = coords[p*ctx->dim+0];
721:     xi[1]  = coords[p*ctx->dim+1];
722:     VecRestoreArray(real, &xi);
723:     SNESSolve(snes, real, ref);
724:     VecGetArray(ref, &xi);
725:     xir[0] = PetscRealPart(xi[0]);
726:     xir[1] = PetscRealPart(xi[1]);
727:     if (4*dof == xSize) {
728:       for (comp = 0; comp < dof; ++comp)
729:         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];
730:     } else if (dof == xSize) {
731:       for (comp = 0; comp < dof; ++comp) a[p*dof+comp] = x[0*dof+comp];
732:     } else {
733:       PetscInt d;

736:       xir[0] = 2.0*xir[0] - 1.0; xir[1] = 2.0*xir[1] - 1.0;
737:       PetscFEComputeTabulation(fem, 1, xir, 0, T);
738:       for (comp = 0; comp < dof; ++comp) {
739:         a[p*dof+comp] = 0.0;
740:         for (d = 0; d < xSize/dof; ++d) {
741:           a[p*dof+comp] += x[d*dof+comp]*T->T[0][d*dof+comp];
742:         }
743:       }
744:     }
745:     VecRestoreArray(ref, &xi);
746:     DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
747:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
748:   }
749:   PetscTabulationDestroy(&T);
750:   VecRestoreArray(v, &a);
751:   VecRestoreArrayRead(ctx->coords, &coords);

753:   SNESDestroy(&snes);
754:   VecDestroy(&r);
755:   VecDestroy(&ref);
756:   VecDestroy(&real);
757:   MatDestroy(&J);
758:   return 0;
759: }

761: static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
762: {
763:   const PetscScalar *vertices = (const PetscScalar*) ctx;
764:   const PetscScalar x0        = vertices[0];
765:   const PetscScalar y0        = vertices[1];
766:   const PetscScalar z0        = vertices[2];
767:   const PetscScalar x1        = vertices[9];
768:   const PetscScalar y1        = vertices[10];
769:   const PetscScalar z1        = vertices[11];
770:   const PetscScalar x2        = vertices[6];
771:   const PetscScalar y2        = vertices[7];
772:   const PetscScalar z2        = vertices[8];
773:   const PetscScalar x3        = vertices[3];
774:   const PetscScalar y3        = vertices[4];
775:   const PetscScalar z3        = vertices[5];
776:   const PetscScalar x4        = vertices[12];
777:   const PetscScalar y4        = vertices[13];
778:   const PetscScalar z4        = vertices[14];
779:   const PetscScalar x5        = vertices[15];
780:   const PetscScalar y5        = vertices[16];
781:   const PetscScalar z5        = vertices[17];
782:   const PetscScalar x6        = vertices[18];
783:   const PetscScalar y6        = vertices[19];
784:   const PetscScalar z6        = vertices[20];
785:   const PetscScalar x7        = vertices[21];
786:   const PetscScalar y7        = vertices[22];
787:   const PetscScalar z7        = vertices[23];
788:   const PetscScalar f_1       = x1 - x0;
789:   const PetscScalar g_1       = y1 - y0;
790:   const PetscScalar h_1       = z1 - z0;
791:   const PetscScalar f_3       = x3 - x0;
792:   const PetscScalar g_3       = y3 - y0;
793:   const PetscScalar h_3       = z3 - z0;
794:   const PetscScalar f_4       = x4 - x0;
795:   const PetscScalar g_4       = y4 - y0;
796:   const PetscScalar h_4       = z4 - z0;
797:   const PetscScalar f_01      = x2 - x1 - x3 + x0;
798:   const PetscScalar g_01      = y2 - y1 - y3 + y0;
799:   const PetscScalar h_01      = z2 - z1 - z3 + z0;
800:   const PetscScalar f_12      = x7 - x3 - x4 + x0;
801:   const PetscScalar g_12      = y7 - y3 - y4 + y0;
802:   const PetscScalar h_12      = z7 - z3 - z4 + z0;
803:   const PetscScalar f_02      = x5 - x1 - x4 + x0;
804:   const PetscScalar g_02      = y5 - y1 - y4 + y0;
805:   const PetscScalar h_02      = z5 - z1 - z4 + z0;
806:   const PetscScalar f_012     = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
807:   const PetscScalar g_012     = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
808:   const PetscScalar h_012     = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
809:   const PetscScalar *ref;
810:   PetscScalar       *real;

812:   VecGetArrayRead(Xref,  &ref);
813:   VecGetArray(Xreal, &real);
814:   {
815:     const PetscScalar p0 = ref[0];
816:     const PetscScalar p1 = ref[1];
817:     const PetscScalar p2 = ref[2];

819:     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;
820:     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;
821:     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;
822:   }
823:   PetscLogFlops(114);
824:   VecRestoreArrayRead(Xref,  &ref);
825:   VecRestoreArray(Xreal, &real);
826:   return 0;
827: }

829: static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
830: {
831:   const PetscScalar *vertices = (const PetscScalar*) ctx;
832:   const PetscScalar x0        = vertices[0];
833:   const PetscScalar y0        = vertices[1];
834:   const PetscScalar z0        = vertices[2];
835:   const PetscScalar x1        = vertices[9];
836:   const PetscScalar y1        = vertices[10];
837:   const PetscScalar z1        = vertices[11];
838:   const PetscScalar x2        = vertices[6];
839:   const PetscScalar y2        = vertices[7];
840:   const PetscScalar z2        = vertices[8];
841:   const PetscScalar x3        = vertices[3];
842:   const PetscScalar y3        = vertices[4];
843:   const PetscScalar z3        = vertices[5];
844:   const PetscScalar x4        = vertices[12];
845:   const PetscScalar y4        = vertices[13];
846:   const PetscScalar z4        = vertices[14];
847:   const PetscScalar x5        = vertices[15];
848:   const PetscScalar y5        = vertices[16];
849:   const PetscScalar z5        = vertices[17];
850:   const PetscScalar x6        = vertices[18];
851:   const PetscScalar y6        = vertices[19];
852:   const PetscScalar z6        = vertices[20];
853:   const PetscScalar x7        = vertices[21];
854:   const PetscScalar y7        = vertices[22];
855:   const PetscScalar z7        = vertices[23];
856:   const PetscScalar f_xy      = x2 - x1 - x3 + x0;
857:   const PetscScalar g_xy      = y2 - y1 - y3 + y0;
858:   const PetscScalar h_xy      = z2 - z1 - z3 + z0;
859:   const PetscScalar f_yz      = x7 - x3 - x4 + x0;
860:   const PetscScalar g_yz      = y7 - y3 - y4 + y0;
861:   const PetscScalar h_yz      = z7 - z3 - z4 + z0;
862:   const PetscScalar f_xz      = x5 - x1 - x4 + x0;
863:   const PetscScalar g_xz      = y5 - y1 - y4 + y0;
864:   const PetscScalar h_xz      = z5 - z1 - z4 + z0;
865:   const PetscScalar f_xyz     = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
866:   const PetscScalar g_xyz     = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
867:   const PetscScalar h_xyz     = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
868:   const PetscScalar *ref;

870:   VecGetArrayRead(Xref,  &ref);
871:   {
872:     const PetscScalar x       = ref[0];
873:     const PetscScalar y       = ref[1];
874:     const PetscScalar z       = ref[2];
875:     const PetscInt    rows[3] = {0, 1, 2};
876:     PetscScalar       values[9];

878:     values[0] = (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0;
879:     values[1] = (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0;
880:     values[2] = (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0;
881:     values[3] = (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0;
882:     values[4] = (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0;
883:     values[5] = (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0;
884:     values[6] = (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0;
885:     values[7] = (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0;
886:     values[8] = (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0;

888:     MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);
889:   }
890:   PetscLogFlops(152);
891:   VecRestoreArrayRead(Xref,  &ref);
892:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
893:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
894:   return 0;
895: }

897: static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
898: {
899:   DM             dmCoord;
900:   SNES           snes;
901:   KSP            ksp;
902:   PC             pc;
903:   Vec            coordsLocal, r, ref, real;
904:   Mat            J;
905:   const PetscScalar *coords;
906:   PetscScalar    *a;
907:   PetscInt       p;

909:   DMGetCoordinatesLocal(dm, &coordsLocal);
910:   DMGetCoordinateDM(dm, &dmCoord);
911:   SNESCreate(PETSC_COMM_SELF, &snes);
912:   SNESSetOptionsPrefix(snes, "hex_interp_");
913:   VecCreate(PETSC_COMM_SELF, &r);
914:   VecSetSizes(r, 3, 3);
915:   VecSetType(r,dm->vectype);
916:   VecDuplicate(r, &ref);
917:   VecDuplicate(r, &real);
918:   MatCreate(PETSC_COMM_SELF, &J);
919:   MatSetSizes(J, 3, 3, 3, 3);
920:   MatSetType(J, MATSEQDENSE);
921:   MatSetUp(J);
922:   SNESSetFunction(snes, r, HexMap_Private, NULL);
923:   SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);
924:   SNESGetKSP(snes, &ksp);
925:   KSPGetPC(ksp, &pc);
926:   PCSetType(pc, PCLU);
927:   SNESSetFromOptions(snes);

929:   VecGetArrayRead(ctx->coords, &coords);
930:   VecGetArray(v, &a);
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;

937:     /* Can make this do all points at once */
938:     DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
940:     DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
942:     SNESSetFunction(snes, NULL, NULL, vertices);
943:     SNESSetJacobian(snes, NULL, NULL, NULL, vertices);
944:     VecGetArray(real, &xi);
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:     VecRestoreArray(real, &xi);
949:     SNESSolve(snes, real, ref);
950:     VecGetArray(ref, &xi);
951:     xir[0] = PetscRealPart(xi[0]);
952:     xir[1] = PetscRealPart(xi[1]);
953:     xir[2] = PetscRealPart(xi[2]);
954:     if (8*ctx->dof == xSize) {
955:       for (comp = 0; comp < ctx->dof; ++comp) {
956:         a[p*ctx->dof+comp] =
957:           x[0*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*(1-xir[2]) +
958:           x[3*ctx->dof+comp]*    xir[0]*(1-xir[1])*(1-xir[2]) +
959:           x[2*ctx->dof+comp]*    xir[0]*    xir[1]*(1-xir[2]) +
960:           x[1*ctx->dof+comp]*(1-xir[0])*    xir[1]*(1-xir[2]) +
961:           x[4*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*   xir[2] +
962:           x[5*ctx->dof+comp]*    xir[0]*(1-xir[1])*   xir[2] +
963:           x[6*ctx->dof+comp]*    xir[0]*    xir[1]*   xir[2] +
964:           x[7*ctx->dof+comp]*(1-xir[0])*    xir[1]*   xir[2];
965:       }
966:     } else {
967:       for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
968:     }
969:     VecRestoreArray(ref, &xi);
970:     DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
971:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
972:   }
973:   VecRestoreArray(v, &a);
974:   VecRestoreArrayRead(ctx->coords, &coords);

976:   SNESDestroy(&snes);
977:   VecDestroy(&r);
978:   VecDestroy(&ref);
979:   VecDestroy(&real);
980:   MatDestroy(&J);
981:   return 0;
982: }

984: /*@C
985:   DMInterpolationEvaluate - Using the input from dm and x, calculates interpolated field values at the interpolation points.

987:   Input Parameters:
988: + ctx - The DMInterpolationInfo context
989: . dm  - The DM
990: - x   - The local vector containing the field to be interpolated

992:   Output Parameters:
993: . v   - The vector containing the interpolated values

995:   Note: A suitable v can be obtained using DMInterpolationGetVector().

997:   Level: beginner

999: .seealso: DMInterpolationGetVector(), DMInterpolationAddPoints(), DMInterpolationCreate()
1000: @*/
1001: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
1002: {
1003:   PetscDS        ds;
1004:   PetscInt       n, p, Nf, field;
1005:   PetscBool      useDS = PETSC_FALSE;

1010:   VecGetLocalSize(v, &n);
1012:   if (!n) return 0;
1013:   DMGetDS(dm, &ds);
1014:   if (ds) {
1015:     useDS = PETSC_TRUE;
1016:     PetscDSGetNumFields(ds, &Nf);
1017:     for (field = 0; field < Nf; ++field) {
1018:       PetscObject  obj;
1019:       PetscClassId id;

1021:       PetscDSGetDiscretization(ds, field, &obj);
1022:       PetscObjectGetClassId(obj, &id);
1023:       if (id != PETSCFE_CLASSID) {useDS = PETSC_FALSE; break;}
1024:     }
1025:   }
1026:   if (useDS) {
1027:     const PetscScalar *coords;
1028:     PetscScalar       *interpolant;
1029:     PetscInt           cdim, d;

1031:     DMGetCoordinateDim(dm, &cdim);
1032:     VecGetArrayRead(ctx->coords, &coords);
1033:     VecGetArrayWrite(v, &interpolant);
1034:     for (p = 0; p < ctx->n; ++p) {
1035:       PetscReal    pcoords[3], xi[3];
1036:       PetscScalar *xa   = NULL;
1037:       PetscInt     coff = 0, foff = 0, clSize;

1039:       if (ctx->cells[p] < 0) continue;
1040:       for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p*cdim+d]);
1041:       DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi);
1042:       DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa);
1043:       for (field = 0; field < Nf; ++field) {
1044:         PetscTabulation T;
1045:         PetscFE         fe;

1047:         PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
1048:         PetscFECreateTabulation(fe, 1, 1, xi, 0, &T);
1049:         {
1050:           const PetscReal *basis = T->T[0];
1051:           const PetscInt   Nb    = T->Nb;
1052:           const PetscInt   Nc    = T->Nc;
1053:           PetscInt         f, fc;

1055:           for (fc = 0; fc < Nc; ++fc) {
1056:             interpolant[p*ctx->dof+coff+fc] = 0.0;
1057:             for (f = 0; f < Nb; ++f) {
1058:               interpolant[p*ctx->dof+coff+fc] += xa[foff+f]*basis[(0*Nb + f)*Nc + fc];
1059:             }
1060:           }
1061:           coff += Nc;
1062:           foff += Nb;
1063:         }
1064:         PetscTabulationDestroy(&T);
1065:       }
1066:       DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa);
1069:     }
1070:     VecRestoreArrayRead(ctx->coords, &coords);
1071:     VecRestoreArrayWrite(v, &interpolant);
1072:   } else {
1073:     DMPolytopeType ct;

1075:     /* TODO Check each cell individually */
1076:     DMPlexGetCellType(dm, ctx->cells[0], &ct);
1077:     switch (ct) {
1078:       case DM_POLYTOPE_SEGMENT:       DMInterpolate_Segment_Private(ctx, dm, x, v);break;
1079:       case DM_POLYTOPE_TRIANGLE:      DMInterpolate_Triangle_Private(ctx, dm, x, v);break;
1080:       case DM_POLYTOPE_QUADRILATERAL: DMInterpolate_Quad_Private(ctx, dm, x, v);break;
1081:       case DM_POLYTOPE_TETRAHEDRON:   DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);break;
1082:       case DM_POLYTOPE_HEXAHEDRON:    DMInterpolate_Hex_Private(ctx, dm, x, v);break;
1083:       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for cell type %s", DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]);
1084:     }
1085:   }
1086:   return 0;
1087: }

1089: /*@C
1090:   DMInterpolationDestroy - Destroys a DMInterpolationInfo context

1092:   Collective on ctx

1094:   Input Parameter:
1095: . ctx - the context

1097:   Level: beginner

1099: .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
1100: @*/
1101: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
1102: {
1104:   VecDestroy(&(*ctx)->coords);
1105:   PetscFree((*ctx)->points);
1106:   PetscFree((*ctx)->cells);
1107:   PetscFree(*ctx);
1108:   *ctx = NULL;
1109:   return 0;
1110: }

1112: /*@C
1113:   SNESMonitorFields - Monitors the residual for each field separately

1115:   Collective on SNES

1117:   Input Parameters:
1118: + snes   - the SNES context
1119: . its    - iteration number
1120: . fgnorm - 2-norm of residual
1121: - vf  - PetscViewerAndFormat of type ASCII

1123:   Notes:
1124:   This routine prints the residual norm at each iteration.

1126:   Level: intermediate

1128: .seealso: SNESMonitorSet(), SNESMonitorDefault()
1129: @*/
1130: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
1131: {
1132:   PetscViewer        viewer = vf->viewer;
1133:   Vec                res;
1134:   DM                 dm;
1135:   PetscSection       s;
1136:   const PetscScalar *r;
1137:   PetscReal         *lnorms, *norms;
1138:   PetscInt           numFields, f, pStart, pEnd, p;

1141:   SNESGetFunction(snes, &res, NULL, NULL);
1142:   SNESGetDM(snes, &dm);
1143:   DMGetLocalSection(dm, &s);
1144:   PetscSectionGetNumFields(s, &numFields);
1145:   PetscSectionGetChart(s, &pStart, &pEnd);
1146:   PetscCalloc2(numFields, &lnorms, numFields, &norms);
1147:   VecGetArrayRead(res, &r);
1148:   for (p = pStart; p < pEnd; ++p) {
1149:     for (f = 0; f < numFields; ++f) {
1150:       PetscInt fdof, foff, d;

1152:       PetscSectionGetFieldDof(s, p, f, &fdof);
1153:       PetscSectionGetFieldOffset(s, p, f, &foff);
1154:       for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff+d]));
1155:     }
1156:   }
1157:   VecRestoreArrayRead(res, &r);
1158:   MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1159:   PetscViewerPushFormat(viewer,vf->format);
1160:   PetscViewerASCIIAddTab(viewer, ((PetscObject) snes)->tablevel);
1161:   PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e [", its, (double) fgnorm);
1162:   for (f = 0; f < numFields; ++f) {
1163:     if (f > 0) PetscViewerASCIIPrintf(viewer, ", ");
1164:     PetscViewerASCIIPrintf(viewer, "%14.12e", (double) PetscSqrtReal(norms[f]));
1165:   }
1166:   PetscViewerASCIIPrintf(viewer, "]\n");
1167:   PetscViewerASCIISubtractTab(viewer, ((PetscObject) snes)->tablevel);
1168:   PetscViewerPopFormat(viewer);
1169:   PetscFree2(lnorms, norms);
1170:   return 0;
1171: }

1173: /********************* Residual Computation **************************/

1175: PetscErrorCode DMPlexGetAllCells_Internal(DM plex, IS *cellIS)
1176: {
1177:   PetscInt       depth;

1179:   DMPlexGetDepth(plex, &depth);
1180:   DMGetStratumIS(plex, "dim", depth, cellIS);
1181:   if (!*cellIS) DMGetStratumIS(plex, "depth", depth, cellIS);
1182:   return 0;
1183: }

1185: /*@
1186:   DMPlexSNESComputeResidualFEM - Sums the local residual into vector F from the local input X using pointwise functions specified by the user

1188:   Input Parameters:
1189: + dm - The mesh
1190: . X  - Local solution
1191: - user - The user context

1193:   Output Parameter:
1194: . F  - Local output vector

1196:   Notes:
1197:   The residual is summed into F; the caller is responsible for using VecZeroEntries() or otherwise ensuring that any data in F is intentional.

1199:   Level: developer

1201: .seealso: DMPlexComputeJacobianAction()
1202: @*/
1203: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1204: {
1205:   DM             plex;
1206:   IS             allcellIS;
1207:   PetscInt       Nds, s;

1209:   DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1210:   DMPlexGetAllCells_Internal(plex, &allcellIS);
1211:   DMGetNumDS(dm, &Nds);
1212:   for (s = 0; s < Nds; ++s) {
1213:     PetscDS          ds;
1214:     IS               cellIS;
1215:     PetscFormKey key;

1217:     DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
1218:     key.value = 0;
1219:     key.field = 0;
1220:     key.part  = 0;
1221:     if (!key.label) {
1222:       PetscObjectReference((PetscObject) allcellIS);
1223:       cellIS = allcellIS;
1224:     } else {
1225:       IS pointIS;

1227:       key.value = 1;
1228:       DMLabelGetStratumIS(key.label, key.value, &pointIS);
1229:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1230:       ISDestroy(&pointIS);
1231:     }
1232:     DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user);
1233:     ISDestroy(&cellIS);
1234:   }
1235:   ISDestroy(&allcellIS);
1236:   DMDestroy(&plex);
1237:   return 0;
1238: }

1240: PetscErrorCode DMSNESComputeResidual(DM dm, Vec X, Vec F, void *user)
1241: {
1242:   DM             plex;
1243:   IS             allcellIS;
1244:   PetscInt       Nds, s;

1246:   DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1247:   DMPlexGetAllCells_Internal(plex, &allcellIS);
1248:   DMGetNumDS(dm, &Nds);
1249:   for (s = 0; s < Nds; ++s) {
1250:     PetscDS ds;
1251:     DMLabel label;
1252:     IS      cellIS;

1254:     DMGetRegionNumDS(dm, s, &label, NULL, &ds);
1255:     {
1256:       PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1};
1257:       PetscWeakForm     wf;
1258:       PetscInt          Nm = 2, m, Nk = 0, k, kp, off = 0;
1259:       PetscFormKey *reskeys;

1261:       /* Get unique residual keys */
1262:       for (m = 0; m < Nm; ++m) {
1263:         PetscInt Nkm;
1264:         PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm);
1265:         Nk  += Nkm;
1266:       }
1267:       PetscMalloc1(Nk, &reskeys);
1268:       for (m = 0; m < Nm; ++m) {
1269:         PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys);
1270:       }
1272:       PetscFormKeySort(Nk, reskeys);
1273:       for (k = 0, kp = 1; kp < Nk; ++kp) {
1274:         if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) {
1275:           ++k;
1276:           if (kp != k) reskeys[k] = reskeys[kp];
1277:         }
1278:       }
1279:       Nk = k;

1281:       PetscDSGetWeakForm(ds, &wf);
1282:       for (k = 0; k < Nk; ++k) {
1283:         DMLabel  label = reskeys[k].label;
1284:         PetscInt val   = reskeys[k].value;

1286:         if (!label) {
1287:           PetscObjectReference((PetscObject) allcellIS);
1288:           cellIS = allcellIS;
1289:         } else {
1290:           IS pointIS;

1292:           DMLabelGetStratumIS(label, val, &pointIS);
1293:           ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1294:           ISDestroy(&pointIS);
1295:         }
1296:         DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user);
1297:         ISDestroy(&cellIS);
1298:       }
1299:       PetscFree(reskeys);
1300:     }
1301:   }
1302:   ISDestroy(&allcellIS);
1303:   DMDestroy(&plex);
1304:   return 0;
1305: }

1307: /*@
1308:   DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input X

1310:   Input Parameters:
1311: + dm - The mesh
1312: - user - The user context

1314:   Output Parameter:
1315: . X  - Local solution

1317:   Level: developer

1319: .seealso: DMPlexComputeJacobianAction()
1320: @*/
1321: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
1322: {
1323:   DM             plex;

1325:   DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
1326:   DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL);
1327:   DMDestroy(&plex);
1328:   return 0;
1329: }

1331: /*@
1332:   DMSNESComputeJacobianAction - Compute the action of the Jacobian J(X) on Y

1334:   Input Parameters:
1335: + dm   - The DM
1336: . X    - Local solution vector
1337: . Y    - Local input vector
1338: - user - The user context

1340:   Output Parameter:
1341: . F    - lcoal output vector

1343:   Level: developer

1345:   Notes:
1346:   Users will typically use DMSNESCreateJacobianMF() followed by MatMult() instead of calling this routine directly.

1348: .seealso: DMSNESCreateJacobianMF(), DMPlexSNESComputeResidualFEM()
1349: @*/
1350: PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, void *user)
1351: {
1352:   DM             plex;
1353:   IS             allcellIS;
1354:   PetscInt       Nds, s;

1356:   DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1357:   DMPlexGetAllCells_Internal(plex, &allcellIS);
1358:   DMGetNumDS(dm, &Nds);
1359:   for (s = 0; s < Nds; ++s) {
1360:     PetscDS ds;
1361:     DMLabel label;
1362:     IS      cellIS;

1364:     DMGetRegionNumDS(dm, s, &label, NULL, &ds);
1365:     {
1366:       PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3};
1367:       PetscWeakForm     wf;
1368:       PetscInt          Nm = 4, m, Nk = 0, k, kp, off = 0;
1369:       PetscFormKey *jackeys;

1371:       /* Get unique Jacobian keys */
1372:       for (m = 0; m < Nm; ++m) {
1373:         PetscInt Nkm;
1374:         PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm);
1375:         Nk  += Nkm;
1376:       }
1377:       PetscMalloc1(Nk, &jackeys);
1378:       for (m = 0; m < Nm; ++m) {
1379:         PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys);
1380:       }
1382:       PetscFormKeySort(Nk, jackeys);
1383:       for (k = 0, kp = 1; kp < Nk; ++kp) {
1384:         if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) {
1385:           ++k;
1386:           if (kp != k) jackeys[k] = jackeys[kp];
1387:         }
1388:       }
1389:       Nk = k;

1391:       PetscDSGetWeakForm(ds, &wf);
1392:       for (k = 0; k < Nk; ++k) {
1393:         DMLabel  label = jackeys[k].label;
1394:         PetscInt val   = jackeys[k].value;

1396:         if (!label) {
1397:           PetscObjectReference((PetscObject) allcellIS);
1398:           cellIS = allcellIS;
1399:         } else {
1400:           IS pointIS;

1402:           DMLabelGetStratumIS(label, val, &pointIS);
1403:           ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1404:           ISDestroy(&pointIS);
1405:         }
1406:         DMPlexComputeJacobian_Action_Internal(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, user);
1407:         ISDestroy(&cellIS);
1408:       }
1409:       PetscFree(jackeys);
1410:     }
1411:   }
1412:   ISDestroy(&allcellIS);
1413:   DMDestroy(&plex);
1414:   return 0;
1415: }

1417: /*@
1418:   DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.

1420:   Input Parameters:
1421: + dm - The mesh
1422: . X  - Local input vector
1423: - user - The user context

1425:   Output Parameter:
1426: . Jac  - Jacobian matrix

1428:   Note:
1429:   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1430:   like a GPU, or vectorize on a multicore machine.

1432:   Level: developer

1434: .seealso: FormFunctionLocal()
1435: @*/
1436: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1437: {
1438:   DM             plex;
1439:   IS             allcellIS;
1440:   PetscBool      hasJac, hasPrec;
1441:   PetscInt       Nds, s;

1443:   DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1444:   DMPlexGetAllCells_Internal(plex, &allcellIS);
1445:   DMGetNumDS(dm, &Nds);
1446:   for (s = 0; s < Nds; ++s) {
1447:     PetscDS          ds;
1448:     IS               cellIS;
1449:     PetscFormKey key;

1451:     DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
1452:     key.value = 0;
1453:     key.field = 0;
1454:     key.part  = 0;
1455:     if (!key.label) {
1456:       PetscObjectReference((PetscObject) allcellIS);
1457:       cellIS = allcellIS;
1458:     } else {
1459:       IS pointIS;

1461:       key.value = 1;
1462:       DMLabelGetStratumIS(key.label, key.value, &pointIS);
1463:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1464:       ISDestroy(&pointIS);
1465:     }
1466:     if (!s) {
1467:       PetscDSHasJacobian(ds, &hasJac);
1468:       PetscDSHasJacobianPreconditioner(ds, &hasPrec);
1469:       if (hasJac && hasPrec) MatZeroEntries(Jac);
1470:       MatZeroEntries(JacP);
1471:     }
1472:     DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user);
1473:     ISDestroy(&cellIS);
1474:   }
1475:   ISDestroy(&allcellIS);
1476:   DMDestroy(&plex);
1477:   return 0;
1478: }

1480: struct _DMSNESJacobianMFCtx
1481: {
1482:   DM    dm;
1483:   Vec   X;
1484:   void *ctx;
1485: };

1487: static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A)
1488: {
1489:   struct _DMSNESJacobianMFCtx *ctx;

1491:   MatShellGetContext(A, &ctx);
1492:   MatShellSetContext(A, NULL);
1493:   DMDestroy(&ctx->dm);
1494:   VecDestroy(&ctx->X);
1495:   PetscFree(ctx);
1496:   return 0;
1497: }

1499: static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z)
1500: {
1501:   struct _DMSNESJacobianMFCtx *ctx;

1503:   MatShellGetContext(A, &ctx);
1504:   DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx);
1505:   return 0;
1506: }

1508: /*@
1509:   DMSNESCreateJacobianMF - Create a Mat which computes the action of the Jacobian matrix-free

1511:   Collective on dm

1513:   Input Parameters:
1514: + dm   - The DM
1515: . X    - The evaluation point for the Jacobian
1516: - user - A user context, or NULL

1518:   Output Parameter:
1519: . J    - The Mat

1521:   Level: advanced

1523:   Notes:
1524:   Vec X is kept in Mat J, so updating X then updates the evaluation point.

1526: .seealso: DMSNESComputeJacobianAction()
1527: @*/
1528: PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, void *user, Mat *J)
1529: {
1530:   struct _DMSNESJacobianMFCtx *ctx;
1531:   PetscInt                     n, N;

1533:   MatCreate(PetscObjectComm((PetscObject) dm), J);
1534:   MatSetType(*J, MATSHELL);
1535:   VecGetLocalSize(X, &n);
1536:   VecGetSize(X, &N);
1537:   MatSetSizes(*J, n, n, N, N);
1538:   PetscObjectReference((PetscObject) dm);
1539:   PetscObjectReference((PetscObject) X);
1540:   PetscMalloc1(1, &ctx);
1541:   ctx->dm  = dm;
1542:   ctx->X   = X;
1543:   ctx->ctx = user;
1544:   MatShellSetContext(*J, ctx);
1545:   MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void)) DMSNESJacobianMF_Destroy_Private);
1546:   MatShellSetOperation(*J, MATOP_MULT,    (void (*)(void)) DMSNESJacobianMF_Mult_Private);
1547:   return 0;
1548: }

1550: /*
1551:      MatComputeNeumannOverlap - Computes an unassembled (Neumann) local overlapping Mat in nonlinear context.

1553:    Input Parameters:
1554: +     X - SNES linearization point
1555: .     ovl - index set of overlapping subdomains

1557:    Output Parameter:
1558: .     J - unassembled (Neumann) local matrix

1560:    Level: intermediate

1562: .seealso: DMCreateNeumannOverlap(), MATIS, PCHPDDMSetAuxiliaryMat()
1563: */
1564: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx)
1565: {
1566:   SNES           snes;
1567:   Mat            pJ;
1568:   DM             ovldm,origdm;
1569:   DMSNES         sdm;
1570:   PetscErrorCode (*bfun)(DM,Vec,void*);
1571:   PetscErrorCode (*jfun)(DM,Vec,Mat,Mat,void*);
1572:   void           *bctx,*jctx;

1574:   PetscObjectQuery((PetscObject)ovl,"_DM_Overlap_HPDDM_MATIS",(PetscObject*)&pJ);
1576:   PetscObjectQuery((PetscObject)ovl,"_DM_Original_HPDDM",(PetscObject*)&origdm);
1578:   MatGetDM(pJ,&ovldm);
1579:   DMSNESGetBoundaryLocal(origdm,&bfun,&bctx);
1580:   DMSNESSetBoundaryLocal(ovldm,bfun,bctx);
1581:   DMSNESGetJacobianLocal(origdm,&jfun,&jctx);
1582:   DMSNESSetJacobianLocal(ovldm,jfun,jctx);
1583:   PetscObjectQuery((PetscObject)ovl,"_DM_Overlap_HPDDM_SNES",(PetscObject*)&snes);
1584:   if (!snes) {
1585:     SNESCreate(PetscObjectComm((PetscObject)ovl),&snes);
1586:     SNESSetDM(snes,ovldm);
1587:     PetscObjectCompose((PetscObject)ovl,"_DM_Overlap_HPDDM_SNES",(PetscObject)snes);
1588:     PetscObjectDereference((PetscObject)snes);
1589:   }
1590:   DMGetDMSNES(ovldm,&sdm);
1591:   VecLockReadPush(X);
1592:   PetscStackPush("SNES user Jacobian function");
1593:   (*sdm->ops->computejacobian)(snes,X,pJ,pJ,sdm->jacobianctx);
1594:   PetscStackPop;
1595:   VecLockReadPop(X);
1596:   /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
1597:   {
1598:     Mat locpJ;

1600:     MatISGetLocalMat(pJ,&locpJ);
1601:     MatCopy(locpJ,J,SAME_NONZERO_PATTERN);
1602:   }
1603:   return 0;
1604: }

1606: /*@
1607:   DMPlexSetSNESLocalFEM - Use DMPlex's internal FEM routines to compute SNES boundary values, residual, and Jacobian.

1609:   Input Parameters:
1610: + dm - The DM object
1611: . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see PetscDSAddBoundary())
1612: . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see PetscDSSetResidual())
1613: - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see PetscDSSetJacobian())

1615:   Level: developer
1616: @*/
1617: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
1618: {
1619:   DMSNESSetBoundaryLocal(dm,DMPlexSNESComputeBoundaryFEM,boundaryctx);
1620:   DMSNESSetFunctionLocal(dm,DMPlexSNESComputeResidualFEM,residualctx);
1621:   DMSNESSetJacobianLocal(dm,DMPlexSNESComputeJacobianFEM,jacobianctx);
1622:   PetscObjectComposeFunction((PetscObject)dm,"MatComputeNeumannOverlap_C",MatComputeNeumannOverlap_Plex);
1623:   return 0;
1624: }

1626: /*@C
1627:   DMSNESCheckDiscretization - Check the discretization error of the exact solution

1629:   Input Parameters:
1630: + snes - the SNES object
1631: . dm   - the DM
1632: . t    - the time
1633: . u    - a DM vector
1634: - tol  - A tolerance for the check, or -1 to print the results instead

1636:   Output Parameters:
1637: . error - An array which holds the discretization error in each field, or NULL

1639:   Note: The user must call PetscDSSetExactSolution() beforehand

1641:   Level: developer

1643: .seealso: DNSNESCheckFromOptions(), DMSNESCheckResidual(), DMSNESCheckJacobian(), PetscDSSetExactSolution()
1644: @*/
1645: PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
1646: {
1647:   PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
1648:   void            **ectxs;
1649:   PetscReal        *err;
1650:   MPI_Comm          comm;
1651:   PetscInt          Nf, f;


1658:   DMComputeExactSolution(dm, t, u, NULL);
1659:   VecViewFromOptions(u, NULL, "-vec_view");

1661:   PetscObjectGetComm((PetscObject) snes, &comm);
1662:   DMGetNumFields(dm, &Nf);
1663:   PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err);
1664:   {
1665:     PetscInt Nds, s;

1667:     DMGetNumDS(dm, &Nds);
1668:     for (s = 0; s < Nds; ++s) {
1669:       PetscDS         ds;
1670:       DMLabel         label;
1671:       IS              fieldIS;
1672:       const PetscInt *fields;
1673:       PetscInt        dsNf, f;

1675:       DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds);
1676:       PetscDSGetNumFields(ds, &dsNf);
1677:       ISGetIndices(fieldIS, &fields);
1678:       for (f = 0; f < dsNf; ++f) {
1679:         const PetscInt field = fields[f];
1680:         PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]);
1681:       }
1682:       ISRestoreIndices(fieldIS, &fields);
1683:     }
1684:   }
1685:   if (Nf > 1) {
1686:     DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err);
1687:     if (tol >= 0.0) {
1688:       for (f = 0; f < Nf; ++f) {
1690:       }
1691:     } else if (error) {
1692:       for (f = 0; f < Nf; ++f) error[f] = err[f];
1693:     } else {
1694:       PetscPrintf(comm, "L_2 Error: [");
1695:       for (f = 0; f < Nf; ++f) {
1696:         if (f) PetscPrintf(comm, ", ");
1697:         PetscPrintf(comm, "%g", (double)err[f]);
1698:       }
1699:       PetscPrintf(comm, "]\n");
1700:     }
1701:   } else {
1702:     DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]);
1703:     if (tol >= 0.0) {
1705:     } else if (error) {
1706:       error[0] = err[0];
1707:     } else {
1708:       PetscPrintf(comm, "L_2 Error: %g\n", (double) err[0]);
1709:     }
1710:   }
1711:   PetscFree3(exacts, ectxs, err);
1712:   return 0;
1713: }

1715: /*@C
1716:   DMSNESCheckResidual - Check the residual of the exact solution

1718:   Input Parameters:
1719: + snes - the SNES object
1720: . dm   - the DM
1721: . u    - a DM vector
1722: - tol  - A tolerance for the check, or -1 to print the results instead

1724:   Output Parameters:
1725: . residual - The residual norm of the exact solution, or NULL

1727:   Level: developer

1729: .seealso: DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckJacobian()
1730: @*/
1731: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
1732: {
1733:   MPI_Comm       comm;
1734:   Vec            r;
1735:   PetscReal      res;

1741:   PetscObjectGetComm((PetscObject) snes, &comm);
1742:   DMComputeExactSolution(dm, 0.0, u, NULL);
1743:   VecDuplicate(u, &r);
1744:   SNESComputeFunction(snes, u, r);
1745:   VecNorm(r, NORM_2, &res);
1746:   if (tol >= 0.0) {
1748:   } else if (residual) {
1749:     *residual = res;
1750:   } else {
1751:     PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);
1752:     VecChop(r, 1.0e-10);
1753:     PetscObjectSetName((PetscObject) r, "Initial Residual");
1754:     PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
1755:     VecViewFromOptions(r, NULL, "-vec_view");
1756:   }
1757:   VecDestroy(&r);
1758:   return 0;
1759: }

1761: /*@C
1762:   DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test

1764:   Input Parameters:
1765: + snes - the SNES object
1766: . dm   - the DM
1767: . u    - a DM vector
1768: - tol  - A tolerance for the check, or -1 to print the results instead

1770:   Output Parameters:
1771: + isLinear - Flag indicaing that the function looks linear, or NULL
1772: - convRate - The rate of convergence of the linear model, or NULL

1774:   Level: developer

1776: .seealso: DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckResidual()
1777: @*/
1778: PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
1779: {
1780:   MPI_Comm       comm;
1781:   PetscDS        ds;
1782:   Mat            J, M;
1783:   MatNullSpace   nullspace;
1784:   PetscReal      slope, intercept;
1785:   PetscBool      hasJac, hasPrec, isLin = PETSC_FALSE;

1792:   PetscObjectGetComm((PetscObject) snes, &comm);
1793:   DMComputeExactSolution(dm, 0.0, u, NULL);
1794:   /* Create and view matrices */
1795:   DMCreateMatrix(dm, &J);
1796:   DMGetDS(dm, &ds);
1797:   PetscDSHasJacobian(ds, &hasJac);
1798:   PetscDSHasJacobianPreconditioner(ds, &hasPrec);
1799:   if (hasJac && hasPrec) {
1800:     DMCreateMatrix(dm, &M);
1801:     SNESComputeJacobian(snes, u, J, M);
1802:     PetscObjectSetName((PetscObject) M, "Preconditioning Matrix");
1803:     PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");
1804:     MatViewFromOptions(M, NULL, "-mat_view");
1805:     MatDestroy(&M);
1806:   } else {
1807:     SNESComputeJacobian(snes, u, J, J);
1808:   }
1809:   PetscObjectSetName((PetscObject) J, "Jacobian");
1810:   PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");
1811:   MatViewFromOptions(J, NULL, "-mat_view");
1812:   /* Check nullspace */
1813:   MatGetNullSpace(J, &nullspace);
1814:   if (nullspace) {
1815:     PetscBool isNull;
1816:     MatNullSpaceTest(nullspace, J, &isNull);
1818:   }
1819:   /* Taylor test */
1820:   {
1821:     PetscRandom rand;
1822:     Vec         du, uhat, r, rhat, df;
1823:     PetscReal   h;
1824:     PetscReal  *es, *hs, *errors;
1825:     PetscReal   hMax = 1.0, hMin = 1e-6, hMult = 0.1;
1826:     PetscInt    Nv, v;

1828:     /* Choose a perturbation direction */
1829:     PetscRandomCreate(comm, &rand);
1830:     VecDuplicate(u, &du);
1831:     VecSetRandom(du, rand);
1832:     PetscRandomDestroy(&rand);
1833:     VecDuplicate(u, &df);
1834:     MatMult(J, du, df);
1835:     /* Evaluate residual at u, F(u), save in vector r */
1836:     VecDuplicate(u, &r);
1837:     SNESComputeFunction(snes, u, r);
1838:     /* Look at the convergence of our Taylor approximation as we approach u */
1839:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
1840:     PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
1841:     VecDuplicate(u, &uhat);
1842:     VecDuplicate(u, &rhat);
1843:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1844:       VecWAXPY(uhat, h, du, u);
1845:       /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1846:       SNESComputeFunction(snes, uhat, rhat);
1847:       VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
1848:       VecNorm(rhat, NORM_2, &errors[Nv]);

1850:       es[Nv] = PetscLog10Real(errors[Nv]);
1851:       hs[Nv] = PetscLog10Real(h);
1852:     }
1853:     VecDestroy(&uhat);
1854:     VecDestroy(&rhat);
1855:     VecDestroy(&df);
1856:     VecDestroy(&r);
1857:     VecDestroy(&du);
1858:     for (v = 0; v < Nv; ++v) {
1859:       if ((tol >= 0) && (errors[v] > tol)) break;
1860:       else if (errors[v] > PETSC_SMALL)    break;
1861:     }
1862:     if (v == Nv) isLin = PETSC_TRUE;
1863:     PetscLinearRegression(Nv, hs, es, &slope, &intercept);
1864:     PetscFree3(es, hs, errors);
1865:     /* Slope should be about 2 */
1866:     if (tol >= 0) {
1868:     } else if (isLinear || convRate) {
1869:       if (isLinear) *isLinear = isLin;
1870:       if (convRate) *convRate = slope;
1871:     } else {
1872:       if (!isLin) PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double) slope);
1873:       else        PetscPrintf(comm, "Function appears to be linear\n");
1874:     }
1875:   }
1876:   MatDestroy(&J);
1877:   return 0;
1878: }

1880: PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
1881: {
1882:   DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL);
1883:   DMSNESCheckResidual(snes, dm, u, -1.0, NULL);
1884:   DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL);
1885:   return 0;
1886: }

1888: /*@C
1889:   DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information

1891:   Input Parameters:
1892: + snes - the SNES object
1893: - u    - representative SNES vector

1895:   Note: The user must call PetscDSSetExactSolution() beforehand

1897:   Level: developer
1898: @*/
1899: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1900: {
1901:   DM             dm;
1902:   Vec            sol;
1903:   PetscBool      check;

1905:   PetscOptionsHasName(((PetscObject)snes)->options,((PetscObject)snes)->prefix, "-dmsnes_check", &check);
1906:   if (!check) return 0;
1907:   SNESGetDM(snes, &dm);
1908:   VecDuplicate(u, &sol);
1909:   SNESSetSolution(snes, sol);
1910:   DMSNESCheck_Internal(snes, dm, sol);
1911:   VecDestroy(&sol);
1912:   return 0;
1913: }