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: {

 44:   DMCreate(comm, dm);
 45:   DMSetType(*dm, DMPLEX);
 46:   DMSetFromOptions(*dm);
 47:   DMViewFromOptions(*dm, NULL, "-dm_view");
 48:   DMGetDimension(*dm, &user->dim);
 49:   return(0);
 50: }

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

 63:   PetscRandomCreate(PetscObjectComm((PetscObject) dmSw), &rnd);
 64:   PetscRandomSetInterval(rnd, -1.0, 1.0);
 65:   PetscRandomSetFromOptions(rnd);

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

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

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

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

120:       initialConditions[(n*2 + 0)*dim + 0] = n+1;
121:       initialConditions[(n*2 + 0)*dim + 1] = 0;
122:       initialConditions[(n*2 + 1)*dim + 0] = 0;
123:       initialConditions[(n*2 + 1)*dim + 1] = PetscSqrtReal(1000./(n+1.));
124:     }
125:   }
126:   VecRestoreArray(u, &initialConditions);
127:   return(0);
128: }

130: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
131: {
132:   PetscInt      *cellid;
133:   PetscInt       dim, cStart, cEnd, c, Np = user->particlesPerCell, p;

137:   DMGetDimension(dm, &dim);
138:   DMCreate(PetscObjectComm((PetscObject) dm), sw);
139:   DMSetType(*sw, DMSWARM);
140:   DMSetDimension(*sw, dim);

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

154:       cellid[n] = c;
155:     }
156:   }
157:   DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
158:   PetscObjectSetName((PetscObject) *sw, "Particles");
159:   DMViewFromOptions(*sw, NULL, "-sw_view");
160:   return(0);
161: }

163: /* Create particle RHS Functions for gravity with G = 1 for simplification */
164: static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
165: {
166:   const PetscScalar *v;
167:   PetscScalar       *xres;
168:   PetscInt          Np, p, dim, d;
169:   PetscErrorCode    ierr;

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

188: static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
189: {
190:   const PetscScalar *x;
191:   PetscScalar       *vres;
192:   PetscInt          Np, p, dim, d;
193:   PetscErrorCode    ierr;

196:   /* The DM is not currently pushed down to the splits */
197:   dim  = ((AppCtx *) ctx)->dim;
198:   VecGetLocalSize(Vres, &Np);
199:   Np  /= dim;
200:   VecGetArray(Vres, &vres);
201:   VecGetArrayRead(X, &x);
202:   for (p = 0; p < Np; ++p) {
203:     const PetscScalar rsqr = DMPlex_NormD_Internal(dim, &x[p*dim]);

205:     for (d = 0; d < dim; ++d) {
206:       vres[p*dim+d] = -(1000./(p+1.)) * x[p*dim+d]/PetscSqr(rsqr);
207:     }
208:   }
209:   VecRestoreArray(Vres, &vres);
210:   VecRestoreArrayRead(X, &x);
211:   return(0);
212: }

214: static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t , Vec U, Vec R, void *ctx)
215: {
216:   DM                dm;
217:   const PetscScalar *u;
218:   PetscScalar       *r;
219:   PetscInt          Np, p, dim, d;
220:   PetscErrorCode    ierr;

223:   TSGetDM(ts, &dm);
224:   DMGetDimension(dm, &dim);
225:   VecGetLocalSize(U, &Np);
226:   Np  /= 2*dim;
227:   VecGetArrayRead(U, &u);
228:   VecGetArray(R, &r);
229:   for (p = 0; p < Np; ++p) {
230:     const PetscScalar rsqr = DMPlex_NormD_Internal(dim, &u[p*2*dim]);

232:     for (d = 0; d < dim; ++d) {
233:         r[p*2*dim+d]   = u[p*2*dim+d+2];
234:         r[p*2*dim+d+2] = -(1000./(1.+p)) * u[p*2*dim+d]/PetscSqr(rsqr);
235:     }
236:   }
237:   VecRestoreArrayRead(U, &u);
238:   VecRestoreArray(R, &r);
239:   return(0);
240: }

242: static PetscErrorCode InitializeSolve(TS ts, Vec u)
243: {
244:   DM             dm;
245:   AppCtx        *user;

249:   TSGetDM(ts, &dm);
250:   DMGetApplicationContext(dm, &user);
251:   SetInitialCoordinates(dm);
252:   SetInitialConditions(dm, u);
253:   return(0);
254: }

256: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
257: {
258:   MPI_Comm           comm;
259:   DM                 sdm;
260:   AppCtx            *user;
261:   const PetscScalar *u, *coords;
262:   PetscScalar       *e;
263:   PetscReal          t;
264:   PetscInt           dim, Np, p;
265:   PetscErrorCode     ierr;

268:   PetscObjectGetComm((PetscObject) ts, &comm);
269:   TSGetDM(ts, &sdm);
270:   DMGetApplicationContext(sdm, &user);
271:   DMGetDimension(sdm, &dim);
272:   TSGetSolveTime(ts, &t);
273:   VecGetArray(E, &e);
274:   VecGetArrayRead(U, &u);
275:   VecGetLocalSize(U, &Np);
276:   Np  /= 2*dim;
277:   DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
278:   for (p = 0; p < Np; ++p) {
279:     const PetscScalar *x     = &u[(p*2+0)*dim];
280:     const PetscScalar *v     = &u[(p*2+1)*dim];
281:     const PetscReal    x0    = p+1.;
282:     const PetscReal    omega = PetscSqrtReal(1000./(p+1.))/x0;
283:     const PetscReal    xe[3] = { x0*PetscCosReal(omega*t),       x0*PetscSinReal(omega*t),       0.0};
284:     const PetscReal    ve[3] = {-x0*omega*PetscSinReal(omega*t), x0*omega*PetscCosReal(omega*t), 0.0};
285:     PetscInt           d;

287:     for (d = 0; d < dim; ++d) {
288:       e[(p*2+0)*dim+d] = x[d] - xe[d];
289:       e[(p*2+1)*dim+d] = v[d] - ve[d];
290:     }
291:     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))));}
292:   }
293:   DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
294:   VecRestoreArrayRead(U, &u);
295:   VecRestoreArray(E, &e);
296:   return(0);
297: }

