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_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>
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: PetscInt c;
104: for (c = 0; c < dim; ++c) g0[c * dim + c] += 1.0;
105: }
107: 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[])
108: {
109: for (PetscInt d = 0; d < dim; ++d) g2[d * dim + d] += 1.0;
110: }
112: 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[])
113: {
114: for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] += 1.0;
115: }
117: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
118: {
119: PetscFE feu, feq;
120: PetscDS ds;
121: PetscBool simplex;
122: PetscInt dim;
124: PetscFunctionBeginUser;
125: PetscCall(DMGetDimension(dm, &dim));
126: PetscCall(DMPlexIsSimplex(dm, &simplex));
127: if (user->em == EM_MIXED) {
128: DMLabel label;
130: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "field_", PETSC_DETERMINE, &feq));
131: PetscCall(PetscObjectSetName((PetscObject)feq, "field"));
132: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "potential_", PETSC_DETERMINE, &feu));
133: PetscCall(PetscObjectSetName((PetscObject)feu, "potential"));
134: PetscCall(PetscFECopyQuadrature(feq, feu));
135: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feq));
136: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)feu));
137: PetscCall(DMCreateDS(dm));
138: PetscCall(PetscFEDestroy(&feu));
139: PetscCall(PetscFEDestroy(&feq));
141: PetscCall(DMGetLabel(dm, "marker", &label));
142: PetscCall(DMGetDS(dm, &ds));
143: PetscCall(PetscDSSetResidual(ds, 0, f0_q, f1_q));
144: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_qq, NULL, NULL, NULL));
145: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_qu, NULL));
146: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_uq, NULL, NULL));
147: } else if (user->em == EM_PRIMAL) {
148: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &feu));
149: PetscCall(PetscObjectSetName((PetscObject)feu, "potential"));
150: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feu));
151: PetscCall(DMCreateDS(dm));
152: PetscCall(PetscFEDestroy(&feu));
153: PetscCall(DMGetDS(dm, &ds));
154: PetscCall(PetscDSSetResidual(ds, 0, NULL, laplacian_f1));
155: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, laplacian_g3));
156: }
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: static PetscErrorCode CreatePoisson(DM dm, AppCtx *user)
161: {
162: SNES snes;
163: Mat J;
164: MatNullSpace nullSpace;
166: PetscFunctionBeginUser;
167: PetscCall(CreateFEM(dm, user));
168: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
169: PetscCall(SNESSetOptionsPrefix(snes, "em_"));
170: PetscCall(SNESSetDM(snes, dm));
171: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, user));
172: PetscCall(SNESSetFromOptions(snes));
174: PetscCall(DMCreateMatrix(dm, &J));
175: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_TRUE, 0, NULL, &nullSpace));
176: PetscCall(MatSetNullSpace(J, nullSpace));
177: PetscCall(MatNullSpaceDestroy(&nullSpace));
178: PetscCall(SNESSetJacobian(snes, J, J, NULL, NULL));
179: PetscCall(MatDestroy(&J));
180: user->snes = snes;
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
185: {
186: PetscReal v0[1] = {1.};
187: PetscInt dim;
189: PetscFunctionBeginUser;
190: PetscCall(DMGetDimension(dm, &dim));
191: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
192: PetscCall(DMSetType(*sw, DMSWARM));
193: PetscCall(DMSetDimension(*sw, dim));
194: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
195: PetscCall(DMSwarmSetCellDM(*sw, dm));
196: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
197: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
198: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT));
199: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "initCoordinates", dim, PETSC_REAL));
200: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "initVelocity", dim, PETSC_REAL));
201: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "E_field", dim, PETSC_REAL));
202: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
203: PetscCall(DMSwarmComputeLocalSizeFromOptions(*sw));
204: PetscCall(DMSwarmInitializeCoordinates(*sw));
205: PetscCall(DMSwarmInitializeVelocitiesFromOptions(*sw, v0));
206: PetscCall(DMSetFromOptions(*sw));
207: PetscCall(DMSetApplicationContext(*sw, user));
208: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
209: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
210: {
211: Vec gc, gc0, gv, gv0;
213: PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc));
214: PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "initCoordinates", &gc0));
215: PetscCall(VecCopy(gc, gc0));
216: PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc));
217: PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "initCoordinates", &gc0));
218: PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "velocity", &gv));
219: PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "initVelocity", &gv0));
220: PetscCall(VecCopy(gv, gv0));
221: PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "velocity", &gv));
222: PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "initVelocity", &gv0));
223: }
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: static PetscErrorCode ComputeFieldAtParticles_Coulomb(SNES snes, DM sw, PetscReal E[])
228: {
229: PetscReal *coords;
230: PetscInt dim, d, Np, p, q;
231: PetscMPIInt size;
233: PetscFunctionBegin;
234: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)snes), &size));
235: PetscCheck(size == 1, PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Coulomb code only works in serial");
236: PetscCall(DMGetDimension(sw, &dim));
237: PetscCall(DMSwarmGetLocalSize(sw, &Np));
239: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
240: for (p = 0; p < Np; ++p) {
241: PetscReal *pcoord = &coords[p * dim];
242: PetscReal *pE = &E[p * dim];
243: /* Calculate field at particle p due to particle q */
244: for (q = 0; q < Np; ++q) {
245: PetscReal *qcoord = &coords[q * dim];
246: PetscReal rpq[3], r;
248: if (p == q) continue;
249: for (d = 0; d < dim; ++d) rpq[d] = pcoord[d] - qcoord[d];
250: r = DMPlex_NormD_Internal(dim, rpq);
251: for (d = 0; d < dim; ++d) pE[d] += rpq[d] / PetscPowRealInt(r, 3);
252: }
253: }
254: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: static PetscErrorCode ComputeFieldAtParticles_Primal(SNES snes, DM sw, PetscReal E[])
259: {
260: DM dm;
261: PetscDS ds;
262: PetscFE fe;
263: Mat M_p;
264: Vec phi, locPhi, rho, f;
265: PetscReal *coords, chargeTol = 1e-13;
266: PetscInt dim, d, cStart, cEnd, c, Np;
267: char oldField[PETSC_MAX_PATH_LEN];
268: const char *tmp;
270: PetscFunctionBegin;
271: PetscCall(DMGetDimension(sw, &dim));
272: PetscCall(DMSwarmGetLocalSize(sw, &Np));
273: PetscCall(SNESGetDM(snes, &dm));
275: PetscCall(DMSwarmVectorGetField(sw, &tmp));
276: PetscCall(PetscStrncpy(oldField, tmp, PETSC_MAX_PATH_LEN));
277: PetscCall(DMSwarmVectorDefineField(sw, "w_q"));
278: PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
279: PetscCall(DMSwarmVectorDefineField(sw, oldField));
281: /* Create the charges rho */
282: PetscCall(DMGetGlobalVector(dm, &rho));
283: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
284: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
285: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
286: PetscCall(MatMultTranspose(M_p, f, rho));
287: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
288: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
289: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
290: PetscCall(MatDestroy(&M_p));
291: {
292: PetscScalar sum;
293: PetscInt n;
294: PetscReal phi_0 = 1.; /* (sigma*sigma*sigma)*(timeScale*timeScale)/(m_e*q_e*epsi_0)*/
296: /* Remove constant from rho */
297: PetscCall(VecGetSize(rho, &n));
298: PetscCall(VecSum(rho, &sum));
299: PetscCall(VecShift(rho, -sum / n));
300: PetscCall(VecSum(rho, &sum));
301: PetscCheck(PetscAbsScalar(sum) < chargeTol, PetscObjectComm((PetscObject)sw), PETSC_ERR_PLIB, "Charge should have no DC component: %g", (double)PetscRealPart(sum));
302: /* Nondimensionalize rho */
303: PetscCall(VecScale(rho, phi_0));
304: }
305: PetscCall(VecViewFromOptions(rho, NULL, "-poisson_rho_view"));
307: PetscCall(DMGetGlobalVector(dm, &phi));
308: PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
309: PetscCall(VecSet(phi, 0.0));
310: PetscCall(SNESSolve(snes, rho, phi));
311: PetscCall(DMRestoreGlobalVector(dm, &rho));
312: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
314: PetscCall(DMGetLocalVector(dm, &locPhi));
315: PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
316: PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
317: PetscCall(DMRestoreGlobalVector(dm, &phi));
319: PetscCall(DMGetDS(dm, &ds));
320: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
321: PetscCall(DMSwarmSortGetAccess(sw));
322: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
323: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
324: for (c = cStart; c < cEnd; ++c) {
325: PetscTabulation tab;
326: PetscScalar *clPhi = NULL;
327: PetscReal *pcoord, *refcoord;
328: PetscReal v[3], J[9], invJ[9], detJ;
329: PetscInt *points;
330: PetscInt Ncp, cp;
332: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
333: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
334: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
335: for (cp = 0; cp < Ncp; ++cp)
336: for (d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
337: PetscCall(DMPlexCoordinatesToReference(dm, c, Ncp, pcoord, refcoord));
338: PetscCall(PetscFECreateTabulation(fe, 1, Ncp, refcoord, 1, &tab));
339: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v, J, invJ, &detJ));
340: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
341: for (cp = 0; cp < Ncp; ++cp) {
342: const PetscReal *basisDer = tab->T[1];
343: const PetscInt p = points[cp];
345: for (d = 0; d < dim; ++d) E[p * dim + d] = 0.;
346: PetscCall(PetscFEFreeInterpolateGradient_Static(fe, basisDer, clPhi, dim, invJ, NULL, cp, &E[p * dim]));
347: for (d = 0; d < dim; ++d) E[p * dim + d] *= -1.0;
348: }
349: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
350: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
351: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
352: PetscCall(PetscTabulationDestroy(&tab));
353: PetscCall(PetscFree(points));
354: }
355: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
356: PetscCall(DMSwarmSortRestoreAccess(sw));
357: PetscCall(DMRestoreLocalVector(dm, &locPhi));
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: static PetscErrorCode ComputeFieldAtParticles_Mixed(SNES snes, DM sw, PetscReal E[])
362: {
363: DM dm, potential_dm;
364: IS potential_IS;
365: PetscDS ds;
366: PetscFE fe;
367: PetscFEGeom feGeometry;
368: Mat M_p;
369: Vec phi, locPhi, rho, f, temp_rho;
370: PetscQuadrature q;
371: PetscReal *coords, chargeTol = 1e-13;
372: PetscInt dim, d, cStart, cEnd, c, Np, fields = 1;
373: char oldField[PETSC_MAX_PATH_LEN];
374: const char *tmp;
376: PetscFunctionBegin;
377: PetscCall(DMGetDimension(sw, &dim));
378: PetscCall(DMSwarmGetLocalSize(sw, &Np));
380: /* Create the charges rho */
381: PetscCall(SNESGetDM(snes, &dm));
382: PetscCall(DMGetGlobalVector(dm, &rho));
383: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
384: PetscCall(DMCreateSubDM(dm, 1, &fields, &potential_IS, &potential_dm));
386: PetscCall(DMSwarmVectorGetField(sw, &tmp));
387: PetscCall(PetscStrncpy(oldField, tmp, PETSC_MAX_PATH_LEN));
388: PetscCall(DMSwarmVectorDefineField(sw, "w_q"));
389: PetscCall(DMCreateMassMatrix(sw, potential_dm, &M_p));
390: PetscCall(DMSwarmVectorDefineField(sw, oldField));
392: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
393: PetscCall(DMGetGlobalVector(potential_dm, &temp_rho));
394: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
395: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
396: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
397: PetscCall(MatMultTranspose(M_p, f, temp_rho));
398: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
399: PetscCall(MatDestroy(&M_p));
400: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
401: PetscCall(VecViewFromOptions(rho, NULL, "-poisson_rho_view"));
402: PetscCall(VecISCopy(rho, potential_IS, SCATTER_FORWARD, temp_rho));
403: PetscCall(DMRestoreGlobalVector(potential_dm, &temp_rho));
404: PetscCall(DMDestroy(&potential_dm));
405: PetscCall(ISDestroy(&potential_IS));
406: {
407: PetscScalar sum;
408: PetscInt n;
409: PetscReal phi_0 = 1.; /*(sigma*sigma*sigma)*(timeScale*timeScale)/(m_e*q_e*epsi_0);*/
411: /*Remove constant from rho*/
412: PetscCall(VecGetSize(rho, &n));
413: PetscCall(VecSum(rho, &sum));
414: PetscCall(VecShift(rho, -sum / n));
415: PetscCall(VecSum(rho, &sum));
416: PetscCheck(PetscAbsScalar(sum) < chargeTol, PetscObjectComm((PetscObject)sw), PETSC_ERR_PLIB, "Charge should have no DC component: %g", (double)PetscRealPart(sum));
417: /* Nondimensionalize rho */
418: PetscCall(VecScale(rho, phi_0));
419: }
420: PetscCall(DMGetGlobalVector(dm, &phi));
421: PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
422: PetscCall(VecSet(phi, 0.0));
423: PetscCall(SNESSolve(snes, NULL, phi));
424: PetscCall(DMRestoreGlobalVector(dm, &rho));
425: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
427: PetscCall(DMGetLocalVector(dm, &locPhi));
428: PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
429: PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
430: PetscCall(DMRestoreGlobalVector(dm, &phi));
432: PetscCall(DMGetDS(dm, &ds));
433: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
434: PetscCall(DMSwarmSortGetAccess(sw));
435: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
436: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
437: for (c = cStart; c < cEnd; ++c) {
438: PetscTabulation tab;
439: PetscScalar *clPhi = NULL;
440: PetscReal *pcoord, *refcoord;
441: PetscReal v[3], J[9], invJ[9], detJ;
442: PetscInt *points;
443: PetscInt Ncp, cp;
445: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
446: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
447: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
448: for (cp = 0; cp < Ncp; ++cp)
449: for (d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
450: PetscCall(DMPlexCoordinatesToReference(dm, c, Ncp, pcoord, refcoord));
451: PetscCall(PetscFECreateTabulation(fe, 1, Ncp, refcoord, 1, &tab));
452: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v, J, invJ, &detJ));
453: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
454: for (cp = 0; cp < Ncp; ++cp) {
455: const PetscInt p = points[cp];
457: for (d = 0; d < dim; ++d) E[p * dim + d] = 0.;
458: PetscCall(PetscFEGetQuadrature(fe, &q));
459: PetscCall(PetscFECreateCellGeometry(fe, q, &feGeometry));
460: PetscCall(PetscFEInterpolateAtPoints_Static(fe, tab, clPhi, &feGeometry, cp, &E[p * dim]));
461: PetscCall(PetscFEDestroyCellGeometry(fe, &feGeometry));
462: }
463: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
464: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
465: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
466: PetscCall(PetscTabulationDestroy(&tab));
467: PetscCall(PetscFree(points));
468: }
469: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
470: PetscCall(DMSwarmSortRestoreAccess(sw));
471: PetscCall(DMRestoreLocalVector(dm, &locPhi));
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: static PetscErrorCode ComputeFieldAtParticles(SNES snes, DM sw, PetscReal E[])
476: {
477: AppCtx *ctx;
478: PetscInt dim, Np;
480: PetscFunctionBegin;
483: PetscAssertPointer(E, 3);
484: PetscCall(DMGetDimension(sw, &dim));
485: PetscCall(DMSwarmGetLocalSize(sw, &Np));
486: PetscCall(DMGetApplicationContext(sw, &ctx));
487: PetscCall(PetscArrayzero(E, Np * dim));
489: switch (ctx->em) {
490: case EM_PRIMAL:
491: PetscCall(ComputeFieldAtParticles_Primal(snes, sw, E));
492: break;
493: case EM_COULOMB:
494: PetscCall(ComputeFieldAtParticles_Coulomb(snes, sw, E));
495: break;
496: case EM_MIXED:
497: PetscCall(ComputeFieldAtParticles_Mixed(snes, sw, E));
498: break;
499: case EM_NONE:
500: break;
501: default:
502: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No solver for electrostatic model %s", EMTypes[ctx->em]);
503: }
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx)
508: {
509: DM sw;
510: SNES snes = ((AppCtx *)ctx)->snes;
511: const PetscReal *coords, *vel;
512: const PetscScalar *u;
513: PetscScalar *g;
514: PetscReal *E;
515: PetscInt dim, d, Np, p;
517: PetscFunctionBeginUser;
518: PetscCall(TSGetDM(ts, &sw));
519: PetscCall(DMGetDimension(sw, &dim));
520: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
521: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
522: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
523: PetscCall(VecGetLocalSize(U, &Np));
524: PetscCall(VecGetArrayRead(U, &u));
525: PetscCall(VecGetArray(G, &g));
527: PetscCall(ComputeFieldAtParticles(snes, sw, E));
529: Np /= 2 * dim;
530: for (p = 0; p < Np; ++p) {
531: const PetscReal x0 = coords[p * dim + 0];
532: const PetscReal vy0 = vel[p * dim + 1];
533: const PetscReal omega = vy0 / x0;
535: for (d = 0; d < dim; ++d) {
536: g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d];
537: g[(p * 2 + 1) * dim + d] = E[p * dim + d] - PetscSqr(omega) * u[(p * 2 + 0) * dim + d];
538: }
539: }
540: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
541: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
542: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
543: PetscCall(VecRestoreArrayRead(U, &u));
544: PetscCall(VecRestoreArray(G, &g));
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: /* J_{ij} = dF_i/dx_j
549: J_p = ( 0 1)
550: (-w^2 0)
551: TODO Now there is another term with w^2 from the electric field. I think we will need to invert the operator.
552: Perhaps we can approximate the Jacobian using only the cellwise P-P gradient from Coulomb
553: */
554: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, void *ctx)
555: {
556: DM sw;
557: const PetscReal *coords, *vel;
558: PetscInt dim, d, Np, p, rStart;
560: PetscFunctionBeginUser;
561: PetscCall(TSGetDM(ts, &sw));
562: PetscCall(DMGetDimension(sw, &dim));
563: PetscCall(VecGetLocalSize(U, &Np));
564: PetscCall(MatGetOwnershipRange(J, &rStart, NULL));
565: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
566: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
567: Np /= 2 * dim;
568: for (p = 0; p < Np; ++p) {
569: const PetscReal x0 = coords[p * dim + 0];
570: const PetscReal vy0 = vel[p * dim + 1];
571: const PetscReal omega = vy0 / x0;
572: PetscScalar vals[4] = {0., 1., -PetscSqr(omega), 0.};
574: for (d = 0; d < dim; ++d) {
575: const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
576: PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES));
577: }
578: }
579: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
580: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
581: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
582: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
586: static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
587: {
588: DM sw;
589: const PetscScalar *v;
590: PetscScalar *xres;
591: PetscInt Np, p, dim, d;
593: PetscFunctionBeginUser;
594: PetscCall(TSGetDM(ts, &sw));
595: PetscCall(DMGetDimension(sw, &dim));
596: PetscCall(VecGetLocalSize(Xres, &Np));
597: Np /= dim;
598: PetscCall(VecGetArrayRead(V, &v));
599: PetscCall(VecGetArray(Xres, &xres));
600: for (p = 0; p < Np; ++p) {
601: for (d = 0; d < dim; ++d) xres[p * dim + d] = v[p * dim + d];
602: }
603: PetscCall(VecRestoreArrayRead(V, &v));
604: PetscCall(VecRestoreArray(Xres, &xres));
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
609: {
610: DM sw;
611: SNES snes = ((AppCtx *)ctx)->snes;
612: const PetscScalar *x;
613: const PetscReal *coords, *vel;
614: PetscScalar *vres;
615: PetscReal *E;
616: PetscInt Np, p, dim, d;
618: PetscFunctionBeginUser;
619: PetscCall(TSGetDM(ts, &sw));
620: PetscCall(DMGetDimension(sw, &dim));
621: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
622: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
623: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
624: PetscCall(VecGetLocalSize(Vres, &Np));
625: PetscCall(VecGetArrayRead(X, &x));
626: PetscCall(VecGetArray(Vres, &vres));
627: PetscCheck(dim == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension must be 2");
629: PetscCall(ComputeFieldAtParticles(snes, sw, E));
631: Np /= dim;
632: for (p = 0; p < Np; ++p) {
633: const PetscReal x0 = coords[p * dim + 0];
634: const PetscReal vy0 = vel[p * dim + 1];
635: const PetscReal omega = vy0 / x0;
637: for (d = 0; d < dim; ++d) vres[p * dim + d] = E[p * dim + d] - PetscSqr(omega) * x[p * dim + d];
638: }
639: PetscCall(VecRestoreArrayRead(X, &x));
640: PetscCall(VecRestoreArray(Vres, &vres));
641: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
642: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
643: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
644: PetscFunctionReturn(PETSC_SUCCESS);
645: }
647: static PetscErrorCode CreateSolution(TS ts)
648: {
649: DM sw;
650: Vec u;
651: PetscInt dim, Np;
653: PetscFunctionBegin;
654: PetscCall(TSGetDM(ts, &sw));
655: PetscCall(DMGetDimension(sw, &dim));
656: PetscCall(DMSwarmGetLocalSize(sw, &Np));
657: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
658: PetscCall(VecSetBlockSize(u, dim));
659: PetscCall(VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE));
660: PetscCall(VecSetUp(u));
661: PetscCall(TSSetSolution(ts, u));
662: PetscCall(VecDestroy(&u));
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: static PetscErrorCode SetProblem(TS ts)
667: {
668: AppCtx *user;
669: DM sw;
671: PetscFunctionBegin;
672: PetscCall(TSGetDM(ts, &sw));
673: PetscCall(DMGetApplicationContext(sw, (void **)&user));
674: // Define unified system for (X, V)
675: {
676: Mat J;
677: PetscInt dim, Np;
679: PetscCall(DMGetDimension(sw, &dim));
680: PetscCall(DMSwarmGetLocalSize(sw, &Np));
681: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
682: PetscCall(MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE));
683: PetscCall(MatSetBlockSize(J, 2 * dim));
684: PetscCall(MatSetFromOptions(J));
685: PetscCall(MatSetUp(J));
686: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, user));
687: PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, user));
688: PetscCall(MatDestroy(&J));
689: }
690: /* Define split system for X and V */
691: {
692: Vec u;
693: IS isx, isv, istmp;
694: const PetscInt *idx;
695: PetscInt dim, Np, rstart;
697: PetscCall(TSGetSolution(ts, &u));
698: PetscCall(DMGetDimension(sw, &dim));
699: PetscCall(DMSwarmGetLocalSize(sw, &Np));
700: PetscCall(VecGetOwnershipRange(u, &rstart, NULL));
701: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 0, 2, &istmp));
702: PetscCall(ISGetIndices(istmp, &idx));
703: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isx));
704: PetscCall(ISRestoreIndices(istmp, &idx));
705: PetscCall(ISDestroy(&istmp));
706: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 1, 2, &istmp));
707: PetscCall(ISGetIndices(istmp, &idx));
708: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isv));
709: PetscCall(ISRestoreIndices(istmp, &idx));
710: PetscCall(ISDestroy(&istmp));
711: PetscCall(TSRHSSplitSetIS(ts, "position", isx));
712: PetscCall(TSRHSSplitSetIS(ts, "momentum", isv));
713: PetscCall(ISDestroy(&isx));
714: PetscCall(ISDestroy(&isv));
715: PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, user));
716: PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, user));
717: }
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: PetscErrorCode circleSingleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
722: {
723: x[0] = p + 1.;
724: x[1] = 0.;
725: return PETSC_SUCCESS;
726: }
728: PetscErrorCode circleSingleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *ctx)
729: {
730: v[0] = 0.;
731: v[1] = PetscSqrtReal(1000. / (p + 1.));
732: return PETSC_SUCCESS;
733: }
735: /* Put 5 particles into each circle */
736: PetscErrorCode circleMultipleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
737: {
738: const PetscInt n = 5;
739: const PetscReal r0 = (p / n) + 1.;
740: const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
742: x[0] = r0 * PetscCosReal(th0);
743: x[1] = r0 * PetscSinReal(th0);
744: return PETSC_SUCCESS;
745: }
747: /* Put 5 particles into each circle */
748: PetscErrorCode circleMultipleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *ctx)
749: {
750: const PetscInt n = 5;
751: const PetscReal r0 = (p / n) + 1.;
752: const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
753: const PetscReal omega = PetscSqrtReal(1000. / r0) / r0;
755: v[0] = -r0 * omega * PetscSinReal(th0);
756: v[1] = r0 * omega * PetscCosReal(th0);
757: return PETSC_SUCCESS;
758: }
760: /*
761: InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values.
763: Input Parameters:
764: + ts - The TS
765: - useInitial - Flag to also set the initial conditions to the current coordinates and velocities and setup the problem
767: Output Parameter:
768: . u - The initialized solution vector
770: Level: advanced
772: .seealso: InitializeSolve()
773: */
774: static PetscErrorCode InitializeSolveAndSwarm(TS ts, PetscBool useInitial)
775: {
776: DM sw;
777: Vec u, gc, gv, gc0, gv0;
778: IS isx, isv;
779: AppCtx *user;
781: PetscFunctionBeginUser;
782: PetscCall(TSGetDM(ts, &sw));
783: PetscCall(DMGetApplicationContext(sw, &user));
784: if (useInitial) {
785: PetscReal v0[1] = {1.};
787: PetscCall(DMSwarmInitializeCoordinates(sw));
788: PetscCall(DMSwarmInitializeVelocitiesFromOptions(sw, v0));
789: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
790: PetscCall(TSReset(ts));
791: PetscCall(CreateSolution(ts));
792: PetscCall(SetProblem(ts));
793: }
794: PetscCall(TSGetSolution(ts, &u));
795: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
796: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
797: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
798: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initCoordinates", &gc0));
799: if (useInitial) PetscCall(VecCopy(gc, gc0));
800: PetscCall(VecISCopy(u, isx, SCATTER_FORWARD, gc));
801: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
802: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initCoordinates", &gc0));
803: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
804: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initVelocity", &gv0));
805: if (useInitial) PetscCall(VecCopy(gv, gv0));
806: PetscCall(VecISCopy(u, isv, SCATTER_FORWARD, gv));
807: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
808: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initVelocity", &gv0));
809: PetscFunctionReturn(PETSC_SUCCESS);
810: }
812: static PetscErrorCode InitializeSolve(TS ts, Vec u)
813: {
814: PetscFunctionBegin;
815: PetscCall(TSSetSolution(ts, u));
816: PetscCall(InitializeSolveAndSwarm(ts, PETSC_TRUE));
817: PetscFunctionReturn(PETSC_SUCCESS);
818: }
820: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
821: {
822: MPI_Comm comm;
823: DM sw;
824: AppCtx *user;
825: const PetscScalar *u;
826: const PetscReal *coords, *vel;
827: PetscScalar *e;
828: PetscReal t;
829: PetscInt dim, Np, p;
831: PetscFunctionBeginUser;
832: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
833: PetscCall(TSGetDM(ts, &sw));
834: PetscCall(DMGetApplicationContext(sw, &user));
835: PetscCall(DMGetDimension(sw, &dim));
836: PetscCall(TSGetSolveTime(ts, &t));
837: PetscCall(VecGetArray(E, &e));
838: PetscCall(VecGetArrayRead(U, &u));
839: PetscCall(VecGetLocalSize(U, &Np));
840: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
841: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
842: Np /= 2 * dim;
843: for (p = 0; p < Np; ++p) {
844: /* TODO generalize initial conditions and project into plane instead of assuming x-y */
845: const PetscReal r0 = DMPlex_NormD_Internal(dim, &coords[p * dim]);
846: const PetscReal th0 = PetscAtan2Real(coords[p * dim + 1], coords[p * dim + 0]);
847: const PetscReal v0 = DMPlex_NormD_Internal(dim, &vel[p * dim]);
848: const PetscReal omega = v0 / r0;
849: const PetscReal ct = PetscCosReal(omega * t + th0);
850: const PetscReal st = PetscSinReal(omega * t + th0);
851: const PetscScalar *x = &u[(p * 2 + 0) * dim];
852: const PetscScalar *v = &u[(p * 2 + 1) * dim];
853: const PetscReal xe[3] = {r0 * ct, r0 * st, 0.0};
854: const PetscReal ve[3] = {-v0 * st, v0 * ct, 0.0};
855: PetscInt d;
857: for (d = 0; d < dim; ++d) {
858: e[(p * 2 + 0) * dim + d] = x[d] - xe[d];
859: e[(p * 2 + 1) * dim + d] = v[d] - ve[d];
860: }
861: if (user->error) {
862: const PetscReal en = 0.5 * DMPlex_DotRealD_Internal(dim, v, v);
863: const PetscReal exen = 0.5 * PetscSqr(v0);
864: 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)));
865: }
866: }
867: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
868: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
869: PetscCall(VecRestoreArrayRead(U, &u));
870: PetscCall(VecRestoreArray(E, &e));
871: PetscFunctionReturn(PETSC_SUCCESS);
872: }
874: static PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
875: {
876: const PetscInt ostep = ((AppCtx *)ctx)->ostep;
877: const EMType em = ((AppCtx *)ctx)->em;
878: DM sw;
879: const PetscScalar *u;
880: PetscReal *coords, *E;
881: PetscReal enKin = 0., enEM = 0.;
882: PetscInt dim, d, Np, p, q;
884: PetscFunctionBeginUser;
885: if (step % ostep == 0) {
886: PetscCall(TSGetDM(ts, &sw));
887: PetscCall(DMGetDimension(sw, &dim));
888: PetscCall(VecGetArrayRead(U, &u));
889: PetscCall(VecGetLocalSize(U, &Np));
890: Np /= 2 * dim;
891: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
892: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
893: if (!step) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Time Step Part Energy\n"));
894: for (p = 0; p < Np; ++p) {
895: const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);
896: PetscReal *pcoord = &coords[p * dim];
898: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " %5" PetscInt_FMT " %10.4lf\n", (double)t, step, p, (double)(0.5 * v2)));
899: enKin += 0.5 * v2;
900: if (em == EM_NONE) {
901: continue;
902: } else if (em == EM_COULOMB) {
903: for (q = p + 1; q < Np; ++q) {
904: PetscReal *qcoord = &coords[q * dim];
905: PetscReal rpq[3], r;
906: for (d = 0; d < dim; ++d) rpq[d] = pcoord[d] - qcoord[d];
907: r = DMPlex_NormD_Internal(dim, rpq);
908: enEM += 1. / r;
909: }
910: } else if (em == EM_PRIMAL || em == EM_MIXED) {
911: for (d = 0; d < dim; ++d) enEM += E[p * dim + d];
912: }
913: }
914: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " KE\t %10.4lf\n", (double)t, step, (double)enKin));
915: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " PE\t %1.10g\n", (double)t, step, (double)enEM));
916: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " E\t %10.4lf\n", (double)t, step, (double)(enKin + enEM)));
917: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
918: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
919: PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)ts), NULL));
920: PetscCall(VecRestoreArrayRead(U, &u));
921: }
922: PetscFunctionReturn(PETSC_SUCCESS);
923: }
925: static PetscErrorCode SetUpMigrateParticles(TS ts, PetscInt n, PetscReal t, Vec x, PetscBool *flg, void *ctx)
926: {
927: DM sw;
929: PetscFunctionBeginUser;
930: *flg = PETSC_TRUE;
931: PetscCall(TSGetDM(ts, &sw));
932: PetscCall(DMViewFromOptions(sw, NULL, "-migrate_view_pre"));
933: {
934: Vec u, gc, gv;
935: IS isx, isv;
937: PetscCall(TSGetSolution(ts, &u));
938: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
939: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
940: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
941: PetscCall(VecISCopy(u, isx, SCATTER_REVERSE, gc));
942: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
943: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
944: PetscCall(VecISCopy(u, isv, SCATTER_REVERSE, gv));
945: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
946: }
947: PetscFunctionReturn(PETSC_SUCCESS);
948: }
950: static PetscErrorCode MigrateParticles(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx)
951: {
952: DM sw;
954: PetscFunctionBeginUser;
955: PetscCall(TSGetDM(ts, &sw));
956: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
957: PetscCall(CreateSolution(ts));
958: PetscCall(SetProblem(ts));
959: PetscCall(InitializeSolveAndSwarm(ts, PETSC_FALSE));
960: PetscFunctionReturn(PETSC_SUCCESS);
961: }
963: int main(int argc, char **argv)
964: {
965: DM dm, sw;
966: TS ts;
967: Vec u;
968: AppCtx user;
970: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
971: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
972: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
973: PetscCall(CreatePoisson(dm, &user));
974: PetscCall(CreateSwarm(dm, &user, &sw));
975: PetscCall(DMSetApplicationContext(sw, &user));
977: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
978: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
979: PetscCall(TSSetDM(ts, sw));
980: PetscCall(TSSetMaxTime(ts, 0.1));
981: PetscCall(TSSetTimeStep(ts, 0.00001));
982: PetscCall(TSSetMaxSteps(ts, 100));
983: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
984: PetscCall(TSMonitorSet(ts, EnergyMonitor, &user, NULL));
985: PetscCall(TSSetFromOptions(ts));
986: PetscCall(DMSwarmVectorDefineField(sw, "velocity"));
987: PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
988: PetscCall(TSSetComputeExactError(ts, ComputeError));
989: PetscCall(TSSetResize(ts, PETSC_FALSE, SetUpMigrateParticles, MigrateParticles, NULL));
991: PetscCall(CreateSolution(ts));
992: PetscCall(TSGetSolution(ts, &u));
993: PetscCall(TSComputeInitialCondition(ts, u));
994: PetscCall(TSSolve(ts, NULL));
996: PetscCall(SNESDestroy(&user.snes));
997: PetscCall(TSDestroy(&ts));
998: PetscCall(DMDestroy(&sw));
999: PetscCall(DMDestroy(&dm));
1000: PetscCall(PetscFinalize());
1001: return 0;
1002: }
1004: /*TEST
1006: build:
1007: requires: double !complex
1009: testset:
1010: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1011: 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 \
1012: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1013: -ts_type basicsymplectic\
1014: -dm_view -output_step 50 -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10
1015: test:
1016: suffix: none_bsi_2d_1
1017: args: -ts_basicsymplectic_type 1 -em_type none -error
1018: test:
1019: suffix: none_bsi_2d_2
1020: args: -ts_basicsymplectic_type 2 -em_type none -error
1021: test:
1022: suffix: none_bsi_2d_3
1023: args: -ts_basicsymplectic_type 3 -em_type none -error
1024: test:
1025: suffix: none_bsi_2d_4
1026: args: -ts_basicsymplectic_type 4 -em_type none -error
1027: test:
1028: suffix: coulomb_bsi_2d_1
1029: args: -ts_basicsymplectic_type 1
1030: test:
1031: suffix: coulomb_bsi_2d_2
1032: args: -ts_basicsymplectic_type 2
1033: test:
1034: suffix: coulomb_bsi_2d_3
1035: args: -ts_basicsymplectic_type 3
1036: test:
1037: suffix: coulomb_bsi_2d_4
1038: args: -ts_basicsymplectic_type 4
1040: testset:
1041: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1042: 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 \
1043: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1044: -ts_type basicsymplectic\
1045: -em_type primal -em_pc_type svd\
1046: -dm_view -output_step 50 -error -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10\
1047: -petscspace_degree 2 -petscfe_default_quadrature_order 3 -sigma 1.0e-8 -timeScale 2.0e-14
1048: test:
1049: suffix: poisson_bsi_2d_1
1050: args: -ts_basicsymplectic_type 1
1051: test:
1052: suffix: poisson_bsi_2d_2
1053: args: -ts_basicsymplectic_type 2
1054: test:
1055: suffix: poisson_bsi_2d_3
1056: args: -ts_basicsymplectic_type 3
1057: test:
1058: suffix: poisson_bsi_2d_4
1059: args: -ts_basicsymplectic_type 4
1061: testset:
1062: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1063: args: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1064: -ts_type theta -ts_theta_theta 0.5 -ts_convergence_estimate -convest_num_refine 2 \
1065: -mat_type baij -ksp_error_if_not_converged -em_pc_type svd\
1066: -dm_view -output_step 50 -error -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10\
1067: -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14
1068: test:
1069: suffix: im_2d_0
1070: 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
1072: testset:
1073: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1074: 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 \
1075: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV -dm_swarm_num_species 1\
1076: -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
1077: -em_snes_type ksponly -em_pc_type svd -em_type primal -petscspace_degree 1\
1078: -dm_view -output_step 50\
1079: -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_dt 0.01 -ts_max_time 1.0 -ts_max_steps 10
1080: test:
1081: suffix: bsi_2d_mesh_1
1082: args: -ts_basicsymplectic_type 4
1083: test:
1084: suffix: bsi_2d_mesh_1_par_2
1085: nsize: 2
1086: args: -ts_basicsymplectic_type 4
1087: test:
1088: suffix: bsi_2d_mesh_1_par_3
1089: nsize: 3
1090: args: -ts_basicsymplectic_type 4
1091: test:
1092: suffix: bsi_2d_mesh_1_par_4
1093: nsize: 4
1094: args: -ts_basicsymplectic_type 4 -dm_swarm_num_particles 0,0,2,0
1096: testset:
1097: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1098: 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 \
1099: -dm_swarm_num_particles 10 -dm_swarm_coordinate_function circleMultipleX -dm_swarm_velocity_function circleMultipleV \
1100: -ts_convergence_estimate -convest_num_refine 2 \
1101: -em_pc_type lu\
1102: -dm_view -output_step 50 -error\
1103: -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10
1104: test:
1105: suffix: bsi_2d_multiple_1
1106: args: -ts_type basicsymplectic -ts_basicsymplectic_type 1
1107: test:
1108: suffix: bsi_2d_multiple_2
1109: args: -ts_type basicsymplectic -ts_basicsymplectic_type 2
1110: test:
1111: suffix: bsi_2d_multiple_3
1112: args: -ts_type basicsymplectic -ts_basicsymplectic_type 3 -ts_dt 0.001
1113: test:
1114: suffix: im_2d_multiple_0
1115: args: -ts_type theta -ts_theta_theta 0.5 \
1116: -mat_type baij -ksp_error_if_not_converged -em_pc_type lu
1118: testset:
1119: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1120: 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 \
1121: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1122: -em_pc_type fieldsplit -ksp_rtol 1e-10 -em_ksp_type preonly -em_type mixed -em_ksp_error_if_not_converged\
1123: -dm_view -output_step 50 -error -dm_refine 0\
1124: -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10
1125: test:
1126: suffix: bsi_4_rt_1
1127: args: -ts_type basicsymplectic -ts_basicsymplectic_type 4\
1128: -pc_fieldsplit_detect_saddle_point\
1129: -pc_type fieldsplit\
1130: -pc_fieldsplit_type schur\
1131: -pc_fieldsplit_schur_precondition full \
1132: -field_petscspace_degree 2\
1133: -field_petscfe_default_quadrature_order 1\
1134: -field_petscspace_type sum \
1135: -field_petscspace_variables 2 \
1136: -field_petscspace_components 2 \
1137: -field_petscspace_sum_spaces 2 \
1138: -field_petscspace_sum_concatenate true \
1139: -field_sumcomp_0_petscspace_variables 2 \
1140: -field_sumcomp_0_petscspace_type tensor \
1141: -field_sumcomp_0_petscspace_tensor_spaces 2 \
1142: -field_sumcomp_0_petscspace_tensor_uniform false \
1143: -field_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
1144: -field_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
1145: -field_sumcomp_1_petscspace_variables 2 \
1146: -field_sumcomp_1_petscspace_type tensor \
1147: -field_sumcomp_1_petscspace_tensor_spaces 2 \
1148: -field_sumcomp_1_petscspace_tensor_uniform false \
1149: -field_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
1150: -field_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
1151: -field_petscdualspace_form_degree -1 \
1152: -field_petscdualspace_order 1 \
1153: -field_petscdualspace_lagrange_trimmed true\
1154: -ksp_gmres_restart 500
1156: TEST*/