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 usePV; /* Use Pauli-Villars preconditioning */ 15 } AppCtx; 16 17 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 18 { 19 PetscFunctionBegin; 20 options->usePV = PETSC_TRUE; 21 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(0); 26 } 27 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(0); 36 } 37 38 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx) 39 { 40 PetscSection s; 41 PetscInt vStart, vEnd, v; 42 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(0); 55 } 56 57 static PetscErrorCode SetupAuxDiscretization(DM dm, AppCtx *user) 58 { 59 DM dmAux, coordDM; 60 PetscSection s; 61 Vec gauge; 62 PetscInt eStart, eEnd, e; 63 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(0); 84 } 85 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; 92 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(0); 106 } 107 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]}; 112 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(0); 143 } 144 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]}; 150 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(0); 185 } 186 187 #include <petsc/private/dmpleximpl.h> 188 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]; 193 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(0); 205 } 206 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 PetscSection s; 213 const PetscScalar *ua; 214 PetscScalar *fa; 215 PetscScalar id[9] = {1., 0., 0., 0., 1., 0., 0., 0., 1.}; 216 PetscInt dim, vStart, vEnd; 217 218 PetscFunctionBeginUser; 219 PetscCall(DMGetDimension(dm, &dim)); 220 PetscCall(DMGetLocalSection(dm, &s)); 221 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 222 PetscCall(VecGetArrayRead(u, &ua)); 223 PetscCall(VecGetArray(f, &fa)); 224 // Loop over y 225 for (PetscInt v = vStart; v < vEnd; ++v) { 226 const PetscInt *supp; 227 PetscInt xdof, xoff; 228 229 PetscCall(DMPlexGetSupport(dm, v, &supp)); 230 PetscCall(PetscSectionGetDof(s, v, &xdof)); 231 PetscCall(PetscSectionGetOffset(s, v, &xoff)); 232 // Diagonal 233 for (PetscInt i = 0; i < xdof; ++i) fa[xoff + i] += (M + 4) * ua[xoff + i]; 234 // Loop over mu 235 for (PetscInt d = 0; d < dim; ++d) { 236 const PetscInt *cone; 237 PetscInt yoff; 238 239 // Left action -(1 + \gamma_\mu)/2 \otimes U^\dagger_\mu(y) \delta_{x - \mu,y} \psi(y) 240 PetscCall(DMPlexGetCone(dm, supp[2 * d + 0], &cone)); 241 PetscCall(PetscSectionGetOffset(s, cone[0], &yoff)); 242 PetscCall(ComputeAction(d, PETSC_FALSE, id, &ua[yoff], &fa[xoff])); 243 // Right edge -(1 - \gamma_\mu)/2 \otimes U_\mu(x) \delta_{x + \mu,y} \psi(y) 244 PetscCall(DMPlexGetCone(dm, supp[2 * d + 1], &cone)); 245 PetscCall(PetscSectionGetOffset(s, cone[1], &yoff)); 246 PetscCall(ComputeAction(d, PETSC_TRUE, id, &ua[yoff], &fa[xoff])); 247 } 248 } 249 PetscCall(VecRestoreArray(f, &fa)); 250 PetscCall(VecRestoreArrayRead(u, &ua)); 251 PetscFunctionReturn(0); 252 } 253 254 static PetscErrorCode CreateJacobian(DM dm, Mat *J) 255 { 256 PetscFunctionBegin; 257 PetscFunctionReturn(0); 258 } 259 260 static PetscErrorCode ComputeJacobian(DM dm, Vec u, Mat J) 261 { 262 PetscFunctionBegin; 263 PetscFunctionReturn(0); 264 } 265 266 static PetscErrorCode PrintTraversal(DM dm) 267 { 268 MPI_Comm comm; 269 PetscInt vStart, vEnd; 270 271 PetscFunctionBeginUser; 272 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 273 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 274 for (PetscInt v = vStart; v < vEnd; ++v) { 275 const PetscInt *supp; 276 PetscInt Ns; 277 278 PetscCall(DMPlexGetSupportSize(dm, v, &Ns)); 279 PetscCall(DMPlexGetSupport(dm, v, &supp)); 280 PetscCall(PrintVertex(dm, v)); 281 PetscCall(PetscPrintf(comm, "\n")); 282 for (PetscInt s = 0; s < Ns; ++s) { 283 const PetscInt *cone; 284 285 PetscCall(DMPlexGetCone(dm, supp[s], &cone)); 286 PetscCall(PetscPrintf(comm, " Edge %" PetscInt_FMT ": ", supp[s])); 287 PetscCall(PrintVertex(dm, cone[0])); 288 PetscCall(PetscPrintf(comm, " -- ")); 289 PetscCall(PrintVertex(dm, cone[1])); 290 PetscCall(PetscPrintf(comm, "\n")); 291 } 292 } 293 PetscFunctionReturn(0); 294 } 295 296 static PetscErrorCode ComputeFFT(Mat FT, PetscInt Nc, Vec x, Vec p) 297 { 298 Vec *xComp, *pComp; 299 PetscInt n, N; 300 301 PetscFunctionBeginUser; 302 PetscCall(PetscMalloc2(Nc, &xComp, Nc, &pComp)); 303 PetscCall(VecGetLocalSize(x, &n)); 304 PetscCall(VecGetSize(x, &N)); 305 for (PetscInt i = 0; i < Nc; ++i) { 306 const char *vtype; 307 308 // HACK: Make these from another DM up front 309 PetscCall(VecCreate(PetscObjectComm((PetscObject)x), &xComp[i])); 310 PetscCall(VecGetType(x, &vtype)); 311 PetscCall(VecSetType(xComp[i], vtype)); 312 PetscCall(VecSetSizes(xComp[i], n / Nc, N / Nc)); 313 PetscCall(VecDuplicate(xComp[i], &pComp[i])); 314 } 315 PetscCall(VecStrideGatherAll(x, xComp, INSERT_VALUES)); 316 for (PetscInt i = 0; i < Nc; ++i) PetscCall(MatMult(FT, xComp[i], pComp[i])); 317 PetscCall(VecStrideScatterAll(pComp, p, INSERT_VALUES)); 318 for (PetscInt i = 0; i < Nc; ++i) { 319 PetscCall(VecDestroy(&xComp[i])); 320 PetscCall(VecDestroy(&pComp[i])); 321 } 322 PetscCall(PetscFree2(xComp, pComp)); 323 PetscFunctionReturn(0); 324 } 325 326 /* 327 Test the action of the Wilson operator in the free field case U = I, 328 329 \eta(x) = D_W(x - y) \psi(y) 330 331 The Wilson operator is a convolution for the free field, so we can check that by the convolution theorem 332 333 \hat\eta(x) = \mathcal{F}(D_W(x - y) \psi(y)) 334 = \hat D_W(p) \mathcal{F}\psi(p) 335 336 The Fourier series for the Wilson operator is 337 338 M + \sum_\mu 2 \sin^2(p_\mu / 2) + i \gamma_\mu \sin(p_\mu) 339 */ 340 static PetscErrorCode TestFreeField(DM dm) 341 { 342 PetscSection s; 343 Mat FT; 344 Vec psi, psiHat; 345 Vec eta, etaHat; 346 Vec DHat; // The product \hat D_w \hat psi 347 PetscRandom r; 348 const PetscScalar *psih; 349 PetscScalar *dh; 350 PetscReal *coef, nrm; 351 const PetscInt *extent, Nc = 12; 352 PetscInt dim, V = 1, vStart, vEnd; 353 PetscContainer c; 354 PetscBool constRhs = PETSC_FALSE; 355 356 PetscFunctionBeginUser; 357 PetscCall(PetscOptionsGetBool(NULL, NULL, "-const_rhs", &constRhs, NULL)); 358 359 PetscCall(DMGetLocalSection(dm, &s)); 360 PetscCall(DMGetGlobalVector(dm, &psi)); 361 PetscCall(PetscObjectSetName((PetscObject)psi, "psi")); 362 PetscCall(DMGetGlobalVector(dm, &psiHat)); 363 PetscCall(PetscObjectSetName((PetscObject)psiHat, "psihat")); 364 PetscCall(DMGetGlobalVector(dm, &eta)); 365 PetscCall(PetscObjectSetName((PetscObject)eta, "eta")); 366 PetscCall(DMGetGlobalVector(dm, &etaHat)); 367 PetscCall(PetscObjectSetName((PetscObject)etaHat, "etahat")); 368 PetscCall(DMGetGlobalVector(dm, &DHat)); 369 PetscCall(PetscObjectSetName((PetscObject)DHat, "Dhat")); 370 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r)); 371 PetscCall(PetscRandomSetType(r, PETSCRAND48)); 372 if (constRhs) PetscCall(VecSet(psi, 1.)); 373 else PetscCall(VecSetRandom(psi, r)); 374 PetscCall(PetscRandomDestroy(&r)); 375 376 PetscCall(DMGetDimension(dm, &dim)); 377 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 378 PetscCall(PetscObjectQuery((PetscObject)dm, "_extent", (PetscObject *)&c)); 379 PetscCall(PetscContainerGetPointer(c, (void **)&extent)); 380 PetscCall(MatCreateFFT(PetscObjectComm((PetscObject)dm), dim, extent, MATFFTW, &FT)); 381 382 PetscCall(PetscMalloc1(dim, &coef)); 383 for (PetscInt d = 0; d < dim; ++d) { 384 coef[d] = 2. * PETSC_PI / extent[d]; 385 V *= extent[d]; 386 } 387 PetscCall(ComputeResidual(dm, psi, eta)); 388 PetscCall(VecViewFromOptions(eta, NULL, "-psi_view")); 389 PetscCall(VecViewFromOptions(eta, NULL, "-eta_view")); 390 PetscCall(ComputeFFT(FT, Nc, psi, psiHat)); 391 PetscCall(VecScale(psiHat, 1. / V)); 392 PetscCall(ComputeFFT(FT, Nc, eta, etaHat)); 393 PetscCall(VecScale(etaHat, 1. / V)); 394 PetscCall(VecGetArrayRead(psiHat, &psih)); 395 PetscCall(VecGetArray(DHat, &dh)); 396 for (PetscInt v = vStart; v < vEnd; ++v) { 397 PetscScalar tmp[12], tmp1 = 0.; 398 PetscInt dof, off; 399 400 PetscCall(PetscSectionGetDof(s, v, &dof)); 401 PetscCall(PetscSectionGetOffset(s, v, &off)); 402 for (PetscInt d = 0, prod = 1; d < dim; prod *= extent[d], ++d) { 403 const PetscInt idx = (v / prod) % extent[d]; 404 405 tmp1 += 2. * PetscSqr(PetscSinReal(0.5 * coef[d] * idx)); 406 for (PetscInt i = 0; i < dof; ++i) tmp[i] = PETSC_i * PetscSinReal(coef[d] * idx) * psih[off + i]; 407 for (PetscInt c = 0; c < 3; ++c) ComputeGamma(d, 3, &tmp[c]); 408 for (PetscInt i = 0; i < dof; ++i) dh[off + i] += tmp[i]; 409 } 410 for (PetscInt i = 0; i < dof; ++i) dh[off + i] += (M + tmp1) * psih[off + i]; 411 } 412 PetscCall(VecRestoreArrayRead(psiHat, &psih)); 413 PetscCall(VecRestoreArray(DHat, &dh)); 414 415 { 416 Vec *etaComp, *DComp; 417 PetscInt n, N; 418 419 PetscCall(PetscMalloc2(Nc, &etaComp, Nc, &DComp)); 420 PetscCall(VecGetLocalSize(etaHat, &n)); 421 PetscCall(VecGetSize(etaHat, &N)); 422 for (PetscInt i = 0; i < Nc; ++i) { 423 const char *vtype; 424 425 // HACK: Make these from another DM up front 426 PetscCall(VecCreate(PetscObjectComm((PetscObject)etaHat), &etaComp[i])); 427 PetscCall(VecGetType(etaHat, &vtype)); 428 PetscCall(VecSetType(etaComp[i], vtype)); 429 PetscCall(VecSetSizes(etaComp[i], n / Nc, N / Nc)); 430 PetscCall(VecDuplicate(etaComp[i], &DComp[i])); 431 } 432 PetscCall(VecStrideGatherAll(etaHat, etaComp, INSERT_VALUES)); 433 PetscCall(VecStrideGatherAll(DHat, DComp, INSERT_VALUES)); 434 for (PetscInt i = 0; i < Nc; ++i) { 435 if (!i) { 436 PetscCall(VecViewFromOptions(etaComp[i], NULL, "-etahat_view")); 437 PetscCall(VecViewFromOptions(DComp[i], NULL, "-dhat_view")); 438 } 439 PetscCall(VecAXPY(etaComp[i], -1., DComp[i])); 440 PetscCall(VecNorm(etaComp[i], NORM_INFINITY, &nrm)); 441 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Slice %" PetscInt_FMT ": %g\n", i, (double)nrm)); 442 } 443 PetscCall(VecStrideScatterAll(etaComp, etaHat, INSERT_VALUES)); 444 for (PetscInt i = 0; i < Nc; ++i) { 445 PetscCall(VecDestroy(&etaComp[i])); 446 PetscCall(VecDestroy(&DComp[i])); 447 } 448 PetscCall(PetscFree2(etaComp, DComp)); 449 PetscCall(VecNorm(etaHat, NORM_INFINITY, &nrm)); 450 PetscCheck(nrm < PETSC_SMALL, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Free field test failed: %g", (double)nrm); 451 } 452 453 PetscCall(PetscFree(coef)); 454 PetscCall(MatDestroy(&FT)); 455 PetscCall(DMRestoreGlobalVector(dm, &psi)); 456 PetscCall(DMRestoreGlobalVector(dm, &psiHat)); 457 PetscCall(DMRestoreGlobalVector(dm, &eta)); 458 PetscCall(DMRestoreGlobalVector(dm, &etaHat)); 459 PetscCall(DMRestoreGlobalVector(dm, &DHat)); 460 PetscFunctionReturn(0); 461 } 462 463 int main(int argc, char **argv) 464 { 465 DM dm; 466 Vec u, f; 467 Mat J; 468 AppCtx user; 469 470 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 471 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 472 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 473 PetscCall(SetupDiscretization(dm, &user)); 474 PetscCall(SetupAuxDiscretization(dm, &user)); 475 PetscCall(DMCreateGlobalVector(dm, &u)); 476 PetscCall(DMCreateGlobalVector(dm, &f)); 477 PetscCall(PrintTraversal(dm)); 478 PetscCall(ComputeResidual(dm, u, f)); 479 PetscCall(VecViewFromOptions(f, NULL, "-res_view")); 480 PetscCall(CreateJacobian(dm, &J)); 481 PetscCall(ComputeJacobian(dm, u, J)); 482 PetscCall(VecViewFromOptions(u, NULL, "-sol_view")); 483 PetscCall(TestFreeField(dm)); 484 PetscCall(VecDestroy(&f)); 485 PetscCall(VecDestroy(&u)); 486 PetscCall(DMDestroy(&dm)); 487 PetscCall(PetscFinalize()); 488 return 0; 489 } 490 491 /*TEST 492 493 build: 494 requires: complex 495 496 test: 497 requires: fftw 498 suffix: dirac_free_field 499 args: -dm_plex_dim 4 -dm_plex_shape hypercubic -dm_plex_box_faces 3,3,3,3 -dm_view -sol_view \ 500 -dm_plex_check_symmetry -dm_plex_check_skeleton -dm_plex_check_faces -dm_plex_check_pointsf 501 502 TEST*/ 503