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, user, user, 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;
268: PetscFunctionBegin;
269: PetscCall(DMGetDimension(sw, &dim));
270: PetscCall(DMSwarmGetLocalSize(sw, &Np));
272: /* Create the charges rho */
273: PetscCall(SNESGetDM(snes, &dm));
274: PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
275: PetscCall(DMGetGlobalVector(dm, &rho));
276: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
277: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
278: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
279: PetscCall(MatMultTranspose(M_p, f, rho));
280: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
281: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
282: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
283: PetscCall(MatDestroy(&M_p));
284: {
285: PetscScalar sum;
286: PetscInt n;
287: PetscReal phi_0 = 1.; /* (sigma*sigma*sigma)*(timeScale*timeScale)/(m_e*q_e*epsi_0)*/
289: /* Remove constant from rho */
290: PetscCall(VecGetSize(rho, &n));
291: PetscCall(VecSum(rho, &sum));
292: PetscCall(VecShift(rho, -sum / n));
293: PetscCall(VecSum(rho, &sum));
294: PetscCheck(PetscAbsScalar(sum) < chargeTol, PetscObjectComm((PetscObject)sw), PETSC_ERR_PLIB, "Charge should have no DC component: %g", (double)PetscRealPart(sum));
295: /* Nondimensionalize rho */
296: PetscCall(VecScale(rho, phi_0));
297: }
298: PetscCall(VecViewFromOptions(rho, NULL, "-poisson_rho_view"));
300: PetscCall(DMGetGlobalVector(dm, &phi));
301: PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
302: PetscCall(VecSet(phi, 0.0));
303: PetscCall(SNESSolve(snes, rho, phi));
304: PetscCall(DMRestoreGlobalVector(dm, &rho));
305: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
307: PetscCall(DMGetLocalVector(dm, &locPhi));
308: PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
309: PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
310: PetscCall(DMRestoreGlobalVector(dm, &phi));
312: PetscCall(DMGetDS(dm, &ds));
313: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
314: PetscCall(DMSwarmSortGetAccess(sw));
315: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
316: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
317: for (c = cStart; c < cEnd; ++c) {
318: PetscTabulation tab;
319: PetscScalar *clPhi = NULL;
320: PetscReal *pcoord, *refcoord;
321: PetscReal v[3], J[9], invJ[9], detJ;
322: PetscInt *points;
323: PetscInt Ncp, cp;
325: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
326: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
327: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
328: for (cp = 0; cp < Ncp; ++cp)
329: for (d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
330: PetscCall(DMPlexCoordinatesToReference(dm, c, Ncp, pcoord, refcoord));
331: PetscCall(PetscFECreateTabulation(fe, 1, Ncp, refcoord, 1, &tab));
332: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v, J, invJ, &detJ));
333: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
334: for (cp = 0; cp < Ncp; ++cp) {
335: const PetscReal *basisDer = tab->T[1];
336: const PetscInt p = points[cp];
338: for (d = 0; d < dim; ++d) E[p * dim + d] = 0.;
339: PetscCall(PetscFEFreeInterpolateGradient_Static(fe, basisDer, clPhi, dim, invJ, NULL, cp, &E[p * dim]));
340: for (d = 0; d < dim; ++d) E[p * dim + d] *= -1.0;
341: }
342: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
343: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
344: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
345: PetscCall(PetscTabulationDestroy(&tab));
346: PetscCall(PetscFree(points));
347: }
348: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
349: PetscCall(DMSwarmSortRestoreAccess(sw));
350: PetscCall(DMRestoreLocalVector(dm, &locPhi));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode ComputeFieldAtParticles_Mixed(SNES snes, DM sw, PetscReal E[])
355: {
356: DM dm, potential_dm;
357: IS potential_IS;
358: PetscDS ds;
359: PetscFE fe;
360: PetscFEGeom feGeometry;
361: Mat M_p;
362: Vec phi, locPhi, rho, f, temp_rho;
363: PetscQuadrature q;
364: PetscReal *coords, chargeTol = 1e-13;
365: PetscInt dim, d, cStart, cEnd, c, Np, fields = 1;
367: PetscFunctionBegin;
368: PetscCall(DMGetDimension(sw, &dim));
369: PetscCall(DMSwarmGetLocalSize(sw, &Np));
371: /* Create the charges rho */
372: PetscCall(SNESGetDM(snes, &dm));
373: PetscCall(DMGetGlobalVector(dm, &rho));
374: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
376: PetscCall(DMCreateSubDM(dm, 1, &fields, &potential_IS, &potential_dm));
377: PetscCall(DMCreateMassMatrix(sw, potential_dm, &M_p));
378: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
379: PetscCall(DMGetGlobalVector(potential_dm, &temp_rho));
380: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
381: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
382: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
383: PetscCall(MatMultTranspose(M_p, f, temp_rho));
384: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
385: PetscCall(MatDestroy(&M_p));
386: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
387: PetscCall(VecViewFromOptions(rho, NULL, "-poisson_rho_view"));
388: PetscCall(VecISCopy(rho, potential_IS, SCATTER_FORWARD, temp_rho));
389: PetscCall(DMRestoreGlobalVector(potential_dm, &temp_rho));
390: PetscCall(DMDestroy(&potential_dm));
391: PetscCall(ISDestroy(&potential_IS));
392: {
393: PetscScalar sum;
394: PetscInt n;
395: PetscReal phi_0 = 1.; /*(sigma*sigma*sigma)*(timeScale*timeScale)/(m_e*q_e*epsi_0);*/
397: /*Remove constant from rho*/
398: PetscCall(VecGetSize(rho, &n));
399: PetscCall(VecSum(rho, &sum));
400: PetscCall(VecShift(rho, -sum / n));
401: PetscCall(VecSum(rho, &sum));
402: PetscCheck(PetscAbsScalar(sum) < chargeTol, PetscObjectComm((PetscObject)sw), PETSC_ERR_PLIB, "Charge should have no DC component: %g", (double)PetscRealPart(sum));
403: /* Nondimensionalize rho */
404: PetscCall(VecScale(rho, phi_0));
405: }
406: PetscCall(DMGetGlobalVector(dm, &phi));
407: PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
408: PetscCall(VecSet(phi, 0.0));
409: PetscCall(SNESSolve(snes, NULL, phi));
410: PetscCall(DMRestoreGlobalVector(dm, &rho));
411: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
413: PetscCall(DMGetLocalVector(dm, &locPhi));
414: PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
415: PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
416: PetscCall(DMRestoreGlobalVector(dm, &phi));
418: PetscCall(DMGetDS(dm, &ds));
419: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
420: PetscCall(DMSwarmSortGetAccess(sw));
421: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
422: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
423: for (c = cStart; c < cEnd; ++c) {
424: PetscTabulation tab;
425: PetscScalar *clPhi = NULL;
426: PetscReal *pcoord, *refcoord;
427: PetscReal v[3], J[9], invJ[9], detJ;
428: PetscInt *points;
429: PetscInt Ncp, cp;
431: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
432: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
433: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
434: for (cp = 0; cp < Ncp; ++cp)
435: for (d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
436: PetscCall(DMPlexCoordinatesToReference(dm, c, Ncp, pcoord, refcoord));
437: PetscCall(PetscFECreateTabulation(fe, 1, Ncp, refcoord, 1, &tab));
438: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v, J, invJ, &detJ));
439: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
440: for (cp = 0; cp < Ncp; ++cp) {
441: const PetscInt p = points[cp];
443: for (d = 0; d < dim; ++d) E[p * dim + d] = 0.;
444: PetscCall(PetscFEGetQuadrature(fe, &q));
445: PetscCall(PetscFECreateCellGeometry(fe, q, &feGeometry));
446: PetscCall(PetscFEInterpolateAtPoints_Static(fe, tab, clPhi, &feGeometry, cp, &E[p * dim]));
447: PetscCall(PetscFEDestroyCellGeometry(fe, &feGeometry));
448: }
449: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
450: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
451: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
452: PetscCall(PetscTabulationDestroy(&tab));
453: PetscCall(PetscFree(points));
454: }
455: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
456: PetscCall(DMSwarmSortRestoreAccess(sw));
457: PetscCall(DMRestoreLocalVector(dm, &locPhi));
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: static PetscErrorCode ComputeFieldAtParticles(SNES snes, DM sw, PetscReal E[])
462: {
463: AppCtx *ctx;
464: PetscInt dim, Np;
466: PetscFunctionBegin;
469: PetscAssertPointer(E, 3);
470: PetscCall(DMGetDimension(sw, &dim));
471: PetscCall(DMSwarmGetLocalSize(sw, &Np));
472: PetscCall(DMGetApplicationContext(sw, &ctx));
473: PetscCall(PetscArrayzero(E, Np * dim));
475: switch (ctx->em) {
476: case EM_PRIMAL:
477: PetscCall(ComputeFieldAtParticles_Primal(snes, sw, E));
478: break;
479: case EM_COULOMB:
480: PetscCall(ComputeFieldAtParticles_Coulomb(snes, sw, E));
481: break;
482: case EM_MIXED:
483: PetscCall(ComputeFieldAtParticles_Mixed(snes, sw, E));
484: break;
485: case EM_NONE:
486: break;
487: default:
488: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No solver for electrostatic model %s", EMTypes[ctx->em]);
489: }
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx)
495: {
496: DM sw;
497: SNES snes = ((AppCtx *)ctx)->snes;
498: const PetscReal *coords, *vel;
499: const PetscScalar *u;
500: PetscScalar *g;
501: PetscReal *E;
502: PetscInt dim, d, Np, p;
504: PetscFunctionBeginUser;
505: PetscCall(TSGetDM(ts, &sw));
506: PetscCall(DMGetDimension(sw, &dim));
507: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
508: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
509: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
510: PetscCall(VecGetLocalSize(U, &Np));
511: PetscCall(VecGetArrayRead(U, &u));
512: PetscCall(VecGetArray(G, &g));
514: PetscCall(ComputeFieldAtParticles(snes, sw, E));
516: Np /= 2 * dim;
517: for (p = 0; p < Np; ++p) {
518: const PetscReal x0 = coords[p * dim + 0];
519: const PetscReal vy0 = vel[p * dim + 1];
520: const PetscReal omega = vy0 / x0;
522: for (d = 0; d < dim; ++d) {
523: g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d];
524: g[(p * 2 + 1) * dim + d] = E[p * dim + d] - PetscSqr(omega) * u[(p * 2 + 0) * dim + d];
525: }
526: }
527: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
528: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
529: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
530: PetscCall(VecRestoreArrayRead(U, &u));
531: PetscCall(VecRestoreArray(G, &g));
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: /* J_{ij} = dF_i/dx_j
536: J_p = ( 0 1)
537: (-w^2 0)
538: TODO Now there is another term with w^2 from the electric field. I think we will need to invert the operator.
539: Perhaps we can approximate the Jacobian using only the cellwise P-P gradient from Coulomb
540: */
541: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, void *ctx)
542: {
543: DM sw;
544: const PetscReal *coords, *vel;
545: PetscInt dim, d, Np, p, rStart;
547: PetscFunctionBeginUser;
548: PetscCall(TSGetDM(ts, &sw));
549: PetscCall(DMGetDimension(sw, &dim));
550: PetscCall(VecGetLocalSize(U, &Np));
551: PetscCall(MatGetOwnershipRange(J, &rStart, NULL));
552: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
553: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
554: Np /= 2 * dim;
555: for (p = 0; p < Np; ++p) {
556: const PetscReal x0 = coords[p * dim + 0];
557: const PetscReal vy0 = vel[p * dim + 1];
558: const PetscReal omega = vy0 / x0;
559: PetscScalar vals[4] = {0., 1., -PetscSqr(omega), 0.};
561: for (d = 0; d < dim; ++d) {
562: const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
563: PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES));
564: }
565: }
566: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
567: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
568: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
569: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
574: {
575: DM sw;
576: const PetscScalar *v;
577: PetscScalar *xres;
578: PetscInt Np, p, dim, d;
580: PetscFunctionBeginUser;
581: PetscCall(TSGetDM(ts, &sw));
582: PetscCall(DMGetDimension(sw, &dim));
583: PetscCall(VecGetLocalSize(Xres, &Np));
584: Np /= dim;
585: PetscCall(VecGetArrayRead(V, &v));
586: PetscCall(VecGetArray(Xres, &xres));
587: for (p = 0; p < Np; ++p) {
588: for (d = 0; d < dim; ++d) xres[p * dim + d] = v[p * dim + d];
589: }
590: PetscCall(VecRestoreArrayRead(V, &v));
591: PetscCall(VecRestoreArray(Xres, &xres));
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
596: {
597: DM sw;
598: SNES snes = ((AppCtx *)ctx)->snes;
599: const PetscScalar *x;
600: const PetscReal *coords, *vel;
601: PetscScalar *vres;
602: PetscReal *E;
603: PetscInt Np, p, dim, d;
605: PetscFunctionBeginUser;
606: PetscCall(TSGetDM(ts, &sw));
607: PetscCall(DMGetDimension(sw, &dim));
608: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
609: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
610: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
611: PetscCall(VecGetLocalSize(Vres, &Np));
612: PetscCall(VecGetArrayRead(X, &x));
613: PetscCall(VecGetArray(Vres, &vres));
614: PetscCheck(dim == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension must be 2");
616: PetscCall(ComputeFieldAtParticles(snes, sw, E));
618: Np /= dim;
619: for (p = 0; p < Np; ++p) {
620: const PetscReal x0 = coords[p * dim + 0];
621: const PetscReal vy0 = vel[p * dim + 1];
622: const PetscReal omega = vy0 / x0;
624: for (d = 0; d < dim; ++d) vres[p * dim + d] = E[p * dim + d] - PetscSqr(omega) * x[p * dim + d];
625: }
626: PetscCall(VecRestoreArrayRead(X, &x));
627: PetscCall(VecRestoreArray(Vres, &vres));
628: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
629: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
630: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: static PetscErrorCode CreateSolution(TS ts)
635: {
636: DM sw;
637: Vec u;
638: PetscInt dim, Np;
640: PetscFunctionBegin;
641: PetscCall(TSGetDM(ts, &sw));
642: PetscCall(DMGetDimension(sw, &dim));
643: PetscCall(DMSwarmGetLocalSize(sw, &Np));
644: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
645: PetscCall(VecSetBlockSize(u, dim));
646: PetscCall(VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE));
647: PetscCall(VecSetUp(u));
648: PetscCall(TSSetSolution(ts, u));
649: PetscCall(VecDestroy(&u));
650: PetscFunctionReturn(PETSC_SUCCESS);
651: }
653: static PetscErrorCode SetProblem(TS ts)
654: {
655: AppCtx *user;
656: DM sw;
658: PetscFunctionBegin;
659: PetscCall(TSGetDM(ts, &sw));
660: PetscCall(DMGetApplicationContext(sw, (void **)&user));
661: // Define unified system for (X, V)
662: {
663: Mat J;
664: PetscInt dim, Np;
666: PetscCall(DMGetDimension(sw, &dim));
667: PetscCall(DMSwarmGetLocalSize(sw, &Np));
668: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
669: PetscCall(MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE));
670: PetscCall(MatSetBlockSize(J, 2 * dim));
671: PetscCall(MatSetFromOptions(J));
672: PetscCall(MatSetUp(J));
673: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, user));
674: PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, user));
675: PetscCall(MatDestroy(&J));
676: }
677: /* Define split system for X and V */
678: {
679: Vec u;
680: IS isx, isv, istmp;
681: const PetscInt *idx;
682: PetscInt dim, Np, rstart;
684: PetscCall(TSGetSolution(ts, &u));
685: PetscCall(DMGetDimension(sw, &dim));
686: PetscCall(DMSwarmGetLocalSize(sw, &Np));
687: PetscCall(VecGetOwnershipRange(u, &rstart, NULL));
688: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 0, 2, &istmp));
689: PetscCall(ISGetIndices(istmp, &idx));
690: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isx));
691: PetscCall(ISRestoreIndices(istmp, &idx));
692: PetscCall(ISDestroy(&istmp));
693: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 1, 2, &istmp));
694: PetscCall(ISGetIndices(istmp, &idx));
695: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isv));
696: PetscCall(ISRestoreIndices(istmp, &idx));
697: PetscCall(ISDestroy(&istmp));
698: PetscCall(TSRHSSplitSetIS(ts, "position", isx));
699: PetscCall(TSRHSSplitSetIS(ts, "momentum", isv));
700: PetscCall(ISDestroy(&isx));
701: PetscCall(ISDestroy(&isv));
702: PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, user));
703: PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, user));
704: }
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: PetscErrorCode circleSingleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
709: {
710: x[0] = p + 1.;
711: x[1] = 0.;
712: return PETSC_SUCCESS;
713: }
715: PetscErrorCode circleSingleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *ctx)
716: {
717: v[0] = 0.;
718: v[1] = PetscSqrtReal(1000. / (p + 1.));
719: return PETSC_SUCCESS;
720: }
722: /* Put 5 particles into each circle */
723: PetscErrorCode circleMultipleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
724: {
725: const PetscInt n = 5;
726: const PetscReal r0 = (p / n) + 1.;
727: const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
729: x[0] = r0 * PetscCosReal(th0);
730: x[1] = r0 * PetscSinReal(th0);
731: return PETSC_SUCCESS;
732: }
734: /* Put 5 particles into each circle */
735: PetscErrorCode circleMultipleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *ctx)
736: {
737: const PetscInt n = 5;
738: const PetscReal r0 = (p / n) + 1.;
739: const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
740: const PetscReal omega = PetscSqrtReal(1000. / r0) / r0;
742: v[0] = -r0 * omega * PetscSinReal(th0);
743: v[1] = r0 * omega * PetscCosReal(th0);
744: return PETSC_SUCCESS;
745: }
747: /*
748: InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values.
750: Input Parameters:
751: + ts - The TS
752: - useInitial - Flag to also set the initial conditions to the current coordinates and velocities and setup the problem
754: Output Parameter:
755: . u - The initialized solution vector
757: Level: advanced
759: .seealso: InitializeSolve()
760: */
761: static PetscErrorCode InitializeSolveAndSwarm(TS ts, PetscBool useInitial)
762: {
763: DM sw;
764: Vec u, gc, gv, gc0, gv0;
765: IS isx, isv;
766: AppCtx *user;
768: PetscFunctionBeginUser;
769: PetscCall(TSGetDM(ts, &sw));
770: PetscCall(DMGetApplicationContext(sw, &user));
771: if (useInitial) {
772: PetscReal v0[1] = {1.};
774: PetscCall(DMSwarmInitializeCoordinates(sw));
775: PetscCall(DMSwarmInitializeVelocitiesFromOptions(sw, v0));
776: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
777: PetscCall(TSReset(ts));
778: PetscCall(CreateSolution(ts));
779: PetscCall(SetProblem(ts));
780: }
781: PetscCall(TSGetSolution(ts, &u));
782: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
783: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
784: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
785: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initCoordinates", &gc0));
786: if (useInitial) PetscCall(VecCopy(gc, gc0));
787: PetscCall(VecISCopy(u, isx, SCATTER_FORWARD, gc));
788: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
789: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initCoordinates", &gc0));
790: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
791: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initVelocity", &gv0));
792: if (useInitial) PetscCall(VecCopy(gv, gv0));
793: PetscCall(VecISCopy(u, isv, SCATTER_FORWARD, gv));
794: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
795: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initVelocity", &gv0));
796: PetscFunctionReturn(PETSC_SUCCESS);
797: }
799: static PetscErrorCode InitializeSolve(TS ts, Vec u)
800: {
801: PetscFunctionBegin;
802: PetscCall(TSSetSolution(ts, u));
803: PetscCall(InitializeSolveAndSwarm(ts, PETSC_TRUE));
804: PetscFunctionReturn(PETSC_SUCCESS);
805: }
807: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
808: {
809: MPI_Comm comm;
810: DM sw;
811: AppCtx *user;
812: const PetscScalar *u;
813: const PetscReal *coords, *vel;
814: PetscScalar *e;
815: PetscReal t;
816: PetscInt dim, Np, p;
818: PetscFunctionBeginUser;
819: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
820: PetscCall(TSGetDM(ts, &sw));
821: PetscCall(DMGetApplicationContext(sw, &user));
822: PetscCall(DMGetDimension(sw, &dim));
823: PetscCall(TSGetSolveTime(ts, &t));
824: PetscCall(VecGetArray(E, &e));
825: PetscCall(VecGetArrayRead(U, &u));
826: PetscCall(VecGetLocalSize(U, &Np));
827: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
828: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
829: Np /= 2 * dim;
830: for (p = 0; p < Np; ++p) {
831: /* TODO generalize initial conditions and project into plane instead of assuming x-y */
832: const PetscReal r0 = DMPlex_NormD_Internal(dim, &coords[p * dim]);
833: const PetscReal th0 = PetscAtan2Real(coords[p * dim + 1], coords[p * dim + 0]);
834: const PetscReal v0 = DMPlex_NormD_Internal(dim, &vel[p * dim]);
835: const PetscReal omega = v0 / r0;
836: const PetscReal ct = PetscCosReal(omega * t + th0);
837: const PetscReal st = PetscSinReal(omega * t + th0);
838: const PetscScalar *x = &u[(p * 2 + 0) * dim];
839: const PetscScalar *v = &u[(p * 2 + 1) * dim];
840: const PetscReal xe[3] = {r0 * ct, r0 * st, 0.0};
841: const PetscReal ve[3] = {-v0 * st, v0 * ct, 0.0};
842: PetscInt d;
844: for (d = 0; d < dim; ++d) {
845: e[(p * 2 + 0) * dim + d] = x[d] - xe[d];
846: e[(p * 2 + 1) * dim + d] = v[d] - ve[d];
847: }
848: if (user->error) {
849: const PetscReal en = 0.5 * DMPlex_DotRealD_Internal(dim, v, v);
850: const PetscReal exen = 0.5 * PetscSqr(v0);
851: 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)));
852: }
853: }
854: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
855: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
856: PetscCall(VecRestoreArrayRead(U, &u));
857: PetscCall(VecRestoreArray(E, &e));
858: PetscFunctionReturn(PETSC_SUCCESS);
859: }
861: static PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
862: {
863: const PetscInt ostep = ((AppCtx *)ctx)->ostep;
864: const EMType em = ((AppCtx *)ctx)->em;
865: DM sw;
866: const PetscScalar *u;
867: PetscReal *coords, *E;
868: PetscReal enKin = 0., enEM = 0.;
869: PetscInt dim, d, Np, p, q;
871: PetscFunctionBeginUser;
872: if (step % ostep == 0) {
873: PetscCall(TSGetDM(ts, &sw));
874: PetscCall(DMGetDimension(sw, &dim));
875: PetscCall(VecGetArrayRead(U, &u));
876: PetscCall(VecGetLocalSize(U, &Np));
877: Np /= 2 * dim;
878: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
879: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
880: if (!step) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Time Step Part Energy\n"));
881: for (p = 0; p < Np; ++p) {
882: const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);
883: PetscReal *pcoord = &coords[p * dim];
885: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " %5" PetscInt_FMT " %10.4lf\n", (double)t, step, p, (double)(0.5 * v2)));
886: enKin += 0.5 * v2;
887: if (em == EM_NONE) {
888: continue;
889: } else if (em == EM_COULOMB) {
890: for (q = p + 1; q < Np; ++q) {
891: PetscReal *qcoord = &coords[q * dim];
892: PetscReal rpq[3], r;
893: for (d = 0; d < dim; ++d) rpq[d] = pcoord[d] - qcoord[d];
894: r = DMPlex_NormD_Internal(dim, rpq);
895: enEM += 1. / r;
896: }
897: } else if (em == EM_PRIMAL || em == EM_MIXED) {
898: for (d = 0; d < dim; ++d) enEM += E[p * dim + d];
899: }
900: }
901: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " KE\t %10.4lf\n", (double)t, step, (double)enKin));
902: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " PE\t %1.10g\n", (double)t, step, (double)enEM));
903: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " E\t %10.4lf\n", (double)t, step, (double)(enKin + enEM)));
904: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
905: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
906: PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)ts), NULL));
907: PetscCall(VecRestoreArrayRead(U, &u));
908: }
909: PetscFunctionReturn(PETSC_SUCCESS);
910: }
912: static PetscErrorCode SetUpMigrateParticles(TS ts, PetscInt n, PetscReal t, Vec x, PetscBool *flg, void *ctx)
913: {
914: DM sw;
916: PetscFunctionBeginUser;
917: *flg = PETSC_TRUE;
918: PetscCall(TSGetDM(ts, &sw));
919: PetscCall(DMViewFromOptions(sw, NULL, "-migrate_view_pre"));
920: {
921: Vec u, gc, gv;
922: IS isx, isv;
924: PetscCall(TSGetSolution(ts, &u));
925: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
926: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
927: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
928: PetscCall(VecISCopy(u, isx, SCATTER_REVERSE, gc));
929: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
930: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
931: PetscCall(VecISCopy(u, isv, SCATTER_REVERSE, gv));
932: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
933: }
934: PetscFunctionReturn(PETSC_SUCCESS);
935: }
937: static PetscErrorCode MigrateParticles(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx)
938: {
939: DM sw;
941: PetscFunctionBeginUser;
942: PetscCall(TSGetDM(ts, &sw));
943: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
944: PetscCall(CreateSolution(ts));
945: PetscCall(SetProblem(ts));
946: PetscCall(InitializeSolveAndSwarm(ts, PETSC_FALSE));
947: PetscFunctionReturn(PETSC_SUCCESS);
948: }
950: int main(int argc, char **argv)
951: {
952: DM dm, sw;
953: TS ts;
954: Vec u;
955: AppCtx user;
957: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
958: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
959: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
960: PetscCall(CreatePoisson(dm, &user));
961: PetscCall(CreateSwarm(dm, &user, &sw));
962: PetscCall(DMSetApplicationContext(sw, &user));
964: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
965: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
966: PetscCall(TSSetDM(ts, sw));
967: PetscCall(TSSetMaxTime(ts, 0.1));
968: PetscCall(TSSetTimeStep(ts, 0.00001));
969: PetscCall(TSSetMaxSteps(ts, 100));
970: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
971: PetscCall(TSMonitorSet(ts, EnergyMonitor, &user, NULL));
972: PetscCall(TSSetFromOptions(ts));
973: PetscCall(DMSwarmVectorDefineField(sw, "velocity"));
974: PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
975: PetscCall(TSSetComputeExactError(ts, ComputeError));
976: PetscCall(TSSetResize(ts, SetUpMigrateParticles, MigrateParticles, NULL));
978: PetscCall(CreateSolution(ts));
979: PetscCall(TSGetSolution(ts, &u));
980: PetscCall(TSComputeInitialCondition(ts, u));
981: PetscCall(TSSolve(ts, NULL));
983: PetscCall(SNESDestroy(&user.snes));
984: PetscCall(TSDestroy(&ts));
985: PetscCall(DMDestroy(&sw));
986: PetscCall(DMDestroy(&dm));
987: PetscCall(PetscFinalize());
988: return 0;
989: }
991: /*TEST
993: build:
994: requires: double !complex
996: testset:
997: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
998: 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 \
999: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1000: -ts_type basicsymplectic\
1001: -dm_view -output_step 50 -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10
1002: test:
1003: suffix: none_bsi_2d_1
1004: args: -ts_basicsymplectic_type 1 -em_type none -error
1005: test:
1006: suffix: none_bsi_2d_2
1007: args: -ts_basicsymplectic_type 2 -em_type none -error
1008: test:
1009: suffix: none_bsi_2d_3
1010: args: -ts_basicsymplectic_type 3 -em_type none -error
1011: test:
1012: suffix: none_bsi_2d_4
1013: args: -ts_basicsymplectic_type 4 -em_type none -error
1014: test:
1015: suffix: coulomb_bsi_2d_1
1016: args: -ts_basicsymplectic_type 1
1017: test:
1018: suffix: coulomb_bsi_2d_2
1019: args: -ts_basicsymplectic_type 2
1020: test:
1021: suffix: coulomb_bsi_2d_3
1022: args: -ts_basicsymplectic_type 3
1023: test:
1024: suffix: coulomb_bsi_2d_4
1025: args: -ts_basicsymplectic_type 4
1027: testset:
1028: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1029: 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 \
1030: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1031: -ts_type basicsymplectic\
1032: -em_type primal -em_pc_type svd\
1033: -dm_view -output_step 50 -error -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10\
1034: -petscspace_degree 2 -petscfe_default_quadrature_order 3 -sigma 1.0e-8 -timeScale 2.0e-14
1035: test:
1036: suffix: poisson_bsi_2d_1
1037: args: -ts_basicsymplectic_type 1
1038: test:
1039: suffix: poisson_bsi_2d_2
1040: args: -ts_basicsymplectic_type 2
1041: test:
1042: suffix: poisson_bsi_2d_3
1043: args: -ts_basicsymplectic_type 3
1044: test:
1045: suffix: poisson_bsi_2d_4
1046: args: -ts_basicsymplectic_type 4
1048: testset:
1049: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1050: args: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1051: -ts_type theta -ts_theta_theta 0.5 -ts_convergence_estimate -convest_num_refine 2 \
1052: -mat_type baij -ksp_error_if_not_converged -em_pc_type svd\
1053: -dm_view -output_step 50 -error -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10\
1054: -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14
1055: test:
1056: suffix: im_2d_0
1057: 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
1059: testset:
1060: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1061: 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 \
1062: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV -dm_swarm_num_species 1\
1063: -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
1064: -em_snes_type ksponly -em_pc_type svd -em_type primal -petscspace_degree 1\
1065: -dm_view -output_step 50\
1066: -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_dt 0.01 -ts_max_time 1.0 -ts_max_steps 10
1067: test:
1068: suffix: bsi_2d_mesh_1
1069: args: -ts_basicsymplectic_type 4
1070: test:
1071: suffix: bsi_2d_mesh_1_par_2
1072: nsize: 2
1073: args: -ts_basicsymplectic_type 4
1074: test:
1075: suffix: bsi_2d_mesh_1_par_3
1076: nsize: 3
1077: args: -ts_basicsymplectic_type 4
1078: test:
1079: suffix: bsi_2d_mesh_1_par_4
1080: nsize: 4
1081: args: -ts_basicsymplectic_type 4 -dm_swarm_num_particles 0,0,2,0
1083: testset:
1084: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1085: 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 \
1086: -dm_swarm_num_particles 10 -dm_swarm_coordinate_function circleMultipleX -dm_swarm_velocity_function circleMultipleV \
1087: -ts_convergence_estimate -convest_num_refine 2 \
1088: -em_pc_type lu\
1089: -dm_view -output_step 50 -error\
1090: -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10
1091: test:
1092: suffix: bsi_2d_multiple_1
1093: args: -ts_type basicsymplectic -ts_basicsymplectic_type 1
1094: test:
1095: suffix: bsi_2d_multiple_2
1096: args: -ts_type basicsymplectic -ts_basicsymplectic_type 2
1097: test:
1098: suffix: bsi_2d_multiple_3
1099: args: -ts_type basicsymplectic -ts_basicsymplectic_type 3 -ts_dt 0.001
1100: test:
1101: suffix: im_2d_multiple_0
1102: args: -ts_type theta -ts_theta_theta 0.5 \
1103: -mat_type baij -ksp_error_if_not_converged -em_pc_type lu
1105: testset:
1106: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1107: 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 \
1108: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1109: -em_pc_type fieldsplit -ksp_rtol 1e-10 -em_ksp_type preonly -em_type mixed -em_ksp_error_if_not_converged\
1110: -dm_view -output_step 50 -error -dm_refine 0\
1111: -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10
1112: test:
1113: suffix: bsi_4_rt_1
1114: args: -ts_type basicsymplectic -ts_basicsymplectic_type 4\
1115: -pc_fieldsplit_detect_saddle_point\
1116: -pc_type fieldsplit\
1117: -pc_fieldsplit_type schur\
1118: -pc_fieldsplit_schur_precondition full \
1119: -field_petscspace_degree 2\
1120: -field_petscfe_default_quadrature_order 1\
1121: -field_petscspace_type sum \
1122: -field_petscspace_variables 2 \
1123: -field_petscspace_components 2 \
1124: -field_petscspace_sum_spaces 2 \
1125: -field_petscspace_sum_concatenate true \
1126: -field_sumcomp_0_petscspace_variables 2 \
1127: -field_sumcomp_0_petscspace_type tensor \
1128: -field_sumcomp_0_petscspace_tensor_spaces 2 \
1129: -field_sumcomp_0_petscspace_tensor_uniform false \
1130: -field_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
1131: -field_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
1132: -field_sumcomp_1_petscspace_variables 2 \
1133: -field_sumcomp_1_petscspace_type tensor \
1134: -field_sumcomp_1_petscspace_tensor_spaces 2 \
1135: -field_sumcomp_1_petscspace_tensor_uniform false \
1136: -field_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
1137: -field_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
1138: -field_petscdualspace_form_degree -1 \
1139: -field_petscdualspace_order 1 \
1140: -field_petscdualspace_lagrange_trimmed true\
1141: -ksp_gmres_restart 500
1143: TEST*/