Actual source code: ex2.c
1: static char help[] = "Landau Damping/Two Stream instability test using Vlasov-Poisson equations\n";
3: /*
4: To run the code with particles sinusoidally perturbed in x space use the test "pp_poisson_bsi_1d_4" or "pp_poisson_bsi_2d_4"
5: According to Lukas, good damping results come at ~16k particles per cell
7: To visualize the efield use
9: -monitor_efield
11: To visualize the swarm distribution use
13: -ts_monitor_hg_swarm
15: To visualize the particles, we can use
17: -ts_monitor_sp_swarm -ts_monitor_sp_swarm_retain 0 -ts_monitor_sp_swarm_phase 1 -draw_size 500,500
19: */
20: #include <petscts.h>
21: #include <petscdmplex.h>
22: #include <petscdmswarm.h>
23: #include <petscfe.h>
24: #include <petscds.h>
25: #include <petscbag.h>
26: #include <petscdraw.h>
27: #include <petsc/private/dmpleximpl.h>
28: #include <petsc/private/petscfeimpl.h>
29: #include "petscdm.h"
30: #include "petscdmlabel.h"
32: PETSC_EXTERN PetscErrorCode stream(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
33: PETSC_EXTERN PetscErrorCode line(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
35: const char *EMTypes[] = {"primal", "mixed", "coulomb", "none", "EMType", "EM_", NULL};
36: typedef enum {
37: EM_PRIMAL,
38: EM_MIXED,
39: EM_COULOMB,
40: EM_NONE
41: } EMType;
43: typedef enum {
44: V0,
45: X0,
46: T0,
47: M0,
48: Q0,
49: PHI0,
50: POISSON,
51: VLASOV,
52: SIGMA,
53: NUM_CONSTANTS
54: } ConstantType;
55: typedef struct {
56: PetscScalar v0; /* Velocity scale, often the thermal velocity */
57: PetscScalar t0; /* Time scale */
58: PetscScalar x0; /* Space scale */
59: PetscScalar m0; /* Mass scale */
60: PetscScalar q0; /* Charge scale */
61: PetscScalar kb;
62: PetscScalar epsi0;
63: PetscScalar phi0; /* Potential scale */
64: PetscScalar poissonNumber; /* Non-Dimensional Poisson Number */
65: PetscScalar vlasovNumber; /* Non-Dimensional Vlasov Number */
66: PetscReal sigma; /* Nondimensional charge per length in x */
67: } Parameter;
69: typedef struct {
70: PetscBag bag; /* Problem parameters */
71: PetscBool error; /* Flag for printing the error */
72: PetscBool efield_monitor; /* Flag to show electric field monitor */
73: PetscBool initial_monitor;
74: PetscBool fake_1D; /* Run simulation in 2D but zeroing second dimension */
75: PetscBool perturbed_weights; /* Uniformly sample x,v space with gaussian weights */
76: PetscBool poisson_monitor;
77: PetscInt ostep; /* print the energy at each ostep time steps */
78: PetscInt numParticles;
79: PetscReal timeScale; /* Nondimensionalizing time scale */
80: PetscReal charges[2]; /* The charges of each species */
81: PetscReal masses[2]; /* The masses of each species */
82: PetscReal thermal_energy[2]; /* Thermal Energy (used to get other constants)*/
83: PetscReal cosine_coefficients[2]; /*(alpha, k)*/
84: PetscReal totalWeight;
85: PetscReal stepSize;
86: PetscInt steps;
87: PetscReal initVel;
88: EMType em; /* Type of electrostatic model */
89: SNES snes;
90: PetscDraw drawef;
91: PetscDrawLG drawlg_ef;
92: PetscDraw drawic_x;
93: PetscDraw drawic_v;
94: PetscDraw drawic_w;
95: PetscDrawHG drawhgic_x;
96: PetscDrawHG drawhgic_v;
97: PetscDrawHG drawhgic_w;
98: PetscDraw EDraw;
99: PetscDraw RhoDraw;
100: PetscDraw PotDraw;
101: PetscDrawSP EDrawSP;
102: PetscDrawSP RhoDrawSP;
103: PetscDrawSP PotDrawSP;
104: PetscBool monitor_positions; /* Flag to show particle positins at each time step */
105: PetscDraw positionDraw;
106: PetscDrawSP positionDrawSP;
107: DM swarm;
108: PetscRandom random;
109: PetscBool twostream;
110: PetscBool checkweights;
111: PetscInt checkVRes; /* Flag to check/output velocity residuals for nightly tests */
112: } AppCtx;
114: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
115: {
116: PetscFunctionBeginUser;
117: PetscInt d = 2;
118: PetscInt maxSpecies = 2;
119: options->error = PETSC_FALSE;
120: options->efield_monitor = PETSC_FALSE;
121: options->initial_monitor = PETSC_FALSE;
122: options->fake_1D = PETSC_FALSE;
123: options->perturbed_weights = PETSC_FALSE;
124: options->poisson_monitor = PETSC_FALSE;
125: options->ostep = 100;
126: options->timeScale = 2.0e-14;
127: options->charges[0] = -1.0;
128: options->charges[1] = 1.0;
129: options->masses[0] = 1.0;
130: options->masses[1] = 1000.0;
131: options->thermal_energy[0] = 1.0;
132: options->thermal_energy[1] = 1.0;
133: options->cosine_coefficients[0] = 0.01;
134: options->cosine_coefficients[1] = 0.5;
135: options->initVel = 1;
136: options->totalWeight = 1.0;
137: options->drawef = NULL;
138: options->drawlg_ef = NULL;
139: options->drawic_x = NULL;
140: options->drawic_v = NULL;
141: options->drawic_w = NULL;
142: options->drawhgic_x = NULL;
143: options->drawhgic_v = NULL;
144: options->drawhgic_w = NULL;
145: options->EDraw = NULL;
146: options->RhoDraw = NULL;
147: options->PotDraw = NULL;
148: options->EDrawSP = NULL;
149: options->RhoDrawSP = NULL;
150: options->PotDrawSP = NULL;
151: options->em = EM_COULOMB;
152: options->numParticles = 32768;
153: options->monitor_positions = PETSC_FALSE;
154: options->positionDraw = NULL;
155: options->positionDrawSP = NULL;
156: options->twostream = PETSC_FALSE;
157: options->checkweights = PETSC_FALSE;
158: options->checkVRes = 0;
160: PetscOptionsBegin(comm, "", "Landau Damping and Two Stream options", "DMSWARM");
161: PetscCall(PetscOptionsBool("-error", "Flag to print the error", "ex2.c", options->error, &options->error, NULL));
162: PetscCall(PetscOptionsBool("-monitor_efield", "Flag to show efield plot", "ex2.c", options->efield_monitor, &options->efield_monitor, NULL));
163: PetscCall(PetscOptionsBool("-monitor_ics", "Flag to show initial condition histograms", "ex2.c", options->initial_monitor, &options->initial_monitor, NULL));
164: PetscCall(PetscOptionsBool("-monitor_positions", "The flag to show particle positions", "ex2.c", options->monitor_positions, &options->monitor_positions, NULL));
165: PetscCall(PetscOptionsBool("-monitor_poisson", "The flag to show charges, Efield and potential solve", "ex2.c", options->poisson_monitor, &options->poisson_monitor, NULL));
166: PetscCall(PetscOptionsBool("-fake_1D", "Flag to run a 1D simulation (but really in 2D)", "ex2.c", options->fake_1D, &options->fake_1D, NULL));
167: PetscCall(PetscOptionsBool("-twostream", "Run two stream instability", "ex2.c", options->twostream, &options->twostream, NULL));
168: PetscCall(PetscOptionsBool("-perturbed_weights", "Flag to run uniform sampling with perturbed weights", "ex2.c", options->perturbed_weights, &options->perturbed_weights, NULL));
169: PetscCall(PetscOptionsBool("-check_weights", "Ensure all particle weights are positive", "ex2.c", options->checkweights, &options->checkweights, NULL));
170: PetscCall(PetscOptionsInt("-check_vel_res", "Check particle velocity residuals for nightly tests", "ex2.c", options->checkVRes, &options->checkVRes, NULL));
171: PetscCall(PetscOptionsInt("-output_step", "Number of time steps between output", "ex2.c", options->ostep, &options->ostep, NULL));
172: PetscCall(PetscOptionsReal("-timeScale", "Nondimensionalizing time scale", "ex2.c", options->timeScale, &options->timeScale, NULL));
173: PetscCall(PetscOptionsReal("-initial_velocity", "Initial velocity of perturbed particle", "ex2.c", options->initVel, &options->initVel, NULL));
174: PetscCall(PetscOptionsReal("-total_weight", "Total weight of all particles", "ex2.c", options->totalWeight, &options->totalWeight, NULL));
175: PetscCall(PetscOptionsRealArray("-cosine_coefficients", "Amplitude and frequency of cosine equation used in initialization", "ex2.c", options->cosine_coefficients, &d, NULL));
176: PetscCall(PetscOptionsRealArray("-charges", "Species charges", "ex2.c", options->charges, &maxSpecies, NULL));
177: PetscCall(PetscOptionsEnum("-em_type", "Type of electrostatic solver", "ex2.c", EMTypes, (PetscEnum)options->em, (PetscEnum *)&options->em, NULL));
178: PetscOptionsEnd();
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: static PetscErrorCode SetupContext(DM dm, DM sw, AppCtx *user)
183: {
184: PetscFunctionBeginUser;
185: if (user->efield_monitor) {
186: PetscDrawAxis axis_ef;
187: PetscCall(PetscDrawCreate(PETSC_COMM_WORLD, NULL, "monitor_efield", 0, 300, 400, 300, &user->drawef));
188: PetscCall(PetscDrawSetSave(user->drawef, "ex9_Efield.png"));
189: PetscCall(PetscDrawSetFromOptions(user->drawef));
190: PetscCall(PetscDrawLGCreate(user->drawef, 1, &user->drawlg_ef));
191: PetscCall(PetscDrawLGGetAxis(user->drawlg_ef, &axis_ef));
192: PetscCall(PetscDrawAxisSetLabels(axis_ef, "Electron Electric Field", "time", "E_max"));
193: PetscCall(PetscDrawLGSetLimits(user->drawlg_ef, 0., user->steps * user->stepSize, -10., 0.));
194: PetscCall(PetscDrawAxisSetLimits(axis_ef, 0., user->steps * user->stepSize, -10., 0.));
195: }
196: if (user->initial_monitor) {
197: PetscDrawAxis axis1, axis2, axis3;
198: PetscReal dmboxlower[2], dmboxupper[2];
199: PetscInt dim, cStart, cEnd;
200: PetscCall(DMGetDimension(sw, &dim));
201: PetscCall(DMGetBoundingBox(dm, dmboxlower, dmboxupper));
202: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
204: PetscCall(PetscDrawCreate(PETSC_COMM_WORLD, NULL, "monitor_initial_conditions_x", 0, 300, 400, 300, &user->drawic_x));
205: PetscCall(PetscDrawSetSave(user->drawic_x, "ex9_ic_x.png"));
206: PetscCall(PetscDrawSetFromOptions(user->drawic_x));
207: PetscCall(PetscDrawHGCreate(user->drawic_x, (int)dim, &user->drawhgic_x));
208: PetscCall(PetscDrawHGGetAxis(user->drawhgic_x, &axis1));
209: PetscCall(PetscDrawHGSetNumberBins(user->drawhgic_x, (int)(cEnd - cStart)));
210: PetscCall(PetscDrawAxisSetLabels(axis1, "Initial X Distribution", "X", "counts"));
211: PetscCall(PetscDrawAxisSetLimits(axis1, dmboxlower[0], dmboxupper[0], 0, 1500));
213: PetscCall(PetscDrawCreate(PETSC_COMM_WORLD, NULL, "monitor_initial_conditions_v", 400, 300, 400, 300, &user->drawic_v));
214: PetscCall(PetscDrawSetSave(user->drawic_v, "ex9_ic_v.png"));
215: PetscCall(PetscDrawSetFromOptions(user->drawic_v));
216: PetscCall(PetscDrawHGCreate(user->drawic_v, (int)dim, &user->drawhgic_v));
217: PetscCall(PetscDrawHGGetAxis(user->drawhgic_v, &axis2));
218: PetscCall(PetscDrawHGSetNumberBins(user->drawhgic_v, 1000));
219: PetscCall(PetscDrawAxisSetLabels(axis2, "Initial V_x Distribution", "V", "counts"));
220: PetscCall(PetscDrawAxisSetLimits(axis2, -1, 1, 0, 1500));
222: PetscCall(PetscDrawCreate(PETSC_COMM_WORLD, NULL, "monitor_initial_conditions_w", 800, 300, 400, 300, &user->drawic_w));
223: PetscCall(PetscDrawSetSave(user->drawic_w, "ex9_ic_w.png"));
224: PetscCall(PetscDrawSetFromOptions(user->drawic_w));
225: PetscCall(PetscDrawHGCreate(user->drawic_w, (int)dim, &user->drawhgic_w));
226: PetscCall(PetscDrawHGGetAxis(user->drawhgic_w, &axis3));
227: PetscCall(PetscDrawHGSetNumberBins(user->drawhgic_w, 10));
228: PetscCall(PetscDrawAxisSetLabels(axis3, "Initial W Distribution", "weight", "counts"));
229: PetscCall(PetscDrawAxisSetLimits(axis3, 0, 0.01, 0, 5000));
230: }
231: if (user->monitor_positions) {
232: PetscDrawAxis axis;
234: PetscCall(PetscDrawCreate(PETSC_COMM_WORLD, NULL, "position_monitor_species1", 0, 0, 400, 300, &user->positionDraw));
235: PetscCall(PetscDrawSetFromOptions(user->positionDraw));
236: PetscCall(PetscDrawSPCreate(user->positionDraw, 10, &user->positionDrawSP));
237: PetscCall(PetscDrawSPSetDimension(user->positionDrawSP, 1));
238: PetscCall(PetscDrawSPGetAxis(user->positionDrawSP, &axis));
239: PetscCall(PetscDrawSPReset(user->positionDrawSP));
240: PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "x", "v"));
241: PetscCall(PetscDrawSetSave(user->positionDraw, "ex9_pos.png"));
242: }
243: if (user->poisson_monitor) {
244: PetscDrawAxis axis_E, axis_Rho, axis_Pot;
246: PetscCall(PetscDrawCreate(PETSC_COMM_WORLD, NULL, "Efield_monitor", 0, 0, 400, 300, &user->EDraw));
247: PetscCall(PetscDrawSetFromOptions(user->EDraw));
248: PetscCall(PetscDrawSPCreate(user->EDraw, 10, &user->EDrawSP));
249: PetscCall(PetscDrawSPSetDimension(user->EDrawSP, 1));
250: PetscCall(PetscDrawSPGetAxis(user->EDrawSP, &axis_E));
251: PetscCall(PetscDrawSPReset(user->EDrawSP));
252: PetscCall(PetscDrawAxisSetLabels(axis_E, "Particles", "x", "E"));
253: PetscCall(PetscDrawSetSave(user->EDraw, "ex9_E_spatial.png"));
255: PetscCall(PetscDrawCreate(PETSC_COMM_WORLD, NULL, "rho_monitor", 0, 0, 400, 300, &user->RhoDraw));
256: PetscCall(PetscDrawSetFromOptions(user->RhoDraw));
257: PetscCall(PetscDrawSPCreate(user->RhoDraw, 10, &user->RhoDrawSP));
258: PetscCall(PetscDrawSPSetDimension(user->RhoDrawSP, 1));
259: PetscCall(PetscDrawSPGetAxis(user->RhoDrawSP, &axis_Rho));
260: PetscCall(PetscDrawSPReset(user->RhoDrawSP));
261: PetscCall(PetscDrawAxisSetLabels(axis_Rho, "Particles", "x", "rho"));
262: PetscCall(PetscDrawSetSave(user->RhoDraw, "ex9_rho_spatial.png"));
264: PetscCall(PetscDrawCreate(PETSC_COMM_WORLD, NULL, "potential_monitor", 0, 0, 400, 300, &user->PotDraw));
265: PetscCall(PetscDrawSetFromOptions(user->PotDraw));
266: PetscCall(PetscDrawSPCreate(user->PotDraw, 10, &user->PotDrawSP));
267: PetscCall(PetscDrawSPSetDimension(user->PotDrawSP, 1));
268: PetscCall(PetscDrawSPGetAxis(user->PotDrawSP, &axis_Pot));
269: PetscCall(PetscDrawSPReset(user->PotDrawSP));
270: PetscCall(PetscDrawAxisSetLabels(axis_Pot, "Particles", "x", "potential"));
271: PetscCall(PetscDrawSetSave(user->PotDraw, "ex9_phi_spatial.png"));
272: }
273: PetscFunctionReturn(PETSC_SUCCESS);
274: }
276: static PetscErrorCode DestroyContext(AppCtx *user)
277: {
278: PetscFunctionBeginUser;
279: PetscCall(PetscDrawLGDestroy(&user->drawlg_ef));
280: PetscCall(PetscDrawDestroy(&user->drawef));
281: PetscCall(PetscDrawHGDestroy(&user->drawhgic_x));
282: PetscCall(PetscDrawDestroy(&user->drawic_x));
283: PetscCall(PetscDrawHGDestroy(&user->drawhgic_v));
284: PetscCall(PetscDrawDestroy(&user->drawic_v));
285: PetscCall(PetscDrawHGDestroy(&user->drawhgic_w));
286: PetscCall(PetscDrawDestroy(&user->drawic_w));
287: PetscCall(PetscDrawSPDestroy(&user->positionDrawSP));
288: PetscCall(PetscDrawDestroy(&user->positionDraw));
290: PetscCall(PetscDrawSPDestroy(&user->EDrawSP));
291: PetscCall(PetscDrawDestroy(&user->EDraw));
292: PetscCall(PetscDrawSPDestroy(&user->RhoDrawSP));
293: PetscCall(PetscDrawDestroy(&user->RhoDraw));
294: PetscCall(PetscDrawSPDestroy(&user->PotDrawSP));
295: PetscCall(PetscDrawDestroy(&user->PotDraw));
297: PetscCall(PetscBagDestroy(&user->bag));
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
301: static PetscErrorCode CheckNonNegativeWeights(DM sw, AppCtx *user)
302: {
303: const PetscScalar *w;
304: PetscInt Np;
306: PetscFunctionBeginUser;
307: if (!user->checkweights) PetscFunctionReturn(PETSC_SUCCESS);
308: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&w));
309: PetscCall(DMSwarmGetLocalSize(sw, &Np));
310: 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]);
311: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: static PetscErrorCode computeParticleMoments(DM sw, PetscReal moments[3], AppCtx *user)
316: {
317: DM dm;
318: const PetscReal *coords;
319: const PetscScalar *w;
320: PetscReal mom[3] = {0.0, 0.0, 0.0};
321: PetscInt cell, cStart, cEnd, dim;
323: PetscFunctionBeginUser;
324: PetscCall(DMGetDimension(sw, &dim));
325: PetscCall(DMSwarmGetCellDM(sw, &dm));
326: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
327: PetscCall(DMSwarmSortGetAccess(sw));
328: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&coords));
329: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&w));
330: for (cell = cStart; cell < cEnd; ++cell) {
331: PetscInt *pidx;
332: PetscInt Np, p, d;
334: PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx));
335: for (p = 0; p < Np; ++p) {
336: const PetscInt idx = pidx[p];
337: const PetscReal *c = &coords[idx * dim];
339: mom[0] += PetscRealPart(w[idx]);
340: mom[1] += PetscRealPart(w[idx]) * c[0];
341: for (d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx]) * c[d] * c[d];
342: //if (w[idx] < 0. ) PetscPrintf(PETSC_COMM_WORLD, "error, negative weight %" PetscInt_FMT " \n", idx);
343: }
344: PetscCall(PetscFree(pidx));
345: }
346: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&coords));
347: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
348: PetscCall(DMSwarmSortRestoreAccess(sw));
349: PetscCallMPI(MPIU_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: static void f0_1(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[])
354: {
355: f0[0] = u[0];
356: }
358: static void f0_x(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[])
359: {
360: f0[0] = x[0] * u[0];
361: }
363: static void f0_r2(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[])
364: {
365: PetscInt d;
367: f0[0] = 0.0;
368: for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d]) * u[0];
369: }
371: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
372: {
373: PetscDS prob;
374: PetscScalar mom;
375: PetscInt field = 0;
377: PetscFunctionBeginUser;
378: PetscCall(DMGetDS(dm, &prob));
379: PetscCall(PetscDSSetObjective(prob, field, &f0_1));
380: PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
381: moments[0] = PetscRealPart(mom);
382: PetscCall(PetscDSSetObjective(prob, field, &f0_x));
383: PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
384: moments[1] = PetscRealPart(mom);
385: PetscCall(PetscDSSetObjective(prob, field, &f0_r2));
386: PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
387: moments[2] = PetscRealPart(mom);
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: static PetscErrorCode MonitorEField(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
392: {
393: AppCtx *user = (AppCtx *)ctx;
394: DM dm, sw;
395: PetscReal *E;
396: PetscReal Enorm = 0., lgEnorm, lgEmax, sum = 0., Emax = 0., temp = 0., *weight, chargesum = 0.;
397: PetscReal *x, *v;
398: PetscInt *species, dim, p, d, Np, cStart, cEnd;
399: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
400: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
401: Vec rho;
403: PetscFunctionBeginUser;
404: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS);
405: PetscCall(TSGetDM(ts, &sw));
406: PetscCall(DMSwarmGetCellDM(sw, &dm));
407: PetscCall(DMGetDimension(sw, &dim));
408: PetscCall(DMSwarmGetLocalSize(sw, &Np));
409: PetscCall(DMSwarmSortGetAccess(sw));
410: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
411: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
412: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
413: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
414: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
415: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
417: for (p = 0; p < Np; ++p) {
418: for (d = 0; d < 1; ++d) {
419: temp = PetscAbsReal(E[p * dim + d]);
420: if (temp > Emax) Emax = temp;
421: }
422: Enorm += PetscSqrtReal(E[p * dim] * E[p * dim]);
423: sum += E[p * dim];
424: chargesum += user->charges[0] * weight[p];
425: }
426: lgEnorm = Enorm != 0 ? PetscLog10Real(Enorm) : -16.;
427: lgEmax = Emax != 0 ? PetscLog10Real(Emax) : -16.;
429: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
430: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
431: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
432: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
433: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
435: Parameter *param;
436: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
437: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "charges", &rho));
438: if (user->em == EM_PRIMAL) {
439: PetscCall(computeParticleMoments(sw, pmoments, user));
440: PetscCall(computeFEMMoments(dm, rho, fmoments, user));
441: } else if (user->em == EM_MIXED) {
442: DM potential_dm;
443: IS potential_IS;
444: PetscInt fields = 1;
445: PetscCall(DMCreateSubDM(dm, 1, &fields, &potential_IS, &potential_dm));
447: PetscCall(computeParticleMoments(sw, pmoments, user));
448: PetscCall(computeFEMMoments(potential_dm, rho, fmoments, user));
449: PetscCall(DMDestroy(&potential_dm));
450: PetscCall(ISDestroy(&potential_IS));
451: }
452: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "charges", &rho));
454: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%f\t%+e\t%e\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", (double)t, (double)sum, (double)Enorm, (double)lgEnorm, (double)Emax, (double)lgEmax, (double)chargesum, (double)pmoments[0], (double)pmoments[1], (double)pmoments[2], (double)fmoments[0], (double)fmoments[1], (double)fmoments[2]));
455: PetscCall(PetscDrawLGAddPoint(user->drawlg_ef, &t, &lgEmax));
456: PetscCall(PetscDrawLGDraw(user->drawlg_ef));
457: PetscCall(PetscDrawSave(user->drawef));
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: PetscErrorCode MonitorInitialConditions(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
462: {
463: AppCtx *user = (AppCtx *)ctx;
464: DM dm, sw;
465: const PetscScalar *u;
466: PetscReal *weight, *pos, *vel;
467: PetscInt dim, p, Np, cStart, cEnd;
469: PetscFunctionBegin;
470: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* -1 indicates interpolated solution */
471: PetscCall(TSGetDM(ts, &sw));
472: PetscCall(DMSwarmGetCellDM(sw, &dm));
473: PetscCall(DMGetDimension(sw, &dim));
474: PetscCall(DMSwarmGetLocalSize(sw, &Np));
475: PetscCall(DMSwarmSortGetAccess(sw));
476: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
478: if (step == 0) {
479: PetscCall(PetscDrawHGReset(user->drawhgic_x));
480: PetscCall(PetscDrawHGGetDraw(user->drawhgic_x, &user->drawic_x));
481: PetscCall(PetscDrawClear(user->drawic_x));
482: PetscCall(PetscDrawFlush(user->drawic_x));
484: PetscCall(PetscDrawHGReset(user->drawhgic_v));
485: PetscCall(PetscDrawHGGetDraw(user->drawhgic_v, &user->drawic_v));
486: PetscCall(PetscDrawClear(user->drawic_v));
487: PetscCall(PetscDrawFlush(user->drawic_v));
489: PetscCall(PetscDrawHGReset(user->drawhgic_w));
490: PetscCall(PetscDrawHGGetDraw(user->drawhgic_w, &user->drawic_w));
491: PetscCall(PetscDrawClear(user->drawic_w));
492: PetscCall(PetscDrawFlush(user->drawic_w));
494: PetscCall(VecGetArrayRead(U, &u));
495: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&vel));
496: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
497: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&pos));
499: PetscCall(VecGetLocalSize(U, &Np));
500: Np /= dim * 2;
501: for (p = 0; p < Np; ++p) {
502: PetscCall(PetscDrawHGAddValue(user->drawhgic_x, pos[p * dim]));
503: PetscCall(PetscDrawHGAddValue(user->drawhgic_v, vel[p * dim]));
504: PetscCall(PetscDrawHGAddValue(user->drawhgic_w, weight[p]));
505: }
507: PetscCall(VecRestoreArrayRead(U, &u));
508: PetscCall(PetscDrawHGDraw(user->drawhgic_x));
509: PetscCall(PetscDrawHGSave(user->drawhgic_x));
511: PetscCall(PetscDrawHGDraw(user->drawhgic_v));
512: PetscCall(PetscDrawHGSave(user->drawhgic_v));
514: PetscCall(PetscDrawHGDraw(user->drawhgic_w));
515: PetscCall(PetscDrawHGSave(user->drawhgic_w));
517: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&pos));
518: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&vel));
519: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
520: }
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: static PetscErrorCode MonitorPositions_2D(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
525: {
526: AppCtx *user = (AppCtx *)ctx;
527: DM dm, sw;
528: PetscScalar *x, *v, *weight;
529: PetscReal lower[3], upper[3], speed;
530: const PetscInt *s;
531: PetscInt dim, cStart, cEnd, c;
533: PetscFunctionBeginUser;
534: if (step > 0 && step % user->ostep == 0) {
535: PetscCall(TSGetDM(ts, &sw));
536: PetscCall(DMSwarmGetCellDM(sw, &dm));
537: PetscCall(DMGetDimension(dm, &dim));
538: PetscCall(DMGetBoundingBox(dm, lower, upper));
539: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
540: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
541: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
542: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
543: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&s));
544: PetscCall(DMSwarmSortGetAccess(sw));
545: PetscCall(PetscDrawSPReset(user->positionDrawSP));
546: PetscCall(PetscDrawSPSetLimits(user->positionDrawSP, lower[0], upper[0], lower[1], upper[1]));
547: PetscCall(PetscDrawSPSetLimits(user->positionDrawSP, lower[0], upper[0], -12, 12));
548: for (c = 0; c < cEnd - cStart; ++c) {
549: PetscInt *pidx, Npc, q;
550: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
551: for (q = 0; q < Npc; ++q) {
552: const PetscInt p = pidx[q];
553: if (s[p] == 0) {
554: speed = PetscSqrtReal(PetscSqr(v[p * dim]) + PetscSqr(v[p * dim + 1]));
555: if (dim == 1 || user->fake_1D) {
556: PetscCall(PetscDrawSPAddPointColorized(user->positionDrawSP, &x[p * dim], &v[p * dim], &speed));
557: } else {
558: PetscCall(PetscDrawSPAddPointColorized(user->positionDrawSP, &x[p * dim], &x[p * dim + 1], &speed));
559: }
560: } else if (s[p] == 1) {
561: PetscCall(PetscDrawSPAddPoint(user->positionDrawSP, &x[p * dim], &v[p * dim]));
562: }
563: }
564: PetscCall(PetscFree(pidx));
565: }
566: PetscCall(PetscDrawSPDraw(user->positionDrawSP, PETSC_TRUE));
567: PetscCall(PetscDrawSave(user->positionDraw));
568: PetscCall(DMSwarmSortRestoreAccess(sw));
569: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
570: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
571: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
572: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&s));
573: }
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: static PetscErrorCode MonitorPoisson(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
578: {
579: AppCtx *user = (AppCtx *)ctx;
580: DM dm, sw;
581: PetscScalar *x, *E, *weight, *pot, *charges;
582: PetscReal lower[3], upper[3], xval;
583: PetscInt dim, cStart, cEnd, c;
585: PetscFunctionBeginUser;
586: if (step > 0 && step % user->ostep == 0) {
587: PetscCall(TSGetDM(ts, &sw));
588: PetscCall(DMSwarmGetCellDM(sw, &dm));
589: PetscCall(DMGetDimension(dm, &dim));
590: PetscCall(DMGetBoundingBox(dm, lower, upper));
591: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
593: PetscCall(PetscDrawSPReset(user->RhoDrawSP));
594: PetscCall(PetscDrawSPReset(user->EDrawSP));
595: PetscCall(PetscDrawSPReset(user->PotDrawSP));
596: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
597: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
598: PetscCall(DMSwarmGetField(sw, "potential", NULL, NULL, (void **)&pot));
599: PetscCall(DMSwarmGetField(sw, "charges", NULL, NULL, (void **)&charges));
600: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
602: PetscCall(DMSwarmSortGetAccess(sw));
603: for (c = 0; c < cEnd - cStart; ++c) {
604: PetscReal Esum = 0.0;
605: PetscInt *pidx, Npc, q;
606: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
607: for (q = 0; q < Npc; ++q) {
608: const PetscInt p = pidx[q];
609: Esum += E[p * dim];
610: }
611: xval = (c + 0.5) * ((upper[0] - lower[0]) / (cEnd - cStart));
612: PetscCall(PetscDrawSPAddPoint(user->EDrawSP, &xval, &Esum));
613: PetscCall(PetscFree(pidx));
614: }
615: for (c = 0; c < (cEnd - cStart); ++c) {
616: xval = (c + 0.5) * ((upper[0] - lower[0]) / (cEnd - cStart));
617: PetscCall(PetscDrawSPAddPoint(user->RhoDrawSP, &xval, &charges[c]));
618: PetscCall(PetscDrawSPAddPoint(user->PotDrawSP, &xval, &pot[c]));
619: }
620: PetscCall(PetscDrawSPDraw(user->RhoDrawSP, PETSC_TRUE));
621: PetscCall(PetscDrawSave(user->RhoDraw));
622: PetscCall(PetscDrawSPDraw(user->EDrawSP, PETSC_TRUE));
623: PetscCall(PetscDrawSave(user->EDraw));
624: PetscCall(PetscDrawSPDraw(user->PotDrawSP, PETSC_TRUE));
625: PetscCall(PetscDrawSave(user->PotDraw));
626: PetscCall(DMSwarmSortRestoreAccess(sw));
627: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
628: PetscCall(DMSwarmRestoreField(sw, "potential", NULL, NULL, (void **)&pot));
629: PetscCall(DMSwarmRestoreField(sw, "charges", NULL, NULL, (void **)&charges));
630: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
631: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
632: }
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
637: {
638: PetscBag bag;
639: Parameter *p;
641: PetscFunctionBeginUser;
642: /* setup PETSc parameter bag */
643: PetscCall(PetscBagGetData(ctx->bag, (void **)&p));
644: PetscCall(PetscBagSetName(ctx->bag, "par", "Vlasov-Poisson Parameters"));
645: bag = ctx->bag;
646: PetscCall(PetscBagRegisterScalar(bag, &p->v0, 1.0, "v0", "Velocity scale, m/s"));
647: PetscCall(PetscBagRegisterScalar(bag, &p->t0, 1.0, "t0", "Time scale, s"));
648: PetscCall(PetscBagRegisterScalar(bag, &p->x0, 1.0, "x0", "Space scale, m"));
649: PetscCall(PetscBagRegisterScalar(bag, &p->v0, 1.0, "phi0", "Potential scale, kg*m^2/A*s^3"));
650: PetscCall(PetscBagRegisterScalar(bag, &p->q0, 1.0, "q0", "Charge Scale, A*s"));
651: PetscCall(PetscBagRegisterScalar(bag, &p->m0, 1.0, "m0", "Mass Scale, kg"));
652: PetscCall(PetscBagRegisterScalar(bag, &p->epsi0, 1.0, "epsi0", "Permittivity of Free Space, kg"));
653: PetscCall(PetscBagRegisterScalar(bag, &p->kb, 1.0, "kb", "Boltzmann Constant, m^2 kg/s^2 K^1"));
655: PetscCall(PetscBagRegisterScalar(bag, &p->sigma, 1.0, "sigma", "Charge per unit area, C/m^3"));
656: PetscCall(PetscBagRegisterScalar(bag, &p->poissonNumber, 1.0, "poissonNumber", "Non-Dimensional Poisson Number"));
657: PetscCall(PetscBagRegisterScalar(bag, &p->vlasovNumber, 1.0, "vlasovNumber", "Non-Dimensional Vlasov Number"));
658: PetscCall(PetscBagSetFromOptions(bag));
659: {
660: PetscViewer viewer;
661: PetscViewerFormat format;
662: PetscBool flg;
664: PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg));
665: if (flg) {
666: PetscCall(PetscViewerPushFormat(viewer, format));
667: PetscCall(PetscBagView(bag, viewer));
668: PetscCall(PetscViewerFlush(viewer));
669: PetscCall(PetscViewerPopFormat(viewer));
670: PetscCall(PetscViewerDestroy(&viewer));
671: }
672: }
673: PetscFunctionReturn(PETSC_SUCCESS);
674: }
676: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
677: {
678: PetscFunctionBeginUser;
679: PetscCall(DMCreate(comm, dm));
680: PetscCall(DMSetType(*dm, DMPLEX));
681: PetscCall(DMSetFromOptions(*dm));
682: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: 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[])
687: {
688: f0[0] = -constants[SIGMA];
689: }
691: 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[])
692: {
693: PetscInt d;
694: for (d = 0; d < dim; ++d) f1[d] = u_x[d];
695: }
697: 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[])
698: {
699: PetscInt d;
700: for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
701: }
703: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
704: {
705: *u = 0.0;
706: return PETSC_SUCCESS;
707: }
709: /*
710: / I -grad\ / q \ = /0\
711: \-div 0 / \phi/ \f/
712: */
713: 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[])
714: {
715: for (PetscInt d = 0; d < dim; ++d) f0[d] += u[uOff[0] + d];
716: }
718: 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[])
719: {
720: for (PetscInt d = 0; d < dim; ++d) f1[d * dim + d] = u[uOff[1]];
721: }
723: static void f0_phi_backgroundCharge(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[])
724: {
725: f0[0] += constants[SIGMA];
726: for (PetscInt d = 0; d < dim; ++d) f0[0] += u_x[uOff_x[0] + d * dim + d];
727: }
729: /* Boundary residual. Dirichlet boundary for u means u_bdy=p*n */
730: 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[])
731: {
732: for (PetscInt d = 0; d < dim; ++d) g0[d * dim + d] = 1.0;
733: }
735: static void g2_qphi(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[])
736: {
737: for (PetscInt d = 0; d < dim; ++d) g2[d * dim + d] = 1.0;
738: }
740: static void g1_phiq(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[])
741: {
742: for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
743: }
745: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
746: {
747: PetscFE fephi, feq;
748: PetscDS ds;
749: PetscBool simplex;
750: PetscInt dim;
752: PetscFunctionBeginUser;
753: PetscCall(DMGetDimension(dm, &dim));
754: PetscCall(DMPlexIsSimplex(dm, &simplex));
755: if (user->em == EM_MIXED) {
756: DMLabel label;
757: const PetscInt id = 1;
759: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "field_", PETSC_DETERMINE, &feq));
760: PetscCall(PetscObjectSetName((PetscObject)feq, "field"));
761: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "potential_", PETSC_DETERMINE, &fephi));
762: PetscCall(PetscObjectSetName((PetscObject)fephi, "potential"));
763: PetscCall(PetscFECopyQuadrature(feq, fephi));
764: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feq));
765: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fephi));
766: PetscCall(DMCreateDS(dm));
767: PetscCall(PetscFEDestroy(&fephi));
768: PetscCall(PetscFEDestroy(&feq));
770: PetscCall(DMGetLabel(dm, "marker", &label));
771: PetscCall(DMGetDS(dm, &ds));
773: PetscCall(PetscDSSetResidual(ds, 0, f0_q, f1_q));
774: PetscCall(PetscDSSetResidual(ds, 1, f0_phi_backgroundCharge, NULL));
775: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_qq, NULL, NULL, NULL));
776: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_qphi, NULL));
777: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_phiq, NULL, NULL));
779: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))zero, NULL, NULL, NULL));
781: } else if (user->em == EM_PRIMAL) {
782: MatNullSpace nullsp;
783: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &fephi));
784: PetscCall(PetscObjectSetName((PetscObject)fephi, "potential"));
785: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fephi));
786: PetscCall(DMCreateDS(dm));
787: PetscCall(DMGetDS(dm, &ds));
788: PetscCall(PetscDSSetResidual(ds, 0, ion_f0, laplacian_f1));
789: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, laplacian_g3));
790: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_TRUE, 0, NULL, &nullsp));
791: PetscCall(PetscObjectCompose((PetscObject)fephi, "nullspace", (PetscObject)nullsp));
792: PetscCall(MatNullSpaceDestroy(&nullsp));
793: PetscCall(PetscFEDestroy(&fephi));
794: }
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
798: static PetscErrorCode CreatePoisson(DM dm, AppCtx *user)
799: {
800: SNES snes;
801: Mat J;
802: MatNullSpace nullSpace;
804: PetscFunctionBeginUser;
805: PetscCall(CreateFEM(dm, user));
806: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
807: PetscCall(SNESSetOptionsPrefix(snes, "em_"));
808: PetscCall(SNESSetDM(snes, dm));
809: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, user));
810: PetscCall(SNESSetFromOptions(snes));
812: PetscCall(DMCreateMatrix(dm, &J));
813: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_TRUE, 0, NULL, &nullSpace));
814: PetscCall(MatSetNullSpace(J, nullSpace));
815: PetscCall(MatNullSpaceDestroy(&nullSpace));
816: PetscCall(SNESSetJacobian(snes, J, J, NULL, NULL));
817: PetscCall(MatDestroy(&J));
818: user->snes = snes;
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }
822: PetscErrorCode PetscPDFPertubedConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
823: {
824: p[0] = (1 + 0.01 * PetscCosReal(0.5 * x[0])) / (2 * PETSC_PI);
825: p[1] = (1 + 0.01 * PetscCosReal(0.5 * x[1])) / (2 * PETSC_PI);
826: return PETSC_SUCCESS;
827: }
828: PetscErrorCode PetscPDFPertubedConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
829: {
830: p[0] = (1. + 0.01 * PetscCosReal(0.5 * x[0])) / (2 * PETSC_PI);
831: return PETSC_SUCCESS;
832: }
834: PetscErrorCode PetscPDFCosine1D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
835: {
836: const PetscReal alpha = scale ? scale[0] : 0.0;
837: const PetscReal k = scale ? scale[1] : 1.;
838: p[0] = (1 + alpha * PetscCosReal(k * x[0]));
839: return PETSC_SUCCESS;
840: }
842: PetscErrorCode PetscPDFCosine2D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
843: {
844: const PetscReal alpha = scale ? scale[0] : 0.;
845: const PetscReal k = scale ? scale[0] : 1.;
846: p[0] = (1 + alpha * PetscCosReal(k * (x[0] + x[1])));
847: return PETSC_SUCCESS;
848: }
850: PetscErrorCode PetscPDFCosine1D_TwoStream(const PetscReal x[], const PetscReal scale[], PetscReal p[])
851: {
852: const PetscReal alpha = scale ? scale[0] : 0.0;
853: const PetscReal k = scale ? scale[1] : 1.;
854: p[0] = (1. + alpha * PetscCosReal(k * x[0]));
855: return PETSC_SUCCESS;
856: }
858: static PetscErrorCode InitializeParticles_PerturbedWeights(DM sw, AppCtx *user)
859: {
860: DM vdm, dm;
861: PetscScalar *weight;
862: PetscReal *x, *v, vmin[3], vmax[3], gmin[3], gmax[3], xi0[3];
863: PetscInt *N, Ns, dim, *cellid, *species, Np, cStart, cEnd, Npc, n;
864: PetscInt Np_global, p, q, s, c, d, cv;
865: PetscBool flg;
866: PetscMPIInt size, rank;
867: Parameter *param;
869: PetscFunctionBegin;
870: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sw), &size));
871: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
872: PetscOptionsBegin(PetscObjectComm((PetscObject)sw), "", "DMSwarm Options", "DMSWARM");
873: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
874: PetscCall(PetscOptionsInt("-dm_swarm_num_species", "The number of species", "DMSwarmSetNumSpecies", Ns, &Ns, &flg));
875: if (flg) PetscCall(DMSwarmSetNumSpecies(sw, Ns));
876: PetscCall(PetscCalloc1(Ns, &N));
877: n = Ns;
878: PetscCall(PetscOptionsIntArray("-dm_swarm_num_particles", "The target number of particles", "", N, &n, NULL));
879: PetscOptionsEnd();
881: PetscCall(DMGetDimension(sw, &dim));
882: PetscCall(DMSwarmGetCellDM(sw, &dm));
883: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
885: PetscCall(DMCreate(PETSC_COMM_SELF, &vdm));
886: PetscCall(DMSetType(vdm, DMPLEX));
887: PetscCall(DMPlexSetOptionsPrefix(vdm, "v"));
888: PetscCall(DMSetFromOptions(vdm));
889: PetscCall(DMViewFromOptions(vdm, NULL, "-vdm_view"));
891: PetscInt vStart, vEnd;
892: PetscCall(DMPlexGetHeightStratum(vdm, 0, &vStart, &vEnd));
893: PetscCall(DMGetBoundingBox(vdm, vmin, vmax));
895: PetscCall(DMGetBoundingBox(dm, gmin, gmax));
896: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
897: Np = (cEnd - cStart) * (vEnd - vStart);
898: PetscCallMPI(MPIU_Allreduce(&Np, &Np_global, 1, MPIU_INT, MPIU_SUM, PETSC_COMM_WORLD));
899: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Global Np = %" PetscInt_FMT "\n", Np_global));
900: PetscCall(DMSwarmSetLocalSizes(sw, Np, 0));
901: Npc = Np / (cEnd - cStart);
902: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
903: for (c = 0, p = 0; c < cEnd - cStart; ++c) {
904: for (s = 0; s < Ns; ++s) {
905: for (q = 0; q < Npc; ++q, ++p) cellid[p] = c;
906: }
907: }
908: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
909: PetscCall(PetscFree(N));
911: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
912: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
913: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
914: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
916: PetscCall(DMSwarmSortGetAccess(sw));
917: for (c = 0; c < cEnd - cStart; ++c) {
918: const PetscInt cell = c + cStart;
919: PetscInt *pidx, Npc;
920: PetscReal centroid[3], volume;
922: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
923: PetscCall(DMPlexComputeCellGeometryFVM(dm, cell, &volume, centroid, NULL));
924: for (q = 0; q < Npc; ++q) {
925: const PetscInt p = pidx[q];
927: for (d = 0; d < dim; ++d) {
928: x[p * dim + d] = centroid[d];
929: v[p * dim + d] = vmin[0] + (q + 0.5) * (vmax[0] - vmin[0]) / Npc;
930: if (user->fake_1D && d > 0) v[p * dim + d] = 0;
931: }
932: }
933: PetscCall(PetscFree(pidx));
934: }
935: PetscCall(DMGetCoordinatesLocalSetUp(vdm));
937: /* Setup Quadrature for spatial and velocity weight calculations*/
938: PetscQuadrature quad_x;
939: PetscInt Nq_x;
940: const PetscReal *wq_x, *xq_x;
941: PetscReal *xq_x_extended;
942: PetscReal weightsum = 0., totalcellweight = 0., *weight_x, *weight_v;
943: PetscReal scale[2] = {user->cosine_coefficients[0], user->cosine_coefficients[1]};
945: PetscCall(PetscCalloc2(cEnd - cStart, &weight_x, Np, &weight_v));
946: if (user->fake_1D) PetscCall(PetscDTGaussTensorQuadrature(1, 1, 5, -1.0, 1.0, &quad_x));
947: else PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 5, -1.0, 1.0, &quad_x));
948: PetscCall(PetscQuadratureGetData(quad_x, NULL, NULL, &Nq_x, &xq_x, &wq_x));
949: if (user->fake_1D) {
950: PetscCall(PetscCalloc1(Nq_x * dim, &xq_x_extended));
951: for (PetscInt i = 0; i < Nq_x; ++i) xq_x_extended[i * dim] = xq_x[i];
952: }
953: /* Integrate the density function to get the weights of particles in each cell */
954: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
955: for (c = cStart; c < cEnd; ++c) {
956: PetscReal v0_x[3], J_x[9], invJ_x[9], detJ_x, xr_x[3], den_x;
957: PetscInt *pidx, Npc, q;
958: PetscInt Ncx;
959: const PetscScalar *array_x;
960: PetscScalar *coords_x = NULL;
961: PetscBool isDGx;
962: weight_x[c] = 0.;
964: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDGx, &Ncx, &array_x, &coords_x));
965: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
966: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0_x, J_x, invJ_x, &detJ_x));
967: for (q = 0; q < Nq_x; ++q) {
968: /*Transform quadrature points from ref space to real space (0,12.5664)*/
969: if (user->fake_1D) CoordinatesRefToReal(dim, dim, xi0, v0_x, J_x, &xq_x_extended[q * dim], xr_x);
970: else CoordinatesRefToReal(dim, dim, xi0, v0_x, J_x, &xq_x[q * dim], xr_x);
972: /*Transform quadrature points from real space to ideal real space (0, 2PI/k)*/
973: if (user->fake_1D) {
974: if (user->twostream) PetscCall(PetscPDFCosine1D_TwoStream(xr_x, scale, &den_x));
975: else PetscCall(PetscPDFCosine1D(xr_x, scale, &den_x));
976: detJ_x = J_x[0];
977: } else PetscCall(PetscPDFCosine2D(xr_x, scale, &den_x));
978: /*We have to transform the quadrature weights as well*/
979: weight_x[c] += den_x * (wq_x[q] * detJ_x);
980: }
981: // Get the cell numbering for consistent output between sequential and distributed tests
982: IS globalOrdering;
983: const PetscInt *ordering;
984: PetscCall(DMPlexGetCellNumbering(dm, &globalOrdering));
985: PetscCall(ISGetIndices(globalOrdering, &ordering));
986: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "c:%" PetscInt_FMT " [x_a,x_b] = %1.15f,%1.15f -> cell weight = %1.15f\n", ordering[c], (double)PetscRealPart(coords_x[0]), (double)PetscRealPart(coords_x[2]), (double)weight_x[c]));
987: PetscCall(ISRestoreIndices(globalOrdering, &ordering));
988: totalcellweight += weight_x[c];
989: // Confirm the number of particles per spatial cell conforms to the size of the velocity grid
990: PetscCheck(Npc == vEnd - vStart, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of particles %" PetscInt_FMT " in cell (rank %d/%d) != %" PetscInt_FMT " number of velocity vertices", Npc, rank, size, vEnd - vStart);
992: /* Set weights to be gaussian in velocity cells (using exact solution) */
993: for (cv = 0; cv < vEnd - vStart; ++cv) {
994: PetscInt Nc;
995: const PetscScalar *array_v;
996: PetscScalar *coords_v = NULL;
997: PetscBool isDG;
998: PetscCall(DMPlexGetCellCoordinates(vdm, cv, &isDG, &Nc, &array_v, &coords_v));
1000: const PetscInt p = pidx[cv];
1001: // Two stream function from 1/2pi v^2 e^(-v^2/2)
1002: if (user->twostream)
1003: weight_v[p] = 1. / (PetscSqrtReal(2 * PETSC_PI)) * ((coords_v[0] * PetscExpReal(-PetscSqr(coords_v[0]) / 2.)) - (coords_v[1] * PetscExpReal(-PetscSqr(coords_v[1]) / 2.))) - 0.5 * PetscErfReal(coords_v[0] / PetscSqrtReal(2.)) + 0.5 * (PetscErfReal(coords_v[1] / PetscSqrtReal(2.)));
1004: else weight_v[p] = 0.5 * (PetscErfReal(coords_v[1] / PetscSqrtReal(2.)) - PetscErfReal(coords_v[0] / PetscSqrtReal(2.)));
1006: weight[p] = user->totalWeight * weight_v[p] * weight_x[c];
1007: if (weight[p] > 1.) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "weights: %g, %g, %g\n", user->totalWeight, weight_v[p], weight_x[c]));
1008: //PetscPrintf(PETSC_COMM_WORLD, "particle %"PetscInt_FMT": %g, weight_v: %g weight_x: %g\n", p, weight[p], weight_v[p], weight_x[p]);
1009: weightsum += weight[p];
1011: PetscCall(DMPlexRestoreCellCoordinates(vdm, cv, &isDG, &Nc, &array_v, &coords_v));
1012: }
1013: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDGx, &Ncx, &array_x, &coords_x));
1014: PetscCall(PetscFree(pidx));
1015: }
1016: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
1017: PetscReal global_cellweight, global_weightsum;
1018: PetscCallMPI(MPIU_Allreduce(&totalcellweight, &global_cellweight, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
1019: PetscCallMPI(MPIU_Allreduce(&weightsum, &global_weightsum, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
1020: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "particle weight sum = %1.10f cell weight sum = %1.10f\n", (double)global_cellweight, (double)global_weightsum));
1021: if (user->fake_1D) PetscCall(PetscFree(xq_x_extended));
1022: PetscCall(PetscFree2(weight_x, weight_v));
1023: PetscCall(PetscQuadratureDestroy(&quad_x));
1024: PetscCall(DMSwarmSortRestoreAccess(sw));
1025: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
1026: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
1027: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1028: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
1029: PetscCall(DMDestroy(&vdm));
1030: PetscFunctionReturn(PETSC_SUCCESS);
1031: }
1033: static PetscErrorCode InitializeConstants(DM sw, AppCtx *user)
1034: {
1035: DM dm;
1036: PetscInt *species;
1037: PetscReal *weight, totalCharge = 0., totalWeight = 0., gmin[3], gmax[3], global_charge, global_weight;
1038: PetscInt Np, dim;
1040: PetscFunctionBegin;
1041: PetscCall(DMSwarmGetCellDM(sw, &dm));
1042: PetscCall(DMGetDimension(sw, &dim));
1043: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1044: PetscCall(DMGetBoundingBox(dm, gmin, gmax));
1045: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
1046: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1047: for (PetscInt p = 0; p < Np; ++p) {
1048: totalWeight += weight[p];
1049: totalCharge += user->charges[species[p]] * weight[p];
1050: }
1051: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
1052: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1053: {
1054: Parameter *param;
1055: PetscReal Area;
1057: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1058: switch (dim) {
1059: case 1:
1060: Area = (gmax[0] - gmin[0]);
1061: break;
1062: case 2:
1063: if (user->fake_1D) {
1064: Area = (gmax[0] - gmin[0]);
1065: } else {
1066: Area = (gmax[0] - gmin[0]) * (gmax[1] - gmin[1]);
1067: }
1068: break;
1069: case 3:
1070: Area = (gmax[0] - gmin[0]) * (gmax[1] - gmin[1]) * (gmax[2] - gmin[2]);
1071: break;
1072: default:
1073: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", dim);
1074: }
1075: PetscCallMPI(MPIU_Allreduce(&totalWeight, &global_weight, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
1076: PetscCallMPI(MPIU_Allreduce(&totalCharge, &global_charge, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
1077: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "dim = %" PetscInt_FMT "\ttotalWeight = %f, user->charges[species[0]] = %f\ttotalCharge = %f, Total Area = %f\n", dim, (double)global_weight, (double)user->charges[0], (double)global_charge, (double)Area));
1078: param->sigma = PetscAbsReal(global_charge / (Area));
1080: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sigma: %g\n", (double)param->sigma));
1081: 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,
1082: (double)param->vlasovNumber));
1083: }
1084: /* Setup Constants */
1085: {
1086: PetscDS ds;
1087: Parameter *param;
1088: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1089: PetscScalar constants[NUM_CONSTANTS];
1090: constants[SIGMA] = param->sigma;
1091: constants[V0] = param->v0;
1092: constants[T0] = param->t0;
1093: constants[X0] = param->x0;
1094: constants[M0] = param->m0;
1095: constants[Q0] = param->q0;
1096: constants[PHI0] = param->phi0;
1097: constants[POISSON] = param->poissonNumber;
1098: constants[VLASOV] = param->vlasovNumber;
1099: PetscCall(DMGetDS(dm, &ds));
1100: PetscCall(PetscDSSetConstants(ds, NUM_CONSTANTS, constants));
1101: }
1102: PetscFunctionReturn(PETSC_SUCCESS);
1103: }
1105: static PetscErrorCode InitializeVelocities_Fake1D(DM sw, AppCtx *user)
1106: {
1107: DM dm;
1108: PetscReal *v;
1109: PetscInt *species, cStart, cEnd;
1110: PetscInt dim, Np;
1112: PetscFunctionBegin;
1113: PetscCall(DMGetDimension(sw, &dim));
1114: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1115: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
1116: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1117: PetscCall(DMSwarmGetCellDM(sw, &dm));
1118: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1119: PetscRandom rnd;
1120: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)sw), &rnd));
1121: PetscCall(PetscRandomSetInterval(rnd, 0, 1.));
1122: PetscCall(PetscRandomSetFromOptions(rnd));
1124: for (PetscInt p = 0; p < Np; ++p) {
1125: PetscReal a[3] = {0., 0., 0.}, vel[3] = {0., 0., 0.};
1127: PetscCall(PetscRandomGetValueReal(rnd, &a[0]));
1128: if (user->perturbed_weights) {
1129: PetscCall(PetscPDFSampleConstant1D(a, NULL, vel));
1130: } else {
1131: PetscCall(PetscPDFSampleGaussian1D(a, NULL, vel));
1132: }
1133: v[p * dim] = vel[0];
1134: }
1135: PetscCall(PetscRandomDestroy(&rnd));
1136: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
1137: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1138: PetscFunctionReturn(PETSC_SUCCESS);
1139: }
1141: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
1142: {
1143: PetscReal v0[2] = {1., 0.};
1144: PetscInt dim;
1146: PetscFunctionBeginUser;
1147: PetscCall(DMGetDimension(dm, &dim));
1148: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
1149: PetscCall(DMSetType(*sw, DMSWARM));
1150: PetscCall(DMSetDimension(*sw, dim));
1151: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
1152: PetscCall(DMSwarmSetCellDM(*sw, dm));
1153: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
1154: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
1155: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT));
1156: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "initCoordinates", dim, PETSC_REAL));
1157: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "initVelocity", dim, PETSC_REAL));
1158: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "E_field", dim, PETSC_REAL));
1159: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "potential", dim, PETSC_REAL));
1160: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "charges", dim, PETSC_REAL));
1161: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
1162: PetscCall(DMSetApplicationContext(*sw, user));
1163: PetscCall(DMSetFromOptions(*sw));
1164: user->swarm = *sw;
1165: if (user->perturbed_weights) {
1166: PetscCall(InitializeParticles_PerturbedWeights(*sw, user));
1167: } else {
1168: PetscCall(DMSwarmComputeLocalSizeFromOptions(*sw));
1169: PetscCall(DMSwarmInitializeCoordinates(*sw));
1170: if (user->fake_1D) {
1171: PetscCall(InitializeVelocities_Fake1D(*sw, user));
1172: } else {
1173: PetscCall(DMSwarmInitializeVelocitiesFromOptions(*sw, v0));
1174: }
1175: }
1176: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
1177: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
1178: {
1179: Vec gc, gc0, gv, gv0;
1181: PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc));
1182: PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "initCoordinates", &gc0));
1183: PetscCall(VecCopy(gc, gc0));
1184: PetscCall(VecViewFromOptions(gc, NULL, "-ic_x_view"));
1185: PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc));
1186: PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "initCoordinates", &gc0));
1187: PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "velocity", &gv));
1188: PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "initVelocity", &gv0));
1189: PetscCall(VecCopy(gv, gv0));
1190: PetscCall(VecViewFromOptions(gv, NULL, "-ic_v_view"));
1191: PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "velocity", &gv));
1192: PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "initVelocity", &gv0));
1193: }
1194: PetscFunctionReturn(PETSC_SUCCESS);
1195: }
1197: static PetscErrorCode ComputeFieldAtParticles_Coulomb(SNES snes, DM sw, PetscReal E[])
1198: {
1199: AppCtx *user;
1200: PetscReal *coords;
1201: PetscInt *species, dim, Np, Ns;
1202: PetscMPIInt size;
1204: PetscFunctionBegin;
1205: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)snes), &size));
1206: PetscCheck(size == 1, PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Coulomb code only works in serial");
1207: PetscCall(DMGetDimension(sw, &dim));
1208: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1209: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1210: PetscCall(DMGetApplicationContext(sw, (void *)&user));
1212: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1213: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1214: for (PetscInt p = 0; p < Np; ++p) {
1215: PetscReal *pcoord = &coords[p * dim];
1216: PetscReal pE[3] = {0., 0., 0.};
1218: /* Calculate field at particle p due to particle q */
1219: for (PetscInt q = 0; q < Np; ++q) {
1220: PetscReal *qcoord = &coords[q * dim];
1221: PetscReal rpq[3], r, r3, q_q;
1223: if (p == q) continue;
1224: q_q = user->charges[species[q]] * 1.;
1225: for (PetscInt d = 0; d < dim; ++d) rpq[d] = pcoord[d] - qcoord[d];
1226: r = DMPlex_NormD_Internal(dim, rpq);
1227: if (r < PETSC_SQRT_MACHINE_EPSILON) continue;
1228: r3 = PetscPowRealInt(r, 3);
1229: for (PetscInt d = 0; d < dim; ++d) pE[d] += q_q * rpq[d] / r3;
1230: }
1231: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] = pE[d];
1232: }
1233: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1234: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1235: PetscFunctionReturn(PETSC_SUCCESS);
1236: }
1238: static PetscErrorCode ComputeFieldAtParticles_Primal(SNES snes, DM sw, PetscReal E[])
1239: {
1240: DM dm;
1241: AppCtx *user;
1242: PetscDS ds;
1243: PetscFE fe;
1244: Mat M_p, M;
1245: Vec phi, locPhi, rho, f;
1246: PetscReal *coords;
1247: PetscInt dim, cStart, cEnd, Np;
1248: PetscQuadrature q;
1250: PetscFunctionBegin;
1251: PetscCall(DMGetDimension(sw, &dim));
1252: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1253: PetscCall(DMGetApplicationContext(sw, (void *)&user));
1255: KSP ksp;
1256: Vec rho0;
1257: char oldField[PETSC_MAX_PATH_LEN];
1258: const char *tmp;
1260: /* Create the charges rho */
1261: PetscCall(SNESGetDM(snes, &dm));
1262: PetscCall(DMSwarmVectorGetField(sw, &tmp));
1263: PetscCall(PetscStrncpy(oldField, tmp, PETSC_MAX_PATH_LEN));
1264: PetscCall(DMSwarmVectorDefineField(sw, "w_q"));
1265: PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
1266: PetscCall(DMSwarmVectorDefineField(sw, oldField));
1268: PetscCall(DMCreateMassMatrix(dm, dm, &M));
1269: PetscCall(DMGetGlobalVector(dm, &rho0));
1270: PetscCall(PetscObjectSetName((PetscObject)rho0, "Charge density (rho0) from Primal Compute"));
1271: PetscCall(DMGetGlobalVector(dm, &rho));
1272: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
1273: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
1275: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
1276: PetscCall(MatMultTranspose(M_p, f, rho));
1277: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
1278: PetscCall(MatViewFromOptions(M, NULL, "-m_view"));
1279: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
1280: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
1282: PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &ksp));
1283: PetscCall(KSPSetOptionsPrefix(ksp, "em_proj_"));
1284: PetscCall(KSPSetOperators(ksp, M, M));
1285: PetscCall(KSPSetFromOptions(ksp));
1286: PetscCall(KSPSolve(ksp, rho, rho0));
1287: PetscCall(VecViewFromOptions(rho0, NULL, "-rho0_view"));
1289: PetscInt rhosize;
1290: PetscReal *charges;
1291: const PetscScalar *rho_vals;
1292: PetscCall(DMSwarmGetField(sw, "charges", NULL, NULL, (void **)&charges));
1293: PetscCall(VecGetLocalSize(rho0, &rhosize));
1294: PetscCall(VecGetArrayRead(rho0, &rho_vals));
1295: for (PetscInt c = 0; c < rhosize; ++c) charges[c] = rho_vals[c];
1296: PetscCall(VecRestoreArrayRead(rho0, &rho_vals));
1297: PetscCall(DMSwarmRestoreField(sw, "charges", NULL, NULL, (void **)&charges));
1299: PetscCall(VecScale(rho, -1.0));
1301: PetscCall(VecViewFromOptions(rho0, NULL, "-rho0_view"));
1302: PetscCall(VecViewFromOptions(rho, NULL, "-rho_view"));
1303: PetscCall(DMRestoreGlobalVector(dm, &rho0));
1304: PetscCall(KSPDestroy(&ksp));
1305: PetscCall(MatDestroy(&M_p));
1306: PetscCall(MatDestroy(&M));
1308: PetscCall(DMGetGlobalVector(dm, &phi));
1309: PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
1310: PetscCall(VecSet(phi, 0.0));
1311: PetscCall(SNESSolve(snes, rho, phi));
1312: PetscCall(DMRestoreGlobalVector(dm, &rho));
1313: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
1315: PetscInt phisize;
1316: PetscReal *pot;
1317: const PetscScalar *phi_vals;
1318: PetscCall(DMSwarmGetField(sw, "potential", NULL, NULL, (void **)&pot));
1319: PetscCall(VecGetLocalSize(phi, &phisize));
1320: PetscCall(VecGetArrayRead(phi, &phi_vals));
1321: for (PetscInt c = 0; c < phisize; ++c) pot[c] = phi_vals[c];
1322: PetscCall(VecRestoreArrayRead(phi, &phi_vals));
1323: PetscCall(DMSwarmRestoreField(sw, "potential", NULL, NULL, (void **)&pot));
1325: PetscCall(DMGetLocalVector(dm, &locPhi));
1326: PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
1327: PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
1328: PetscCall(DMRestoreGlobalVector(dm, &phi));
1330: PetscCall(DMGetDS(dm, &ds));
1331: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
1332: PetscCall(DMSwarmSortGetAccess(sw));
1333: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1334: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1336: for (PetscInt c = cStart; c < cEnd; ++c) {
1337: PetscTabulation tab;
1338: PetscScalar *clPhi = NULL;
1339: PetscReal *pcoord, *refcoord;
1340: PetscReal v[3], J[9], invJ[9], detJ;
1341: PetscInt *points;
1342: PetscInt Ncp;
1344: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
1345: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
1346: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
1347: for (PetscInt cp = 0; cp < Ncp; ++cp)
1348: for (PetscInt d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
1349: PetscCall(DMPlexCoordinatesToReference(dm, c, Ncp, pcoord, refcoord));
1350: PetscCall(PetscFECreateTabulation(fe, 1, Ncp, refcoord, 1, &tab));
1351: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v, J, invJ, &detJ));
1352: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
1353: for (PetscInt cp = 0; cp < Ncp; ++cp) {
1354: const PetscReal *basisDer = tab->T[1];
1355: const PetscInt p = points[cp];
1357: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] = 0.;
1358: PetscCall(PetscFEGetQuadrature(fe, &q));
1359: PetscCall(PetscFEFreeInterpolateGradient_Static(fe, basisDer, clPhi, dim, invJ, NULL, cp, &E[p * dim]));
1360: for (PetscInt d = 0; d < dim; ++d) {
1361: E[p * dim + d] *= -1.0;
1362: if (user->fake_1D && d > 0) E[p * dim + d] = 0;
1363: }
1364: }
1365: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
1366: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
1367: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
1368: PetscCall(PetscTabulationDestroy(&tab));
1369: PetscCall(PetscFree(points));
1370: }
1371: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1372: PetscCall(DMSwarmSortRestoreAccess(sw));
1373: PetscCall(DMRestoreLocalVector(dm, &locPhi));
1374: PetscFunctionReturn(PETSC_SUCCESS);
1375: }
1377: static PetscErrorCode ComputeFieldAtParticles_Mixed(SNES snes, DM sw, PetscReal E[])
1378: {
1379: AppCtx *user;
1380: DM dm, potential_dm;
1381: KSP ksp;
1382: IS potential_IS;
1383: PetscDS ds;
1384: PetscFE fe;
1385: PetscFEGeom feGeometry;
1386: Mat M_p, M;
1387: Vec phi, locPhi, rho, f, temp_rho, rho0;
1388: PetscQuadrature q;
1389: PetscReal *coords, *pot;
1390: PetscInt dim, cStart, cEnd, Np, fields = 1;
1391: char oldField[PETSC_MAX_PATH_LEN];
1392: const char *tmp;
1394: PetscFunctionBegin;
1395: PetscCall(DMGetDimension(sw, &dim));
1396: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1397: PetscCall(DMGetApplicationContext(sw, &user));
1399: /* Create the charges rho */
1400: PetscCall(SNESGetDM(snes, &dm));
1401: PetscCall(DMGetGlobalVector(dm, &rho));
1402: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
1404: PetscCall(DMCreateSubDM(dm, 1, &fields, &potential_IS, &potential_dm));
1406: PetscCall(DMSwarmVectorGetField(sw, &tmp));
1407: PetscCall(PetscStrncpy(oldField, tmp, PETSC_MAX_PATH_LEN));
1408: PetscCall(DMSwarmVectorDefineField(sw, "w_q"));
1409: PetscCall(DMCreateMassMatrix(sw, potential_dm, &M_p));
1410: PetscCall(DMSwarmVectorDefineField(sw, oldField));
1412: PetscCall(DMCreateMassMatrix(potential_dm, potential_dm, &M));
1413: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
1414: PetscCall(MatViewFromOptions(M, NULL, "-m_view"));
1415: PetscCall(DMGetGlobalVector(potential_dm, &temp_rho));
1416: PetscCall(PetscObjectSetName((PetscObject)temp_rho, "Mf"));
1417: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
1418: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
1419: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
1420: PetscCall(MatMultTranspose(M_p, f, temp_rho));
1421: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
1422: PetscCall(DMGetGlobalVector(potential_dm, &rho0));
1423: PetscCall(PetscObjectSetName((PetscObject)rho0, "Charge density (rho0) from Mixed Compute"));
1425: PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &ksp));
1426: PetscCall(KSPSetOptionsPrefix(ksp, "em_proj"));
1427: PetscCall(KSPSetOperators(ksp, M, M));
1428: PetscCall(KSPSetFromOptions(ksp));
1429: PetscCall(KSPSolve(ksp, temp_rho, rho0));
1430: PetscCall(VecViewFromOptions(rho0, NULL, "-rho0_view"));
1432: PetscInt rhosize;
1433: PetscReal *charges;
1434: const PetscScalar *rho_vals;
1435: Parameter *param;
1436: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1437: PetscCall(DMSwarmGetField(sw, "charges", NULL, NULL, (void **)&charges));
1438: PetscCall(VecGetLocalSize(rho0, &rhosize));
1440: /* Integral over reference element is size 1. Reference element area is 4. Scale rho0 by 1/4 because the basis function is 1/4 */
1441: PetscCall(VecScale(rho0, 0.25));
1442: PetscCall(VecGetArrayRead(rho0, &rho_vals));
1443: for (PetscInt c = 0; c < rhosize; ++c) charges[c] = rho_vals[c];
1444: PetscCall(VecRestoreArrayRead(rho0, &rho_vals));
1445: PetscCall(DMSwarmRestoreField(sw, "charges", NULL, NULL, (void **)&charges));
1447: PetscCall(VecISCopy(rho, potential_IS, SCATTER_FORWARD, temp_rho));
1448: PetscCall(VecScale(rho, 0.25));
1449: PetscCall(VecViewFromOptions(rho0, NULL, "-rho0_view"));
1450: PetscCall(VecViewFromOptions(temp_rho, NULL, "-temprho_view"));
1451: PetscCall(VecViewFromOptions(rho, NULL, "-rho_view"));
1452: PetscCall(DMRestoreGlobalVector(potential_dm, &temp_rho));
1453: PetscCall(DMRestoreGlobalVector(potential_dm, &rho0));
1455: PetscCall(MatDestroy(&M_p));
1456: PetscCall(MatDestroy(&M));
1457: PetscCall(KSPDestroy(&ksp));
1458: PetscCall(DMDestroy(&potential_dm));
1459: PetscCall(ISDestroy(&potential_IS));
1461: PetscCall(DMGetGlobalVector(dm, &phi));
1462: PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
1463: PetscCall(VecSet(phi, 0.0));
1464: PetscCall(SNESSolve(snes, rho, phi));
1465: PetscCall(DMRestoreGlobalVector(dm, &rho));
1467: PetscInt phisize;
1468: const PetscScalar *phi_vals;
1469: PetscCall(DMSwarmGetField(sw, "potential", NULL, NULL, (void **)&pot));
1470: PetscCall(VecGetLocalSize(phi, &phisize));
1471: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
1472: PetscCall(VecGetArrayRead(phi, &phi_vals));
1473: for (PetscInt c = 0; c < phisize; ++c) pot[c] = phi_vals[c];
1474: PetscCall(VecRestoreArrayRead(phi, &phi_vals));
1475: PetscCall(DMSwarmRestoreField(sw, "potential", NULL, NULL, (void **)&pot));
1477: PetscCall(DMGetLocalVector(dm, &locPhi));
1478: PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
1479: PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
1480: PetscCall(DMRestoreGlobalVector(dm, &phi));
1482: PetscCall(DMGetDS(dm, &ds));
1483: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
1484: PetscCall(DMSwarmSortGetAccess(sw));
1485: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1486: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1487: PetscCall(PetscFEGetQuadrature(fe, &q));
1488: PetscCall(PetscFECreateCellGeometry(fe, q, &feGeometry));
1489: for (PetscInt c = cStart; c < cEnd; ++c) {
1490: PetscTabulation tab;
1491: PetscScalar *clPhi = NULL;
1492: PetscReal *pcoord, *refcoord;
1493: PetscInt *points;
1494: PetscInt Ncp;
1496: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
1497: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
1498: PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
1499: for (PetscInt cp = 0; cp < Ncp; ++cp)
1500: for (PetscInt d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
1501: PetscCall(DMPlexCoordinatesToReference(dm, c, Ncp, pcoord, refcoord));
1502: PetscCall(PetscFECreateTabulation(fe, 1, Ncp, refcoord, 1, &tab));
1503: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, q, feGeometry.v, feGeometry.J, feGeometry.invJ, feGeometry.detJ));
1504: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
1506: for (PetscInt cp = 0; cp < Ncp; ++cp) {
1507: const PetscInt p = points[cp];
1509: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] = 0.;
1510: PetscCall(PetscFEInterpolateAtPoints_Static(fe, tab, clPhi, &feGeometry, cp, &E[p * dim]));
1511: PetscCall(PetscFEPushforward(fe, &feGeometry, 1, &E[p * dim]));
1512: for (PetscInt d = 0; d < dim; ++d) {
1513: E[p * dim + d] *= -2.0;
1514: if (user->fake_1D && d > 0) E[p * dim + d] = 0;
1515: }
1516: }
1517: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
1518: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
1519: PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
1520: PetscCall(PetscTabulationDestroy(&tab));
1521: PetscCall(PetscFree(points));
1522: }
1523: PetscCall(PetscFEDestroyCellGeometry(fe, &feGeometry));
1524: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1525: PetscCall(DMSwarmSortRestoreAccess(sw));
1526: PetscCall(DMRestoreLocalVector(dm, &locPhi));
1527: PetscFunctionReturn(PETSC_SUCCESS);
1528: }
1530: static PetscErrorCode ComputeFieldAtParticles(SNES snes, DM sw, PetscReal E[])
1531: {
1532: AppCtx *ctx;
1533: PetscInt dim, Np;
1535: PetscFunctionBegin;
1538: PetscAssertPointer(E, 3);
1539: PetscCall(DMGetDimension(sw, &dim));
1540: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1541: PetscCall(DMGetApplicationContext(sw, &ctx));
1542: PetscCall(PetscArrayzero(E, Np * dim));
1544: switch (ctx->em) {
1545: case EM_PRIMAL:
1546: PetscCall(ComputeFieldAtParticles_Primal(snes, sw, E));
1547: break;
1548: case EM_COULOMB:
1549: PetscCall(ComputeFieldAtParticles_Coulomb(snes, sw, E));
1550: break;
1551: case EM_MIXED:
1552: PetscCall(ComputeFieldAtParticles_Mixed(snes, sw, E));
1553: break;
1554: case EM_NONE:
1555: break;
1556: default:
1557: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No solver for electrostatic model %s", EMTypes[ctx->em]);
1558: }
1559: PetscFunctionReturn(PETSC_SUCCESS);
1560: }
1562: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx)
1563: {
1564: DM sw;
1565: SNES snes = ((AppCtx *)ctx)->snes;
1566: const PetscReal *coords, *vel;
1567: const PetscScalar *u;
1568: PetscScalar *g;
1569: PetscReal *E, m_p = 1., q_p = -1.;
1570: PetscInt dim, d, Np, p;
1572: PetscFunctionBeginUser;
1573: PetscCall(TSGetDM(ts, &sw));
1574: PetscCall(DMGetDimension(sw, &dim));
1575: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
1576: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
1577: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
1578: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1579: PetscCall(VecGetArrayRead(U, &u));
1580: PetscCall(VecGetArray(G, &g));
1582: PetscCall(ComputeFieldAtParticles(snes, sw, E));
1584: Np /= 2 * dim;
1585: for (p = 0; p < Np; ++p) {
1586: for (d = 0; d < dim; ++d) {
1587: g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d];
1588: g[(p * 2 + 1) * dim + d] = q_p * E[p * dim + d] / m_p;
1589: }
1590: }
1591: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
1592: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
1593: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
1594: PetscCall(VecRestoreArrayRead(U, &u));
1595: PetscCall(VecRestoreArray(G, &g));
1596: PetscFunctionReturn(PETSC_SUCCESS);
1597: }
1599: /* J_{ij} = dF_i/dx_j
1600: J_p = ( 0 1)
1601: (-w^2 0)
1602: TODO Now there is another term with w^2 from the electric field. I think we will need to invert the operator.
1603: Perhaps we can approximate the Jacobian using only the cellwise P-P gradient from Coulomb
1604: */
1605: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, void *ctx)
1606: {
1607: DM sw;
1608: const PetscReal *coords, *vel;
1609: PetscInt dim, d, Np, p, rStart;
1611: PetscFunctionBeginUser;
1612: PetscCall(TSGetDM(ts, &sw));
1613: PetscCall(DMGetDimension(sw, &dim));
1614: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1615: PetscCall(MatGetOwnershipRange(J, &rStart, NULL));
1616: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
1617: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
1618: Np /= 2 * dim;
1619: for (p = 0; p < Np; ++p) {
1620: const PetscReal x0 = coords[p * dim + 0];
1621: const PetscReal vy0 = vel[p * dim + 1];
1622: const PetscReal omega = vy0 / x0;
1623: PetscScalar vals[4] = {0., 1., -PetscSqr(omega), 0.};
1625: for (d = 0; d < dim; ++d) {
1626: const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
1627: PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES));
1628: }
1629: }
1630: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
1631: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
1632: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1633: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1634: PetscFunctionReturn(PETSC_SUCCESS);
1635: }
1637: static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
1638: {
1639: AppCtx *user = (AppCtx *)ctx;
1640: DM sw;
1641: const PetscScalar *v;
1642: PetscScalar *xres;
1643: PetscInt Np, p, d, dim;
1645: PetscFunctionBeginUser;
1646: PetscCall(TSGetDM(ts, &sw));
1647: PetscCall(DMGetDimension(sw, &dim));
1648: PetscCall(VecGetLocalSize(Xres, &Np));
1649: PetscCall(VecGetArrayRead(V, &v));
1650: PetscCall(VecGetArray(Xres, &xres));
1651: Np /= dim;
1652: for (p = 0; p < Np; ++p) {
1653: for (d = 0; d < dim; ++d) {
1654: xres[p * dim + d] = v[p * dim + d];
1655: if (user->fake_1D && d > 0) xres[p * dim + d] = 0;
1656: }
1657: }
1658: PetscCall(VecRestoreArrayRead(V, &v));
1659: PetscCall(VecRestoreArray(Xres, &xres));
1660: PetscFunctionReturn(PETSC_SUCCESS);
1661: }
1663: static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
1664: {
1665: DM sw;
1666: AppCtx *user = (AppCtx *)ctx;
1667: SNES snes = ((AppCtx *)ctx)->snes;
1668: const PetscScalar *x;
1669: const PetscReal *coords, *vel;
1670: PetscReal *E, m_p, q_p;
1671: PetscScalar *vres;
1672: PetscInt Np, p, dim, d;
1673: Parameter *param;
1675: PetscFunctionBeginUser;
1676: PetscCall(TSGetDM(ts, &sw));
1677: PetscCall(DMGetDimension(sw, &dim));
1678: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
1679: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
1680: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
1681: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1682: m_p = user->masses[0] * param->m0;
1683: q_p = user->charges[0] * param->q0;
1684: PetscCall(VecGetLocalSize(Vres, &Np));
1685: PetscCall(VecGetArrayRead(X, &x));
1686: PetscCall(VecGetArray(Vres, &vres));
1687: PetscCheck(dim == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension must be 2");
1688: PetscCall(ComputeFieldAtParticles(snes, sw, E));
1690: Np /= dim;
1691: for (p = 0; p < Np; ++p) {
1692: for (d = 0; d < dim; ++d) {
1693: vres[p * dim + d] = q_p * E[p * dim + d] / m_p;
1694: if (user->fake_1D && d > 0) vres[p * dim + d] = 0.;
1695: }
1696: }
1697: PetscCall(VecRestoreArrayRead(X, &x));
1698: /*
1699: Synchronized, ordered output for parallel/sequential test cases.
1700: In the 1D (on the 2D mesh) case, every y component should be zero.
1701: */
1702: if (user->checkVRes) {
1703: PetscBool pr = user->checkVRes > 1 ? PETSC_TRUE : PETSC_FALSE;
1704: PetscInt step;
1706: PetscCall(TSGetStepNumber(ts, &step));
1707: if (pr) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "step: %" PetscInt_FMT "\n", step));
1708: for (PetscInt p = 0; p < Np; ++p) {
1709: if (pr) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Residual: %.12g %.12g\n", (double)PetscRealPart(vres[p * dim + 0]), (double)PetscRealPart(vres[p * dim + 1])));
1710: 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]));
1711: }
1712: if (pr) PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
1713: }
1714: PetscCall(VecRestoreArray(Vres, &vres));
1715: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
1716: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
1717: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
1718: PetscFunctionReturn(PETSC_SUCCESS);
1719: }
1721: static PetscErrorCode CreateSolution(TS ts)
1722: {
1723: DM sw;
1724: Vec u;
1725: PetscInt dim, Np;
1727: PetscFunctionBegin;
1728: PetscCall(TSGetDM(ts, &sw));
1729: PetscCall(DMGetDimension(sw, &dim));
1730: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1731: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
1732: PetscCall(VecSetBlockSize(u, dim));
1733: PetscCall(VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE));
1734: PetscCall(VecSetUp(u));
1735: PetscCall(TSSetSolution(ts, u));
1736: PetscCall(VecDestroy(&u));
1737: PetscFunctionReturn(PETSC_SUCCESS);
1738: }
1740: static PetscErrorCode SetProblem(TS ts)
1741: {
1742: AppCtx *user;
1743: DM sw;
1745: PetscFunctionBegin;
1746: PetscCall(TSGetDM(ts, &sw));
1747: PetscCall(DMGetApplicationContext(sw, (void **)&user));
1748: // Define unified system for (X, V)
1749: {
1750: Mat J;
1751: PetscInt dim, Np;
1753: PetscCall(DMGetDimension(sw, &dim));
1754: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1755: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
1756: PetscCall(MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE));
1757: PetscCall(MatSetBlockSize(J, 2 * dim));
1758: PetscCall(MatSetFromOptions(J));
1759: PetscCall(MatSetUp(J));
1760: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, user));
1761: PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, user));
1762: PetscCall(MatDestroy(&J));
1763: }
1764: /* Define split system for X and V */
1765: {
1766: Vec u;
1767: IS isx, isv, istmp;
1768: const PetscInt *idx;
1769: PetscInt dim, Np, rstart;
1771: PetscCall(TSGetSolution(ts, &u));
1772: PetscCall(DMGetDimension(sw, &dim));
1773: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1774: PetscCall(VecGetOwnershipRange(u, &rstart, NULL));
1775: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 0, 2, &istmp));
1776: PetscCall(ISGetIndices(istmp, &idx));
1777: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isx));
1778: PetscCall(ISRestoreIndices(istmp, &idx));
1779: PetscCall(ISDestroy(&istmp));
1780: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 1, 2, &istmp));
1781: PetscCall(ISGetIndices(istmp, &idx));
1782: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isv));
1783: PetscCall(ISRestoreIndices(istmp, &idx));
1784: PetscCall(ISDestroy(&istmp));
1785: PetscCall(TSRHSSplitSetIS(ts, "position", isx));
1786: PetscCall(TSRHSSplitSetIS(ts, "momentum", isv));
1787: PetscCall(ISDestroy(&isx));
1788: PetscCall(ISDestroy(&isv));
1789: PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, user));
1790: PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, user));
1791: }
1792: PetscFunctionReturn(PETSC_SUCCESS);
1793: }
1795: static PetscErrorCode DMSwarmTSRedistribute(TS ts)
1796: {
1797: DM sw;
1798: Vec u;
1799: PetscReal t, maxt, dt;
1800: PetscInt n, maxn;
1802: PetscFunctionBegin;
1803: PetscCall(TSGetDM(ts, &sw));
1804: PetscCall(TSGetTime(ts, &t));
1805: PetscCall(TSGetMaxTime(ts, &maxt));
1806: PetscCall(TSGetTimeStep(ts, &dt));
1807: PetscCall(TSGetStepNumber(ts, &n));
1808: PetscCall(TSGetMaxSteps(ts, &maxn));
1810: PetscCall(TSReset(ts));
1811: PetscCall(TSSetDM(ts, sw));
1812: PetscCall(TSSetFromOptions(ts));
1813: PetscCall(TSSetTime(ts, t));
1814: PetscCall(TSSetMaxTime(ts, maxt));
1815: PetscCall(TSSetTimeStep(ts, dt));
1816: PetscCall(TSSetStepNumber(ts, n));
1817: PetscCall(TSSetMaxSteps(ts, maxn));
1819: PetscCall(CreateSolution(ts));
1820: PetscCall(SetProblem(ts));
1821: PetscCall(TSGetSolution(ts, &u));
1822: PetscFunctionReturn(PETSC_SUCCESS);
1823: }
1825: PetscErrorCode line(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
1826: {
1827: DM sw, cdm;
1828: PetscInt Np;
1829: PetscReal low[2], high[2];
1830: AppCtx *user = (AppCtx *)ctx;
1832: sw = user->swarm;
1833: PetscCall(DMSwarmGetCellDM(sw, &cdm));
1834: // Get the bounding box so we can equally space the particles
1835: PetscCall(DMGetLocalBoundingBox(cdm, low, high));
1836: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1837: // shift it by h/2 so nothing is initialized directly on a boundary
1838: x[0] = ((high[0] - low[0]) / Np) * (p + 0.5);
1839: x[1] = 0.;
1840: return PETSC_SUCCESS;
1841: }
1843: /*
1844: InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values.
1846: Input Parameters:
1847: + ts - The TS
1848: - useInitial - Flag to also set the initial conditions to the current coordinates and velocities and setup the problem
1850: Output Parameters:
1851: . u - The initialized solution vector
1853: Level: advanced
1855: .seealso: InitializeSolve()
1856: */
1857: static PetscErrorCode InitializeSolveAndSwarm(TS ts, PetscBool useInitial)
1858: {
1859: DM sw;
1860: Vec u, gc, gv, gc0, gv0;
1861: IS isx, isv;
1862: PetscInt dim;
1863: AppCtx *user;
1865: PetscFunctionBeginUser;
1866: PetscCall(TSGetDM(ts, &sw));
1867: PetscCall(DMGetApplicationContext(sw, &user));
1868: PetscCall(DMGetDimension(sw, &dim));
1869: if (useInitial) {
1870: PetscReal v0[2] = {1., 0.};
1871: if (user->perturbed_weights) {
1872: PetscCall(InitializeParticles_PerturbedWeights(sw, user));
1873: } else {
1874: PetscCall(DMSwarmComputeLocalSizeFromOptions(sw));
1875: PetscCall(DMSwarmInitializeCoordinates(sw));
1876: if (user->fake_1D) {
1877: PetscCall(InitializeVelocities_Fake1D(sw, user));
1878: } else {
1879: PetscCall(DMSwarmInitializeVelocitiesFromOptions(sw, v0));
1880: }
1881: }
1882: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
1883: PetscCall(DMSwarmTSRedistribute(ts));
1884: }
1885: PetscCall(TSGetSolution(ts, &u));
1886: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
1887: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
1888: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
1889: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initCoordinates", &gc0));
1890: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
1891: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initVelocity", &gv0));
1892: if (useInitial) {
1893: PetscCall(VecCopy(gc, gc0));
1894: PetscCall(VecCopy(gv, gv0));
1895: }
1896: PetscCall(VecISCopy(u, isx, SCATTER_FORWARD, gc));
1897: PetscCall(VecISCopy(u, isv, SCATTER_FORWARD, gv));
1898: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
1899: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initCoordinates", &gc0));
1900: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
1901: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initVelocity", &gv0));
1902: PetscFunctionReturn(PETSC_SUCCESS);
1903: }
1905: static PetscErrorCode InitializeSolve(TS ts, Vec u)
1906: {
1907: PetscFunctionBegin;
1908: PetscCall(TSSetSolution(ts, u));
1909: PetscCall(InitializeSolveAndSwarm(ts, PETSC_TRUE));
1910: PetscFunctionReturn(PETSC_SUCCESS);
1911: }
1913: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
1914: {
1915: MPI_Comm comm;
1916: DM sw;
1917: AppCtx *user;
1918: const PetscScalar *u;
1919: const PetscReal *coords, *vel;
1920: PetscScalar *e;
1921: PetscReal t;
1922: PetscInt dim, Np, p;
1924: PetscFunctionBeginUser;
1925: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
1926: PetscCall(TSGetDM(ts, &sw));
1927: PetscCall(DMGetApplicationContext(sw, &user));
1928: PetscCall(DMGetDimension(sw, &dim));
1929: PetscCall(TSGetSolveTime(ts, &t));
1930: PetscCall(VecGetArray(E, &e));
1931: PetscCall(VecGetArrayRead(U, &u));
1932: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1933: PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
1934: PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
1935: Np /= 2 * dim;
1936: for (p = 0; p < Np; ++p) {
1937: /* TODO generalize initial conditions and project into plane instead of assuming x-y */
1938: const PetscReal r0 = DMPlex_NormD_Internal(dim, &coords[p * dim]);
1939: const PetscReal th0 = PetscAtan2Real(coords[p * dim + 1], coords[p * dim + 0]);
1940: const PetscReal v0 = DMPlex_NormD_Internal(dim, &vel[p * dim]);
1941: const PetscReal omega = v0 / r0;
1942: const PetscReal ct = PetscCosReal(omega * t + th0);
1943: const PetscReal st = PetscSinReal(omega * t + th0);
1944: const PetscScalar *x = &u[(p * 2 + 0) * dim];
1945: const PetscScalar *v = &u[(p * 2 + 1) * dim];
1946: const PetscReal xe[3] = {r0 * ct, r0 * st, 0.0};
1947: const PetscReal ve[3] = {-v0 * st, v0 * ct, 0.0};
1948: PetscInt d;
1950: for (d = 0; d < dim; ++d) {
1951: e[(p * 2 + 0) * dim + d] = x[d] - xe[d];
1952: e[(p * 2 + 1) * dim + d] = v[d] - ve[d];
1953: }
1954: if (user->error) {
1955: const PetscReal en = 0.5 * DMPlex_DotRealD_Internal(dim, v, v);
1956: const PetscReal exen = 0.5 * PetscSqr(v0);
1957: PetscCall(PetscPrintf(comm, "t %.4g: p%" PetscInt_FMT " error [%.2g %.2g] sol [(%.6lf %.6lf) (%.6lf %.6lf)] exact [(%.6lf %.6lf) (%.6lf %.6lf)] energy/exact energy %g / %g (%.10lf%%)\n", (double)t, p, (double)DMPlex_NormD_Internal(dim, &e[(p * 2 + 0) * dim]), (double)DMPlex_NormD_Internal(dim, &e[(p * 2 + 1) * dim]), (double)x[0], (double)x[1], (double)v[0], (double)v[1], (double)xe[0], (double)xe[1], (double)ve[0], (double)ve[1], (double)en, (double)exen, (double)(PetscAbsReal(exen - en) * 100. / exen)));
1958: }
1959: }
1960: PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
1961: PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
1962: PetscCall(VecRestoreArrayRead(U, &u));
1963: PetscCall(VecRestoreArray(E, &e));
1964: PetscFunctionReturn(PETSC_SUCCESS);
1965: }
1967: static PetscErrorCode MigrateParticles(TS ts)
1968: {
1969: DM sw, cdm;
1970: const PetscReal *L;
1972: PetscFunctionBeginUser;
1973: PetscCall(TSGetDM(ts, &sw));
1974: PetscCall(DMViewFromOptions(sw, NULL, "-migrate_view_pre"));
1975: {
1976: Vec u, gc, gv, position, momentum;
1977: IS isx, isv;
1978: PetscReal *pos, *mom;
1980: PetscCall(TSGetSolution(ts, &u));
1981: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
1982: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
1983: PetscCall(VecGetSubVector(u, isx, &position));
1984: PetscCall(VecGetSubVector(u, isv, &momentum));
1985: PetscCall(VecGetArray(position, &pos));
1986: PetscCall(VecGetArray(momentum, &mom));
1987: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
1988: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
1989: PetscCall(VecISCopy(u, isx, SCATTER_REVERSE, gc));
1990: PetscCall(VecISCopy(u, isv, SCATTER_REVERSE, gv));
1992: PetscCall(DMSwarmGetCellDM(sw, &cdm));
1993: PetscCall(DMGetPeriodicity(cdm, NULL, NULL, &L));
1994: if ((L[0] || L[1]) >= 0.) {
1995: PetscReal *x, *v, upper[3], lower[3];
1996: PetscInt Np, dim;
1998: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1999: PetscCall(DMGetDimension(cdm, &dim));
2000: PetscCall(DMGetBoundingBox(cdm, lower, upper));
2001: PetscCall(VecGetArray(gc, &x));
2002: PetscCall(VecGetArray(gv, &v));
2003: for (PetscInt p = 0; p < Np; ++p) {
2004: for (PetscInt d = 0; d < dim; ++d) {
2005: if (pos[p * dim + d] < lower[d]) {
2006: x[p * dim + d] = pos[p * dim + d] + (upper[d] - lower[d]);
2007: } else if (pos[p * dim + d] > upper[d]) {
2008: x[p * dim + d] = pos[p * dim + d] - (upper[d] - lower[d]);
2009: } else {
2010: x[p * dim + d] = pos[p * dim + d];
2011: }
2012: 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]);
2013: v[p * dim + d] = mom[p * dim + d];
2014: }
2015: }
2016: PetscCall(VecRestoreArray(gc, &x));
2017: PetscCall(VecRestoreArray(gv, &v));
2018: }
2019: PetscCall(VecRestoreArray(position, &pos));
2020: PetscCall(VecRestoreArray(momentum, &mom));
2021: PetscCall(VecRestoreSubVector(u, isx, &position));
2022: PetscCall(VecRestoreSubVector(u, isv, &momentum));
2023: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
2024: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
2025: }
2026: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
2027: PetscCall(DMSwarmTSRedistribute(ts));
2028: PetscCall(InitializeSolveAndSwarm(ts, PETSC_FALSE));
2029: PetscFunctionReturn(PETSC_SUCCESS);
2030: }
2032: int main(int argc, char **argv)
2033: {
2034: DM dm, sw;
2035: TS ts;
2036: Vec u;
2037: PetscReal dt;
2038: PetscInt maxn;
2039: AppCtx user;
2041: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
2042: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
2043: PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &user.bag));
2044: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
2045: PetscCall(CreatePoisson(dm, &user));
2046: PetscCall(CreateSwarm(dm, &user, &sw));
2047: PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
2048: PetscCall(InitializeConstants(sw, &user));
2049: PetscCall(DMSetApplicationContext(sw, &user));
2051: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
2052: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
2053: PetscCall(TSSetDM(ts, sw));
2054: PetscCall(TSSetMaxTime(ts, 0.1));
2055: PetscCall(TSSetTimeStep(ts, 0.00001));
2056: PetscCall(TSSetMaxSteps(ts, 100));
2057: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
2059: if (user.efield_monitor) PetscCall(TSMonitorSet(ts, MonitorEField, &user, NULL));
2060: if (user.initial_monitor) PetscCall(TSMonitorSet(ts, MonitorInitialConditions, &user, NULL));
2061: if (user.monitor_positions) PetscCall(TSMonitorSet(ts, MonitorPositions_2D, &user, NULL));
2062: if (user.poisson_monitor) PetscCall(TSMonitorSet(ts, MonitorPoisson, &user, NULL));
2064: PetscCall(TSSetFromOptions(ts));
2065: PetscCall(TSGetTimeStep(ts, &dt));
2066: PetscCall(TSGetMaxSteps(ts, &maxn));
2067: user.steps = maxn;
2068: user.stepSize = dt;
2069: PetscCall(SetupContext(dm, sw, &user));
2070: PetscCall(DMSwarmVectorDefineField(sw, "velocity"));
2071: PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
2072: PetscCall(TSSetComputeExactError(ts, ComputeError));
2073: PetscCall(TSSetPostStep(ts, MigrateParticles));
2074: PetscCall(CreateSolution(ts));
2075: PetscCall(TSGetSolution(ts, &u));
2076: PetscCall(TSComputeInitialCondition(ts, u));
2077: PetscCall(CheckNonNegativeWeights(sw, &user));
2078: PetscCall(TSSolve(ts, NULL));
2080: PetscCall(SNESDestroy(&user.snes));
2081: PetscCall(TSDestroy(&ts));
2082: PetscCall(DMDestroy(&sw));
2083: PetscCall(DMDestroy(&dm));
2084: PetscCall(DestroyContext(&user));
2085: PetscCall(PetscFinalize());
2086: return 0;
2087: }
2089: /*TEST
2091: build:
2092: requires: !complex double
2094: # This tests that we can put particles in a box and compute the Coulomb force
2095: # Recommend -draw_size 500,500
2096: testset:
2097: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
2098: args: -dm_plex_dim 2 -fake_1D -dm_plex_simplex 0 -dm_plex_box_faces 20,1 \
2099: -dm_plex_box_lower 0,-1 -dm_plex_box_upper 12.5664,1 \
2100: -dm_plex_box_bd periodic,none \
2101: -dm_swarm_coordinate_density constant -dm_swarm_num_particles 100 \
2102: -sigma 1.0e-8 -timeScale 2.0e-14 \
2103: -ts_type basicsymplectic -ts_basicsymplectic_type 1 \
2104: -ts_monitor_sp_swarm -ts_monitor_sp_swarm_retain 0 -ts_monitor_sp_swarm_phase 0 \
2105: -output_step 50 -check_vel_res
2106: test:
2107: suffix: none_1d
2108: args: -em_type none -error
2109: test:
2110: suffix: coulomb_1d
2111: args: -em_type coulomb
2113: # for viewers
2114: #-ts_monitor_sp_swarm_phase -ts_monitor_sp_swarm -em_snes_monitor -ts_monitor_sp_swarm_multi_species 0 -ts_monitor_sp_swarm_retain 0
2115: testset:
2116: nsize: {{1 2}}
2117: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
2118: args: -dm_plex_dim 2 -fake_1D -dm_plex_simplex 0 -dm_plex_box_faces 36,1 \
2119: -dm_plex_box_lower 0.,-0.5 -dm_plex_box_upper 12.5664,0.5 \
2120: -dm_plex_box_bd periodic,none \
2121: -vdm_plex_dim 1 -vdm_plex_simplex 0 -vdm_plex_box_faces 10 \
2122: -vdm_plex_box_lower -3 -vdm_plex_box_upper 3 \
2123: -dm_swarm_num_species 1 -dm_swarm_num_particles 360 \
2124: -twostream -charges -1.,1. -sigma 1.0e-8 \
2125: -cosine_coefficients 0.01,0.5 -perturbed_weights -total_weight 1. \
2126: -ts_type basicsymplectic -ts_basicsymplectic_type 2 \
2127: -ts_dt 0.01 -ts_max_time 5 -ts_max_steps 10 \
2128: -em_snes_atol 1.e-15 -em_snes_error_if_not_converged -em_ksp_error_if_not_converged \
2129: -output_step 1 -check_vel_res
2130: test:
2131: suffix: two_stream_c0
2132: args: -em_type primal -petscfe_default_quadrature_order 2 -petscspace_degree 2 -em_pc_type svd
2133: test:
2134: suffix: two_stream_rt
2135: requires: superlu_dist
2136: args: -em_type mixed \
2137: -potential_petscspace_degree 0 \
2138: -potential_petscdualspace_lagrange_use_moments \
2139: -potential_petscdualspace_lagrange_moment_order 2 \
2140: -field_petscspace_degree 2 -field_petscfe_default_quadrature_order 1 \
2141: -field_petscspace_type sum \
2142: -field_petscspace_variables 2 \
2143: -field_petscspace_components 2 \
2144: -field_petscspace_sum_spaces 2 \
2145: -field_petscspace_sum_concatenate true \
2146: -field_sumcomp_0_petscspace_variables 2 \
2147: -field_sumcomp_0_petscspace_type tensor \
2148: -field_sumcomp_0_petscspace_tensor_spaces 2 \
2149: -field_sumcomp_0_petscspace_tensor_uniform false \
2150: -field_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
2151: -field_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
2152: -field_sumcomp_1_petscspace_variables 2 \
2153: -field_sumcomp_1_petscspace_type tensor \
2154: -field_sumcomp_1_petscspace_tensor_spaces 2 \
2155: -field_sumcomp_1_petscspace_tensor_uniform false \
2156: -field_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
2157: -field_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
2158: -field_petscdualspace_form_degree -1 \
2159: -field_petscdualspace_order 1 \
2160: -field_petscdualspace_lagrange_trimmed true \
2161: -em_snes_error_if_not_converged \
2162: -em_ksp_type preonly -em_ksp_error_if_not_converged \
2163: -em_pc_type fieldsplit -em_pc_fieldsplit_type schur \
2164: -em_pc_fieldsplit_schur_fact_type full -em_pc_fieldsplit_schur_precondition full \
2165: -em_fieldsplit_field_pc_type lu \
2166: -em_fieldsplit_field_pc_factor_mat_solver_type superlu_dist \
2167: -em_fieldsplit_potential_pc_type svd
2169: # For verification, we use
2170: # -dm_plex_box_faces 100,1 -vdm_plex_box_faces 8000 -dm_swarm_num_particles 800000
2171: # -ts_monitor_sp_swarm_multi_species 0 -ts_monitor_sp_swarm_retain 0 -ts_monitor_sp_swarm_phase 1 -draw_size 500,500
2172: testset:
2173: nsize: {{1 2}}
2174: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
2175: args: -dm_plex_dim 2 -fake_1D -dm_plex_simplex 0 -dm_plex_box_faces 10,1 \
2176: -dm_plex_box_lower 0.,-0.5 -dm_plex_box_upper 12.5664,0.5 \
2177: -dm_plex_box_bd periodic,none \
2178: -vdm_plex_dim 1 -vdm_plex_simplex 0 -vdm_plex_box_faces 10 \
2179: -vdm_plex_box_lower -10 -vdm_plex_box_upper 10 \
2180: -dm_swarm_num_species 1 -dm_swarm_num_particles 100 \
2181: -charges -1.,1. \
2182: -cosine_coefficients 0.01,0.5 -perturbed_weights -total_weight 1. \
2183: -ts_type basicsymplectic -ts_basicsymplectic_type 1 \
2184: -ts_dt 0.03 -ts_max_time 500 -ts_max_steps 1 \
2185: -em_snes_atol 1.e-12 -em_snes_error_if_not_converged -em_ksp_error_if_not_converged \
2186: -output_step 1 -check_vel_res
2188: test:
2189: suffix: uniform_equilibrium_1d
2190: args: -cosine_coefficients 0.0,0.5 -em_type primal -petscspace_degree 1 -em_pc_type svd
2191: test:
2192: suffix: uniform_primal_1d
2193: args: -em_type primal -petscspace_degree 1 -em_pc_type svd
2194: test:
2195: requires: superlu_dist
2196: suffix: uniform_mixed_1d
2197: args: -em_type mixed \
2198: -potential_petscspace_degree 0 \
2199: -potential_petscdualspace_lagrange_use_moments \
2200: -potential_petscdualspace_lagrange_moment_order 2 \
2201: -field_petscspace_degree 2 -field_petscfe_default_quadrature_order 1 \
2202: -field_petscspace_type sum \
2203: -field_petscspace_variables 2 \
2204: -field_petscspace_components 2 \
2205: -field_petscspace_sum_spaces 2 \
2206: -field_petscspace_sum_concatenate true \
2207: -field_sumcomp_0_petscspace_variables 2 \
2208: -field_sumcomp_0_petscspace_type tensor \
2209: -field_sumcomp_0_petscspace_tensor_spaces 2 \
2210: -field_sumcomp_0_petscspace_tensor_uniform false \
2211: -field_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
2212: -field_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
2213: -field_sumcomp_1_petscspace_variables 2 \
2214: -field_sumcomp_1_petscspace_type tensor \
2215: -field_sumcomp_1_petscspace_tensor_spaces 2 \
2216: -field_sumcomp_1_petscspace_tensor_uniform false \
2217: -field_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
2218: -field_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
2219: -field_petscdualspace_form_degree -1 \
2220: -field_petscdualspace_order 1 \
2221: -field_petscdualspace_lagrange_trimmed true \
2222: -em_snes_error_if_not_converged \
2223: -em_ksp_type preonly -em_ksp_error_if_not_converged \
2224: -em_pc_type fieldsplit -em_pc_fieldsplit_type schur \
2225: -em_pc_fieldsplit_schur_fact_type full -em_pc_fieldsplit_schur_precondition full \
2226: -em_fieldsplit_field_pc_type lu \
2227: -em_fieldsplit_field_pc_factor_mat_solver_type superlu_dist \
2228: -em_fieldsplit_potential_pc_type svd
2230: TEST*/