Actual source code: ex11.c
1: static char help[] = "Tests multifield and multicomponent L2 projection.\n";
3: #include <petscdmswarm.h>
4: #include <petscksp.h>
5: #include <petscdmplex.h>
6: #include <petscds.h>
8: typedef struct {
9: PetscBool pass; // Don't fail when moments are wrong
10: PetscBool fv; // Use an FV discretization, instead of FE
11: PetscInt Npc; // The number of partices per cell
12: PetscInt field; // The field to project
13: PetscInt Nm; // The number of moments to match
14: PetscReal mtol; // Tolerance for checking moment conservation
15: PetscSimplePointFn *func; // Function used to set particle weights
16: } AppCtx;
18: typedef enum {
19: FUNCTION_CONSTANT,
20: FUNCTION_LINEAR,
21: FUNCTION_SIN,
22: FUNCTION_X2_X4,
23: FUNCTION_UNKNOWN,
24: NUM_FUNCTIONS
25: } FunctionType;
26: const char *const FunctionTypes[] = {"constant", "linear", "sin", "x2_x4", "unknown", "FunctionType", "FUNCTION_", NULL};
28: static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
29: {
30: u[0] = 1.0;
31: return PETSC_SUCCESS;
32: }
34: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
35: {
36: u[0] = 0.0;
37: for (PetscInt d = 0; d < dim; ++d) u[0] += x[d];
38: return PETSC_SUCCESS;
39: }
41: static PetscErrorCode sinx(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
42: {
43: u[0] = 1;
44: for (PetscInt d = 0; d < dim; ++d) u[0] *= PetscSinReal(2. * PETSC_PI * x[d]);
45: return PETSC_SUCCESS;
46: }
48: static PetscErrorCode x2_x4(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
49: {
50: u[0] = 1;
51: for (PetscInt d = 0; d < dim; ++d) u[0] *= PetscSqr(x[d]) - PetscSqr(PetscSqr(x[d]));
52: return PETSC_SUCCESS;
53: }
55: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
56: {
57: FunctionType func = FUNCTION_LINEAR;
58: PetscBool flg;
60: PetscFunctionBeginUser;
61: options->pass = PETSC_FALSE;
62: options->fv = PETSC_FALSE;
63: options->Npc = 1;
64: options->field = 0;
65: options->Nm = 1;
66: options->mtol = 100. * PETSC_MACHINE_EPSILON;
68: PetscOptionsBegin(comm, "", "L2 Projection Options", "DMPLEX");
69: PetscCall(PetscOptionsBool("-pass", "Don't fail when moments are wrong", __FILE__, options->pass, &options->pass, NULL));
70: PetscCall(PetscOptionsBool("-fv", "Use FV instead of FE", __FILE__, options->fv, &options->fv, NULL));
71: PetscCall(PetscOptionsBoundedInt("-npc", "Number of particles per cell", __FILE__, options->Npc, &options->Npc, NULL, 1));
72: PetscCall(PetscOptionsBoundedInt("-field", "The field to project", __FILE__, options->field, &options->field, NULL, 0));
73: PetscCall(PetscOptionsBoundedInt("-moments", "Number of moments to match", __FILE__, options->Nm, &options->Nm, NULL, 0));
74: PetscCheck(options->Nm < 4, comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot match %" PetscInt_FMT " > 3 moments", options->Nm);
75: PetscCall(PetscOptionsReal("-mtol", "Tolerance for moment checks", "ex2.c", options->mtol, &options->mtol, NULL));
76: PetscCall(PetscOptionsEnum("-func", "Type of particle weight function", __FILE__, FunctionTypes, (PetscEnum)func, (PetscEnum *)&func, &flg));
77: switch (func) {
78: case FUNCTION_CONSTANT:
79: options->func = constant;
80: break;
81: case FUNCTION_LINEAR:
82: options->func = linear;
83: break;
84: case FUNCTION_SIN:
85: options->func = sinx;
86: break;
87: case FUNCTION_X2_X4:
88: options->func = x2_x4;
89: break;
90: default:
91: PetscCheck(flg, comm, PETSC_ERR_ARG_WRONG, "Cannot handle function \"%s\"", FunctionTypes[func]);
92: }
93: PetscOptionsEnd();
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
98: {
99: PetscFunctionBeginUser;
100: PetscCall(DMCreate(comm, dm));
101: PetscCall(DMSetType(*dm, DMPLEX));
102: PetscCall(DMSetFromOptions(*dm));
103: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user)
108: {
109: PetscFE fe;
110: PetscFV fv;
111: DMPolytopeType ct;
112: PetscInt dim, cStart;
114: PetscFunctionBeginUser;
115: PetscCall(DMGetDimension(dm, &dim));
116: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
117: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
118: if (user->fv) {
119: PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv));
120: PetscCall(PetscObjectSetName((PetscObject)fv, "fv"));
121: PetscCall(PetscFVSetNumComponents(fv, 1));
122: PetscCall(PetscFVSetSpatialDimension(fv, dim));
123: PetscCall(PetscFVCreateDualSpace(fv, ct));
124: PetscCall(PetscFVSetFromOptions(fv));
125: PetscCall(DMAddField(dm, NULL, (PetscObject)fv));
126: PetscCall(PetscFVDestroy(&fv));
127: PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv));
128: PetscCall(PetscObjectSetName((PetscObject)fv, "fv2"));
129: PetscCall(PetscFVSetNumComponents(fv, 2));
130: PetscCall(PetscFVSetSpatialDimension(fv, dim));
131: PetscCall(PetscFVCreateDualSpace(fv, ct));
132: PetscCall(PetscFVSetFromOptions(fv));
133: PetscCall(DMAddField(dm, NULL, (PetscObject)fv));
134: PetscCall(PetscFVDestroy(&fv));
135: } else {
136: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, -1, &fe));
137: PetscCall(PetscObjectSetName((PetscObject)fe, "fe"));
138: PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
139: PetscCall(PetscFEDestroy(&fe));
140: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 2, ct, NULL, -1, &fe));
141: PetscCall(PetscObjectSetName((PetscObject)fe, "fe2"));
142: PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
143: PetscCall(PetscFEDestroy(&fe));
144: }
145: PetscCall(DMCreateDS(dm));
146: if (user->fv) {
147: DMLabel label;
148: PetscInt values[1] = {1};
150: PetscCall(DMCreateLabel(dm, "ghost"));
151: PetscCall(DMGetLabel(dm, "marker", &label));
152: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "dummy", label, 1, values, 0, 0, NULL, NULL, NULL, NULL, NULL));
153: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "dummy", label, 1, values, 1, 0, NULL, NULL, NULL, NULL, NULL));
154: }
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: static PetscErrorCode CreateSwarm(DM dm, DM *sw, AppCtx *user)
159: {
160: PetscScalar *coords, *wvals, *xvals;
161: PetscInt Npc = user->Npc, dim, Np;
163: PetscFunctionBeginUser;
164: PetscCall(DMGetDimension(dm, &dim));
166: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
167: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
168: PetscCall(DMSetType(*sw, DMSWARM));
169: PetscCall(DMSetDimension(*sw, dim));
170: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
171: PetscCall(DMSwarmSetCellDM(*sw, dm));
172: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
173: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "x_q", 2, PETSC_SCALAR));
174: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
175: PetscCall(DMSwarmInsertPointsUsingCellDM(*sw, DMSWARMPIC_LAYOUT_GAUSS, Npc));
176: PetscCall(DMSetFromOptions(*sw));
178: PetscCall(DMSwarmGetLocalSize(*sw, &Np));
179: PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
180: PetscCall(DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **)&wvals));
181: PetscCall(DMSwarmGetField(*sw, "x_q", NULL, NULL, (void **)&xvals));
182: for (PetscInt p = 0; p < Np; ++p) {
183: PetscCall(user->func(dim, 0., &coords[p * dim], 1, &wvals[p], user));
184: for (PetscInt c = 0; c < 2; ++c) PetscCall(user->func(dim, 0., &coords[p * dim], 1, &xvals[p * 2 + c], user));
185: }
186: PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
187: PetscCall(DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **)&wvals));
188: PetscCall(DMSwarmRestoreField(*sw, "x_q", NULL, NULL, (void **)&xvals));
190: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: static PetscErrorCode computeParticleMoments(DM sw, Vec u, PetscReal moments[3], AppCtx *user)
195: {
196: DM dm;
197: const PetscReal *coords;
198: const PetscScalar *w;
199: PetscReal mom[3] = {0.0, 0.0, 0.0};
200: PetscInt dim, cStart, cEnd, Nc;
202: PetscFunctionBeginUser;
203: PetscCall(DMGetDimension(sw, &dim));
204: PetscCall(DMSwarmGetCellDM(sw, &dm));
205: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
206: PetscCall(DMSwarmSortGetAccess(sw));
207: PetscCall(DMSwarmGetFieldInfo(sw, user->field ? "x_q" : "w_q", &Nc, NULL));
208: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
209: PetscCall(VecGetArrayRead(u, &w));
210: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
211: PetscInt *pidx;
212: PetscInt Np;
214: PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx));
215: for (PetscInt p = 0; p < Np; ++p) {
216: const PetscInt idx = pidx[p];
217: const PetscReal *x = &coords[idx * dim];
219: for (PetscInt c = 0; c < Nc; ++c) {
220: mom[0] += PetscRealPart(w[idx * Nc + c]);
221: mom[1] += PetscRealPart(w[idx * Nc + c]) * x[0];
222: for (PetscInt d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx * Nc + c]) * PetscSqr(x[d]);
223: }
224: }
225: PetscCall(PetscFree(pidx));
226: }
227: PetscCall(VecRestoreArrayRead(u, &w));
228: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
229: PetscCall(DMSwarmSortRestoreAccess(sw));
230: PetscCallMPI(MPIU_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: static void f0_1(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
235: {
236: const PetscInt Nc = uOff[1] - uOff[0];
238: for (PetscInt c = 0; c < Nc; ++c) f0[0] += u[c];
239: }
241: static void f0_x(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
242: {
243: const PetscInt Nc = uOff[1] - uOff[0];
245: for (PetscInt c = 0; c < Nc; ++c) f0[0] += x[0] * u[c];
246: }
248: static void f0_r2(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
249: {
250: const PetscInt Nc = uOff[1] - uOff[0];
252: for (PetscInt c = 0; c < Nc; ++c)
253: for (PetscInt d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d]) * u[c];
254: }
256: static PetscErrorCode computeFieldMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
257: {
258: PetscDS ds;
259: PetscScalar mom;
261: PetscFunctionBeginUser;
262: PetscCall(DMGetDS(dm, &ds));
263: PetscCall(PetscDSSetObjective(ds, 0, &f0_1));
264: mom = 0.;
265: PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
266: moments[0] = PetscRealPart(mom);
267: PetscCall(PetscDSSetObjective(ds, 0, &f0_x));
268: mom = 0.;
269: PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
270: moments[1] = PetscRealPart(mom);
271: PetscCall(PetscDSSetObjective(ds, 0, &f0_r2));
272: mom = 0.;
273: PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
274: moments[2] = PetscRealPart(mom);
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: static PetscErrorCode TestParticlesToField(DM sw, DM dm, Vec fhat, AppCtx *user)
279: {
280: const char *fieldnames[1] = {user->field ? "x_q" : "w_q"};
281: Vec fields[1] = {fhat}, f;
282: PetscReal pmoments[3]; // \int f, \int x f, \int r^2 f
283: PetscReal fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f
285: PetscFunctionBeginUser;
286: PetscCall(DMSwarmProjectFields(sw, dm, 1, fieldnames, fields, SCATTER_FORWARD));
288: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, fieldnames[0], &f));
289: PetscCall(computeParticleMoments(sw, f, pmoments, user));
290: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, fieldnames[0], &f));
291: PetscCall(computeFieldMoments(dm, fhat, fmoments, user));
292: PetscCall(VecViewFromOptions(fhat, NULL, "-fhat_view"));
293: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
294: for (PetscInt m = 0; m < user->Nm; ++m) {
295: if (user->pass) {
296: if (PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) > user->mtol) {
297: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "p projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", pmoments[0], pmoments[1], pmoments[2]));
298: }
299: } else {
300: PetscCheck(PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) <= user->mtol, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Moment %" PetscInt_FMT " error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]),
301: user->mtol);
302: }
303: }
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: static PetscErrorCode TestFieldToParticles(DM sw, DM dm, Vec fhat, AppCtx *user)
308: {
309: const char *fieldnames[1] = {user->field ? "x_q" : "w_q"};
310: Vec fields[1] = {fhat}, f;
311: PetscReal pmoments[3]; // \int f, \int x f, \int r^2 f
312: PetscReal fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f
314: PetscFunctionBeginUser;
315: PetscCall(DMSwarmProjectFields(sw, dm, 1, fieldnames, fields, SCATTER_REVERSE));
317: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, fieldnames[0], &f));
318: PetscCall(computeParticleMoments(sw, f, pmoments, user));
319: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, fieldnames[0], &f));
320: PetscCall(computeFieldMoments(dm, fhat, fmoments, user));
321: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
322: for (PetscInt m = 0; m < user->Nm; ++m) {
323: if (user->pass) {
324: if (PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) > user->mtol) {
325: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "p projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", pmoments[0], pmoments[1], pmoments[2]));
326: }
327: } else {
328: PetscCheck(PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) <= user->mtol, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Moment %" PetscInt_FMT " error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]),
329: user->mtol);
330: }
331: }
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: int main(int argc, char *argv[])
336: {
337: DM dm, subdm, sw;
338: Vec fhat;
339: IS subis;
340: AppCtx user;
342: PetscFunctionBeginUser;
343: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
344: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
345: PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm, &user));
346: PetscCall(CreateDiscretization(dm, &user));
347: PetscCall(CreateSwarm(dm, &sw, &user));
349: PetscCall(DMCreateSubDM(dm, 1, &user.field, &subis, &subdm));
351: PetscCall(DMGetGlobalVector(subdm, &fhat));
352: PetscCall(PetscObjectSetName((PetscObject)fhat, "FEM f"));
353: PetscCall(TestParticlesToField(sw, subdm, fhat, &user));
354: PetscCall(TestFieldToParticles(sw, subdm, fhat, &user));
355: PetscCall(DMRestoreGlobalVector(subdm, &fhat));
357: PetscCall(ISDestroy(&subis));
358: PetscCall(DMDestroy(&subdm));
359: PetscCall(DMDestroy(&dm));
360: PetscCall(DMDestroy(&sw));
361: PetscCall(PetscFinalize());
362: return 0;
363: }
365: /*TEST
367: # Swarm does not handle complex or quad
368: build:
369: requires: !complex double
371: testset:
372: requires: triangle
373: args: -dm_refine 1 -petscspace_degree 2 -moments 3 \
374: -ptof_pc_type lu \
375: -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
377: test:
378: suffix: tri_fe_f0
379: args: -field 0
381: test:
382: suffix: tri_fe_f1
383: args: -field 1
385: test:
386: suffix: quad_fe_f0
387: args: -dm_plex_simplex 0 -field 0
389: test:
390: suffix: quad_fe_f1
391: args: -dm_plex_simplex 0 -field 1
393: testset:
394: requires: triangle
395: args: -dm_refine 1 -moments 1 -fv \
396: -ptof_pc_type lu \
397: -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
399: test:
400: suffix: tri_fv_f0
401: args: -field 0
403: test:
404: suffix: tri_fv_f1
405: args: -field 1
407: test:
408: suffix: quad_fv_f0
409: args: -dm_plex_simplex 0 -field 0
411: test:
412: suffix: quad_fv_f1
413: args: -dm_plex_simplex 0 -field 1
415: TEST*/