Actual source code: ex27.c
1: static char help[] = "Particle Basis Landau Example using nonlinear solve + Implicit Midpoint-like time stepping.";
3: /* TODO
5: 1) SNES is sensitive to epsilon (but not to h). Should we do continuation in it?
7: 2) Put this timestepper in library, maybe by changing DG
9: 3) Add monitor to visualize distributions
11: */
13: /* References
14: [1] https://arxiv.org/abs/1910.03080v2
15: */
17: #include <petscdmplex.h>
18: #include <petscdmswarm.h>
19: #include <petscts.h>
20: #include <petscviewer.h>
21: #include <petscmath.h>
23: typedef struct {
24: /* Velocity space grid and functions */
25: PetscInt N; /* The number of partices per spatial cell */
26: PetscReal L; /* Velocity space is [-L, L]^d */
27: PetscReal h; /* Spacing for grid 2L / N^{1/d} */
28: PetscReal epsilon; /* gaussian regularization parameter */
29: PetscReal momentTol; /* Tolerance for checking moment conservation */
30: } AppCtx;
32: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
33: {
34: PetscFunctionBeginUser;
35: options->N = 1;
36: options->momentTol = 100.0 * PETSC_MACHINE_EPSILON;
37: options->L = 1.0;
38: options->h = -1.0;
39: options->epsilon = -1.0;
41: PetscOptionsBegin(comm, "", "Collision Options", "DMPLEX");
42: PetscCall(PetscOptionsInt("-N", "Number of particles per spatial cell", "ex27.c", options->N, &options->N, NULL));
43: PetscCall(PetscOptionsReal("-L", "Velocity-space extent", "ex27.c", options->L, &options->L, NULL));
44: PetscCall(PetscOptionsReal("-h", "Velocity-space resolution", "ex27.c", options->h, &options->h, NULL));
45: PetscCall(PetscOptionsReal("-epsilon", "Mollifier regularization parameter", "ex27.c", options->epsilon, &options->epsilon, NULL));
46: PetscOptionsEnd();
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
51: {
52: PetscFunctionBeginUser;
53: PetscCall(DMCreate(comm, dm));
54: PetscCall(DMSetType(*dm, DMPLEX));
55: PetscCall(DMSetFromOptions(*dm));
56: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: static PetscErrorCode SetInitialCoordinates(DM sw)
61: {
62: AppCtx *user;
63: PetscRandom rnd, rndv;
64: DM dm;
65: DMPolytopeType ct;
66: PetscBool simplex;
67: PetscReal *centroid, *coords, *velocity, *xi0, *v0, *J, *invJ, detJ, *vals;
68: PetscInt dim, d, cStart, cEnd, c, Np, p;
70: PetscFunctionBeginUser;
71: /* Randomization for coordinates */
72: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)sw), &rnd));
73: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
74: PetscCall(PetscRandomSetFromOptions(rnd));
75: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)sw), &rndv));
76: PetscCall(PetscRandomSetInterval(rndv, -1., 1.));
77: PetscCall(PetscRandomSetFromOptions(rndv));
78: PetscCall(DMGetApplicationContext(sw, &user));
79: Np = user->N;
80: PetscCall(DMGetDimension(sw, &dim));
81: PetscCall(DMSwarmGetCellDM(sw, &dm));
82: PetscCall(DMGetCoordinatesLocalSetUp(dm));
83: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
84: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
85: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
86: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
87: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
88: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
89: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&velocity));
90: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&vals));
91: for (c = cStart; c < cEnd; ++c) {
92: if (Np == 1) {
93: PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL));
94: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
95: vals[c] = 1.0;
96: } else {
97: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); /* affine */
98: for (p = 0; p < Np; ++p) {
99: const PetscInt n = c * Np + p;
100: PetscReal sum = 0.0, refcoords[3];
102: for (d = 0; d < dim; ++d) {
103: PetscCall(PetscRandomGetValueReal(rnd, &refcoords[d]));
104: sum += refcoords[d];
105: }
106: if (simplex && sum > 0.0)
107: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
108: vals[n] = 1.0;
109: PetscCall(DMPlexReferenceToCoordinates(dm, c, 1, refcoords, &coords[n * dim]));
110: }
111: }
112: }
113: /* Random velocity IC */
114: for (c = cStart; c < cEnd; ++c) {
115: for (p = 0; p < Np; ++p) {
116: for (d = 0; d < dim; ++d) {
117: PetscReal v_val;
119: PetscCall(PetscRandomGetValueReal(rndv, &v_val));
120: velocity[p * dim + d] = v_val;
121: }
122: }
123: }
124: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
125: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&velocity));
126: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&vals));
127: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
128: PetscCall(PetscRandomDestroy(&rnd));
129: PetscCall(PetscRandomDestroy(&rndv));
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: /* Get velocities from swarm and place in solution vector */
134: static PetscErrorCode SetInitialConditions(DM dmSw, Vec u)
135: {
136: DM dm;
137: AppCtx *user;
138: PetscReal *velocity;
139: PetscScalar *initialConditions;
140: PetscInt dim, d, cStart, cEnd, c, Np, p, n;
142: PetscFunctionBeginUser;
143: PetscCall(VecGetLocalSize(u, &n));
144: PetscCall(DMGetApplicationContext(dmSw, &user));
145: Np = user->N;
146: PetscCall(DMSwarmGetCellDM(dmSw, &dm));
147: PetscCall(DMGetDimension(dm, &dim));
148: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
149: PetscCall(DMSwarmGetField(dmSw, "velocity", NULL, NULL, (void **)&velocity));
150: PetscCall(VecGetArray(u, &initialConditions));
151: for (c = cStart; c < cEnd; ++c) {
152: for (p = 0; p < Np; ++p) {
153: const PetscInt n = c * Np + p;
154: for (d = 0; d < dim; d++) initialConditions[n * dim + d] = velocity[n * dim + d];
155: }
156: }
157: PetscCall(VecRestoreArray(u, &initialConditions));
158: PetscCall(DMSwarmRestoreField(dmSw, "velocity", NULL, NULL, (void **)&velocity));
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
163: {
164: DMSwarmCellDM celldm;
165: PetscInt *swarm_cellid;
166: PetscInt dim, cStart, cEnd, c, Np = user->N, p;
167: PetscBool view = PETSC_FALSE;
168: const char *cellid;
170: PetscFunctionBeginUser;
171: PetscCall(DMGetDimension(dm, &dim));
172: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
173: PetscCall(DMSetType(*sw, DMSWARM));
174: PetscCall(DMSetDimension(*sw, dim));
175: /* h = 2L/n and N = n^d */
176: if (user->h < 0.) user->h = 2. * user->L / PetscPowReal(user->N, 1. / dim);
177: /* From Section 4 in [1], \epsilon = 0.64 h^.98 */
178: if (user->epsilon < 0.) user->epsilon = 0.64 * pow(user->h, 1.98);
179: PetscCall(PetscOptionsGetBool(NULL, NULL, "-param_view", &view, NULL));
180: if (view) PetscCall(PetscPrintf(PETSC_COMM_SELF, "N: %" PetscInt_FMT " L: %g h: %g eps: %g\n", user->N, (double)user->L, (double)user->h, (double)user->epsilon));
181: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
182: PetscCall(DMSwarmSetCellDM(*sw, dm));
183: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
184: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
185: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
186: PetscCall(DMSwarmGetCellDMActive(*sw, &celldm));
187: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
188: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
189: PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0));
190: PetscCall(DMSetFromOptions(*sw));
191: PetscCall(DMSwarmGetField(*sw, cellid, NULL, NULL, (void **)&swarm_cellid));
192: for (c = cStart; c < cEnd; ++c) {
193: for (p = 0; p < Np; ++p) {
194: const PetscInt n = c * Np + p;
195: swarm_cellid[n] = c;
196: }
197: }
198: PetscCall(DMSwarmRestoreField(*sw, cellid, NULL, NULL, (void **)&swarm_cellid));
199: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
200: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /* Internal dmplex function, same as found in dmpleximpl.h */
205: static void DMPlex_WaxpyD_Internal(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w)
206: {
207: PetscInt d;
209: for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
210: }
212: /* Internal dmplex function, same as found in dmpleximpl.h */
213: static PetscReal DMPlex_DotD_Internal(PetscInt dim, const PetscScalar *x, const PetscReal *y)
214: {
215: PetscReal sum = 0.0;
216: PetscInt d;
218: for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d]) * y[d];
219: return sum;
220: }
222: /* Internal dmplex function, same as found in dmpleximpl.h */
223: static void DMPlex_MultAdd2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
224: {
225: PetscScalar z[2];
226: z[0] = x[0];
227: z[1] = x[ldx];
228: y[0] += A[0] * z[0] + A[1] * z[1];
229: y[ldx] += A[2] * z[0] + A[3] * z[1];
230: (void)PetscLogFlops(6.0);
231: }
233: /* Internal dmplex function, same as found in dmpleximpl.h to avoid private includes. */
234: static void DMPlex_MultAdd3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
235: {
236: PetscScalar z[3];
237: z[0] = x[0];
238: z[1] = x[ldx];
239: z[2] = x[ldx * 2];
240: y[0] += A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
241: y[ldx] += A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
242: y[ldx * 2] += A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
243: (void)PetscLogFlops(15.0);
244: }
246: /*
247: Gaussian - The Gaussian function G(x)
249: Input Parameters:
250: + dim - The number of dimensions, or size of x
251: . mu - The mean, or center
252: . sigma - The standard deviation, or width
253: - x - The evaluation point of the function
255: Output Parameter:
256: . ret - The value G(x)
257: */
258: static PetscReal Gaussian(PetscInt dim, const PetscReal mu[], PetscReal sigma, const PetscReal x[])
259: {
260: PetscReal arg = 0.0;
261: PetscInt d;
263: for (d = 0; d < dim; ++d) arg += PetscSqr(x[d] - mu[d]);
264: return PetscPowReal(2.0 * PETSC_PI * sigma, -dim / 2.0) * PetscExpReal(-arg / (2.0 * sigma));
265: }
267: /*
268: ComputeGradS - Compute grad_v dS_eps/df
270: Input Parameters:
271: + dim - The dimension
272: . Np - The number of particles
273: . vp - The velocity v_p of the particle at which we evaluate
274: . velocity - The velocity field for all particles
275: . epsilon - The regularization strength
277: Output Parameter:
278: . integral - The output grad_v dS_eps/df (v_p)
280: Note:
281: This comes from (3.6) in [1], and we are computing
282: $ \nabla_v S_p = \grad \psi_\epsilon(v_p - v) log \sum_q \psi_\epsilon(v - v_q)
283: which is discretized by using a one-point quadrature in each box l at its center v^c_l
284: $ \sum_l h^d \nabla\psi_\epsilon(v_p - v^c_l) \log\left( \sum_q w_q \psi_\epsilon(v^c_l - v_q) \right)
285: where h^d is the volume of each box.
286: */
287: static PetscErrorCode ComputeGradS(PetscInt dim, PetscInt Np, const PetscReal vp[], const PetscReal velocity[], PetscReal integral[], AppCtx *ctx)
288: {
289: PetscReal vc_l[3], L = ctx->L, h = ctx->h, epsilon = ctx->epsilon, init = 0.5 * h - L;
290: PetscInt nx = roundf(2. * L / h);
291: PetscInt ny = dim > 1 ? nx : 1;
292: PetscInt nz = dim > 2 ? nx : 1;
293: PetscInt i, j, k, d, q, dbg = 0;
295: PetscFunctionBeginHot;
296: for (d = 0; d < dim; ++d) integral[d] = 0.0;
297: for (k = 0, vc_l[2] = init; k < nz; ++k, vc_l[2] += h) {
298: for (j = 0, vc_l[1] = init; j < ny; ++j, vc_l[1] += h) {
299: for (i = 0, vc_l[0] = init; i < nx; ++i, vc_l[0] += h) {
300: PetscReal sum = 0.0;
302: if (dbg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "(%" PetscInt_FMT " %" PetscInt_FMT ") vc_l: %g %g\n", i, j, (double)vc_l[0], (double)vc_l[1]));
303: /* \log \sum_k \psi(v - v_k) */
304: for (q = 0; q < Np; ++q) sum += Gaussian(dim, &velocity[q * dim], epsilon, vc_l);
305: sum = PetscLogReal(sum);
306: for (d = 0; d < dim; ++d) integral[d] += (-1. / epsilon) * PetscAbsReal(vp[d] - vc_l[d]) * Gaussian(dim, vp, epsilon, vc_l) * sum;
307: }
308: }
309: }
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: /* Q = 1/|xi| (I - xi xi^T / |xi|^2), xi = vp - vq */
314: static PetscErrorCode QCompute(PetscInt dim, const PetscReal vp[], const PetscReal vq[], PetscReal Q[])
315: {
316: PetscReal xi[3], xi2, xi3, mag;
317: PetscInt d, e;
319: PetscFunctionBeginHot;
320: DMPlex_WaxpyD_Internal(dim, -1.0, vq, vp, xi);
321: xi2 = DMPlex_DotD_Internal(dim, xi, xi);
322: mag = PetscSqrtReal(xi2);
323: xi3 = xi2 * mag;
324: for (d = 0; d < dim; ++d) {
325: for (e = 0; e < dim; ++e) Q[d * dim + e] = -xi[d] * xi[e] / xi3;
326: Q[d * dim + d] += 1. / mag;
327: }
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t, Vec U, Vec R, void *ctx)
332: {
333: AppCtx *user = (AppCtx *)ctx;
334: PetscInt dbg = 0;
335: DM sw; /* Particles */
336: Vec sol; /* Solution vector at current time */
337: const PetscScalar *u; /* input solution vector */
338: PetscScalar *r;
339: PetscReal *velocity;
340: PetscInt dim, Np, p, q;
342: PetscFunctionBeginUser;
343: PetscCall(VecZeroEntries(R));
344: PetscCall(TSGetDM(ts, &sw));
345: PetscCall(DMGetDimension(sw, &dim));
346: PetscCall(VecGetLocalSize(U, &Np));
347: PetscCall(TSGetSolution(ts, &sol));
348: PetscCall(VecGetArray(sol, &velocity));
349: PetscCall(VecGetArray(R, &r));
350: PetscCall(VecGetArrayRead(U, &u));
351: Np /= dim;
352: if (dbg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Part ppr x y\n"));
353: for (p = 0; p < Np; ++p) {
354: PetscReal gradS_p[3] = {0., 0., 0.};
356: PetscCall(ComputeGradS(dim, Np, &velocity[p * dim], velocity, gradS_p, user));
357: for (q = 0; q < Np; ++q) {
358: PetscReal gradS_q[3] = {0., 0., 0.}, GammaS[3] = {0., 0., 0.}, Q[9];
360: if (q == p) continue;
361: PetscCall(ComputeGradS(dim, Np, &velocity[q * dim], velocity, gradS_q, user));
362: DMPlex_WaxpyD_Internal(dim, -1.0, gradS_q, gradS_p, GammaS);
363: PetscCall(QCompute(dim, &u[p * dim], &u[q * dim], Q));
364: switch (dim) {
365: case 2:
366: DMPlex_MultAdd2DReal_Internal(Q, 1, GammaS, &r[p * dim]);
367: break;
368: case 3:
369: DMPlex_MultAdd3DReal_Internal(Q, 1, GammaS, &r[p * dim]);
370: break;
371: default:
372: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Do not support dimension %" PetscInt_FMT, dim);
373: }
374: }
375: if (dbg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Final %4" PetscInt_FMT " %10.8lf %10.8lf\n", p, r[p * dim + 0], r[p * dim + 1]));
376: }
377: PetscCall(VecRestoreArrayRead(U, &u));
378: PetscCall(VecRestoreArray(R, &r));
379: PetscCall(VecRestoreArray(sol, &velocity));
380: PetscCall(VecViewFromOptions(R, NULL, "-residual_view"));
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: /*
385: TS Post Step Function. Copy the solution back into the swarm for migration. We may also need to reform
386: the solution vector in cases of particle migration, but we forgo that here since there is no velocity space grid
387: to migrate between.
388: */
389: static PetscErrorCode UpdateSwarm(TS ts)
390: {
391: PetscInt idx, n;
392: const PetscScalar *u;
393: PetscScalar *velocity;
394: DM sw;
395: Vec sol;
397: PetscFunctionBeginUser;
398: PetscCall(TSGetDM(ts, &sw));
399: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&velocity));
400: PetscCall(TSGetSolution(ts, &sol));
401: PetscCall(VecGetArrayRead(sol, &u));
402: PetscCall(VecGetLocalSize(sol, &n));
403: for (idx = 0; idx < n; ++idx) velocity[idx] = u[idx];
404: PetscCall(VecRestoreArrayRead(sol, &u));
405: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&velocity));
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: static PetscErrorCode InitializeSolve(TS ts, Vec u)
410: {
411: DM dm;
412: AppCtx *user;
414: PetscFunctionBeginUser;
415: PetscCall(TSGetDM(ts, &dm));
416: PetscCall(DMGetApplicationContext(dm, &user));
417: PetscCall(SetInitialCoordinates(dm));
418: PetscCall(SetInitialConditions(dm, u));
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: int main(int argc, char **argv)
423: {
424: TS ts; /* nonlinear solver */
425: DM dm, sw; /* Velocity space mesh and Particle Swarm */
426: Vec u, v; /* problem vector */
427: MPI_Comm comm;
428: AppCtx user;
430: PetscFunctionBeginUser;
431: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
432: comm = PETSC_COMM_WORLD;
433: PetscCall(ProcessOptions(comm, &user));
434: /* Initialize objects and set initial conditions */
435: PetscCall(CreateMesh(comm, &dm, &user));
436: PetscCall(CreateParticles(dm, &sw, &user));
437: PetscCall(DMSetApplicationContext(sw, &user));
438: PetscCall(DMSwarmVectorDefineField(sw, "velocity"));
439: PetscCall(TSCreate(comm, &ts));
440: PetscCall(TSSetDM(ts, sw));
441: PetscCall(TSSetMaxTime(ts, 10.0));
442: PetscCall(TSSetTimeStep(ts, 0.1));
443: PetscCall(TSSetMaxSteps(ts, 1));
444: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
445: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user));
446: PetscCall(TSSetFromOptions(ts));
447: PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
448: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &v));
449: PetscCall(VecDuplicate(v, &u));
450: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &v));
451: PetscCall(TSComputeInitialCondition(ts, u));
452: PetscCall(TSSetPostStep(ts, UpdateSwarm));
453: PetscCall(TSSolve(ts, u));
454: PetscCall(VecDestroy(&u));
455: PetscCall(TSDestroy(&ts));
456: PetscCall(DMDestroy(&sw));
457: PetscCall(DMDestroy(&dm));
458: PetscCall(PetscFinalize());
459: return 0;
460: }
462: /*TEST
463: build:
464: requires: triangle !single !complex
465: test:
466: suffix: midpoint
467: args: -N 3 -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1,-1 -dm_plex_box_upper 1,1 -dm_view \
468: -ts_type theta -ts_theta_theta 0.5 -ts_dmswarm_monitor_moments -ts_monitor_frequency 1 -snes_fd
469: TEST*/