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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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 DM dmAux; 213 Vec gauge; 214 PetscSection s, sGauge; 215 const PetscScalar *ua; 216 PetscScalar *fa, *link; 217 PetscInt dim, vStart, vEnd; 218 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)); 225 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; 234 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; 244 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 } 262 263 static PetscErrorCode CreateJacobian(DM dm, Mat *J) 264 { 265 PetscFunctionBegin; 266 PetscFunctionReturn(PETSC_SUCCESS); 267 } 268 269 static PetscErrorCode ComputeJacobian(DM dm, Vec u, Mat J) 270 { 271 PetscFunctionBegin; 272 PetscFunctionReturn(PETSC_SUCCESS); 273 } 274 275 static PetscErrorCode PrintTraversal(DM dm) 276 { 277 MPI_Comm comm; 278 PetscInt vStart, vEnd; 279 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; 286 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; 293 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 } 304 305 static PetscErrorCode ComputeFFT(Mat FT, PetscInt Nc, Vec x, Vec p) 306 { 307 Vec *xComp, *pComp; 308 PetscInt n, N; 309 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; 316 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 } 334 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; 343 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 } 353 354 /* 355 Test the action of the Wilson operator in the free field case U = I, 356 357 \eta(x) = D_W(x - y) \psi(y) 358 359 The Wilson operator is a convolution for the free field, so we can check that by the convolution theorem 360 361 \hat\eta(x) = \mathcal{F}(D_W(x - y) \psi(y)) 362 = \hat D_W(p) \mathcal{F}\psi(p) 363 364 The Fourier series for the Wilson operator is 365 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; 383 384 PetscFunctionBeginUser; 385 PetscCall(PetscOptionsGetBool(NULL, NULL, "-const_rhs", &constRhs, NULL)); 386 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)); 404 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)); 410 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; 428 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]; 433 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)); 443 444 { 445 Vec *etaComp, *DComp; 446 PetscInt n, N; 447 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; 453 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 } 481 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 } 491 492 int main(int argc, char **argv) 493 { 494 DM dm; 495 Vec u, f; 496 Mat J; 497 AppCtx user; 498 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 } 519 520 /*TEST 521 522 build: 523 requires: complex 524 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 530 531 TEST*/ 532