299: int main(int argc, char **argv)
300: {
301:   TS                 ts;
302:   TSConvergedReason  reason;
303:   DM                 dm, sw;
304:   Vec                u;
305:   IS                 is1, is2;
306:   PetscInt          *idx1, *idx2;
307:   MPI_Comm           comm;
308:   AppCtx             user;
309:   const PetscScalar *endVals;
310:   PetscReal          ftime   = .1;
311:   PetscInt           locSize, dim, d, Np, p, steps;
312:   PetscErrorCode     ierr;

314:   PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
315:   comm = PETSC_COMM_WORLD;
316:   ProcessOptions(comm, &user);

318:   CreateMesh(comm, &dm, &user);
319:   CreateParticles(dm, &sw, &user);
320:   DMSetApplicationContext(sw, &user);

322:   TSCreate(comm, &ts);
323:   TSSetType(ts, TSBASICSYMPLECTIC);
324:   TSSetDM(ts, sw);
325:   TSSetMaxTime(ts, ftime);
326:   TSSetTimeStep(ts, 0.0001);
327:   TSSetMaxSteps(ts, 10);
328:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
329:   TSSetTime(ts, 0.0);
330:   TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user);

332:   DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u);
333:   DMGetDimension(sw, &dim);
334:   VecGetLocalSize(u, &locSize);
335:   Np   = locSize/(2*dim);
336:   PetscMalloc1(locSize/2, &idx1);
337:   PetscMalloc1(locSize/2, &idx2);
338:   for (p = 0; p < Np; ++p) {
339:     for (d = 0; d < dim; ++d) {
340:       idx1[p*dim+d] = (p*2+0)*dim + d;
341:       idx2[p*dim+d] = (p*2+1)*dim + d;
342:     }
343:   }
344:   ISCreateGeneral(comm, locSize/2, idx1, PETSC_OWN_POINTER, &is1);
345:   ISCreateGeneral(comm, locSize/2, idx2, PETSC_OWN_POINTER, &is2);
346:   TSRHSSplitSetIS(ts, "position", is1);
347:   TSRHSSplitSetIS(ts, "momentum", is2);
348:   ISDestroy(&is1);
349:   ISDestroy(&is2);
350:   TSRHSSplitSetRHSFunction(ts,"position",NULL,RHSFunction1,&user);
351:   TSRHSSplitSetRHSFunction(ts,"momentum",NULL,RHSFunction2,&user);

353:   TSSetFromOptions(ts);
354:   TSSetComputeInitialCondition(ts, InitializeSolve);
355:   TSSetComputeExactError(ts, ComputeError);
356:   TSComputeInitialCondition(ts, u);
357:   VecViewFromOptions(u, NULL, "-init_view");
358:   TSSolve(ts, u);
359:   TSGetSolveTime(ts, &ftime);
360:   TSGetConvergedReason(ts, &reason);
361:   TSGetStepNumber(ts, &steps);
362:   PetscPrintf(comm,"%s at time %g after %D steps\n", TSConvergedReasons[reason], (double) ftime, steps);

364:   VecGetArrayRead(u, &endVals);
365:   for (p = 0; p < Np; ++p) {
366:     const PetscReal norm = DMPlex_NormD_Internal(dim, &endVals[(p*2 + 1)*dim]);
367:     PetscPrintf(comm, "Particle %D initial Energy: %g  Final Energy: %g\n", p, (double) (0.5*(1000./(p+1))), (double) 0.5*PetscSqr(norm));
368:   }
369:   VecRestoreArrayRead(u, &endVals);
370:   DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u);
371:   TSDestroy(&ts);
372:   DMDestroy(&sw);
373:   DMDestroy(&dm);
374:   PetscFinalize();
375:   return ierr;
376: }

378: /*TEST

380:    build:
381:      requires: triangle !single !complex
382:    test:
383:      suffix: bsi1
384:      args: -dm_plex_box_faces 1,1 -dm_view -sw_view -ts_basicsymplectic_type 1 -ts_max_time 0.1 -ts_monitor_sp_swarm -ts_convergence_estimate -convest_num_refine 2
385:    test:
386:      suffix: bsi2
387:      args: -dm_plex_box_faces 1,1 -dm_view -sw_view -ts_basicsymplectic_type 2 -ts_max_time 0.1 -ts_monitor_sp_swarm -ts_convergence_estimate -convest_num_refine 2
388:    test:
389:      suffix: euler
390:      args: -dm_plex_box_faces 1,1 -dm_view -sw_view -ts_type euler -ts_max_time 0.1 -ts_monitor_sp_swarm -ts_convergence_estimate -convest_num_refine 2

392: TEST*/