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