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, &centroid, 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*/