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();

 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

 46: /* Create the mesh for velocity space */
 47: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
 48: {
 49:   PetscFunctionBeginUser;
 50:   PetscCall(DMCreate(comm, dm));
 51:   PetscCall(DMSetType(*dm, DMPLEX));
 52:   PetscCall(DMSetFromOptions(*dm));
 53:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 54:   PetscFunctionReturn(PETSC_SUCCESS);
 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;

 68:   PetscFunctionBeginUser;
 69:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)sw), &rnd));
 70:   PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
 71:   PetscCall(PetscRandomSetFromOptions(rnd));

 73:   PetscCall(DMGetApplicationContext(sw, &user));
 74:   Np = user->particlesPerCell;
 75:   PetscCall(DMGetDimension(sw, &dim));
 76:   PetscCall(DMSwarmGetCellDM(sw, &dm));
 77:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
 78:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 79:   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
 80:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
 81:   PetscCall(PetscMalloc5(dim, &centroid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
 82:   for (d = 0; d < dim; ++d) xi0[d] = -1.0;
 83:   PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
 84:   PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&vals));
 85:   for (c = cStart; c < cEnd; ++c) {
 86:     if (Np == 1) {
 87:       PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL));
 88:       for (d = 0; d < dim; ++d) {
 89:         coords[c * dim + d] = centroid[d];
 90:         if ((coords[c * dim + d] >= -1) && (coords[c * dim + d] <= 1)) {
 91:           vals[c] = 1.0;
 92:         } else {
 93:           vals[c] = 0.;
 94:         }
 95:       }
 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:   PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
114:   PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&vals));
115:   PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
116:   PetscCall(PetscRandomDestroy(&rnd));
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

120: /* The initial conditions are just the initial particle weights */
121: static PetscErrorCode SetInitialConditions(DM dmSw, Vec u)
122: {
123:   DM           dm;
124:   AppCtx      *user;
125:   PetscReal   *vals;
126:   PetscScalar *initialConditions;
127:   PetscInt     dim, d, cStart, cEnd, c, Np, p, n;

129:   PetscFunctionBeginUser;
130:   PetscCall(VecGetLocalSize(u, &n));
131:   PetscCall(DMGetApplicationContext(dmSw, &user));
132:   Np = user->particlesPerCell;
133:   PetscCall(DMSwarmGetCellDM(dmSw, &dm));
134:   PetscCall(DMGetDimension(dm, &dim));
135:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
136:   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);
137:   PetscCall(DMSwarmGetField(dmSw, "w_q", NULL, NULL, (void **)&vals));
138:   PetscCall(VecGetArray(u, &initialConditions));
139:   for (c = cStart; c < cEnd; ++c) {
140:     for (p = 0; p < Np; ++p) {
141:       const PetscInt n = c * Np + p;
142:       for (d = 0; d < dim; d++) initialConditions[n] = vals[n];
143:     }
144:   }
145:   PetscCall(VecRestoreArray(u, &initialConditions));
146:   PetscCall(DMSwarmRestoreField(dmSw, "w_q", NULL, NULL, (void **)&vals));
147:   PetscFunctionReturn(PETSC_SUCCESS);
148: }

150: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
151: {
152:   PetscInt *cellid;
153:   PetscInt  dim, cStart, cEnd, c, Np = user->particlesPerCell, p;

155:   PetscFunctionBeginUser;
156:   PetscCall(DMGetDimension(dm, &dim));
157:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
158:   PetscCall(DMSetType(*sw, DMSWARM));
159:   PetscCall(DMSetDimension(*sw, dim));
160:   PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
161:   PetscCall(DMSwarmSetCellDM(*sw, dm));
162:   PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", dim, PETSC_REAL));
163:   PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
164:   PetscCall(DMSwarmFinalizeFieldRegister(*sw));
165:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
166:   PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0));
167:   PetscCall(DMSetFromOptions(*sw));
168:   PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
169:   for (c = cStart; c < cEnd; ++c) {
170:     for (p = 0; p < Np; ++p) {
171:       const PetscInt n = c * Np + p;
172:       cellid[n]        = c;
173:     }
174:   }
175:   PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
176:   PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
177:   PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
178:   PetscFunctionReturn(PETSC_SUCCESS);
179: }

