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*/