Actual source code: ex7.c
1: static char help[] = "Fermions on a hypercubic lattice.\n\n";
3: #include <petscdmplex.h>
4: #include <petscsnes.h>
6: /* Common operations:
8: - View the input \psi as ASCII in lexicographic order: -psi_view
9: */
11: const PetscReal M = 1.0;
13: typedef struct {
14: PetscBool usePV; /* Use Pauli-Villars preconditioning */
15: } AppCtx;
17: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
18: {
19: PetscFunctionBegin;
20: options->usePV = PETSC_TRUE;
22: PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
23: PetscCall(PetscOptionsBool("-use_pv", "Use Pauli-Villars preconditioning", "ex1.c", options->usePV, &options->usePV, NULL));
24: PetscOptionsEnd();
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
29: {
30: PetscFunctionBegin;
31: PetscCall(DMCreate(comm, dm));
32: PetscCall(DMSetType(*dm, DMPLEX));
33: PetscCall(DMSetFromOptions(*dm));
34: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
39: {
40: PetscSection s;
41: PetscInt vStart, vEnd, v;
43: PetscFunctionBegin;
44: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
45: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
46: PetscCall(PetscSectionSetChart(s, vStart, vEnd));
47: for (v = vStart; v < vEnd; ++v) {
48: PetscCall(PetscSectionSetDof(s, v, 12));
49: /* TODO Divide the values into fields/components */
50: }
51: PetscCall(PetscSectionSetUp(s));
52: PetscCall(DMSetLocalSection(dm, s));
53: PetscCall(PetscSectionDestroy(&s));
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: static PetscErrorCode SetupAuxDiscretization(DM dm, AppCtx *user)
58: {
59: DM dmAux, coordDM;
60: PetscSection s;
61: Vec gauge;
62: PetscInt eStart, eEnd, e;
64: PetscFunctionBegin;
65: /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
66: PetscCall(DMGetCoordinateDM(dm, &coordDM));
67: PetscCall(DMClone(dm, &dmAux));
68: PetscCall(DMSetCoordinateDM(dmAux, coordDM));
69: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
70: PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
71: PetscCall(PetscSectionSetChart(s, eStart, eEnd));
72: for (e = eStart; e < eEnd; ++e) {
73: /* TODO Should we store the whole SU(3) matrix, or the symmetric part? */
74: PetscCall(PetscSectionSetDof(s, e, 9));
75: }
76: PetscCall(PetscSectionSetUp(s));
77: PetscCall(DMSetLocalSection(dmAux, s));
78: PetscCall(PetscSectionDestroy(&s));
79: PetscCall(DMCreateLocalVector(dmAux, &gauge));
80: PetscCall(DMDestroy(&dmAux));
81: PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, gauge));
82: PetscCall(VecDestroy(&gauge));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: static PetscErrorCode PrintVertex(DM dm, PetscInt v)
87: {
88: MPI_Comm comm;
89: PetscContainer c;
90: PetscInt *extent;
91: PetscInt dim, cStart, cEnd, sum;
93: PetscFunctionBeginUser;
94: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
95: PetscCall(DMGetDimension(dm, &dim));
96: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
97: PetscCall(PetscObjectQuery((PetscObject)dm, "_extent", (PetscObject *)&c));
98: PetscCall(PetscContainerGetPointer(c, (void **)&extent));
99: sum = 1;
100: PetscCall(PetscPrintf(comm, "Vertex %" PetscInt_FMT ":", v));
101: for (PetscInt d = 0; d < dim; ++d) {
102: PetscCall(PetscPrintf(comm, " %" PetscInt_FMT, (v / sum) % extent[d]));
103: if (d < dim) sum *= extent[d];
104: }
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: // Apply \gamma_\mu
109: static PetscErrorCode ComputeGamma(PetscInt d, PetscInt ldx, PetscScalar f[])
110: {
111: const PetscScalar fin[4] = {f[0 * ldx], f[1 * ldx], f[2 * ldx], f[3 * ldx]};
113: PetscFunctionBeginHot;
114: switch (d) {
115: case 0:
116: f[0 * ldx] = PETSC_i * fin[3];
117: f[1 * ldx] = PETSC_i * fin[2];
118: f[2 * ldx] = -PETSC_i * fin[1];
119: f[3 * ldx] = -PETSC_i * fin[0];
120: break;
121: case 1:
122: f[0 * ldx] = -fin[3];
123: f[1 * ldx] = fin[2];
124: f[2 * ldx] = fin[1];
125: f[3 * ldx] = -fin[0];
126: break;
127: case 2:
128: f[0 * ldx] = PETSC_i * fin[2];
129: f[1 * ldx] = -PETSC_i * fin[3];
130: f[2 * ldx] = -PETSC_i * fin[0];
131: f[3 * ldx] = PETSC_i * fin[1];
132: break;
133: case 3:
134: f[0 * ldx] = fin[2];
135: f[1 * ldx] = fin[3];
136: f[2 * ldx] = fin[0];
137: f[3 * ldx] = fin[1];
138: break;
139: default:
140: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Direction for gamma %" PetscInt_FMT " not in [0, 4)", d);
141: }
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: // Apply (1 \pm \gamma_\mu)/2
146: static PetscErrorCode ComputeGammaFactor(PetscInt d, PetscBool forward, PetscInt ldx, PetscScalar f[])
147: {
148: const PetscReal sign = forward ? -1. : 1.;
149: const PetscScalar fin[4] = {f[0 * ldx], f[1 * ldx], f[2 * ldx], f[3 * ldx]};
151: PetscFunctionBeginHot;
152: switch (d) {
153: case 0:
154: f[0 * ldx] += sign * PETSC_i * fin[3];
155: f[1 * ldx] += sign * PETSC_i * fin[2];
156: f[2 * ldx] += sign * -PETSC_i * fin[1];
157: f[3 * ldx] += sign * -PETSC_i * fin[0];
158: break;
159: case 1:
160: f[0 * ldx] += -sign * fin[3];
161: f[1 * ldx] += sign * fin[2];
162: f[2 * ldx] += sign * fin[1];
163: f[3 * ldx] += -sign * fin[0];
164: break;
165: case 2:
166: f[0 * ldx] += sign * PETSC_i * fin[2];
167: f[1 * ldx] += sign * -PETSC_i * fin[3];
168: f[2 * ldx] += sign * -PETSC_i * fin[0];
169: f[3 * ldx] += sign * PETSC_i * fin[1];
170: break;
171: case 3:
172: f[0 * ldx] += sign * fin[2];
173: f[1 * ldx] += sign * fin[3];
174: f[2 * ldx] += sign * fin[0];
175: f[3 * ldx] += sign * fin[1];
176: break;
177: default:
178: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Direction for gamma %" PetscInt_FMT " not in [0, 4)", d);
179: }
180: f[0 * ldx] *= 0.5;
181: f[1 * ldx] *= 0.5;
182: f[2 * ldx] *= 0.5;
183: f[3 * ldx] *= 0.5;
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: #include <petsc/private/dmpleximpl.h>
189: // ComputeAction() sums the action of 1/2 (1 \pm \gamma_\mu) U \psi into f
190: static PetscErrorCode ComputeAction(PetscInt d, PetscBool forward, const PetscScalar U[], const PetscScalar psi[], PetscScalar f[])
191: {
192: PetscScalar tmp[12];
194: PetscFunctionBeginHot;
195: // Apply U
196: for (PetscInt beta = 0; beta < 4; ++beta) {
197: if (forward) DMPlex_Mult3D_Internal(U, 1, &psi[beta * 3], &tmp[beta * 3]);
198: else DMPlex_MultTranspose3D_Internal(U, 1, &psi[beta * 3], &tmp[beta * 3]);
199: }
200: // Apply (1 \pm \gamma_\mu)/2 to each color
201: for (PetscInt c = 0; c < 3; ++c) PetscCall(ComputeGammaFactor(d, forward, 3, &tmp[c]));
202: // Note that we are subtracting this contribution
203: for (PetscInt i = 0; i < 12; ++i) f[i] -= tmp[i];
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: /*
208: The assembly loop runs over vertices. Each vertex has 2d edges in its support. The edges are ordered first by the dimension along which they run, and second from smaller to larger index, expect for edges which loop periodically. The vertices on edges are also ordered from smaller to larger index except for periodic edges.
209: */
210: static PetscErrorCode ComputeResidual(DM dm, Vec u, Vec f)
211: {
212: DM dmAux;
213: Vec gauge;
214: PetscSection s, sGauge;
215: const PetscScalar *ua;
216: PetscScalar *fa, *link;
217: PetscInt dim, vStart, vEnd;
219: PetscFunctionBeginUser;
220: PetscCall(DMGetDimension(dm, &dim));
221: PetscCall(DMGetLocalSection(dm, &s));
222: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
223: PetscCall(VecGetArrayRead(u, &ua));
224: PetscCall(VecGetArray(f, &fa));
226: PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &gauge));
227: PetscCall(VecGetDM(gauge, &dmAux));
228: PetscCall(DMGetLocalSection(dmAux, &sGauge));
229: PetscCall(VecGetArray(gauge, &link));
230: // Loop over y
231: for (PetscInt v = vStart; v < vEnd; ++v) {
232: const PetscInt *supp;
233: PetscInt xdof, xoff;
235: PetscCall(DMPlexGetSupport(dm, v, &supp));
236: PetscCall(PetscSectionGetDof(s, v, &xdof));
237: PetscCall(PetscSectionGetOffset(s, v, &xoff));
238: // Diagonal
239: for (PetscInt i = 0; i < xdof; ++i) fa[xoff + i] += (M + 4) * ua[xoff + i];
240: // Loop over mu
241: for (PetscInt d = 0; d < dim; ++d) {
242: const PetscInt *cone;
243: PetscInt yoff, goff;
245: // Left action -(1 + \gamma_\mu)/2 \otimes U^\dagger_\mu(y) \delta_{x - \mu,y} \psi(y)
246: PetscCall(DMPlexGetCone(dm, supp[2 * d + 0], &cone));
247: PetscCall(PetscSectionGetOffset(s, cone[0], &yoff));
248: PetscCall(PetscSectionGetOffset(sGauge, supp[2 * d + 0], &goff));
249: PetscCall(ComputeAction(d, PETSC_FALSE, &link[goff], &ua[yoff], &fa[xoff]));
250: // Right edge -(1 - \gamma_\mu)/2 \otimes U_\mu(x) \delta_{x + \mu,y} \psi(y)
251: PetscCall(DMPlexGetCone(dm, supp[2 * d + 1], &cone));
252: PetscCall(PetscSectionGetOffset(s, cone[1], &yoff));
253: PetscCall(PetscSectionGetOffset(sGauge, supp[2 * d + 1], &goff));
254: PetscCall(ComputeAction(d, PETSC_TRUE, &link[goff], &ua[yoff], &fa[xoff]));
255: }
256: }
257: PetscCall(VecRestoreArray(f, &fa));
258: PetscCall(VecRestoreArray(gauge, &link));
259: PetscCall(VecRestoreArrayRead(u, &ua));
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: static PetscErrorCode CreateJacobian(DM dm, Mat *J)
264: {
265: PetscFunctionBegin;
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: static PetscErrorCode ComputeJacobian(DM dm, Vec u, Mat J)
270: {
271: PetscFunctionBegin;
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: static PetscErrorCode PrintTraversal(DM dm)
276: {
277: MPI_Comm comm;
278: PetscInt vStart, vEnd;
280: PetscFunctionBeginUser;
281: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
282: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
283: for (PetscInt v = vStart; v < vEnd; ++v) {
284: const PetscInt *supp;
285: PetscInt Ns;
287: PetscCall(DMPlexGetSupportSize(dm, v, &Ns));
288: PetscCall(DMPlexGetSupport(dm, v, &supp));
289: PetscCall(PrintVertex(dm, v));
290: PetscCall(PetscPrintf(comm, "\n"));
291: for (PetscInt s = 0; s < Ns; ++s) {
292: const PetscInt *cone;
294: PetscCall(DMPlexGetCone(dm, supp[s], &cone));
295: PetscCall(PetscPrintf(comm, " Edge %" PetscInt_FMT ": ", supp[s]));
296: PetscCall(PrintVertex(dm, cone[0]));
297: PetscCall(PetscPrintf(comm, " -- "));
298: PetscCall(PrintVertex(dm, cone[1]));
299: PetscCall(PetscPrintf(comm, "\n"));
300: }
301: }
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: static PetscErrorCode ComputeFFT(Mat FT, PetscInt Nc, Vec x, Vec p)
306: {
307: Vec *xComp, *pComp;
308: PetscInt n, N;
310: PetscFunctionBeginUser;
311: PetscCall(PetscMalloc2(Nc, &xComp, Nc, &pComp));
312: PetscCall(VecGetLocalSize(x, &n));
313: PetscCall(VecGetSize(x, &N));
314: for (PetscInt i = 0; i < Nc; ++i) {
315: const char *vtype;
317: // HACK: Make these from another DM up front
318: PetscCall(VecCreate(PetscObjectComm((PetscObject)x), &xComp[i]));
319: PetscCall(VecGetType(x, &vtype));
320: PetscCall(VecSetType(xComp[i], vtype));
321: PetscCall(VecSetSizes(xComp[i], n / Nc, N / Nc));
322: PetscCall(VecDuplicate(xComp[i], &pComp[i]));
323: }
324: PetscCall(VecStrideGatherAll(x, xComp, INSERT_VALUES));
325: for (PetscInt i = 0; i < Nc; ++i) PetscCall(MatMult(FT, xComp[i], pComp[i]));
326: PetscCall(VecStrideScatterAll(pComp, p, INSERT_VALUES));
327: for (PetscInt i = 0; i < Nc; ++i) {
328: PetscCall(VecDestroy(&xComp[i]));
329: PetscCall(VecDestroy(&pComp[i]));
330: }
331: PetscCall(PetscFree2(xComp, pComp));
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: // Sets each link to be the identity for the free field test
336: static PetscErrorCode SetGauge_Identity(DM dm)
337: {
338: DM auxDM;
339: Vec auxVec;
340: PetscSection s;
341: PetscScalar id[9] = {1., 0., 0., 0., 1., 0., 0., 0., 1.};
342: PetscInt eStart, eEnd;
344: PetscFunctionBegin;
345: PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &auxVec));
346: PetscCall(VecGetDM(auxVec, &auxDM));
347: PetscCall(DMGetLocalSection(auxDM, &s));
348: PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
349: for (PetscInt i = eStart; i < eEnd; ++i) PetscCall(VecSetValuesSection(auxVec, s, i, id, INSERT_VALUES));
350: PetscCall(VecViewFromOptions(auxVec, NULL, "-gauge_view"));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: /*
355: Test the action of the Wilson operator in the free field case U = I,
357: \eta(x) = D_W(x - y) \psi(y)
359: The Wilson operator is a convolution for the free field, so we can check that by the convolution theorem
361: \hat\eta(x) = \mathcal{F}(D_W(x - y) \psi(y))
362: = \hat D_W(p) \mathcal{F}\psi(p)
364: The Fourier series for the Wilson operator is
366: M + \sum_\mu 2 \sin^2(p_\mu / 2) + i \gamma_\mu \sin(p_\mu)
367: */
368: static PetscErrorCode TestFreeField(DM dm)
369: {
370: PetscSection s;
371: Mat FT;
372: Vec psi, psiHat;
373: Vec eta, etaHat;
374: Vec DHat; // The product \hat D_w \hat psi
375: PetscRandom r;
376: const PetscScalar *psih;
377: PetscScalar *dh;
378: PetscReal *coef, nrm;
379: const PetscInt *extent, Nc = 12;
380: PetscInt dim, V = 1, vStart, vEnd;
381: PetscContainer c;
382: PetscBool constRhs = PETSC_FALSE;
384: PetscFunctionBeginUser;
385: PetscCall(PetscOptionsGetBool(NULL, NULL, "-const_rhs", &constRhs, NULL));
387: PetscCall(SetGauge_Identity(dm));
388: PetscCall(DMGetLocalSection(dm, &s));
389: PetscCall(DMGetGlobalVector(dm, &psi));
390: PetscCall(PetscObjectSetName((PetscObject)psi, "psi"));
391: PetscCall(DMGetGlobalVector(dm, &psiHat));
392: PetscCall(PetscObjectSetName((PetscObject)psiHat, "psihat"));
393: PetscCall(DMGetGlobalVector(dm, &eta));
394: PetscCall(PetscObjectSetName((PetscObject)eta, "eta"));
395: PetscCall(DMGetGlobalVector(dm, &etaHat));
396: PetscCall(PetscObjectSetName((PetscObject)etaHat, "etahat"));
397: PetscCall(DMGetGlobalVector(dm, &DHat));
398: PetscCall(PetscObjectSetName((PetscObject)DHat, "Dhat"));
399: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
400: PetscCall(PetscRandomSetType(r, PETSCRAND48));
401: if (constRhs) PetscCall(VecSet(psi, 1.));
402: else PetscCall(VecSetRandom(psi, r));
403: PetscCall(PetscRandomDestroy(&r));
405: PetscCall(DMGetDimension(dm, &dim));
406: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
407: PetscCall(PetscObjectQuery((PetscObject)dm, "_extent", (PetscObject *)&c));
408: PetscCall(PetscContainerGetPointer(c, (void **)&extent));
409: PetscCall(MatCreateFFT(PetscObjectComm((PetscObject)dm), dim, extent, MATFFTW, &FT));
411: PetscCall(PetscMalloc1(dim, &coef));
412: for (PetscInt d = 0; d < dim; ++d) {
413: coef[d] = 2. * PETSC_PI / extent[d];
414: V *= extent[d];
415: }
416: PetscCall(ComputeResidual(dm, psi, eta));
417: PetscCall(VecViewFromOptions(eta, NULL, "-psi_view"));
418: PetscCall(VecViewFromOptions(eta, NULL, "-eta_view"));
419: PetscCall(ComputeFFT(FT, Nc, psi, psiHat));
420: PetscCall(VecScale(psiHat, 1. / V));
421: PetscCall(ComputeFFT(FT, Nc, eta, etaHat));
422: PetscCall(VecScale(etaHat, 1. / V));
423: PetscCall(VecGetArrayRead(psiHat, &psih));
424: PetscCall(VecGetArray(DHat, &dh));
425: for (PetscInt v = vStart; v < vEnd; ++v) {
426: PetscScalar tmp[12], tmp1 = 0.;
427: PetscInt dof, off;
429: PetscCall(PetscSectionGetDof(s, v, &dof));
430: PetscCall(PetscSectionGetOffset(s, v, &off));
431: for (PetscInt d = 0, prod = 1; d < dim; prod *= extent[d], ++d) {
432: const PetscInt idx = (v / prod) % extent[d];
434: tmp1 += 2. * PetscSqr(PetscSinReal(0.5 * coef[d] * idx));
435: for (PetscInt i = 0; i < dof; ++i) tmp[i] = PETSC_i * PetscSinReal(coef[d] * idx) * psih[off + i];
436: for (PetscInt c = 0; c < 3; ++c) PetscCall(ComputeGamma(d, 3, &tmp[c]));
437: for (PetscInt i = 0; i < dof; ++i) dh[off + i] += tmp[i];
438: }
439: for (PetscInt i = 0; i < dof; ++i) dh[off + i] += (M + tmp1) * psih[off + i];
440: }
441: PetscCall(VecRestoreArrayRead(psiHat, &psih));
442: PetscCall(VecRestoreArray(DHat, &dh));
444: {
445: Vec *etaComp, *DComp;
446: PetscInt n, N;
448: PetscCall(PetscMalloc2(Nc, &etaComp, Nc, &DComp));
449: PetscCall(VecGetLocalSize(etaHat, &n));
450: PetscCall(VecGetSize(etaHat, &N));
451: for (PetscInt i = 0; i < Nc; ++i) {
452: const char *vtype;
454: // HACK: Make these from another DM up front
455: PetscCall(VecCreate(PetscObjectComm((PetscObject)etaHat), &etaComp[i]));
456: PetscCall(VecGetType(etaHat, &vtype));
457: PetscCall(VecSetType(etaComp[i], vtype));
458: PetscCall(VecSetSizes(etaComp[i], n / Nc, N / Nc));
459: PetscCall(VecDuplicate(etaComp[i], &DComp[i]));
460: }
461: PetscCall(VecStrideGatherAll(etaHat, etaComp, INSERT_VALUES));
462: PetscCall(VecStrideGatherAll(DHat, DComp, INSERT_VALUES));
463: for (PetscInt i = 0; i < Nc; ++i) {
464: if (!i) {
465: PetscCall(VecViewFromOptions(etaComp[i], NULL, "-etahat_view"));
466: PetscCall(VecViewFromOptions(DComp[i], NULL, "-dhat_view"));
467: }
468: PetscCall(VecAXPY(etaComp[i], -1., DComp[i]));
469: PetscCall(VecNorm(etaComp[i], NORM_INFINITY, &nrm));
470: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Slice %" PetscInt_FMT ": %g\n", i, (double)nrm));
471: }
472: PetscCall(VecStrideScatterAll(etaComp, etaHat, INSERT_VALUES));
473: for (PetscInt i = 0; i < Nc; ++i) {
474: PetscCall(VecDestroy(&etaComp[i]));
475: PetscCall(VecDestroy(&DComp[i]));
476: }
477: PetscCall(PetscFree2(etaComp, DComp));
478: PetscCall(VecNorm(etaHat, NORM_INFINITY, &nrm));
479: PetscCheck(nrm < PETSC_SMALL, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Free field test failed: %g", (double)nrm);
480: }
482: PetscCall(PetscFree(coef));
483: PetscCall(MatDestroy(&FT));
484: PetscCall(DMRestoreGlobalVector(dm, &psi));
485: PetscCall(DMRestoreGlobalVector(dm, &psiHat));
486: PetscCall(DMRestoreGlobalVector(dm, &eta));
487: PetscCall(DMRestoreGlobalVector(dm, &etaHat));
488: PetscCall(DMRestoreGlobalVector(dm, &DHat));
489: PetscFunctionReturn(PETSC_SUCCESS);
490: }
492: int main(int argc, char **argv)
493: {
494: DM dm;
495: Vec u, f;
496: Mat J;
497: AppCtx user;
499: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
500: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
501: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
502: PetscCall(SetupDiscretization(dm, &user));
503: PetscCall(SetupAuxDiscretization(dm, &user));
504: PetscCall(DMCreateGlobalVector(dm, &u));
505: PetscCall(DMCreateGlobalVector(dm, &f));
506: PetscCall(PrintTraversal(dm));
507: PetscCall(ComputeResidual(dm, u, f));
508: PetscCall(VecViewFromOptions(f, NULL, "-res_view"));
509: PetscCall(CreateJacobian(dm, &J));
510: PetscCall(ComputeJacobian(dm, u, J));
511: PetscCall(VecViewFromOptions(u, NULL, "-sol_view"));
512: PetscCall(TestFreeField(dm));
513: PetscCall(VecDestroy(&f));
514: PetscCall(VecDestroy(&u));
515: PetscCall(DMDestroy(&dm));
516: PetscCall(PetscFinalize());
517: return 0;
518: }
520: /*TEST
522: build:
523: requires: complex
525: test:
526: requires: fftw
527: suffix: dirac_free_field
528: args: -dm_plex_dim 4 -dm_plex_shape hypercubic -dm_plex_box_faces 3,3,3,3 -dm_view -sol_view \
529: -dm_plex_check_symmetry -dm_plex_check_skeleton -dm_plex_check_faces -dm_plex_check_pointsf
531: TEST*/