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