Actual source code: ex5.c

  1: static char help[] = "Vlasov example of particles orbiting around a central massive point.\n";

  3: #include <petscdmplex.h>
  4: #include <petsc/private/dmpleximpl.h>
  5: #include <petsc/private/petscfeimpl.h>
  6: #include <petscdmswarm.h>
  7: #include <petscts.h>

  9: typedef struct {
 10:   PetscInt  dim;
 11:   PetscInt  particlesPerCell; /* The number of partices per cell */
 12:   PetscReal momentTol;        /* Tolerance for checking moment conservation */
 13:   PetscBool monitor;          /* Flag for using the TS monitor */
 14:   PetscBool error;            /* Flag for printing the error */
 15:   PetscInt  ostep;            /* print the energy at each ostep time steps */
 16: } AppCtx;

 18: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 19: {

 23:   options->monitor          = PETSC_FALSE;
 24:   options->error            = PETSC_FALSE;
 25:   options->particlesPerCell = 1;
 26:   options->momentTol        = 100.0*PETSC_MACHINE_EPSILON;
 27:   options->ostep            = 100;

 29:   PetscOptionsBegin(comm, "", "Vlasov Options", "DMPLEX");
 30:   PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL);
 31:   PetscOptionsBool("-monitor", "Flag to use the TS monitor", "ex4.c", options->monitor, &options->monitor, NULL);
 32:   PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL);
 33:   PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex4.c", options->particlesPerCell, &options->particlesPerCell, NULL);
 34:   PetscOptionsEnd();

 36:   return 0;
 37: }

 39: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
 40: {
 42:   DMCreate(comm, dm);
 43:   DMSetType(*dm, DMPLEX);
 44:   DMSetFromOptions(*dm);
 45:   DMViewFromOptions(*dm, NULL, "-dm_view");
 46:   DMGetDimension(*dm, &user->dim);
 47:   return 0;
 48: }

 50: static PetscErrorCode SetInitialCoordinates(DM dmSw)
 51: {
 52:   DM             dm;
 53:   AppCtx        *user;
 54:   PetscRandom    rnd;
 55:   PetscBool      simplex;
 56:   PetscReal     *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
 57:   PetscInt       dim, d, cStart, cEnd, c, Np, p;

 60:   PetscRandomCreate(PetscObjectComm((PetscObject) dmSw), &rnd);
 61:   PetscRandomSetInterval(rnd, -1.0, 1.0);
 62:   PetscRandomSetFromOptions(rnd);

 64:   DMGetApplicationContext(dmSw, &user);
 65:   Np   = user->particlesPerCell;
 66:   DMSwarmGetCellDM(dmSw, &dm);
 67:   DMGetDimension(dm, &dim);
 68:   DMPlexIsSimplex(dm, &simplex);
 69:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
 70:   PetscMalloc5(dim, &centroid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
 71:   for (d = 0; d < dim; ++d) xi0[d] = -1.0;
 72:   DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
 73:   for (c = cStart; c < cEnd; ++c) {
 74:     if (Np == 1) {
 75:       DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
 76:       for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
 77:     } else {
 78:       DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
 79:       for (p = 0; p < Np; ++p) {
 80:         const PetscInt n   = c*Np + p;
 81:         PetscReal      sum = 0.0, refcoords[3];

 83:         for (d = 0; d < dim; ++d) {
 84:           PetscRandomGetValueReal(rnd, &refcoords[d]);
 85:           sum += refcoords[d];
 86:         }
 87:         if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
 88:         CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
 89:       }
 90:     }
 91:   }
 92:   DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
 93:   PetscFree5(centroid, xi0, v0, J, invJ);
 94:   PetscRandomDestroy(&rnd);
 95:   return 0;
 96: }

 98: static PetscErrorCode SetInitialConditions(DM dmSw, Vec u)
 99: {
100:   DM             dm;
101:   AppCtx         *user;
102:   PetscScalar    *initialConditions;
103:   PetscInt       dim, cStart, cEnd, c, Np, p;

106:   DMGetApplicationContext(dmSw, &user);
107:   Np   = user->particlesPerCell;
108:   DMSwarmGetCellDM(dmSw, &dm);
109:   DMGetDimension(dm, &dim);
110:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
111:   VecGetArray(u, &initialConditions);
112:   for (c = cStart; c < cEnd; ++c) {
113:     for (p = 0; p < Np; ++p) {
114:       const PetscInt n = c*Np + p;

116:       initialConditions[(n*2 + 0)*dim + 0] = n+1;
117:       initialConditions[(n*2 + 0)*dim + 1] = 0;
118:       initialConditions[(n*2 + 1)*dim + 0] = 0;
119:       initialConditions[(n*2 + 1)*dim + 1] = PetscSqrtReal(1000./(n+1.));
120:     }
121:   }
122:   VecRestoreArray(u, &initialConditions);
123:   return 0;
124: }

126: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
127: {
128:   PetscInt      *cellid;
129:   PetscInt       dim, cStart, cEnd, c, Np = user->particlesPerCell, p;

132:   DMGetDimension(dm, &dim);
133:   DMCreate(PetscObjectComm((PetscObject) dm), sw);
134:   DMSetType(*sw, DMSWARM);
135:   DMSetDimension(*sw, dim);

137:   DMSwarmSetType(*sw, DMSWARM_PIC);
138:   DMSwarmSetCellDM(*sw, dm);
139:   DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2*dim, PETSC_REAL);
140:   DMSwarmFinalizeFieldRegister(*sw);
141:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
142:   DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0);
143:   DMSetFromOptions(*sw);
144:   DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
145:   for (c = cStart; c < cEnd; ++c) {
146:     for (p = 0; p < Np; ++p) {
147:       const PetscInt n = c*Np + p;

149:       cellid[n] = c;
150:     }
151:   }
152:   DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
153:   PetscObjectSetName((PetscObject) *sw, "Particles");
154:   DMViewFromOptions(*sw, NULL, "-sw_view");
155:   return 0;
156: }

