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 printTr; // Print the traversal
15: PetscBool usePV; // Use Pauli-Villars preconditioning
16: } AppCtx;
18: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
19: {
20: PetscFunctionBegin;
21: options->printTr = PETSC_FALSE;
22: options->usePV = PETSC_TRUE;
24: PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
25: PetscCall(PetscOptionsBool("-print_traversal", "Print the mesh traversal", "ex1.c", options->printTr, &options->printTr, NULL));
26: PetscCall(PetscOptionsBool("-use_pv", "Use Pauli-Villars preconditioning", "ex1.c", options->usePV, &options->usePV, NULL));
27: PetscOptionsEnd();
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
32: {
33: PetscFunctionBegin;
34: PetscCall(DMCreate(comm, dm));
35: PetscCall(DMSetType(*dm, DMPLEX));
36: PetscCall(DMSetFromOptions(*dm));
37: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
42: {
43: PetscSection s;
44: PetscInt vStart, vEnd, v;
46: PetscFunctionBegin;
47: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
48: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
49: PetscCall(PetscSectionSetChart(s, vStart, vEnd));
50: for (v = vStart; v < vEnd; ++v) {
51: PetscCall(PetscSectionSetDof(s, v, 12));
52: /* TODO Divide the values into fields/components */
53: }
54: PetscCall(PetscSectionSetUp(s));
55: PetscCall(DMSetLocalSection(dm, s));
56: PetscCall(PetscSectionDestroy(&s));
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: static PetscErrorCode SetupAuxDiscretization(DM dm, AppCtx *user)
61: {
62: DM dmAux, coordDM;
63: PetscSection s;
64: Vec gauge;
65: PetscInt eStart, eEnd, e;
67: PetscFunctionBegin;
68: /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
69: PetscCall(DMGetCoordinateDM(dm, &coordDM));
70: PetscCall(DMClone(dm, &dmAux));
71: PetscCall(DMSetCoordinateDM(dmAux, coordDM));
72: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
73: PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
74: PetscCall(PetscSectionSetChart(s, eStart, eEnd));
75: for (e = eStart; e < eEnd; ++e) {
76: /* TODO Should we store the whole SU(3) matrix, or the symmetric part? */
77: PetscCall(PetscSectionSetDof(s, e, 9));
78: }
79: PetscCall(PetscSectionSetUp(s));
80: PetscCall(DMSetLocalSection(dmAux, s));
81: PetscCall(PetscSectionDestroy(&s));
82: PetscCall(DMCreateLocalVector(dmAux, &gauge));
83: PetscCall(DMDestroy(&dmAux));
84: PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, gauge));
85: PetscCall(VecDestroy(&gauge));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: static PetscErrorCode PrintVertex(DM dm, PetscInt v)
90: {
91: MPI_Comm comm;
92: PetscContainer c;
93: PetscInt *extent;
94: PetscInt dim, cStart, cEnd, sum;
96: PetscFunctionBeginUser;
97: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
98: PetscCall(DMGetDimension(dm, &dim));
99: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
100: PetscCall(PetscObjectQuery((PetscObject)dm, "_extent", (PetscObject *)&c));
101: PetscCall(PetscContainerGetPointer(c, (void **)&extent));
102: sum = 1;
103: PetscCall(PetscPrintf(comm, "Vertex %" PetscInt_FMT ":", v));
104: for (PetscInt d = 0; d < dim; ++d) {
105: PetscCall(PetscPrintf(comm, " %" PetscInt_FMT, (v / sum) % extent[d]));
106: if (d < dim) sum *= extent[d];
107: }
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: // Apply \gamma_\mu
112: static PetscErrorCode ComputeGamma(PetscInt d, PetscInt ldx, PetscScalar f[])
113: {
114: const PetscScalar fin[4] = {f[0 * ldx], f[1 * ldx], f[2 * ldx], f[3 * ldx]};
116: PetscFunctionBeginHot;
117: switch (d) {
118: case 0:
119: f[0 * ldx] = PETSC_i * fin[3];
120: f[1 * ldx] = PETSC_i * fin[2];
121: f[2 * ldx] = -PETSC_i * fin[1];
122: f[3 * ldx] = -PETSC_i * fin[0];
123: break;
124: case 1:
125: f[0 * ldx] = -fin[3];
126: f[1 * ldx] = fin[2];
127: f[2 * ldx] = fin[1];
128: f[3 * ldx] = -fin[0];
129: break;
130: case 2:
131: f[0 * ldx] = PETSC_i * fin[2];
132: f[1 * ldx] = -PETSC_i * fin[3];
133: f[2 * ldx] = -PETSC_i * fin[0];
134: f[3 * ldx] = PETSC_i * fin[1];
135: break;
136: case 3:
137: f[0 * ldx] = fin[2];
138: f[1 * ldx] = fin[3];
139: f[2 * ldx] = fin[0];
140: f[3 * ldx] = fin[1];
141: break;
142: default:
143: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Direction for gamma %" PetscInt_FMT " not in [0, 4)", d);
144: }
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: // Apply (1 \pm \gamma_\mu)/2
149: static PetscErrorCode ComputeGammaFactor(PetscInt d, PetscBool forward, PetscInt ldx, PetscScalar f[])
150: {
151: const PetscReal sign = forward ? -1. : 1.;
152: const PetscScalar fin[4] = {f[0 * ldx], f[1 * ldx], f[2 * ldx], f[3 * ldx]};
154: PetscFunctionBeginHot;
155: switch (d) {
156: case 0:
157: f[0 * ldx] += sign * PETSC_i * fin[3];
158: f[1 * ldx] += sign * PETSC_i * fin[2];
159: f[2 * ldx] += sign * -PETSC_i * fin[1];
160: f[3 * ldx] += sign * -PETSC_i * fin[0];
161: break;
162: case 1:
163: f[0 * ldx] += -sign * fin[3];
164: f[1 * ldx] += sign * fin[2];
165: f[2 * ldx] += sign * fin[1];
166: f[3 * ldx] += -sign * fin[0];
167: break;
168: case 2:
169: f[0 * ldx] += sign * PETSC_i * fin[2];
170: f[1 * ldx] += sign * -PETSC_i * fin[3];
171: f[2 * ldx] += sign * -PETSC_i * fin[0];
172: f[3 * ldx] += sign * PETSC_i * fin[1];
173: break;
174: case 3:
175: f[0 * ldx] += sign * fin[2];
176: f[1 * ldx] += sign * fin[3];
177: f[2 * ldx] += sign * fin[0];
178: f[3 * ldx] += sign * fin[1];
179: break;
180: default:
181: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Direction for gamma %" PetscInt_FMT " not in [0, 4)", d);
182: }
183: f[0 * ldx] *= 0.5;
184: f[1 * ldx] *= 0.5;
185: f[2 * ldx] *= 0.5;
186: f[3 * ldx] *= 0.5;
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: #include <petsc/private/dmpleximpl.h>
192: // ComputeAction() sums the action of 1/2 (1 \pm \gamma_\mu) U \psi into f
193: static PetscErrorCode ComputeAction(PetscInt d, PetscBool forward, const PetscScalar U[], const PetscScalar psi[], PetscScalar f[])
194: {
195: PetscScalar tmp[12];
197: PetscFunctionBeginHot;
198: // Apply U
199: for (PetscInt beta = 0; beta < 4; ++beta) {
200: if (forward) DMPlex_Mult3D_Internal(U, 1, &psi[beta * 3], &tmp[beta * 3]);
201: else DMPlex_MultTranspose3D_Internal(U, 1, &psi[beta * 3], &tmp[beta * 3]);
202: }
203: // Apply (1 \pm \gamma_\mu)/2 to each color
204: for (PetscInt c = 0; c < 3; ++c) PetscCall(ComputeGammaFactor(d, forward, 3, &tmp[c]));
205: // Note that we are subtracting this contribution
206: for (PetscInt i = 0; i < 12; ++i) f[i] -= tmp[i];
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: /*
211: 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.
212: */
213: static PetscErrorCode ComputeResidualLocal(DM dm, Vec u, Vec f)
214: {
215: DM dmAux;
216: Vec gauge;
217: PetscSection s, sGauge;
218: const PetscScalar *ua;
219: PetscScalar *fa, *link;
220: PetscInt dim, vStart, vEnd;
222: PetscFunctionBeginUser;
223: PetscCall(DMGetDimension(dm, &dim));
224: PetscCall(DMGetLocalSection(dm, &s));
225: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
226: PetscCall(VecGetArrayRead(u, &ua));
227: PetscCall(VecGetArray(f, &fa));
229: PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &gauge));
230: PetscCall(VecGetDM(gauge, &dmAux));
231: PetscCall(DMGetLocalSection(dmAux, &sGauge));
232: PetscCall(VecGetArray(gauge, &link));
233: // Loop over y
234: for (PetscInt v = vStart; v < vEnd; ++v) {
235: const PetscInt *supp;
236: PetscInt xdof, xoff;
238: PetscCall(DMPlexGetSupport(dm, v, &supp));
239: PetscCall(PetscSectionGetDof(s, v, &xdof));
240: PetscCall(PetscSectionGetOffset(s, v, &xoff));
241: // Diagonal
242: for (PetscInt i = 0; i < xdof; ++i) fa[xoff + i] += (M + 4) * ua[xoff + i];
243: // Loop over mu
244: for (PetscInt d = 0; d < dim; ++d) {
245: const PetscInt *cone;
246: PetscInt yoff, goff;
248: // Left action -(1 + \gamma_\mu)/2 \otimes U^\dagger_\mu(y) \delta_{x - \mu,y} \psi(y)
249: PetscCall(DMPlexGetCone(dm, supp[2 * d + 0], &cone));
250: PetscCall(PetscSectionGetOffset(s, cone[0], &yoff));
251: PetscCall(PetscSectionGetOffset(sGauge, supp[2 * d + 0], &goff));
252: PetscCall(ComputeAction(d, PETSC_FALSE, &link[goff], &ua[yoff], &fa[xoff]));
253: // Right edge -(1 - \gamma_\mu)/2 \otimes U_\mu(x) \delta_{x + \mu,y} \psi(y)
254: PetscCall(DMPlexGetCone(dm, supp[2 * d + 1], &cone));
255: PetscCall(PetscSectionGetOffset(s, cone[1], &yoff));
256: PetscCall(PetscSectionGetOffset(sGauge, supp[2 * d + 1], &goff));
257: PetscCall(ComputeAction(d, PETSC_TRUE, &link[goff], &ua[yoff], &fa[xoff]));
258: }
259: }
260: PetscCall(VecRestoreArray(f, &fa));
261: PetscCall(VecRestoreArray(gauge, &link));
262: PetscCall(VecRestoreArrayRead(u, &ua));
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: static PetscErrorCode ComputeResidual(DM dm, Vec u, Vec f)
267: {
268: Vec lu, lf;
270: PetscFunctionBegin;
271: PetscCall(DMGetLocalVector(dm, &lu));
272: PetscCall(DMGetLocalVector(dm, &lf));
273: PetscCall(DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lu));
274: PetscCall(DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lu));
275: PetscCall(VecSet(lf, 0.));
276: PetscCall(ComputeResidualLocal(dm, lu, lf));
277: PetscCall(DMLocalToGlobalBegin(dm, lf, INSERT_VALUES, f));
278: PetscCall(DMLocalToGlobalEnd(dm, lf, INSERT_VALUES, f));
279: PetscCall(DMRestoreLocalVector(dm, &lu));
280: PetscCall(DMRestoreLocalVector(dm, &lf));
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: static PetscErrorCode CreateJacobian(DM dm, Mat *J)
285: {
286: PetscFunctionBegin;
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: static PetscErrorCode ComputeJacobian(DM dm, Vec u, Mat J)
291: {
292: PetscFunctionBegin;
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: static PetscErrorCode PrintTraversal(DM dm)
297: {
298: MPI_Comm comm;
299: PetscInt vStart, vEnd;
301: PetscFunctionBeginUser;
302: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
303: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
304: for (PetscInt v = vStart; v < vEnd; ++v) {
305: const PetscInt *supp;
306: PetscInt Ns;
308: PetscCall(DMPlexGetSupportSize(dm, v, &Ns));
309: PetscCall(DMPlexGetSupport(dm, v, &supp));
310: PetscCall(PrintVertex(dm, v));
311: PetscCall(PetscPrintf(comm, "\n"));
312: for (PetscInt s = 0; s < Ns; ++s) {
313: const PetscInt *cone;
315: PetscCall(DMPlexGetCone(dm, supp[s], &cone));
316: PetscCall(PetscPrintf(comm, " Edge %" PetscInt_FMT ": ", supp[s]));
317: PetscCall(PrintVertex(dm, cone[0]));
318: PetscCall(PetscPrintf(comm, " -- "));
319: PetscCall(PrintVertex(dm, cone[1]));
320: PetscCall(PetscPrintf(comm, "\n"));
321: }
322: }
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: static PetscErrorCode ComputeFFT(Mat FT, PetscInt Nc, Vec x, Vec p)
327: {
328: Vec *xComp, *pComp, xTmp, pTmp;
329: PetscInt n, N;
331: PetscFunctionBeginUser;
332: PetscCall(PetscMalloc2(Nc, &xComp, Nc, &pComp));
333: PetscCall(VecGetLocalSize(x, &n));
334: PetscCall(VecGetSize(x, &N));
335: for (PetscInt i = 0; i < Nc; ++i) {
336: const char *vtype;
338: // HACK: Make these from another DM up front
339: PetscCall(VecCreate(PetscObjectComm((PetscObject)x), &xComp[i]));
340: PetscCall(VecGetType(x, &vtype));
341: PetscCall(VecSetType(xComp[i], vtype));
342: PetscCall(VecSetSizes(xComp[i], n / Nc, N / Nc));
343: PetscCall(VecDuplicate(xComp[i], &pComp[i]));
344: }
345: PetscCall(MatCreateVecsFFTW(FT, &xTmp, &pTmp, NULL));
346: PetscCall(VecStrideGatherAll(x, xComp, INSERT_VALUES));
347: for (PetscInt i = 0; i < Nc; ++i) {
348: PetscCall(VecScatterPetscToFFTW(FT, xComp[i], xTmp));
349: PetscCall(MatMult(FT, xTmp, pTmp));
350: PetscCall(VecScatterFFTWToPetsc(FT, pTmp, pComp[i]));
351: }
352: PetscCall(VecStrideScatterAll(pComp, p, INSERT_VALUES));
353: for (PetscInt i = 0; i < Nc; ++i) {
354: PetscCall(VecDestroy(&xComp[i]));
355: PetscCall(VecDestroy(&pComp[i]));
356: }
357: PetscCall(VecDestroy(&xTmp));
358: PetscCall(VecDestroy(&pTmp));
359: PetscCall(PetscFree2(xComp, pComp));
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: // Sets each link to be the identity for the free field test
364: static PetscErrorCode SetGauge_Identity(DM dm)
365: {
366: DM auxDM;
367: Vec auxVec;
368: PetscSection s;
369: PetscScalar id[9] = {1., 0., 0., 0., 1., 0., 0., 0., 1.};
370: PetscInt eStart, eEnd;
372: PetscFunctionBegin;
373: PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &auxVec));
374: PetscCall(VecGetDM(auxVec, &auxDM));
375: PetscCall(DMGetLocalSection(auxDM, &s));
376: PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
377: for (PetscInt i = eStart; i < eEnd; ++i) PetscCall(VecSetValuesSection(auxVec, s, i, id, INSERT_VALUES));
378: PetscCall(VecViewFromOptions(auxVec, NULL, "-gauge_view"));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: /*
383: Test the action of the Wilson operator in the free field case U = I,
385: \eta(x) = D_W(x - y) \psi(y)
387: The Wilson operator is a convolution for the free field, so we can check that by the convolution theorem
389: \hat\eta(x) = \mathcal{F}(D_W(x - y) \psi(y))
390: = \hat D_W(p) \mathcal{F}\psi(p)
392: The Fourier series for the Wilson operator is
394: M + \sum_\mu 2 \sin^2(p_\mu / 2) + i \gamma_\mu \sin(p_\mu)
395: */
396: static PetscErrorCode TestFreeField(DM dm)
397: {
398: PetscSection s;
399: Mat FT;
400: Vec psi, psiHat;
401: Vec eta, etaHat;
402: Vec DHat; // The product \hat D_w \hat psi
403: PetscRandom r;
404: const PetscScalar *psih;
405: PetscScalar *dh;
406: PetscReal *coef, nrm;
407: const PetscInt *extent, Nc = 12;
408: PetscInt dim, V = 1, vStart, vEnd;
409: PetscContainer c;
410: PetscBool constRhs = PETSC_FALSE;
411: PetscMPIInt size;
413: PetscFunctionBeginUser;
414: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
415: if (size > 1) PetscFunctionReturn(PETSC_SUCCESS);
416: PetscCall(PetscOptionsGetBool(NULL, NULL, "-const_rhs", &constRhs, NULL));
418: PetscCall(SetGauge_Identity(dm));
419: PetscCall(DMGetLocalSection(dm, &s));
420: PetscCall(DMGetGlobalVector(dm, &psi));
421: PetscCall(PetscObjectSetName((PetscObject)psi, "psi"));
422: PetscCall(DMGetGlobalVector(dm, &psiHat));
423: PetscCall(PetscObjectSetName((PetscObject)psiHat, "psihat"));
424: PetscCall(DMGetGlobalVector(dm, &eta));
425: PetscCall(PetscObjectSetName((PetscObject)eta, "eta"));
426: PetscCall(DMGetGlobalVector(dm, &etaHat));
427: PetscCall(PetscObjectSetName((PetscObject)etaHat, "etahat"));
428: PetscCall(DMGetGlobalVector(dm, &DHat));
429: PetscCall(PetscObjectSetName((PetscObject)DHat, "Dhat"));
430: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
431: PetscCall(PetscRandomSetType(r, PETSCRAND48));
432: if (constRhs) PetscCall(VecSet(psi, 1.));
433: else PetscCall(VecSetRandom(psi, r));
434: PetscCall(PetscRandomDestroy(&r));
436: PetscCall(DMGetDimension(dm, &dim));
437: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
438: PetscCall(PetscObjectQuery((PetscObject)dm, "_extent", (PetscObject *)&c));
439: PetscCall(PetscContainerGetPointer(c, (void **)&extent));
440: PetscCall(MatCreateFFT(PetscObjectComm((PetscObject)dm), dim, extent, MATFFTW, &FT));
442: PetscCall(PetscMalloc1(dim, &coef));
443: for (PetscInt d = 0; d < dim; ++d) {
444: coef[d] = 2. * PETSC_PI / extent[d];
445: V *= extent[d];
446: }
447: PetscCall(ComputeResidual(dm, psi, eta));
448: PetscCall(VecViewFromOptions(psi, NULL, "-psi_view"));
449: PetscCall(VecViewFromOptions(eta, NULL, "-eta_view"));
450: PetscCall(ComputeFFT(FT, Nc, psi, psiHat));
451: PetscCall(VecScale(psiHat, 1. / V));
452: PetscCall(ComputeFFT(FT, Nc, eta, etaHat));
453: PetscCall(VecScale(etaHat, 1. / V));
454: PetscCall(VecGetArrayRead(psiHat, &psih));
455: PetscCall(VecGetArray(DHat, &dh));
456: for (PetscInt v = vStart; v < vEnd; ++v) {
457: PetscScalar tmp[12], tmp1 = 0.;
458: PetscInt dof, off;
460: PetscCall(PetscSectionGetDof(s, v, &dof));
461: PetscCall(PetscSectionGetOffset(s, v, &off));
462: for (PetscInt d = 0, prod = 1; d < dim; prod *= extent[d], ++d) {
463: const PetscInt idx = (v / prod) % extent[d];
465: tmp1 += 2. * PetscSqr(PetscSinReal(0.5 * coef[d] * idx));
466: for (PetscInt i = 0; i < dof; ++i) tmp[i] = PETSC_i * PetscSinReal(coef[d] * idx) * psih[off + i];
467: for (PetscInt c = 0; c < 3; ++c) PetscCall(ComputeGamma(d, 3, &tmp[c]));
468: for (PetscInt i = 0; i < dof; ++i) dh[off + i] += tmp[i];
469: }
470: for (PetscInt i = 0; i < dof; ++i) dh[off + i] += (M + tmp1) * psih[off + i];
471: }
472: PetscCall(VecRestoreArrayRead(psiHat, &psih));
473: PetscCall(VecRestoreArray(DHat, &dh));
475: {
476: Vec *etaComp, *DComp;
477: PetscInt n, N;
479: PetscCall(PetscMalloc2(Nc, &etaComp, Nc, &DComp));
480: PetscCall(VecGetLocalSize(etaHat, &n));
481: PetscCall(VecGetSize(etaHat, &N));
482: for (PetscInt i = 0; i < Nc; ++i) {
483: const char *vtype;
485: // HACK: Make these from another DM up front
486: PetscCall(VecCreate(PetscObjectComm((PetscObject)etaHat), &etaComp[i]));
487: PetscCall(VecGetType(etaHat, &vtype));
488: PetscCall(VecSetType(etaComp[i], vtype));
489: PetscCall(VecSetSizes(etaComp[i], n / Nc, N / Nc));
490: PetscCall(VecDuplicate(etaComp[i], &DComp[i]));
491: }
492: PetscCall(VecStrideGatherAll(etaHat, etaComp, INSERT_VALUES));
493: PetscCall(VecStrideGatherAll(DHat, DComp, INSERT_VALUES));
494: for (PetscInt i = 0; i < Nc; ++i) {
495: if (!i) {
496: PetscCall(VecViewFromOptions(etaComp[i], NULL, "-etahat_view"));
497: PetscCall(VecViewFromOptions(DComp[i], NULL, "-dhat_view"));
498: }
499: PetscCall(VecAXPY(etaComp[i], -1., DComp[i]));
500: PetscCall(VecNorm(etaComp[i], NORM_INFINITY, &nrm));
501: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Slice %" PetscInt_FMT ": %g\n", i, (double)nrm));
502: }
503: PetscCall(VecStrideScatterAll(etaComp, etaHat, INSERT_VALUES));
504: for (PetscInt i = 0; i < Nc; ++i) {
505: PetscCall(VecDestroy(&etaComp[i]));
506: PetscCall(VecDestroy(&DComp[i]));
507: }
508: PetscCall(PetscFree2(etaComp, DComp));
509: PetscCall(VecNorm(etaHat, NORM_INFINITY, &nrm));
510: PetscCheck(nrm < PETSC_SMALL, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Free field test failed: %g", (double)nrm);
511: }
513: PetscCall(PetscFree(coef));
514: PetscCall(MatDestroy(&FT));
515: PetscCall(DMRestoreGlobalVector(dm, &psi));
516: PetscCall(DMRestoreGlobalVector(dm, &psiHat));
517: PetscCall(DMRestoreGlobalVector(dm, &eta));
518: PetscCall(DMRestoreGlobalVector(dm, &etaHat));
519: PetscCall(DMRestoreGlobalVector(dm, &DHat));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: int main(int argc, char **argv)
524: {
525: DM dm;
526: Vec u, f;
527: Mat J;
528: AppCtx user;
530: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
531: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
532: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
533: PetscCall(SetupDiscretization(dm, &user));
534: PetscCall(SetupAuxDiscretization(dm, &user));
535: PetscCall(DMCreateGlobalVector(dm, &u));
536: PetscCall(DMCreateGlobalVector(dm, &f));
537: if (user.printTr) PetscCall(PrintTraversal(dm));
538: PetscCall(ComputeResidual(dm, u, f));
539: PetscCall(VecViewFromOptions(f, NULL, "-res_view"));
540: PetscCall(CreateJacobian(dm, &J));
541: PetscCall(ComputeJacobian(dm, u, J));
542: PetscCall(VecViewFromOptions(u, NULL, "-sol_view"));
543: PetscCall(TestFreeField(dm));
544: PetscCall(VecDestroy(&f));
545: PetscCall(VecDestroy(&u));
546: PetscCall(DMDestroy(&dm));
547: PetscCall(PetscFinalize());
548: return 0;
549: }
551: /*TEST
553: build:
554: requires: complex
556: test:
557: requires: fftw
558: suffix: dirac_free_field
559: args: -dm_plex_dim 4 -dm_plex_shape hypercubic -dm_plex_box_faces 3,3,3,3 -dm_view -sol_view \
560: -dm_plex_check_symmetry -dm_plex_check_skeleton -dm_plex_check_faces -dm_plex_check_pointsf
562: # Problem on rank 3 on MacOSX
563: test:
564: requires: fftw
565: suffix: dirac_free_field_par
566: nsize: 16
567: args: -dm_plex_dim 4 -dm_plex_shape hypercubic -dm_plex_box_faces 4,4,4,4 -dm_view -sol_view \
568: -dm_plex_check_symmetry -dm_plex_check_skeleton -dm_plex_check_faces -dm_plex_check_pointsf
570: TEST*/