Actual source code: ex28.c
1: static char help[] = "Example application of the Bhatnagar-Gross-Krook (BGK) collision operator.\n\
2: This example is a 0D-1V setting for the kinetic equation\n\
3: https://en.wikipedia.org/wiki/Bhatnagar%E2%80%93Gross%E2%80%93Krook_operator\n";
5: #include <petscdmplex.h>
6: #include <petscdmswarm.h>
7: #include <petscts.h>
8: #include <petscdraw.h>
9: #include <petscviewer.h>
11: typedef struct {
12: PetscInt particlesPerCell; /* The number of partices per cell */
13: PetscReal momentTol; /* Tolerance for checking moment conservation */
14: PetscBool monitorhg; /* Flag for using the TS histogram monitor */
15: PetscBool monitorsp; /* Flag for using the TS scatter monitor */
16: PetscBool monitorks; /* Monitor to perform KS test to the maxwellian */
17: PetscBool error; /* Flag for printing the error */
18: PetscInt ostep; /* print the energy at each ostep time steps */
19: PetscDraw draw; /* The draw object for histogram monitoring */
20: PetscDrawHG drawhg; /* The histogram draw context for monitoring */
21: PetscDrawSP drawsp; /* The scatter plot draw context for the monitor */
22: PetscDrawSP drawks; /* Scatterplot draw context for KS test */
23: } AppCtx;
25: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
26: {
27: PetscFunctionBeginUser;
28: options->monitorhg = PETSC_FALSE;
29: options->monitorsp = PETSC_FALSE;
30: options->monitorks = PETSC_FALSE;
31: options->particlesPerCell = 1;
32: options->momentTol = 100.0 * PETSC_MACHINE_EPSILON;
33: options->ostep = 100;
35: PetscOptionsBegin(comm, "", "Collision Options", "DMPLEX");
36: PetscCall(PetscOptionsBool("-monitorhg", "Flag to use the TS histogram monitor", "ex28.c", options->monitorhg, &options->monitorhg, NULL));
37: PetscCall(PetscOptionsBool("-monitorsp", "Flag to use the TS scatter plot monitor", "ex28.c", options->monitorsp, &options->monitorsp, NULL));
38: PetscCall(PetscOptionsBool("-monitorks", "Flag to plot KS test results", "ex28.c", options->monitorks, &options->monitorks, NULL));
39: PetscCall(PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex28.c", options->particlesPerCell, &options->particlesPerCell, NULL));
40: PetscCall(PetscOptionsInt("-output_step", "Number of time steps between output", "ex28.c", options->ostep, &options->ostep, NULL));
41: PetscOptionsEnd();
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: /* Create the mesh for velocity space */
46: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
47: {
48: PetscFunctionBeginUser;
49: PetscCall(DMCreate(comm, dm));
50: PetscCall(DMSetType(*dm, DMPLEX));
51: PetscCall(DMSetFromOptions(*dm));
52: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
53: PetscFunctionReturn(PETSC_SUCCESS);
54: }
56: /* Since we are putting the same number of particles in each cell, this amounts to a uniform distribution of v */
57: static PetscErrorCode SetInitialCoordinates(DM sw)
58: {
59: AppCtx *user;
60: PetscRandom rnd;
61: DM dm;
62: DMPolytopeType ct;
63: PetscBool simplex;
64: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ, *vals;
65: PetscInt dim, d, cStart, cEnd, c, Np, p;
67: PetscFunctionBeginUser;
68: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)sw), &rnd));
69: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
70: PetscCall(PetscRandomSetFromOptions(rnd));
72: PetscCall(DMGetApplicationContext(sw, &user));
73: Np = user->particlesPerCell;
74: PetscCall(DMGetDimension(sw, &dim));
75: PetscCall(DMSwarmGetCellDM(sw, &dm));
76: PetscCall(DMGetCoordinatesLocalSetUp(dm));
77: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
78: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
79: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
80: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
81: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
82: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
83: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&vals));
84: for (c = cStart; c < cEnd; ++c) {
85: if (Np == 1) {
86: PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL));
87: for (d = 0; d < dim; ++d) {
88: coords[c * dim + d] = centroid[d];
89: if ((coords[c * dim + d] >= -1) && (coords[c * dim + d] <= 1)) {
90: vals[c] = 1.0;
91: } else {
92: vals[c] = 0.;
93: }
94: }
95: } else {
96: PetscCall(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: PetscCall(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: PetscCall(DMPlexReferenceToCoordinates(dm, c, 1, refcoords, &coords[n * dim]));
109: }
110: }
111: }
112: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
113: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&vals));
114: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
115: PetscCall(PetscRandomDestroy(&rnd));
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: /* The initial conditions are just the initial particle weights */
120: static PetscErrorCode SetInitialConditions(DM dmSw, Vec u)
121: {
122: DM dm;
123: AppCtx *user;
124: PetscReal *vals;
125: PetscScalar *initialConditions;
126: PetscInt dim, d, cStart, cEnd, c, Np, p, n;
128: PetscFunctionBeginUser;
129: PetscCall(VecGetLocalSize(u, &n));
130: PetscCall(DMGetApplicationContext(dmSw, &user));
131: Np = user->particlesPerCell;
132: PetscCall(DMSwarmGetCellDM(dmSw, &dm));
133: PetscCall(DMGetDimension(dm, &dim));
134: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
135: PetscCheck(n == (cEnd - cStart) * Np, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "TS solution local size %" PetscInt_FMT " != %" PetscInt_FMT " nm particles", n, (cEnd - cStart) * Np);
136: PetscCall(DMSwarmGetField(dmSw, "w_q", NULL, NULL, (void **)&vals));
137: PetscCall(VecGetArray(u, &initialConditions));
138: for (c = cStart; c < cEnd; ++c) {
139: for (p = 0; p < Np; ++p) {
140: const PetscInt n = c * Np + p;
141: for (d = 0; d < dim; d++) initialConditions[n] = vals[n];
142: }
143: }
144: PetscCall(VecRestoreArray(u, &initialConditions));
145: PetscCall(DMSwarmRestoreField(dmSw, "w_q", NULL, NULL, (void **)&vals));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
150: {
151: DMSwarmCellDM celldm;
152: PetscInt *swarm_cellid;
153: PetscInt dim, cStart, cEnd, c, Np = user->particlesPerCell, p;
154: const char *cellid;
156: PetscFunctionBeginUser;
157: PetscCall(DMGetDimension(dm, &dim));
158: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
159: PetscCall(DMSetType(*sw, DMSWARM));
160: PetscCall(DMSetDimension(*sw, dim));
161: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
162: PetscCall(DMSwarmSetCellDM(*sw, dm));
163: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", dim, PETSC_REAL));
164: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
165: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
166: PetscCall(DMSwarmGetCellDMActive(*sw, &celldm));
167: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
168: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
169: PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0));
170: PetscCall(DMSetFromOptions(*sw));
171: PetscCall(DMSwarmGetField(*sw, cellid, NULL, NULL, (void **)&swarm_cellid));
172: for (c = cStart; c < cEnd; ++c) {
173: for (p = 0; p < Np; ++p) {
174: const PetscInt n = c * Np + p;
175: swarm_cellid[n] = c;
176: }
177: }
178: PetscCall(DMSwarmRestoreField(*sw, cellid, NULL, NULL, (void **)&swarm_cellid));
179: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
180: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: /*
185: f_t = 1/\tau \left( f_eq - f \right)
186: n_t = 1/\tau \left( \int f_eq - \int f \right) = 1/\tau (n - n) = 0
187: v_t = 1/\tau \left( \int v f_eq - \int v f \right) = 1/\tau (v - v) = 0
188: E_t = 1/\tau \left( \int v^2 f_eq - \int v^2 f \right) = 1/\tau (T - T) = 0
190: Let's look at a single cell:
192: \int_C f_t = 1/\tau \left( \int_C f_eq - \int_C f \right)
193: \sum_{x_i \in C} w^i_t = 1/\tau \left( F_eq(C) - \sum_{x_i \in C} w_i \right)
194: */
196: /* This computes the 1D Maxwellian distribution for given mass n, velocity v, and temperature T */
197: static PetscReal ComputePDF(PetscReal m, PetscReal n, PetscReal T, PetscReal v[])
198: {
199: return (n / PetscSqrtReal(2.0 * PETSC_PI * T / m)) * PetscExpReal(-0.5 * m * PetscSqr(v[0]) / T);
200: }
202: /*
203: erf z = \frac{2}{\sqrt\pi} \int^z_0 dt e^{-t^2} and erf \infty = 1
205: We begin with our distribution
207: \sqrt{\frac{m}{2 \pi T}} e^{-m v^2/2T}
209: Let t = \sqrt{m/2T} v, z = \sqrt{m/2T} w, so that we now have
211: \sqrt{\frac{m}{2 \pi T}} \int^w_0 dv e^{-m v^2/2T}
212: = \sqrt{\frac{m}{2 \pi T}} \int^{\sqrt{m/2T} w}_0 \sqrt{2T/m} dt e^{-t^2}
213: = 1/\sqrt{\pi} \int^{\sqrt{m/2T} w}_0 dt e^{-t^2}
214: = 1/2 erf(\sqrt{m/2T} w)
215: */
216: static PetscReal ComputeCDF(PetscReal m, PetscReal n, PetscReal T, PetscReal va, PetscReal vb)
217: {
218: PetscReal alpha = PetscSqrtReal(0.5 * m / T);
219: PetscReal za = alpha * va;
220: PetscReal zb = alpha * vb;
221: PetscReal sum = 0.0;
223: sum += zb >= 0 ? erf(zb) : -erf(-zb);
224: sum -= za >= 0 ? erf(za) : -erf(-za);
225: return 0.5 * n * sum;
226: }
228: static PetscErrorCode CheckDistribution(DM dm, PetscReal m, PetscReal n, PetscReal T, PetscReal v[])
229: {
230: PetscSection coordSection;
231: Vec coordsLocal;
232: PetscReal *xq, *wq;
233: PetscReal vmin, vmax, neq, veq, Teq;
234: PetscInt Nq = 100, q, cStart, cEnd, c;
236: PetscFunctionBeginUser;
237: PetscCall(DMGetBoundingBox(dm, &vmin, &vmax));
238: /* Check analytic over entire line */
239: neq = ComputeCDF(m, n, T, vmin, vmax);
240: PetscCheck(PetscAbsReal(neq - n) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Int f %g != %g mass (%g)", (double)neq, (double)n, (double)(neq - n));
241: /* Check analytic over cells */
242: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
243: PetscCall(DMGetCoordinateSection(dm, &coordSection));
244: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
245: neq = 0.0;
246: for (c = cStart; c < cEnd; ++c) {
247: PetscScalar *vcoords = NULL;
249: PetscCall(DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &vcoords));
250: neq += ComputeCDF(m, n, T, vcoords[0], vcoords[1]);
251: PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &vcoords));
252: }
253: PetscCheck(PetscAbsReal(neq - n) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell Int f %g != %g mass (%g)", (double)neq, (double)n, (double)(neq - n));
254: /* Check quadrature over entire line */
255: PetscCall(PetscMalloc2(Nq, &xq, Nq, &wq));
256: PetscCall(PetscDTGaussQuadrature(100, vmin, vmax, xq, wq));
257: neq = 0.0;
258: for (q = 0; q < Nq; ++q) neq += ComputePDF(m, n, T, &xq[q]) * wq[q];
259: PetscCheck(PetscAbsReal(neq - n) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Int f %g != %g mass (%g)", (double)neq, (double)n, (double)(neq - n));
260: /* Check omemnts with quadrature */
261: veq = 0.0;
262: for (q = 0; q < Nq; ++q) veq += xq[q] * ComputePDF(m, n, T, &xq[q]) * wq[q];
263: veq /= neq;
264: PetscCheck(PetscAbsReal(veq - v[0]) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Int v f %g != %g velocity (%g)", (double)veq, (double)v[0], (double)(veq - v[0]));
265: Teq = 0.0;
266: for (q = 0; q < Nq; ++q) Teq += PetscSqr(xq[q]) * ComputePDF(m, n, T, &xq[q]) * wq[q];
267: Teq = Teq * m / neq - PetscSqr(veq);
268: PetscCheck(PetscAbsReal(Teq - T) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Int v^2 f %g != %g temperature (%g)", (double)Teq, (double)T, (double)(Teq - T));
269: PetscCall(PetscFree2(xq, wq));
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t, Vec U, Vec R, void *ctx)
274: {
275: const PetscScalar *u;
276: PetscSection coordSection;
277: Vec coordsLocal;
278: PetscScalar *r, *coords;
279: PetscReal n = 0.0, v = 0.0, E = 0.0, T = 0.0, m = 1.0, cn = 0.0, cv = 0.0, cE = 0.0, pE = 0.0, eqE = 0.0;
280: PetscInt dim, d, Np, Ncp, p, cStart, cEnd, c;
281: DM dmSw, plex;
283: PetscFunctionBeginUser;
284: PetscCall(VecGetLocalSize(U, &Np));
285: PetscCall(VecGetArrayRead(U, &u));
286: PetscCall(VecGetArray(R, &r));
287: PetscCall(TSGetDM(ts, &dmSw));
288: PetscCall(DMSwarmGetCellDM(dmSw, &plex));
289: PetscCall(DMGetDimension(dmSw, &dim));
290: PetscCall(DMGetCoordinatesLocal(plex, &coordsLocal));
291: PetscCall(DMGetCoordinateSection(plex, &coordSection));
292: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
293: Np /= dim;
294: Ncp = Np / (cEnd - cStart);
295: /* Calculate moments of particle distribution, note that velocity is in the coordinate */
296: PetscCall(DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
297: for (p = 0; p < Np; ++p) {
298: PetscReal m1 = 0.0, m2 = 0.0;
300: for (d = 0; d < dim; ++d) {
301: m1 += PetscRealPart(coords[p * dim + d]);
302: m2 += PetscSqr(PetscRealPart(coords[p * dim + d]));
303: }
304: n += u[p];
305: v += u[p] * m1;
306: E += u[p] * m2;
307: }
308: v /= n;
309: T = E * m / n - v * v;
310: PetscCall(PetscInfo(ts, "Time %.2f: mass %.4f velocity: %+.4f temperature: %.4f\n", (double)t, (double)n, (double)v, (double)T));
311: PetscCall(CheckDistribution(plex, m, n, T, &v));
312: /*
313: Begin cellwise evaluation of the collision operator. Essentially, penalize the weights of the particles
314: in that cell against the maxwellian for the number of particles expected to be in that cell
315: */
316: for (c = cStart; c < cEnd; ++c) {
317: PetscScalar *vcoords = NULL;
318: PetscReal relaxation = 1.0, neq;
319: PetscInt sp = c * Ncp, q;
321: /* Calculate equilibrium occupation for this velocity cell */
322: PetscCall(DMPlexVecGetClosure(plex, coordSection, coordsLocal, c, NULL, &vcoords));
323: neq = ComputeCDF(m, n, T, vcoords[0], vcoords[1]);
324: PetscCall(DMPlexVecRestoreClosure(plex, coordSection, coordsLocal, c, NULL, &vcoords));
325: for (q = 0; q < Ncp; ++q) r[sp + q] = (1.0 / relaxation) * (neq - u[sp + q]);
326: }
327: /* Check update */
328: for (p = 0; p < Np; ++p) {
329: PetscReal m1 = 0.0, m2 = 0.0;
330: PetscScalar *vcoords = NULL;
332: for (d = 0; d < dim; ++d) {
333: m1 += PetscRealPart(coords[p * dim + d]);
334: m2 += PetscSqr(PetscRealPart(coords[p * dim + d]));
335: }
336: cn += r[p];
337: cv += r[p] * m1;
338: cE += r[p] * m2;
339: pE += u[p] * m2;
340: PetscCall(DMPlexVecGetClosure(plex, coordSection, coordsLocal, p, NULL, &vcoords));
341: eqE += ComputeCDF(m, n, T, vcoords[0], vcoords[1]) * m2;
342: PetscCall(DMPlexVecRestoreClosure(plex, coordSection, coordsLocal, p, NULL, &vcoords));
343: }
344: PetscCall(PetscInfo(ts, "Time %.2f: mass update %.8f velocity update: %+.8f energy update: %.8f (%.8f, %.8f)\n", (double)t, (double)cn, (double)cv, (double)cE, (double)pE, (double)eqE));
345: PetscCall(DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
346: PetscCall(VecRestoreArrayRead(U, &u));
347: PetscCall(VecRestoreArray(R, &r));
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: static PetscErrorCode HGMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
352: {
353: AppCtx *user = (AppCtx *)ctx;
354: const PetscScalar *u;
355: DM sw, dm;
356: PetscInt dim, Np, p;
358: PetscFunctionBeginUser;
359: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS);
360: if ((user->ostep > 0) && (!(step % user->ostep))) {
361: PetscDrawAxis axis;
363: PetscCall(TSGetDM(ts, &sw));
364: PetscCall(DMSwarmGetCellDM(sw, &dm));
365: PetscCall(DMGetDimension(dm, &dim));
366: PetscCall(PetscDrawHGReset(user->drawhg));
367: PetscCall(PetscDrawHGGetAxis(user->drawhg, &axis));
368: PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "V", "f(V)"));
369: PetscCall(PetscDrawAxisSetLimits(axis, -3, 3, 0, 100));
370: PetscCall(PetscDrawHGSetLimits(user->drawhg, -3.0, 3.0, 0, 10));
371: PetscCall(VecGetLocalSize(U, &Np));
372: Np /= dim;
373: PetscCall(VecGetArrayRead(U, &u));
374: /* get points from solution vector */
375: for (p = 0; p < Np; ++p) PetscCall(PetscDrawHGAddValue(user->drawhg, u[p]));
376: PetscCall(PetscDrawHGDraw(user->drawhg));
377: PetscCall(VecRestoreArrayRead(U, &u));
378: }
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: static PetscErrorCode SPMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
383: {
384: AppCtx *user = (AppCtx *)ctx;
385: const PetscScalar *u;
386: PetscReal *v, *coords;
387: PetscInt Np, p;
388: DM dmSw;
390: PetscFunctionBeginUser;
391: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS);
392: if ((user->ostep > 0) && (!(step % user->ostep))) {
393: PetscDrawAxis axis;
395: PetscCall(TSGetDM(ts, &dmSw));
396: PetscCall(PetscDrawSPReset(user->drawsp));
397: PetscCall(PetscDrawSPGetAxis(user->drawsp, &axis));
398: PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "V", "w"));
399: PetscCall(VecGetLocalSize(U, &Np));
400: PetscCall(VecGetArrayRead(U, &u));
401: /* get points from solution vector */
402: PetscCall(PetscMalloc1(Np, &v));
403: for (p = 0; p < Np; ++p) v[p] = PetscRealPart(u[p]);
404: PetscCall(DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
405: for (p = 0; p < Np - 1; ++p) PetscCall(PetscDrawSPAddPoint(user->drawsp, &coords[p], &v[p]));
406: PetscCall(PetscDrawSPDraw(user->drawsp, PETSC_TRUE));
407: PetscCall(VecRestoreArrayRead(U, &u));
408: PetscCall(DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
409: PetscCall(PetscFree(v));
410: }
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: static PetscErrorCode KSConv(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
415: {
416: AppCtx *user = (AppCtx *)ctx;
417: const PetscScalar *u;
418: PetscScalar sup;
419: PetscReal *v, *coords, T = 0., vel = 0., step_cast, w_sum;
420: PetscInt dim, Np, p, cStart, cEnd;
421: DM sw, plex;
423: PetscFunctionBeginUser;
424: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS);
425: if ((user->ostep > 0) && (!(step % user->ostep))) {
426: PetscDrawAxis axis;
427: PetscCall(PetscDrawSPGetAxis(user->drawks, &axis));
428: PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "t", "D_n"));
429: PetscCall(PetscDrawSPSetLimits(user->drawks, 0., 100, 0., 3.5));
430: PetscCall(TSGetDM(ts, &sw));
431: PetscCall(VecGetLocalSize(U, &Np));
432: PetscCall(VecGetArrayRead(U, &u));
433: /* get points from solution vector */
434: PetscCall(PetscMalloc1(Np, &v));
435: PetscCall(DMSwarmGetCellDM(sw, &plex));
436: PetscCall(DMGetDimension(sw, &dim));
437: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
438: for (p = 0; p < Np; ++p) v[p] = PetscRealPart(u[p]);
439: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
440: w_sum = 0.;
441: for (p = 0; p < Np; ++p) {
442: w_sum += u[p];
443: T += u[p] * coords[p] * coords[p];
444: vel += u[p] * coords[p];
445: }
446: vel /= w_sum;
447: T = T / w_sum - vel * vel;
448: sup = 0.0;
449: for (p = 0; p < Np; ++p) {
450: PetscReal tmp = 0.;
451: tmp = PetscAbs(u[p] - ComputePDF(1.0, w_sum, T, &coords[p * dim]));
452: if (tmp > sup) sup = tmp;
453: }
454: step_cast = (PetscReal)step;
455: PetscCall(PetscDrawSPAddPoint(user->drawks, &step_cast, &sup));
456: PetscCall(PetscDrawSPDraw(user->drawks, PETSC_TRUE));
457: PetscCall(VecRestoreArrayRead(U, &u));
458: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
459: PetscCall(PetscFree(v));
460: }
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: static PetscErrorCode InitializeSolve(TS ts, Vec u)
465: {
466: DM dm;
467: AppCtx *user;
469: PetscFunctionBeginUser;
470: PetscCall(TSGetDM(ts, &dm));
471: PetscCall(DMGetApplicationContext(dm, &user));
472: PetscCall(SetInitialCoordinates(dm));
473: PetscCall(SetInitialConditions(dm, u));
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
476: /*
477: A single particle is generated for each velocity space cell using the dmswarmpicfield_coor and is used to evaluate collisions in that cell.
478: 0 weight ghost particles are initialized outside of a small velocity domain to ensure the tails of the amxwellian are resolved.
479: */
480: int main(int argc, char **argv)
481: {
482: TS ts; /* nonlinear solver */
483: DM dm, sw; /* Velocity space mesh and Particle Swarm */
484: Vec u, w; /* swarm vector */
485: MPI_Comm comm;
486: AppCtx user;
488: PetscFunctionBeginUser;
489: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
490: comm = PETSC_COMM_WORLD;
491: PetscCall(ProcessOptions(comm, &user));
492: PetscCall(CreateMesh(comm, &dm, &user));
493: PetscCall(CreateParticles(dm, &sw, &user));
494: PetscCall(DMSetApplicationContext(sw, &user));
495: PetscCall(TSCreate(comm, &ts));
496: PetscCall(TSSetDM(ts, sw));
497: PetscCall(TSSetMaxTime(ts, 10.0));
498: PetscCall(TSSetTimeStep(ts, 0.01));
499: PetscCall(TSSetMaxSteps(ts, 100000));
500: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
501: if (user.monitorhg) {
502: PetscCall(PetscDrawCreate(comm, NULL, "monitor", 0, 0, 400, 300, &user.draw));
503: PetscCall(PetscDrawSetFromOptions(user.draw));
504: PetscCall(PetscDrawHGCreate(user.draw, 20, &user.drawhg));
505: PetscCall(PetscDrawHGSetColor(user.drawhg, 3));
506: PetscCall(TSMonitorSet(ts, HGMonitor, &user, NULL));
507: } else if (user.monitorsp) {
508: PetscCall(PetscDrawCreate(comm, NULL, "monitor", 0, 0, 400, 300, &user.draw));
509: PetscCall(PetscDrawSetFromOptions(user.draw));
510: PetscCall(PetscDrawSPCreate(user.draw, 1, &user.drawsp));
511: PetscCall(TSMonitorSet(ts, SPMonitor, &user, NULL));
512: } else if (user.monitorks) {
513: PetscCall(PetscDrawCreate(comm, NULL, "monitor", 0, 0, 400, 300, &user.draw));
514: PetscCall(PetscDrawSetFromOptions(user.draw));
515: PetscCall(PetscDrawSPCreate(user.draw, 1, &user.drawks));
516: PetscCall(TSMonitorSet(ts, KSConv, &user, NULL));
517: }
518: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user));
519: PetscCall(TSSetFromOptions(ts));
520: PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
521: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &w));
522: PetscCall(VecDuplicate(w, &u));
523: PetscCall(VecCopy(w, u));
524: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &w));
525: PetscCall(TSComputeInitialCondition(ts, u));
526: PetscCall(TSSolve(ts, u));
527: if (user.monitorhg) {
528: PetscCall(PetscDrawSave(user.draw));
529: PetscCall(PetscDrawHGDestroy(&user.drawhg));
530: PetscCall(PetscDrawDestroy(&user.draw));
531: }
532: if (user.monitorsp) {
533: PetscCall(PetscDrawSave(user.draw));
534: PetscCall(PetscDrawSPDestroy(&user.drawsp));
535: PetscCall(PetscDrawDestroy(&user.draw));
536: }
537: if (user.monitorks) {
538: PetscCall(PetscDrawSave(user.draw));
539: PetscCall(PetscDrawSPDestroy(&user.drawks));
540: PetscCall(PetscDrawDestroy(&user.draw));
541: }
542: PetscCall(VecDestroy(&u));
543: PetscCall(TSDestroy(&ts));
544: PetscCall(DMDestroy(&sw));
545: PetscCall(DMDestroy(&dm));
546: PetscCall(PetscFinalize());
547: return 0;
548: }
550: /*TEST
551: build:
552: requires: double !complex
553: test:
554: suffix: 1
555: args: -particles_per_cell 1 -output_step 10 -ts_type euler -dm_plex_dim 1 -dm_plex_box_faces 200 -dm_plex_box_lower -10 -dm_plex_box_upper 10 -dm_view -monitorsp
556: test:
557: suffix: 2
558: args: -particles_per_cell 1 -output_step 50 -ts_type euler -dm_plex_dim 1 -dm_plex_box_faces 200 -dm_plex_box_lower -10 -dm_plex_box_upper 10 -dm_view -monitorks
559: TEST*/