158: /* Create particle RHS Functions for gravity with G = 1 for simplification */
159: static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
160: {
161:   const PetscScalar *v;
162:   PetscScalar       *xres;
163:   PetscInt          Np, p, dim, d;

166:   /* The DM is not currently pushed down to the splits */
167:   dim  = ((AppCtx *) ctx)->dim;
168:   VecGetLocalSize(Xres, &Np);
169:   Np  /= dim;
170:   VecGetArray(Xres, &xres);
171:   VecGetArrayRead(V, &v);
172:   for (p = 0; p < Np; ++p) {
173:      for (d = 0; d < dim; ++d) {
174:        xres[p*dim+d] = v[p*dim+d];
175:      }
176:   }
177:   VecRestoreArrayRead(V,& v);
178:   VecRestoreArray(Xres, &xres);
179:   return 0;
180: }

182: static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
183: {
184:   const PetscScalar *x;
185:   PetscScalar       *vres;
186:   PetscInt          Np, p, dim, d;

189:   /* The DM is not currently pushed down to the splits */
190:   dim  = ((AppCtx *) ctx)->dim;
191:   VecGetLocalSize(Vres, &Np);
192:   Np  /= dim;
193:   VecGetArray(Vres, &vres);
194:   VecGetArrayRead(X, &x);
195:   for (p = 0; p < Np; ++p) {
196:     const PetscScalar rsqr = DMPlex_NormD_Internal(dim, &x[p*dim]);

198:     for (d = 0; d < dim; ++d) {
199:       vres[p*dim+d] = -(1000./(p+1.)) * x[p*dim+d]/PetscSqr(rsqr);
200:     }
201:   }
202:   VecRestoreArray(Vres, &vres);
203:   VecRestoreArrayRead(X, &x);
204:   return 0;
205: }

207: static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t , Vec U, Vec R, void *ctx)
208: {
209:   DM                dm;
210:   const PetscScalar *u;
211:   PetscScalar       *r;
212:   PetscInt          Np, p, dim, d;

215:   TSGetDM(ts, &dm);
216:   DMGetDimension(dm, &dim);
217:   VecGetLocalSize(U, &Np);
218:   Np  /= 2*dim;
219:   VecGetArrayRead(U, &u);
220:   VecGetArray(R, &r);
221:   for (p = 0; p < Np; ++p) {
222:     const PetscScalar rsqr = DMPlex_NormD_Internal(dim, &u[p*2*dim]);

224:     for (d = 0; d < dim; ++d) {
225:         r[p*2*dim+d]   = u[p*2*dim+d+2];
226:         r[p*2*dim+d+2] = -(1000./(1.+p)) * u[p*2*dim+d]/PetscSqr(rsqr);
227:     }
228:   }
229:   VecRestoreArrayRead(U, &u);
230:   VecRestoreArray(R, &r);
231:   return 0;
232: }