181: /*
182:   f_t = 1/\tau \left( f_eq - f \right)
183:   n_t = 1/\tau \left( \int f_eq - \int f \right) = 1/\tau (n - n) = 0
184:   v_t = 1/\tau \left( \int v f_eq - \int v f \right) = 1/\tau (v - v) = 0
185:   E_t = 1/\tau \left( \int v^2 f_eq - \int v^2 f \right) = 1/\tau (T - T) = 0

187:   Let's look at a single cell:

189:     \int_C f_t             = 1/\tau \left( \int_C f_eq - \int_C f \right)
190:     \sum_{x_i \in C} w^i_t = 1/\tau \left( F_eq(C) - \sum_{x_i \in C} w_i \right)
191: */

193: /* This computes the 1D Maxwellian distribution for given mass n, velocity v, and temperature T */
194: static PetscReal ComputePDF(PetscReal m, PetscReal n, PetscReal T, PetscReal v[])
195: {
196:   return (n / PetscSqrtReal(2.0 * PETSC_PI * T / m)) * PetscExpReal(-0.5 * m * PetscSqr(v[0]) / T);
197: }

199: /*
200:   erf z = \frac{2}{\sqrt\pi} \int^z_0 dt e^{-t^2} and erf \infty = 1

202:   We begin with our distribution

204:     \sqrt{\frac{m}{2 \pi T}} e^{-m v^2/2T}

206:   Let t = \sqrt{m/2T} v, z = \sqrt{m/2T} w, so that we now have

208:       \sqrt{\frac{m}{2 \pi T}} \int^w_0 dv e^{-m v^2/2T}
209:     = \sqrt{\frac{m}{2 \pi T}} \int^{\sqrt{m/2T} w}_0 \sqrt{2T/m} dt e^{-t^2}
210:     = 1/\sqrt{\pi} \int^{\sqrt{m/2T} w}_0 dt e^{-t^2}
211:     = 1/2 erf(\sqrt{m/2T} w)
212: */
213: static PetscReal ComputeCDF(PetscReal m, PetscReal n, PetscReal T, PetscReal va, PetscReal vb)
214: {
215:   PetscReal alpha = PetscSqrtReal(0.5 * m / T);
216:   PetscReal za    = alpha * va;
217:   PetscReal zb    = alpha * vb;
218:   PetscReal sum   = 0.0;

220:   sum += zb >= 0 ? erf(zb) : -erf(-zb);
221:   sum -= za >= 0 ? erf(za) : -erf(-za);
222:   return 0.5 * n * sum;
223: }

