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, ¢roid, 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*/