234: static PetscErrorCode InitializeSolve(TS ts, Vec u)
235: {
236:   DM             dm;
237:   AppCtx        *user;

240:   TSGetDM(ts, &dm);
241:   DMGetApplicationContext(dm, &user);
242:   SetInitialCoordinates(dm);
243:   SetInitialConditions(dm, u);
244:   return 0;
245: }

247: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
248: {
249:   MPI_Comm          comm;
250:   DM                sdm;
251:   AppCtx            *user;
252:   const PetscScalar *u, *coords;
253:   PetscScalar       *e;
254:   PetscReal         t;
255:   PetscInt          dim, Np, p;

258:   PetscObjectGetComm((PetscObject) ts, &comm);
259:   TSGetDM(ts, &sdm);
260:   DMGetApplicationContext(sdm, &user);
261:   DMGetDimension(sdm, &dim);
262:   TSGetSolveTime(ts, &t);
263:   VecGetArray(E, &e);
264:   VecGetArrayRead(U, &u);
265:   VecGetLocalSize(U, &Np);
266:   Np  /= 2*dim;
267:   DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
268:   for (p = 0; p < Np; ++p) {
269:     const PetscScalar *x     = &u[(p*2+0)*dim];
270:     const PetscScalar *v     = &u[(p*2+1)*dim];
271:     const PetscReal   x0    = p+1.;
272:     const PetscReal   omega = PetscSqrtReal(1000./(p+1.))/x0;
273:     const PetscReal   xe[3] = { x0*PetscCosReal(omega*t),       x0*PetscSinReal(omega*t),       0.0};
274:     const PetscReal   ve[3] = {-x0*omega*PetscSinReal(omega*t), x0*omega*PetscCosReal(omega*t), 0.0};
275:     PetscInt          d;

277:     for (d = 0; d < dim; ++d) {
278:       e[(p*2+0)*dim+d] = x[d] - xe[d];
279:       e[(p*2+1)*dim+d] = v[d] - ve[d];
280:     }
281:     if (user->error) PetscPrintf(comm, "t %.4g: p%D error [%.2g %.2g] sol [(%.6lf %.6lf) (%.6lf %.6lf)] exact [(%.6lf %.6lf) (%.6lf %.6lf)] energy/exact energy %g / %g\n", t, p, (double) DMPlex_NormD_Internal(dim, &e[(p*2+0)*dim]), (double) DMPlex_NormD_Internal(dim, &e[(p*2+1)*dim]), (double) x[0], (double) x[1], (double) v[0], (double) v[1], (double) xe[0], (double) xe[1], (double) ve[0], (double) ve[1], 0.5*DMPlex_NormD_Internal(dim, v), (double) (0.5*(1000./(p+1))));
282:   }
283:   DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
284:   VecRestoreArrayRead(U, &u);
285:   VecRestoreArray(E, &e);
286:   return 0;
287: }

289: int main(int argc, char **argv)
290: {
291:   TS                 ts;
292:   TSConvergedReason  reason;
293:   DM                 dm, sw;
294:   Vec                u;
295:   IS                 is1, is2;
296:   PetscInt          *idx1, *idx2;
297:   MPI_Comm           comm;
298:   AppCtx             user;
299:   const PetscScalar *endVals;
300:   PetscReal          ftime   = .1;
301:   PetscInt           locSize, dim, d, Np, p, steps;

303:   PetscInitialize(&argc, &argv, NULL, help);
304:   comm = PETSC_COMM_WORLD;
305:   ProcessOptions(comm, &user);

307:   CreateMesh(comm, &dm, &user);
308:   CreateParticles(dm, &sw, &user);
309:   DMSetApplicationContext(sw, &user);

311:   TSCreate(comm, &ts);
312:   TSSetType(ts, TSBASICSYMPLECTIC);
313:   TSSetDM(ts, sw);
314:   TSSetMaxTime(ts, ftime);
315:   TSSetTimeStep(ts, 0.0001);
316:   TSSetMaxSteps(ts, 10);
317:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
318:   TSSetTime(ts, 0.0);
319:   TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user);

