1 static char help[] = "Fermions on a hypercubic lattice.\n\n";
2
3 #include <petscdmplex.h>
4 #include <petscsnes.h>
5
6 /* Common operations:
7
8 - View the input \psi as ASCII in lexicographic order: -psi_view
9 */
10
11 const PetscReal M = 1.0;
12
13 typedef struct {
14 PetscBool printTr; // Print the traversal
15 PetscBool usePV; // Use Pauli-Villars preconditioning
16 } AppCtx;
17
ProcessOptions(MPI_Comm comm,AppCtx * options)18 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
19 {
20 PetscFunctionBegin;
21 options->printTr = PETSC_FALSE;
22 options->usePV = PETSC_TRUE;
23
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 }
30
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)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 }
40
SetupDiscretization(DM dm,AppCtx * ctx)41 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
42 {
43 PetscSection s;
44 PetscInt vStart, vEnd, v;
45
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 }
59
SetupAuxDiscretization(DM dm,AppCtx * user)60 static PetscErrorCode SetupAuxDiscretization(DM dm, AppCtx *user)
61 {
62 DM dmAux, coordDM;
63 PetscSection s;
64 Vec gauge;
65 PetscInt eStart, eEnd, e;
66
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 }
88
PrintVertex(DM dm,PetscInt v)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;
95
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, &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 }
110
111 // Apply \gamma_\mu
ComputeGamma(PetscInt d,PetscInt ldx,PetscScalar f[])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]};
115
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 }
147
148 // Apply (1 \pm \gamma_\mu)/2
ComputeGammaFactor(PetscInt d,PetscBool forward,PetscInt ldx,PetscScalar f[])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]};
153
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 }
189
190 #include <petsc/private/dmpleximpl.h>
191
192 // ComputeAction() sums the action of 1/2 (1 \pm \gamma_\mu) U \psi into f
ComputeAction(PetscInt d,PetscBool forward,const PetscScalar U[],const PetscScalar psi[],PetscScalar f[])193 static PetscErrorCode ComputeAction(PetscInt d, PetscBool forward, const PetscScalar U[], const PetscScalar psi[], PetscScalar f[])
194 {
195 PetscScalar tmp[12];
196
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 }
209
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 */
ComputeResidualLocal(DM dm,Vec u,Vec f)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;
221
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));
228
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;
237
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;
247
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 }
265
ComputeResidual(DM dm,Vec u,Vec f)266 static PetscErrorCode ComputeResidual(DM dm, Vec u, Vec f)
267 {
268 Vec lu, lf;
269
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 }
283
CreateJacobian(DM dm,Mat * J)284 static PetscErrorCode CreateJacobian(DM dm, Mat *J)
285 {
286 PetscFunctionBegin;
287 PetscFunctionReturn(PETSC_SUCCESS);
288 }
289
ComputeJacobian(DM dm,Vec u,Mat J)290 static PetscErrorCode ComputeJacobian(DM dm, Vec u, Mat J)
291 {
292 PetscFunctionBegin;
293 PetscFunctionReturn(PETSC_SUCCESS);
294 }
295
PrintTraversal(DM dm)296 static PetscErrorCode PrintTraversal(DM dm)
297 {
298 MPI_Comm comm;
299 PetscInt vStart, vEnd;
300
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;
307
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;
314
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 }
325
ComputeFFT(Mat FT,PetscInt Nc,Vec x,Vec p)326 static PetscErrorCode ComputeFFT(Mat FT, PetscInt Nc, Vec x, Vec p)
327 {
328 Vec *xComp, *pComp, xTmp, pTmp;
329 PetscInt n, N;
330
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;
337
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 }
362
363 // Sets each link to be the identity for the free field test
SetGauge_Identity(DM dm)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;
371
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 }
381
382 /*
383 Test the action of the Wilson operator in the free field case U = I,
384
385 \eta(x) = D_W(x - y) \psi(y)
386
387 The Wilson operator is a convolution for the free field, so we can check that by the convolution theorem
388
389 \hat\eta(x) = \mathcal{F}(D_W(x - y) \psi(y))
390 = \hat D_W(p) \mathcal{F}\psi(p)
391
392 The Fourier series for the Wilson operator is
393
394 M + \sum_\mu 2 \sin^2(p_\mu / 2) + i \gamma_\mu \sin(p_\mu)
395 */
TestFreeField(DM dm)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;
412
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));
417
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));
435
436 PetscCall(DMGetDimension(dm, &dim));
437 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
438 PetscCall(PetscObjectQuery((PetscObject)dm, "_extent", (PetscObject *)&c));
439 PetscCall(PetscContainerGetPointer(c, &extent));
440 PetscCall(MatCreateFFT(PetscObjectComm((PetscObject)dm), dim, extent, MATFFTW, &FT));
441
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;
459
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];
464
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));
474
475 {
476 Vec *etaComp, *DComp;
477 PetscInt n, N;
478
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;
484
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 }
512
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 }
522
main(int argc,char ** argv)523 int main(int argc, char **argv)
524 {
525 DM dm;
526 Vec u, f;
527 Mat J;
528 AppCtx user;
529
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 }
550
551 /*TEST
552
553 build:
554 requires: complex
555
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
561
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
569
570 TEST*/
571