225: static PetscErrorCode CheckDistribution(DM dm, PetscReal m, PetscReal n, PetscReal T, PetscReal v[])
226: {
227:   PetscSection coordSection;
228:   Vec          coordsLocal;
229:   PetscReal   *xq, *wq;
230:   PetscReal    vmin, vmax, neq, veq, Teq;
231:   PetscInt     Nq = 100, q, cStart, cEnd, c;

233:   PetscFunctionBeginUser;
234:   PetscCall(DMGetBoundingBox(dm, &vmin, &vmax));
235:   /* Check analytic over entire line */
236:   neq = ComputeCDF(m, n, T, vmin, vmax);
237:   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));
238:   /* Check analytic over cells */
239:   PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
240:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
241:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
242:   neq = 0.0;
243:   for (c = cStart; c < cEnd; ++c) {
244:     PetscScalar *vcoords = NULL;

246:     PetscCall(DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &vcoords));
247:     neq += ComputeCDF(m, n, T, vcoords[0], vcoords[1]);
248:     PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &vcoords));
249:   }
250:   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));
251:   /* Check quadrature over entire line */
252:   PetscCall(PetscMalloc2(Nq, &xq, Nq, &wq));
253:   PetscCall(PetscDTGaussQuadrature(100, vmin, vmax, xq, wq));
254:   neq = 0.0;
255:   for (q = 0; q < Nq; ++q) neq += ComputePDF(m, n, T, &xq[q]) * wq[q];
256:   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));
257:   /* Check omemnts with quadrature */
258:   veq = 0.0;
259:   for (q = 0; q < Nq; ++q) veq += xq[q] * ComputePDF(m, n, T, &xq[q]) * wq[q];
260:   veq /= neq;
261:   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]));
262:   Teq = 0.0;
263:   for (q = 0; q < Nq; ++q) Teq += PetscSqr(xq[q]) * ComputePDF(m, n, T, &xq[q]) * wq[q];
264:   Teq = Teq * m / neq - PetscSqr(veq);
265:   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));
266:   PetscCall(PetscFree2(xq, wq));
267:   PetscFunctionReturn(PETSC_SUCCESS);
268: }

270: static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t, Vec U, Vec R, void *ctx)
271: {
272:   const PetscScalar *u;
273:   PetscSection       coordSection;
274:   Vec                coordsLocal;
275:   PetscScalar       *r, *coords;
276:   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;
277:   PetscInt           dim, d, Np, Ncp, p, cStart, cEnd, c;
278:   DM                 dmSw, plex;

280:   PetscFunctionBeginUser;
281:   PetscCall(VecGetLocalSize(U, &Np));
282:   PetscCall(VecGetArrayRead(U, &u));
283:   PetscCall(VecGetArray(R, &r));
284:   PetscCall(TSGetDM(ts, &dmSw));
285:   PetscCall(DMSwarmGetCellDM(dmSw, &plex));
286:   PetscCall(DMGetDimension(dmSw, &dim));
287:   PetscCall(DMGetCoordinatesLocal(plex, &coordsLocal));
288:   PetscCall(DMGetCoordinateSection(plex, &coordSection));
289:   PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
290:   Np /= dim;
291:   Ncp = Np / (cEnd - cStart);
292:   /* Calculate moments of particle distribution, note that velocity is in the coordinate */
293:   PetscCall(DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
294:   for (p = 0; p < Np; ++p) {
295:     PetscReal m1 = 0.0, m2 = 0.0;

297:     for (d = 0; d < dim; ++d) {
298:       m1 += PetscRealPart(coords[p * dim + d]);
299:       m2 += PetscSqr(PetscRealPart(coords[p * dim + d]));
300:     }
301:     n += u[p];
302:     v += u[p] * m1;
303:     E += u[p] * m2;
304:   }
305:   v /= n;
306:   T = E * m / n - v * v;
307:   PetscCall(PetscInfo(ts, "Time %.2f: mass %.4f velocity: %+.4f temperature: %.4f\n", (double)t, (double)n, (double)v, (double)T));
308:   PetscCall(CheckDistribution(plex, m, n, T, &v));
309:   /*
310:      Begin cellwise evaluation of the collision operator. Essentially, penalize the weights of the particles
311:      in that cell against the maxwellian for the number of particles expected to be in that cell
312:   */
313:   for (c = cStart; c < cEnd; ++c) {
314:     PetscScalar *vcoords    = NULL;
315:     PetscReal    relaxation = 1.0, neq;
316:     PetscInt     sp         = c * Ncp, q;

318:     /* Calculate equilibrium occupation for this velocity cell */
319:     PetscCall(DMPlexVecGetClosure(plex, coordSection, coordsLocal, c, NULL, &vcoords));
320:     neq = ComputeCDF(m, n, T, vcoords[0], vcoords[1]);
321:     PetscCall(DMPlexVecRestoreClosure(plex, coordSection, coordsLocal, c, NULL, &vcoords));
322:     for (q = 0; q < Ncp; ++q) r[sp + q] = (1.0 / relaxation) * (neq - u[sp + q]);
323:   }
324:   /* Check update */
325:   for (p = 0; p < Np; ++p) {
326:     PetscReal    m1 = 0.0, m2 = 0.0;
327:     PetscScalar *vcoords = NULL;

329:     for (d = 0; d < dim; ++d) {
330:       m1 += PetscRealPart(coords[p * dim + d]);
331:       m2 += PetscSqr(PetscRealPart(coords[p * dim + d]));
332:     }
333:     cn += r[p];
334:     cv += r[p] * m1;
335:     cE += r[p] * m2;
336:     pE += u[p] * m2;
337:     PetscCall(DMPlexVecGetClosure(plex, coordSection, coordsLocal, p, NULL, &vcoords));
338:     eqE += ComputeCDF(m, n, T, vcoords[0], vcoords[1]) * m2;
339:     PetscCall(DMPlexVecRestoreClosure(plex, coordSection, coordsLocal, p, NULL, &vcoords));
340:   }
341:   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));
342:   PetscCall(DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
343:   PetscCall(VecRestoreArrayRead(U, &u));
344:   PetscCall(VecRestoreArray(R, &r));
345:   PetscFunctionReturn(PETSC_SUCCESS);
346: }

