Actual source code: ex5.c
petsc-3.14.6 2021-03-30
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; /* The topological mesh dimension */
11: PetscBool simplex; /* Flag for simplices or tensor cells */
12: char filename[PETSC_MAX_PATH_LEN]; /* Name of the mesh filename if any */
13: PetscInt particlesPerCell; /* The number of partices per cell */
14: PetscReal momentTol; /* Tolerance for checking moment conservation */
15: PetscBool monitor; /* Flag for using the TS monitor */
16: PetscBool error; /* Flag for printing the error */
17: PetscInt ostep; /* print the energy at each ostep time steps */
18: } AppCtx;
20: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
21: {
25: options->dim = 2;
26: options->simplex = PETSC_TRUE;
27: options->monitor = PETSC_FALSE;
28: options->error = PETSC_FALSE;
29: options->particlesPerCell = 1;
30: options->momentTol = 100.0*PETSC_MACHINE_EPSILON;
31: options->ostep = 100;
33: PetscStrcpy(options->filename, "");
35: PetscOptionsBegin(comm, "", "Vlasov Options", "DMPLEX");
36: PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL);
37: PetscOptionsInt("-dim", "The topological mesh dimension", "ex4.c", options->dim, &options->dim, NULL);
38: PetscOptionsBool("-monitor", "Flag to use the TS monitor", "ex4.c", options->monitor, &options->monitor, NULL);
39: PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL);
40: PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex4.c", options->simplex, &options->simplex, NULL);
41: PetscOptionsString("-mesh", "Name of the mesh filename if any", "ex4.c", options->filename, options->filename, sizeof(options->filename), NULL);
42: PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex4.c", options->particlesPerCell, &options->particlesPerCell, NULL);
43: PetscOptionsEnd();
45: return(0);
46: }
48: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
49: {
50: PetscBool flg;
54: PetscStrcmp(user->filename, "", &flg);
55: if (flg) {
56: DMPlexCreateBoxMesh(comm, user->dim, user->simplex, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);
57: } else {
58: DMPlexCreateFromFile(comm, user->filename, PETSC_TRUE, dm);
59: DMGetDimension(*dm, &user->dim);
60: }
61: {
62: DM distributedMesh = NULL;
64: DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
65: if (distributedMesh) {
66: DMDestroy(dm);
67: *dm = distributedMesh;
68: }
69: }
70: DMLocalizeCoordinates(*dm); /* needed for periodic */
71: DMSetFromOptions(*dm);
72: PetscObjectSetName((PetscObject) *dm, "Mesh");
73: DMViewFromOptions(*dm, NULL, "-dm_view");
74: return(0);
75: }
77: static PetscErrorCode SetInitialCoordinates(DM dmSw)
78: {
79: DM dm;
80: AppCtx *user;
81: PetscRandom rnd;
82: PetscBool simplex;
83: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
84: PetscInt dim, d, cStart, cEnd, c, Np, p;
88: PetscRandomCreate(PetscObjectComm((PetscObject) dmSw), &rnd);
89: PetscRandomSetInterval(rnd, -1.0, 1.0);
90: PetscRandomSetFromOptions(rnd);
92: DMGetApplicationContext(dmSw, (void **) &user);
93: simplex = user->simplex;
94: Np = user->particlesPerCell;
95: DMSwarmGetCellDM(dmSw, &dm);
96: DMGetDimension(dm, &dim);
97: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
98: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
99: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
100: DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
101: for (c = cStart; c < cEnd; ++c) {
102: if (Np == 1) {
103: DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
104: for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
105: } else {
106: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
107: for (p = 0; p < Np; ++p) {
108: const PetscInt n = c*Np + p;
109: PetscReal sum = 0.0, refcoords[3];
111: for (d = 0; d < dim; ++d) {
112: PetscRandomGetValueReal(rnd, &refcoords[d]);
113: sum += refcoords[d];
114: }
115: if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
116: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
117: }
118: }
119: }
120: DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
121: PetscFree5(centroid, xi0, v0, J, invJ);
122: PetscRandomDestroy(&rnd);
123: return(0);
124: }
126: static PetscErrorCode SetInitialConditions(DM dmSw, Vec u)
127: {
128: DM dm;
129: AppCtx *user;
130: PetscScalar *initialConditions;
131: PetscInt dim, cStart, cEnd, c, Np, p;
135: DMGetApplicationContext(dmSw, (void **) &user);
136: Np = user->particlesPerCell;
137: DMSwarmGetCellDM(dmSw, &dm);
138: DMGetDimension(dm, &dim);
139: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
140: VecGetArray(u, &initialConditions);
141: for (c = cStart; c < cEnd; ++c) {
142: for (p = 0; p < Np; ++p) {
143: const PetscInt n = c*Np + p;
145: initialConditions[(n*2 + 0)*dim + 0] = n+1;
146: initialConditions[(n*2 + 0)*dim + 1] = 0;
147: initialConditions[(n*2 + 1)*dim + 0] = 0;
148: initialConditions[(n*2 + 1)*dim + 1] = PetscSqrtReal(1000./(n+1.));
149: }
150: }
151: VecRestoreArray(u, &initialConditions);
152: return(0);
153: }
155: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
156: {
157: PetscInt *cellid;
158: PetscInt dim, cStart, cEnd, c, Np = user->particlesPerCell, p;
162: DMGetDimension(dm, &dim);
163: DMCreate(PetscObjectComm((PetscObject) dm), sw);
164: DMSetType(*sw, DMSWARM);
165: DMSetDimension(*sw, dim);
167: DMSwarmSetType(*sw, DMSWARM_PIC);
168: DMSwarmSetCellDM(*sw, dm);
169: DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2*dim, PETSC_REAL);
170: DMSwarmFinalizeFieldRegister(*sw);
171: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
172: DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0);
173: DMSetFromOptions(*sw);
174: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
175: for (c = cStart; c < cEnd; ++c) {
176: for (p = 0; p < Np; ++p) {
177: const PetscInt n = c*Np + p;
179: cellid[n] = c;
180: }
181: }
182: DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
183: PetscObjectSetName((PetscObject) *sw, "Particles");
184: DMViewFromOptions(*sw, NULL, "-sw_view");
185: return(0);
186: }
188: /* Create particle RHS Functions for gravity with G = 1 for simplification */
189: static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
190: {
191: DM dm;
192: const PetscScalar *v;
193: PetscScalar *xres;
194: PetscInt Np, p, dim, d;
195: PetscErrorCode ierr;
198: TSGetDM(ts, &dm);
199: DMGetDimension(dm, &dim);
200: VecGetLocalSize(Xres, &Np);
201: Np /= dim;
202: VecGetArray(Xres, &xres);
203: VecGetArrayRead(V, &v);
204: for (p = 0; p < Np; ++p) {
205: for (d = 0; d < dim; ++d) {
206: xres[p*dim+d] = v[p*dim+d];
207: }
208: }
209: VecRestoreArrayRead(V,& v);
210: VecRestoreArray(Xres, &xres);
211: return(0);
212: }
214: static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
215: {
216: DM dm;
217: const PetscScalar *x;
218: PetscScalar *vres;
219: PetscInt Np, p, dim, d;
220: PetscErrorCode ierr;
224: TSGetDM(ts, &dm);
225: DMGetDimension(dm, &dim);
226: VecGetLocalSize(Vres, &Np);
227: Np /= dim;
228: VecGetArray(Vres, &vres);
229: VecGetArrayRead(X, &x);
230: for (p = 0; p < Np; ++p) {
231: const PetscScalar rsqr = DMPlex_NormD_Internal(dim, &x[p*dim]);
233: for (d = 0; d < dim; ++d) {
234: vres[p*dim+d] = -(1000./(p+1.)) * x[p*dim+d]/PetscSqr(rsqr);
235: }
236: }
237: VecRestoreArray(Vres, &vres);
238: VecRestoreArrayRead(X, &x);
239: return(0);
240: }
242: static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t , Vec U, Vec R, void *ctx)
243: {
244: DM dm;
245: const PetscScalar *u;
246: PetscScalar *r;
247: PetscInt Np, p, dim, d;
248: PetscErrorCode ierr;
251: TSGetDM(ts, &dm);
252: DMGetDimension(dm, &dim);
253: VecGetLocalSize(U, &Np);
254: Np /= 2*dim;
255: VecGetArrayRead(U, &u);
256: VecGetArray(R, &r);
257: for (p = 0; p < Np; ++p) {
258: const PetscScalar rsqr = DMPlex_NormD_Internal(dim, &u[p*2*dim]);
260: for (d = 0; d < dim; ++d) {
261: r[p*2*dim+d] = u[p*2*dim+d+2];
262: r[p*2*dim+d+2] = -(1000./(1.+p)) * u[p*2*dim+d]/PetscSqr(rsqr);
263: }
264: }
265: VecRestoreArrayRead(U, &u);
266: VecRestoreArray(R, &r);
267: return(0);
268: }
270: static PetscErrorCode InitializeSolve(TS ts, Vec u)
271: {
272: DM dm;
273: AppCtx *user;
277: TSGetDM(ts, &dm);
278: DMGetApplicationContext(dm, (void **) &user);
279: SetInitialCoordinates(dm);
280: SetInitialConditions(dm, u);
281: return(0);
282: }
284: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
285: {
286: MPI_Comm comm;
287: DM sdm;
288: AppCtx *user;
289: const PetscScalar *u, *coords;
290: PetscScalar *e;
291: PetscReal t;
292: PetscInt dim, Np, p;
293: PetscErrorCode ierr;
296: PetscObjectGetComm((PetscObject) ts, &comm);
297: TSGetDM(ts, &sdm);
298: DMGetApplicationContext(sdm, (void **) &user);
299: DMGetDimension(sdm, &dim);
300: TSGetSolveTime(ts, &t);
301: VecGetArray(E, &e);
302: VecGetArrayRead(U, &u);
303: VecGetLocalSize(U, &Np);
304: Np /= 2*dim;
305: DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
306: for (p = 0; p < Np; ++p) {
307: const PetscScalar *x = &u[(p*2+0)*dim];
308: const PetscScalar *v = &u[(p*2+1)*dim];
309: const PetscReal x0 = p+1.;
310: const PetscReal omega = PetscSqrtReal(1000./(p+1.))/x0;
311: const PetscReal xe[3] = { x0*PetscCosReal(omega*t), x0*PetscSinReal(omega*t), 0.0};
312: const PetscReal ve[3] = {-x0*omega*PetscSinReal(omega*t), x0*omega*PetscCosReal(omega*t), 0.0};
313: PetscInt d;
315: for (d = 0; d < dim; ++d) {
316: e[(p*2+0)*dim+d] = x[d] - xe[d];
317: e[(p*2+1)*dim+d] = v[d] - ve[d];
318: }
319: if (user->error) {PetscPrintf(comm, "p%D error [%.2g %.2g] sol [(%.6lf %.6lf) (%.6lf %.6lf)] exact [(%.6lf %.6lf) (%.6lf %.6lf)] energy/exact energy %g / %g\n", 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))));}
320: }
321: DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
322: VecRestoreArrayRead(U, &u);
323: VecRestoreArray(E, &e);
324: return(0);
325: }
327: int main(int argc, char **argv)
328: {
329: TS ts;
330: TSConvergedReason reason;
331: DM dm, sw;
332: Vec u;
333: IS is1, is2;
334: PetscInt *idx1, *idx2;
335: MPI_Comm comm;
336: AppCtx user;
337: const PetscScalar *endVals;
338: PetscReal ftime = .1;
339: PetscInt locSize, dim, d, Np, p, steps;
340: PetscErrorCode ierr;
342: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
343: comm = PETSC_COMM_WORLD;
344: ProcessOptions(comm, &user);
346: CreateMesh(comm, &dm, &user);
347: CreateParticles(dm, &sw, &user);
348: DMSetApplicationContext(sw, &user);
350: TSCreate(comm, &ts);
351: TSSetType(ts, TSBASICSYMPLECTIC);
352: TSSetDM(ts, sw);
353: TSSetMaxTime(ts, ftime);
354: TSSetTimeStep(ts, 0.0001);
355: TSSetMaxSteps(ts, 10);
356: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
357: TSSetTime(ts, 0.0);
358: TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user);
360: DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u);
361: DMGetDimension(sw, &dim);
362: VecGetLocalSize(u, &locSize);
363: Np = locSize/(2*dim);
364: PetscMalloc1(locSize/2, &idx1);
365: PetscMalloc1(locSize/2, &idx2);
366: for (p = 0; p < Np; ++p) {
367: for (d = 0; d < dim; ++d) {
368: idx1[p*dim+d] = (p*2+0)*dim + d;
369: idx2[p*dim+d] = (p*2+1)*dim + d;
370: }
371: }
372: ISCreateGeneral(comm, locSize/2, idx1, PETSC_OWN_POINTER, &is1);
373: ISCreateGeneral(comm, locSize/2, idx2, PETSC_OWN_POINTER, &is2);
374: TSRHSSplitSetIS(ts, "position", is1);
375: TSRHSSplitSetIS(ts, "momentum", is2);
376: ISDestroy(&is1);
377: ISDestroy(&is2);
378: TSRHSSplitSetRHSFunction(ts,"position",NULL,RHSFunction1,&user);
379: TSRHSSplitSetRHSFunction(ts,"momentum",NULL,RHSFunction2,&user);
381: TSSetFromOptions(ts);
382: TSSetComputeInitialCondition(ts, InitializeSolve);
383: TSSetComputeExactError(ts, ComputeError);
384: TSComputeInitialCondition(ts, u);
385: VecViewFromOptions(u, NULL, "-init_view");
386: TSSolve(ts, u);
387: TSGetSolveTime(ts, &ftime);
388: TSGetConvergedReason(ts, &reason);
389: TSGetStepNumber(ts, &steps);
390: PetscPrintf(comm,"%s at time %g after %D steps\n", TSConvergedReasons[reason], (double) ftime, steps);
392: VecGetArrayRead(u, &endVals);
393: for (p = 0; p < Np; ++p) {
394: const PetscReal norm = DMPlex_NormD_Internal(dim, &endVals[(p*2 + 1)*dim]);
395: PetscPrintf(comm, "Particle %D initial Energy: %g Final Energy: %g\n", p, (double) (0.5*(1000./(p+1))), (double) 0.5*PetscSqr(norm));
396: }
397: VecRestoreArrayRead(u, &endVals);
398: DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u);
399: TSDestroy(&ts);
400: DMDestroy(&sw);
401: DMDestroy(&dm);
402: PetscFinalize();
403: return ierr;
404: }
407: /*TEST
409: build:
410: requires: triangle !single !complex
411: test:
412: suffix: bsi1
413: args: -dim 2 -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
414: test:
415: suffix: bsi2
416: args: -dim 2 -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
417: test:
418: suffix: euler
419: args: -dim 2 -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
421: TEST*/