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