Actual source code: ex1.c
1: static char help[] = "Tests projection with DMSwarm using general particle shapes.\n";
3: #include <petsc/private/petscfeimpl.h>
4: #include <petscdmswarm.h>
5: #include <petscdmplex.h>
6: #include <petscds.h>
7: #include <petscksp.h>
9: typedef struct {
10: PetscReal L[3]; /* Dimensions of the mesh bounding box */
11: } AppCtx;
13: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
14: {
15: PetscReal low[3], high[3];
16: PetscInt dim;
18: PetscFunctionBeginUser;
19: PetscCall(DMCreate(comm, dm));
20: PetscCall(DMSetType(*dm, DMPLEX));
21: PetscCall(DMSetFromOptions(*dm));
22: PetscCall(DMGetBoundingBox(*dm, low, high));
23: PetscCall(DMGetDimension(*dm, &dim));
24: for (PetscInt d = 0; d < dim; ++d) user->L[d] = high[d] - low[d];
25: PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
26: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
31: {
32: AppCtx *ctx = (AppCtx *)a_ctx;
33: PetscInt d;
34: u[0] = 0.0;
35: for (d = 0; d < dim; ++d) u[0] += x[d] / ctx->L[d];
36: return PETSC_SUCCESS;
37: }
39: 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[])
40: {
41: g0[0] = 1.0;
42: }
44: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
45: {
46: PetscFE fe;
47: PetscDS ds;
48: PetscBool simplex;
49: PetscInt dim;
51: PetscFunctionBeginUser;
52: PetscCall(DMGetDimension(dm, &dim));
53: PetscCall(DMPlexIsSimplex(dm, &simplex));
54: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, NULL, -1, &fe));
55: PetscCall(PetscObjectSetName((PetscObject)fe, "fe"));
56: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
57: PetscCall(DMCreateDS(dm));
58: PetscCall(PetscFEDestroy(&fe));
59: /* Setup to form mass matrix */
60: PetscCall(DMGetDS(dm, &ds));
61: PetscCall(PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL));
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
66: {
67: PetscDS prob;
68: PetscFE fe;
69: PetscQuadrature quad;
70: PetscScalar *vals;
71: PetscReal *v0, *J, *invJ, detJ, *coords, *xi0;
72: PetscInt *cellid;
73: const PetscReal *qpoints, *qweights;
74: PetscInt cStart, cEnd, c, Nq, q, dim;
76: PetscFunctionBeginUser;
77: PetscCall(DMGetDimension(dm, &dim));
78: PetscCall(DMGetDS(dm, &prob));
79: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
80: PetscCall(PetscDSGetDiscretization(prob, 0, (PetscObject *)&fe));
81: PetscCall(PetscFEGetQuadrature(fe, &quad));
82: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, &qweights));
83: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
84: PetscCall(DMSetType(*sw, DMSWARM));
85: PetscCall(DMSetDimension(*sw, dim));
86: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
87: PetscCall(DMSwarmSetCellDM(*sw, dm));
88: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "f_q", 1, PETSC_REAL));
89: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
90: PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Nq, 0));
91: PetscCall(DMSetFromOptions(*sw));
92: PetscCall(PetscMalloc4(dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
93: for (c = 0; c < dim; c++) xi0[c] = -1.;
94: PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
95: PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
96: PetscCall(DMSwarmGetField(*sw, "f_q", NULL, NULL, (void **)&vals));
97: for (c = cStart; c < cEnd; ++c) {
98: for (q = 0; q < Nq; ++q) {
99: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
100: cellid[c * Nq + q] = c;
101: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q * dim], &coords[(c * Nq + q) * dim]);
102: PetscCall(linear(dim, 0.0, &coords[(c * Nq + q) * dim], 1, &vals[c * Nq + q], (void *)user));
103: vals[c * Nq + q] *= qweights[q] * detJ;
104: }
105: }
106: PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
107: PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
108: PetscCall(DMSwarmRestoreField(*sw, "f_q", NULL, NULL, (void **)&vals));
109: PetscCall(PetscFree4(xi0, v0, J, invJ));
110: PetscCall(DMSwarmMigrate(*sw, PETSC_FALSE));
111: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
112: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode TestL2Projection(DM dm, DM sw, AppCtx *user)
117: {
118: PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
119: KSP ksp;
120: Mat mass;
121: Vec u, rhs, uproj;
122: void *ctxs[1];
123: PetscReal error;
125: PetscFunctionBeginUser;
126: ctxs[0] = (void *)user;
127: funcs[0] = linear;
129: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "f_q", &u));
130: PetscCall(VecViewFromOptions(u, NULL, "-f_view"));
131: PetscCall(DMGetGlobalVector(dm, &rhs));
132: PetscCall(DMCreateMassMatrix(sw, dm, &mass));
133: PetscCall(MatMultTranspose(mass, u, rhs));
134: PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));
135: PetscCall(MatDestroy(&mass));
136: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "f_q", &u));
137: PetscCall(DMGetGlobalVector(dm, &uproj));
138: PetscCall(DMCreateMatrix(dm, &mass));
139: PetscCall(DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user));
140: PetscCall(MatViewFromOptions(mass, NULL, "-mass_mat_view"));
141: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
142: PetscCall(KSPSetOperators(ksp, mass, mass));
143: PetscCall(KSPSetFromOptions(ksp));
144: PetscCall(KSPSolve(ksp, rhs, uproj));
145: PetscCall(KSPDestroy(&ksp));
146: PetscCall(PetscObjectSetName((PetscObject)uproj, "Full Projection"));
147: PetscCall(VecViewFromOptions(uproj, NULL, "-proj_vec_view"));
148: PetscCall(DMComputeL2Diff(dm, 0.0, funcs, ctxs, uproj, &error));
149: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error));
150: PetscCall(MatDestroy(&mass));
151: PetscCall(DMRestoreGlobalVector(dm, &rhs));
152: PetscCall(DMRestoreGlobalVector(dm, &uproj));
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: int main(int argc, char *argv[])
157: {
158: MPI_Comm comm;
159: DM dm, sw;
160: AppCtx user;
162: PetscFunctionBeginUser;
163: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
164: comm = PETSC_COMM_WORLD;
165: PetscCall(CreateMesh(comm, &dm, &user));
166: PetscCall(CreateFEM(dm, &user));
167: PetscCall(CreateParticles(dm, &sw, &user));
168: PetscCall(TestL2Projection(dm, sw, &user));
169: PetscCall(DMDestroy(&dm));
170: PetscCall(DMDestroy(&sw));
171: PetscCall(PetscFinalize());
172: return 0;
173: }
175: /*TEST
177: build:
178: requires: double !complex
179: test:
180: suffix: proj_0
181: requires: triangle
182: args: -dm_plex_box_faces 2,2\
183: -dm_view -sw_view\
184: -petscspace_degree 1 -petscfe_default_quadrature_order 1\
185: -pc_type lu
186: test:
187: suffix: proj_1
188: requires: triangle
189: args: -dm_plex_simplex 0 -dm_plex_box_faces 2,2\
190: -dm_view -sw_view\
191: -petscspace_degree 1 -petscfe_default_quadrature_order 1\
192: -pc_type lu
193: test:
194: suffix: proj_2
195: requires: ctetgen
196: args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2\
197: -dm_view -sw_view\
198: -petscspace_degree 1 -petscfe_default_quadrature_order 1\
199: -pc_type lu
200: test:
201: suffix: proj_3
202: requires: ctetgen
203: args: -dm_plex_dim 3 -dm_plex_simplex 0\
204: -dm_plex_box_faces 2,2,2\
205: -dm_view -sw_view\
206: -petscspace_degree 1 -petscfe_default_quadrature_order 1\
207: -pc_type lu
209: TEST*/