348: static PetscErrorCode HGMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
349: {
350:   AppCtx            *user = (AppCtx *)ctx;
351:   const PetscScalar *u;
352:   DM                 sw, dm;
353:   PetscInt           dim, Np, p;

355:   PetscFunctionBeginUser;
356:   if (step < 0) PetscFunctionReturn(PETSC_SUCCESS);
357:   if (((user->ostep > 0) && (!(step % user->ostep)))) {
358:     PetscDrawAxis axis;

360:     PetscCall(TSGetDM(ts, &sw));
361:     PetscCall(DMSwarmGetCellDM(sw, &dm));
362:     PetscCall(DMGetDimension(dm, &dim));
363:     PetscCall(PetscDrawHGReset(user->drawhg));
364:     PetscCall(PetscDrawHGGetAxis(user->drawhg, &axis));
365:     PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "V", "f(V)"));
366:     PetscCall(PetscDrawAxisSetLimits(axis, -3, 3, 0, 100));
367:     PetscCall(PetscDrawHGSetLimits(user->drawhg, -3.0, 3.0, 0, 10));
368:     PetscCall(VecGetLocalSize(U, &Np));
369:     Np /= dim;
370:     PetscCall(VecGetArrayRead(U, &u));
371:     /* get points from solution vector */
372:     for (p = 0; p < Np; ++p) PetscCall(PetscDrawHGAddValue(user->drawhg, u[p]));
373:     PetscCall(PetscDrawHGDraw(user->drawhg));
374:     PetscCall(VecRestoreArrayRead(U, &u));
375:   }
376:   PetscFunctionReturn(PETSC_SUCCESS);
377: }

379: static PetscErrorCode SPMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
380: {
381:   AppCtx            *user = (AppCtx *)ctx;
382:   const PetscScalar *u;
383:   PetscReal         *v, *coords;
384:   PetscInt           Np, p;
385:   DM                 dmSw;

387:   PetscFunctionBeginUser;

389:   if (step < 0) PetscFunctionReturn(PETSC_SUCCESS);
390:   if (((user->ostep > 0) && (!(step % user->ostep)))) {
391:     PetscDrawAxis axis;

393:     PetscCall(TSGetDM(ts, &dmSw));
394:     PetscCall(PetscDrawSPReset(user->drawsp));
395:     PetscCall(PetscDrawSPGetAxis(user->drawsp, &axis));
396:     PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "V", "w"));
397:     PetscCall(VecGetLocalSize(U, &Np));
398:     PetscCall(VecGetArrayRead(U, &u));
399:     /* get points from solution vector */
400:     PetscCall(PetscMalloc1(Np, &v));
401:     for (p = 0; p < Np; ++p) v[p] = PetscRealPart(u[p]);
402:     PetscCall(DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
403:     for (p = 0; p < Np - 1; ++p) PetscCall(PetscDrawSPAddPoint(user->drawsp, &coords[p], &v[p]));
404:     PetscCall(PetscDrawSPDraw(user->drawsp, PETSC_TRUE));
405:     PetscCall(VecRestoreArrayRead(U, &u));
406:     PetscCall(DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
407:     PetscCall(PetscFree(v));
408:   }
409:   PetscFunctionReturn(PETSC_SUCCESS);
410: }

