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: {
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: PetscOptionsBool("-monitorhg", "Flag to use the TS histogram monitor", "ex28.c", options->monitorhg, &options->monitorhg, NULL);
37: PetscOptionsBool("-monitorsp", "Flag to use the TS scatter plot monitor", "ex28.c", options->monitorsp, &options->monitorsp, NULL);
38: PetscOptionsBool("-monitorks", "Flag to plot KS test results", "ex28.c", options->monitorks, &options->monitorks, NULL);
39: PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex28.c", options->particlesPerCell, &options->particlesPerCell, NULL);
40: PetscOptionsInt("-output_step", "Number of time steps between output", "ex28.c", options->ostep, &options->ostep, PETSC_NULL);
41: PetscOptionsEnd();
43: return 0;
44: }
46: /* Create the mesh for velocity space */
47: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
48: {
50: DMCreate(comm, dm);
51: DMSetType(*dm, DMPLEX);
52: DMSetFromOptions(*dm);
53: DMViewFromOptions(*dm, NULL, "-dm_view");
54: return 0;
55: }
57: /* Since we are putting the same number of particles in each cell, this amounts to a uniform distribution of v */
58: static PetscErrorCode SetInitialCoordinates(DM sw)
59: {
60: AppCtx *user;
61: PetscRandom rnd;
62: DM dm;
63: DMPolytopeType ct;
64: PetscBool simplex;
65: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ, *vals;
66: PetscInt dim, d, cStart, cEnd, c, Np, p;
69: PetscRandomCreate(PetscObjectComm((PetscObject)sw), &rnd);
70: PetscRandomSetInterval(rnd, -1.0, 1.0);
71: PetscRandomSetFromOptions(rnd);
73: DMGetApplicationContext(sw, &user);
74: Np = user->particlesPerCell;
75: DMGetDimension(sw, &dim);
76: DMSwarmGetCellDM(sw, &dm);
77: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
78: DMPlexGetCellType(dm, cStart, &ct);
79: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
80: 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: DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
83: DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&vals);
84: for (c = cStart; c < cEnd; ++c) {
85: if (Np == 1) {
86: 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: 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: DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
113: DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&vals);
114: PetscFree5(centroid, xi0, v0, J, invJ);
115: PetscRandomDestroy(&rnd);
116: return 0;
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;
129: VecGetLocalSize(u, &n);
130: DMGetApplicationContext(dmSw, &user);
131: Np = user->particlesPerCell;
132: DMSwarmGetCellDM(dmSw, &dm);
133: DMGetDimension(dm, &dim);
134: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
136: DMSwarmGetField(dmSw, "w_q", NULL, NULL, (void **)&vals);
137: 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: VecRestoreArray(u, &initialConditions);
145: DMSwarmRestoreField(dmSw, "w_q", NULL, NULL, (void **)&vals);
146: return 0;
147: }
149: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
150: {
151: PetscInt *cellid;
152: PetscInt dim, cStart, cEnd, c, Np = user->particlesPerCell, p;
155: DMGetDimension(dm, &dim);
156: DMCreate(PetscObjectComm((PetscObject)dm), sw);
157: DMSetType(*sw, DMSWARM);
158: DMSetDimension(*sw, dim);
159: DMSwarmSetType(*sw, DMSWARM_PIC);
160: DMSwarmSetCellDM(*sw, dm);
161: DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", dim, PETSC_REAL);
162: DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR);
163: DMSwarmFinalizeFieldRegister(*sw);
164: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
165: DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0);
166: DMSetFromOptions(*sw);
167: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
168: for (c = cStart; c < cEnd; ++c) {
169: for (p = 0; p < Np; ++p) {
170: const PetscInt n = c * Np + p;
171: cellid[n] = c;
172: }
173: }
174: DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
175: PetscObjectSetName((PetscObject)*sw, "Particles");
176: DMViewFromOptions(*sw, NULL, "-sw_view");
177: return 0;
178: }
180: /*
181: f_t = 1/\tau \left( f_eq - f \right)
182: n_t = 1/\tau \left( \int f_eq - \int f \right) = 1/\tau (n - n) = 0
183: v_t = 1/\tau \left( \int v f_eq - \int v f \right) = 1/\tau (v - v) = 0
184: E_t = 1/\tau \left( \int v^2 f_eq - \int v^2 f \right) = 1/\tau (T - T) = 0
186: Let's look at a single cell:
188: \int_C f_t = 1/\tau \left( \int_C f_eq - \int_C f \right)
189: \sum_{x_i \in C} w^i_t = 1/\tau \left( F_eq(C) - \sum_{x_i \in C} w_i \right)
190: */
192: /* This computes the 1D Maxwellian distribution for given mass n, velocity v, and temperature T */
193: static PetscReal ComputePDF(PetscReal m, PetscReal n, PetscReal T, PetscReal v[])
194: {
195: return (n / PetscSqrtReal(2.0 * PETSC_PI * T / m)) * PetscExpReal(-0.5 * m * PetscSqr(v[0]) / T);
196: }
198: /*
199: erf z = \frac{2}{\sqrt\pi} \int^z_0 dt e^{-t^2} and erf \infty = 1
201: We begin with our distribution
203: \sqrt{\frac{m}{2 \pi T}} e^{-m v^2/2T}
205: Let t = \sqrt{m/2T} v, z = \sqrt{m/2T} w, so that we now have
207: \sqrt{\frac{m}{2 \pi T}} \int^w_0 dv e^{-m v^2/2T}
208: = \sqrt{\frac{m}{2 \pi T}} \int^{\sqrt{m/2T} w}_0 \sqrt{2T/m} dt e^{-t^2}
209: = 1/\sqrt{\pi} \int^{\sqrt{m/2T} w}_0 dt e^{-t^2}
210: = 1/2 erf(\sqrt{m/2T} w)
211: */
212: static PetscReal ComputeCDF(PetscReal m, PetscReal n, PetscReal T, PetscReal va, PetscReal vb)
213: {
214: PetscReal alpha = PetscSqrtReal(0.5 * m / T);
215: PetscReal za = alpha * va;
216: PetscReal zb = alpha * vb;
217: PetscReal sum = 0.0;
219: sum += zb >= 0 ? erf(zb) : -erf(-zb);
220: sum -= za >= 0 ? erf(za) : -erf(-za);
221: return 0.5 * n * sum;
222: }
224: static PetscErrorCode CheckDistribution(DM dm, PetscReal m, PetscReal n, PetscReal T, PetscReal v[])
225: {
226: PetscSection coordSection;
227: Vec coordsLocal;
228: PetscReal *xq, *wq;
229: PetscReal vmin, vmax, neq, veq, Teq;
230: PetscInt Nq = 100, q, cStart, cEnd, c;
233: DMGetBoundingBox(dm, &vmin, &vmax);
234: /* Check analytic over entire line */
235: neq = ComputeCDF(m, n, T, vmin, vmax);
237: /* Check analytic over cells */
238: DMGetCoordinatesLocal(dm, &coordsLocal);
239: DMGetCoordinateSection(dm, &coordSection);
240: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
241: neq = 0.0;
242: for (c = cStart; c < cEnd; ++c) {
243: PetscScalar *vcoords = NULL;
245: DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &vcoords);
246: neq += ComputeCDF(m, n, T, vcoords[0], vcoords[1]);
247: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &vcoords);
248: }
250: /* Check quadrature over entire line */
251: PetscMalloc2(Nq, &xq, Nq, &wq);
252: PetscDTGaussQuadrature(100, vmin, vmax, xq, wq);
253: neq = 0.0;
254: for (q = 0; q < Nq; ++q) neq += ComputePDF(m, n, T, &xq[q]) * wq[q];
256: /* Check omemnts with quadrature */
257: veq = 0.0;
258: for (q = 0; q < Nq; ++q) veq += xq[q] * ComputePDF(m, n, T, &xq[q]) * wq[q];
259: veq /= neq;
261: Teq = 0.0;
262: for (q = 0; q < Nq; ++q) Teq += PetscSqr(xq[q]) * ComputePDF(m, n, T, &xq[q]) * wq[q];
263: Teq = Teq * m / neq - PetscSqr(veq);
265: PetscFree2(xq, wq);
266: return 0;
267: }
269: static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t, Vec U, Vec R, void *ctx)
270: {
271: const PetscScalar *u;
272: PetscSection coordSection;
273: Vec coordsLocal;
274: PetscScalar *r, *coords;
275: 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;
276: PetscInt dim, d, Np, Ncp, p, cStart, cEnd, c;
277: DM dmSw, plex;
280: VecGetLocalSize(U, &Np);
281: VecGetArrayRead(U, &u);
282: VecGetArray(R, &r);
283: TSGetDM(ts, &dmSw);
284: DMSwarmGetCellDM(dmSw, &plex);
285: DMGetDimension(dmSw, &dim);
286: DMGetCoordinatesLocal(plex, &coordsLocal);
287: DMGetCoordinateSection(plex, &coordSection);
288: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
289: Np /= dim;
290: Ncp = Np / (cEnd - cStart);
291: /* Calculate moments of particle distribution, note that velocity is in the coordinate */
292: DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
293: for (p = 0; p < Np; ++p) {
294: PetscReal m1 = 0.0, m2 = 0.0;
296: for (d = 0; d < dim; ++d) {
297: m1 += PetscRealPart(coords[p * dim + d]);
298: m2 += PetscSqr(PetscRealPart(coords[p * dim + d]));
299: }
300: n += u[p];
301: v += u[p] * m1;
302: E += u[p] * m2;
303: }
304: v /= n;
305: T = E * m / n - v * v;
306: PetscInfo(ts, "Time %.2f: mass %.4f velocity: %+.4f temperature: %.4f\n", (double)t, (double)n, (double)v, (double)T);
307: CheckDistribution(plex, m, n, T, &v);
308: /*
309: Begin cellwise evaluation of the collision operator. Essentially, penalize the weights of the particles
310: in that cell against the maxwellian for the number of particles expected to be in that cell
311: */
312: for (c = cStart; c < cEnd; ++c) {
313: PetscScalar *vcoords = NULL;
314: PetscReal relaxation = 1.0, neq;
315: PetscInt sp = c * Ncp, q;
317: /* Calculate equilibrium occupation for this velocity cell */
318: DMPlexVecGetClosure(plex, coordSection, coordsLocal, c, NULL, &vcoords);
319: neq = ComputeCDF(m, n, T, vcoords[0], vcoords[1]);
320: DMPlexVecRestoreClosure(plex, coordSection, coordsLocal, c, NULL, &vcoords);
321: for (q = 0; q < Ncp; ++q) r[sp + q] = (1.0 / relaxation) * (neq - u[sp + q]);
322: }
323: /* Check update */
324: for (p = 0; p < Np; ++p) {
325: PetscReal m1 = 0.0, m2 = 0.0;
326: PetscScalar *vcoords = NULL;
328: for (d = 0; d < dim; ++d) {
329: m1 += PetscRealPart(coords[p * dim + d]);
330: m2 += PetscSqr(PetscRealPart(coords[p * dim + d]));
331: }
332: cn += r[p];
333: cv += r[p] * m1;
334: cE += r[p] * m2;
335: pE += u[p] * m2;
336: DMPlexVecGetClosure(plex, coordSection, coordsLocal, p, NULL, &vcoords);
337: eqE += ComputeCDF(m, n, T, vcoords[0], vcoords[1]) * m2;
338: DMPlexVecRestoreClosure(plex, coordSection, coordsLocal, p, NULL, &vcoords);
339: }
340: 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);
341: DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
342: VecRestoreArrayRead(U, &u);
343: VecRestoreArray(R, &r);
344: return 0;
345: }
347: static PetscErrorCode HGMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
348: {
349: AppCtx *user = (AppCtx *)ctx;
350: const PetscScalar *u;
351: DM sw, dm;
352: PetscInt dim, Np, p;
355: if (step < 0) return 0;
356: if (((user->ostep > 0) && (!(step % user->ostep)))) {
357: PetscDrawAxis axis;
359: TSGetDM(ts, &sw);
360: DMSwarmGetCellDM(sw, &dm);
361: DMGetDimension(dm, &dim);
362: PetscDrawHGReset(user->drawhg);
363: PetscDrawHGGetAxis(user->drawhg, &axis);
364: PetscDrawAxisSetLabels(axis, "Particles", "V", "f(V)");
365: PetscDrawAxisSetLimits(axis, -3, 3, 0, 100);
366: PetscDrawHGSetLimits(user->drawhg, -3.0, 3.0, 0, 10);
367: VecGetLocalSize(U, &Np);
368: Np /= dim;
369: VecGetArrayRead(U, &u);
370: /* get points from solution vector */
371: for (p = 0; p < Np; ++p) PetscDrawHGAddValue(user->drawhg, u[p]);
372: PetscDrawHGDraw(user->drawhg);
373: VecRestoreArrayRead(U, &u);
374: }
375: return 0;
376: }
378: static PetscErrorCode SPMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
379: {
380: AppCtx *user = (AppCtx *)ctx;
381: const PetscScalar *u;
382: PetscReal *v, *coords;
383: PetscInt Np, p;
384: DM dmSw;
388: if (step < 0) return 0;
389: if (((user->ostep > 0) && (!(step % user->ostep)))) {
390: PetscDrawAxis axis;
392: TSGetDM(ts, &dmSw);
393: PetscDrawSPReset(user->drawsp);
394: PetscDrawSPGetAxis(user->drawsp, &axis);
395: PetscDrawAxisSetLabels(axis, "Particles", "V", "w");
396: VecGetLocalSize(U, &Np);
397: VecGetArrayRead(U, &u);
398: /* get points from solution vector */
399: PetscMalloc1(Np, &v);
400: for (p = 0; p < Np; ++p) v[p] = PetscRealPart(u[p]);
401: DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
402: for (p = 0; p < Np - 1; ++p) PetscDrawSPAddPoint(user->drawsp, &coords[p], &v[p]);
403: PetscDrawSPDraw(user->drawsp, PETSC_TRUE);
404: VecRestoreArrayRead(U, &u);
405: DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
406: PetscFree(v);
407: }
408: return 0;
409: }
411: static PetscErrorCode KSConv(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
412: {
413: AppCtx *user = (AppCtx *)ctx;
414: const PetscScalar *u;
415: PetscScalar sup;
416: PetscReal *v, *coords, T = 0., vel = 0., step_cast, w_sum;
417: PetscInt dim, Np, p, cStart, cEnd;
418: DM sw, plex;
421: if (step < 0) return 0;
422: if (((user->ostep > 0) && (!(step % user->ostep)))) {
423: PetscDrawAxis axis;
424: PetscDrawSPGetAxis(user->drawks, &axis);
425: PetscDrawAxisSetLabels(axis, "Particles", "t", "D_n");
426: PetscDrawSPSetLimits(user->drawks, 0., 100, 0., 3.5);
427: TSGetDM(ts, &sw);
428: VecGetLocalSize(U, &Np);
429: VecGetArrayRead(U, &u);
430: /* get points from solution vector */
431: PetscMalloc1(Np, &v);
432: DMSwarmGetCellDM(sw, &plex);
433: DMGetDimension(sw, &dim);
434: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
435: for (p = 0; p < Np; ++p) v[p] = PetscRealPart(u[p]);
436: DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
437: w_sum = 0.;
438: for (p = 0; p < Np; ++p) {
439: w_sum += u[p];
440: T += u[p] * coords[p] * coords[p];
441: vel += u[p] * coords[p];
442: }
443: vel /= w_sum;
444: T = T / w_sum - vel * vel;
445: sup = 0.0;
446: for (p = 0; p < Np; ++p) {
447: PetscReal tmp = 0.;
448: tmp = PetscAbs(u[p] - ComputePDF(1.0, w_sum, T, &coords[p * dim]));
449: if (tmp > sup) sup = tmp;
450: }
451: step_cast = (PetscReal)step;
452: PetscDrawSPAddPoint(user->drawks, &step_cast, &sup);
453: PetscDrawSPDraw(user->drawks, PETSC_TRUE);
454: VecRestoreArrayRead(U, &u);
455: DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
456: PetscFree(v);
457: }
458: return 0;
459: }
461: static PetscErrorCode InitializeSolve(TS ts, Vec u)
462: {
463: DM dm;
464: AppCtx *user;
467: TSGetDM(ts, &dm);
468: DMGetApplicationContext(dm, &user);
469: SetInitialCoordinates(dm);
470: SetInitialConditions(dm, u);
471: return 0;
472: }
473: /*
474: A single particle is generated for each velocity space cell using the dmswarmpicfield_coor and is used to evaluate collisions in that cell.
475: 0 weight ghost particles are initialized outside of a small velocity domain to ensure the tails of the amxwellian are resolved.
476: */
477: int main(int argc, char **argv)
478: {
479: TS ts; /* nonlinear solver */
480: DM dm, sw; /* Velocity space mesh and Particle Swarm */
481: Vec u, w; /* swarm vector */
482: MPI_Comm comm;
483: AppCtx user;
486: PetscInitialize(&argc, &argv, NULL, help);
487: comm = PETSC_COMM_WORLD;
488: ProcessOptions(comm, &user);
489: CreateMesh(comm, &dm, &user);
490: CreateParticles(dm, &sw, &user);
491: DMSetApplicationContext(sw, &user);
492: TSCreate(comm, &ts);
493: TSSetDM(ts, sw);
494: TSSetMaxTime(ts, 10.0);
495: TSSetTimeStep(ts, 0.01);
496: TSSetMaxSteps(ts, 100000);
497: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
498: if (user.monitorhg) {
499: PetscDrawCreate(comm, NULL, "monitor", 0, 0, 400, 300, &user.draw);
500: PetscDrawSetFromOptions(user.draw);
501: PetscDrawHGCreate(user.draw, 20, &user.drawhg);
502: PetscDrawHGSetColor(user.drawhg, 3);
503: TSMonitorSet(ts, HGMonitor, &user, NULL);
504: } else if (user.monitorsp) {
505: PetscDrawCreate(comm, NULL, "monitor", 0, 0, 400, 300, &user.draw);
506: PetscDrawSetFromOptions(user.draw);
507: PetscDrawSPCreate(user.draw, 1, &user.drawsp);
508: TSMonitorSet(ts, SPMonitor, &user, NULL);
509: } else if (user.monitorks) {
510: PetscDrawCreate(comm, NULL, "monitor", 0, 0, 400, 300, &user.draw);
511: PetscDrawSetFromOptions(user.draw);
512: PetscDrawSPCreate(user.draw, 1, &user.drawks);
513: TSMonitorSet(ts, KSConv, &user, NULL);
514: }
515: TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user);
516: TSSetFromOptions(ts);
517: TSSetComputeInitialCondition(ts, InitializeSolve);
518: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &w);
519: VecDuplicate(w, &u);
520: VecCopy(w, u);
521: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &w);
522: TSComputeInitialCondition(ts, u);
523: TSSolve(ts, u);
524: if (user.monitorhg) {
525: PetscDrawSave(user.draw);
526: PetscDrawHGDestroy(&user.drawhg);
527: PetscDrawDestroy(&user.draw);
528: }
529: if (user.monitorsp) {
530: PetscDrawSave(user.draw);
531: PetscDrawSPDestroy(&user.drawsp);
532: PetscDrawDestroy(&user.draw);
533: }
534: if (user.monitorks) {
535: PetscDrawSave(user.draw);
536: PetscDrawSPDestroy(&user.drawks);
537: PetscDrawDestroy(&user.draw);
538: }
539: VecDestroy(&u);
540: TSDestroy(&ts);
541: DMDestroy(&sw);
542: DMDestroy(&dm);
543: PetscFinalize();
544: return 0;
545: }
547: /*TEST
548: build:
549: requires: double !complex
550: test:
551: suffix: 1
552: 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
553: test:
554: suffix: 2
555: 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
556: TEST*/