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: {
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: PetscOptionsInt("-N", "Number of particles per spatial cell", "ex27.c", options->N, &options->N, NULL);
43: PetscOptionsReal("-L", "Velocity-space extent", "ex27.c", options->L, &options->L, NULL);
44: PetscOptionsReal("-h", "Velocity-space resolution", "ex27.c", options->h, &options->h, NULL);
45: PetscOptionsReal("-epsilon", "Mollifier regularization parameter", "ex27.c", options->epsilon, &options->epsilon, NULL);
46: PetscOptionsEnd();
47: return 0;
48: }
50: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
51: {
53: DMCreate(comm, dm);
54: DMSetType(*dm, DMPLEX);
55: DMSetFromOptions(*dm);
56: DMViewFromOptions(*dm, NULL, "-dm_view");
57: return 0;
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;
71: /* Randomization for coordinates */
72: PetscRandomCreate(PetscObjectComm((PetscObject)sw), &rnd);
73: PetscRandomSetInterval(rnd, -1.0, 1.0);
74: PetscRandomSetFromOptions(rnd);
75: PetscRandomCreate(PetscObjectComm((PetscObject)sw), &rndv);
76: PetscRandomSetInterval(rndv, -1., 1.);
77: PetscRandomSetFromOptions(rndv);
78: DMGetApplicationContext(sw, &user);
79: Np = user->N;
80: DMGetDimension(sw, &dim);
81: DMSwarmGetCellDM(sw, &dm);
82: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
83: DMPlexGetCellType(dm, cStart, &ct);
84: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
85: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ);
86: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
87: DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
88: DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&velocity);
89: DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&vals);
90: for (c = cStart; c < cEnd; ++c) {
91: if (Np == 1) {
92: DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
93: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
94: vals[c] = 1.0;
95: } else {
96: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
97: for (p = 0; p < Np; ++p) {
98: const PetscInt n = c * Np + p;
99: PetscReal sum = 0.0, refcoords[3];
101: for (d = 0; d < dim; ++d) {
102: PetscRandomGetValueReal(rnd, &refcoords[d]);
103: sum += refcoords[d];
104: }
105: if (simplex && sum > 0.0)
106: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
107: vals[n] = 1.0;
108: DMPlexReferenceToCoordinates(dm, c, 1, refcoords, &coords[n * dim]);
109: }
110: }
111: }
112: /* Random velocity IC */
113: for (c = cStart; c < cEnd; ++c) {
114: for (p = 0; p < Np; ++p) {
115: for (d = 0; d < dim; ++d) {
116: PetscReal v_val;
118: PetscRandomGetValueReal(rndv, &v_val);
119: velocity[p * dim + d] = v_val;
120: }
121: }
122: }
123: DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
124: DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&velocity);
125: DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&vals);
126: PetscFree5(centroid, xi0, v0, J, invJ);
127: PetscRandomDestroy(&rnd);
128: PetscRandomDestroy(&rndv);
129: return 0;
130: }
132: /* Get velocities from swarm and place in solution vector */
133: static PetscErrorCode SetInitialConditions(DM dmSw, Vec u)
134: {
135: DM dm;
136: AppCtx *user;
137: PetscReal *velocity;
138: PetscScalar *initialConditions;
139: PetscInt dim, d, cStart, cEnd, c, Np, p, n;
142: VecGetLocalSize(u, &n);
143: DMGetApplicationContext(dmSw, &user);
144: Np = user->N;
145: DMSwarmGetCellDM(dmSw, &dm);
146: DMGetDimension(dm, &dim);
147: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
148: DMSwarmGetField(dmSw, "velocity", NULL, NULL, (void **)&velocity);
149: VecGetArray(u, &initialConditions);
150: for (c = cStart; c < cEnd; ++c) {
151: for (p = 0; p < Np; ++p) {
152: const PetscInt n = c * Np + p;
153: for (d = 0; d < dim; d++) initialConditions[n * dim + d] = velocity[n * dim + d];
154: }
155: }
156: VecRestoreArray(u, &initialConditions);
157: DMSwarmRestoreField(dmSw, "velocity", NULL, NULL, (void **)&velocity);
158: return 0;
159: }
161: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
162: {
163: PetscInt *cellid;
164: PetscInt dim, cStart, cEnd, c, Np = user->N, p;
165: PetscBool view = PETSC_FALSE;
168: DMGetDimension(dm, &dim);
169: DMCreate(PetscObjectComm((PetscObject)dm), sw);
170: DMSetType(*sw, DMSWARM);
171: DMSetDimension(*sw, dim);
172: /* h = 2L/n and N = n^d */
173: if (user->h < 0.) user->h = 2. * user->L / PetscPowReal(user->N, 1. / dim);
174: /* From Section 4 in [1], \epsilon = 0.64 h^.98 */
175: if (user->epsilon < 0.) user->epsilon = 0.64 * pow(user->h, 1.98);
176: PetscOptionsGetBool(NULL, NULL, "-param_view", &view, NULL);
177: if (view) 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);
178: DMSwarmSetType(*sw, DMSWARM_PIC);
179: DMSwarmSetCellDM(*sw, dm);
180: DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL);
181: DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR);
182: DMSwarmFinalizeFieldRegister(*sw);
183: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
184: DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0);
185: DMSetFromOptions(*sw);
186: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
187: for (c = cStart; c < cEnd; ++c) {
188: for (p = 0; p < Np; ++p) {
189: const PetscInt n = c * Np + p;
190: cellid[n] = c;
191: }
192: }
193: DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
194: PetscObjectSetName((PetscObject)*sw, "Particles");
195: DMViewFromOptions(*sw, NULL, "-sw_view");
196: return 0;
197: }
199: /* Internal dmplex function, same as found in dmpleximpl.h */
200: static void DMPlex_WaxpyD_Internal(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w)
201: {
202: PetscInt d;
204: for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
205: }
207: /* Internal dmplex function, same as found in dmpleximpl.h */
208: static PetscReal DMPlex_DotD_Internal(PetscInt dim, const PetscScalar *x, const PetscReal *y)
209: {
210: PetscReal sum = 0.0;
211: PetscInt d;
213: for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d]) * y[d];
214: return sum;
215: }
217: /* Internal dmplex function, same as found in dmpleximpl.h */
218: static void DMPlex_MultAdd2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
219: {
220: PetscScalar z[2];
221: z[0] = x[0];
222: z[1] = x[ldx];
223: y[0] += A[0] * z[0] + A[1] * z[1];
224: y[ldx] += A[2] * z[0] + A[3] * z[1];
225: (void)PetscLogFlops(6.0);
226: }
228: /* Internal dmplex function, same as found in dmpleximpl.h to avoid private includes. */
229: static void DMPlex_MultAdd3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
230: {
231: PetscScalar z[3];
232: z[0] = x[0];
233: z[1] = x[ldx];
234: z[2] = x[ldx * 2];
235: y[0] += A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
236: y[ldx] += A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
237: y[ldx * 2] += A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
238: (void)PetscLogFlops(15.0);
239: }
241: /*
242: Gaussian - The Gaussian function G(x)
244: Input Parameters:
245: + dim - The number of dimensions, or size of x
246: . mu - The mean, or center
247: . sigma - The standard deviation, or width
248: - x - The evaluation point of the function
250: Output Parameter:
251: . ret - The value G(x)
252: */
253: static PetscReal Gaussian(PetscInt dim, const PetscReal mu[], PetscReal sigma, const PetscReal x[])
254: {
255: PetscReal arg = 0.0;
256: PetscInt d;
258: for (d = 0; d < dim; ++d) arg += PetscSqr(x[d] - mu[d]);
259: return PetscPowReal(2.0 * PETSC_PI * sigma, -dim / 2.0) * PetscExpReal(-arg / (2.0 * sigma));
260: }
262: /*
263: ComputeGradS - Compute grad_v dS_eps/df
265: Input Parameters:
266: + dim - The dimension
267: . Np - The number of particles
268: . vp - The velocity v_p of the particle at which we evaluate
269: . velocity - The velocity field for all particles
270: . epsilon - The regularization strength
272: Output Parameter:
273: . integral - The output grad_v dS_eps/df (v_p)
275: Note:
276: This comes from (3.6) in [1], and we are computing
277: $ \nabla_v S_p = \grad \psi_\epsilon(v_p - v) log \sum_q \psi_\epsilon(v - v_q)
278: which is discretized by using a one-point quadrature in each box l at its center v^c_l
279: $ \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)
280: where h^d is the volume of each box.
281: */
282: static PetscErrorCode ComputeGradS(PetscInt dim, PetscInt Np, const PetscReal vp[], const PetscReal velocity[], PetscReal integral[], AppCtx *ctx)
283: {
284: PetscReal vc_l[3], L = ctx->L, h = ctx->h, epsilon = ctx->epsilon, init = 0.5 * h - L;
285: PetscInt nx = roundf(2. * L / h);
286: PetscInt ny = dim > 1 ? nx : 1;
287: PetscInt nz = dim > 2 ? nx : 1;
288: PetscInt i, j, k, d, q, dbg = 0;
291: for (d = 0; d < dim; ++d) integral[d] = 0.0;
292: for (k = 0, vc_l[2] = init; k < nz; ++k, vc_l[2] += h) {
293: for (j = 0, vc_l[1] = init; j < ny; ++j, vc_l[1] += h) {
294: for (i = 0, vc_l[0] = init; i < nx; ++i, vc_l[0] += h) {
295: PetscReal sum = 0.0;
297: if (dbg) PetscPrintf(PETSC_COMM_SELF, "(%" PetscInt_FMT " %" PetscInt_FMT ") vc_l: %g %g\n", i, j, (double)vc_l[0], (double)vc_l[1]);
298: /* \log \sum_k \psi(v - v_k) */
299: for (q = 0; q < Np; ++q) sum += Gaussian(dim, &velocity[q * dim], epsilon, vc_l);
300: sum = PetscLogReal(sum);
301: for (d = 0; d < dim; ++d) integral[d] += (-1. / (epsilon)) * PetscAbsReal(vp[d] - vc_l[d]) * (Gaussian(dim, vp, epsilon, vc_l)) * sum;
302: }
303: }
304: }
305: return 0;
306: }
308: /* Q = 1/|xi| (I - xi xi^T / |xi|^2), xi = vp - vq */
309: static PetscErrorCode QCompute(PetscInt dim, const PetscReal vp[], const PetscReal vq[], PetscReal Q[])
310: {
311: PetscReal xi[3], xi2, xi3, mag;
312: PetscInt d, e;
315: DMPlex_WaxpyD_Internal(dim, -1.0, vq, vp, xi);
316: xi2 = DMPlex_DotD_Internal(dim, xi, xi);
317: mag = PetscSqrtReal(xi2);
318: xi3 = xi2 * mag;
319: for (d = 0; d < dim; ++d) {
320: for (e = 0; e < dim; ++e) Q[d * dim + e] = -xi[d] * xi[e] / xi3;
321: Q[d * dim + d] += 1. / mag;
322: }
323: return 0;
324: }
326: static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t, Vec U, Vec R, void *ctx)
327: {
328: AppCtx *user = (AppCtx *)ctx;
329: PetscInt dbg = 0;
330: DM sw; /* Particles */
331: Vec sol; /* Solution vector at current time */
332: const PetscScalar *u; /* input solution vector */
333: PetscScalar *r;
334: PetscReal *velocity;
335: PetscInt dim, Np, p, q;
338: VecZeroEntries(R);
339: TSGetDM(ts, &sw);
340: DMGetDimension(sw, &dim);
341: VecGetLocalSize(U, &Np);
342: TSGetSolution(ts, &sol);
343: VecGetArray(sol, &velocity);
344: VecGetArray(R, &r);
345: VecGetArrayRead(U, &u);
346: Np /= dim;
347: if (dbg) PetscPrintf(PETSC_COMM_WORLD, "Part ppr x y\n");
348: for (p = 0; p < Np; ++p) {
349: PetscReal gradS_p[3] = {0., 0., 0.};
351: ComputeGradS(dim, Np, &velocity[p * dim], velocity, gradS_p, user);
352: for (q = 0; q < Np; ++q) {
353: PetscReal gradS_q[3] = {0., 0., 0.}, GammaS[3] = {0., 0., 0.}, Q[9];
355: if (q == p) continue;
356: ComputeGradS(dim, Np, &velocity[q * dim], velocity, gradS_q, user);
357: DMPlex_WaxpyD_Internal(dim, -1.0, gradS_q, gradS_p, GammaS);
358: QCompute(dim, &u[p * dim], &u[q * dim], Q);
359: switch (dim) {
360: case 2:
361: DMPlex_MultAdd2DReal_Internal(Q, 1, GammaS, &r[p * dim]);
362: break;
363: case 3:
364: DMPlex_MultAdd3DReal_Internal(Q, 1, GammaS, &r[p * dim]);
365: break;
366: default:
367: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Do not support dimension %" PetscInt_FMT, dim);
368: }
369: }
370: if (dbg) PetscPrintf(PETSC_COMM_WORLD, "Final %4" PetscInt_FMT " %10.8lf %10.8lf\n", p, r[p * dim + 0], r[p * dim + 1]);
371: }
372: VecRestoreArrayRead(U, &u);
373: VecRestoreArray(R, &r);
374: VecRestoreArray(sol, &velocity);
375: VecViewFromOptions(R, NULL, "-residual_view");
376: return 0;
377: }
379: /*
380: TS Post Step Function. Copy the solution back into the swarm for migration. We may also need to reform
381: the solution vector in cases of particle migration, but we forgo that here since there is no velocity space grid
382: to migrate between.
383: */
384: static PetscErrorCode UpdateSwarm(TS ts)
385: {
386: PetscInt idx, n;
387: const PetscScalar *u;
388: PetscScalar *velocity;
389: DM sw;
390: Vec sol;
393: TSGetDM(ts, &sw);
394: DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&velocity);
395: TSGetSolution(ts, &sol);
396: VecGetArrayRead(sol, &u);
397: VecGetLocalSize(sol, &n);
398: for (idx = 0; idx < n; ++idx) velocity[idx] = u[idx];
399: VecRestoreArrayRead(sol, &u);
400: DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&velocity);
401: return 0;
402: }
404: static PetscErrorCode InitializeSolve(TS ts, Vec u)
405: {
406: DM dm;
407: AppCtx *user;
410: TSGetDM(ts, &dm);
411: DMGetApplicationContext(dm, &user);
412: SetInitialCoordinates(dm);
413: SetInitialConditions(dm, u);
414: return 0;
415: }
417: int main(int argc, char **argv)
418: {
419: TS ts; /* nonlinear solver */
420: DM dm, sw; /* Velocity space mesh and Particle Swarm */
421: Vec u, v; /* problem vector */
422: MPI_Comm comm;
423: AppCtx user;
426: PetscInitialize(&argc, &argv, NULL, help);
427: comm = PETSC_COMM_WORLD;
428: ProcessOptions(comm, &user);
429: /* Initialize objects and set initial conditions */
430: CreateMesh(comm, &dm, &user);
431: CreateParticles(dm, &sw, &user);
432: DMSetApplicationContext(sw, &user);
433: DMSwarmVectorDefineField(sw, "velocity");
434: TSCreate(comm, &ts);
435: TSSetDM(ts, sw);
436: TSSetMaxTime(ts, 10.0);
437: TSSetTimeStep(ts, 0.1);
438: TSSetMaxSteps(ts, 1);
439: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
440: TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user);
441: TSSetFromOptions(ts);
442: TSSetComputeInitialCondition(ts, InitializeSolve);
443: DMSwarmCreateGlobalVectorFromField(sw, "velocity", &v);
444: VecDuplicate(v, &u);
445: DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &v);
446: TSComputeInitialCondition(ts, u);
447: TSSetPostStep(ts, UpdateSwarm);
448: TSSolve(ts, u);
449: VecDestroy(&u);
450: TSDestroy(&ts);
451: DMDestroy(&sw);
452: DMDestroy(&dm);
453: PetscFinalize();
454: return 0;
455: }
457: /*TEST
458: build:
459: requires: triangle !single !complex
460: test:
461: suffix: midpoint
462: 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 \
463: -ts_type theta -ts_theta_theta 0.5 -ts_dmswarm_monitor_moments -ts_monitor_frequency 1 -snes_fd
464: TEST*/