Actual source code: ex5.c
1: static char help[] = "Vlasov example of central orbits\n";
3: /*
4: To visualize the orbit, we can used
6: -ts_monitor_sp_swarm -ts_monitor_sp_swarm_retain -1 -ts_monitor_sp_swarm_phase 0 -draw_size 500,500
8: and we probably want it to run fast and not check convergence
10: -convest_num_refine 0 -ts_dt 0.01 -ts_max_steps 100 -ts_max_time 100 -output_step 10
11: */
13: #include <petscts.h>
14: #include <petscdmplex.h>
15: #include <petscdmswarm.h>
16: #include <petsc/private/dmpleximpl.h>
18: PETSC_EXTERN PetscErrorCode circleSingleX(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
19: PETSC_EXTERN PetscErrorCode circleSingleV(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
20: PETSC_EXTERN PetscErrorCode circleMultipleX(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
21: PETSC_EXTERN PetscErrorCode circleMultipleV(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
23: typedef struct {
24: PetscBool error; /* Flag for printing the error */
25: PetscInt ostep; /* print the energy at each ostep time steps */
26: } AppCtx;
28: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
29: {
31: options->error = PETSC_FALSE;
32: options->ostep = 100;
34: PetscOptionsBegin(comm, "", "Central Orbit Options", "DMSWARM");
35: PetscOptionsBool("-error", "Flag to print the error", "ex5.c", options->error, &options->error, NULL);
36: PetscOptionsInt("-output_step", "Number of time steps between output", "ex5.c", options->ostep, &options->ostep, PETSC_NULL);
37: PetscOptionsEnd();
38: return 0;
39: }
41: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
42: {
44: DMCreate(comm, dm);
45: DMSetType(*dm, DMPLEX);
46: DMSetFromOptions(*dm);
47: DMViewFromOptions(*dm, NULL, "-dm_view");
48: return 0;
49: }
51: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
52: {
53: PetscReal v0[1] = {1.};
54: PetscInt dim;
57: DMGetDimension(dm, &dim);
58: DMCreate(PetscObjectComm((PetscObject)dm), sw);
59: DMSetType(*sw, DMSWARM);
60: DMSetDimension(*sw, dim);
61: DMSwarmSetType(*sw, DMSWARM_PIC);
62: DMSwarmSetCellDM(*sw, dm);
63: DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR);
64: DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL);
65: DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT);
66: DMSwarmRegisterPetscDatatypeField(*sw, "initCoordinates", dim, PETSC_REAL);
67: DMSwarmRegisterPetscDatatypeField(*sw, "initVelocity", dim, PETSC_REAL);
68: DMSwarmFinalizeFieldRegister(*sw);
69: DMSwarmComputeLocalSizeFromOptions(*sw);
70: DMSwarmInitializeCoordinates(*sw);
71: DMSwarmInitializeVelocitiesFromOptions(*sw, v0);
72: DMSetFromOptions(*sw);
73: DMSetApplicationContext(*sw, user);
74: PetscObjectSetName((PetscObject)*sw, "Particles");
75: DMViewFromOptions(*sw, NULL, "-sw_view");
76: {
77: Vec gc, gc0, gv, gv0;
79: DMSwarmCreateGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc);
80: DMSwarmCreateGlobalVectorFromField(*sw, "initCoordinates", &gc0);
81: VecCopy(gc, gc0);
82: DMSwarmDestroyGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc);
83: DMSwarmDestroyGlobalVectorFromField(*sw, "initCoordinates", &gc0);
84: DMSwarmCreateGlobalVectorFromField(*sw, "velocity", &gv);
85: DMSwarmCreateGlobalVectorFromField(*sw, "initVelocity", &gv0);
86: VecCopy(gv, gv0);
87: DMSwarmDestroyGlobalVectorFromField(*sw, "velocity", &gv);
88: DMSwarmDestroyGlobalVectorFromField(*sw, "initVelocity", &gv0);
89: }
90: return 0;
91: }
93: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx)
94: {
95: DM sw;
96: const PetscReal *coords, *vel;
97: const PetscScalar *u;
98: PetscScalar *g;
99: PetscInt dim, d, Np, p;
102: TSGetDM(ts, &sw);
103: DMGetDimension(sw, &dim);
104: DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
105: DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel);
106: VecGetLocalSize(U, &Np);
107: VecGetArrayRead(U, &u);
108: VecGetArray(G, &g);
109: Np /= 2 * dim;
110: for (p = 0; p < Np; ++p) {
111: const PetscReal x0 = coords[p * dim + 0];
112: const PetscReal vy0 = vel[p * dim + 1];
113: const PetscReal omega = vy0 / x0;
115: for (d = 0; d < dim; ++d) {
116: g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d];
117: g[(p * 2 + 1) * dim + d] = -PetscSqr(omega) * u[(p * 2 + 0) * dim + d];
118: }
119: }
120: DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
121: DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel);
122: VecRestoreArrayRead(U, &u);
123: VecRestoreArray(G, &g);
124: return 0;
125: }
127: /* J_{ij} = dF_i/dx_j
128: J_p = ( 0 1)
129: (-w^2 0)
130: */
131: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, void *ctx)
132: {
133: DM sw;
134: const PetscReal *coords, *vel;
135: PetscInt dim, d, Np, p, rStart;
138: TSGetDM(ts, &sw);
139: DMGetDimension(sw, &dim);
140: VecGetLocalSize(U, &Np);
141: MatGetOwnershipRange(J, &rStart, NULL);
142: DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
143: DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel);
144: Np /= 2 * dim;
145: for (p = 0; p < Np; ++p) {
146: const PetscReal x0 = coords[p * dim + 0];
147: const PetscReal vy0 = vel[p * dim + 1];
148: const PetscReal omega = vy0 / x0;
149: PetscScalar vals[4] = {0., 1., -PetscSqr(omega), 0.};
151: for (d = 0; d < dim; ++d) {
152: const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
153: MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES);
154: }
155: }
156: DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
157: DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel);
158: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
159: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
160: return 0;
161: }
163: static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
164: {
165: const PetscScalar *v;
166: PetscScalar *xres;
167: PetscInt Np, p;
170: VecGetLocalSize(Xres, &Np);
171: VecGetArrayRead(V, &v);
172: VecGetArray(Xres, &xres);
173: for (p = 0; p < Np; ++p) xres[p] = v[p];
174: VecRestoreArrayRead(V, &v);
175: VecRestoreArray(Xres, &xres);
176: return 0;
177: }
179: static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
180: {
181: DM sw;
182: const PetscScalar *x;
183: const PetscReal *coords, *vel;
184: PetscScalar *vres;
185: PetscInt Np, p, dim, d;
188: TSGetDM(ts, &sw);
189: DMGetDimension(sw, &dim);
190: DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
191: DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel);
192: VecGetLocalSize(Vres, &Np);
193: VecGetArrayRead(X, &x);
194: VecGetArray(Vres, &vres);
196: Np /= dim;
197: for (p = 0; p < Np; ++p) {
198: const PetscReal x0 = coords[p * dim + 0];
199: const PetscReal vy0 = vel[p * dim + 1];
200: const PetscReal omega = vy0 / x0;
202: for (d = 0; d < dim; ++d) vres[p * dim + d] = -PetscSqr(omega) * x[p * dim + d];
203: }
204: VecRestoreArrayRead(X, &x);
205: VecRestoreArray(Vres, &vres);
206: DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
207: DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel);
208: return 0;
209: }
211: static PetscErrorCode CreateSolution(TS ts)
212: {
213: DM sw;
214: Vec u;
215: PetscInt dim, Np;
217: TSGetDM(ts, &sw);
218: DMGetDimension(sw, &dim);
219: DMSwarmGetLocalSize(sw, &Np);
220: VecCreate(PETSC_COMM_WORLD, &u);
221: VecSetBlockSize(u, dim);
222: VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE);
223: VecSetUp(u);
224: TSSetSolution(ts, u);
225: VecDestroy(&u);
226: return 0;
227: }
229: static PetscErrorCode SetProblem(TS ts)
230: {
231: AppCtx *user;
232: DM sw;
234: TSGetDM(ts, &sw);
235: DMGetApplicationContext(sw, (void **)&user);
236: // Define unified system for (X, V)
237: {
238: Mat J;
239: PetscInt dim, Np;
241: DMGetDimension(sw, &dim);
242: DMSwarmGetLocalSize(sw, &Np);
243: MatCreate(PETSC_COMM_WORLD, &J);
244: MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE);
245: MatSetBlockSize(J, 2 * dim);
246: MatSetFromOptions(J);
247: MatSetUp(J);
248: TSSetRHSFunction(ts, NULL, RHSFunction, &user);
249: TSSetRHSJacobian(ts, J, J, RHSJacobian, &user);
250: MatDestroy(&J);
251: }
252: // Define split system for X and V
253: {
254: Vec u;
255: IS isx, isv, istmp;
256: const PetscInt *idx;
257: PetscInt dim, Np, rstart;
259: TSGetSolution(ts, &u);
260: DMGetDimension(sw, &dim);
261: DMSwarmGetLocalSize(sw, &Np);
262: VecGetOwnershipRange(u, &rstart, NULL);
263: ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 0, 2, &istmp);
264: ISGetIndices(istmp, &idx);
265: ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isx);
266: ISRestoreIndices(istmp, &idx);
267: ISDestroy(&istmp);
268: ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 1, 2, &istmp);
269: ISGetIndices(istmp, &idx);
270: ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isv);
271: ISRestoreIndices(istmp, &idx);
272: ISDestroy(&istmp);
273: TSRHSSplitSetIS(ts, "position", isx);
274: TSRHSSplitSetIS(ts, "momentum", isv);
275: ISDestroy(&isx);
276: ISDestroy(&isv);
277: TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, &user);
278: TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, &user);
279: }
280: // Define symplectic formulation U_t = S . G, where G = grad F
281: {
282: //TSDiscGradSetFormulation(ts, RHSJacobianS, RHSObjectiveF, RHSFunctionG, &user);
283: }
284: return 0;
285: }
287: static PetscErrorCode DMSwarmTSRedistribute(TS ts)
288: {
289: DM sw;
290: Vec u;
291: PetscReal t, maxt, dt;
292: PetscInt n, maxn;
294: TSGetDM(ts, &sw);
295: TSGetTime(ts, &t);
296: TSGetMaxTime(ts, &maxt);
297: TSGetTimeStep(ts, &dt);
298: TSGetStepNumber(ts, &n);
299: TSGetMaxSteps(ts, &maxn);
301: TSReset(ts);
302: TSSetDM(ts, sw);
303: /* TODO Check whether TS was set from options */
304: TSSetFromOptions(ts);
305: TSSetTime(ts, t);
306: TSSetMaxTime(ts, maxt);
307: TSSetTimeStep(ts, dt);
308: TSSetStepNumber(ts, n);
309: TSSetMaxSteps(ts, maxn);
311: CreateSolution(ts);
312: SetProblem(ts);
313: TSGetSolution(ts, &u);
314: return 0;
315: }
317: PetscErrorCode circleSingleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
318: {
319: x[0] = p + 1.;
320: x[1] = 0.;
321: return 0;
322: }
324: PetscErrorCode circleSingleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *ctx)
325: {
326: v[0] = 0.;
327: v[1] = PetscSqrtReal(1000. / (p + 1.));
328: return 0;
329: }
331: /* Put 5 particles into each circle */
332: PetscErrorCode circleMultipleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
333: {
334: const PetscInt n = 5;
335: const PetscReal r0 = (p / n) + 1.;
336: const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
338: x[0] = r0 * PetscCosReal(th0);
339: x[1] = r0 * PetscSinReal(th0);
340: return 0;
341: }
343: /* Put 5 particles into each circle */
344: PetscErrorCode circleMultipleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *ctx)
345: {
346: const PetscInt n = 5;
347: const PetscReal r0 = (p / n) + 1.;
348: const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
349: const PetscReal omega = PetscSqrtReal(1000. / r0) / r0;
351: v[0] = -r0 * omega * PetscSinReal(th0);
352: v[1] = r0 * omega * PetscCosReal(th0);
353: return 0;
354: }
356: /*
357: InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values.
359: Input Parameters:
360: + ts - The TS
361: - useInitial - Flag to also set the initial conditions to the current coordinates and velocities and setup the problem
363: Output Parameters:
364: . u - The initialized solution vector
366: Level: advanced
368: .seealso: InitializeSolve()
369: */
370: static PetscErrorCode InitializeSolveAndSwarm(TS ts, PetscBool useInitial)
371: {
372: DM sw;
373: Vec u, gc, gv, gc0, gv0;
374: IS isx, isv;
375: AppCtx *user;
378: TSGetDM(ts, &sw);
379: DMGetApplicationContext(sw, &user);
380: if (useInitial) {
381: PetscReal v0[1] = {1.};
383: DMSwarmInitializeCoordinates(sw);
384: DMSwarmInitializeVelocitiesFromOptions(sw, v0);
385: DMSwarmMigrate(sw, PETSC_TRUE);
386: DMSwarmTSRedistribute(ts);
387: }
388: TSGetSolution(ts, &u);
389: TSRHSSplitGetIS(ts, "position", &isx);
390: TSRHSSplitGetIS(ts, "momentum", &isv);
391: DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc);
392: DMSwarmCreateGlobalVectorFromField(sw, "initCoordinates", &gc0);
393: if (useInitial) VecCopy(gc, gc0);
394: VecISCopy(u, isx, SCATTER_FORWARD, gc);
395: DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc);
396: DMSwarmDestroyGlobalVectorFromField(sw, "initCoordinates", &gc0);
397: DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv);
398: DMSwarmCreateGlobalVectorFromField(sw, "initVelocity", &gv0);
399: if (useInitial) VecCopy(gv, gv0);
400: VecISCopy(u, isv, SCATTER_FORWARD, gv);
401: DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv);
402: DMSwarmDestroyGlobalVectorFromField(sw, "initVelocity", &gv0);
403: return 0;
404: }
406: static PetscErrorCode InitializeSolve(TS ts, Vec u)
407: {
408: TSSetSolution(ts, u);
409: InitializeSolveAndSwarm(ts, PETSC_TRUE);
410: return 0;
411: }
413: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
414: {
415: MPI_Comm comm;
416: DM sw;
417: AppCtx *user;
418: const PetscScalar *u;
419: const PetscReal *coords, *vel;
420: PetscScalar *e;
421: PetscReal t;
422: PetscInt dim, Np, p;
425: PetscObjectGetComm((PetscObject)ts, &comm);
426: TSGetDM(ts, &sw);
427: DMGetApplicationContext(sw, &user);
428: DMGetDimension(sw, &dim);
429: TSGetSolveTime(ts, &t);
430: VecGetArray(E, &e);
431: VecGetArrayRead(U, &u);
432: VecGetLocalSize(U, &Np);
433: DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
434: DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel);
435: Np /= 2 * dim;
436: for (p = 0; p < Np; ++p) {
437: /* TODO generalize initial conditions and project into plane instead of assuming x-y */
438: const PetscReal r0 = DMPlex_NormD_Internal(dim, &coords[p * dim]);
439: const PetscReal th0 = PetscAtan2Real(coords[p * dim + 1], coords[p * dim + 0]);
440: const PetscReal v0 = DMPlex_NormD_Internal(dim, &vel[p * dim]);
441: const PetscReal omega = v0 / r0;
442: const PetscReal ct = PetscCosReal(omega * t + th0);
443: const PetscReal st = PetscSinReal(omega * t + th0);
444: const PetscScalar *x = &u[(p * 2 + 0) * dim];
445: const PetscScalar *v = &u[(p * 2 + 1) * dim];
446: const PetscReal xe[3] = {r0 * ct, r0 * st, 0.0};
447: const PetscReal ve[3] = {-v0 * st, v0 * ct, 0.0};
448: PetscInt d;
450: for (d = 0; d < dim; ++d) {
451: e[(p * 2 + 0) * dim + d] = x[d] - xe[d];
452: e[(p * 2 + 1) * dim + d] = v[d] - ve[d];
453: }
454: if (user->error) {
455: const PetscReal en = 0.5 * DMPlex_DotRealD_Internal(dim, v, v);
456: const PetscReal exen = 0.5 * PetscSqr(v0);
457: PetscPrintf(comm, "t %.4g: p%" PetscInt_FMT " error [%.2g %.2g] sol [(%.6lf %.6lf) (%.6lf %.6lf)] exact [(%.6lf %.6lf) (%.6lf %.6lf)] energy/exact energy %g / %g (%.10lf%%)\n", (double)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], (double)en, (double)exen, (double)(PetscAbsReal(exen - en) * 100. / exen));
458: }
459: }
460: DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
461: DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel);
462: VecRestoreArrayRead(U, &u);
463: VecRestoreArray(E, &e);
464: return 0;
465: }
467: static PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
468: {
469: const PetscInt ostep = ((AppCtx *)ctx)->ostep;
470: DM sw;
471: const PetscScalar *u;
472: PetscInt dim, Np, p;
475: if (step % ostep == 0) {
476: TSGetDM(ts, &sw);
477: DMGetDimension(sw, &dim);
478: VecGetArrayRead(U, &u);
479: VecGetLocalSize(U, &Np);
480: Np /= 2 * dim;
481: if (!step) PetscPrintf(PetscObjectComm((PetscObject)ts), "Time Step Part Energy\n");
482: for (p = 0; p < Np; ++p) {
483: const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);
485: PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " %4" PetscInt_FMT " %10.4lf\n", (double)t, step, p, (double)(0.5 * v2));
486: }
487: PetscSynchronizedFlush(PetscObjectComm((PetscObject)ts), NULL);
488: VecRestoreArrayRead(U, &u);
489: }
490: return 0;
491: }
493: static PetscErrorCode MigrateParticles(TS ts)
494: {
495: DM sw;
498: TSGetDM(ts, &sw);
499: DMViewFromOptions(sw, NULL, "-migrate_view_pre");
500: {
501: Vec u, gc, gv;
502: IS isx, isv;
504: TSGetSolution(ts, &u);
505: TSRHSSplitGetIS(ts, "position", &isx);
506: TSRHSSplitGetIS(ts, "momentum", &isv);
507: DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc);
508: VecISCopy(u, isx, SCATTER_REVERSE, gc);
509: DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc);
510: DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv);
511: VecISCopy(u, isv, SCATTER_REVERSE, gv);
512: DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv);
513: }
514: DMSwarmMigrate(sw, PETSC_TRUE);
515: DMSwarmTSRedistribute(ts);
516: InitializeSolveAndSwarm(ts, PETSC_FALSE);
517: //VecViewFromOptions(u, NULL, "-sol_migrate_view");
518: return 0;
519: }
521: int main(int argc, char **argv)
522: {
523: DM dm, sw;
524: TS ts;
525: Vec u;
526: AppCtx user;
529: PetscInitialize(&argc, &argv, NULL, help);
530: ProcessOptions(PETSC_COMM_WORLD, &user);
531: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
532: CreateSwarm(dm, &user, &sw);
533: DMSetApplicationContext(sw, &user);
535: TSCreate(PETSC_COMM_WORLD, &ts);
536: TSSetProblemType(ts, TS_NONLINEAR);
537: TSSetDM(ts, sw);
538: TSSetMaxTime(ts, 0.1);
539: TSSetTimeStep(ts, 0.00001);
540: TSSetMaxSteps(ts, 100);
541: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
542: TSMonitorSet(ts, EnergyMonitor, &user, NULL);
543: TSSetFromOptions(ts);
544: TSSetComputeInitialCondition(ts, InitializeSolve);
545: TSSetComputeExactError(ts, ComputeError);
546: TSSetPostStep(ts, MigrateParticles);
548: CreateSolution(ts);
549: TSGetSolution(ts, &u);
550: TSComputeInitialCondition(ts, u);
551: TSSolve(ts, NULL);
553: TSDestroy(&ts);
554: DMDestroy(&sw);
555: DMDestroy(&dm);
556: PetscFinalize();
557: return 0;
558: }
560: /*TEST
562: build:
563: requires: double !complex
565: testset:
566: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
567: args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
568: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
569: -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
570: -dm_view -output_step 50 -error
571: test:
572: suffix: bsi_2d_1
573: args: -ts_basicsymplectic_type 1
574: test:
575: suffix: bsi_2d_2
576: args: -ts_basicsymplectic_type 2
577: test:
578: suffix: bsi_2d_3
579: args: -ts_basicsymplectic_type 3
580: test:
581: suffix: bsi_2d_4
582: args: -ts_basicsymplectic_type 4 -ts_dt 0.0001
584: testset:
585: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
586: args: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
587: -ts_type theta -ts_theta_theta 0.5 -ts_convergence_estimate -convest_num_refine 2 \
588: -mat_type baij -ksp_error_if_not_converged -pc_type lu \
589: -dm_view -output_step 50 -error
590: test:
591: suffix: im_2d_0
592: args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5
594: testset:
595: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
596: args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 10,10 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 -petscpartitioner_type simple \
597: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
598: -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
599: -dm_view -output_step 50
600: test:
601: suffix: bsi_2d_mesh_1
602: args: -ts_basicsymplectic_type 1 -error
603: test:
604: suffix: bsi_2d_mesh_1_par_2
605: nsize: 2
606: args: -ts_basicsymplectic_type 1
607: test:
608: suffix: bsi_2d_mesh_1_par_3
609: nsize: 3
610: args: -ts_basicsymplectic_type 1
611: test:
612: suffix: bsi_2d_mesh_1_par_4
613: nsize: 4
614: args: -ts_basicsymplectic_type 1 -dm_swarm_num_particles 0,0,2,0
616: testset:
617: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
618: args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 10,10 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
619: -dm_swarm_num_particles 10 -dm_swarm_coordinate_function circleMultipleX -dm_swarm_velocity_function circleMultipleV \
620: -ts_convergence_estimate -convest_num_refine 2 \
621: -dm_view -output_step 50 -error
622: test:
623: suffix: bsi_2d_multiple_1
624: args: -ts_type basicsymplectic -ts_basicsymplectic_type 1
625: test:
626: suffix: bsi_2d_multiple_2
627: args: -ts_type basicsymplectic -ts_basicsymplectic_type 2
628: test:
629: suffix: bsi_2d_multiple_3
630: args: -ts_type basicsymplectic -ts_basicsymplectic_type 3 -ts_dt 0.001
631: test:
632: suffix: im_2d_multiple_0
633: args: -ts_type theta -ts_theta_theta 0.5 \
634: -mat_type baij -ksp_error_if_not_converged -pc_type lu
636: TEST*/