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: PetscInt *cellid;
165: PetscInt dim, cStart, cEnd, c, Np = user->N, p;
166: PetscBool view = PETSC_FALSE;
168: PetscFunctionBeginUser;
169: PetscCall(DMGetDimension(dm, &dim));
170: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
171: PetscCall(DMSetType(*sw, DMSWARM));
172: PetscCall(DMSetDimension(*sw, dim));
173: /* h = 2L/n and N = n^d */
174: if (user->h < 0.) user->h = 2. * user->L / PetscPowReal(user->N, 1. / dim);
175: /* From Section 4 in [1], \epsilon = 0.64 h^.98 */
176: if (user->epsilon < 0.) user->epsilon = 0.64 * pow(user->h, 1.98);
177: PetscCall(PetscOptionsGetBool(NULL, NULL, "-param_view", &view, NULL));
178: 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));
179: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
180: PetscCall(DMSwarmSetCellDM(*sw, dm));
181: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
182: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
183: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
184: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
185: PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0));
186: PetscCall(DMSetFromOptions(*sw));
187: PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
188: for (c = cStart; c < cEnd; ++c) {
189: for (p = 0; p < Np; ++p) {
190: const PetscInt n = c * Np + p;
191: cellid[n] = c;
192: }
193: }
194: PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
195: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
196: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: /* Internal dmplex function, same as found in dmpleximpl.h */
201: static void DMPlex_WaxpyD_Internal(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w)
202: {
203: PetscInt d;
205: for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
206: }
208: /* Internal dmplex function, same as found in dmpleximpl.h */
209: static PetscReal DMPlex_DotD_Internal(PetscInt dim, const PetscScalar *x, const PetscReal *y)
210: {
211: PetscReal sum = 0.0;
212: PetscInt d;
214: for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d]) * y[d];
215: return sum;
216: }
218: /* Internal dmplex function, same as found in dmpleximpl.h */
219: static void DMPlex_MultAdd2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
220: {
221: PetscScalar z[2];
222: z[0] = x[0];
223: z[1] = x[ldx];
224: y[0] += A[0] * z[0] + A[1] * z[1];
225: y[ldx] += A[2] * z[0] + A[3] * z[1];
226: (void)PetscLogFlops(6.0);
227: }
229: /* Internal dmplex function, same as found in dmpleximpl.h to avoid private includes. */
230: static void DMPlex_MultAdd3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
231: {
232: PetscScalar z[3];
233: z[0] = x[0];
234: z[1] = x[ldx];
235: z[2] = x[ldx * 2];
236: y[0] += A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
237: y[ldx] += A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
238: y[ldx * 2] += A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
239: (void)PetscLogFlops(15.0);
240: }
242: /*
243: Gaussian - The Gaussian function G(x)
245: Input Parameters:
246: + dim - The number of dimensions, or size of x
247: . mu - The mean, or center
248: . sigma - The standard deviation, or width
249: - x - The evaluation point of the function
251: Output Parameter:
252: . ret - The value G(x)
253: */
254: static PetscReal Gaussian(PetscInt dim, const PetscReal mu[], PetscReal sigma, const PetscReal x[])
255: {
256: PetscReal arg = 0.0;
257: PetscInt d;
259: for (d = 0; d < dim; ++d) arg += PetscSqr(x[d] - mu[d]);
260: return PetscPowReal(2.0 * PETSC_PI * sigma, -dim / 2.0) * PetscExpReal(-arg / (2.0 * sigma));
261: }
263: /*
264: ComputeGradS - Compute grad_v dS_eps/df
266: Input Parameters:
267: + dim - The dimension
268: . Np - The number of particles
269: . vp - The velocity v_p of the particle at which we evaluate
270: . velocity - The velocity field for all particles
271: . epsilon - The regularization strength
273: Output Parameter:
274: . integral - The output grad_v dS_eps/df (v_p)
276: Note:
277: This comes from (3.6) in [1], and we are computing
278: $ \nabla_v S_p = \grad \psi_\epsilon(v_p - v) log \sum_q \psi_\epsilon(v - v_q)
279: which is discretized by using a one-point quadrature in each box l at its center v^c_l
280: $ \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)
281: where h^d is the volume of each box.
282: */
283: static PetscErrorCode ComputeGradS(PetscInt dim, PetscInt Np, const PetscReal vp[], const PetscReal velocity[], PetscReal integral[], AppCtx *ctx)
284: {
285: PetscReal vc_l[3], L = ctx->L, h = ctx->h, epsilon = ctx->epsilon, init = 0.5 * h - L;
286: PetscInt nx = roundf(2. * L / h);
287: PetscInt ny = dim > 1 ? nx : 1;
288: PetscInt nz = dim > 2 ? nx : 1;
289: PetscInt i, j, k, d, q, dbg = 0;
291: PetscFunctionBeginHot;
292: for (d = 0; d < dim; ++d) integral[d] = 0.0;
293: for (k = 0, vc_l[2] = init; k < nz; ++k, vc_l[2] += h) {
294: for (j = 0, vc_l[1] = init; j < ny; ++j, vc_l[1] += h) {
295: for (i = 0, vc_l[0] = init; i < nx; ++i, vc_l[0] += h) {
296: PetscReal sum = 0.0;
298: 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]));
299: /* \log \sum_k \psi(v - v_k) */
300: for (q = 0; q < Np; ++q) sum += Gaussian(dim, &velocity[q * dim], epsilon, vc_l);
301: sum = PetscLogReal(sum);
302: for (d = 0; d < dim; ++d) integral[d] += (-1. / epsilon) * PetscAbsReal(vp[d] - vc_l[d]) * Gaussian(dim, vp, epsilon, vc_l) * sum;
303: }
304: }
305: }
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /* Q = 1/|xi| (I - xi xi^T / |xi|^2), xi = vp - vq */
310: static PetscErrorCode QCompute(PetscInt dim, const PetscReal vp[], const PetscReal vq[], PetscReal Q[])
311: {
312: PetscReal xi[3], xi2, xi3, mag;
313: PetscInt d, e;
315: PetscFunctionBeginHot;
316: DMPlex_WaxpyD_Internal(dim, -1.0, vq, vp, xi);
317: xi2 = DMPlex_DotD_Internal(dim, xi, xi);
318: mag = PetscSqrtReal(xi2);
319: xi3 = xi2 * mag;
320: for (d = 0; d < dim; ++d) {
321: for (e = 0; e < dim; ++e) Q[d * dim + e] = -xi[d] * xi[e] / xi3;
322: Q[d * dim + d] += 1. / mag;
323: }
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t, Vec U, Vec R, void *ctx)
328: {
329: AppCtx *user = (AppCtx *)ctx;
330: PetscInt dbg = 0;
331: DM sw; /* Particles */
332: Vec sol; /* Solution vector at current time */
333: const PetscScalar *u; /* input solution vector */
334: PetscScalar *r;
335: PetscReal *velocity;
336: PetscInt dim, Np, p, q;
338: PetscFunctionBeginUser;
339: PetscCall(VecZeroEntries(R));
340: PetscCall(TSGetDM(ts, &sw));
341: PetscCall(DMGetDimension(sw, &dim));
342: PetscCall(VecGetLocalSize(U, &Np));
343: PetscCall(TSGetSolution(ts, &sol));
344: PetscCall(VecGetArray(sol, &velocity));
345: PetscCall(VecGetArray(R, &r));
346: PetscCall(VecGetArrayRead(U, &u));
347: Np /= dim;
348: if (dbg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Part ppr x y\n"));
349: for (p = 0; p < Np; ++p) {
350: PetscReal gradS_p[3] = {0., 0., 0.};
352: PetscCall(ComputeGradS(dim, Np, &velocity[p * dim], velocity, gradS_p, user));
353: for (q = 0; q < Np; ++q) {
354: PetscReal gradS_q[3] = {0., 0., 0.}, GammaS[3] = {0., 0., 0.}, Q[9];
356: if (q == p) continue;
357: PetscCall(ComputeGradS(dim, Np, &velocity[q * dim], velocity, gradS_q, user));
358: DMPlex_WaxpyD_Internal(dim, -1.0, gradS_q, gradS_p, GammaS);
359: PetscCall(QCompute(dim, &u[p * dim], &u[q * dim], Q));
360: switch (dim) {
361: case 2:
362: DMPlex_MultAdd2DReal_Internal(Q, 1, GammaS, &r[p * dim]);
363: break;
364: case 3:
365: DMPlex_MultAdd3DReal_Internal(Q, 1, GammaS, &r[p * dim]);
366: break;
367: default:
368: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Do not support dimension %" PetscInt_FMT, dim);
369: }
370: }
371: 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]));
372: }
373: PetscCall(VecRestoreArrayRead(U, &u));
374: PetscCall(VecRestoreArray(R, &r));
375: PetscCall(VecRestoreArray(sol, &velocity));
376: PetscCall(VecViewFromOptions(R, NULL, "-residual_view"));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: /*
381: TS Post Step Function. Copy the solution back into the swarm for migration. We may also need to reform
382: the solution vector in cases of particle migration, but we forgo that here since there is no velocity space grid
383: to migrate between.
384: */
385: static PetscErrorCode UpdateSwarm(TS ts)
386: {
387: PetscInt idx, n;
388: const PetscScalar *u;
389: PetscScalar *velocity;
390: DM sw;
391: Vec sol;
393: PetscFunctionBeginUser;
394: PetscCall(TSGetDM(ts, &sw));
395: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&velocity));
396: PetscCall(TSGetSolution(ts, &sol));
397: PetscCall(VecGetArrayRead(sol, &u));
398: PetscCall(VecGetLocalSize(sol, &n));
399: for (idx = 0; idx < n; ++idx) velocity[idx] = u[idx];
400: PetscCall(VecRestoreArrayRead(sol, &u));
401: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&velocity));
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: static PetscErrorCode InitializeSolve(TS ts, Vec u)
406: {
407: DM dm;
408: AppCtx *user;
410: PetscFunctionBeginUser;
411: PetscCall(TSGetDM(ts, &dm));
412: PetscCall(DMGetApplicationContext(dm, &user));
413: PetscCall(SetInitialCoordinates(dm));
414: PetscCall(SetInitialConditions(dm, u));
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: int main(int argc, char **argv)
419: {
420: TS ts; /* nonlinear solver */
421: DM dm, sw; /* Velocity space mesh and Particle Swarm */
422: Vec u, v; /* problem vector */
423: MPI_Comm comm;
424: AppCtx user;
426: PetscFunctionBeginUser;
427: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
428: comm = PETSC_COMM_WORLD;
429: PetscCall(ProcessOptions(comm, &user));
430: /* Initialize objects and set initial conditions */
431: PetscCall(CreateMesh(comm, &dm, &user));
432: PetscCall(CreateParticles(dm, &sw, &user));
433: PetscCall(DMSetApplicationContext(sw, &user));
434: PetscCall(DMSwarmVectorDefineField(sw, "velocity"));
435: PetscCall(TSCreate(comm, &ts));
436: PetscCall(TSSetDM(ts, sw));
437: PetscCall(TSSetMaxTime(ts, 10.0));
438: PetscCall(TSSetTimeStep(ts, 0.1));
439: PetscCall(TSSetMaxSteps(ts, 1));
440: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
441: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user));
442: PetscCall(TSSetFromOptions(ts));
443: PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
444: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &v));
445: PetscCall(VecDuplicate(v, &u));
446: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &v));
447: PetscCall(TSComputeInitialCondition(ts, u));
448: PetscCall(TSSetPostStep(ts, UpdateSwarm));
449: PetscCall(TSSolve(ts, u));
450: PetscCall(VecDestroy(&u));
451: PetscCall(TSDestroy(&ts));
452: PetscCall(DMDestroy(&sw));
453: PetscCall(DMDestroy(&dm));
454: PetscCall(PetscFinalize());
455: return 0;
456: }
458: /*TEST
459: build:
460: requires: triangle !single !complex
461: test:
462: suffix: midpoint
463: 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 \
464: -ts_type theta -ts_theta_theta 0.5 -ts_dmswarm_monitor_moments -ts_monitor_frequency 1 -snes_fd
465: TEST*/