412: static PetscErrorCode KSConv(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
413: {
414:   AppCtx            *user = (AppCtx *)ctx;
415:   const PetscScalar *u;
416:   PetscScalar        sup;
417:   PetscReal         *v, *coords, T = 0., vel = 0., step_cast, w_sum;
418:   PetscInt           dim, Np, p, cStart, cEnd;
419:   DM                 sw, plex;

421:   PetscFunctionBeginUser;
422:   if (step < 0) PetscFunctionReturn(PETSC_SUCCESS);
423:   if (((user->ostep > 0) && (!(step % user->ostep)))) {
424:     PetscDrawAxis axis;
425:     PetscCall(PetscDrawSPGetAxis(user->drawks, &axis));
426:     PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "t", "D_n"));
427:     PetscCall(PetscDrawSPSetLimits(user->drawks, 0., 100, 0., 3.5));
428:     PetscCall(TSGetDM(ts, &sw));
429:     PetscCall(VecGetLocalSize(U, &Np));
430:     PetscCall(VecGetArrayRead(U, &u));
431:     /* get points from solution vector */
432:     PetscCall(PetscMalloc1(Np, &v));
433:     PetscCall(DMSwarmGetCellDM(sw, &plex));
434:     PetscCall(DMGetDimension(sw, &dim));
435:     PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
436:     for (p = 0; p < Np; ++p) v[p] = PetscRealPart(u[p]);
437:     PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
438:     w_sum = 0.;
439:     for (p = 0; p < Np; ++p) {
440:       w_sum += u[p];
441:       T += u[p] * coords[p] * coords[p];
442:       vel += u[p] * coords[p];
443:     }
444:     vel /= w_sum;
445:     T   = T / w_sum - vel * vel;
446:     sup = 0.0;
447:     for (p = 0; p < Np; ++p) {
448:       PetscReal tmp = 0.;
449:       tmp           = PetscAbs(u[p] - ComputePDF(1.0, w_sum, T, &coords[p * dim]));
450:       if (tmp > sup) sup = tmp;
451:     }
452:     step_cast = (PetscReal)step;
453:     PetscCall(PetscDrawSPAddPoint(user->drawks, &step_cast, &sup));
454:     PetscCall(PetscDrawSPDraw(user->drawks, PETSC_TRUE));
455:     PetscCall(VecRestoreArrayRead(U, &u));
456:     PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
457:     PetscCall(PetscFree(v));
458:   }
459:   PetscFunctionReturn(PETSC_SUCCESS);
460: }

