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*/