Actual source code: ex1.c
1: static char help[] = "Tests projection with DMSwarm using general particle shapes.\n";
3: #include <petsc/private/dmswarmimpl.h>
4: #include <petsc/private/petscfeimpl.h>
6: #include <petscdmplex.h>
7: #include <petscds.h>
8: #include <petscksp.h>
10: typedef struct {
11: PetscInt dummy;
12: } AppCtx;
14: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
15: {
17: DMCreate(comm, dm);
18: DMSetType(*dm, DMPLEX);
19: DMSetFromOptions(*dm);
20: DMViewFromOptions(*dm, NULL, "-dm_view");
21: return 0;
22: }
24: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
25: {
26: PetscInt d;
27: u[0] = 0.0;
28: for (d = 0; d < dim; ++d) u[0] += x[d];
29: return 0;
30: }
32: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
33: {
34: g0[0] = 1.0;
35: }
37: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
38: {
39: PetscDS prob;
40: PetscFE fe;
41: PetscQuadrature quad;
42: PetscScalar *vals;
43: PetscReal *v0, *J, *invJ, detJ, *coords, *xi0;
44: PetscInt *cellid;
45: const PetscReal *qpoints;
46: PetscInt Ncell, c, Nq, q, dim;
47: PetscBool simplex;
50: DMGetDimension(dm, &dim);
51: DMPlexIsSimplex(dm, &simplex);
52: PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, NULL, -1, &fe);
53: DMGetDS(dm, &prob);
54: PetscDSSetDiscretization(prob, 0, (PetscObject)fe);
55: PetscFEDestroy(&fe);
56: PetscDSSetJacobian(prob, 0, 0, identity, NULL, NULL, NULL);
57: DMPlexGetHeightStratum(dm, 0, NULL, &Ncell);
58: PetscDSGetDiscretization(prob, 0, (PetscObject *)&fe);
59: PetscFEGetQuadrature(fe, &quad);
60: PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, NULL);
62: DMCreate(PetscObjectComm((PetscObject)dm), sw);
63: DMSetType(*sw, DMSWARM);
64: DMSetDimension(*sw, dim);
66: DMSwarmSetType(*sw, DMSWARM_PIC);
67: DMSwarmSetCellDM(*sw, dm);
68: DMSwarmRegisterPetscDatatypeField(*sw, "f_q", 1, PETSC_SCALAR);
69: DMSwarmFinalizeFieldRegister(*sw);
70: DMSwarmSetLocalSizes(*sw, Ncell * Nq, 0);
71: DMSetFromOptions(*sw);
73: PetscMalloc4(dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ);
74: for (c = 0; c < dim; c++) xi0[c] = -1.;
75: DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
76: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
77: DMSwarmGetField(*sw, "f_q", NULL, NULL, (void **)&vals);
78: for (c = 0; c < Ncell; ++c) {
79: for (q = 0; q < Nq; ++q) {
80: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
81: cellid[c * Nq + q] = c;
82: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q * dim], &coords[(c * Nq + q) * dim]);
83: linear(dim, 0.0, &coords[(c * Nq + q) * dim], 1, &vals[c * Nq + q], NULL);
84: }
85: }
86: DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
87: DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
88: DMSwarmRestoreField(*sw, "f_q", NULL, NULL, (void **)&vals);
89: PetscFree4(xi0, v0, J, invJ);
90: return 0;
91: }
93: static PetscErrorCode TestL2Projection(DM dm, DM sw, AppCtx *user)
94: {
95: PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
96: KSP ksp;
97: Mat mass;
98: Vec u, rhs, uproj;
99: PetscReal error;
102: funcs[0] = linear;
104: DMSwarmCreateGlobalVectorFromField(sw, "f_q", &u);
105: VecViewFromOptions(u, NULL, "-f_view");
106: DMGetGlobalVector(dm, &rhs);
107: DMCreateMassMatrix(sw, dm, &mass);
108: MatMult(mass, u, rhs);
109: MatDestroy(&mass);
110: VecDestroy(&u);
112: DMGetGlobalVector(dm, &uproj);
113: DMCreateMatrix(dm, &mass);
114: DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user);
115: MatViewFromOptions(mass, NULL, "-mass_mat_view");
116: KSPCreate(PETSC_COMM_WORLD, &ksp);
117: KSPSetOperators(ksp, mass, mass);
118: KSPSetFromOptions(ksp);
119: KSPSolve(ksp, rhs, uproj);
120: KSPDestroy(&ksp);
121: PetscObjectSetName((PetscObject)uproj, "Full Projection");
122: VecViewFromOptions(uproj, NULL, "-proj_vec_view");
123: DMComputeL2Diff(dm, 0.0, funcs, NULL, uproj, &error);
124: PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error);
125: MatDestroy(&mass);
126: DMRestoreGlobalVector(dm, &rhs);
127: DMRestoreGlobalVector(dm, &uproj);
128: return 0;
129: }
131: int main(int argc, char *argv[])
132: {
133: MPI_Comm comm;
134: DM dm, sw;
135: AppCtx user;
138: PetscInitialize(&argc, &argv, NULL, help);
139: comm = PETSC_COMM_WORLD;
140: CreateMesh(comm, &dm, &user);
141: CreateParticles(dm, &sw, &user);
142: PetscObjectSetName((PetscObject)dm, "Mesh");
143: DMViewFromOptions(dm, NULL, "-dm_view");
144: DMViewFromOptions(sw, NULL, "-sw_view");
145: TestL2Projection(dm, sw, &user);
146: DMDestroy(&dm);
147: DMDestroy(&sw);
148: PetscFinalize();
149: return 0;
150: }
152: /*TEST
154: test:
155: suffix: proj_0
156: requires: pragmatic
157: TODO: broken
158: args: -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic
159: test:
160: suffix: proj_1
161: requires: pragmatic
162: TODO: broken
163: args: -dm_plex_simplex 0 -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic
164: test:
165: suffix: proj_2
166: requires: pragmatic
167: TODO: broken
168: args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic
169: test:
170: suffix: proj_3
171: requires: pragmatic
172: TODO: broken
173: args: -dm_plex_simplex 0 -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic
175: TEST*/