Actual source code: ex4.c
petsc-3.13.4 2020-08-01
1: static char help[] = "Example of simple hamiltonian system (harmonic oscillator) with particles and a basic symplectic integrator\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: PetscReal omega; /* Oscillation frequency omega */
14: PetscInt particlesPerCell; /* The number of partices per cell */
15: PetscReal momentTol; /* Tolerance for checking moment conservation */
16: PetscBool monitor; /* Flag for using the TS monitor */
17: PetscBool error; /* Flag for printing the error */
18: PetscInt ostep; /* print the energy at each ostep time steps */
19: } AppCtx;
21: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
22: {
26: options->dim = 2;
27: options->simplex = PETSC_TRUE;
28: options->monitor = PETSC_FALSE;
29: options->error = PETSC_FALSE;
30: options->particlesPerCell = 1;
31: options->momentTol = 100.0*PETSC_MACHINE_EPSILON;
32: options->omega = 64.0;
33: options->ostep = 100;
35: PetscStrcpy(options->filename, "");
37: PetscOptionsBegin(comm, "", "Harmonic Oscillator Options", "DMPLEX");
38: PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL);
39: PetscOptionsInt("-dim", "The topological mesh dimension", "ex4.c", options->dim, &options->dim, NULL);
40: PetscOptionsBool("-monitor", "Flag to use the TS monitor", "ex4.c", options->monitor, &options->monitor, NULL);
41: PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL);
42: PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex4.c", options->simplex, &options->simplex, NULL);
43: PetscOptionsString("-mesh", "Name of the mesh filename if any", "ex4.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);
44: PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex4.c", options->particlesPerCell, &options->particlesPerCell, NULL);
45: PetscOptionsReal("-omega", "Oscillator frequency", "ex4.c", options->omega, &options->omega, PETSC_NULL);
46: PetscOptionsEnd();
48: return(0);
49: }
51: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
52: {
53: PetscBool flg;
57: PetscStrcmp(user->filename, "", &flg);
58: if (flg) {
59: DMPlexCreateBoxMesh(comm, user->dim, user->simplex, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);
60: } else {
61: DMPlexCreateFromFile(comm, user->filename, PETSC_TRUE, dm);
62: DMGetDimension(*dm, &user->dim);
63: }
64: {
65: DM distributedMesh = NULL;
67: DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
68: if (distributedMesh) {
69: DMDestroy(dm);
70: *dm = distributedMesh;
71: }
72: }
73: DMLocalizeCoordinates(*dm); /* needed for periodic */
74: DMSetFromOptions(*dm);
75: PetscObjectSetName((PetscObject) *dm, "Mesh");
76: DMViewFromOptions(*dm, NULL, "-dm_view");
77: return(0);
78: }
80: static PetscErrorCode SetInitialCoordinates(DM dmSw)
81: {
82: DM dm;
83: AppCtx *user;
84: PetscRandom rnd;
85: PetscBool simplex;
86: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
87: PetscInt dim, d, cStart, cEnd, c, Np, p;
91: PetscRandomCreate(PetscObjectComm((PetscObject) dmSw), &rnd);
92: PetscRandomSetInterval(rnd, -1.0, 1.0);
93: PetscRandomSetFromOptions(rnd);
95: DMGetApplicationContext(dmSw, (void **) &user);
96: simplex = user->simplex;
97: Np = user->particlesPerCell;
98: DMSwarmGetCellDM(dmSw, &dm);
99: DMGetDimension(dm, &dim);
100: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
101: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
102: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
103: DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
104: for (c = cStart; c < cEnd; ++c) {
105: if (Np == 1) {
106: DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
107: for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
108: } else {
109: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
110: for (p = 0; p < Np; ++p) {
111: const PetscInt n = c*Np + p;
112: PetscReal sum = 0.0, refcoords[3];
114: for (d = 0; d < dim; ++d) {
115: PetscRandomGetValueReal(rnd, &refcoords[d]);
116: sum += refcoords[d];
117: }
118: if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
119: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
120: }
121: }
122: }
123: DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
124: PetscFree5(centroid, xi0, v0, J, invJ);
125: PetscRandomDestroy(&rnd);
126: return(0);
127: }
129: static PetscErrorCode SetInitialConditions(DM dmSw, Vec u)
130: {
131: DM dm;
132: AppCtx *user;
133: PetscReal *coords;
134: PetscScalar *initialConditions;
135: PetscInt dim, cStart, cEnd, c, Np, p;
139: DMGetApplicationContext(dmSw, (void **) &user);
140: Np = user->particlesPerCell;
141: DMSwarmGetCellDM(dmSw, &dm);
142: DMGetDimension(dm, &dim);
143: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
144: DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
145: VecGetArray(u, &initialConditions);
146: for (c = cStart; c < cEnd; ++c) {
147: for (p = 0; p < Np; ++p) {
148: const PetscInt n = c*Np + p;
150: initialConditions[n*2+0] = DMPlex_NormD_Internal(dim, &coords[n*dim]);
151: initialConditions[n*2+1] = 0.0;
152: }
153: }
154: VecRestoreArray(u, &initialConditions);
155: DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
156: return(0);
157: }
159: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
160: {
161: PetscInt *cellid;
162: PetscInt dim, cStart, cEnd, c, Np = user->particlesPerCell, p;
166: DMGetDimension(dm, &dim);
167: DMCreate(PetscObjectComm((PetscObject) dm), sw);
168: DMSetType(*sw, DMSWARM);
169: DMSetDimension(*sw, dim);
171: DMSwarmSetType(*sw, DMSWARM_PIC);
172: DMSwarmSetCellDM(*sw, dm);
173: DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2, PETSC_REAL);
174: DMSwarmFinalizeFieldRegister(*sw);
175: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
176: DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0);
177: DMSetFromOptions(*sw);
178: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
179: for (c = cStart; c < cEnd; ++c) {
180: for (p = 0; p < Np; ++p) {
181: const PetscInt n = c*Np + p;
183: cellid[n] = c;
184: }
185: }
186: DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
187: PetscObjectSetName((PetscObject) *sw, "Particles");
188: DMViewFromOptions(*sw, NULL, "-sw_view");
189: return(0);
190: }
192: /* Create particle RHS Functions */
193: static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
194: {
195: const PetscScalar *v;
196: PetscScalar *xres;
197: PetscInt Np, p;
198: PetscErrorCode ierr;
201: VecGetArray(Xres, &xres);
202: VecGetArrayRead(V, &v);
203: VecGetLocalSize(Xres, &Np);
204: for (p = 0; p < Np; ++p) {
205: xres[p] = v[p];
206: }
207: VecRestoreArrayRead(V, &v);
208: VecRestoreArray(Xres, &xres);
209: return(0);
210: }
212: static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
213: {
214: AppCtx *user = (AppCtx *)ctx;
215: const PetscScalar *x;
216: PetscScalar *vres;
217: PetscInt Np, p;
218: PetscErrorCode ierr;
221: VecGetArray(Vres, &vres);
222: VecGetArrayRead(X, &x);
223: VecGetLocalSize(Vres, &Np);
224: for (p = 0; p < Np; ++p) {
225: vres[p] = -PetscSqr(user->omega)*x[p];
226: }
227: VecRestoreArrayRead(X, &x);
228: VecRestoreArray(Vres, &vres);
229: return(0);
230: }
232: static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t, Vec U, Vec R, void *ctx)
233: {
234: AppCtx *user = (AppCtx *) ctx;
235: DM dm;
236: const PetscScalar *u;
237: PetscScalar *r;
238: PetscInt Np, p;
239: PetscErrorCode ierr;
242: TSGetDM(ts, &dm);
243: VecGetArray(R, &r);
244: VecGetArrayRead(U, &u);
245: VecGetLocalSize(U, &Np);
246: Np /= 2;
247: for (p = 0; p < Np; ++p) {
248: r[p*2+0] = u[p*2+1];
249: r[p*2+1] = -PetscSqr(user->omega)*u[p*2+0];
250: }
251: VecRestoreArrayRead(U, &u);
252: VecRestoreArray(R, &r);
253: return(0);
254: }
256: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
257: {
258: AppCtx *user = (AppCtx *) ctx;
259: const PetscReal omega = user->omega;
260: const PetscScalar *u;
261: MPI_Comm comm;
262: PetscReal dt;
263: PetscInt Np, p;
264: PetscErrorCode ierr;
267: if (step%user->ostep == 0) {
268: PetscObjectGetComm((PetscObject) ts, &comm);
269: if (!step) {PetscPrintf(comm, "Time Step Part Energy Mod Energy\n");}
270: TSGetTimeStep(ts, &dt);
271: VecGetArrayRead(U, &u);
272: VecGetLocalSize(U, &Np);
273: Np /= 2;
274: for (p = 0; p < Np; ++p) {
275: const PetscReal x = PetscRealPart(u[p*2+0]);
276: const PetscReal v = PetscRealPart(u[p*2+1]);
277: const PetscReal E = 0.5*(v*v + PetscSqr(omega)*x*x);
278: const PetscReal mE = 0.5*(v*v + PetscSqr(omega)*x*x - PetscSqr(omega)*dt*x*v);
280: PetscPrintf(comm, "%.6lf %4D %4D %10.4lf %10.4lf\n", t, step, p, (double) E, (double) mE);
281: }
282: VecRestoreArrayRead(U, &u);
283: }
284: return(0);
285: }
287: static PetscErrorCode InitializeSolve(TS ts, Vec u)
288: {
289: DM dm;
290: AppCtx *user;
294: TSGetDM(ts, &dm);
295: DMGetApplicationContext(dm, (void **) &user);
296: SetInitialCoordinates(dm);
297: SetInitialConditions(dm, u);
298: return(0);
299: }
301: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
302: {
303: MPI_Comm comm;
304: DM sdm;
305: AppCtx *user;
306: const PetscScalar *u, *coords;
307: PetscScalar *e;
308: PetscReal t, omega;
309: PetscInt dim, Np, p;
310: PetscErrorCode ierr;
313: PetscObjectGetComm((PetscObject) ts, &comm);
314: TSGetDM(ts, &sdm);
315: DMGetApplicationContext(sdm, (void **) &user);
316: omega = user->omega;
317: DMGetDimension(sdm, &dim);
318: TSGetSolveTime(ts, &t);
319: VecGetArray(E, &e);
320: VecGetArrayRead(U, &u);
321: VecGetLocalSize(U, &Np);
322: Np /= 2;
323: DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
324: for (p = 0; p < Np; ++p) {
325: const PetscReal x = PetscRealPart(u[p*2+0]);
326: const PetscReal v = PetscRealPart(u[p*2+1]);
327: const PetscReal x0 = DMPlex_NormD_Internal(dim, &coords[p*dim]);
328: const PetscReal ex = x0*PetscCosReal(omega*t);
329: const PetscReal ev = -x0*omega*PetscSinReal(omega*t);
331: if (user->error) {PetscPrintf(comm, "p%D error [%.2g %.2g] sol [%.6lf %.6lf] exact [%.6lf %.6lf] energy/exact energy %g / %g\n", p, (double) PetscAbsReal(x-ex), (double) PetscAbsReal(v-ev), (double) x, (double) v, (double) ex, (double) ev, 0.5*(v*v + PetscSqr(omega)*x*x), (double) 0.5*PetscSqr(omega*x0));}
332: e[p*2+0] = x - ex;
333: e[p*2+1] = v - ev;
334: }
335: DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
336: VecRestoreArrayRead(U, &u);
337: VecRestoreArray(E, &e);
338: return(0);
339: }
341: int main(int argc,char **argv)
342: {
343: TS ts; /* nonlinear solver */
344: DM dm, sw; /* Mesh and particle managers */
345: Vec u; /* swarm vector */
346: IS is1, is2;
347: PetscInt n;
348: MPI_Comm comm;
349: AppCtx user;
352: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
353: comm = PETSC_COMM_WORLD;
354: ProcessOptions(comm, &user);
356: CreateMesh(comm, &dm, &user);
357: CreateParticles(dm, &sw, &user);
358: DMSetApplicationContext(sw, &user);
360: TSCreate(comm, &ts);
361: TSSetType(ts, TSBASICSYMPLECTIC);
362: TSSetDM(ts, sw);
363: TSSetMaxTime(ts, 0.1);
364: TSSetTimeStep(ts, 0.00001);
365: TSSetMaxSteps(ts, 100);
366: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
367: if (user.monitor) {TSMonitorSet(ts, Monitor, &user, NULL);}
368: TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user);
370: DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u);
371: VecGetLocalSize(u, &n);
372: ISCreateStride(comm, n/2, 0, 2, &is1);
373: ISCreateStride(comm, n/2, 1, 2, &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: TSSolve(ts, u);
386: DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u);
388: TSDestroy(&ts);
389: DMDestroy(&sw);
390: DMDestroy(&dm);
391: PetscFinalize();
392: return ierr;
393: }
395: /*TEST
397: build:
398: requires: triangle !single !complex
399: test:
400: suffix: 1
401: args: -dm_plex_box_faces 1,1 -ts_basicsymplectic_type 1 -ts_convergence_estimate -convest_num_refine 2 -dm_view -sw_view -monitor -output_step 50 -error
402: test:
403: suffix: 2
404: args: -dm_plex_box_faces 1,1 -ts_basicsymplectic_type 2 -ts_convergence_estimate -convest_num_refine 2 -dm_view -sw_view -monitor -output_step 50 -error
405: test:
406: suffix: 3
407: args: -dm_plex_box_faces 1,1 -ts_basicsymplectic_type 3 -ts_convergence_estimate -convest_num_refine 2 -ts_dt 0.0001 -dm_view -sw_view -monitor -output_step 50 -error
408: test:
409: suffix: 4
410: args: -dm_plex_box_faces 1,1 -ts_basicsymplectic_type 4 -ts_convergence_estimate -convest_num_refine 2 -ts_dt 0.0001 -dm_view -sw_view -monitor -output_step 50 -error
412: TEST*/