462: static PetscErrorCode InitializeSolve(TS ts, Vec u)
463: {
464:   DM      dm;
465:   AppCtx *user;

467:   PetscFunctionBeginUser;
468:   PetscCall(TSGetDM(ts, &dm));
469:   PetscCall(DMGetApplicationContext(dm, &user));
470:   PetscCall(SetInitialCoordinates(dm));
471:   PetscCall(SetInitialConditions(dm, u));
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }
474: /*
475:      A single particle is generated for each velocity space cell using the dmswarmpicfield_coor and is used to evaluate collisions in that cell.
476:      0 weight ghost particles are initialized outside of a small velocity domain to ensure the tails of the amxwellian are resolved.
477:    */
478: int main(int argc, char **argv)
479: {
480:   TS       ts;     /* nonlinear solver */
481:   DM       dm, sw; /* Velocity space mesh and Particle Swarm */
482:   Vec      u, w;   /* swarm vector */
483:   MPI_Comm comm;
484:   AppCtx   user;

486:   PetscFunctionBeginUser;
487:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
488:   comm = PETSC_COMM_WORLD;
489:   PetscCall(ProcessOptions(comm, &user));
490:   PetscCall(CreateMesh(comm, &dm, &user));
491:   PetscCall(CreateParticles(dm, &sw, &user));
492:   PetscCall(DMSetApplicationContext(sw, &user));
493:   PetscCall(TSCreate(comm, &ts));
494:   PetscCall(TSSetDM(ts, sw));
495:   PetscCall(TSSetMaxTime(ts, 10.0));
496:   PetscCall(TSSetTimeStep(ts, 0.01));
497:   PetscCall(TSSetMaxSteps(ts, 100000));
498:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
499:   if (user.monitorhg) {
500:     PetscCall(PetscDrawCreate(comm, NULL, "monitor", 0, 0, 400, 300, &user.draw));
501:     PetscCall(PetscDrawSetFromOptions(user.draw));
502:     PetscCall(PetscDrawHGCreate(user.draw, 20, &user.drawhg));
503:     PetscCall(PetscDrawHGSetColor(user.drawhg, 3));
504:     PetscCall(TSMonitorSet(ts, HGMonitor, &user, NULL));
505:   } else if (user.monitorsp) {
506:     PetscCall(PetscDrawCreate(comm, NULL, "monitor", 0, 0, 400, 300, &user.draw));
507:     PetscCall(PetscDrawSetFromOptions(user.draw));
508:     PetscCall(PetscDrawSPCreate(user.draw, 1, &user.drawsp));
509:     PetscCall(TSMonitorSet(ts, SPMonitor, &user, NULL));
510:   } else if (user.monitorks) {
511:     PetscCall(PetscDrawCreate(comm, NULL, "monitor", 0, 0, 400, 300, &user.draw));
512:     PetscCall(PetscDrawSetFromOptions(user.draw));
513:     PetscCall(PetscDrawSPCreate(user.draw, 1, &user.drawks));
514:     PetscCall(TSMonitorSet(ts, KSConv, &user, NULL));
515:   }
516:   PetscCall(TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user));
517:   PetscCall(TSSetFromOptions(ts));
518:   PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
519:   PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &w));
520:   PetscCall(VecDuplicate(w, &u));
521:   PetscCall(VecCopy(w, u));
522:   PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &w));
523:   PetscCall(TSComputeInitialCondition(ts, u));
524:   PetscCall(TSSolve(ts, u));
525:   if (user.monitorhg) {
526:     PetscCall(PetscDrawSave(user.draw));
527:     PetscCall(PetscDrawHGDestroy(&user.drawhg));
528:     PetscCall(PetscDrawDestroy(&user.draw));
529:   }
530:   if (user.monitorsp) {
531:     PetscCall(PetscDrawSave(user.draw));
532:     PetscCall(PetscDrawSPDestroy(&user.drawsp));
533:     PetscCall(PetscDrawDestroy(&user.draw));
534:   }
535:   if (user.monitorks) {
536:     PetscCall(PetscDrawSave(user.draw));
537:     PetscCall(PetscDrawSPDestroy(&user.drawks));
538:     PetscCall(PetscDrawDestroy(&user.draw));
539:   }
540:   PetscCall(VecDestroy(&u));
541:   PetscCall(TSDestroy(&ts));
542:   PetscCall(DMDestroy(&sw));
543:   PetscCall(DMDestroy(&dm));
544:   PetscCall(PetscFinalize());
545:   return 0;
546: }

548: /*TEST
549:    build:
550:      requires: double !complex
551:    test:
552:      suffix: 1
553:      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
554:    test:
555:      suffix: 2
556:      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
557: TEST*/