Actual source code: ex3.c
1: static char help[] = "Landau Damping test using Vlasov-Poisson equations\n";
3: /*
4: This example solves the Vlasov-Poisson system for Landau damping (6X + 6V).
5: The system is solved using a Particle-In-Cell (PIC) method with DMSwarm for particles and DMPlex/PetscFE for the field.
6: This particular example uses the velocity mesh from DMPlexLandauCreateVelocitySpace for 3D velocity space.
8: Options:
9: -particle_monitor [prefix] : Output particle data (x, v, w, E) to binary files at each output step.
10: Optional prefix for filenames (default: "particles").
12: */
13: #include <petscts.h>
14: #include <petscdmplex.h>
15: #include <petscdmswarm.h>
16: #include <petscfe.h>
17: #include <petscds.h>
18: #include <petscbag.h>
19: #include <petscdraw.h>
20: #include <petsclandau.h>
21: #include <petscdmcomposite.h>
22: #include <petsc/private/dmpleximpl.h>
23: #include <petsc/private/petscfeimpl.h>
24: #include <petsc/private/dmswarmimpl.h>
25: #include "petscdm.h"
26: #include "petscdmlabel.h"
28: PETSC_EXTERN PetscErrorCode stream(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
29: PETSC_EXTERN PetscErrorCode line(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
31: typedef enum {
32: V0,
33: X0,
34: T0,
35: M0,
36: Q0,
37: PHI0,
38: POISSON,
39: VLASOV,
40: SIGMA,
41: NUM_CONSTANTS
42: } ConstantType;
43: typedef struct {
44: PetscScalar v0; /* Velocity scale, often the thermal velocity */
45: PetscScalar t0; /* Time scale */
46: PetscScalar x0; /* Space scale */
47: PetscScalar m0; /* Mass scale */
48: PetscScalar q0; /* Charge scale */
49: PetscScalar kb;
50: PetscScalar epsi0;
51: PetscScalar phi0; /* Potential scale */
52: PetscScalar poissonNumber; /* Non-Dimensional Poisson Number */
53: PetscScalar vlasovNumber; /* Non-Dimensional Vlasov Number */
54: PetscReal sigma; /* Nondimensional charge per length in x */
55: } Parameter;
57: typedef struct {
58: PetscBag bag; // Problem parameters
59: PetscBool error; // Flag for printing the error
60: PetscBool efield_monitor; // Flag to show electric field monitor
61: PetscBool moment_monitor; // Flag to show distribution moment monitor
62: char particle_monitor_prefix[PETSC_MAX_PATH_LEN];
63: PetscBool particle_monitor; // Flag to output particle data
64: PetscInt ostep; // Print the energy at each ostep time steps
65: PetscInt numParticles;
66: PetscReal timeScale; /* Nondimensionalizing time scale */
67: PetscReal charges[2]; /* The charges of each species */
68: PetscReal masses[2]; /* The masses of each species */
69: PetscReal thermal_energy[2]; /* Thermal Energy (used to get other constants)*/
70: PetscReal cosine_coefficients[2]; /*(alpha, k)*/
71: PetscReal totalWeight;
72: PetscReal stepSize;
73: PetscInt steps;
74: PetscReal initVel;
75: SNES snes; // EM solver
76: DM dmPot; // The DM for potential
77: Mat M; // The finite element mass matrix for potential
78: PetscFEGeom *fegeom; // Geometric information for the DM cells
79: PetscBool validE; // Flag to indicate E-field in swarm is valid
80: PetscReal drawlgEmin; // The minimum lg(E) to plot
81: PetscDrawLG drawlgE; // Logarithm of maximum electric field
82: DM swarm;
83: PetscBool checkweights;
84: PetscInt checkVRes; // Flag to check/output velocity residuals for nightly tests
85: DM landau_pack;
86: PetscBool use_landau_velocity_space;
87: PetscLogEvent RhsXEvent, RhsVEvent, ESolveEvent, ETabEvent;
88: } AppCtx;
90: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
91: {
92: PetscFunctionBeginUser;
93: PetscInt d = 2;
94: PetscInt maxSpecies = 2;
95: options->error = PETSC_FALSE;
96: options->efield_monitor = PETSC_FALSE;
97: options->moment_monitor = PETSC_FALSE;
98: options->particle_monitor = PETSC_FALSE;
99: options->ostep = 100;
100: options->timeScale = 2.0e-14;
101: options->charges[0] = -1.0;
102: options->charges[1] = 1.0;
103: options->masses[0] = 1.0;
104: options->masses[1] = 1000.0;
105: options->thermal_energy[0] = 1.0;
106: options->thermal_energy[1] = 1.0;
107: options->cosine_coefficients[0] = 0.01;
108: options->cosine_coefficients[1] = 0.5;
109: options->initVel = 1;
110: options->totalWeight = 1.0;
111: options->validE = PETSC_FALSE;
112: options->drawlgEmin = -6;
113: options->drawlgE = NULL;
114: options->snes = NULL;
115: options->dmPot = NULL;
116: options->M = NULL;
117: options->numParticles = 32768;
118: options->checkweights = PETSC_FALSE;
119: options->checkVRes = 0;
120: options->landau_pack = NULL;
121: options->use_landau_velocity_space = PETSC_FALSE;
123: PetscOptionsBegin(comm, "", "Landau Damping options", "DMSWARM");
124: PetscCall(PetscOptionsBool("-error", "Flag to print the error", "ex3.c", options->error, &options->error, NULL));
125: PetscCall(PetscOptionsBool("-efield_monitor", "Flag to plot log(max E) over time", "ex3.c", options->efield_monitor, &options->efield_monitor, NULL));
126: PetscCall(PetscOptionsReal("-efield_min_monitor", "Minimum E field to plot", "ex3.c", options->drawlgEmin, &options->drawlgEmin, NULL));
127: PetscCall(PetscOptionsBool("-moments_monitor", "Flag to show moments table", "ex3.c", options->moment_monitor, &options->moment_monitor, NULL));
128: PetscCall(PetscOptionsString("-particle_monitor", "Prefix for particle data files", "ex3.c", options->particle_monitor_prefix, options->particle_monitor_prefix, sizeof(options->particle_monitor_prefix), &options->particle_monitor));
129: PetscCall(PetscOptionsBool("-check_weights", "Ensure all particle weights are positive", "ex3.c", options->checkweights, &options->checkweights, NULL));
130: PetscCall(PetscOptionsInt("-check_vel_res", "Check particle velocity residuals for nightly tests", "ex3.c", options->checkVRes, &options->checkVRes, NULL));
131: PetscCall(PetscOptionsInt("-output_step", "Number of time steps between output", "ex3.c", options->ostep, &options->ostep, NULL));
132: PetscCall(PetscOptionsReal("-timeScale", "Nondimensionalizing time scale", "ex3.c", options->timeScale, &options->timeScale, NULL));
133: PetscCall(PetscOptionsReal("-initial_velocity", "Initial velocity of perturbed particle", "ex3.c", options->initVel, &options->initVel, NULL));
134: PetscCall(PetscOptionsReal("-total_weight", "Total weight of all particles", "ex3.c", options->totalWeight, &options->totalWeight, NULL));
135: PetscCall(PetscOptionsRealArray("-cosine_coefficients", "Amplitude and frequency of cosine equation used in initialization", "ex3.c", options->cosine_coefficients, &d, NULL));
136: PetscCall(PetscOptionsRealArray("-charges", "Species charges", "ex3.c", options->charges, &maxSpecies, NULL));
137: PetscCall(PetscOptionsBool("-use_landau_velocity_space", "Use Landau velocity space", "ex3.c", options->use_landau_velocity_space, &options->use_landau_velocity_space, NULL));
138: PetscOptionsEnd();
140: PetscCall(PetscLogEventRegister("RhsX", TS_CLASSID, &options->RhsXEvent));
141: PetscCall(PetscLogEventRegister("RhsV", TS_CLASSID, &options->RhsVEvent));
142: PetscCall(PetscLogEventRegister("ESolve", TS_CLASSID, &options->ESolveEvent));
143: PetscCall(PetscLogEventRegister("ETab", TS_CLASSID, &options->ETabEvent));
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: static PetscErrorCode SetupContext(DM dm, DM sw, AppCtx *user)
148: {
149: MPI_Comm comm;
151: PetscFunctionBeginUser;
152: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
153: if (user->efield_monitor) {
154: PetscDraw draw;
155: PetscDrawAxis axis;
157: PetscCall(PetscDrawCreate(comm, NULL, "Max Electric Field", 0, 300, 400, 300, &draw));
158: PetscCall(PetscDrawSetSave(draw, "ex3_Efield"));
159: PetscCall(PetscDrawSetFromOptions(draw));
160: PetscCall(PetscDrawLGCreate(draw, 1, &user->drawlgE));
161: PetscCall(PetscDrawDestroy(&draw));
162: PetscCall(PetscDrawLGGetAxis(user->drawlgE, &axis));
163: PetscCall(PetscDrawAxisSetLabels(axis, "Electron Electric Field", "time", "E_max"));
164: PetscCall(PetscDrawLGSetLimits(user->drawlgE, 0., user->steps * user->stepSize, user->drawlgEmin, 0.));
165: PetscCall(PetscDrawLGSetUseMarkers(user->drawlgE, PETSC_FALSE));
166: }
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: static PetscErrorCode DestroyContext(AppCtx *user)
171: {
172: PetscFunctionBeginUser;
173: if (user->landau_pack) PetscCall(DMPlexLandauDestroyVelocitySpace(&user->landau_pack));
174: PetscCall(PetscDrawLGDestroy(&user->drawlgE));
175: PetscCall(PetscBagDestroy(&user->bag));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: static PetscErrorCode CheckNonNegativeWeights(DM sw, AppCtx *user)
180: {
181: const PetscScalar *w;
182: PetscInt Np;
184: PetscFunctionBeginUser;
185: if (!user->checkweights) PetscFunctionReturn(PETSC_SUCCESS);
186: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&w));
187: PetscCall(DMSwarmGetLocalSize(sw, &Np));
188: for (PetscInt p = 0; p < Np; ++p) PetscCheck(w[p] >= 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Particle %" PetscInt_FMT " has negative weight %g", p, w[p]);
189: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: static void f0_Dirichlet(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[])
194: {
195: for (PetscInt d = 0; d < dim; ++d) f0[0] += 0.5 * PetscSqr(u_x[d]);
196: }
198: static PetscErrorCode computeFieldEnergy(DM dm, Vec u, PetscReal *En)
199: {
200: PetscDS ds;
201: const PetscInt field = 0;
202: PetscInt Nf;
203: void *ctx;
205: PetscFunctionBegin;
206: PetscCall(DMGetApplicationContext(dm, &ctx));
207: PetscCall(DMGetDS(dm, &ds));
208: PetscCall(PetscDSGetNumFields(ds, &Nf));
209: PetscCheck(Nf == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "We currently only support 1 field, not %" PetscInt_FMT, Nf);
210: PetscCall(PetscDSSetObjective(ds, field, &f0_Dirichlet));
211: PetscCall(DMPlexComputeIntegralFEM(dm, u, En, ctx));
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: static void f0_grad_phi2(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[])
216: {
217: f0[0] = 0.;
218: for (PetscInt d = 0; d < dim; ++d) f0[0] += PetscSqr(u_x[uOff_x[0] + d]); // + d * dim cause segfault
219: }
221: static PetscErrorCode MonitorEField(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
222: {
223: AppCtx *user = (AppCtx *)ctx;
224: DM sw;
225: PetscScalar intESq;
226: PetscReal *E, *x, *weight;
227: PetscReal Enorm = 0., lgEnorm, lgEmax, sum = 0., Emax = 0., chargesum = 0.;
228: PetscReal pmoments[16]; /* \int f, \int v f, \int v^2 f */
229: PetscInt *species, dim, Np, gNp;
230: MPI_Comm comm;
232: PetscFunctionBeginUser;
233: if (step < 0 || !user->validE) PetscFunctionReturn(PETSC_SUCCESS);
234: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
235: PetscCall(TSGetDM(ts, &sw));
236: PetscCall(DMGetDimension(sw, &dim));
237: PetscCall(DMSwarmGetLocalSize(sw, &Np));
238: PetscCall(DMSwarmGetSize(sw, &gNp));
239: PetscCall(DMSwarmSortGetAccess(sw));
240: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
241: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
242: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
243: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
245: for (PetscInt p = 0; p < Np; ++p) {
246: PetscReal Emag = 0.;
247: for (PetscInt d = 0; d < dim; ++d) {
248: PetscReal temp = PetscAbsReal(E[p * dim + d]);
249: if (temp > Emax) Emax = temp;
250: Emag += PetscSqr(E[p * dim + d]);
251: }
252: Enorm += PetscSqrtReal(Emag);
253: sum += E[p * dim];
254: chargesum += user->charges[0] * weight[p];
255: }
256: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &Emax, 1, MPIU_REAL, MPIU_MAX, comm));
257: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &Enorm, 1, MPIU_REAL, MPIU_SUM, comm));
258: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &sum, 1, MPIU_REAL, MPIU_SUM, comm));
259: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &chargesum, 1, MPIU_REAL, MPIU_SUM, comm));
260: lgEnorm = Enorm != 0 ? PetscLog10Real(Enorm) : -16.;
261: lgEmax = Emax != 0 ? PetscLog10Real(Emax) : user->drawlgEmin;
262: if (lgEmax < user->drawlgEmin) lgEmax = user->drawlgEmin;
264: PetscDS ds;
265: Vec phi;
267: PetscCall(DMGetNamedGlobalVector(user->dmPot, "phi", &phi));
268: PetscCall(DMGetDS(user->dmPot, &ds));
269: PetscCall(PetscDSSetObjective(ds, 0, &f0_grad_phi2));
270: PetscCall(DMPlexComputeIntegralFEM(user->dmPot, phi, &intESq, user));
271: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "phi", &phi));
273: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
274: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
275: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
276: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
277: PetscCall(PetscDrawLGAddPoint(user->drawlgE, &t, &lgEmax));
278: PetscCall(PetscDrawLGDraw(user->drawlgE));
279: PetscDraw draw;
280: PetscCall(PetscDrawLGGetDraw(user->drawlgE, &draw));
281: PetscCall(PetscDrawSave(draw));
283: PetscCall(DMSwarmComputeMoments(sw, "velocity", "w_q", pmoments));
284: PetscCall(PetscPrintf(comm, "E: %f\t%+e\t%e\t%f\t%20.15e\t%f\t%f\t%f\t%20.15e\t%20.15e\t%20.15e\t%" PetscInt_FMT "\t(%" PetscInt_FMT ")\n", (double)t, (double)sum, (double)Enorm, (double)lgEnorm, (double)Emax, (double)lgEmax, (double)chargesum, (double)pmoments[0], (double)pmoments[1], (double)pmoments[1 + dim], (double)PetscSqrtReal(intESq), gNp, step));
285: PetscCall(DMViewFromOptions(sw, NULL, "-sw_efield_view"));
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: static PetscErrorCode MonitorMoments(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
290: {
291: DM sw;
292: PetscReal pmoments[16]; /* \int f, \int v f, \int v^2 f */
294: PetscFunctionBeginUser;
295: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS);
296: PetscCall(TSGetDM(ts, &sw));
298: PetscCall(DMSwarmComputeMoments(sw, "velocity", "w_q", pmoments));
300: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%f\t%f\t%f\t%f\n", (double)t, (double)pmoments[0], (double)pmoments[1], (double)pmoments[3]));
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: static PetscErrorCode MonitorParticles(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
305: {
306: AppCtx *user = (AppCtx *)ctx;
307: DM sw;
308: PetscInt Np, dim;
309: const PetscReal *x, *v, *E;
310: const PetscScalar *w;
311: PetscViewer viewer;
312: char filename[64];
313: MPI_Comm comm;
315: PetscFunctionBeginUser;
316: if (step < 0 || step % user->ostep != 0) PetscFunctionReturn(PETSC_SUCCESS);
317: PetscCall(TSGetDM(ts, &sw));
318: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
319: PetscCall(DMGetDimension(sw, &dim));
320: PetscCall(DMSwarmGetLocalSize(sw, &Np));
322: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
323: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
324: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&w));
325: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
327: if (user->particle_monitor_prefix[0]) {
328: PetscCall(PetscSNPrintf(filename, sizeof(filename), "%s_step_%d.bin", user->particle_monitor_prefix, (int)step));
329: } else {
330: PetscCall(PetscSNPrintf(filename, sizeof(filename), "particles_step_%d.bin", (int)step));
331: }
332: PetscCall(PetscViewerBinaryOpen(comm, filename, FILE_MODE_WRITE, &viewer));
334: {
335: Vec vx, vv, vw, vE;
336: PetscCall(VecCreateMPIWithArray(comm, dim, Np * dim, PETSC_DECIDE, x, &vx));
337: PetscCall(VecCreateMPIWithArray(comm, dim, Np * dim, PETSC_DECIDE, v, &vv));
338: PetscCall(VecCreateMPIWithArray(comm, 1, Np, PETSC_DECIDE, w, &vw));
339: PetscCall(VecCreateMPIWithArray(comm, dim, Np * dim, PETSC_DECIDE, E, &vE));
341: PetscCall(VecView(vx, viewer));
342: PetscCall(VecView(vv, viewer));
343: PetscCall(VecView(vw, viewer));
344: PetscCall(VecView(vE, viewer));
346: PetscCall(VecDestroy(&vx));
347: PetscCall(VecDestroy(&vv));
348: PetscCall(VecDestroy(&vw));
349: PetscCall(VecDestroy(&vE));
350: }
351: PetscCall(PetscViewerDestroy(&viewer));
353: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
354: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
355: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
356: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: PetscErrorCode PetscPDFPertubedConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
361: {
362: p[0] = (1. + 0.01 * PetscCosReal(0.5 * x[0])) / (4. * PETSC_PI * 4. * PETSC_PI);
363: return PETSC_SUCCESS;
364: }
365: PetscErrorCode PetscPDFPertubedConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
366: {
367: p[0] = (1. + 0.01 * PetscCosReal(0.5 * x[0])) / (2 * PETSC_PI);
368: return PETSC_SUCCESS;
369: }
371: PetscErrorCode PetscPDFCosine1D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
372: {
373: const PetscReal alpha = scale ? scale[0] : 0.0;
374: const PetscReal k = scale ? scale[1] : 1.;
375: p[0] = (1 + alpha * PetscCosReal(k * x[0]));
376: return PETSC_SUCCESS;
377: }
379: PetscErrorCode PetscPDFCosine2D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
380: {
381: const PetscReal alpha = scale ? scale[0] : 0.;
382: const PetscReal k = scale ? scale[1] : 1.;
383: p[0] = (1 + alpha * PetscCosReal(k * x[0]));
384: return PETSC_SUCCESS;
385: }
387: PetscErrorCode PetscPDFCosine3D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
388: {
389: const PetscReal alpha = scale ? scale[0] : 0.;
390: const PetscReal k = scale ? scale[1] : 1.;
391: p[0] = (1 + alpha * PetscCosReal(k * x[0]));
392: return PETSC_SUCCESS;
393: }
395: static PetscErrorCode CreateVelocityDM(DM sw, DM *vdm)
396: {
397: PetscFE fe;
398: DMPolytopeType ct;
399: PetscInt dim, cStart;
400: const char *prefix = "v";
401: AppCtx *user;
403: PetscFunctionBegin;
404: PetscCall(DMGetDimension(sw, &dim));
405: PetscCall(DMGetApplicationContext(sw, &user));
406: if (dim == 3 && user->use_landau_velocity_space) {
407: LandauCtx *ctx;
408: Vec X;
409: Mat J;
411: PetscCall(DMPlexLandauCreateVelocitySpace(PETSC_COMM_SELF, dim, prefix, &X, &J, &user->landau_pack));
412: PetscCall(DMGetApplicationContext(user->landau_pack, &ctx));
413: *vdm = ctx->plex[0];
414: PetscCall(PetscObjectReference((PetscObject)*vdm));
415: PetscCall(VecDestroy(&X));
416: PetscCall(PetscObjectSetName((PetscObject)*vdm, "velocity"));
417: } else {
418: PetscCall(DMCreate(PETSC_COMM_SELF, vdm));
419: PetscCall(DMSetType(*vdm, DMPLEX));
420: PetscCall(DMPlexSetOptionsPrefix(*vdm, prefix));
421: PetscCall(DMSetFromOptions(*vdm));
422: PetscCall(PetscObjectSetName((PetscObject)*vdm, "velocity"));
423: }
424: PetscCall(DMViewFromOptions(*vdm, NULL, "-dm_view"));
426: PetscCall(DMGetDimension(*vdm, &dim));
427: PetscCall(DMPlexGetHeightStratum(*vdm, 0, &cStart, NULL));
428: PetscCall(DMPlexGetCellType(*vdm, cStart, &ct));
429: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, prefix, PETSC_DETERMINE, &fe));
430: PetscCall(PetscObjectSetName((PetscObject)fe, "distribution"));
431: PetscCall(DMSetField(*vdm, 0, NULL, (PetscObject)fe));
432: PetscCall(DMCreateDS(*vdm));
433: PetscCall(PetscFEDestroy(&fe));
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: /*
438: InitializeParticles_Centroid - Initialize a regular grid of particles.
440: Input Parameters:
441: + sw - The `DMSWARM`
442: - force1D - Treat the spatial domain as 1D
444: Notes:
445: This functions sets the species, cellid, spatial coordinate, and velocity fields for all particles.
447: It places one particle in the centroid of each cell in the implicit tensor product of the spatial
448: and velocity meshes.
449: */
450: static PetscErrorCode InitializeParticles_Centroid(DM sw)
451: {
452: DM_Swarm *swarm = (DM_Swarm *)sw->data;
453: DMSwarmCellDM celldm;
454: DM xdm, vdm;
455: PetscReal vmin[3], vmax[3];
456: PetscReal *x, *v;
457: PetscInt *species, *cellid;
458: PetscInt dim, xcStart, xcEnd, vcStart, vcEnd, Ns, Np, Npc, debug;
459: PetscBool flg;
460: MPI_Comm comm;
461: const char *cellidname;
463: PetscFunctionBegin;
464: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
466: PetscOptionsBegin(comm, "", "DMSwarm Options", "DMSWARM");
467: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
468: PetscCall(PetscOptionsInt("-dm_swarm_num_species", "The number of species", "DMSwarmSetNumSpecies", Ns, &Ns, &flg));
469: if (flg) PetscCall(DMSwarmSetNumSpecies(sw, Ns));
470: PetscCall(PetscOptionsBoundedInt("-dm_swarm_print_coords", "Debug output level for particle coordinate computations", "InitializeParticles", 0, &swarm->printCoords, NULL, 0));
471: PetscCall(PetscOptionsBoundedInt("-dm_swarm_print_weights", "Debug output level for particle weight computations", "InitializeWeights", 0, &swarm->printWeights, NULL, 0));
472: PetscOptionsEnd();
473: debug = swarm->printCoords;
475: PetscCall(DMGetDimension(sw, &dim));
476: PetscCall(DMSwarmGetCellDM(sw, &xdm));
477: PetscCall(DMPlexGetHeightStratum(xdm, 0, &xcStart, &xcEnd));
479: PetscCall(DMSwarmGetCellDMByName(sw, "velocity", &celldm));
480: PetscCall(DMSwarmCellDMGetDM(celldm, &vdm));
481: PetscCall(DMPlexGetHeightStratum(vdm, 0, &vcStart, &vcEnd));
483: // One particle per centroid on the tensor product grid
484: Npc = (vcEnd - vcStart) * Ns;
485: Np = (xcEnd - xcStart) * Npc;
486: PetscCall(DMSwarmSetLocalSizes(sw, Np, 0));
487: if (debug) {
488: PetscInt gNp, gNc = xcEnd - xcStart;
489: PetscCallMPI(MPIU_Allreduce(&Np, &gNp, 1, MPIU_INT, MPIU_SUM, comm));
490: PetscCall(PetscPrintf(comm, "Global Np = %" PetscInt_FMT "\n", gNp));
491: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &gNc, 1, MPIU_INT, MPIU_SUM, comm));
492: PetscCall(PetscPrintf(comm, "Global X-cells = %" PetscInt_FMT "\n", gNc));
493: PetscCall(PetscPrintf(comm, "Global V-cells = %" PetscInt_FMT "\n", vcEnd - vcStart));
494: }
496: // Set species and cellid
497: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
498: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellidname));
499: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
500: PetscCall(DMSwarmGetField(sw, cellidname, NULL, NULL, (void **)&cellid));
501: for (PetscInt c = 0, p = 0; c < xcEnd - xcStart; ++c) {
502: for (PetscInt s = 0; s < Ns; ++s) {
503: for (PetscInt q = 0; q < Npc / Ns; ++q, ++p) {
504: species[p] = s;
505: cellid[p] = c;
506: }
507: }
508: }
509: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
510: PetscCall(DMSwarmRestoreField(sw, cellidname, NULL, NULL, (void **)&cellid));
512: // Set particle coordinates
513: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
514: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
515: PetscCall(DMSwarmSortGetAccess(sw));
516: PetscCall(DMGetBoundingBox(vdm, vmin, vmax));
517: PetscCall(DMGetCoordinatesLocalSetUp(xdm));
518: PetscCall(DMGetCoordinatesLocalSetUp(vdm));
519: for (PetscInt c = 0; c < xcEnd - xcStart; ++c) {
520: const PetscInt cell = c + xcStart;
521: PetscInt *pidx, Npc;
522: PetscReal centroid[3], volume;
524: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
525: PetscCall(DMPlexComputeCellGeometryFVM(xdm, cell, &volume, centroid, NULL));
526: for (PetscInt s = 0; s < Ns; ++s) {
527: for (PetscInt q = 0; q < Npc / Ns; ++q) {
528: const PetscInt p = pidx[q * Ns + s];
529: const PetscInt vc = vcStart + q;
530: PetscReal vcentroid[3], vvolume;
532: PetscCall(DMPlexComputeCellGeometryFVM(vdm, vc, &vvolume, vcentroid, NULL));
533: for (PetscInt d = 0; d < dim; ++d) {
534: x[p * dim + d] = centroid[d];
535: v[p * dim + d] = vcentroid[d];
536: }
538: if (debug > 1) {
539: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Particle %4" PetscInt_FMT " ", p));
540: PetscCall(PetscPrintf(PETSC_COMM_SELF, " x: ("));
541: for (PetscInt d = 0; d < dim; ++d) {
542: if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
543: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", x[p * dim + d]));
544: }
545: PetscCall(PetscPrintf(PETSC_COMM_SELF, ") v:("));
546: for (PetscInt d = 0; d < dim; ++d) {
547: if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
548: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", v[p * dim + d]));
549: }
550: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
551: }
552: }
553: }
554: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
555: }
556: PetscCall(DMSwarmSortRestoreAccess(sw));
557: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
558: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
562: /*
563: InitializeWeights - Compute weight for each local particle
565: Input Parameters:
566: + sw - The `DMSwarm`
567: . totalWeight - The sum of all particle weights
568: . func - The PDF for the particle spatial distribution
569: - param - The PDF parameters
571: Notes:
572: The PDF for velocity is assumed to be a Gaussian
574: The particle weights are returned in the `w_q` field of `sw`.
575: */
576: static PetscErrorCode InitializeWeights(DM sw, PetscReal totalWeight, PetscProbFn *func, const PetscReal param[])
577: {
578: DM xdm, vdm;
579: DMSwarmCellDM celldm;
580: PetscScalar *weight;
581: PetscQuadrature xquad;
582: const PetscReal *xq, *xwq;
583: const PetscInt order = 5;
584: PetscReal xi0[3];
585: PetscReal xwtot = 0., pwtot = 0.;
586: PetscInt xNq;
587: PetscInt dim, Ns, xcStart, xcEnd, vcStart, vcEnd, debug = ((DM_Swarm *)sw->data)->printWeights;
588: MPI_Comm comm;
589: PetscMPIInt rank;
591: PetscFunctionBegin;
592: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
593: PetscCallMPI(MPI_Comm_rank(comm, &rank));
594: PetscCall(DMGetDimension(sw, &dim));
595: PetscCall(DMSwarmGetCellDM(sw, &xdm));
596: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
597: PetscCall(DMPlexGetHeightStratum(xdm, 0, &xcStart, &xcEnd));
598: PetscCall(DMSwarmGetCellDMByName(sw, "velocity", &celldm));
599: PetscCall(DMSwarmCellDMGetDM(celldm, &vdm));
600: PetscCall(DMPlexGetHeightStratum(vdm, 0, &vcStart, &vcEnd));
602: // Setup Quadrature for spatial and velocity weight calculations
603: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, order, -1.0, 1.0, &xquad));
604: PetscCall(PetscQuadratureGetData(xquad, NULL, NULL, &xNq, &xq, &xwq));
605: for (PetscInt d = 0; d < dim; ++d) xi0[d] = -1.0;
607: // Integrate the density function to get the weights of particles in each cell
608: PetscCall(DMGetCoordinatesLocalSetUp(vdm));
609: PetscCall(DMSwarmSortGetAccess(sw));
610: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
611: for (PetscInt c = xcStart; c < xcEnd; ++c) {
612: PetscReal xv0[3], xJ[9], xinvJ[9], xdetJ, xqr[3], xden, xw = 0.;
613: PetscInt *pidx, Npc;
614: PetscInt xNc;
615: const PetscScalar *xarray;
616: PetscScalar *xcoords = NULL;
617: PetscBool xisDG;
619: PetscCall(DMPlexGetCellCoordinates(xdm, c, &xisDG, &xNc, &xarray, &xcoords));
620: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
621: PetscCheck(Npc == (vcEnd - vcStart) * Ns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of particles %" PetscInt_FMT " in cell (rank %d) != %" PetscInt_FMT " number of velocity vertices", Npc, rank, (vcEnd - vcStart) * Ns);
622: PetscCall(DMPlexComputeCellGeometryFEM(xdm, c, NULL, xv0, xJ, xinvJ, &xdetJ));
623: for (PetscInt q = 0; q < xNq; ++q) {
624: // Transform quadrature points from ref space to real space
625: CoordinatesRefToReal(dim, dim, xi0, xv0, xJ, &xq[q * dim], xqr);
626: // Get probability density at quad point
627: // No need to scale xqr since PDF will be periodic
628: PetscCall((*func)(xqr, param, &xden));
629: xw += xden * (xwq[q] * xdetJ);
630: }
631: xwtot += xw;
632: if (debug) {
633: IS globalOrdering;
634: const PetscInt *ordering;
636: PetscCall(DMPlexGetCellNumbering(xdm, &globalOrdering));
637: PetscCall(ISGetIndices(globalOrdering, &ordering));
638: PetscCall(PetscSynchronizedPrintf(comm, "c:%" PetscInt_FMT " [x_a,x_b] = %1.15f,%1.15f -> cell weight = %1.15f\n", ordering[c], (double)PetscRealPart(xcoords[0]), (double)PetscRealPart(xcoords[0 + dim]), (double)xw));
639: PetscCall(ISRestoreIndices(globalOrdering, &ordering));
640: }
641: // Set weights to be Gaussian in velocity cells
642: for (PetscInt vc = vcStart; vc < vcEnd; ++vc) {
643: const PetscInt p = pidx[vc * Ns + 0];
644: PetscReal vw = 0.;
645: PetscInt vNc;
646: const PetscScalar *varray;
647: PetscScalar *vcoords = NULL;
648: PetscBool visDG;
650: PetscCall(DMPlexGetCellCoordinates(vdm, vc, &visDG, &vNc, &varray, &vcoords));
651: PetscCheck(vNc > 0 && vNc % dim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Velocity cell %" PetscInt_FMT " has invalid coordinates (vNc=%" PetscInt_FMT ", dim=%" PetscInt_FMT ")", vc, vNc, dim);
652: {
653: PetscReal vmin[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
654: PetscReal vmax[3] = {-PETSC_MAX_REAL, -PETSC_MAX_REAL, -PETSC_MAX_REAL};
655: PetscInt numVert = vNc / dim;
656: for (PetscInt i = 0; i < numVert; ++i) {
657: for (PetscInt d = 0; d < dim; ++d) {
658: PetscReal coord = PetscRealPart(vcoords[i * dim + d]);
659: vmin[d] = PetscMin(vmin[d], coord);
660: vmax[d] = PetscMax(vmax[d], coord);
661: }
662: }
663: vw = 1.0;
664: for (PetscInt d = 0; d < dim; ++d) vw *= 0.5 * (PetscErfReal(vmax[d] / PetscSqrtReal(2.)) - PetscErfReal(vmin[d] / PetscSqrtReal(2.)));
665: PetscCheck(PetscIsNormalReal(vw), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Particle %" PetscInt_FMT " velocity weight is not normal: vw=%g, vmin=(%g,%g,%g), vmax=(%g,%g,%g)", p, vw, vmin[0], vmin[1], vmin[2], vmax[0], vmax[1], vmax[2]);
666: }
668: weight[p] = totalWeight * vw * xw;
669: pwtot += weight[p];
670: PetscCheck(weight[p] <= 10., PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Particle %" PetscInt_FMT " weight exceeded 10: weight=%g, xw=%g, vw=%g, totalWeight=%g", p, weight[p], xw, vw, totalWeight);
671: PetscCall(DMPlexRestoreCellCoordinates(vdm, vc, &visDG, &vNc, &varray, &vcoords));
672: if (debug > 1) PetscCall(PetscPrintf(comm, "particle %" PetscInt_FMT ": %g, vw: %g xw: %g\n", p, weight[p], vw, xw));
673: }
674: PetscCall(DMPlexRestoreCellCoordinates(xdm, c, &xisDG, &xNc, &xarray, &xcoords));
675: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
676: }
677: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
678: PetscCall(DMSwarmSortRestoreAccess(sw));
679: PetscCall(PetscQuadratureDestroy(&xquad));
681: if (debug) {
682: PetscReal wtot[2] = {pwtot, xwtot};
684: PetscCall(PetscSynchronizedFlush(comm, NULL));
685: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, wtot, 2, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
686: PetscCall(PetscPrintf(comm, "particle weight sum = %1.10f cell weight sum = %1.10f\n", (double)wtot[0], (double)wtot[1]));
687: }
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: static PetscErrorCode InitializeParticles_PerturbedWeights(DM sw, AppCtx *user)
692: {
693: PetscReal scale[2] = {user->cosine_coefficients[0], user->cosine_coefficients[1]};
694: PetscInt dim;
696: PetscFunctionBegin;
697: PetscCall(DMGetDimension(sw, &dim));
698: PetscCall(InitializeParticles_Centroid(sw));
699: PetscCall(InitializeWeights(sw, user->totalWeight, dim == 1 ? PetscPDFCosine1D : (dim == 2 ? PetscPDFCosine2D : PetscPDFCosine3D), scale));
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: static PetscErrorCode InitializeConstants(DM sw, AppCtx *user)
704: {
705: DM dm;
706: PetscInt *species;
707: PetscReal *weight, totalCharge = 0., totalWeight = 0., gmin[3], gmax[3];
708: PetscInt Np, dim;
710: PetscFunctionBegin;
711: PetscCall(DMSwarmGetCellDM(sw, &dm));
712: PetscCall(DMGetDimension(sw, &dim));
713: PetscCall(DMSwarmGetLocalSize(sw, &Np));
714: PetscCall(DMGetBoundingBox(dm, gmin, gmax));
715: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
716: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
717: for (PetscInt p = 0; p < Np; ++p) {
718: totalWeight += weight[p];
719: totalCharge += user->charges[species[p]] * weight[p];
720: }
721: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
722: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
723: {
724: Parameter *param;
725: PetscReal Area;
727: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
728: switch (dim) {
729: case 1:
730: Area = (gmax[0] - gmin[0]);
731: break;
732: case 2:
733: Area = (gmax[0] - gmin[0]) * (gmax[1] - gmin[1]);
734: break;
735: case 3:
736: Area = (gmax[0] - gmin[0]) * (gmax[1] - gmin[1]) * (gmax[2] - gmin[2]);
737: break;
738: default:
739: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", dim);
740: }
741: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &totalWeight, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
742: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &totalCharge, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
743: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "dim = %" PetscInt_FMT "\ttotalWeight = %f, user->charges[species[0]] = %f\ttotalCharge = %f, Total Area = %f\n", dim, (double)totalWeight, (double)user->charges[0], (double)totalCharge, (double)Area));
744: param->sigma = PetscAbsReal(totalCharge / (Area));
746: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sigma: %g\n", (double)param->sigma));
747: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(x0,v0,t0,m0,q0,phi0): (%e, %e, %e, %e, %e, %e) - (P, V) = (%e, %e)\n", (double)param->x0, (double)param->v0, (double)param->t0, (double)param->m0, (double)param->q0, (double)param->phi0, (double)param->poissonNumber,
748: (double)param->vlasovNumber));
749: }
750: /* Setup Constants */
751: {
752: PetscDS ds;
753: Parameter *param;
754: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
755: PetscScalar constants[NUM_CONSTANTS];
756: constants[SIGMA] = param->sigma;
757: constants[V0] = param->v0;
758: constants[T0] = param->t0;
759: constants[X0] = param->x0;
760: constants[M0] = param->m0;
761: constants[Q0] = param->q0;
762: constants[PHI0] = param->phi0;
763: constants[POISSON] = param->poissonNumber;
764: constants[VLASOV] = param->vlasovNumber;
765: PetscCall(DMGetDS(dm, &ds));
766: PetscCall(PetscDSSetConstants(ds, NUM_CONSTANTS, constants));
767: }
768: PetscFunctionReturn(PETSC_SUCCESS);
769: }
771: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
772: {
773: PetscBag bag;
774: Parameter *p;
776: PetscFunctionBeginUser;
777: /* setup PETSc parameter bag */
778: PetscCall(PetscBagGetData(ctx->bag, (void **)&p));
779: PetscCall(PetscBagSetName(ctx->bag, "par", "Vlasov-Poisson Parameters"));
780: bag = ctx->bag;
781: PetscCall(PetscBagRegisterScalar(bag, &p->v0, 1.0, "v0", "Velocity scale, m/s"));
782: PetscCall(PetscBagRegisterScalar(bag, &p->t0, 1.0, "t0", "Time scale, s"));
783: PetscCall(PetscBagRegisterScalar(bag, &p->x0, 1.0, "x0", "Space scale, m"));
784: PetscCall(PetscBagRegisterScalar(bag, &p->phi0, 1.0, "phi0", "Potential scale, kg*m^2/A*s^3"));
785: PetscCall(PetscBagRegisterScalar(bag, &p->q0, 1.0, "q0", "Charge Scale, A*s"));
786: PetscCall(PetscBagRegisterScalar(bag, &p->m0, 1.0, "m0", "Mass Scale, kg"));
787: PetscCall(PetscBagRegisterScalar(bag, &p->epsi0, 1.0, "epsi0", "Permittivity of Free Space, kg"));
788: PetscCall(PetscBagRegisterScalar(bag, &p->kb, 1.0, "kb", "Boltzmann Constant, m^2 kg/s^2 K^1"));
790: PetscCall(PetscBagRegisterScalar(bag, &p->sigma, 1.0, "sigma", "Charge per unit area, C/m^3"));
791: PetscCall(PetscBagRegisterScalar(bag, &p->poissonNumber, 1.0, "poissonNumber", "Non-Dimensional Poisson Number"));
792: PetscCall(PetscBagRegisterScalar(bag, &p->vlasovNumber, 1.0, "vlasovNumber", "Non-Dimensional Vlasov Number"));
793: PetscCall(PetscBagSetFromOptions(bag));
794: {
795: PetscViewer viewer;
796: PetscViewerFormat format;
797: PetscBool flg;
799: PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg));
800: if (flg) {
801: PetscCall(PetscViewerPushFormat(viewer, format));
802: PetscCall(PetscBagView(bag, viewer));
803: PetscCall(PetscViewerFlush(viewer));
804: PetscCall(PetscViewerPopFormat(viewer));
805: PetscCall(PetscViewerDestroy(&viewer));
806: }
807: }
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
812: {
813: const char *prefix = "x";
815: PetscFunctionBeginUser;
816: PetscCall(DMCreate(comm, dm));
817: PetscCall(DMSetType(*dm, DMPLEX));
818: PetscCall(DMPlexSetOptionsPrefix(*dm, prefix));
819: PetscCall(DMSetFromOptions(*dm));
820: PetscCall(PetscObjectSetName((PetscObject)*dm, "space"));
821: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
823: // Cache the mesh geometry
824: DMField coordField;
825: IS cellIS;
826: PetscQuadrature quad;
827: PetscReal *wt, *pt;
828: PetscInt cdim, cStart, cEnd;
830: PetscCall(DMGetCoordinateField(*dm, &coordField));
831: PetscCheck(coordField, comm, PETSC_ERR_USER, "DM must have a coordinate field");
832: PetscCall(DMGetCoordinateDim(*dm, &cdim));
833: PetscCall(DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd));
834: PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cStart, cStart, 1, &cellIS));
835: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
836: PetscCall(PetscMalloc1(1, &wt));
837: PetscCall(PetscMalloc1(cdim, &pt));
838: wt[0] = 1.;
839: for (PetscInt d = 0; d < cdim; ++d) pt[d] = -1.;
840: PetscCall(PetscQuadratureSetData(quad, cdim, 1, 1, pt, wt));
841: PetscCall(DMFieldCreateFEGeom(coordField, cellIS, quad, PETSC_FEGEOM_BASIC, &user->fegeom));
842: PetscCall(PetscQuadratureDestroy(&quad));
843: PetscCall(ISDestroy(&cellIS));
844: PetscFunctionReturn(PETSC_SUCCESS);
845: }
847: static void ion_f0(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[])
848: {
849: f0[0] = -constants[SIGMA];
850: }
852: 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[])
853: {
854: PetscInt d;
855: for (d = 0; d < dim; ++d) f1[d] = u_x[d];
856: }
858: 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[])
859: {
860: for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
861: }
863: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
864: {
865: PetscFE fephi;
866: PetscDS ds;
867: PetscBool simplex;
868: PetscInt dim;
869: MatNullSpace nullsp;
871: PetscFunctionBeginUser;
872: PetscCall(DMGetDimension(dm, &dim));
873: PetscCall(DMPlexIsSimplex(dm, &simplex));
875: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &fephi));
876: PetscCall(PetscObjectSetName((PetscObject)fephi, "potential"));
877: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fephi));
878: PetscCall(DMCreateDS(dm));
879: PetscCall(DMGetDS(dm, &ds));
880: PetscCall(PetscDSSetResidual(ds, 0, ion_f0, laplacian_f1));
881: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, laplacian_g3));
882: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_TRUE, 0, NULL, &nullsp));
883: PetscCall(PetscObjectCompose((PetscObject)fephi, "nullspace", (PetscObject)nullsp));
884: PetscCall(MatNullSpaceDestroy(&nullsp));
885: PetscCall(PetscFEDestroy(&fephi));
886: PetscFunctionReturn(PETSC_SUCCESS);
887: }
889: static PetscErrorCode CreatePoisson(DM dm, AppCtx *user)
890: {
891: SNES snes;
892: Mat J;
893: MatNullSpace nullSpace;
895: PetscFunctionBeginUser;
896: PetscCall(CreateFEM(dm, user));
897: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
898: PetscCall(SNESSetOptionsPrefix(snes, "em_"));
899: PetscCall(SNESSetDM(snes, dm));
900: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, user));
901: PetscCall(SNESSetFromOptions(snes));
903: PetscCall(DMCreateMatrix(dm, &J));
904: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_TRUE, 0, NULL, &nullSpace));
905: PetscCall(MatSetNullSpace(J, nullSpace));
906: PetscCall(MatNullSpaceDestroy(&nullSpace));
907: PetscCall(SNESSetJacobian(snes, J, J, NULL, NULL));
908: PetscCall(MatDestroy(&J));
910: user->dmPot = dm;
911: PetscCall(PetscObjectReference((PetscObject)user->dmPot));
913: PetscCall(DMCreateMassMatrix(user->dmPot, user->dmPot, &user->M));
914: PetscCall(DMPlexCreateClosureIndex(dm, NULL));
915: user->snes = snes;
916: PetscFunctionReturn(PETSC_SUCCESS);
917: }
919: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
920: {
921: DMSwarmCellDM celldm;
922: PetscInt dim;
924: PetscFunctionBeginUser;
925: PetscCall(DMGetDimension(dm, &dim));
926: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
927: PetscCall(DMSetType(*sw, DMSWARM));
928: PetscCall(DMSetDimension(*sw, dim));
929: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
930: PetscCall(DMSetApplicationContext(*sw, user));
932: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
933: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
934: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT));
935: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "E_field", dim, PETSC_REAL));
937: const char *fieldnames[2] = {DMSwarmPICField_coor, "velocity"};
939: PetscCall(DMSwarmCellDMCreate(dm, 2, fieldnames, 1, fieldnames, &celldm));
940: PetscCall(DMSwarmAddCellDM(*sw, celldm));
941: PetscCall(DMSwarmCellDMDestroy(&celldm));
943: const char *vfieldnames[2] = {"w_q"};
944: DM vdm;
946: PetscCall(CreateVelocityDM(*sw, &vdm));
947: PetscCall(DMSwarmCellDMCreate(vdm, 1, vfieldnames, 1, &fieldnames[1], &celldm));
948: PetscCall(DMSwarmAddCellDM(*sw, celldm));
949: PetscCall(DMSwarmCellDMDestroy(&celldm));
950: PetscCall(DMDestroy(&vdm));
952: DM mdm;
954: PetscCall(DMClone(dm, &mdm));
955: PetscCall(PetscObjectSetName((PetscObject)mdm, "moments"));
956: PetscCall(DMCopyDisc(dm, mdm));
957: PetscCall(DMSwarmCellDMCreate(mdm, 1, vfieldnames, 1, fieldnames, &celldm));
958: PetscCall(DMDestroy(&mdm));
959: PetscCall(DMSwarmAddCellDM(*sw, celldm));
960: PetscCall(DMSwarmCellDMDestroy(&celldm));
962: PetscCall(DMSetFromOptions(*sw));
963: PetscCall(DMSetUp(*sw));
965: PetscCall(DMSwarmSetCellDMActive(*sw, "space"));
966: user->swarm = *sw;
967: // TODO: This is redundant init since it is done in InitializeSolveAndSwarm, however DMSetUp() requires the local size be set
968: PetscCall(InitializeParticles_PerturbedWeights(*sw, user));
969: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
970: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
971: PetscFunctionReturn(PETSC_SUCCESS);
972: }
974: static PetscErrorCode ComputeFieldAtParticles_Primal(SNES snes, DM sw, Mat M_p, PetscReal E[])
975: {
976: DM dm;
977: AppCtx *user;
978: PetscDS ds;
979: PetscFE fe;
980: KSP ksp;
981: Vec rhoRhs; // Weak charge density, \int phi_i rho
982: Vec rho; // Charge density, M^{-1} rhoRhs
983: Vec phi, locPhi; // Potential
984: Vec f; // Particle weights
985: PetscReal *coords;
986: PetscInt dim, cStart, cEnd, Np;
988: PetscFunctionBegin;
989: PetscCall(DMGetApplicationContext(sw, (void *)&user));
990: PetscCall(PetscLogEventBegin(user->ESolveEvent, snes, sw, 0, 0));
991: PetscCall(DMGetDimension(sw, &dim));
992: PetscCall(DMSwarmGetLocalSize(sw, &Np));
994: PetscCall(SNESGetDM(snes, &dm));
995: PetscCall(DMGetGlobalVector(dm, &rhoRhs));
996: PetscCall(PetscObjectSetName((PetscObject)rhoRhs, "Weak charge density"));
997: PetscCall(DMGetNamedGlobalVector(user->dmPot, "rho", &rho));
998: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
999: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
1001: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
1002: PetscCall(MatViewFromOptions(user->M, NULL, "-m_view"));
1003: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
1005: PetscCall(MatMultTranspose(M_p, f, rhoRhs));
1006: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
1008: // Low-pass filter rhoRhs
1009: PetscInt window = 0;
1010: PetscCall(PetscOptionsGetInt(NULL, NULL, "-rho_average", &window, NULL));
1011: if (window) {
1012: PetscScalar *a;
1013: PetscInt n;
1014: PetscReal width = 2. * window + 1.;
1016: // This will only work for P_1
1017: // This works because integration against a test function is linear
1018: // Do a real integral against weight function for higher order
1019: PetscCall(VecGetLocalSize(rhoRhs, &n));
1020: PetscCall(VecGetArrayWrite(rhoRhs, &a));
1021: for (PetscInt i = 0; i < n; ++i) {
1022: PetscScalar avg = a[i];
1023: for (PetscInt j = 1; j <= window; ++j) avg += a[(i - j + n) % n] + a[(i + j) % n];
1024: a[i] = avg / width;
1025: //a[i] = (a[(i - 1 + n) % n] + a[i] + a[(i + 1) % n]) / 3.;
1026: }
1027: PetscCall(VecRestoreArrayWrite(rhoRhs, &a));
1028: }
1030: PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &ksp));
1031: PetscCall(KSPSetOptionsPrefix(ksp, "em_proj_"));
1032: PetscCall(KSPSetOperators(ksp, user->M, user->M));
1033: PetscCall(KSPSetFromOptions(ksp));
1034: PetscCall(KSPSolve(ksp, rhoRhs, rho));
1036: PetscCall(VecScale(rhoRhs, -1.0));
1038: PetscCall(VecViewFromOptions(rhoRhs, NULL, "-rho_view"));
1039: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "rho", &rho));
1040: PetscCall(KSPDestroy(&ksp));
1042: PetscCall(DMGetNamedGlobalVector(user->dmPot, "phi", &phi));
1043: PetscCall(VecSet(phi, 0.0));
1044: PetscCall(SNESSolve(snes, rhoRhs, phi));
1045: PetscCall(DMRestoreGlobalVector(dm, &rhoRhs));
1046: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
1048: PetscCall(DMGetLocalVector(dm, &locPhi));
1049: PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
1050: PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
1051: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "phi", &phi));
1052: PetscCall(PetscLogEventEnd(user->ESolveEvent, snes, sw, 0, 0));
1054: PetscCall(DMGetDS(dm, &ds));
1055: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
1056: PetscCall(DMSwarmSortGetAccess(sw));
1057: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1058: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1060: PetscCall(PetscLogEventBegin(user->ETabEvent, snes, sw, 0, 0));
1061: PetscTabulation tab;
1062: PetscReal *pcoord, *refcoord;
1063: PetscFEGeom *chunkgeom = NULL;
1064: PetscInt maxNcp = 0;
1066: for (PetscInt c = cStart; c < cEnd; ++c) {
1067: PetscInt Ncp;
1069: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, c, &Ncp));
1070: maxNcp = PetscMax(maxNcp, Ncp);
1071: }
1072: PetscCall(DMGetWorkArray(dm, maxNcp * dim, MPIU_REAL, &refcoord));
1073: PetscCall(PetscArrayzero(refcoord, maxNcp * dim));
1074: PetscCall(DMGetWorkArray(dm, maxNcp * dim, MPIU_REAL, &pcoord));
1075: PetscCall(PetscFECreateTabulation(fe, 1, maxNcp, refcoord, 1, &tab));
1076: for (PetscInt c = cStart; c < cEnd; ++c) {
1077: PetscScalar *clPhi = NULL;
1078: PetscInt *points;
1079: PetscInt Ncp;
1081: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
1082: for (PetscInt cp = 0; cp < Ncp; ++cp)
1083: for (PetscInt d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
1084: {
1085: PetscCall(PetscFEGeomGetChunk(user->fegeom, c - cStart, c - cStart + 1, &chunkgeom));
1086: for (PetscInt i = 0; i < Ncp; ++i) {
1087: const PetscReal x0[3] = {-1., -1., -1.};
1088: CoordinatesRealToRef(dim, dim, x0, chunkgeom->v, chunkgeom->invJ, &pcoord[dim * i], &refcoord[dim * i]);
1089: }
1090: }
1091: PetscCall(PetscFEComputeTabulation(fe, Ncp, refcoord, 1, tab));
1092: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
1093: for (PetscInt cp = 0; cp < Ncp; ++cp) {
1094: const PetscReal *basisDer = tab->T[1];
1095: const PetscInt p = points[cp];
1097: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] = 0.;
1098: PetscCall(PetscFEFreeInterpolateGradient_Static(fe, basisDer, clPhi, dim, chunkgeom->invJ, NULL, cp, &E[p * dim]));
1099: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] *= -1.0;
1100: }
1101: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
1102: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
1103: }
1104: PetscCall(DMRestoreWorkArray(dm, maxNcp * dim, MPIU_REAL, &pcoord));
1105: PetscCall(DMRestoreWorkArray(dm, maxNcp * dim, MPIU_REAL, &refcoord));
1106: PetscCall(PetscTabulationDestroy(&tab));
1107: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1108: PetscCall(DMSwarmSortRestoreAccess(sw));
1109: PetscCall(DMRestoreLocalVector(dm, &locPhi));
1110: PetscCall(PetscFEGeomRestoreChunk(user->fegeom, 0, 1, &chunkgeom));
1111: PetscCall(PetscLogEventEnd(user->ETabEvent, snes, sw, 0, 0));
1112: PetscFunctionReturn(PETSC_SUCCESS);
1113: }
1115: static PetscErrorCode ComputeFieldAtParticles(SNES snes, DM sw)
1116: {
1117: AppCtx *user;
1118: Mat M_p;
1119: PetscReal *E;
1120: PetscInt dim, Np;
1122: PetscFunctionBegin;
1125: PetscCall(DMGetDimension(sw, &dim));
1126: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1127: PetscCall(DMGetApplicationContext(sw, &user));
1129: PetscCall(DMSwarmSetCellDMActive(sw, "moments"));
1130: // TODO: Could share sort context with space cellDM
1131: PetscCall(DMSwarmMigrate(sw, PETSC_FALSE));
1132: PetscCall(DMCreateMassMatrix(sw, user->dmPot, &M_p));
1133: PetscCall(DMSwarmSetCellDMActive(sw, "space"));
1135: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
1136: PetscCall(PetscArrayzero(E, Np * dim));
1137: user->validE = PETSC_TRUE;
1139: PetscCall(ComputeFieldAtParticles_Primal(snes, sw, M_p, E));
1140: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
1141: PetscCall(MatDestroy(&M_p));
1142: PetscFunctionReturn(PETSC_SUCCESS);
1143: }
1145: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx)
1146: {
1147: DM sw;
1148: SNES snes = ((AppCtx *)ctx)->snes;
1149: const PetscScalar *u;
1150: PetscScalar *g;
1151: PetscReal *E, m_p = 1., q_p = -1.;
1152: PetscInt dim, d, Np, p;
1154: PetscFunctionBeginUser;
1155: PetscCall(TSGetDM(ts, &sw));
1156: PetscCall(ComputeFieldAtParticles(snes, sw));
1158: PetscCall(DMGetDimension(sw, &dim));
1159: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1160: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
1161: PetscCall(VecGetArrayRead(U, &u));
1162: PetscCall(VecGetArray(G, &g));
1163: Np /= 2 * dim;
1164: for (p = 0; p < Np; ++p) {
1165: for (d = 0; d < dim; ++d) {
1166: g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d];
1167: g[(p * 2 + 1) * dim + d] = q_p * E[p * dim + d] / m_p;
1168: }
1169: }
1170: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
1171: PetscCall(VecRestoreArrayRead(U, &u));
1172: PetscCall(VecRestoreArray(G, &g));
1173: PetscFunctionReturn(PETSC_SUCCESS);
1174: }
1176: /* J_{ij} = dF_i/dx_j
1177: J_p = ( 0 1)
1178: (-w^2 0)
1179: TODO Now there is another term with w^2 from the electric field. I think we will need to invert the operator.
1180: Perhaps we can approximate the Jacobian using only the cellwise P-P gradient from Coulomb
1181: */
1182: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, void *ctx)
1183: {
1184: DM sw;
1185: const PetscReal *coords, *vel;
1186: PetscInt dim, d, Np, p, rStart;
1188: PetscFunctionBeginUser;
1189: PetscCall(TSGetDM(ts, &sw));
1190: PetscCall(DMGetDimension(sw, &dim));
1191: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1192: PetscCall(MatGetOwnershipRange(J, &rStart, NULL));
1193: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1194: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&vel));
1195: Np /= 2 * dim;
1196: for (p = 0; p < Np; ++p) {
1197: // TODO This is not right because dv/dx has the electric field in it
1198: PetscScalar vals[4] = {0., 1., -1, 0.};
1200: for (d = 0; d < dim; ++d) {
1201: const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
1202: PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES));
1203: }
1204: }
1205: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1206: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&vel));
1207: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1208: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1209: PetscFunctionReturn(PETSC_SUCCESS);
1210: }
1212: static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
1213: {
1214: AppCtx *user = (AppCtx *)ctx;
1215: DM sw;
1216: const PetscScalar *v;
1217: PetscScalar *xres;
1218: PetscInt Np, p, d, dim;
1220: PetscFunctionBeginUser;
1221: PetscCall(PetscLogEventBegin(user->RhsXEvent, ts, 0, 0, 0));
1222: PetscCall(TSGetDM(ts, &sw));
1223: PetscCall(DMGetDimension(sw, &dim));
1224: PetscCall(VecGetLocalSize(Xres, &Np));
1225: PetscCall(VecGetArrayRead(V, &v));
1226: PetscCall(VecGetArray(Xres, &xres));
1227: Np /= dim;
1228: for (p = 0; p < Np; ++p) {
1229: for (d = 0; d < dim; ++d) xres[p * dim + d] = v[p * dim + d];
1230: }
1231: PetscCall(VecRestoreArrayRead(V, &v));
1232: PetscCall(VecRestoreArray(Xres, &xres));
1233: PetscCall(PetscLogEventEnd(user->RhsXEvent, ts, 0, 0, 0));
1234: PetscFunctionReturn(PETSC_SUCCESS);
1235: }
1237: static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
1238: {
1239: DM sw;
1240: AppCtx *user = (AppCtx *)ctx;
1241: SNES snes = ((AppCtx *)ctx)->snes;
1242: const PetscScalar *x;
1243: PetscScalar *vres;
1244: PetscReal *E, m_p, q_p;
1245: PetscInt Np, p, dim, d;
1246: Parameter *param;
1248: PetscFunctionBeginUser;
1249: PetscCall(PetscLogEventBegin(user->RhsVEvent, ts, 0, 0, 0));
1250: PetscCall(TSGetDM(ts, &sw));
1251: PetscCall(ComputeFieldAtParticles(snes, sw));
1253: PetscCall(DMGetDimension(sw, &dim));
1254: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
1255: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1256: m_p = user->masses[0] * param->m0;
1257: q_p = user->charges[0] * param->q0;
1258: PetscCall(VecGetLocalSize(Vres, &Np));
1259: PetscCall(VecGetArrayRead(X, &x));
1260: PetscCall(VecGetArray(Vres, &vres));
1261: Np /= dim;
1262: for (p = 0; p < Np; ++p) {
1263: for (d = 0; d < dim; ++d) vres[p * dim + d] = q_p * E[p * dim + d] / m_p;
1264: }
1265: PetscCall(VecRestoreArrayRead(X, &x));
1266: /*
1267: Synchronized, ordered output for parallel/sequential test cases.
1268: In the 1D (on the 2D mesh) case, every y component should be zero.
1269: */
1270: if (user->checkVRes) {
1271: PetscBool pr = user->checkVRes > 1 ? PETSC_TRUE : PETSC_FALSE;
1272: PetscInt step;
1274: PetscCall(TSGetStepNumber(ts, &step));
1275: if (pr) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "step: %" PetscInt_FMT "\n", step));
1276: for (PetscInt p = 0; p < Np; ++p) {
1277: if (pr) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Residual: %.12g %.12g\n", (double)PetscRealPart(vres[p * dim + 0]), (double)PetscRealPart(vres[p * dim + 1])));
1278: PetscCheck(PetscAbsScalar(vres[p * dim + 1]) < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Y velocity should be 0., not %g", (double)PetscRealPart(vres[p * dim + 1]));
1279: }
1280: if (pr) PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
1281: }
1282: PetscCall(VecRestoreArray(Vres, &vres));
1283: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
1284: PetscCall(PetscLogEventEnd(user->RhsVEvent, ts, 0, 0, 0));
1285: PetscFunctionReturn(PETSC_SUCCESS);
1286: }
1288: /* Discrete Gradients Formulation: S, F, gradF (G) */
1289: PetscErrorCode RHSJacobianS(TS ts, PetscReal t, Vec U, Mat S, void *ctx)
1290: {
1291: PetscScalar vals[4] = {0., 1., -1., 0.};
1292: DM sw;
1293: PetscInt dim, d, Np, p, rStart;
1295: PetscFunctionBeginUser;
1296: PetscCall(TSGetDM(ts, &sw));
1297: PetscCall(DMGetDimension(sw, &dim));
1298: PetscCall(VecGetLocalSize(U, &Np));
1299: PetscCall(MatGetOwnershipRange(S, &rStart, NULL));
1300: Np /= 2 * dim;
1301: for (p = 0; p < Np; ++p) {
1302: for (d = 0; d < dim; ++d) {
1303: const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
1304: PetscCall(MatSetValues(S, 2, rows, 2, rows, vals, INSERT_VALUES));
1305: }
1306: }
1307: PetscCall(MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY));
1308: PetscCall(MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY));
1309: PetscFunctionReturn(PETSC_SUCCESS);
1310: }
1312: PetscErrorCode RHSObjectiveF(TS ts, PetscReal t, Vec U, PetscScalar *F, void *ctx)
1313: {
1314: AppCtx *user = (AppCtx *)ctx;
1315: DM sw;
1316: Vec phi;
1317: const PetscScalar *u;
1318: PetscInt dim, Np, cStart, cEnd;
1319: PetscReal *vel, *coords, m_p = 1.;
1321: PetscFunctionBeginUser;
1322: PetscCall(TSGetDM(ts, &sw));
1323: PetscCall(DMGetDimension(sw, &dim));
1324: PetscCall(DMPlexGetHeightStratum(user->dmPot, 0, &cStart, &cEnd));
1326: PetscCall(DMGetNamedGlobalVector(user->dmPot, "phi", &phi));
1327: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view_dg"));
1328: PetscCall(computeFieldEnergy(user->dmPot, phi, F));
1329: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "phi", &phi));
1331: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1332: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&vel));
1333: PetscCall(DMSwarmSortGetAccess(sw));
1334: PetscCall(VecGetArrayRead(U, &u));
1335: PetscCall(VecGetLocalSize(U, &Np));
1336: Np /= 2 * dim;
1337: for (PetscInt c = cStart; c < cEnd; ++c) {
1338: PetscInt *points;
1339: PetscInt Ncp;
1341: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
1342: for (PetscInt cp = 0; cp < Ncp; ++cp) {
1343: const PetscInt p = points[cp];
1344: const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);
1346: *F += 0.5 * m_p * v2;
1347: }
1348: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
1349: }
1350: PetscCall(VecRestoreArrayRead(U, &u));
1351: PetscCall(DMSwarmSortRestoreAccess(sw));
1352: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1353: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&vel));
1354: PetscFunctionReturn(PETSC_SUCCESS);
1355: }
1357: /* dF/dx = q E dF/dv = v */
1358: PetscErrorCode RHSFunctionG(TS ts, PetscReal t, Vec U, Vec G, void *ctx)
1359: {
1360: DM sw;
1361: SNES snes = ((AppCtx *)ctx)->snes;
1362: const PetscReal *coords, *vel, *E;
1363: const PetscScalar *u;
1364: PetscScalar *g;
1365: PetscReal m_p = 1., q_p = -1.;
1366: PetscInt dim, d, Np, p;
1368: PetscFunctionBeginUser;
1369: PetscCall(TSGetDM(ts, &sw));
1370: PetscCall(DMGetDimension(sw, &dim));
1371: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1372: PetscCall(VecGetArrayRead(U, &u));
1373: PetscCall(VecGetArray(G, &g));
1375: PetscLogEvent COMPUTEFIELD;
1376: PetscCall(PetscLogEventRegister("COMPFIELDATPART", TS_CLASSID, &COMPUTEFIELD));
1377: PetscCall(PetscLogEventBegin(COMPUTEFIELD, 0, 0, 0, 0));
1378: PetscCall(ComputeFieldAtParticles(snes, sw));
1379: PetscCall(PetscLogEventEnd(COMPUTEFIELD, 0, 0, 0, 0));
1380: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1381: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&vel));
1382: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
1383: for (p = 0; p < Np; ++p) {
1384: for (d = 0; d < dim; ++d) {
1385: g[(p * 2 + 0) * dim + d] = -(q_p / m_p) * E[p * dim + d];
1386: g[(p * 2 + 1) * dim + d] = u[(p * 2 + 1) * dim + d];
1387: }
1388: }
1389: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
1390: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1391: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&vel));
1392: PetscCall(VecRestoreArrayRead(U, &u));
1393: PetscCall(VecRestoreArray(G, &g));
1394: PetscFunctionReturn(PETSC_SUCCESS);
1395: }
1397: static PetscErrorCode CreateSolution(TS ts)
1398: {
1399: DM sw;
1400: Vec u;
1401: PetscInt dim, Np;
1403: PetscFunctionBegin;
1404: PetscCall(TSGetDM(ts, &sw));
1405: PetscCall(DMGetDimension(sw, &dim));
1406: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1407: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
1408: PetscCall(VecSetBlockSize(u, dim));
1409: PetscCall(VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE));
1410: PetscCall(VecSetUp(u));
1411: PetscCall(TSSetSolution(ts, u));
1412: PetscCall(VecDestroy(&u));
1413: PetscFunctionReturn(PETSC_SUCCESS);
1414: }
1416: static PetscErrorCode SetProblem(TS ts)
1417: {
1418: AppCtx *user;
1419: DM sw;
1421: PetscFunctionBegin;
1422: PetscCall(TSGetDM(ts, &sw));
1423: PetscCall(DMGetApplicationContext(sw, (void **)&user));
1424: // Define unified system for (X, V)
1425: {
1426: Mat J;
1427: PetscInt dim, Np;
1429: PetscCall(DMGetDimension(sw, &dim));
1430: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1431: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
1432: PetscCall(MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE));
1433: PetscCall(MatSetBlockSize(J, 2 * dim));
1434: PetscCall(MatSetFromOptions(J));
1435: PetscCall(MatSetUp(J));
1436: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, user));
1437: PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, user));
1438: PetscCall(MatDestroy(&J));
1439: }
1440: /* Define split system for X and V */
1441: {
1442: Vec u;
1443: IS isx, isv, istmp;
1444: const PetscInt *idx;
1445: PetscInt dim, Np, rstart;
1447: PetscCall(TSGetSolution(ts, &u));
1448: PetscCall(DMGetDimension(sw, &dim));
1449: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1450: PetscCall(VecGetOwnershipRange(u, &rstart, NULL));
1451: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 0, 2, &istmp));
1452: PetscCall(ISGetIndices(istmp, &idx));
1453: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isx));
1454: PetscCall(ISRestoreIndices(istmp, &idx));
1455: PetscCall(ISDestroy(&istmp));
1456: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 1, 2, &istmp));
1457: PetscCall(ISGetIndices(istmp, &idx));
1458: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isv));
1459: PetscCall(ISRestoreIndices(istmp, &idx));
1460: PetscCall(ISDestroy(&istmp));
1461: PetscCall(TSRHSSplitSetIS(ts, "position", isx));
1462: PetscCall(TSRHSSplitSetIS(ts, "momentum", isv));
1463: PetscCall(ISDestroy(&isx));
1464: PetscCall(ISDestroy(&isv));
1465: PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, user));
1466: PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, user));
1467: }
1468: // Define symplectic formulation U_t = S . G, where G = grad F
1469: {
1470: PetscCall(TSDiscGradSetFormulation(ts, RHSJacobianS, RHSObjectiveF, RHSFunctionG, user));
1471: }
1472: PetscFunctionReturn(PETSC_SUCCESS);
1473: }
1475: static PetscErrorCode DMSwarmTSRedistribute(TS ts)
1476: {
1477: DM sw;
1478: Vec u;
1479: PetscReal t, maxt, dt;
1480: PetscInt n, maxn;
1482: PetscFunctionBegin;
1483: PetscCall(TSGetDM(ts, &sw));
1484: PetscCall(TSGetTime(ts, &t));
1485: PetscCall(TSGetMaxTime(ts, &maxt));
1486: PetscCall(TSGetTimeStep(ts, &dt));
1487: PetscCall(TSGetStepNumber(ts, &n));
1488: PetscCall(TSGetMaxSteps(ts, &maxn));
1490: PetscCall(TSReset(ts));
1491: PetscCall(TSSetDM(ts, sw));
1492: PetscCall(TSSetFromOptions(ts));
1493: PetscCall(TSSetTime(ts, t));
1494: PetscCall(TSSetMaxTime(ts, maxt));
1495: PetscCall(TSSetTimeStep(ts, dt));
1496: PetscCall(TSSetStepNumber(ts, n));
1497: PetscCall(TSSetMaxSteps(ts, maxn));
1499: PetscCall(CreateSolution(ts));
1500: PetscCall(SetProblem(ts));
1501: PetscCall(TSGetSolution(ts, &u));
1502: PetscFunctionReturn(PETSC_SUCCESS);
1503: }
1505: PetscErrorCode line(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
1506: {
1507: DM sw, cdm;
1508: PetscInt Np;
1509: PetscReal low[2], high[2];
1510: AppCtx *user = (AppCtx *)ctx;
1512: sw = user->swarm;
1513: PetscCall(DMSwarmGetCellDM(sw, &cdm));
1514: // Get the bounding box so we can equally space the particles
1515: PetscCall(DMGetLocalBoundingBox(cdm, low, high));
1516: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1517: // shift it by h/2 so nothing is initialized directly on a boundary
1518: x[0] = ((high[0] - low[0]) / Np) * (p + 0.5);
1519: x[1] = 0.;
1520: return PETSC_SUCCESS;
1521: }
1523: /*
1524: InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values.
1526: Input Parameters:
1527: + ts - The TS
1528: - useInitial - Flag to also set the initial conditions to the current coordinates and velocities and setup the problem
1530: Output Parameters:
1531: . u - The initialized solution vector
1533: Level: advanced
1535: .seealso: InitializeSolve()
1536: */
1537: static PetscErrorCode InitializeSolveAndSwarm(TS ts, PetscBool useInitial)
1538: {
1539: DM sw;
1540: Vec u, gc, gv;
1541: IS isx, isv;
1542: PetscInt dim;
1543: AppCtx *user;
1545: PetscFunctionBeginUser;
1546: PetscCall(TSGetDM(ts, &sw));
1547: PetscCall(DMGetApplicationContext(sw, &user));
1548: PetscCall(DMGetDimension(sw, &dim));
1549: if (useInitial) {
1550: PetscCall(InitializeParticles_PerturbedWeights(sw, user));
1551: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
1552: PetscCall(DMSwarmTSRedistribute(ts));
1553: }
1554: PetscCall(DMSetUp(sw));
1555: PetscCall(TSGetSolution(ts, &u));
1556: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
1557: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
1558: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
1559: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
1560: PetscCall(VecISCopy(u, isx, SCATTER_FORWARD, gc));
1561: PetscCall(VecISCopy(u, isv, SCATTER_FORWARD, gv));
1562: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
1563: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
1564: PetscFunctionReturn(PETSC_SUCCESS);
1565: }
1567: static PetscErrorCode InitializeSolve(TS ts, Vec u)
1568: {
1569: PetscFunctionBegin;
1570: PetscCall(TSSetSolution(ts, u));
1571: PetscCall(InitializeSolveAndSwarm(ts, PETSC_TRUE));
1572: PetscFunctionReturn(PETSC_SUCCESS);
1573: }
1575: static PetscErrorCode MigrateParticles(TS ts)
1576: {
1577: DM sw, cdm;
1578: const PetscReal *L;
1579: AppCtx *ctx;
1581: PetscFunctionBeginUser;
1582: PetscCall(TSGetDM(ts, &sw));
1583: PetscCall(DMGetApplicationContext(sw, &ctx));
1584: PetscCall(DMViewFromOptions(sw, NULL, "-migrate_view_pre"));
1585: {
1586: Vec u, gc, gv, position, momentum;
1587: IS isx, isv;
1588: PetscReal *pos, *mom;
1590: PetscCall(TSGetSolution(ts, &u));
1591: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
1592: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
1593: PetscCall(VecGetSubVector(u, isx, &position));
1594: PetscCall(VecGetSubVector(u, isv, &momentum));
1595: PetscCall(VecGetArray(position, &pos));
1596: PetscCall(VecGetArray(momentum, &mom));
1597: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
1598: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
1599: PetscCall(VecISCopy(u, isx, SCATTER_REVERSE, gc));
1600: PetscCall(VecISCopy(u, isv, SCATTER_REVERSE, gv));
1602: PetscCall(DMSwarmGetCellDM(sw, &cdm));
1603: PetscCall(DMGetPeriodicity(cdm, NULL, NULL, &L));
1604: if ((L[0] || L[1]) >= 0.) {
1605: PetscReal *x, *v, upper[3], lower[3];
1606: PetscInt Np, dim;
1608: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1609: PetscCall(DMGetDimension(cdm, &dim));
1610: PetscCall(DMGetBoundingBox(cdm, lower, upper));
1611: PetscCall(VecGetArray(gc, &x));
1612: PetscCall(VecGetArray(gv, &v));
1613: for (PetscInt p = 0; p < Np; ++p) {
1614: for (PetscInt d = 0; d < dim; ++d) {
1615: if (pos[p * dim + d] < lower[d]) {
1616: x[p * dim + d] = pos[p * dim + d] + (upper[d] - lower[d]);
1617: } else if (pos[p * dim + d] > upper[d]) {
1618: x[p * dim + d] = pos[p * dim + d] - (upper[d] - lower[d]);
1619: } else {
1620: x[p * dim + d] = pos[p * dim + d];
1621: }
1622: PetscCheck(x[p * dim + d] >= lower[d] && x[p * dim + d] <= upper[d], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "p: %" PetscInt_FMT "x[%" PetscInt_FMT "] %g", p, d, (double)x[p * dim + d]);
1623: v[p * dim + d] = mom[p * dim + d];
1624: }
1625: }
1626: PetscCall(VecRestoreArray(gc, &x));
1627: PetscCall(VecRestoreArray(gv, &v));
1628: }
1629: PetscCall(VecRestoreArray(position, &pos));
1630: PetscCall(VecRestoreArray(momentum, &mom));
1631: PetscCall(VecRestoreSubVector(u, isx, &position));
1632: PetscCall(VecRestoreSubVector(u, isv, &momentum));
1633: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
1634: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
1635: }
1636: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
1637: // This MUST come last, since it recreates the subswarms and they must DMClone() the new swarm
1638: PetscCall(DMSwarmTSRedistribute(ts));
1639: PetscCall(InitializeSolveAndSwarm(ts, PETSC_FALSE));
1640: {
1641: const char *fieldnames[2] = {DMSwarmPICField_coor, "velocity"};
1642: PetscCall(DMSwarmVectorDefineFields(sw, 2, fieldnames));
1643: }
1644: PetscFunctionReturn(PETSC_SUCCESS);
1645: }
1647: int main(int argc, char **argv)
1648: {
1649: DM dm, sw;
1650: TS ts;
1651: Vec u;
1652: PetscReal dt;
1653: PetscInt maxn;
1654: AppCtx user;
1656: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1657: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
1658: PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &user.bag));
1659: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
1660: PetscCall(CreatePoisson(dm, &user));
1661: PetscCall(CreateSwarm(dm, &user, &sw));
1662: PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
1663: PetscCall(InitializeConstants(sw, &user));
1664: PetscCall(DMSetApplicationContext(sw, &user));
1666: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1667: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
1668: PetscCall(TSSetDM(ts, sw));
1669: PetscCall(TSSetMaxTime(ts, 0.1));
1670: PetscCall(TSSetTimeStep(ts, 0.00001));
1671: PetscCall(TSSetMaxSteps(ts, 100));
1672: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1674: if (user.efield_monitor) PetscCall(TSMonitorSet(ts, MonitorEField, &user, NULL));
1675: if (user.moment_monitor) PetscCall(TSMonitorSet(ts, MonitorMoments, &user, NULL));
1676: if (user.particle_monitor) PetscCall(TSMonitorSet(ts, MonitorParticles, &user, NULL));
1678: PetscCall(TSSetFromOptions(ts));
1679: PetscCall(TSGetTimeStep(ts, &dt));
1680: PetscCall(TSGetMaxSteps(ts, &maxn));
1681: user.steps = maxn;
1682: user.stepSize = dt;
1683: PetscCall(SetupContext(dm, sw, &user));
1684: PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
1685: PetscCall(TSSetPostStep(ts, MigrateParticles));
1686: PetscCall(CreateSolution(ts));
1687: PetscCall(TSGetSolution(ts, &u));
1688: PetscCall(TSComputeInitialCondition(ts, u));
1689: PetscCall(CheckNonNegativeWeights(sw, &user));
1690: PetscCall(TSSolve(ts, NULL));
1692: PetscCall(SNESDestroy(&user.snes));
1693: PetscCall(DMDestroy(&user.dmPot));
1694: PetscCall(MatDestroy(&user.M));
1695: PetscCall(PetscFEGeomDestroy(&user.fegeom));
1696: PetscCall(TSDestroy(&ts));
1697: PetscCall(DMDestroy(&sw));
1698: PetscCall(DMDestroy(&dm));
1699: PetscCall(DestroyContext(&user));
1700: PetscCall(PetscFinalize());
1701: return 0;
1702: }
1704: /*TEST
1706: build:
1707: requires: !complex double
1709: testset:
1710: nsize: 2
1711: args: -cosine_coefficients 0.01 -charges -1. -total_weight 1. -xdm_plex_hash_location -vpetscspace_degree 2 -petscspace_degree 1 -em_snes_atol 1.e-12 \
1712: -em_snes_error_if_not_converged -em_ksp_error_if_not_converged -em_ksp_type cg -em_pc_type gamg -em_mg_coarse_ksp_type preonly -em_mg_coarse_pc_type svd -em_proj_ksp_type cg \
1713: -em_proj_pc_type gamg -em_proj_mg_coarse_ksp_type preonly -em_proj_mg_coarse_pc_type svd -ts_time_step 0.03 -xdm_plex_simplex 0 \
1714: -ts_max_time 100 -output_step 1 -ts_type basicsymplectic -ts_basicsymplectic_type 1 -check_weights -ts_max_steps 60
1716: test:
1717: suffix: landau_damping_1d
1718: args: -xdm_plex_dim 1 -xdm_plex_box_faces 60 -xdm_plex_box_lower 0. -xdm_plex_box_upper 12.5664 -xdm_plex_box_bd periodic -vdm_plex_dim 1 -vdm_plex_box_faces 60 \
1719: -vdm_plex_box_lower -6 -vdm_plex_box_upper 6 -vdm_plex_hash_location -vdm_plex_simplex 0 -vdm_plex_box_bd none
1721: test:
1722: suffix: landau_damping_2d
1723: args: -xdm_plex_dim 2 -xdm_plex_box_bd periodic,periodic -vdm_plex_dim 2 -xdm_plex_box_lower 0.,-.5 -vdm_plex_box_lower -6,-6 -vdm_plex_box_upper 6,6 -xdm_plex_box_faces 6,3 \
1724: -xdm_plex_box_upper 12.5664,.5 -vdm_plex_box_faces 15,15 -vdm_plex_box_bd none,none -vdm_plex_hash_location -vdm_plex_simplex 0
1726: test:
1727: suffix: landau_damping_3d
1728: requires: !defined(PETSCTEST_VALGRIND)
1729: args: -xdm_plex_dim 3 -xdm_plex_box_faces 6,3,3 -xdm_plex_box_lower 0.,-1,-1 -xdm_plex_box_upper 12.5664,1,1 -xdm_plex_box_bd periodic,periodic,periodic -vdm_plex_dim 3 -vdm_plex_box_faces 4,4,4 \
1730: -vdm_plex_box_lower -6,-6,-6 -vdm_plex_box_upper 6,6,6 -vdm_plex_hash_location -vdm_plex_simplex 0 -vdm_plex_box_bd none,none,none
1732: test:
1733: requires: !defined(PETSC_USE_DMLANDAU_2D)
1734: suffix: sphere_3d
1735: nsize: 1
1736: args: -use_landau_velocity_space -xdm_plex_dim 3 -vdm_landau_thermal_temps 1 -vdm_landau_device_type cpu -xdm_plex_box_faces 6,3,3 -xdm_plex_box_lower 0.,-1,-1 \
1737: -xdm_plex_box_upper 12.5664,1,1 -xdm_plex_box_bd periodic,periodic,periodic -vdm_landau_verbose 2 -vdm_landau_sphere -vdm_landau_map_sphere \
1738: -vdm_landau_domain_radius 6,6,6 -vdm_landau_sphere_inner_radius_90degree_scale .35 -vdm_refine 1
1740: TEST*/