321:   DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u);
322:   DMGetDimension(sw, &dim);
323:   VecGetLocalSize(u, &locSize);
324:   Np   = locSize/(2*dim);
325:   PetscMalloc1(locSize/2, &idx1);
326:   PetscMalloc1(locSize/2, &idx2);
327:   for (p = 0; p < Np; ++p) {
328:     for (d = 0; d < dim; ++d) {
329:       idx1[p*dim+d] = (p*2+0)*dim + d;
330:       idx2[p*dim+d] = (p*2+1)*dim + d;
331:     }
332:   }
333:   ISCreateGeneral(comm, locSize/2, idx1, PETSC_OWN_POINTER, &is1);
334:   ISCreateGeneral(comm, locSize/2, idx2, PETSC_OWN_POINTER, &is2);
335:   TSRHSSplitSetIS(ts, "position", is1);
336:   TSRHSSplitSetIS(ts, "momentum", is2);
337:   ISDestroy(&is1);
338:   ISDestroy(&is2);
339:   TSRHSSplitSetRHSFunction(ts,"position",NULL,RHSFunction1,&user);
340:   TSRHSSplitSetRHSFunction(ts,"momentum",NULL,RHSFunction2,&user);

342:   TSSetFromOptions(ts);
343:   TSSetComputeInitialCondition(ts, InitializeSolve);
344:   TSSetComputeExactError(ts, ComputeError);
345:   TSComputeInitialCondition(ts, u);
346:   VecViewFromOptions(u, NULL, "-init_view");
347:   TSSolve(ts, u);
348:   TSGetSolveTime(ts, &ftime);
349:   TSGetConvergedReason(ts, &reason);
350:   TSGetStepNumber(ts, &steps);
351:   PetscPrintf(comm,"%s at time %g after %D steps\n", TSConvergedReasons[reason], (double) ftime, steps);

353:   VecGetArrayRead(u, &endVals);
354:   for (p = 0; p < Np; ++p) {
355:     const PetscReal norm = DMPlex_NormD_Internal(dim, &endVals[(p*2 + 1)*dim]);
356:     PetscPrintf(comm, "Particle %D initial Energy: %g  Final Energy: %g\n", p, (double) (0.5*(1000./(p+1))), (double) 0.5*PetscSqr(norm));
357:   }
358:   VecRestoreArrayRead(u, &endVals);
359:   DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u);
360:   TSDestroy(&ts);
361:   DMDestroy(&sw);
362:   DMDestroy(&dm);
363:   PetscFinalize();
364:   return 0;
365: }

367: /*TEST

369:    build:
370:      requires: triangle !single !complex
371:    test:
372:      suffix: bsi1
373:      args: -dm_plex_box_faces 1,1 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
374:            -ts_basicsymplectic_type 1 -ts_max_time 0.1 -ts_convergence_estimate -convest_num_refine 2 \
375:            -dm_view -sw_view -ts_monitor_sp_swarm -ts_monitor_sp_swarm_retain -ts_monitor_sp_swarm_phase 0
376:    test:
377:      suffix: bsi2
378:      args: -dm_plex_box_faces 1,1 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
379:            -ts_basicsymplectic_type 2 -ts_max_time 0.1 -ts_convergence_estimate -convest_num_refine 2 \
380:            -dm_view -sw_view -ts_monitor_sp_swarm -ts_monitor_sp_swarm_retain -ts_monitor_sp_swarm_phase 0
381:    test:
382:      suffix: euler
383:      args: -dm_plex_box_faces 1,1 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
384:            -ts_type euler -ts_max_time 0.1 -ts_convergence_estimate -convest_num_refine 2 \
385:            -dm_view -sw_view -ts_monitor_sp_swarm -ts_monitor_sp_swarm_retain -ts_monitor_sp_swarm_phase 0

387: TEST*/