Actual source code: ex6.c
1: static char help[] = "Vlasov-Poisson 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_time_step 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>
17: #include <petscfe.h>
18: #include <petscds.h>
19: #include <petsc/private/petscfeimpl.h>
20: #include <petscksp.h>
21: #include <petscsnes.h>
23: PETSC_EXTERN PetscErrorCode circleSingleX(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
24: PETSC_EXTERN PetscErrorCode circleSingleV(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
25: PETSC_EXTERN PetscErrorCode circleMultipleX(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
26: PETSC_EXTERN PetscErrorCode circleMultipleV(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
28: const char *EMTypes[] = {"primal", "mixed", "coulomb", "none", "EMType", "EM_", NULL};
29: typedef enum {
30: EM_PRIMAL,
31: EM_MIXED,
32: EM_COULOMB,
33: EM_NONE
34: } EMType;
36: typedef struct {
37: PetscBool error; /* Flag for printing the error */
38: PetscInt ostep; /* print the energy at each ostep time steps */
39: PetscReal timeScale; /* Nondimensionalizing time scale */
40: PetscReal sigma; /* Linear charge per box length */
41: EMType em; /* Type of electrostatic model */
42: SNES snes;
43: } AppCtx;
45: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
46: {
47: PetscFunctionBeginUser;
48: options->error = PETSC_FALSE;
49: options->ostep = 100;
50: options->timeScale = 1.0e-6;
51: options->sigma = 1.;
52: options->em = EM_COULOMB;
54: PetscOptionsBegin(comm, "", "Central Orbit Options", "DMSWARM");
55: PetscCall(PetscOptionsBool("-error", "Flag to print the error", "ex6.c", options->error, &options->error, NULL));
56: PetscCall(PetscOptionsInt("-output_step", "Number of time steps between output", "ex6.c", options->ostep, &options->ostep, NULL));
57: PetscCall(PetscOptionsReal("-sigma", "Linear charge per box length", "ex6.c", options->sigma, &options->sigma, NULL));
58: PetscCall(PetscOptionsReal("-timeScale", "Nondimensionalizing time scale", "ex6.c", options->timeScale, &options->timeScale, NULL));
59: PetscCall(PetscOptionsEnum("-em_type", "Type of electrostatic solver", "ex6.c", EMTypes, (PetscEnum)options->em, (PetscEnum *)&options->em, NULL));
60: PetscOptionsEnd();
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
65: {
66: PetscFunctionBeginUser;
67: PetscCall(DMCreate(comm, dm));
68: PetscCall(DMSetType(*dm, DMPLEX));
69: PetscCall(DMSetFromOptions(*dm));
70: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: static void laplacian_f1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
75: {
76: PetscInt d;
77: for (d = 0; d < dim; ++d) f1[d] = u_x[d];
78: }
80: static void laplacian_g3(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
81: {
82: PetscInt d;
83: for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
84: }
86: /*
87: / I grad\ /q\ = /0\
88: \-div 0 / \u/ \f/
89: */
90: static void f0_q(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
91: {
92: for (PetscInt c = 0; c < dim; ++c) f0[c] += u[uOff[0] + c];
93: }
95: static void f1_q(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
96: {
97: for (PetscInt d = 0; d < dim; ++d) f1[d * dim + d] += u[uOff[1]];
98: }
100: /* <t, q> */
101: static void g0_qq(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
102: {
103: for (PetscInt c = 0; c < dim; ++c) g0[c * dim + c] += 1.0;
104: }
106: static void g2_qu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
107: {
108: for (PetscInt d = 0; d < dim; ++d) g2[d * dim + d] += 1.0;
109: }
111: static void g1_uq(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
112: {
113: for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] += 1.0;
114: }
116: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
117: {
118: PetscFE feu, feq;
119: PetscDS ds;
120: PetscBool simplex;
121: PetscInt dim;
123: PetscFunctionBeginUser;
124: PetscCall(DMGetDimension(dm, &dim));
125: PetscCall(DMPlexIsSimplex(dm, &simplex));
126: if (user->em == EM_MIXED) {
127: DMLabel label;
129: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "field_", PETSC_DETERMINE, &feq));
130: PetscCall(PetscObjectSetName((PetscObject)feq, "field"));
131: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "potential_", PETSC_DETERMINE, &feu));
132: PetscCall(PetscObjectSetName((PetscObject)feu, "potential"));
133: PetscCall(PetscFECopyQuadrature(feq, feu));
134: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feq));
135: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)feu));
136: PetscCall(DMCreateDS(dm));
137: PetscCall(PetscFEDestroy(&feu));
138: PetscCall(PetscFEDestroy(&feq));
140: PetscCall(DMGetLabel(dm, "marker", &label));
141: PetscCall(DMGetDS(dm, &ds));
142: PetscCall(PetscDSSetResidual(ds, 0, f0_q, f1_q));
143: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_qq, NULL, NULL, NULL));
144: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_qu, NULL));
145: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_uq, NULL, NULL));
146: } else if (user->em == EM_PRIMAL) {
147: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &feu));
148: PetscCall(PetscObjectSetName((PetscObject)feu, "potential"));
149: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feu));
150: PetscCall(DMCreateDS(dm));
151: PetscCall(PetscFEDestroy(&feu));
152: PetscCall(DMGetDS(dm, &ds));
153: PetscCall(PetscDSSetResidual(ds, 0, NULL, laplacian_f1));
154: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, laplacian_g3));
155: }
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: static PetscErrorCode CreatePoisson(DM dm, AppCtx *user)
160: {
161: SNES snes;
162: Mat J;
163: MatNullSpace nullSpace;
165: PetscFunctionBeginUser;
166: PetscCall(CreateFEM(dm, user));
167: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
168: PetscCall(SNESSetOptionsPrefix(snes, "em_"));
169: PetscCall(SNESSetDM(snes, dm));
170: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, user));
171: PetscCall(SNESSetFromOptions(snes));
173: PetscCall(DMCreateMatrix(dm, &J));
174: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_TRUE, 0, NULL, &nullSpace));
175: PetscCall(MatSetNullSpace(J, nullSpace));
176: PetscCall(MatNullSpaceDestroy(&nullSpace));
177: PetscCall(SNESSetJacobian(snes, J, J, NULL, NULL));
178: PetscCall(MatDestroy(&J));
179: user->snes = snes;
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
184: {
185: PetscReal v0[1] = {1.};
186: PetscInt dim;
188: PetscFunctionBeginUser;
189: PetscCall(DMGetDimension(dm, &dim));
190: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
191: PetscCall(DMSetType(*sw, DMSWARM));
192: PetscCall(DMSetDimension(*sw, dim));
193: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
194: PetscCall(DMSwarmSetCellDM(*sw, dm));
195: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
196: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
197: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT));
198: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "initCoordinates", dim, PETSC_REAL));
199: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "initVelocity", dim, PETSC_REAL));
200: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "E_field", dim, PETSC_REAL));
201: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
202: PetscCall(DMSwarmComputeLocalSizeFromOptions(*sw));
203: PetscCall(DMSwarmInitializeCoordinates(*sw));
204: PetscCall(DMSwarmInitializeVelocitiesFromOptions(*sw, v0));
205: PetscCall(DMSetFromOptions(*sw));
206: PetscCall(DMSetApplicationContext(*sw, user));
207: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
208: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
209: {
210: Vec gc, gc0, gv, gv0;
212: PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc));
213: PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "initCoordinates", &gc0));
214: PetscCall(VecCopy(gc, gc0));
215: PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc));
216: PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "initCoordinates", &gc0));
217: PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "velocity", &gv));
218: PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "initVelocity", &gv0));
219: PetscCall(VecCopy(gv, gv0));
220: PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "velocity", &gv));
221: PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "initVelocity", &gv0));
222: }
223: {
224: const char *fieldnames[2] = {DMSwarmPICField_coor, "velocity"};
226: PetscCall(DMSwarmVectorDefineFields(*sw, 2, fieldnames));
227: }
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: static PetscErrorCode ComputeFieldAtParticles_Coulomb(SNES snes, DM sw, PetscReal E[])
232: {
233: PetscReal *coords;
234: PetscInt dim, d, Np, p, q;
235: PetscMPIInt size;
237: PetscFunctionBegin;
238: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)snes), &size));
239: PetscCheck(size == 1, PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Coulomb code only works in serial");
240: PetscCall(DMGetDimension(sw, &dim));
241: PetscCall(DMSwarmGetLocalSize(sw, &Np));
243: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
244: for (p = 0; p < Np; ++p) {
245: PetscReal *pcoord = &coords[p * dim];
246: PetscReal *pE = &E[p * dim];
247: /* Calculate field at particle p due to particle q */
248: for (q = 0; q < Np; ++q) {
249: PetscReal *qcoord = &coords[q * dim];
250: PetscReal rpq[3], r;
252: if (p == q) continue;
253: for (d = 0; d < dim; ++d) rpq[d] = pcoord[d] - qcoord[d];
254: r = DMPlex_NormD_Internal(dim, rpq);
255: for (d = 0; d < dim; ++d) pE[d] += rpq[d] / PetscPowRealInt(r, 3);
256: }
257: }
258: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: static PetscErrorCode ComputeFieldAtParticles_Primal(SNES snes, DM sw, PetscReal E[])
263: {
264: DM dm;
265: PetscDS ds;
266: PetscFE fe;
267: Mat M_p;
268: Vec phi, locPhi, rho, f;
269: PetscReal *coords, chargeTol = 1e-13;
270: PetscInt dim, d, cStart, cEnd, c, Np;
271: const char **oldFields;
272: PetscInt Nf;
273: const char **tmp;
275: PetscFunctionBegin;
276: PetscCall(DMGetDimension(sw, &dim));
277: PetscCall(DMSwarmGetLocalSize(sw, &Np));
278: PetscCall(SNESGetDM(snes, &dm));
280: PetscCall(DMSwarmVectorGetField(sw, &Nf, &tmp));
281: PetscCall(PetscMalloc1(Nf, &oldFields));
282: for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscStrallocpy(tmp[f], (char **)&oldFields[f]));
283: PetscCall(DMSwarmVectorDefineField(sw, "w_q"));
284: PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
285: PetscCall(DMSwarmVectorDefineFields(sw, Nf, oldFields));
286: for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscFree(oldFields[f]));
287: PetscCall(PetscFree(oldFields));
289: /* Create the charges rho */
290: PetscCall(DMGetGlobalVector(dm, &rho));
291: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
292: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
293: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
294: PetscCall(MatMultTranspose(M_p, f, rho));
295: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
296: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
297: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
298: PetscCall(MatDestroy(&M_p));
299: {
300: PetscScalar sum;
301: PetscInt n;
302: PetscReal phi_0 = 1.; /* (sigma*sigma*sigma)*(timeScale*timeScale)/(m_e*q_e*epsi_0)*/
304: /* Remove constant from rho */
305: PetscCall(VecGetSize(rho, &n));
306: PetscCall(VecSum(rho, &sum));
307: PetscCall(VecShift(rho, -sum / n));
308: PetscCall(VecSum(rho, &sum));
309: PetscCheck(PetscAbsScalar(sum) < chargeTol, PetscObjectComm((PetscObject)sw), PETSC_ERR_PLIB, "Charge should have no DC component: %g", (double)PetscRealPart(sum));
310: /* Nondimensionalize rho */
311: PetscCall(VecScale(rho, phi_0));
312: }
313: PetscCall(VecViewFromOptions(rho, NULL, "-poisson_rho_view"));
315: PetscCall(DMGetGlobalVector(dm, &phi));
316: PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
317: PetscCall(VecSet(phi, 0.0));
318: PetscCall(SNESSolve(snes, rho, phi));
319: PetscCall(DMRestoreGlobalVector(dm, &rho));
320: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
322: PetscCall(DMGetLocalVector(dm, &locPhi));
323: PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
324: PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
325: PetscCall(DMRestoreGlobalVector(dm, &phi));
327: PetscCall(DMGetDS(dm, &ds));
328: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
329: PetscCall(DMSwarmSortGetAccess(sw));
330: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
331: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
332: for (c = cStart; c < cEnd; ++c) {
333: PetscTabulation tab;
334: PetscScalar *clPhi = NULL;
335: PetscReal *pcoord, *refcoord;
336: PetscReal v[3], J[9], invJ[9], detJ;
337: PetscInt *points;
338: PetscInt Ncp, cp;
340: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
341: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
342: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
343: for (cp = 0; cp < Ncp; ++cp)
344: for (d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
345: PetscCall(DMPlexCoordinatesToReference(dm, c, Ncp, pcoord, refcoord));
346: PetscCall(PetscFECreateTabulation(fe, 1, Ncp, refcoord, 1, &tab));
347: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v, J, invJ, &detJ));
348: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
349: for (cp = 0; cp < Ncp; ++cp) {
350: const PetscReal *basisDer = tab->T[1];
351: const PetscInt p = points[cp];
353: for (d = 0; d < dim; ++d) E[p * dim + d] = 0.;
354: PetscCall(PetscFEFreeInterpolateGradient_Static(fe, basisDer, clPhi, dim, invJ, NULL, cp, &E[p * dim]));
355: for (d = 0; d < dim; ++d) E[p * dim + d] *= -1.0;
356: }
357: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
358: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
359: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
360: PetscCall(PetscTabulationDestroy(&tab));
361: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
362: }
363: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
364: PetscCall(DMSwarmSortRestoreAccess(sw));
365: PetscCall(DMRestoreLocalVector(dm, &locPhi));
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: static PetscErrorCode ComputeFieldAtParticles_Mixed(SNES snes, DM sw, PetscReal E[])
370: {
371: DM dm, potential_dm;
372: IS potential_IS;
373: PetscDS ds;
374: PetscFE fe;
375: PetscFEGeom feGeometry;
376: Mat M_p;
377: Vec phi, locPhi, rho, f, temp_rho;
378: PetscQuadrature q;
379: PetscReal *coords, chargeTol = 1e-13;
380: PetscInt dim, d, cStart, cEnd, c, Np, pot_field = 1;
381: const char **oldFields;
382: PetscInt Nf;
383: const char **tmp;
385: PetscFunctionBegin;
386: PetscCall(DMGetDimension(sw, &dim));
387: PetscCall(DMSwarmGetLocalSize(sw, &Np));
389: /* Create the charges rho */
390: PetscCall(SNESGetDM(snes, &dm));
391: PetscCall(DMGetGlobalVector(dm, &rho));
392: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
393: PetscCall(DMCreateSubDM(dm, 1, &pot_field, &potential_IS, &potential_dm));
395: PetscCall(DMSwarmVectorGetField(sw, &Nf, &tmp));
396: PetscCall(PetscMalloc1(Nf, &oldFields));
397: for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscStrallocpy(tmp[f], (char **)&oldFields[f]));
398: PetscCall(DMSwarmVectorDefineField(sw, "w_q"));
399: PetscCall(DMCreateMassMatrix(sw, potential_dm, &M_p));
400: PetscCall(DMSwarmVectorDefineFields(sw, Nf, oldFields));
401: for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscFree(oldFields[f]));
402: PetscCall(PetscFree(oldFields));
404: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
405: PetscCall(DMGetGlobalVector(potential_dm, &temp_rho));
406: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
407: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
408: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
409: PetscCall(MatMultTranspose(M_p, f, temp_rho));
410: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
411: PetscCall(MatDestroy(&M_p));
412: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
413: PetscCall(VecViewFromOptions(rho, NULL, "-poisson_rho_view"));
414: PetscCall(VecISCopy(rho, potential_IS, SCATTER_FORWARD, temp_rho));
415: PetscCall(DMRestoreGlobalVector(potential_dm, &temp_rho));
416: PetscCall(DMDestroy(&potential_dm));
417: PetscCall(ISDestroy(&potential_IS));
418: {
419: PetscScalar sum;
420: PetscInt n;
421: PetscReal phi_0 = 1.; /*(sigma*sigma*sigma)*(timeScale*timeScale)/(m_e*q_e*epsi_0);*/
423: /*Remove constant from rho*/
424: PetscCall(VecGetSize(rho, &n));
425: PetscCall(VecSum(rho, &sum));
426: PetscCall(VecShift(rho, -sum / n));
427: PetscCall(VecSum(rho, &sum));
428: PetscCheck(PetscAbsScalar(sum) < chargeTol, PetscObjectComm((PetscObject)sw), PETSC_ERR_PLIB, "Charge should have no DC component: %g", (double)PetscRealPart(sum));
429: /* Nondimensionalize rho */
430: PetscCall(VecScale(rho, phi_0));
431: }
432: PetscCall(DMGetGlobalVector(dm, &phi));
433: PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
434: PetscCall(VecSet(phi, 0.0));
435: PetscCall(SNESSolve(snes, NULL, phi));
436: PetscCall(DMRestoreGlobalVector(dm, &rho));
437: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
439: PetscCall(DMGetLocalVector(dm, &locPhi));
440: PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
441: PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
442: PetscCall(DMRestoreGlobalVector(dm, &phi));
444: PetscCall(DMGetDS(dm, &ds));
445: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
446: PetscCall(DMSwarmSortGetAccess(sw));
447: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
448: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
449: for (c = cStart; c < cEnd; ++c) {
450: PetscTabulation tab;
451: PetscScalar *clPhi = NULL;
452: PetscReal *pcoord, *refcoord;
453: PetscReal v[3], J[9], invJ[9], detJ;
454: PetscInt *points;
455: PetscInt Ncp, cp;
457: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
458: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
459: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
460: for (cp = 0; cp < Ncp; ++cp)
461: for (d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
462: PetscCall(DMPlexCoordinatesToReference(dm, c, Ncp, pcoord, refcoord));
463: PetscCall(PetscFECreateTabulation(fe, 1, Ncp, refcoord, 1, &tab));
464: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v, J, invJ, &detJ));
465: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
466: for (cp = 0; cp < Ncp; ++cp) {
467: const PetscInt p = points[cp];
469: for (d = 0; d < dim; ++d) E[p * dim + d] = 0.;
470: PetscCall(PetscFEGetQuadrature(fe, &q));
471: PetscCall(PetscFECreateCellGeometry(fe, q, &feGeometry));
472: PetscCall(PetscFEInterpolateAtPoints_Static(fe, tab, clPhi, &feGeometry, cp, &E[p * dim]));
473: PetscCall(PetscFEDestroyCellGeometry(fe, &feGeometry));
474: }
475: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
476: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
477: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
478: PetscCall(PetscTabulationDestroy(&tab));
479: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
480: }
481: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
482: PetscCall(DMSwarmSortRestoreAccess(sw));
483: PetscCall(DMRestoreLocalVector(dm, &locPhi));
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }
487: static PetscErrorCode ComputeFieldAtParticles(SNES snes, DM sw, PetscReal E[])
488: {
489: AppCtx *ctx;
490: PetscInt dim, Np;
492: PetscFunctionBegin;
495: PetscAssertPointer(E, 3);
496: PetscCall(DMGetDimension(sw, &dim));
497: PetscCall(DMSwarmGetLocalSize(sw, &Np));
498: PetscCall(DMGetApplicationContext(sw, &ctx));
499: PetscCall(PetscArrayzero(E, Np * dim));
501: switch (ctx->em) {
502: case EM_PRIMAL:
503: PetscCall(ComputeFieldAtParticles_Primal(snes, sw, E));
504: break;
505: case EM_COULOMB:
506: PetscCall(ComputeFieldAtParticles_Coulomb(snes, sw, E));
507: break;
508: case EM_MIXED:
509: PetscCall(ComputeFieldAtParticles_Mixed(snes, sw, E));
510: break;
511: case EM_NONE:
512: break;
513: default:
514: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No solver for electrostatic model %s", EMTypes[ctx->em]);
515: }
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, PetscCtx ctx)
520: {
521: DM sw;
522: SNES snes = ((AppCtx *)ctx)->snes;
523: const PetscReal *coords, *vel;
524: const PetscScalar *u;
525: PetscScalar *g;
526: PetscReal *E;
527: PetscInt dim, d, Np, p;
529: PetscFunctionBeginUser;
530: PetscCall(TSGetDM(ts, &sw));
531: PetscCall(DMGetDimension(sw, &dim));
532: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
533: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
534: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
535: PetscCall(VecGetLocalSize(U, &Np));
536: PetscCall(VecGetArrayRead(U, &u));
537: PetscCall(VecGetArray(G, &g));
539: PetscCall(ComputeFieldAtParticles(snes, sw, E));
541: Np /= 2 * dim;
542: for (p = 0; p < Np; ++p) {
543: const PetscReal x0 = coords[p * dim + 0];
544: const PetscReal vy0 = vel[p * dim + 1];
545: const PetscReal omega = vy0 / x0;
547: for (d = 0; d < dim; ++d) {
548: g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d];
549: g[(p * 2 + 1) * dim + d] = E[p * dim + d] - PetscSqr(omega) * u[(p * 2 + 0) * dim + d];
550: }
551: }
552: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
553: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
554: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
555: PetscCall(VecRestoreArrayRead(U, &u));
556: PetscCall(VecRestoreArray(G, &g));
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: /* J_{ij} = dF_i/dx_j
561: J_p = ( 0 1)
562: (-w^2 0)
563: TODO Now there is another term with w^2 from the electric field. I think we will need to invert the operator.
564: Perhaps we can approximate the Jacobian using only the cellwise P-P gradient from Coulomb
565: */
566: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, PetscCtx ctx)
567: {
568: DM sw;
569: const PetscReal *coords, *vel;
570: PetscInt dim, d, Np, p, rStart;
572: PetscFunctionBeginUser;
573: PetscCall(TSGetDM(ts, &sw));
574: PetscCall(DMGetDimension(sw, &dim));
575: PetscCall(VecGetLocalSize(U, &Np));
576: PetscCall(MatGetOwnershipRange(J, &rStart, NULL));
577: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
578: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
579: Np /= 2 * dim;
580: for (p = 0; p < Np; ++p) {
581: const PetscReal x0 = coords[p * dim + 0];
582: const PetscReal vy0 = vel[p * dim + 1];
583: const PetscReal omega = vy0 / x0;
584: PetscScalar vals[4] = {0., 1., -PetscSqr(omega), 0.};
586: for (d = 0; d < dim; ++d) {
587: const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
588: PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES));
589: }
590: }
591: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
592: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
593: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
594: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
595: PetscFunctionReturn(PETSC_SUCCESS);
596: }
598: static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, PetscCtx ctx)
599: {
600: DM sw;
601: const PetscScalar *v;
602: PetscScalar *xres;
603: PetscInt Np, p, dim, d;
605: PetscFunctionBeginUser;
606: PetscCall(TSGetDM(ts, &sw));
607: PetscCall(DMGetDimension(sw, &dim));
608: PetscCall(VecGetLocalSize(Xres, &Np));
609: Np /= dim;
610: PetscCall(VecGetArrayRead(V, &v));
611: PetscCall(VecGetArray(Xres, &xres));
612: for (p = 0; p < Np; ++p) {
613: for (d = 0; d < dim; ++d) xres[p * dim + d] = v[p * dim + d];
614: }
615: PetscCall(VecRestoreArrayRead(V, &v));
616: PetscCall(VecRestoreArray(Xres, &xres));
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
620: static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, PetscCtx ctx)
621: {
622: DM sw;
623: SNES snes = ((AppCtx *)ctx)->snes;
624: const PetscScalar *x;
625: const PetscReal *coords, *vel;
626: PetscScalar *vres;
627: PetscReal *E;
628: PetscInt Np, p, dim, d;
630: PetscFunctionBeginUser;
631: PetscCall(TSGetDM(ts, &sw));
632: PetscCall(DMGetDimension(sw, &dim));
633: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
634: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
635: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
636: PetscCall(VecGetLocalSize(Vres, &Np));
637: PetscCall(VecGetArrayRead(X, &x));
638: PetscCall(VecGetArray(Vres, &vres));
639: PetscCheck(dim == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension must be 2");
641: PetscCall(ComputeFieldAtParticles(snes, sw, E));
643: Np /= dim;
644: for (p = 0; p < Np; ++p) {
645: const PetscReal x0 = coords[p * dim + 0];
646: const PetscReal vy0 = vel[p * dim + 1];
647: const PetscReal omega = vy0 / x0;
649: for (d = 0; d < dim; ++d) vres[p * dim + d] = E[p * dim + d] - PetscSqr(omega) * x[p * dim + d];
650: }
651: PetscCall(VecRestoreArrayRead(X, &x));
652: PetscCall(VecRestoreArray(Vres, &vres));
653: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
654: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
655: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: /* Discrete Gradients Formulation: S, F, gradF (G) */
660: PetscErrorCode RHSJacobianS(TS ts, PetscReal t, Vec U, Mat S, PetscCtx ctx)
661: {
662: PetscScalar vals[4] = {0., 1., -1., 0.};
663: DM sw;
664: PetscInt dim, d, Np, p, rStart;
666: PetscFunctionBeginUser;
667: PetscCall(TSGetDM(ts, &sw));
668: PetscCall(DMGetDimension(sw, &dim));
669: PetscCall(VecGetLocalSize(U, &Np));
670: PetscCall(MatGetOwnershipRange(S, &rStart, NULL));
671: Np /= 2 * dim;
672: for (p = 0; p < Np; ++p) {
673: for (d = 0; d < dim; ++d) {
674: const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
675: PetscCall(MatSetValues(S, 2, rows, 2, rows, vals, INSERT_VALUES));
676: }
677: }
678: PetscCall(MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY));
679: PetscCall(MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY));
680: PetscFunctionReturn(PETSC_SUCCESS);
681: }
683: PetscErrorCode RHSObjectiveF(TS ts, PetscReal t, Vec U, PetscScalar *F, PetscCtx ctx)
684: {
685: SNES snes = ((AppCtx *)ctx)->snes;
686: DM dm, sw;
687: const PetscScalar *u, *phi_vals;
688: PetscInt dim, Np, cStart, cEnd;
689: PetscReal *vel, *coords, m_p = 1., q_p = -1.;
690: Vec phi;
692: PetscFunctionBeginUser;
693: PetscCall(TSGetDM(ts, &sw));
694: PetscCall(DMGetDimension(sw, &dim));
695: PetscCall(SNESGetDM(snes, &dm));
696: PetscCall(VecGetArrayRead(U, &u));
697: PetscCall(VecGetLocalSize(U, &Np));
698: PetscCall(DMGetGlobalVector(dm, &phi));
699: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view_dg"));
700: PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
701: PetscInt phi_size;
702: PetscCall(VecGetSize(phi, &phi_size));
703: PetscCall(VecGetArrayRead(phi, &phi_vals));
704: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
705: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
707: PetscCall(DMSwarmSortGetAccess(sw));
708: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
709: Np /= 2 * dim;
710: for (PetscInt c = cStart; c < cEnd; ++c) {
711: PetscInt *points;
712: PetscInt Ncp;
713: PetscReal E = phi_vals[c];
715: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
716: for (PetscInt cp = 0; cp < Ncp; ++cp) {
717: const PetscInt p = points[cp];
718: const PetscReal x0 = coords[p * dim + 0];
719: const PetscReal vy0 = vel[p * dim + 1];
720: const PetscReal omega = vy0 / x0;
721: const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);
722: const PetscReal x2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 0) * dim], &u[(p * 2 + 0) * dim]);
723: E += 0.5 * q_p * m_p * (v2) + 0.5 * PetscSqr(omega) * (x2);
725: *F += E;
726: }
727: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
728: }
729: PetscCall(DMSwarmSortRestoreAccess(sw));
730: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
731: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
732: PetscCall(VecRestoreArrayRead(phi, &phi_vals));
733: PetscCall(DMRestoreGlobalVector(dm, &phi));
734: // PetscCall(DMSwarmRestoreField(sw, "potential", NULL, NULL, (void **)&pot));
735: PetscCall(VecRestoreArrayRead(U, &u));
736: PetscFunctionReturn(PETSC_SUCCESS);
737: }
739: /* dF/dx = q E dF/dv = v */
740: PetscErrorCode RHSFunctionG(TS ts, PetscReal t, Vec U, Vec G, PetscCtx ctx)
741: {
742: DM sw;
743: SNES snes = ((AppCtx *)ctx)->snes;
744: const PetscReal *coords, *vel;
745: const PetscScalar *u;
746: PetscScalar *g;
747: PetscReal *E, m_p = 1., q_p = -1.;
748: PetscInt dim, d, Np, p;
750: PetscFunctionBeginUser;
751: PetscCall(TSGetDM(ts, &sw));
752: PetscCall(DMGetDimension(sw, &dim));
753: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
754: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
755: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
756: PetscCall(DMSwarmGetLocalSize(sw, &Np));
757: PetscCall(VecGetArrayRead(U, &u));
758: PetscCall(VecGetArray(G, &g));
760: int COMPUTEFIELD;
761: PetscCall(PetscLogEventRegister("COMPFIELDATPART", TS_CLASSID, &COMPUTEFIELD));
762: PetscCall(PetscLogEventBegin(COMPUTEFIELD, 0, 0, 0, 0));
763: PetscCall(ComputeFieldAtParticles(snes, sw, E));
764: PetscCall(PetscLogEventEnd(COMPUTEFIELD, 0, 0, 0, 0));
765: for (p = 0; p < Np; ++p) {
766: const PetscReal x0 = coords[p * dim + 0];
767: const PetscReal vy0 = vel[p * dim + 1];
768: const PetscReal omega = vy0 / x0;
769: for (d = 0; d < dim; ++d) {
770: g[(p * 2 + 0) * dim + d] = -(q_p / m_p) * E[p * dim + d] + PetscSqr(omega) * u[(p * 2 + 0) * dim + d];
771: g[(p * 2 + 1) * dim + d] = u[(p * 2 + 1) * dim + d];
772: }
773: }
774: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
775: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
776: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
777: PetscCall(VecRestoreArrayRead(U, &u));
778: PetscCall(VecRestoreArray(G, &g));
779: PetscFunctionReturn(PETSC_SUCCESS);
780: }
782: static PetscErrorCode CreateSolution(TS ts)
783: {
784: DM sw;
785: Vec u;
786: PetscInt dim, Np;
788: PetscFunctionBegin;
789: PetscCall(TSGetDM(ts, &sw));
790: PetscCall(DMGetDimension(sw, &dim));
791: PetscCall(DMSwarmGetLocalSize(sw, &Np));
792: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
793: PetscCall(VecSetBlockSize(u, dim));
794: PetscCall(VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE));
795: PetscCall(VecSetUp(u));
796: PetscCall(TSSetSolution(ts, u));
797: PetscCall(VecDestroy(&u));
798: PetscFunctionReturn(PETSC_SUCCESS);
799: }
801: static PetscErrorCode SetProblem(TS ts)
802: {
803: AppCtx *user;
804: DM sw;
806: PetscFunctionBegin;
807: PetscCall(TSGetDM(ts, &sw));
808: PetscCall(DMGetApplicationContext(sw, &user));
809: // Define unified system for (X, V)
810: {
811: Mat J;
812: PetscInt dim, Np;
814: PetscCall(DMGetDimension(sw, &dim));
815: PetscCall(DMSwarmGetLocalSize(sw, &Np));
816: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
817: PetscCall(MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE));
818: PetscCall(MatSetBlockSize(J, 2 * dim));
819: PetscCall(MatSetFromOptions(J));
820: PetscCall(MatSetUp(J));
821: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, user));
822: PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, user));
823: PetscCall(MatDestroy(&J));
824: }
825: /* Define split system for X and V */
826: {
827: Vec u;
828: IS isx, isv, istmp;
829: const PetscInt *idx;
830: PetscInt dim, Np, rstart;
832: PetscCall(TSGetSolution(ts, &u));
833: PetscCall(DMGetDimension(sw, &dim));
834: PetscCall(DMSwarmGetLocalSize(sw, &Np));
835: PetscCall(VecGetOwnershipRange(u, &rstart, NULL));
836: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 0, 2, &istmp));
837: PetscCall(ISGetIndices(istmp, &idx));
838: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isx));
839: PetscCall(ISRestoreIndices(istmp, &idx));
840: PetscCall(ISDestroy(&istmp));
841: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 1, 2, &istmp));
842: PetscCall(ISGetIndices(istmp, &idx));
843: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isv));
844: PetscCall(ISRestoreIndices(istmp, &idx));
845: PetscCall(ISDestroy(&istmp));
846: PetscCall(TSRHSSplitSetIS(ts, "position", isx));
847: PetscCall(TSRHSSplitSetIS(ts, "momentum", isv));
848: PetscCall(ISDestroy(&isx));
849: PetscCall(ISDestroy(&isv));
850: PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, user));
851: PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, user));
852: }
853: // Define symplectic formulation U_t = S . G, where G = grad F
854: {
855: PetscCall(TSDiscGradSetFormulation(ts, RHSJacobianS, RHSObjectiveF, RHSFunctionG, user));
856: }
857: PetscFunctionReturn(PETSC_SUCCESS);
858: }
860: PetscErrorCode circleSingleX(PetscInt dim, PetscReal time, const PetscReal unused[], PetscInt p, PetscScalar x[], PetscCtx ctx)
861: {
862: x[0] = p + 1.;
863: x[1] = 0.;
864: return PETSC_SUCCESS;
865: }
867: PetscErrorCode circleSingleV(PetscInt dim, PetscReal time, const PetscReal unused[], PetscInt p, PetscScalar v[], PetscCtx ctx)
868: {
869: v[0] = 0.;
870: v[1] = PetscSqrtReal(1000. / (p + 1.));
871: return PETSC_SUCCESS;
872: }
874: /* Put 5 particles into each circle */
875: PetscErrorCode circleMultipleX(PetscInt dim, PetscReal time, const PetscReal unused[], PetscInt p, PetscScalar x[], PetscCtx ctx)
876: {
877: const PetscInt n = 5;
878: const PetscReal r0 = (p / n) + 1.;
879: const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
881: x[0] = r0 * PetscCosReal(th0);
882: x[1] = r0 * PetscSinReal(th0);
883: return PETSC_SUCCESS;
884: }
886: /* Put 5 particles into each circle */
887: PetscErrorCode circleMultipleV(PetscInt dim, PetscReal time, const PetscReal unused[], PetscInt p, PetscScalar v[], PetscCtx ctx)
888: {
889: const PetscInt n = 5;
890: const PetscReal r0 = (p / n) + 1.;
891: const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
892: const PetscReal omega = PetscSqrtReal(1000. / r0) / r0;
894: v[0] = -r0 * omega * PetscSinReal(th0);
895: v[1] = r0 * omega * PetscCosReal(th0);
896: return PETSC_SUCCESS;
897: }
899: /*
900: InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values.
902: Input Parameters:
903: + ts - The TS
904: - useInitial - Flag to also set the initial conditions to the current coordinates and velocities and setup the problem
906: Output Parameter:
907: . u - The initialized solution vector
909: Level: advanced
911: .seealso: InitializeSolve()
912: */
913: static PetscErrorCode InitializeSolveAndSwarm(TS ts, PetscBool useInitial)
914: {
915: DM sw;
916: Vec u, gc, gv, gc0, gv0;
917: IS isx, isv;
918: AppCtx *user;
920: PetscFunctionBeginUser;
921: PetscCall(TSGetDM(ts, &sw));
922: PetscCall(DMGetApplicationContext(sw, &user));
923: if (useInitial) {
924: PetscReal v0[1] = {1.};
926: PetscCall(DMSwarmInitializeCoordinates(sw));
927: PetscCall(DMSwarmInitializeVelocitiesFromOptions(sw, v0));
928: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
929: PetscCall(TSReset(ts));
930: PetscCall(CreateSolution(ts));
931: PetscCall(SetProblem(ts));
932: }
933: PetscCall(TSGetSolution(ts, &u));
934: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
935: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
936: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
937: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initCoordinates", &gc0));
938: if (useInitial) PetscCall(VecCopy(gc, gc0));
939: PetscCall(VecISCopy(u, isx, SCATTER_FORWARD, gc));
940: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
941: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initCoordinates", &gc0));
942: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
943: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initVelocity", &gv0));
944: if (useInitial) PetscCall(VecCopy(gv, gv0));
945: PetscCall(VecISCopy(u, isv, SCATTER_FORWARD, gv));
946: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
947: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initVelocity", &gv0));
948: PetscFunctionReturn(PETSC_SUCCESS);
949: }
951: static PetscErrorCode InitializeSolve(TS ts, Vec u)
952: {
953: PetscFunctionBegin;
954: PetscCall(TSSetSolution(ts, u));
955: PetscCall(InitializeSolveAndSwarm(ts, PETSC_TRUE));
956: PetscFunctionReturn(PETSC_SUCCESS);
957: }
959: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
960: {
961: MPI_Comm comm;
962: DM sw;
963: AppCtx *user;
964: const PetscScalar *u;
965: const PetscReal *coords, *vel;
966: PetscScalar *e;
967: PetscReal t;
968: PetscInt dim, Np, p;
970: PetscFunctionBeginUser;
971: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
972: PetscCall(TSGetDM(ts, &sw));
973: PetscCall(DMGetApplicationContext(sw, &user));
974: PetscCall(DMGetDimension(sw, &dim));
975: PetscCall(TSGetSolveTime(ts, &t));
976: PetscCall(VecGetArray(E, &e));
977: PetscCall(VecGetArrayRead(U, &u));
978: PetscCall(VecGetLocalSize(U, &Np));
979: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
980: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
981: Np /= 2 * dim;
982: for (p = 0; p < Np; ++p) {
983: /* TODO generalize initial conditions and project into plane instead of assuming x-y */
984: const PetscReal r0 = DMPlex_NormD_Internal(dim, &coords[p * dim]);
985: const PetscReal th0 = PetscAtan2Real(coords[p * dim + 1], coords[p * dim + 0]);
986: const PetscReal v0 = DMPlex_NormD_Internal(dim, &vel[p * dim]);
987: const PetscReal omega = v0 / r0;
988: const PetscReal ct = PetscCosReal(omega * t + th0);
989: const PetscReal st = PetscSinReal(omega * t + th0);
990: const PetscScalar *x = &u[(p * 2 + 0) * dim];
991: const PetscScalar *v = &u[(p * 2 + 1) * dim];
992: const PetscReal xe[3] = {r0 * ct, r0 * st, 0.0};
993: const PetscReal ve[3] = {-v0 * st, v0 * ct, 0.0};
994: PetscInt d;
996: for (d = 0; d < dim; ++d) {
997: e[(p * 2 + 0) * dim + d] = x[d] - xe[d];
998: e[(p * 2 + 1) * dim + d] = v[d] - ve[d];
999: }
1000: if (user->error) {
1001: const PetscReal en = 0.5 * DMPlex_DotRealD_Internal(dim, v, v);
1002: const PetscReal exen = 0.5 * PetscSqr(v0);
1003: PetscCall(PetscPrintf(comm, "t %.4g: p%" PetscInt_FMT " error [%.2f %.2f] 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)));
1004: }
1005: }
1006: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
1007: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
1008: PetscCall(VecRestoreArrayRead(U, &u));
1009: PetscCall(VecRestoreArray(E, &e));
1010: PetscFunctionReturn(PETSC_SUCCESS);
1011: }
1013: static PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal t, Vec U, PetscCtx ctx)
1014: {
1015: const PetscInt ostep = ((AppCtx *)ctx)->ostep;
1016: const EMType em = ((AppCtx *)ctx)->em;
1017: DM sw;
1018: const PetscScalar *u;
1019: PetscReal *coords, *E;
1020: PetscReal enKin = 0., enEM = 0.;
1021: PetscInt dim, d, Np, p, q;
1023: PetscFunctionBeginUser;
1024: if (step % ostep == 0) {
1025: PetscCall(TSGetDM(ts, &sw));
1026: PetscCall(DMGetDimension(sw, &dim));
1027: PetscCall(VecGetArrayRead(U, &u));
1028: PetscCall(VecGetLocalSize(U, &Np));
1029: Np /= 2 * dim;
1030: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1031: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
1032: if (!step) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Time Step Part Energy\n"));
1033: for (p = 0; p < Np; ++p) {
1034: const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);
1035: PetscReal *pcoord = &coords[p * dim];
1037: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " %5" PetscInt_FMT " %10.4lf\n", (double)t, step, p, (double)(0.5 * v2)));
1038: enKin += 0.5 * v2;
1039: if (em == EM_NONE) {
1040: continue;
1041: } else if (em == EM_COULOMB) {
1042: for (q = p + 1; q < Np; ++q) {
1043: PetscReal *qcoord = &coords[q * dim];
1044: PetscReal rpq[3], r;
1045: for (d = 0; d < dim; ++d) rpq[d] = pcoord[d] - qcoord[d];
1046: r = DMPlex_NormD_Internal(dim, rpq);
1047: enEM += 1. / r;
1048: }
1049: } else if (em == EM_PRIMAL || em == EM_MIXED) {
1050: for (d = 0; d < dim; ++d) enEM += E[p * dim + d];
1051: }
1052: }
1053: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " KE\t %10.4lf\n", (double)t, step, (double)enKin));
1054: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " PE\t %1.10g\n", (double)t, step, (double)enEM));
1055: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " E\t %10.4lf\n", (double)t, step, (double)(enKin + enEM)));
1056: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1057: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
1058: PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)ts), NULL));
1059: PetscCall(VecRestoreArrayRead(U, &u));
1060: }
1061: PetscFunctionReturn(PETSC_SUCCESS);
1062: }
1064: static PetscErrorCode SetUpMigrateParticles(TS ts, PetscInt n, PetscReal t, Vec x, PetscBool *flg, PetscCtx ctx)
1065: {
1066: DM sw;
1068: PetscFunctionBeginUser;
1069: *flg = PETSC_TRUE;
1070: PetscCall(TSGetDM(ts, &sw));
1071: PetscCall(DMViewFromOptions(sw, NULL, "-migrate_view_pre"));
1072: {
1073: Vec u, gc, gv;
1074: IS isx, isv;
1076: PetscCall(TSGetSolution(ts, &u));
1077: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
1078: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
1079: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
1080: PetscCall(VecISCopy(u, isx, SCATTER_REVERSE, gc));
1081: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
1082: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
1083: PetscCall(VecISCopy(u, isv, SCATTER_REVERSE, gv));
1084: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
1085: }
1086: PetscFunctionReturn(PETSC_SUCCESS);
1087: }
1089: static PetscErrorCode MigrateParticles(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], PetscCtx ctx)
1090: {
1091: DM sw;
1093: PetscFunctionBeginUser;
1094: PetscCall(TSGetDM(ts, &sw));
1095: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
1096: PetscCall(CreateSolution(ts));
1097: PetscCall(SetProblem(ts));
1098: PetscCall(InitializeSolveAndSwarm(ts, PETSC_FALSE));
1099: PetscFunctionReturn(PETSC_SUCCESS);
1100: }
1102: int main(int argc, char **argv)
1103: {
1104: DM dm, sw;
1105: TS ts;
1106: Vec u;
1107: AppCtx user;
1109: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1110: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
1111: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
1112: PetscCall(CreatePoisson(dm, &user));
1113: PetscCall(CreateSwarm(dm, &user, &sw));
1114: PetscCall(DMSetApplicationContext(sw, &user));
1116: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1117: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
1118: PetscCall(TSSetDM(ts, sw));
1119: PetscCall(TSSetMaxTime(ts, 0.1));
1120: PetscCall(TSSetTimeStep(ts, 0.00001));
1121: PetscCall(TSSetMaxSteps(ts, 100));
1122: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1123: PetscCall(TSMonitorSet(ts, EnergyMonitor, &user, NULL));
1124: PetscCall(TSSetFromOptions(ts));
1125: PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
1126: PetscCall(TSSetComputeExactError(ts, ComputeError));
1127: PetscCall(TSSetResize(ts, PETSC_FALSE, SetUpMigrateParticles, MigrateParticles, NULL));
1129: PetscCall(CreateSolution(ts));
1130: PetscCall(TSGetSolution(ts, &u));
1131: PetscCall(TSComputeInitialCondition(ts, u));
1132: PetscCall(TSSolve(ts, NULL));
1134: PetscCall(SNESDestroy(&user.snes));
1135: PetscCall(TSDestroy(&ts));
1136: PetscCall(DMDestroy(&sw));
1137: PetscCall(DMDestroy(&dm));
1138: PetscCall(PetscFinalize());
1139: return 0;
1140: }
1142: /*TEST
1144: build:
1145: requires: double !complex
1147: testset:
1148: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1149: 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 \
1150: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1151: -ts_type basicsymplectic\
1152: -dm_view -output_step 50 -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10
1153: test:
1154: suffix: none_bsi_2d_1
1155: args: -ts_basicsymplectic_type 1 -em_type none -error
1156: test:
1157: suffix: none_bsi_2d_2
1158: args: -ts_basicsymplectic_type 2 -em_type none -error
1159: test:
1160: suffix: none_bsi_2d_3
1161: args: -ts_basicsymplectic_type 3 -em_type none -error
1162: test:
1163: suffix: none_bsi_2d_4
1164: args: -ts_basicsymplectic_type 4 -em_type none -error
1165: test:
1166: suffix: coulomb_bsi_2d_1
1167: args: -ts_basicsymplectic_type 1
1168: test:
1169: suffix: coulomb_bsi_2d_2
1170: args: -ts_basicsymplectic_type 2
1171: test:
1172: suffix: coulomb_bsi_2d_3
1173: args: -ts_basicsymplectic_type 3
1174: test:
1175: suffix: coulomb_bsi_2d_4
1176: args: -ts_basicsymplectic_type 4
1178: testset:
1179: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1180: 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 \
1181: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1182: -ts_type basicsymplectic\
1183: -em_type primal -em_pc_type svd\
1184: -dm_view -output_step 50 -error -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10\
1185: -petscspace_degree 2 -petscfe_default_quadrature_order 3 -sigma 1.0e-8 -timeScale 2.0e-14
1186: test:
1187: suffix: poisson_bsi_2d_1
1188: args: -ts_basicsymplectic_type 1
1189: test:
1190: suffix: poisson_bsi_2d_2
1191: args: -ts_basicsymplectic_type 2
1192: test:
1193: suffix: poisson_bsi_2d_3
1194: args: -ts_basicsymplectic_type 3
1195: test:
1196: suffix: poisson_bsi_2d_4
1197: args: -ts_basicsymplectic_type 4
1199: testset:
1200: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1201: 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 \
1202: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1203: -ts_convergence_estimate -convest_num_refine 2 -em_type primal \
1204: -mat_type baij -em_ksp_error_if_not_converged -em_pc_type svd \
1205: -dm_view -output_step 50 -error -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10 \
1206: -sigma 1.0e-8 -timeScale 2.0e-14
1207: test:
1208: suffix: im_2d_0
1209: args: -ts_type theta -ts_theta_theta 0.5
1210: test:
1211: suffix: dg_2d_none
1212: args: -ts_type discgrad -ts_discgrad_type none -snes_type qn
1213: test:
1214: suffix: dg_2d_average
1215: args: -ts_type discgrad -ts_discgrad_type average -snes_type qn
1216: test:
1217: suffix: dg_2d_gonzalez
1218: args: -ts_type discgrad -ts_discgrad_type gonzalez -snes_fd -snes_type newtonls -snes_fd -pc_type lu
1220: testset:
1221: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1222: 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 \
1223: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV -dm_swarm_num_species 1\
1224: -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
1225: -em_snes_type ksponly -em_pc_type svd -em_type primal -petscspace_degree 1\
1226: -dm_view -output_step 50\
1227: -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_time_step 0.01 -ts_max_time 1.0 -ts_max_steps 10
1228: test:
1229: suffix: bsi_2d_mesh_1
1230: args: -ts_basicsymplectic_type 4
1231: test:
1232: suffix: bsi_2d_mesh_1_par_2
1233: nsize: 2
1234: args: -ts_basicsymplectic_type 4
1235: test:
1236: suffix: bsi_2d_mesh_1_par_3
1237: nsize: 3
1238: args: -ts_basicsymplectic_type 4
1239: test:
1240: suffix: bsi_2d_mesh_1_par_4
1241: nsize: 4
1242: args: -ts_basicsymplectic_type 4 -dm_swarm_num_particles 0,0,2,0
1244: testset:
1245: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1246: 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 \
1247: -dm_swarm_num_particles 10 -dm_swarm_coordinate_function circleMultipleX -dm_swarm_velocity_function circleMultipleV \
1248: -ts_convergence_estimate -convest_num_refine 2 \
1249: -em_pc_type lu\
1250: -dm_view -output_step 50 -error\
1251: -sigma 1.0e-8 -timeScale 2.0e-14 -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10
1252: test:
1253: suffix: bsi_2d_multiple_1
1254: args: -ts_type basicsymplectic -ts_basicsymplectic_type 1
1255: test:
1256: suffix: bsi_2d_multiple_2
1257: args: -ts_type basicsymplectic -ts_basicsymplectic_type 2
1258: test:
1259: suffix: bsi_2d_multiple_3
1260: args: -ts_type basicsymplectic -ts_basicsymplectic_type 3 -ts_time_step 0.001
1261: test:
1262: suffix: im_2d_multiple_0
1263: args: -ts_type theta -ts_theta_theta 0.5 \
1264: -mat_type baij -em_ksp_error_if_not_converged -em_pc_type lu
1266: testset:
1267: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1268: args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
1269: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1270: -em_pc_type fieldsplit -ksp_rtol 1e-10 -em_ksp_type preonly -em_type mixed -em_ksp_error_if_not_converged\
1271: -dm_view -output_step 50 -error -dm_refine 0\
1272: -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10
1273: test:
1274: suffix: bsi_4_rt_1
1275: args: -ts_type basicsymplectic -ts_basicsymplectic_type 4\
1276: -pc_fieldsplit_detect_saddle_point\
1277: -pc_type fieldsplit\
1278: -pc_fieldsplit_type schur\
1279: -pc_fieldsplit_schur_precondition full \
1280: -field_petscspace_degree 2\
1281: -field_petscfe_default_quadrature_order 1\
1282: -field_petscspace_type sum \
1283: -field_petscspace_variables 2 \
1284: -field_petscspace_components 2 \
1285: -field_petscspace_sum_spaces 2 \
1286: -field_petscspace_sum_concatenate true \
1287: -field_sumcomp_0_petscspace_variables 2 \
1288: -field_sumcomp_0_petscspace_type tensor \
1289: -field_sumcomp_0_petscspace_tensor_spaces 2 \
1290: -field_sumcomp_0_petscspace_tensor_uniform false \
1291: -field_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
1292: -field_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
1293: -field_sumcomp_1_petscspace_variables 2 \
1294: -field_sumcomp_1_petscspace_type tensor \
1295: -field_sumcomp_1_petscspace_tensor_spaces 2 \
1296: -field_sumcomp_1_petscspace_tensor_uniform false \
1297: -field_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
1298: -field_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
1299: -field_petscdualspace_form_degree -1 \
1300: -field_petscdualspace_order 1 \
1301: -field_petscdualspace_lagrange_trimmed true\
1302: -ksp_gmres_restart 500
1304: TEST*/