Actual source code: ex4.c
1: static char help[] = "Two-level system for Landau Damping using Vlasov-Poisson equations\n";
3: /*
4: Moment Equations:
6: Discretize the moment equations using finite elements and project the moments into the finite element space. Use DMSWARM_REMAP_PFAK, which guarantees that the FE approximation is weakly equivalent to the true moment. The first moment, number density, is given by
8: \int dx \phi_i n_f = \int dx \phi_i n_p
9: \int dx \phi_i n_f = \int dx \phi_i \int dv f
10: \int dx \phi_i n_f = \int dx \phi_i \int dv \sum_p w_p \delta(x - x_p) \delta(v - v_p)
11: \int dx \phi_i n_f = \int dx \phi_i \sum_p w_p \delta(x - x_p)
12: M n_F = M_p w_p
14: where
16: (M_p){ip} = \phi_i(x_p)
18: which is just a scaled version of the charge density. The second moment, momentum density, is given by
20: \int dx \phi_i p_f = m \int dx \phi_i \int dv v f
21: \int dx \phi_i p_f = m \int dx \phi_i \sum_p w_p \delta(x - x_p) v_p
22: M p_F = M_p v_p w_p
24: And finally the third moment, pressure, is given by
26: \int dx \phi_i pr_f = m \int dx \phi_i \int dv (v - u)^2 f
27: \int dx \phi_i pr_f = m \int dx \phi_i \sum_p w_p \delta(x - x_p) (v_p - u)^2
28: M pr_F = M_p (v_p - u)^2 w_p
29: = M_p (v_p - p_F(x_p) / m n_F(x_p))^2 w_p
30: = M_p (v_p - (\sum_j p_F \phi_j(x_p)) / m (\sum_k n_F \phi_k(x_p)))^2 w_p
32: Here we need all FEM basis functions \phi_i that see that particle p.
34: 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"
35: According to Lukas, good damping results come at ~16k particles per cell
37: Swarm CellDMs
38: =============
39: Name: "space"
40: Fields: DMSwarmPICField_coor, "velocity"
41: Coordinates: DMSwarmPICField_coor
43: Name: "velocity"
44: Fields: "w_q"
45: Coordinates: "velocity"
47: Name: "moments"
48: Fields: "w_q"
49: Coordinates: DMSwarmPICField_coor
51: Name: "moment fields"
52: Fields: "velocity"
53: Coordinates: DMSwarmPICField_coor
55: To visualize the maximum electric field use
57: -efield_monitor
59: To monitor velocity moments of the distribution use
61: -ptof_pc_type lu -moments_monitor
63: To monitor the particle positions in phase space use
65: -positions_monitor
67: To monitor the charge density, E field, and potential use
69: -poisson_monitor
71: To monitor the remapping field use
73: -remap_uf_view draw
75: To visualize the swarm distribution use
77: -ts_monitor_hg_swarm
79: To visualize the particles, we can use
81: -ts_monitor_sp_swarm -ts_monitor_sp_swarm_retain 0 -ts_monitor_sp_swarm_phase 1 -draw_size 500,500
82: */
83: #include <petsctao.h>
84: #include <petscts.h>
85: #include <petscdmplex.h>
86: #include <petscdmswarm.h>
87: #include <petscfe.h>
88: #include <petscds.h>
89: #include <petscbag.h>
90: #include <petscdraw.h>
91: #include <petsc/private/petscfeimpl.h>
92: #include <petsc/private/dmswarmimpl.h>
93: #include "petscdm.h"
94: #include "petscdmlabel.h"
96: PETSC_EXTERN PetscErrorCode stream(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
97: PETSC_EXTERN PetscErrorCode line(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
99: const char *EMTypes[] = {"primal", "mixed", "coulomb", "none", "EMType", "EM_", NULL};
100: typedef enum {
101: EM_PRIMAL,
102: EM_MIXED,
103: EM_COULOMB,
104: EM_NONE
105: } EMType;
107: typedef enum {
108: V0,
109: X0,
110: T0,
111: M0,
112: Q0,
113: PHI0,
114: POISSON,
115: VLASOV,
116: SIGMA,
117: NUM_CONSTANTS
118: } ConstantType;
120: typedef enum {
121: E_MONITOR_NONE,
122: E_MONITOR_FULL,
123: E_MONITOR_QUIET
124: } EMonitorType;
125: const char *const EMonitorTypes[] = {"NONE", "FULL", "QUIET", "EMonitorType", "E_MONITOR_", NULL};
127: typedef struct {
128: PetscScalar v0; /* Velocity scale, often the thermal velocity */
129: PetscScalar t0; /* Time scale */
130: PetscScalar x0; /* Space scale */
131: PetscScalar m0; /* Mass scale */
132: PetscScalar q0; /* Charge scale */
133: PetscScalar kb;
134: PetscScalar epsi0;
135: PetscScalar phi0; /* Potential scale */
136: PetscScalar poissonNumber; /* Non-Dimensional Poisson Number */
137: PetscScalar vlasovNumber; /* Non-Dimensional Vlasov Number */
138: PetscReal sigma; /* Nondimensional charge per length in x */
139: } Parameter;
141: typedef struct {
142: PetscInt s; // Starting sample (we ignore some in the beginning)
143: PetscInt e; // Ending sample
144: PetscInt per; // Period of fitting
145: const PetscReal *t; // Time for each sample
146: const PetscReal *Emax; // Emax for each sample
147: } EmaxCtx;
149: typedef struct {
150: PetscBag bag; // Problem parameters
151: PetscBool error; // Flag for printing the error
152: PetscInt remapFreq; // Number of timesteps between remapping
153: EMonitorType efield_monitor; // Flag to show electric field monitor
154: PetscBool moment_monitor; // Flag to show distribution moment monitor
155: PetscBool moment_field_monitor; // Flag to show moment field monitor
156: PetscBool positions_monitor; // Flag to show particle positins at each time step
157: PetscBool poisson_monitor; // Flag to display charge, E field, and potential at each solve
158: PetscBool initial_monitor; // Flag to monitor the initial conditions
159: PetscInt velocity_monitor; // Cell to monitor the velocity distribution for
160: PetscBool perturbed_weights; // Uniformly sample x,v space with gaussian weights
161: PetscInt ostep; // Print the energy at each ostep time steps
162: PetscInt numParticles;
163: PetscReal timeScale; /* Nondimensionalizing time scale */
164: PetscReal charges[2]; /* The charges of each species */
165: PetscReal masses[2]; /* The masses of each species */
166: PetscReal thermal_energy[2]; /* Thermal Energy (used to get other constants)*/
167: PetscReal cosine_coefficients[2]; /*(alpha, k)*/
168: PetscReal totalWeight;
169: PetscReal stepSize;
170: PetscInt steps;
171: PetscReal initVel;
172: EMType em; // Type of electrostatic model
173: SNES snes; // EM solver
174: DM dmMom; // The DM for moment fields
175: DM dmN; // The DM for number density fields
176: IS isN; // The IS mapping dmN into dmMom
177: Mat MN; // The finite element mass matrix for number density
178: DM dmP; // The DM for momentum density fields
179: IS isP; // The IS mapping dmP into dmMom
180: Mat MP; // The finite element mass matrix for momentum density
181: DM dmE; // The DM for energy density (pressure) fields
182: IS isE; // The IS mapping dmE into dmMom
183: Mat ME; // The finite element mass matrix for energy density (pressure)
184: DM dmPot; // The DM for potential
185: Mat fftPot; // Fourier Transform operator for the potential
186: Vec fftX, fftY; // FFT vectors with phases added (complex parts)
187: IS fftReal; // The indices for real parts
188: IS isPot; // The IS for potential, or NULL in primal
189: Mat M; // The finite element mass matrix for potential
190: PetscFEGeom *fegeom; // Geometric information for the DM cells
191: PetscDrawHG drawhgic_x; // Histogram of the particle weight in each X cell
192: PetscDrawHG drawhgic_v; // Histogram of the particle weight in each X cell
193: PetscDrawHG drawhgcell_v; // Histogram of the particle weight in a given cell
194: PetscBool validE; // Flag to indicate E-field in swarm is valid
195: PetscReal drawlgEmin; // The minimum lg(E) to plot
196: PetscDrawLG drawlgE; // Logarithm of maximum electric field
197: PetscDrawSP drawspE; // Electric field at particle positions
198: PetscDrawSP drawspX; // Particle positions
199: PetscViewer viewerRho; // Charge density viewer
200: PetscViewer viewerRhoHat; // Charge density Fourier Transform viewer
201: PetscViewer viewerPhi; // Potential viewer
202: PetscViewer viewerN; // Number density viewer
203: PetscViewer viewerP; // Momentum density viewer
204: PetscViewer viewerE; // Energy density (pressure) viewer
205: PetscViewer viewerNRes; // Number density residual viewer
206: PetscViewer viewerPRes; // Momentum density residual viewer
207: PetscViewer viewerERes; // Energy density (pressure) residual viewer
208: PetscDrawLG drawlgMomRes; // Residuals for the moment equations
209: DM swarm; // The particle swarm
210: PetscRandom random; // Used for particle perturbations
211: PetscBool twostream; // Flag for activating 2-stream setup
212: PetscBool checkweights; // Check weight normalization
213: PetscInt checkVRes; // Flag to check/output velocity residuals for nightly tests
214: PetscBool checkLandau; // Check the Landau damping result
215: EmaxCtx emaxCtx; // Information for fit to decay profile
216: PetscReal gamma; // The damping rate for Landau damping
217: PetscReal omega; // The perturbed oscillation frequency for Landau damping
219: PetscLogEvent RhsXEvent, RhsVEvent, ESolveEvent, ETabEvent;
220: } AppCtx;
222: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
223: {
224: PetscFunctionBeginUser;
225: PetscInt d = 2;
226: PetscInt maxSpecies = 2;
227: options->error = PETSC_FALSE;
228: options->remapFreq = 1;
229: options->efield_monitor = E_MONITOR_NONE;
230: options->moment_monitor = PETSC_FALSE;
231: options->moment_field_monitor = PETSC_FALSE;
232: options->initial_monitor = PETSC_FALSE;
233: options->perturbed_weights = PETSC_FALSE;
234: options->poisson_monitor = PETSC_FALSE;
235: options->positions_monitor = PETSC_FALSE;
236: options->velocity_monitor = -1;
237: options->ostep = 100;
238: options->timeScale = 2.0e-14;
239: options->charges[0] = -1.0;
240: options->charges[1] = 1.0;
241: options->masses[0] = 1.0;
242: options->masses[1] = 1000.0;
243: options->thermal_energy[0] = 1.0;
244: options->thermal_energy[1] = 1.0;
245: options->cosine_coefficients[0] = 0.01;
246: options->cosine_coefficients[1] = 0.5;
247: options->initVel = 1;
248: options->totalWeight = 1.0;
249: options->drawhgic_x = NULL;
250: options->drawhgic_v = NULL;
251: options->drawhgcell_v = NULL;
252: options->validE = PETSC_FALSE;
253: options->drawlgEmin = -6;
254: options->drawlgE = NULL;
255: options->drawspE = NULL;
256: options->drawspX = NULL;
257: options->viewerRho = NULL;
258: options->viewerRhoHat = NULL;
259: options->viewerPhi = NULL;
260: options->viewerN = NULL;
261: options->viewerP = NULL;
262: options->viewerE = NULL;
263: options->viewerNRes = NULL;
264: options->viewerPRes = NULL;
265: options->viewerERes = NULL;
266: options->drawlgMomRes = NULL;
267: options->em = EM_COULOMB;
268: options->snes = NULL;
269: options->dmMom = NULL;
270: options->dmN = NULL;
271: options->MN = NULL;
272: options->dmP = NULL;
273: options->MP = NULL;
274: options->dmE = NULL;
275: options->ME = NULL;
276: options->dmPot = NULL;
277: options->fftPot = NULL;
278: options->fftX = NULL;
279: options->fftY = NULL;
280: options->fftReal = NULL;
281: options->isPot = NULL;
282: options->M = NULL;
283: options->numParticles = 32768;
284: options->twostream = PETSC_FALSE;
285: options->checkweights = PETSC_FALSE;
286: options->checkVRes = 0;
287: options->checkLandau = PETSC_FALSE;
288: options->emaxCtx.s = 50;
289: options->emaxCtx.per = 100;
291: PetscOptionsBegin(comm, "", "Landau Damping and Two Stream options", "DMSWARM");
292: PetscCall(PetscOptionsBool("-error", "Flag to print the error", __FILE__, options->error, &options->error, NULL));
293: PetscCall(PetscOptionsInt("-remap_freq", "Number", __FILE__, options->remapFreq, &options->remapFreq, NULL));
294: PetscCall(PetscOptionsEnum("-efield_monitor", "Flag to record and plot log(max E) over time", __FILE__, EMonitorTypes, (PetscEnum)options->efield_monitor, (PetscEnum *)&options->efield_monitor, NULL));
295: PetscCall(PetscOptionsReal("-efield_min_monitor", "Minimum E field to plot", __FILE__, options->drawlgEmin, &options->drawlgEmin, NULL));
296: PetscCall(PetscOptionsBool("-moments_monitor", "Flag to show moments table", __FILE__, options->moment_monitor, &options->moment_monitor, NULL));
297: PetscCall(PetscOptionsBool("-moment_field_monitor", "Flag to show moment fields", __FILE__, options->moment_field_monitor, &options->moment_field_monitor, NULL));
298: PetscCall(PetscOptionsBool("-ics_monitor", "Flag to show initial condition histograms", __FILE__, options->initial_monitor, &options->initial_monitor, NULL));
299: PetscCall(PetscOptionsBool("-positions_monitor", "The flag to show particle positions", __FILE__, options->positions_monitor, &options->positions_monitor, NULL));
300: PetscCall(PetscOptionsBool("-poisson_monitor", "The flag to show charges, Efield and potential solve", __FILE__, options->poisson_monitor, &options->poisson_monitor, NULL));
301: PetscCall(PetscOptionsInt("-velocity_monitor", "Cell to show velocity histograms", __FILE__, options->velocity_monitor, &options->velocity_monitor, NULL));
302: PetscCall(PetscOptionsBool("-twostream", "Run two stream instability", __FILE__, options->twostream, &options->twostream, NULL));
303: PetscCall(PetscOptionsBool("-perturbed_weights", "Flag to run uniform sampling with perturbed weights", __FILE__, options->perturbed_weights, &options->perturbed_weights, NULL));
304: PetscCall(PetscOptionsBool("-check_weights", "Ensure all particle weights are positive", __FILE__, options->checkweights, &options->checkweights, NULL));
305: PetscCall(PetscOptionsBool("-check_landau", "Check the decay from Landau damping", __FILE__, options->checkLandau, &options->checkLandau, NULL));
306: PetscCall(PetscOptionsInt("-output_step", "Number of time steps between output", __FILE__, options->ostep, &options->ostep, NULL));
307: PetscCall(PetscOptionsReal("-timeScale", "Nondimensionalizing time scale", __FILE__, options->timeScale, &options->timeScale, NULL));
308: PetscCall(PetscOptionsInt("-check_vel_res", "Check particle velocity residuals for nightly tests", __FILE__, options->checkVRes, &options->checkVRes, NULL));
309: PetscCall(PetscOptionsReal("-initial_velocity", "Initial velocity of perturbed particle", __FILE__, options->initVel, &options->initVel, NULL));
310: PetscCall(PetscOptionsReal("-total_weight", "Total weight of all particles", __FILE__, options->totalWeight, &options->totalWeight, NULL));
311: PetscCall(PetscOptionsRealArray("-cosine_coefficients", "Amplitude and frequency of cosine equation used in initialization", __FILE__, options->cosine_coefficients, &d, NULL));
312: PetscCall(PetscOptionsRealArray("-charges", "Species charges", __FILE__, options->charges, &maxSpecies, NULL));
313: PetscCall(PetscOptionsEnum("-em_type", "Type of electrostatic solver", __FILE__, EMTypes, (PetscEnum)options->em, (PetscEnum *)&options->em, NULL));
314: PetscCall(PetscOptionsInt("-emax_start_step", "First time step to use for Emax fits", __FILE__, options->emaxCtx.s, &options->emaxCtx.s, NULL));
315: PetscCall(PetscOptionsInt("-emax_solve_step", "Number of time steps between Emax fits", __FILE__, options->emaxCtx.per, &options->emaxCtx.per, NULL));
316: PetscOptionsEnd();
318: PetscCall(PetscLogEventRegister("RhsX", TS_CLASSID, &options->RhsXEvent));
319: PetscCall(PetscLogEventRegister("RhsV", TS_CLASSID, &options->RhsVEvent));
320: PetscCall(PetscLogEventRegister("ESolve", TS_CLASSID, &options->ESolveEvent));
321: PetscCall(PetscLogEventRegister("ETab", TS_CLASSID, &options->ETabEvent));
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: static PetscErrorCode SetupContext(DM dm, DM sw, AppCtx *ctx)
326: {
327: MPI_Comm comm;
329: PetscFunctionBeginUser;
330: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
331: if (ctx->efield_monitor) {
332: PetscDraw draw;
333: PetscDrawAxis axis;
335: if (ctx->efield_monitor == E_MONITOR_FULL) {
336: PetscCall(PetscDrawCreate(comm, NULL, "Max Electric Field", 0, 0, 400, 300, &draw));
337: PetscCall(PetscDrawSetSave(draw, "ex2_Efield"));
338: PetscCall(PetscDrawSetFromOptions(draw));
339: } else {
340: PetscCall(PetscDrawOpenNull(comm, &draw));
341: }
342: PetscCall(PetscDrawLGCreate(draw, 1, &ctx->drawlgE));
343: PetscCall(PetscDrawDestroy(&draw));
344: PetscCall(PetscDrawLGGetAxis(ctx->drawlgE, &axis));
345: PetscCall(PetscDrawAxisSetLabels(axis, "Max Electric Field", "time", "E_max"));
346: PetscCall(PetscDrawLGSetLimits(ctx->drawlgE, 0., ctx->steps * ctx->stepSize, ctx->drawlgEmin, 0.));
347: }
349: if (ctx->initial_monitor) {
350: PetscDraw drawic_x, drawic_v;
351: PetscDrawAxis axis1, axis2;
352: PetscReal dmboxlower[2], dmboxupper[2];
353: PetscInt dim, cStart, cEnd;
355: PetscCall(DMGetDimension(sw, &dim));
356: PetscCall(DMGetBoundingBox(dm, dmboxlower, dmboxupper));
357: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
359: PetscCall(PetscDrawCreate(comm, NULL, "monitor_initial_conditions_x", 0, 300, 400, 300, &drawic_x));
360: PetscCall(PetscDrawSetSave(drawic_x, "ex2_ic_x"));
361: PetscCall(PetscDrawSetFromOptions(drawic_x));
362: PetscCall(PetscDrawHGCreate(drawic_x, (int)dim, &ctx->drawhgic_x));
363: PetscCall(PetscDrawHGCalcStats(ctx->drawhgic_x, PETSC_TRUE));
364: PetscCall(PetscDrawHGGetAxis(ctx->drawhgic_x, &axis1));
365: PetscCall(PetscDrawHGSetNumberBins(ctx->drawhgic_x, (int)(cEnd - cStart)));
366: PetscCall(PetscDrawAxisSetLabels(axis1, "Initial X Distribution", "X", "weight"));
367: PetscCall(PetscDrawAxisSetLimits(axis1, dmboxlower[0], dmboxupper[0], 0, 0));
368: PetscCall(PetscDrawDestroy(&drawic_x));
370: PetscCall(PetscDrawCreate(comm, NULL, "monitor_initial_conditions_v", 400, 300, 400, 300, &drawic_v));
371: PetscCall(PetscDrawSetSave(drawic_v, "ex9_ic_v"));
372: PetscCall(PetscDrawSetFromOptions(drawic_v));
373: PetscCall(PetscDrawHGCreate(drawic_v, (int)dim, &ctx->drawhgic_v));
374: PetscCall(PetscDrawHGCalcStats(ctx->drawhgic_v, PETSC_TRUE));
375: PetscCall(PetscDrawHGGetAxis(ctx->drawhgic_v, &axis2));
376: PetscCall(PetscDrawHGSetNumberBins(ctx->drawhgic_v, 21));
377: PetscCall(PetscDrawAxisSetLabels(axis2, "Initial V_x Distribution", "V", "weight"));
378: PetscCall(PetscDrawAxisSetLimits(axis2, -6, 6, 0, 0));
379: PetscCall(PetscDrawDestroy(&drawic_v));
380: }
382: if (ctx->velocity_monitor >= 0) {
383: DM vdm;
384: DMSwarmCellDM celldm;
385: PetscDraw drawcell_v;
386: PetscDrawAxis axis;
387: PetscReal dmboxlower[2], dmboxupper[2];
388: PetscInt dim;
389: char title[PETSC_MAX_PATH_LEN];
391: PetscCall(DMSwarmGetCellDMByName(sw, "velocity", &celldm));
392: PetscCall(DMSwarmCellDMGetDM(celldm, &vdm));
393: PetscCall(DMGetDimension(vdm, &dim));
394: PetscCall(DMGetBoundingBox(vdm, dmboxlower, dmboxupper));
396: PetscCall(PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Cell %" PetscInt_FMT ": Velocity Distribution", ctx->velocity_monitor));
397: PetscCall(PetscDrawCreate(comm, NULL, title, 400, 300, 400, 300, &drawcell_v));
398: PetscCall(PetscDrawSetSave(drawcell_v, "ex2_cell_v"));
399: PetscCall(PetscDrawSetFromOptions(drawcell_v));
400: PetscCall(PetscDrawHGCreate(drawcell_v, (int)dim, &ctx->drawhgcell_v));
401: PetscCall(PetscDrawHGCalcStats(ctx->drawhgcell_v, PETSC_TRUE));
402: PetscCall(PetscDrawHGGetAxis(ctx->drawhgcell_v, &axis));
403: PetscCall(PetscDrawHGSetNumberBins(ctx->drawhgcell_v, 21));
404: PetscCall(PetscDrawAxisSetLabels(axis, "V_x Distribution", "V", "weight"));
405: PetscCall(PetscDrawAxisSetLimits(axis, dmboxlower[0], dmboxupper[0], 0, 0));
406: PetscCall(PetscDrawDestroy(&drawcell_v));
407: }
409: if (ctx->positions_monitor) {
410: PetscDraw draw;
411: PetscDrawAxis axis;
413: PetscCall(PetscDrawCreate(comm, NULL, "Particle Position", 0, 0, 400, 300, &draw));
414: PetscCall(PetscDrawSetSave(draw, "ex9_pos"));
415: PetscCall(PetscDrawSetFromOptions(draw));
416: PetscCall(PetscDrawSPCreate(draw, 10, &ctx->drawspX));
417: PetscCall(PetscDrawDestroy(&draw));
418: PetscCall(PetscDrawSPSetDimension(ctx->drawspX, 1));
419: PetscCall(PetscDrawSPGetAxis(ctx->drawspX, &axis));
420: PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "x", "v"));
421: PetscCall(PetscDrawSPReset(ctx->drawspX));
422: }
423: if (ctx->poisson_monitor) {
424: Vec rho, rhohat, phi;
425: PetscDraw draw;
426: PetscDrawAxis axis;
428: PetscCall(PetscDrawCreate(comm, NULL, "Electric_Field", 0, 0, 400, 300, &draw));
429: PetscCall(PetscDrawSetFromOptions(draw));
430: PetscCall(PetscDrawSetSave(draw, "ex9_E_spatial"));
431: PetscCall(PetscDrawSPCreate(draw, 10, &ctx->drawspE));
432: PetscCall(PetscDrawDestroy(&draw));
433: PetscCall(PetscDrawSPSetDimension(ctx->drawspE, 1));
434: PetscCall(PetscDrawSPGetAxis(ctx->drawspE, &axis));
435: PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "x", "E"));
436: PetscCall(PetscDrawSPReset(ctx->drawspE));
438: PetscCall(PetscViewerDrawOpen(comm, NULL, "Charge Density", 0, 0, 400, 300, &ctx->viewerRho));
439: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->viewerRho, "rho_"));
440: PetscCall(PetscViewerDrawGetDraw(ctx->viewerRho, 0, &draw));
441: PetscCall(PetscDrawSetSave(draw, "ex9_rho_spatial"));
442: PetscCall(PetscViewerSetFromOptions(ctx->viewerRho));
443: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "rho", &rho));
444: PetscCall(PetscObjectSetName((PetscObject)rho, "charge_density"));
445: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "rho", &rho));
447: PetscInt dim, N;
449: PetscCall(DMGetDimension(ctx->dmPot, &dim));
450: if (dim == 1) {
451: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "rhohat", &rhohat));
452: PetscCall(VecGetSize(rhohat, &N));
453: PetscCall(MatCreateFFT(comm, dim, &N, MATFFTW, &ctx->fftPot));
454: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "rhohat", &rhohat));
455: PetscCall(MatCreateVecs(ctx->fftPot, &ctx->fftX, &ctx->fftY));
456: PetscCall(ISCreateStride(PETSC_COMM_SELF, N, 0, 1, &ctx->fftReal));
457: }
459: PetscCall(PetscViewerDrawOpen(comm, NULL, "rhohat: Charge Density FT", 0, 0, 400, 300, &ctx->viewerRhoHat));
460: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->viewerRhoHat, "rhohat_"));
461: PetscCall(PetscViewerDrawGetDraw(ctx->viewerRhoHat, 0, &draw));
462: PetscCall(PetscDrawSetSave(draw, "ex9_rho_ft"));
463: PetscCall(PetscViewerSetFromOptions(ctx->viewerRhoHat));
464: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "rhohat", &rhohat));
465: PetscCall(PetscObjectSetName((PetscObject)rhohat, "charge_density_ft"));
466: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "rhohat", &rhohat));
468: PetscCall(PetscViewerDrawOpen(comm, NULL, "Potential", 400, 0, 400, 300, &ctx->viewerPhi));
469: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->viewerPhi, "phi_"));
470: PetscCall(PetscViewerDrawGetDraw(ctx->viewerPhi, 0, &draw));
471: PetscCall(PetscDrawSetSave(draw, "ex9_phi_spatial"));
472: PetscCall(PetscViewerSetFromOptions(ctx->viewerPhi));
473: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "phi", &phi));
474: PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
475: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "phi", &phi));
476: }
477: if (ctx->moment_field_monitor) {
478: Vec n, p, e;
479: Vec nres, pres, eres;
480: PetscDraw draw;
482: PetscCall(PetscViewerDrawOpen(comm, NULL, "Number Density", 400, 0, 400, 300, &ctx->viewerN));
483: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->viewerN, "n_"));
484: PetscCall(PetscViewerDrawGetDraw(ctx->viewerN, 0, &draw));
485: PetscCall(PetscDrawSetSave(draw, "ex4_n_spatial"));
486: PetscCall(PetscViewerSetFromOptions(ctx->viewerN));
487: PetscCall(DMGetNamedGlobalVector(ctx->dmN, "n", &n));
488: PetscCall(PetscObjectSetName((PetscObject)n, "Number Density"));
489: PetscCall(DMRestoreNamedGlobalVector(ctx->dmN, "n", &n));
491: PetscCall(PetscViewerDrawOpen(comm, NULL, "Momentum Density", 800, 0, 400, 300, &ctx->viewerP));
492: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->viewerP, "p_"));
493: PetscCall(PetscViewerDrawGetDraw(ctx->viewerP, 0, &draw));
494: PetscCall(PetscDrawSetSave(draw, "ex4_p_spatial"));
495: PetscCall(PetscViewerSetFromOptions(ctx->viewerP));
496: PetscCall(DMGetNamedGlobalVector(ctx->dmP, "p", &p));
497: PetscCall(PetscObjectSetName((PetscObject)p, "Momentum Density"));
498: PetscCall(DMRestoreNamedGlobalVector(ctx->dmP, "p", &p));
500: PetscCall(PetscViewerDrawOpen(comm, NULL, "Emergy Density (Pressure)", 1200, 0, 400, 300, &ctx->viewerE));
501: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->viewerE, "e_"));
502: PetscCall(PetscViewerDrawGetDraw(ctx->viewerE, 0, &draw));
503: PetscCall(PetscDrawSetSave(draw, "ex4_e_spatial"));
504: PetscCall(PetscViewerSetFromOptions(ctx->viewerE));
505: PetscCall(DMGetNamedGlobalVector(ctx->dmE, "e", &e));
506: PetscCall(PetscObjectSetName((PetscObject)e, "Energy Density (Pressure)"));
507: PetscCall(DMRestoreNamedGlobalVector(ctx->dmE, "e", &e));
509: PetscDrawAxis axis;
511: PetscCall(PetscDrawCreate(comm, NULL, "Moment Residual", 0, 320, 400, 300, &draw));
512: PetscCall(PetscDrawSetSave(draw, "ex4_moment_res"));
513: PetscCall(PetscDrawSetFromOptions(draw));
514: PetscCall(PetscDrawLGCreate(draw, 3, &ctx->drawlgMomRes));
515: PetscCall(PetscDrawDestroy(&draw));
516: PetscCall(PetscDrawLGGetAxis(ctx->drawlgMomRes, &axis));
517: PetscCall(PetscDrawAxisSetLabels(axis, "Moment Residial", "time", "Residual Norm"));
518: PetscCall(PetscDrawLGSetLimits(ctx->drawlgMomRes, 0., ctx->steps * ctx->stepSize, -8, 0));
520: PetscCall(PetscViewerDrawOpen(comm, NULL, "Number Density Residual", 400, 300, 400, 300, &ctx->viewerNRes));
521: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->viewerNRes, "nres_"));
522: PetscCall(PetscViewerDrawGetDraw(ctx->viewerNRes, 0, &draw));
523: PetscCall(PetscDrawSetSave(draw, "ex4_nres_spatial"));
524: PetscCall(PetscViewerSetFromOptions(ctx->viewerNRes));
525: PetscCall(DMGetNamedGlobalVector(ctx->dmN, "nres", &nres));
526: PetscCall(PetscObjectSetName((PetscObject)nres, "Number Density Residual"));
527: PetscCall(DMRestoreNamedGlobalVector(ctx->dmN, "nres", &nres));
529: PetscCall(PetscViewerDrawOpen(comm, NULL, "Momentum Density Residual", 800, 300, 400, 300, &ctx->viewerPRes));
530: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->viewerPRes, "pres_"));
531: PetscCall(PetscViewerDrawGetDraw(ctx->viewerPRes, 0, &draw));
532: PetscCall(PetscDrawSetSave(draw, "ex4_pres_spatial"));
533: PetscCall(PetscViewerSetFromOptions(ctx->viewerPRes));
534: PetscCall(DMGetNamedGlobalVector(ctx->dmP, "pres", &pres));
535: PetscCall(PetscObjectSetName((PetscObject)pres, "Momentum Density Residual"));
536: PetscCall(DMRestoreNamedGlobalVector(ctx->dmP, "pres", &pres));
538: PetscCall(PetscViewerDrawOpen(comm, NULL, "Energy Density Residual", 1200, 300, 400, 300, &ctx->viewerERes));
539: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->viewerERes, "eres_"));
540: PetscCall(PetscViewerDrawGetDraw(ctx->viewerERes, 0, &draw));
541: PetscCall(PetscDrawSetSave(draw, "ex4_eres_spatial"));
542: PetscCall(PetscViewerSetFromOptions(ctx->viewerERes));
543: PetscCall(DMGetNamedGlobalVector(ctx->dmE, "eres", &eres));
544: PetscCall(PetscObjectSetName((PetscObject)eres, "Energy Density Residual"));
545: PetscCall(DMRestoreNamedGlobalVector(ctx->dmE, "eres", &eres));
546: }
547: PetscFunctionReturn(PETSC_SUCCESS);
548: }
550: static PetscErrorCode DestroyContext(AppCtx *ctx)
551: {
552: PetscFunctionBeginUser;
553: PetscCall(PetscDrawHGDestroy(&ctx->drawhgic_x));
554: PetscCall(PetscDrawHGDestroy(&ctx->drawhgic_v));
555: PetscCall(PetscDrawHGDestroy(&ctx->drawhgcell_v));
557: PetscCall(PetscDrawLGDestroy(&ctx->drawlgE));
558: PetscCall(PetscDrawSPDestroy(&ctx->drawspE));
559: PetscCall(PetscDrawSPDestroy(&ctx->drawspX));
560: PetscCall(PetscViewerDestroy(&ctx->viewerRho));
561: PetscCall(PetscViewerDestroy(&ctx->viewerRhoHat));
562: PetscCall(MatDestroy(&ctx->fftPot));
563: PetscCall(VecDestroy(&ctx->fftX));
564: PetscCall(VecDestroy(&ctx->fftY));
565: PetscCall(ISDestroy(&ctx->fftReal));
566: PetscCall(PetscViewerDestroy(&ctx->viewerPhi));
567: PetscCall(PetscViewerDestroy(&ctx->viewerN));
568: PetscCall(PetscViewerDestroy(&ctx->viewerP));
569: PetscCall(PetscViewerDestroy(&ctx->viewerE));
570: PetscCall(PetscViewerDestroy(&ctx->viewerNRes));
571: PetscCall(PetscViewerDestroy(&ctx->viewerPRes));
572: PetscCall(PetscViewerDestroy(&ctx->viewerERes));
573: PetscCall(PetscDrawLGDestroy(&ctx->drawlgMomRes));
575: PetscCall(PetscBagDestroy(&ctx->bag));
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: static PetscErrorCode CheckNonNegativeWeights(DM sw, AppCtx *ctx)
580: {
581: const PetscScalar *w;
582: PetscInt Np;
584: PetscFunctionBeginUser;
585: if (!ctx->checkweights) PetscFunctionReturn(PETSC_SUCCESS);
586: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&w));
587: PetscCall(DMSwarmGetLocalSize(sw, &Np));
588: 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]);
589: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: static void f0_Dirichlet(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
594: {
595: for (PetscInt d = 0; d < dim; ++d) f0[0] += 0.5 * PetscSqr(u_x[d]);
596: }
598: static PetscErrorCode computeFieldEnergy(DM dm, Vec u, PetscReal *En)
599: {
600: PetscDS ds;
601: const PetscInt field = 0;
602: PetscInt Nf;
603: void *ctx;
605: PetscFunctionBegin;
606: PetscCall(DMGetApplicationContext(dm, &ctx));
607: PetscCall(DMGetDS(dm, &ds));
608: PetscCall(PetscDSGetNumFields(ds, &Nf));
609: PetscCheck(Nf == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "We currently only support 1 field, not %" PetscInt_FMT, Nf);
610: PetscCall(PetscDSSetObjective(ds, field, &f0_Dirichlet));
611: PetscCall(DMPlexComputeIntegralFEM(dm, u, En, ctx));
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: static PetscErrorCode computeVelocityFEMMoments(DM sw, PetscReal moments[], AppCtx *ctx)
616: {
617: DMSwarmCellDM celldm;
618: DM vdm;
619: Vec u[1];
620: const char *fields[1] = {"w_q"};
622: PetscFunctionBegin;
623: PetscCall(DMSwarmSetCellDMActive(sw, "velocity"));
624: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
625: PetscCall(DMSwarmCellDMGetDM(celldm, &vdm));
626: PetscCall(DMGetGlobalVector(vdm, &u[0]));
627: PetscCall(DMSwarmProjectFields(sw, vdm, 1, fields, u, SCATTER_FORWARD));
628: PetscCall(DMPlexComputeMoments(vdm, u[0], moments));
629: PetscCall(DMRestoreGlobalVector(vdm, &u[0]));
630: PetscCall(DMSwarmSetCellDMActive(sw, "space"));
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: static void f0_grad_phi2(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
635: {
636: f0[0] = 0.;
637: for (PetscInt d = 0; d < dim; ++d) f0[0] += PetscSqr(u_x[uOff_x[0] + d * dim + d]);
638: }
640: // Our model is E_max(t) = C e^{-gamma t} |cos(omega t - phi)|
641: static PetscErrorCode ComputeEmaxResidual(Tao tao, Vec x, Vec res, void *Ctx)
642: {
643: EmaxCtx *ctx = (EmaxCtx *)Ctx;
644: const PetscScalar *a;
645: PetscScalar *F;
646: PetscReal C, gamma, omega, phi;
648: PetscFunctionBegin;
649: PetscCall(VecGetArrayRead(x, &a));
650: PetscCall(VecGetArray(res, &F));
651: C = PetscRealPart(a[0]);
652: gamma = PetscRealPart(a[1]);
653: omega = PetscRealPart(a[2]);
654: phi = PetscRealPart(a[3]);
655: PetscCall(VecRestoreArrayRead(x, &a));
656: for (PetscInt i = ctx->s; i < ctx->e; ++i) F[i - ctx->s] = PetscPowReal(10., ctx->Emax[i]) - C * PetscExpReal(-gamma * ctx->t[i]) * PetscAbsReal(PetscCosReal(omega * ctx->t[i] - phi));
657: PetscCall(VecRestoreArray(res, &F));
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: // The Jacobian of the residual J = dr(x)/dx
662: static PetscErrorCode ComputeEmaxJacobian(Tao tao, Vec x, Mat J, Mat Jpre, void *Ctx)
663: {
664: EmaxCtx *ctx = (EmaxCtx *)Ctx;
665: const PetscScalar *a;
666: PetscScalar *jac;
667: PetscReal C, gamma, omega, phi;
668: const PetscInt n = ctx->e - ctx->s;
670: PetscFunctionBegin;
671: PetscCall(VecGetArrayRead(x, &a));
672: C = PetscRealPart(a[0]);
673: gamma = PetscRealPart(a[1]);
674: omega = PetscRealPart(a[2]);
675: phi = PetscRealPart(a[3]);
676: PetscCall(VecRestoreArrayRead(x, &a));
677: PetscCall(MatDenseGetArray(J, &jac));
678: for (PetscInt i = 0; i < n; ++i) {
679: const PetscInt k = i + ctx->s;
681: jac[i * 4 + 0] = -PetscExpReal(-gamma * ctx->t[k]) * PetscAbsReal(PetscCosReal(omega * ctx->t[k] - phi));
682: jac[i * 4 + 1] = C * ctx->t[k] * PetscExpReal(-gamma * ctx->t[k]) * PetscAbsReal(PetscCosReal(omega * ctx->t[k] - phi));
683: jac[i * 4 + 2] = C * ctx->t[k] * PetscExpReal(-gamma * ctx->t[k]) * (PetscCosReal(omega * ctx->t[k] - phi) < 0. ? -1. : 1.) * PetscSinReal(omega * ctx->t[k] - phi);
684: jac[i * 4 + 3] = -C * PetscExpReal(-gamma * ctx->t[k]) * (PetscCosReal(omega * ctx->t[k] - phi) < 0. ? -1. : 1.) * PetscSinReal(omega * ctx->t[k] - phi);
685: }
686: PetscCall(MatDenseRestoreArray(J, &jac));
687: PetscFunctionReturn(PETSC_SUCCESS);
688: }
690: // Our model is log_10 E_max(t) = log_10 C - gamma t log_10 e + log_10 |cos(omega t - phi)|
691: static PetscErrorCode ComputeLogEmaxResidual(Tao tao, Vec x, Vec res, void *Ctx)
692: {
693: EmaxCtx *ctx = (EmaxCtx *)Ctx;
694: const PetscScalar *a;
695: PetscScalar *F;
696: PetscReal C, gamma, omega, phi;
698: PetscFunctionBegin;
699: PetscCall(VecGetArrayRead(x, &a));
700: PetscCall(VecGetArray(res, &F));
701: C = PetscRealPart(a[0]);
702: gamma = PetscRealPart(a[1]);
703: omega = PetscRealPart(a[2]);
704: phi = PetscRealPart(a[3]);
705: PetscCall(VecRestoreArrayRead(x, &a));
706: for (PetscInt i = ctx->s; i < ctx->e; ++i) {
707: if (C < 0) {
708: F[i - ctx->s] = 1e10;
709: continue;
710: }
711: F[i - ctx->s] = ctx->Emax[i] - (PetscLog10Real(C) - gamma * ctx->t[i] * PetscLog10Real(PETSC_E) + PetscLog10Real(PetscAbsReal(PetscCosReal(omega * ctx->t[i] - phi))));
712: }
713: PetscCall(VecRestoreArray(res, &F));
714: PetscFunctionReturn(PETSC_SUCCESS);
715: }
717: // The Jacobian of the residual J = dr(x)/dx
718: static PetscErrorCode ComputeLogEmaxJacobian(Tao tao, Vec x, Mat J, Mat Jpre, void *Ctx)
719: {
720: EmaxCtx *ctx = (EmaxCtx *)Ctx;
721: const PetscScalar *a;
722: PetscScalar *jac;
723: PetscReal C, omega, phi;
724: const PetscInt n = ctx->e - ctx->s;
726: PetscFunctionBegin;
727: PetscCall(VecGetArrayRead(x, &a));
728: C = PetscRealPart(a[0]);
729: omega = PetscRealPart(a[2]);
730: phi = PetscRealPart(a[3]);
731: PetscCall(VecRestoreArrayRead(x, &a));
732: PetscCall(MatDenseGetArray(J, &jac));
733: for (PetscInt i = 0; i < n; ++i) {
734: const PetscInt k = i + ctx->s;
736: jac[0 * n + i] = -1. / (PetscLog10Real(PETSC_E) * C);
737: jac[1 * n + i] = ctx->t[k] * PetscLog10Real(PETSC_E);
738: jac[2 * n + i] = (PetscCosReal(omega * ctx->t[k] - phi) < 0. ? -1. : 1.) * ctx->t[k] * PetscSinReal(omega * ctx->t[k] - phi) / (PetscLog10Real(PETSC_E) * PetscAbsReal(PetscCosReal(omega * ctx->t[k] - phi)));
739: jac[3 * n + i] = -(PetscCosReal(omega * ctx->t[k] - phi) < 0. ? -1. : 1.) * PetscSinReal(omega * ctx->t[k] - phi) / (PetscLog10Real(PETSC_E) * PetscAbsReal(PetscCosReal(omega * ctx->t[k] - phi)));
740: }
741: PetscCall(MatDenseRestoreArray(J, &jac));
742: PetscCall(MatViewFromOptions(J, NULL, "-emax_jac_view"));
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: static PetscErrorCode MonitorEField(TS ts, PetscInt step, PetscReal t, Vec U, void *Ctx)
747: {
748: AppCtx *ctx = (AppCtx *)Ctx;
749: DM sw;
750: PetscScalar intESq;
751: PetscReal *E, *x, *weight;
752: PetscReal Enorm = 0., lgEnorm, lgEmax, sum = 0., Emax = 0., chargesum = 0.;
753: PetscReal pmoments[4]; /* \int f, \int v f, \int v^2 f */
754: PetscInt *species, dim, Np, gNp;
755: MPI_Comm comm;
756: PetscMPIInt rank;
758: PetscFunctionBeginUser;
759: if (step < 0 || !ctx->validE) PetscFunctionReturn(PETSC_SUCCESS);
760: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
761: PetscCallMPI(MPI_Comm_rank(comm, &rank));
762: PetscCall(TSGetDM(ts, &sw));
763: PetscCall(DMGetDimension(sw, &dim));
764: PetscCall(DMSwarmGetLocalSize(sw, &Np));
765: PetscCall(DMSwarmGetSize(sw, &gNp));
766: PetscCall(DMSwarmSortGetAccess(sw));
767: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
768: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
769: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
770: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
772: for (PetscInt p = 0; p < Np; ++p) {
773: for (PetscInt d = 0; d < 1; ++d) {
774: PetscReal temp = PetscAbsReal(E[p * dim + d]);
775: if (temp > Emax) Emax = temp;
776: }
777: Enorm += PetscSqrtReal(E[p * dim] * E[p * dim]);
778: sum += E[p * dim];
779: chargesum += ctx->charges[0] * weight[p];
780: }
781: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &Emax, 1, MPIU_REAL, MPIU_MAX, comm));
782: lgEnorm = Enorm != 0 ? PetscLog10Real(Enorm) : -16.;
783: lgEmax = Emax != 0 ? PetscLog10Real(Emax) : ctx->drawlgEmin;
785: PetscDS ds;
786: Vec phi;
788: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "phi", &phi));
789: PetscCall(DMGetDS(ctx->dmPot, &ds));
790: PetscCall(PetscDSSetObjective(ds, 0, &f0_grad_phi2));
791: PetscCall(DMPlexComputeIntegralFEM(ctx->dmPot, phi, &intESq, ctx));
792: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "phi", &phi));
794: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
795: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
796: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
797: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
798: PetscCall(PetscDrawLGAddPoint(ctx->drawlgE, &t, &lgEmax));
799: if (ctx->efield_monitor == E_MONITOR_FULL) {
800: PetscDraw draw;
802: PetscCall(PetscDrawLGDraw(ctx->drawlgE));
803: PetscCall(PetscDrawLGGetDraw(ctx->drawlgE, &draw));
804: PetscCall(PetscDrawSave(draw));
806: PetscCall(DMSwarmComputeMoments(sw, "velocity", "w_q", pmoments));
807: PetscCall(PetscPrintf(comm, "E: %f\t%+e\t%e\t%f\t%20.15e\t%f\t%f\t%f\t%20.15e\t%20.15e\t%20.15e\t%" PetscInt_FMT "\t(%" PetscInt_FMT ")\n", (double)t, (double)sum, (double)Enorm, (double)lgEnorm, (double)Emax, (double)lgEmax, (double)chargesum, (double)pmoments[0], (double)pmoments[1], (double)pmoments[1 + dim], (double)PetscSqrtReal(intESq), gNp, step));
808: PetscCall(DMViewFromOptions(sw, NULL, "-sw_efield_view"));
809: }
811: // Compute decay rate and frequency
812: PetscCall(PetscDrawLGGetData(ctx->drawlgE, NULL, &ctx->emaxCtx.e, &ctx->emaxCtx.t, &ctx->emaxCtx.Emax));
813: if (!rank && !(ctx->emaxCtx.e % ctx->emaxCtx.per)) {
814: Tao tao;
815: Mat J;
816: Vec x, r;
817: PetscScalar *a;
818: PetscBool fitLog = PETSC_TRUE, debug = PETSC_FALSE;
820: PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
821: PetscCall(TaoSetOptionsPrefix(tao, "emax_"));
822: PetscCall(VecCreateSeq(PETSC_COMM_SELF, 4, &x));
823: PetscCall(TaoSetSolution(tao, x));
824: PetscCall(VecCreateSeq(PETSC_COMM_SELF, ctx->emaxCtx.e - ctx->emaxCtx.s, &r));
825: if (fitLog) PetscCall(TaoSetResidualRoutine(tao, r, ComputeLogEmaxResidual, &ctx->emaxCtx));
826: else PetscCall(TaoSetResidualRoutine(tao, r, ComputeEmaxResidual, &ctx->emaxCtx));
827: PetscCall(VecDestroy(&r));
828: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, ctx->emaxCtx.e - ctx->emaxCtx.s, 4, NULL, &J));
829: if (fitLog) PetscCall(TaoSetJacobianResidualRoutine(tao, J, J, ComputeLogEmaxJacobian, &ctx->emaxCtx));
830: else PetscCall(TaoSetJacobianResidualRoutine(tao, J, J, ComputeEmaxJacobian, &ctx->emaxCtx));
831: PetscCall(MatDestroy(&J));
832: PetscCall(TaoSetFromOptions(tao));
833: PetscCall(VecGetArray(x, &a));
834: a[0] = 0.02;
835: a[1] = 0.15;
836: a[2] = 1.4;
837: a[3] = 0.45;
838: PetscCall(VecRestoreArray(x, &a));
839: PetscCall(TaoSolve(tao));
840: if (debug) {
841: PetscCall(PetscPrintf(PETSC_COMM_SELF, "t = ["));
842: for (PetscInt i = 0; i < ctx->emaxCtx.e; ++i) {
843: if (i > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
844: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", ctx->emaxCtx.t[i]));
845: }
846: PetscCall(PetscPrintf(PETSC_COMM_SELF, "]\n"));
847: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Emax = ["));
848: for (PetscInt i = 0; i < ctx->emaxCtx.e; ++i) {
849: if (i > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
850: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", ctx->emaxCtx.Emax[i]));
851: }
852: PetscCall(PetscPrintf(PETSC_COMM_SELF, "]\n"));
853: }
854: PetscDraw draw;
855: PetscDrawAxis axis;
856: char title[PETSC_MAX_PATH_LEN];
858: PetscCall(VecGetArray(x, &a));
859: ctx->gamma = a[1];
860: ctx->omega = a[2];
861: if (ctx->efield_monitor == E_MONITOR_FULL) {
862: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Emax Fit: gamma %g omega %g C %g phi %g\n", a[1], a[2], a[0], a[3]));
863: PetscCall(PetscDrawLGGetDraw(ctx->drawlgE, &draw));
864: PetscCall(PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Max Electric Field gamma %.4g omega %.4g", a[1], a[2]));
865: PetscCall(PetscDrawSetTitle(draw, title));
866: PetscCall(PetscDrawLGGetAxis(ctx->drawlgE, &axis));
867: PetscCall(PetscDrawAxisSetLabels(axis, title, "time", "E_max"));
868: }
869: PetscCall(VecRestoreArray(x, &a));
870: PetscCall(VecDestroy(&x));
871: PetscCall(TaoDestroy(&tao));
872: }
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: static PetscErrorCode MonitorMoments(TS ts, PetscInt step, PetscReal t, Vec U, void *Ctx)
877: {
878: AppCtx *ctx = (AppCtx *)Ctx;
879: DM sw;
880: PetscReal pmoments[4], fmoments[4]; /* \int f, \int v f, \int v^2 f */
882: PetscFunctionBeginUser;
883: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS);
884: PetscCall(TSGetDM(ts, &sw));
886: PetscCall(DMSwarmComputeMoments(sw, "velocity", "w_q", pmoments));
887: PetscCall(computeVelocityFEMMoments(sw, fmoments, ctx));
889: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%f\t%f\t%f\t%f\t%f\t%f\t%f\n", (double)t, (double)pmoments[0], (double)pmoments[1], (double)pmoments[3], (double)fmoments[0], (double)fmoments[1], (double)fmoments[2]));
890: PetscFunctionReturn(PETSC_SUCCESS);
891: }
893: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
894: {
895: u[0] = 0.0;
896: return PETSC_SUCCESS;
897: }
899: /*
900: M_p w_p
901: - Make M_p with "moments"
902: - Get w_p from Swarm
903: M_p v_p w_p
904: - Get v_p from Swarm
905: - pointwise multiply v_p and w_p
906: M_p (v_p - (\sum_j p_F \phi_j(x_p)) / m (\sum_k n_F \phi_k(x_p)))^2 w_p
907: - ProjectField(sw, {n, p} U, {v_p} A, tmp_p)
908: - pointwise multiply tmp_p and w_p
910: Projection works fpr swarms
911: Fields are FE from the CellDM, and aux fields are the swarm fields
912: */
913: static PetscErrorCode ComputeMomentFields(TS ts)
914: {
915: AppCtx *ctx;
916: DM sw;
917: KSP ksp;
918: Mat M_p, D_p;
919: Vec f, v, E, tmpMom;
920: Vec m, mold, mfluxold, mres, n, nrhs, nflux, nres, p, prhs, pflux, pres, e, erhs, eflux, eres;
921: PetscReal dt, t;
922: PetscInt Nts;
924: PetscFunctionBegin;
925: PetscCall(TSGetStepNumber(ts, &Nts));
926: PetscCall(TSGetTimeStep(ts, &dt));
927: PetscCall(TSGetTime(ts, &t));
928: PetscCall(TSGetDM(ts, &sw));
929: PetscCall(DMGetApplicationContext(sw, &ctx));
930: PetscCall(DMSwarmSetCellDMActive(sw, "moment fields"));
931: PetscCall(DMSwarmMigrate(sw, PETSC_FALSE));
932: // TODO In higher dimensions, we will have to create different M_p and D_p for each field
933: PetscCall(DMCreateMassMatrix(sw, ctx->dmN, &M_p));
934: PetscCall(DMCreateGradientMatrix(sw, ctx->dmN, &D_p));
935: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
936: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &v));
937: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "E_field", &E));
938: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
940: PetscCall(MatViewFromOptions(ctx->MN, NULL, "-mn_view"));
941: PetscCall(MatViewFromOptions(ctx->MP, NULL, "-mp_view"));
942: PetscCall(MatViewFromOptions(ctx->ME, NULL, "-me_view"));
943: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
945: PetscCall(DMGetGlobalVector(ctx->dmN, &nrhs));
946: PetscCall(DMGetGlobalVector(ctx->dmN, &nflux));
947: PetscCall(PetscObjectSetName((PetscObject)nrhs, "Weak number density"));
948: PetscCall(DMGetNamedGlobalVector(ctx->dmN, "n", &n));
949: PetscCall(DMGetGlobalVector(ctx->dmP, &prhs));
950: PetscCall(DMGetGlobalVector(ctx->dmP, &pflux));
951: PetscCall(PetscObjectSetName((PetscObject)prhs, "Weak momentum density"));
952: PetscCall(DMGetNamedGlobalVector(ctx->dmP, "p", &p));
953: PetscCall(DMGetGlobalVector(ctx->dmE, &erhs));
954: PetscCall(DMGetGlobalVector(ctx->dmE, &eflux));
955: PetscCall(PetscObjectSetName((PetscObject)erhs, "Weak energy density (pressure)"));
956: PetscCall(DMGetNamedGlobalVector(ctx->dmE, "e", &e));
958: // Compute moments and fluxes
959: PetscCall(VecDuplicate(f, &tmpMom));
961: PetscCall(MatMultTranspose(M_p, f, nrhs));
963: PetscCall(VecPointwiseMult(tmpMom, f, v));
964: PetscCall(MatMultTranspose(M_p, tmpMom, prhs));
965: PetscCall(MatMultTranspose(D_p, tmpMom, nflux));
967: PetscCall(VecPointwiseMult(tmpMom, tmpMom, v));
968: PetscCall(MatMultTranspose(M_p, tmpMom, erhs));
969: PetscCall(MatMultTranspose(D_p, tmpMom, pflux));
971: PetscCall(VecPointwiseMult(tmpMom, tmpMom, v));
972: PetscCall(MatMultTranspose(D_p, tmpMom, eflux));
974: PetscCall(VecPointwiseMult(tmpMom, f, E));
975: PetscCall(MatMultTransposeAdd(M_p, tmpMom, pflux, pflux));
977: PetscCall(VecPointwiseMult(tmpMom, v, E));
978: PetscCall(VecScale(tmpMom, 2.));
979: PetscCall(MatMultTransposeAdd(M_p, tmpMom, eflux, eflux));
981: PetscCall(VecDestroy(&tmpMom));
982: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &v));
983: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
984: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "E_field", &E));
986: PetscCall(MatDestroy(&M_p));
987: PetscCall(MatDestroy(&D_p));
989: PetscCall(KSPCreate(PetscObjectComm((PetscObject)sw), &ksp));
990: PetscCall(KSPSetOptionsPrefix(ksp, "mom_proj_"));
991: PetscCall(KSPSetOperators(ksp, ctx->MN, ctx->MN));
992: PetscCall(KSPSetFromOptions(ksp));
993: PetscCall(KSPSolve(ksp, nrhs, n));
994: PetscCall(KSPSetOperators(ksp, ctx->MP, ctx->MP));
995: PetscCall(KSPSetFromOptions(ksp));
996: PetscCall(KSPSolve(ksp, prhs, p));
997: PetscCall(KSPSetOperators(ksp, ctx->ME, ctx->ME));
998: PetscCall(KSPSetFromOptions(ksp));
999: PetscCall(KSPSolve(ksp, erhs, e));
1000: PetscCall(KSPDestroy(&ksp));
1001: PetscCall(DMRestoreGlobalVector(ctx->dmN, &nrhs));
1002: PetscCall(DMRestoreGlobalVector(ctx->dmP, &prhs));
1003: PetscCall(DMRestoreGlobalVector(ctx->dmE, &erhs));
1005: // Check moment residual
1006: // TODO Fix global2local here
1007: PetscReal res[3], logres[3];
1009: PetscCall(DMGetGlobalVector(ctx->dmMom, &m));
1010: PetscCall(VecISCopy(m, ctx->isN, SCATTER_FORWARD, n));
1011: PetscCall(VecISCopy(m, ctx->isP, SCATTER_FORWARD, p));
1012: PetscCall(VecISCopy(m, ctx->isE, SCATTER_FORWARD, e));
1013: PetscCall(DMGetNamedGlobalVector(ctx->dmMom, "mold", &mold));
1014: PetscCall(DMGetNamedGlobalVector(ctx->dmMom, "mfluxold", &mfluxold));
1015: if (!Nts) goto end;
1017: // e = \Tr{\tau}
1018: // M_p w^{k+1} - M_p w^k - \Delta t D_p (w^k \vb{v}^k) = 0
1019: // M_p \vb{p}^{k+1} - M_p \vb{p}^k - \Delta t D_p \tau - e \Delta t M_p \left( n \vb{E} \right) = 0
1020: // M_p e^{k+1} - M_p e^k - \Delta t D_p \vb{Q} - 2 e \Delta t M_p \left( \vb{p} \cdot \vb{E} \right) = 0
1021: PetscCall(DMGetGlobalVector(ctx->dmMom, &mres));
1022: PetscCall(VecCopy(mfluxold, mres));
1023: PetscCall(VecAXPBYPCZ(mres, 1. / dt, -1. / dt, -1., m, mold));
1025: PetscCall(DMGetNamedGlobalVector(ctx->dmN, "nres", &nres));
1026: PetscCall(DMGetNamedGlobalVector(ctx->dmP, "pres", &pres));
1027: PetscCall(DMGetNamedGlobalVector(ctx->dmE, "eres", &eres));
1028: PetscCall(VecISCopy(mres, ctx->isN, SCATTER_REVERSE, nres));
1029: PetscCall(VecISCopy(mres, ctx->isP, SCATTER_REVERSE, pres));
1030: PetscCall(VecISCopy(mres, ctx->isE, SCATTER_REVERSE, eres));
1031: PetscCall(VecNorm(nres, NORM_2, &res[0]));
1032: PetscCall(VecNorm(pres, NORM_2, &res[1]));
1033: PetscCall(VecNorm(eres, NORM_2, &res[2]));
1034: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sw), "Mass Residual: %g\n", (double)res[0]));
1035: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sw), "Momentum Residual: %g\n", (double)res[1]));
1036: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sw), "Energy Residual: %g\n", (double)res[2]));
1037: PetscCall(DMRestoreNamedGlobalVector(ctx->dmN, "nres", &nres));
1038: PetscCall(DMRestoreNamedGlobalVector(ctx->dmP, "pres", &pres));
1039: PetscCall(DMRestoreNamedGlobalVector(ctx->dmE, "eres", &eres));
1040: PetscCall(DMRestoreGlobalVector(ctx->dmMom, &mres));
1042: for (PetscInt i = 0; i < 3; ++i) logres[i] = PetscLog10Real(res[i]);
1043: PetscCall(PetscDrawLGAddCommonPoint(ctx->drawlgMomRes, t, logres));
1044: PetscCall(PetscDrawLGDraw(ctx->drawlgMomRes));
1045: {
1046: PetscDraw draw;
1048: PetscCall(PetscDrawLGGetDraw(ctx->drawlgMomRes, &draw));
1049: PetscCall(PetscDrawSave(draw));
1050: }
1052: end:
1053: PetscCall(VecCopy(m, mold));
1054: PetscCall(DMRestoreGlobalVector(ctx->dmMom, &m));
1055: PetscCall(DMRestoreNamedGlobalVector(ctx->dmMom, "mold", &mold));
1056: PetscCall(VecISCopy(mfluxold, ctx->isN, SCATTER_FORWARD, nflux));
1057: PetscCall(VecISCopy(mfluxold, ctx->isP, SCATTER_FORWARD, pflux));
1058: PetscCall(VecISCopy(mfluxold, ctx->isE, SCATTER_FORWARD, eflux));
1059: PetscCall(DMRestoreNamedGlobalVector(ctx->dmMom, "mfluxold", &mfluxold));
1061: PetscCall(DMRestoreGlobalVector(ctx->dmN, &nflux));
1062: PetscCall(DMRestoreGlobalVector(ctx->dmP, &pflux));
1063: PetscCall(DMRestoreGlobalVector(ctx->dmE, &eflux));
1064: PetscCall(DMRestoreNamedGlobalVector(ctx->dmN, "n", &n));
1065: PetscCall(DMRestoreNamedGlobalVector(ctx->dmP, "p", &p));
1066: PetscCall(DMRestoreNamedGlobalVector(ctx->dmE, "e", &e));
1067: PetscCall(DMSwarmSetCellDMActive(sw, "space"));
1068: PetscFunctionReturn(PETSC_SUCCESS);
1069: }
1071: static PetscErrorCode MonitorMomentFields(TS ts, PetscInt step, PetscReal t, Vec U, void *Ctx)
1072: {
1073: AppCtx *ctx = (AppCtx *)Ctx;
1074: Vec n, p, e;
1075: Vec nres, pres, eres;
1077: PetscFunctionBeginUser;
1078: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS);
1079: PetscCall(ComputeMomentFields(ts));
1081: PetscCall(DMGetNamedGlobalVector(ctx->dmN, "n", &n));
1082: PetscCall(VecView(n, ctx->viewerN));
1083: PetscCall(DMRestoreNamedGlobalVector(ctx->dmN, "n", &n));
1085: PetscCall(DMGetNamedGlobalVector(ctx->dmP, "p", &p));
1086: PetscCall(VecView(p, ctx->viewerP));
1087: PetscCall(DMRestoreNamedGlobalVector(ctx->dmP, "p", &p));
1089: PetscCall(DMGetNamedGlobalVector(ctx->dmE, "e", &e));
1090: PetscCall(VecView(e, ctx->viewerE));
1091: PetscCall(DMRestoreNamedGlobalVector(ctx->dmE, "e", &e));
1093: PetscCall(DMGetNamedGlobalVector(ctx->dmN, "nres", &nres));
1094: PetscCall(VecView(nres, ctx->viewerNRes));
1095: PetscCall(DMRestoreNamedGlobalVector(ctx->dmN, "nres", &nres));
1097: PetscCall(DMGetNamedGlobalVector(ctx->dmP, "pres", &pres));
1098: PetscCall(VecView(pres, ctx->viewerPRes));
1099: PetscCall(DMRestoreNamedGlobalVector(ctx->dmP, "pres", &pres));
1101: PetscCall(DMGetNamedGlobalVector(ctx->dmE, "eres", &eres));
1102: PetscCall(VecView(eres, ctx->viewerERes));
1103: PetscCall(DMRestoreNamedGlobalVector(ctx->dmE, "eres", &eres));
1104: PetscFunctionReturn(PETSC_SUCCESS);
1105: }
1107: PetscErrorCode MonitorInitialConditions(TS ts, PetscInt step, PetscReal t, Vec U, void *Ctx)
1108: {
1109: AppCtx *ctx = (AppCtx *)Ctx;
1110: DM sw;
1111: PetscDraw drawic_x, drawic_v;
1112: PetscReal *weight, *pos, *vel;
1113: PetscInt dim, Np;
1115: PetscFunctionBegin;
1116: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* -1 indicates interpolated solution */
1117: if (step == 0) {
1118: PetscCall(TSGetDM(ts, &sw));
1119: PetscCall(DMGetDimension(sw, &dim));
1120: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1122: PetscCall(PetscDrawHGReset(ctx->drawhgic_x));
1123: PetscCall(PetscDrawHGGetDraw(ctx->drawhgic_x, &drawic_x));
1124: PetscCall(PetscDrawClear(drawic_x));
1125: PetscCall(PetscDrawFlush(drawic_x));
1127: PetscCall(PetscDrawHGReset(ctx->drawhgic_v));
1128: PetscCall(PetscDrawHGGetDraw(ctx->drawhgic_v, &drawic_v));
1129: PetscCall(PetscDrawClear(drawic_v));
1130: PetscCall(PetscDrawFlush(drawic_v));
1132: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&pos));
1133: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&vel));
1134: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
1135: for (PetscInt p = 0; p < Np; ++p) {
1136: PetscCall(PetscDrawHGAddWeightedValue(ctx->drawhgic_x, pos[p * dim], weight[p]));
1137: PetscCall(PetscDrawHGAddWeightedValue(ctx->drawhgic_v, vel[p * dim], weight[p]));
1138: }
1139: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&pos));
1140: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&vel));
1141: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
1143: PetscCall(PetscDrawHGDraw(ctx->drawhgic_x));
1144: PetscCall(PetscDrawHGSave(ctx->drawhgic_x));
1145: PetscCall(PetscDrawHGDraw(ctx->drawhgic_v));
1146: PetscCall(PetscDrawHGSave(ctx->drawhgic_v));
1147: }
1148: PetscFunctionReturn(PETSC_SUCCESS);
1149: }
1151: // Right now, make the complete velocity histogram
1152: PetscErrorCode MonitorVelocity(TS ts, PetscInt step, PetscReal t, Vec U, void *Ctx)
1153: {
1154: AppCtx *ctx = (AppCtx *)Ctx;
1155: DM sw, dm;
1156: Vec ks;
1157: PetscProbFn *cdf;
1158: PetscDraw drawcell_v;
1159: PetscScalar *ksa;
1160: PetscReal *weight, *vel;
1161: PetscInt *pidx;
1162: PetscInt dim, Npc, cStart, cEnd, cell = ctx->velocity_monitor;
1164: PetscFunctionBegin;
1165: PetscCall(TSGetDM(ts, &sw));
1166: PetscCall(DMGetDimension(sw, &dim));
1168: PetscCall(DMSwarmGetCellDM(sw, &dm));
1169: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1170: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &ks));
1171: PetscCall(PetscObjectSetName((PetscObject)ks, "KS Statistic by Cell"));
1172: PetscCall(VecSetSizes(ks, cEnd - cStart, PETSC_DETERMINE));
1173: PetscCall(VecSetFromOptions(ks));
1174: switch (dim) {
1175: case 1:
1176: //cdf = PetscCDFMaxwellBoltzmann1D;
1177: cdf = PetscCDFGaussian1D;
1178: break;
1179: case 2:
1180: cdf = PetscCDFMaxwellBoltzmann2D;
1181: break;
1182: case 3:
1183: cdf = PetscCDFMaxwellBoltzmann3D;
1184: break;
1185: default:
1186: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported", dim);
1187: }
1189: PetscCall(PetscDrawHGReset(ctx->drawhgcell_v));
1190: PetscCall(PetscDrawHGGetDraw(ctx->drawhgcell_v, &drawcell_v));
1191: PetscCall(PetscDrawClear(drawcell_v));
1192: PetscCall(PetscDrawFlush(drawcell_v));
1194: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&vel));
1195: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
1196: PetscCall(DMSwarmSortGetAccess(sw));
1197: PetscCall(VecGetArrayWrite(ks, &ksa));
1198: for (PetscInt c = cStart; c < cEnd; ++c) {
1199: Vec cellv, cellw;
1200: PetscScalar *cella, *cellaw;
1201: PetscReal totWgt = 0.;
1203: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
1204: PetscCall(VecCreate(PETSC_COMM_SELF, &cellv));
1205: PetscCall(VecSetBlockSize(cellv, dim));
1206: PetscCall(VecSetSizes(cellv, Npc * dim, Npc));
1207: PetscCall(VecSetFromOptions(cellv));
1208: PetscCall(VecCreate(PETSC_COMM_SELF, &cellw));
1209: PetscCall(VecSetSizes(cellw, Npc, Npc));
1210: PetscCall(VecSetFromOptions(cellw));
1211: PetscCall(VecGetArrayWrite(cellv, &cella));
1212: PetscCall(VecGetArrayWrite(cellw, &cellaw));
1213: for (PetscInt q = 0; q < Npc; ++q) {
1214: const PetscInt p = pidx[q];
1215: if (c == cell) PetscCall(PetscDrawHGAddWeightedValue(ctx->drawhgcell_v, vel[p * dim], weight[p]));
1216: for (PetscInt d = 0; d < dim; ++d) cella[q * dim + d] = vel[p * dim + d];
1217: cellaw[q] = weight[p];
1218: totWgt += weight[p];
1219: }
1220: PetscCall(VecRestoreArrayWrite(cellv, &cella));
1221: PetscCall(VecRestoreArrayWrite(cellw, &cellaw));
1222: PetscCall(VecScale(cellw, 1. / totWgt));
1223: PetscCall(PetscProbComputeKSStatisticWeighted(cellv, cellw, cdf, &ksa[c - cStart]));
1224: PetscCall(VecDestroy(&cellv));
1225: PetscCall(VecDestroy(&cellw));
1226: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
1227: }
1228: PetscCall(VecRestoreArrayWrite(ks, &ksa));
1229: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&vel));
1230: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
1231: PetscCall(DMSwarmSortRestoreAccess(sw));
1233: PetscReal minalpha, maxalpha;
1234: PetscInt mincell, maxcell;
1236: PetscCall(VecFilter(ks, PETSC_SMALL));
1237: PetscCall(VecMin(ks, &mincell, &minalpha));
1238: PetscCall(VecMax(ks, &maxcell, &maxalpha));
1239: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Step %" PetscInt_FMT ": Min/Max KS statistic %g/%g in cell %" PetscInt_FMT "/%" PetscInt_FMT "\n", step, minalpha, maxalpha, mincell, maxcell));
1240: PetscCall(VecViewFromOptions(ks, NULL, "-ks_view"));
1241: PetscCall(VecDestroy(&ks));
1243: PetscCall(PetscDrawHGDraw(ctx->drawhgcell_v));
1244: PetscCall(PetscDrawHGSave(ctx->drawhgcell_v));
1245: PetscFunctionReturn(PETSC_SUCCESS);
1246: }
1248: static PetscErrorCode MonitorPositions_2D(TS ts, PetscInt step, PetscReal t, Vec U, void *Ctx)
1249: {
1250: AppCtx *ctx = (AppCtx *)Ctx;
1251: DM dm, sw;
1252: PetscDrawAxis axis;
1253: char title[1024];
1254: PetscScalar *x, *v, *weight;
1255: PetscReal lower[3], upper[3], speed;
1256: const PetscInt *s;
1257: PetscInt dim, cStart, cEnd, c;
1259: PetscFunctionBeginUser;
1260: if (step > 0 && step % ctx->ostep == 0) {
1261: PetscCall(TSGetDM(ts, &sw));
1262: PetscCall(DMSwarmGetCellDM(sw, &dm));
1263: PetscCall(DMGetDimension(dm, &dim));
1264: PetscCall(DMGetBoundingBox(dm, lower, upper));
1265: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1266: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
1267: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
1268: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
1269: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&s));
1270: PetscCall(DMSwarmSortGetAccess(sw));
1271: PetscCall(PetscDrawSPReset(ctx->drawspX));
1272: PetscCall(PetscDrawSPGetAxis(ctx->drawspX, &axis));
1273: PetscCall(PetscSNPrintf(title, 1024, "Step %" PetscInt_FMT " Time: %g", step, (double)t));
1274: PetscCall(PetscDrawAxisSetLabels(axis, title, "x", "v"));
1275: PetscCall(PetscDrawSPSetLimits(ctx->drawspX, lower[0], upper[0], lower[1], upper[1]));
1276: PetscCall(PetscDrawSPSetLimits(ctx->drawspX, lower[0], upper[0], -12, 12));
1277: for (c = 0; c < cEnd - cStart; ++c) {
1278: PetscInt *pidx, Npc, q;
1279: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
1280: for (q = 0; q < Npc; ++q) {
1281: const PetscInt p = pidx[q];
1282: if (s[p] == 0) {
1283: speed = 0.;
1284: for (PetscInt d = 0; d < dim; ++d) speed += PetscSqr(v[p * dim + d]);
1285: speed = PetscSqrtReal(speed);
1286: if (dim == 1) {
1287: PetscCall(PetscDrawSPAddPointColorized(ctx->drawspX, &x[p * dim], &v[p * dim], &speed));
1288: } else {
1289: PetscCall(PetscDrawSPAddPointColorized(ctx->drawspX, &x[p * dim], &x[p * dim + 1], &speed));
1290: }
1291: } else if (s[p] == 1) {
1292: PetscCall(PetscDrawSPAddPoint(ctx->drawspX, &x[p * dim], &v[p * dim]));
1293: }
1294: }
1295: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
1296: }
1297: PetscCall(PetscDrawSPDraw(ctx->drawspX, PETSC_TRUE));
1298: PetscDraw draw;
1299: PetscCall(PetscDrawSPGetDraw(ctx->drawspX, &draw));
1300: PetscCall(PetscDrawSave(draw));
1301: PetscCall(DMSwarmSortRestoreAccess(sw));
1302: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
1303: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
1304: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
1305: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&s));
1306: }
1307: PetscFunctionReturn(PETSC_SUCCESS);
1308: }
1310: static PetscErrorCode MonitorPoisson(TS ts, PetscInt step, PetscReal t, Vec U, void *Ctx)
1311: {
1312: AppCtx *ctx = (AppCtx *)Ctx;
1313: DM dm, sw;
1315: PetscFunctionBeginUser;
1316: if (step > 0 && step % ctx->ostep == 0) {
1317: PetscCall(TSGetDM(ts, &sw));
1318: PetscCall(DMSwarmGetCellDM(sw, &dm));
1320: if (ctx->validE) {
1321: PetscScalar *x, *E, *weight;
1322: PetscReal lower[3], upper[3], xval;
1323: PetscDraw draw;
1324: PetscInt dim, cStart, cEnd;
1326: PetscCall(DMGetDimension(dm, &dim));
1327: PetscCall(DMGetBoundingBox(dm, lower, upper));
1328: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1330: PetscCall(PetscDrawSPReset(ctx->drawspE));
1331: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
1332: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
1333: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
1335: PetscCall(DMSwarmSortGetAccess(sw));
1336: for (PetscInt c = 0; c < cEnd - cStart; ++c) {
1337: PetscReal Eavg = 0.0;
1338: PetscInt *pidx, Npc;
1340: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
1341: for (PetscInt q = 0; q < Npc; ++q) {
1342: const PetscInt p = pidx[q];
1343: Eavg += E[p * dim];
1344: }
1345: Eavg /= Npc;
1346: xval = (c + 0.5) * ((upper[0] - lower[0]) / (cEnd - cStart));
1347: PetscCall(PetscDrawSPAddPoint(ctx->drawspE, &xval, &Eavg));
1348: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
1349: }
1350: PetscCall(PetscDrawSPDraw(ctx->drawspE, PETSC_TRUE));
1351: PetscCall(PetscDrawSPGetDraw(ctx->drawspE, &draw));
1352: PetscCall(PetscDrawSave(draw));
1353: PetscCall(DMSwarmSortRestoreAccess(sw));
1354: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
1355: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
1356: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
1357: }
1359: Vec rho, rhohat, phi;
1361: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "rho", &rho));
1362: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "rhohat", &rhohat));
1363: PetscCall(VecView(rho, ctx->viewerRho));
1364: PetscCall(VecISCopy(ctx->fftX, ctx->fftReal, SCATTER_FORWARD, rho));
1365: PetscCall(MatMult(ctx->fftPot, ctx->fftX, ctx->fftY));
1366: PetscCall(VecFilter(ctx->fftY, PETSC_SMALL));
1367: PetscCall(VecViewFromOptions(ctx->fftX, NULL, "-real_view"));
1368: PetscCall(VecViewFromOptions(ctx->fftY, NULL, "-fft_view"));
1369: PetscCall(VecISCopy(ctx->fftY, ctx->fftReal, SCATTER_REVERSE, rhohat));
1370: PetscCall(VecSetValue(rhohat, 0, 0., INSERT_VALUES)); // Remove large DC component
1371: PetscCall(VecView(rhohat, ctx->viewerRhoHat));
1372: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "rho", &rho));
1373: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "rhohat", &rhohat));
1375: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "phi", &phi));
1376: PetscCall(VecView(phi, ctx->viewerPhi));
1377: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "phi", &phi));
1378: }
1379: PetscFunctionReturn(PETSC_SUCCESS);
1380: }
1382: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
1383: {
1384: PetscBag bag;
1385: Parameter *p;
1387: PetscFunctionBeginUser;
1388: /* setup PETSc parameter bag */
1389: PetscCall(PetscBagGetData(ctx->bag, &p));
1390: PetscCall(PetscBagSetName(ctx->bag, "par", "Vlasov-Poisson Parameters"));
1391: bag = ctx->bag;
1392: PetscCall(PetscBagRegisterScalar(bag, &p->v0, 1.0, "v0", "Velocity scale, m/s"));
1393: PetscCall(PetscBagRegisterScalar(bag, &p->t0, 1.0, "t0", "Time scale, s"));
1394: PetscCall(PetscBagRegisterScalar(bag, &p->x0, 1.0, "x0", "Space scale, m"));
1395: PetscCall(PetscBagRegisterScalar(bag, &p->v0, 1.0, "phi0", "Potential scale, kg*m^2/A*s^3"));
1396: PetscCall(PetscBagRegisterScalar(bag, &p->q0, 1.0, "q0", "Charge Scale, A*s"));
1397: PetscCall(PetscBagRegisterScalar(bag, &p->m0, 1.0, "m0", "Mass Scale, kg"));
1398: PetscCall(PetscBagRegisterScalar(bag, &p->epsi0, 1.0, "epsi0", "Permittivity of Free Space, kg"));
1399: PetscCall(PetscBagRegisterScalar(bag, &p->kb, 1.0, "kb", "Boltzmann Constant, m^2 kg/s^2 K^1"));
1401: PetscCall(PetscBagRegisterScalar(bag, &p->sigma, 1.0, "sigma", "Charge per unit area, C/m^3"));
1402: PetscCall(PetscBagRegisterScalar(bag, &p->poissonNumber, 1.0, "poissonNumber", "Non-Dimensional Poisson Number"));
1403: PetscCall(PetscBagRegisterScalar(bag, &p->vlasovNumber, 1.0, "vlasovNumber", "Non-Dimensional Vlasov Number"));
1404: PetscCall(PetscBagSetFromOptions(bag));
1405: {
1406: PetscViewer viewer;
1407: PetscViewerFormat format;
1408: PetscBool flg;
1410: PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg));
1411: if (flg) {
1412: PetscCall(PetscViewerPushFormat(viewer, format));
1413: PetscCall(PetscBagView(bag, viewer));
1414: PetscCall(PetscViewerFlush(viewer));
1415: PetscCall(PetscViewerPopFormat(viewer));
1416: PetscCall(PetscViewerDestroy(&viewer));
1417: }
1418: }
1419: PetscFunctionReturn(PETSC_SUCCESS);
1420: }
1422: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
1423: {
1424: DMField coordField;
1425: IS cellIS;
1426: PetscQuadrature quad;
1427: PetscReal *wt, *pt;
1428: PetscInt cdim, cStart, cEnd;
1430: PetscFunctionBeginUser;
1431: PetscCall(DMCreate(comm, dm));
1432: PetscCall(DMSetType(*dm, DMPLEX));
1433: PetscCall(DMSetFromOptions(*dm));
1434: PetscCall(PetscObjectSetName((PetscObject)*dm, "space"));
1435: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
1437: // Cache the mesh geometry
1438: PetscCall(DMGetCoordinateField(*dm, &coordField));
1439: PetscCheck(coordField, comm, PETSC_ERR_USER, "DM must have a coordinate field");
1440: PetscCall(DMGetCoordinateDim(*dm, &cdim));
1441: PetscCall(DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd));
1442: PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cStart, cStart, 1, &cellIS));
1443: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
1444: PetscCall(PetscMalloc1(1, &wt));
1445: PetscCall(PetscMalloc1(2, &pt));
1446: wt[0] = 1.;
1447: pt[0] = -1.;
1448: pt[1] = -1.;
1449: PetscCall(PetscQuadratureSetData(quad, cdim, 1, 1, pt, wt));
1450: PetscCall(DMFieldCreateFEGeom(coordField, cellIS, quad, PETSC_FEGEOM_BASIC, &ctx->fegeom));
1451: PetscCall(PetscQuadratureDestroy(&quad));
1452: PetscCall(ISDestroy(&cellIS));
1453: PetscFunctionReturn(PETSC_SUCCESS);
1454: }
1456: 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[])
1457: {
1458: f0[0] = -constants[SIGMA];
1459: }
1461: 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[])
1462: {
1463: for (PetscInt d = 0; d < dim; ++d) f1[d] = u_x[d];
1464: }
1466: 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[])
1467: {
1468: for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
1469: }
1471: /*
1472: / I -grad\ / q \ = /0\
1473: \-div 0 / \phi/ \f/
1474: */
1475: 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[])
1476: {
1477: for (PetscInt d = 0; d < dim; ++d) f0[d] += u[uOff[0] + d];
1478: }
1480: 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[])
1481: {
1482: for (PetscInt d = 0; d < dim; ++d) f1[d * dim + d] = u[uOff[1]];
1483: }
1485: 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[])
1486: {
1487: f0[0] += constants[SIGMA];
1488: for (PetscInt d = 0; d < dim; ++d) f0[0] += u_x[uOff_x[0] + d * dim + d];
1489: }
1491: /* Boundary residual. Dirichlet boundary for u means u_bdy=p*n */
1492: 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[])
1493: {
1494: for (PetscInt d = 0; d < dim; ++d) g0[d * dim + d] = 1.0;
1495: }
1497: 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[])
1498: {
1499: for (PetscInt d = 0; d < dim; ++d) g2[d * dim + d] = 1.0;
1500: }
1502: 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[])
1503: {
1504: for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
1505: }
1507: static PetscErrorCode CreateFEM(DM dm, AppCtx *ctx)
1508: {
1509: PetscFE fephi, feq;
1510: PetscDS ds;
1511: PetscBool simplex;
1512: PetscInt dim;
1514: PetscFunctionBeginUser;
1515: PetscCall(DMGetDimension(dm, &dim));
1516: PetscCall(DMPlexIsSimplex(dm, &simplex));
1517: if (ctx->em == EM_MIXED) {
1518: DMLabel label;
1519: const PetscInt id = 1;
1521: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "field_", PETSC_DETERMINE, &feq));
1522: PetscCall(PetscObjectSetName((PetscObject)feq, "field"));
1523: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "potential_", PETSC_DETERMINE, &fephi));
1524: PetscCall(PetscObjectSetName((PetscObject)fephi, "potential"));
1525: PetscCall(PetscFECopyQuadrature(feq, fephi));
1526: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feq));
1527: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fephi));
1528: PetscCall(DMCreateDS(dm));
1529: PetscCall(PetscFEDestroy(&fephi));
1530: PetscCall(PetscFEDestroy(&feq));
1532: PetscCall(DMGetLabel(dm, "marker", &label));
1533: PetscCall(DMGetDS(dm, &ds));
1535: PetscCall(PetscDSSetResidual(ds, 0, f0_q, f1_q));
1536: PetscCall(PetscDSSetResidual(ds, 1, f0_phi_backgroundCharge, NULL));
1537: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_qq, NULL, NULL, NULL));
1538: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_qphi, NULL));
1539: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_phiq, NULL, NULL));
1541: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)zero, NULL, NULL, NULL));
1543: } else {
1544: MatNullSpace nullsp;
1545: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &fephi));
1546: PetscCall(PetscObjectSetName((PetscObject)fephi, "potential"));
1547: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fephi));
1548: PetscCall(DMCreateDS(dm));
1549: PetscCall(DMGetDS(dm, &ds));
1550: PetscCall(PetscDSSetResidual(ds, 0, ion_f0, laplacian_f1));
1551: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, laplacian_g3));
1552: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_TRUE, 0, NULL, &nullsp));
1553: PetscCall(PetscObjectCompose((PetscObject)fephi, "nullspace", (PetscObject)nullsp));
1554: PetscCall(MatNullSpaceDestroy(&nullsp));
1555: PetscCall(PetscFEDestroy(&fephi));
1556: }
1557: PetscFunctionReturn(PETSC_SUCCESS);
1558: }
1560: static PetscErrorCode CreatePoisson(DM dm, AppCtx *ctx)
1561: {
1562: SNES snes;
1563: Mat J;
1564: MatNullSpace nullSpace;
1566: PetscFunctionBeginUser;
1567: PetscCall(CreateFEM(dm, ctx));
1568: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
1569: PetscCall(SNESSetOptionsPrefix(snes, "em_"));
1570: PetscCall(SNESSetDM(snes, dm));
1571: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, ctx));
1572: PetscCall(SNESSetFromOptions(snes));
1574: PetscCall(DMCreateMatrix(dm, &J));
1575: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_TRUE, 0, NULL, &nullSpace));
1576: PetscCall(MatSetNullSpace(J, nullSpace));
1577: PetscCall(MatNullSpaceDestroy(&nullSpace));
1578: PetscCall(SNESSetJacobian(snes, J, J, NULL, NULL));
1579: PetscCall(MatDestroy(&J));
1580: if (ctx->em == EM_MIXED) {
1581: const PetscInt potential = 1;
1583: PetscCall(DMCreateSubDM(dm, 1, &potential, &ctx->isPot, &ctx->dmPot));
1584: } else {
1585: ctx->dmPot = dm;
1586: PetscCall(PetscObjectReference((PetscObject)ctx->dmPot));
1587: }
1588: PetscCall(DMCreateMassMatrix(ctx->dmPot, ctx->dmPot, &ctx->M));
1589: PetscCall(DMPlexCreateClosureIndex(dm, NULL));
1590: ctx->snes = snes;
1591: PetscFunctionReturn(PETSC_SUCCESS);
1592: }
1594: // Conservation of mass (m = 1.0)
1595: // n_t + 1/ m p_x = 0
1596: static void f0_mass(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[])
1597: {
1598: for (PetscInt d = 0; d < dim; ++d) f0[0] += u_t[uOff[0]] + u_x[uOff_x[1] + d * dim + d];
1599: }
1601: // Conservation of momentum (m = 1, e = 1)
1602: // p_t + (u p)_x = -pr_x + e n E
1603: // p_t + (div u) p + u . grad p = -pr_x + e n E
1604: // p_t + (div p) p / n - (p . grad n) p / n^2 + p / n . grad p = -pr_x + e n E
1605: static void f0_momentum(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[])
1606: {
1607: const PetscScalar n = u[uOff[0]];
1609: for (PetscInt d = 0; d < dim; ++d) {
1610: PetscReal divp = 0.;
1612: f0[d] += u_t[uOff[1] + d];
1613: for (PetscInt e = 0; e < dim; ++e) {
1614: f0[d] += u[uOff[1] + e] * u_x[uOff_x[1] + d * dim + e] / n; // p / n . grad p
1615: f0[d] -= (u[uOff[1] + e] * u_x[uOff_x[0] + e]) * u[uOff[1] + d] / PetscSqr(n); // -(p . grad n) p / n^2
1616: divp += u_x[uOff_x[1] + e * dim + e];
1617: }
1618: f0[d] += divp * u[uOff[1] + d] / n; // (div p) p / n
1619: f0[d] += u_x[uOff_x[2] + d]; // pr_x
1620: f0[d] -= n * a[d]; // -e n E
1621: }
1622: }
1624: // Conservation of energy
1625: // pr_t + (u pr)_x = -3 pr u_x - q_x
1626: // pr_t + (div u) pr + u . grad pr = -3 pr (div u) - q_x
1627: // pr_t + 4 (div u) pr + u . grad pr = -q_x
1628: // pr_t + 4 div p pr / n - 4 (p . grad n) pr / n^2 + p . grad pr / n = -q_x
1629: static void f0_energy(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[])
1630: {
1631: const PetscScalar n = u[uOff[0]];
1632: const PetscScalar pr = u[uOff[2]];
1633: PetscReal divp = 0.;
1635: f0[0] += u_t[uOff[2]];
1636: for (PetscInt d = 0; d < dim; ++d) {
1637: f0[0] += u[uOff[1] + d] * u_x[uOff_x[2] + d] / n; // p . grad pr / n
1638: f0[0] -= 4. * u[uOff[1] + d] * u_x[uOff_x[0] + d] * pr / PetscSqr(n); // -4 (p . grad n) pr / n^2
1639: divp += u_x[uOff_x[1] + d * dim + d];
1640: }
1641: f0[0] += 4. * divp * pr / n; // 4 div p pr / n
1642: }
1644: static PetscErrorCode SetupMomentProblem(DM dm, AppCtx *ctx)
1645: {
1646: PetscDS ds;
1648: PetscFunctionBegin;
1649: PetscCall(DMGetDS(dm, &ds));
1650: PetscCall(PetscDSSetResidual(ds, 0, f0_mass, NULL));
1651: PetscCall(PetscDSSetResidual(ds, 1, f0_momentum, NULL));
1652: PetscCall(PetscDSSetResidual(ds, 2, f0_energy, NULL));
1653: //PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_mass_uu, NULL, NULL, NULL));
1654: PetscFunctionReturn(PETSC_SUCCESS);
1655: }
1657: static PetscErrorCode CreateMomentFields(DM odm, AppCtx *ctx)
1658: {
1659: DM dm;
1660: PetscFE fe;
1661: DMPolytopeType ct;
1662: PetscInt dim, cStart;
1664: PetscFunctionBeginUser;
1665: PetscCall(DMClone(odm, &dm));
1666: PetscCall(DMGetDimension(dm, &dim));
1667: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
1668: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
1669: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, PETSC_DETERMINE, &fe));
1670: PetscCall(PetscObjectSetName((PetscObject)fe, "number density"));
1671: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
1672: PetscCall(PetscFEDestroy(&fe));
1673: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, dim, ct, NULL, PETSC_DETERMINE, &fe));
1674: PetscCall(PetscObjectSetName((PetscObject)fe, "momentum density"));
1675: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe));
1676: PetscCall(PetscFEDestroy(&fe));
1677: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, PETSC_DETERMINE, &fe));
1678: PetscCall(PetscObjectSetName((PetscObject)fe, "energy density"));
1679: PetscCall(DMSetField(dm, 2, NULL, (PetscObject)fe));
1680: PetscCall(PetscFEDestroy(&fe));
1681: PetscCall(DMCreateDS(dm));
1682: PetscCall(SetupMomentProblem(dm, ctx));
1683: ctx->dmMom = dm;
1684: PetscInt field;
1686: field = 0;
1687: PetscCall(DMCreateSubDM(ctx->dmMom, 1, &field, &ctx->isN, &ctx->dmN));
1688: PetscCall(DMCreateMassMatrix(ctx->dmN, ctx->dmN, &ctx->MN));
1689: field = 1;
1690: PetscCall(DMCreateSubDM(ctx->dmMom, 1, &field, &ctx->isP, &ctx->dmP));
1691: PetscCall(DMCreateMassMatrix(ctx->dmP, ctx->dmP, &ctx->MP));
1692: field = 2;
1693: PetscCall(DMCreateSubDM(ctx->dmMom, 1, &field, &ctx->isE, &ctx->dmE));
1694: PetscCall(DMCreateMassMatrix(ctx->dmE, ctx->dmE, &ctx->ME));
1695: PetscFunctionReturn(PETSC_SUCCESS);
1696: }
1698: PetscErrorCode PetscPDFPertubedConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
1699: {
1700: p[0] = (1 + 0.01 * PetscCosReal(0.5 * x[0])) / (2 * PETSC_PI);
1701: p[1] = (1 + 0.01 * PetscCosReal(0.5 * x[1])) / (2 * PETSC_PI);
1702: return PETSC_SUCCESS;
1703: }
1704: PetscErrorCode PetscPDFPertubedConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
1705: {
1706: p[0] = (1. + 0.01 * PetscCosReal(0.5 * x[0])) / (2 * PETSC_PI);
1707: return PETSC_SUCCESS;
1708: }
1710: PetscErrorCode PetscPDFCosine1D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
1711: {
1712: const PetscReal alpha = scale ? scale[0] : 0.0;
1713: const PetscReal k = scale ? scale[1] : 1.;
1714: p[0] = (1 + alpha * PetscCosReal(k * x[0]));
1715: return PETSC_SUCCESS;
1716: }
1718: PetscErrorCode PetscPDFCosine2D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
1719: {
1720: const PetscReal alpha = scale ? scale[0] : 0.;
1721: const PetscReal k = scale ? scale[0] : 1.;
1722: p[0] = (1 + alpha * PetscCosReal(k * (x[0] + x[1])));
1723: return PETSC_SUCCESS;
1724: }
1726: static PetscErrorCode CreateVelocityDM(DM sw, DM *vdm)
1727: {
1728: PetscFE fe;
1729: DMPolytopeType ct;
1730: PetscInt dim, cStart;
1731: const char *prefix = "v";
1733: PetscFunctionBegin;
1734: PetscCall(DMCreate(PETSC_COMM_SELF, vdm));
1735: PetscCall(DMSetType(*vdm, DMPLEX));
1736: PetscCall(DMPlexSetOptionsPrefix(*vdm, prefix));
1737: PetscCall(DMSetFromOptions(*vdm));
1738: PetscCall(PetscObjectSetName((PetscObject)*vdm, "velocity"));
1739: PetscCall(DMViewFromOptions(*vdm, NULL, "-dm_view"));
1741: PetscCall(DMGetDimension(*vdm, &dim));
1742: PetscCall(DMPlexGetHeightStratum(*vdm, 0, &cStart, NULL));
1743: PetscCall(DMPlexGetCellType(*vdm, cStart, &ct));
1744: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, prefix, PETSC_DETERMINE, &fe));
1745: PetscCall(PetscObjectSetName((PetscObject)fe, "distribution"));
1746: PetscCall(DMSetField(*vdm, 0, NULL, (PetscObject)fe));
1747: PetscCall(DMCreateDS(*vdm));
1748: PetscCall(PetscFEDestroy(&fe));
1749: PetscFunctionReturn(PETSC_SUCCESS);
1750: }
1752: /*
1753: InitializeParticles_Centroid - Initialize a regular grid of particles.
1755: Input Parameters:
1756: + sw - The `DMSWARM`
1757: - force1D - Treat the spatial domain as 1D
1759: Notes:
1760: This functions sets the species, cellid, spatial coordinate, and velocity fields for all particles.
1762: It places one particle in the centroid of each cell in the implicit tensor product of the spatial
1763: and velocity meshes.
1764: */
1765: static PetscErrorCode InitializeParticles_Centroid(DM sw)
1766: {
1767: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1768: DMSwarmCellDM celldm;
1769: DM xdm, vdm;
1770: PetscReal vmin[3], vmax[3];
1771: PetscReal *x, *v;
1772: PetscInt *species, *cellid;
1773: PetscInt dim, xcStart, xcEnd, vcStart, vcEnd, Ns, Np, Npc, debug;
1774: PetscBool flg;
1775: MPI_Comm comm;
1776: const char *cellidname;
1778: PetscFunctionBegin;
1779: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1781: PetscOptionsBegin(comm, "", "DMSwarm Options", "DMSWARM");
1782: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1783: PetscCall(PetscOptionsInt("-dm_swarm_num_species", "The number of species", "DMSwarmSetNumSpecies", Ns, &Ns, &flg));
1784: if (flg) PetscCall(DMSwarmSetNumSpecies(sw, Ns));
1785: PetscCall(PetscOptionsBoundedInt("-dm_swarm_print_coords", "Debug output level for particle coordinate computations", "InitializeParticles", 0, &swarm->printCoords, NULL, 0));
1786: PetscCall(PetscOptionsBoundedInt("-dm_swarm_print_weights", "Debug output level for particle weight computations", "InitializeWeights", 0, &swarm->printWeights, NULL, 0));
1787: PetscOptionsEnd();
1788: debug = swarm->printCoords;
1790: PetscCall(DMGetDimension(sw, &dim));
1791: PetscCall(DMSwarmGetCellDM(sw, &xdm));
1792: PetscCall(DMPlexGetHeightStratum(xdm, 0, &xcStart, &xcEnd));
1794: PetscCall(DMSwarmGetCellDMByName(sw, "velocity", &celldm));
1795: PetscCall(DMSwarmCellDMGetDM(celldm, &vdm));
1796: PetscCall(DMPlexGetHeightStratum(vdm, 0, &vcStart, &vcEnd));
1798: // One particle per centroid on the tensor product grid
1799: Npc = (vcEnd - vcStart) * Ns;
1800: Np = (xcEnd - xcStart) * Npc;
1801: PetscCall(DMSwarmSetLocalSizes(sw, Np, 0));
1802: if (debug) {
1803: PetscInt gNp, gNc, Nc = xcEnd - xcStart;
1804: PetscCallMPI(MPIU_Allreduce(&Np, &gNp, 1, MPIU_INT, MPIU_SUM, comm));
1805: PetscCall(PetscPrintf(comm, "Global Np = %" PetscInt_FMT "\n", gNp));
1806: PetscCallMPI(MPIU_Allreduce(&Nc, &gNc, 1, MPIU_INT, MPIU_SUM, comm));
1807: PetscCall(PetscPrintf(comm, "Global X-cells = %" PetscInt_FMT "\n", gNc));
1808: PetscCall(PetscPrintf(comm, "Global V-cells = %" PetscInt_FMT "\n", vcEnd - vcStart));
1809: }
1811: // Set species and cellid
1812: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1813: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellidname));
1814: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1815: PetscCall(DMSwarmGetField(sw, cellidname, NULL, NULL, (void **)&cellid));
1816: for (PetscInt c = 0, p = 0; c < xcEnd - xcStart; ++c) {
1817: for (PetscInt s = 0; s < Ns; ++s) {
1818: for (PetscInt q = 0; q < Npc / Ns; ++q, ++p) {
1819: species[p] = s;
1820: cellid[p] = c;
1821: }
1822: }
1823: }
1824: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1825: PetscCall(DMSwarmRestoreField(sw, cellidname, NULL, NULL, (void **)&cellid));
1827: // Set particle coordinates
1828: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
1829: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
1830: PetscCall(DMSwarmSortGetAccess(sw));
1831: PetscCall(DMGetBoundingBox(vdm, vmin, vmax));
1832: PetscCall(DMGetCoordinatesLocalSetUp(xdm));
1833: for (PetscInt c = 0; c < xcEnd - xcStart; ++c) {
1834: const PetscInt cell = c + xcStart;
1835: PetscInt *pidx, Npc;
1836: PetscReal centroid[3], volume;
1838: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
1839: PetscCall(DMPlexComputeCellGeometryFVM(xdm, cell, &volume, centroid, NULL));
1840: for (PetscInt s = 0; s < Ns; ++s) {
1841: for (PetscInt q = 0; q < Npc / Ns; ++q) {
1842: const PetscInt p = pidx[q * Ns + s];
1844: for (PetscInt d = 0; d < dim; ++d) {
1845: x[p * dim + d] = centroid[d];
1846: v[p * dim + d] = vmin[0] + (q + 0.5) * ((vmax[0] - vmin[0]) / (Npc / Ns));
1847: }
1848: if (debug > 1) {
1849: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Particle %4" PetscInt_FMT " ", p));
1850: PetscCall(PetscPrintf(PETSC_COMM_SELF, " x: ("));
1851: for (PetscInt d = 0; d < dim; ++d) {
1852: if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
1853: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", x[p * dim + d]));
1854: }
1855: PetscCall(PetscPrintf(PETSC_COMM_SELF, ") v:("));
1856: for (PetscInt d = 0; d < dim; ++d) {
1857: if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
1858: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", v[p * dim + d]));
1859: }
1860: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
1861: }
1862: }
1863: }
1864: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
1865: }
1866: PetscCall(DMSwarmSortRestoreAccess(sw));
1867: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
1868: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
1869: PetscFunctionReturn(PETSC_SUCCESS);
1870: }
1872: /*
1873: InitializeWeights - Compute weight for each local particle
1875: Input Parameters:
1876: + sw - The `DMSwarm`
1877: . totalWeight - The sum of all particle weights
1878: . func - The PDF for the particle spatial distribution
1879: - param - The PDF parameters
1881: Notes:
1882: The PDF for velocity is assumed to be a Gaussian
1884: The particle weights are returned in the `w_q` field of `sw`.
1885: */
1886: static PetscErrorCode InitializeWeights(DM sw, PetscReal totalWeight, PetscProbFn *func, const PetscReal param[])
1887: {
1888: DM xdm, vdm;
1889: DMSwarmCellDM celldm;
1890: PetscScalar *weight;
1891: PetscQuadrature xquad;
1892: const PetscReal *xq, *xwq;
1893: const PetscInt order = 5;
1894: PetscReal xi0[3];
1895: PetscReal xwtot = 0., pwtot = 0.;
1896: PetscInt xNq;
1897: PetscInt dim, Ns, xcStart, xcEnd, vcStart, vcEnd, debug = ((DM_Swarm *)sw->data)->printWeights;
1898: MPI_Comm comm;
1899: PetscMPIInt rank;
1901: PetscFunctionBegin;
1902: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1903: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1904: PetscCall(DMGetDimension(sw, &dim));
1905: PetscCall(DMSwarmGetCellDM(sw, &xdm));
1906: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1907: PetscCall(DMPlexGetHeightStratum(xdm, 0, &xcStart, &xcEnd));
1908: PetscCall(DMSwarmGetCellDMByName(sw, "velocity", &celldm));
1909: PetscCall(DMSwarmCellDMGetDM(celldm, &vdm));
1910: PetscCall(DMPlexGetHeightStratum(vdm, 0, &vcStart, &vcEnd));
1912: // Setup Quadrature for spatial and velocity weight calculations
1913: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, order, -1.0, 1.0, &xquad));
1914: PetscCall(PetscQuadratureGetData(xquad, NULL, NULL, &xNq, &xq, &xwq));
1915: for (PetscInt d = 0; d < dim; ++d) xi0[d] = -1.0;
1917: // Integrate the density function to get the weights of particles in each cell
1918: PetscCall(DMGetCoordinatesLocalSetUp(vdm));
1919: PetscCall(DMSwarmSortGetAccess(sw));
1920: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
1921: for (PetscInt c = xcStart; c < xcEnd; ++c) {
1922: PetscReal xv0[3], xJ[9], xinvJ[9], xdetJ, xqr[3], xden, xw = 0.;
1923: PetscInt *pidx, Npc;
1924: PetscInt xNc;
1925: const PetscScalar *xarray;
1926: PetscScalar *xcoords = NULL;
1927: PetscBool xisDG;
1929: PetscCall(DMPlexGetCellCoordinates(xdm, c, &xisDG, &xNc, &xarray, &xcoords));
1930: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
1931: PetscCheck(Npc == (vcEnd - vcStart) * Ns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of particles %" PetscInt_FMT " in cell (rank %d) != %" PetscInt_FMT " number of velocity vertices", Npc, rank, (vcEnd - vcStart) * Ns);
1932: PetscCall(DMPlexComputeCellGeometryFEM(xdm, c, NULL, xv0, xJ, xinvJ, &xdetJ));
1933: for (PetscInt q = 0; q < xNq; ++q) {
1934: // Transform quadrature points from ref space to real space
1935: CoordinatesRefToReal(dim, dim, xi0, xv0, xJ, &xq[q * dim], xqr);
1936: // Get probability density at quad point
1937: // No need to scale xqr since PDF will be periodic
1938: PetscCall((*func)(xqr, param, &xden));
1939: xw += xden * (xwq[q] * xdetJ);
1940: }
1941: xwtot += xw;
1942: if (debug) {
1943: IS globalOrdering;
1944: const PetscInt *ordering;
1946: PetscCall(DMPlexGetCellNumbering(xdm, &globalOrdering));
1947: PetscCall(ISGetIndices(globalOrdering, &ordering));
1948: PetscCall(PetscSynchronizedPrintf(comm, "c:%" PetscInt_FMT " [x_a,x_b] = %1.15f,%1.15f -> cell weight = %1.15f\n", ordering[c], (double)PetscRealPart(xcoords[0]), (double)PetscRealPart(xcoords[0 + dim]), (double)xw));
1949: PetscCall(ISRestoreIndices(globalOrdering, &ordering));
1950: }
1951: // Set weights to be Gaussian in velocity cells
1952: for (PetscInt vc = vcStart; vc < vcEnd; ++vc) {
1953: const PetscInt p = pidx[vc * Ns + 0];
1954: PetscReal vw = 0.;
1955: PetscInt vNc;
1956: const PetscScalar *varray;
1957: PetscScalar *vcoords = NULL;
1958: PetscBool visDG;
1960: PetscCall(DMPlexGetCellCoordinates(vdm, vc, &visDG, &vNc, &varray, &vcoords));
1961: // TODO: Fix 2 stream Ask Joe
1962: // Two stream function from 1/2pi v^2 e^(-v^2/2)
1963: // vw = 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.)));
1964: vw = 0.5 * (PetscErfReal(vcoords[1] / PetscSqrtReal(2.)) - PetscErfReal(vcoords[0] / PetscSqrtReal(2.)));
1966: weight[p] = totalWeight * vw * xw;
1967: pwtot += weight[p];
1968: PetscCheck(weight[p] <= 10., PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Particle %" PetscInt_FMT " weight exceeded 1: %g, %g, %g", p, xw, vw, totalWeight);
1969: PetscCall(DMPlexRestoreCellCoordinates(vdm, vc, &visDG, &vNc, &varray, &vcoords));
1970: if (debug > 1) PetscCall(PetscPrintf(comm, "particle %" PetscInt_FMT ": %g, vw: %g xw: %g\n", p, weight[p], vw, xw));
1971: }
1972: PetscCall(DMPlexRestoreCellCoordinates(xdm, c, &xisDG, &xNc, &xarray, &xcoords));
1973: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
1974: }
1975: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
1976: PetscCall(DMSwarmSortRestoreAccess(sw));
1977: PetscCall(PetscQuadratureDestroy(&xquad));
1979: if (debug) {
1980: PetscReal wtot[2] = {pwtot, xwtot}, gwtot[2];
1982: PetscCall(PetscSynchronizedFlush(comm, NULL));
1983: PetscCallMPI(MPIU_Allreduce(wtot, gwtot, 2, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
1984: PetscCall(PetscPrintf(comm, "particle weight sum = %1.10f cell weight sum = %1.10f\n", (double)gwtot[0], (double)gwtot[1]));
1985: }
1986: PetscFunctionReturn(PETSC_SUCCESS);
1987: }
1989: static PetscErrorCode InitializeParticles_PerturbedWeights(DM sw, AppCtx *ctx)
1990: {
1991: PetscReal scale[2] = {ctx->cosine_coefficients[0], ctx->cosine_coefficients[1]};
1992: PetscInt dim;
1994: PetscFunctionBegin;
1995: PetscCall(DMGetDimension(sw, &dim));
1996: PetscCall(InitializeParticles_Centroid(sw));
1997: PetscCall(InitializeWeights(sw, ctx->totalWeight, dim == 1 ? PetscPDFCosine1D : PetscPDFCosine2D, scale));
1998: PetscFunctionReturn(PETSC_SUCCESS);
1999: }
2001: static PetscErrorCode InitializeConstants(DM sw, AppCtx *ctx)
2002: {
2003: DM dm;
2004: PetscInt *species;
2005: PetscReal *weight, totalCharge = 0., totalWeight = 0., gmin[3], gmax[3], global_charge, global_weight;
2006: PetscInt Np, dim;
2008: PetscFunctionBegin;
2009: PetscCall(DMSwarmGetCellDM(sw, &dm));
2010: PetscCall(DMGetDimension(sw, &dim));
2011: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2012: PetscCall(DMGetBoundingBox(dm, gmin, gmax));
2013: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
2014: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
2015: for (PetscInt p = 0; p < Np; ++p) {
2016: totalWeight += weight[p];
2017: totalCharge += ctx->charges[species[p]] * weight[p];
2018: }
2019: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
2020: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
2021: {
2022: Parameter *param;
2023: PetscReal Area;
2025: PetscCall(PetscBagGetData(ctx->bag, ¶m));
2026: switch (dim) {
2027: case 1:
2028: Area = (gmax[0] - gmin[0]);
2029: break;
2030: case 2:
2031: Area = (gmax[0] - gmin[0]) * (gmax[1] - gmin[1]);
2032: break;
2033: case 3:
2034: Area = (gmax[0] - gmin[0]) * (gmax[1] - gmin[1]) * (gmax[2] - gmin[2]);
2035: break;
2036: default:
2037: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", dim);
2038: }
2039: PetscCallMPI(MPIU_Allreduce(&totalWeight, &global_weight, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
2040: PetscCallMPI(MPIU_Allreduce(&totalCharge, &global_charge, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
2041: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "dim = %" PetscInt_FMT "\ttotalWeight = %f, ctx->charges[species[0]] = %f\ttotalCharge = %f, Total Area = %f\n", dim, (double)global_weight, (double)ctx->charges[0], (double)global_charge, (double)Area));
2042: param->sigma = PetscAbsReal(global_charge / (Area));
2044: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sigma: %g\n", (double)param->sigma));
2045: 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,
2046: (double)param->vlasovNumber));
2047: }
2048: /* Setup Constants */
2049: {
2050: PetscDS ds;
2051: Parameter *param;
2052: PetscCall(PetscBagGetData(ctx->bag, ¶m));
2053: PetscScalar constants[NUM_CONSTANTS];
2054: constants[SIGMA] = param->sigma;
2055: constants[V0] = param->v0;
2056: constants[T0] = param->t0;
2057: constants[X0] = param->x0;
2058: constants[M0] = param->m0;
2059: constants[Q0] = param->q0;
2060: constants[PHI0] = param->phi0;
2061: constants[POISSON] = param->poissonNumber;
2062: constants[VLASOV] = param->vlasovNumber;
2063: PetscCall(DMGetDS(dm, &ds));
2064: PetscCall(PetscDSSetConstants(ds, NUM_CONSTANTS, constants));
2065: }
2066: PetscFunctionReturn(PETSC_SUCCESS);
2067: }
2069: static PetscErrorCode CreateSwarm(DM dm, AppCtx *ctx, DM *sw)
2070: {
2071: DMSwarmCellDM celldm;
2072: DM vdm;
2073: PetscReal v0[2] = {1., 0.};
2074: PetscInt dim;
2076: PetscFunctionBeginUser;
2077: PetscCall(DMGetDimension(dm, &dim));
2078: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
2079: PetscCall(DMSetType(*sw, DMSWARM));
2080: PetscCall(DMSetDimension(*sw, dim));
2081: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
2082: PetscCall(DMSetApplicationContext(*sw, ctx));
2084: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
2085: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
2086: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT));
2087: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "E_field", dim, PETSC_REAL));
2089: const char *fieldnames[2] = {DMSwarmPICField_coor, "velocity"};
2091: PetscCall(DMSwarmCellDMCreate(dm, 2, fieldnames, 1, fieldnames, &celldm));
2092: PetscCall(DMSwarmAddCellDM(*sw, celldm));
2093: PetscCall(DMSwarmCellDMDestroy(&celldm));
2095: const char *vfieldnames[2] = {"w_q"};
2097: PetscCall(CreateVelocityDM(*sw, &vdm));
2098: PetscCall(DMSwarmCellDMCreate(vdm, 1, vfieldnames, 1, &fieldnames[1], &celldm));
2099: PetscCall(DMSwarmAddCellDM(*sw, celldm));
2100: PetscCall(DMSwarmCellDMDestroy(&celldm));
2101: PetscCall(DMDestroy(&vdm));
2103: DM mdm;
2105: PetscCall(DMClone(dm, &mdm));
2106: PetscCall(PetscObjectSetName((PetscObject)mdm, "moments"));
2107: PetscCall(DMCopyDisc(dm, mdm));
2108: PetscCall(DMSwarmCellDMCreate(mdm, 1, vfieldnames, 1, fieldnames, &celldm));
2109: PetscCall(DMDestroy(&mdm));
2110: PetscCall(DMSwarmAddCellDM(*sw, celldm));
2111: PetscCall(DMSwarmCellDMDestroy(&celldm));
2113: DM mfdm;
2115: PetscCall(DMClone(dm, &mfdm));
2116: PetscCall(PetscObjectSetName((PetscObject)mfdm, "moment fields"));
2117: PetscCall(DMCopyDisc(dm, mfdm));
2118: PetscCall(DMSwarmCellDMCreate(mfdm, 1, &fieldnames[1], 1, fieldnames, &celldm));
2119: PetscCall(DMDestroy(&mfdm));
2120: PetscCall(DMSwarmAddCellDM(*sw, celldm));
2121: PetscCall(DMSwarmCellDMDestroy(&celldm));
2123: PetscCall(DMSetFromOptions(*sw));
2124: PetscCall(DMSetUp(*sw));
2126: PetscCall(DMSwarmSetCellDMActive(*sw, "space"));
2127: ctx->swarm = *sw;
2128: // TODO: This is redundant init since it is done in InitializeSolveAndSwarm, however DMSetUp() requires the local size be set
2129: if (ctx->perturbed_weights) {
2130: PetscCall(InitializeParticles_PerturbedWeights(*sw, ctx));
2131: } else {
2132: PetscCall(DMSwarmComputeLocalSizeFromOptions(*sw));
2133: PetscCall(DMSwarmInitializeCoordinates(*sw));
2134: PetscCall(DMSwarmInitializeVelocitiesFromOptions(*sw, v0));
2135: }
2136: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
2137: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
2138: PetscFunctionReturn(PETSC_SUCCESS);
2139: }
2141: static PetscErrorCode ComputeFieldAtParticles_Coulomb(SNES snes, DM sw, PetscReal E[])
2142: {
2143: AppCtx *ctx;
2144: PetscReal *coords;
2145: PetscInt *species, dim, Np, Ns;
2146: PetscMPIInt size;
2148: PetscFunctionBegin;
2149: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)snes), &size));
2150: PetscCheck(size == 1, PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Coulomb code only works in serial");
2151: PetscCall(DMGetDimension(sw, &dim));
2152: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2153: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
2154: PetscCall(DMGetApplicationContext(sw, &ctx));
2156: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
2157: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
2158: for (PetscInt p = 0; p < Np; ++p) {
2159: PetscReal *pcoord = &coords[p * dim];
2160: PetscReal pE[3] = {0., 0., 0.};
2162: /* Calculate field at particle p due to particle q */
2163: for (PetscInt q = 0; q < Np; ++q) {
2164: PetscReal *qcoord = &coords[q * dim];
2165: PetscReal rpq[3], r, r3, q_q;
2167: if (p == q) continue;
2168: q_q = ctx->charges[species[q]] * 1.;
2169: for (PetscInt d = 0; d < dim; ++d) rpq[d] = pcoord[d] - qcoord[d];
2170: r = DMPlex_NormD_Internal(dim, rpq);
2171: if (r < PETSC_SQRT_MACHINE_EPSILON) continue;
2172: r3 = PetscPowRealInt(r, 3);
2173: for (PetscInt d = 0; d < dim; ++d) pE[d] += q_q * rpq[d] / r3;
2174: }
2175: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] = pE[d];
2176: }
2177: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
2178: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
2179: PetscFunctionReturn(PETSC_SUCCESS);
2180: }
2182: static PetscErrorCode ComputeFieldAtParticles_Primal(SNES snes, DM sw, Mat M_p, PetscReal E[])
2183: {
2184: DM dm;
2185: AppCtx *ctx;
2186: PetscDS ds;
2187: PetscFE fe;
2188: KSP ksp;
2189: Vec rhoRhs; // Weak charge density, \int phi_i rho
2190: Vec rho; // Charge density, M^{-1} rhoRhs
2191: Vec phi, locPhi; // Potential
2192: Vec f; // Particle weights
2193: PetscReal *coords;
2194: PetscInt dim, cStart, cEnd, Np;
2196: PetscFunctionBegin;
2197: PetscCall(DMGetApplicationContext(sw, &ctx));
2198: PetscCall(PetscLogEventBegin(ctx->ESolveEvent, snes, sw, 0, 0));
2199: PetscCall(DMGetDimension(sw, &dim));
2200: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2202: PetscCall(SNESGetDM(snes, &dm));
2203: PetscCall(DMGetGlobalVector(dm, &rhoRhs));
2204: PetscCall(PetscObjectSetName((PetscObject)rhoRhs, "Weak charge density"));
2205: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "rho", &rho));
2206: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
2207: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
2209: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
2210: PetscCall(MatViewFromOptions(ctx->M, NULL, "-m_view"));
2211: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
2213: PetscCall(MatMultTranspose(M_p, f, rhoRhs));
2214: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
2216: PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &ksp));
2217: PetscCall(KSPSetOptionsPrefix(ksp, "em_proj_"));
2218: PetscCall(KSPSetOperators(ksp, ctx->M, ctx->M));
2219: PetscCall(KSPSetFromOptions(ksp));
2220: PetscCall(KSPSolve(ksp, rhoRhs, rho));
2222: PetscCall(VecScale(rhoRhs, -1.0));
2224: PetscCall(VecViewFromOptions(rhoRhs, NULL, "-rho_view"));
2225: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "rho", &rho));
2226: PetscCall(KSPDestroy(&ksp));
2228: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "phi", &phi));
2229: PetscCall(VecSet(phi, 0.0));
2230: PetscCall(SNESSolve(snes, rhoRhs, phi));
2231: PetscCall(DMRestoreGlobalVector(dm, &rhoRhs));
2232: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
2234: PetscCall(DMGetLocalVector(dm, &locPhi));
2235: PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
2236: PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
2237: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "phi", &phi));
2238: PetscCall(PetscLogEventEnd(ctx->ESolveEvent, snes, sw, 0, 0));
2240: PetscCall(DMGetDS(dm, &ds));
2241: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
2242: PetscCall(DMSwarmSortGetAccess(sw));
2243: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2244: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
2246: PetscCall(PetscLogEventBegin(ctx->ETabEvent, snes, sw, 0, 0));
2247: PetscTabulation tab;
2248: PetscReal *pcoord, *refcoord;
2249: PetscFEGeom *chunkgeom = NULL;
2250: PetscInt maxNcp = 0;
2252: for (PetscInt c = cStart; c < cEnd; ++c) {
2253: PetscInt Ncp;
2255: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, c, &Ncp));
2256: maxNcp = PetscMax(maxNcp, Ncp);
2257: }
2258: PetscCall(DMGetWorkArray(dm, maxNcp * dim, MPIU_REAL, &refcoord));
2259: PetscCall(DMGetWorkArray(dm, maxNcp * dim, MPIU_REAL, &pcoord));
2260: // This can raise an FP_INEXACT in the dgemm inside
2261: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2262: PetscCall(PetscFECreateTabulation(fe, 1, maxNcp, refcoord, 1, &tab));
2263: PetscCall(PetscFPTrapPop());
2264: for (PetscInt c = cStart; c < cEnd; ++c) {
2265: PetscScalar *clPhi = NULL;
2266: PetscInt *points;
2267: PetscInt Ncp;
2269: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
2270: for (PetscInt cp = 0; cp < Ncp; ++cp)
2271: for (PetscInt d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
2272: {
2273: PetscCall(PetscFEGeomGetChunk(ctx->fegeom, c - cStart, c - cStart + 1, &chunkgeom));
2274: for (PetscInt i = 0; i < Ncp; ++i) {
2275: const PetscReal x0[3] = {-1., -1., -1.};
2276: CoordinatesRealToRef(dim, dim, x0, chunkgeom->v, chunkgeom->invJ, &pcoord[dim * i], &refcoord[dim * i]);
2277: }
2278: }
2279: PetscCall(PetscFEComputeTabulation(fe, Ncp, refcoord, 1, tab));
2280: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
2281: for (PetscInt cp = 0; cp < Ncp; ++cp) {
2282: const PetscReal *basisDer = tab->T[1];
2283: const PetscInt p = points[cp];
2285: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] = 0.;
2286: PetscCall(PetscFEFreeInterpolateGradient_Static(fe, basisDer, clPhi, dim, chunkgeom->invJ, NULL, cp, &E[p * dim]));
2287: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] *= -1.0;
2288: }
2289: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
2290: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
2291: }
2292: PetscCall(DMRestoreWorkArray(dm, maxNcp * dim, MPIU_REAL, &pcoord));
2293: PetscCall(DMRestoreWorkArray(dm, maxNcp * dim, MPIU_REAL, &refcoord));
2294: PetscCall(PetscTabulationDestroy(&tab));
2295: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
2296: PetscCall(DMSwarmSortRestoreAccess(sw));
2297: PetscCall(DMRestoreLocalVector(dm, &locPhi));
2298: PetscCall(PetscFEGeomRestoreChunk(ctx->fegeom, 0, 1, &chunkgeom));
2299: PetscCall(PetscLogEventEnd(ctx->ETabEvent, snes, sw, 0, 0));
2300: PetscFunctionReturn(PETSC_SUCCESS);
2301: }
2303: static PetscErrorCode ComputeFieldAtParticles_Mixed(SNES snes, DM sw, Mat M_p, PetscReal E[])
2304: {
2305: DM dm;
2306: AppCtx *ctx;
2307: PetscDS ds;
2308: PetscFE fe;
2309: KSP ksp;
2310: Vec rhoRhs, rhoRhsFull; // Weak charge density, \int phi_i rho, and embedding in mixed problem
2311: Vec rho; // Charge density, M^{-1} rhoRhs
2312: Vec phi, locPhi, phiFull; // Potential and embedding in mixed problem
2313: Vec f; // Particle weights
2314: PetscReal *coords;
2315: PetscInt dim, cStart, cEnd, Np;
2317: PetscFunctionBegin;
2318: PetscCall(DMGetApplicationContext(sw, &ctx));
2319: PetscCall(PetscLogEventBegin(ctx->ESolveEvent, snes, sw, 0, 0));
2320: PetscCall(DMGetDimension(sw, &dim));
2321: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2323: PetscCall(SNESGetDM(snes, &dm));
2324: PetscCall(DMGetGlobalVector(ctx->dmPot, &rhoRhs));
2325: PetscCall(PetscObjectSetName((PetscObject)rhoRhs, "Weak charge density"));
2326: PetscCall(DMGetGlobalVector(dm, &rhoRhsFull));
2327: PetscCall(PetscObjectSetName((PetscObject)rhoRhsFull, "Weak charge density"));
2328: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "rho", &rho));
2329: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
2330: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
2332: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
2333: PetscCall(MatViewFromOptions(ctx->M, NULL, "-m_view"));
2334: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
2336: PetscCall(MatMultTranspose(M_p, f, rhoRhs));
2337: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
2339: PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &ksp));
2340: PetscCall(KSPSetOptionsPrefix(ksp, "em_proj"));
2341: PetscCall(KSPSetOperators(ksp, ctx->M, ctx->M));
2342: PetscCall(KSPSetFromOptions(ksp));
2343: PetscCall(KSPSolve(ksp, rhoRhs, rho));
2345: PetscCall(VecISCopy(rhoRhsFull, ctx->isPot, SCATTER_FORWARD, rhoRhs));
2346: //PetscCall(VecScale(rhoRhsFull, -1.0));
2348: PetscCall(VecViewFromOptions(rhoRhs, NULL, "-rho_view"));
2349: PetscCall(VecViewFromOptions(rhoRhsFull, NULL, "-rho_full_view"));
2350: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "rho", &rho));
2351: PetscCall(DMRestoreGlobalVector(ctx->dmPot, &rhoRhs));
2352: PetscCall(KSPDestroy(&ksp));
2354: PetscCall(DMGetGlobalVector(dm, &phiFull));
2355: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "phi", &phi));
2356: PetscCall(VecSet(phiFull, 0.0));
2357: PetscCall(SNESSolve(snes, rhoRhsFull, phiFull));
2358: PetscCall(DMRestoreGlobalVector(dm, &rhoRhsFull));
2359: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
2361: PetscCall(VecISCopy(phiFull, ctx->isPot, SCATTER_REVERSE, phi));
2362: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "phi", &phi));
2364: PetscCall(DMGetLocalVector(dm, &locPhi));
2365: PetscCall(DMGlobalToLocalBegin(dm, phiFull, INSERT_VALUES, locPhi));
2366: PetscCall(DMGlobalToLocalEnd(dm, phiFull, INSERT_VALUES, locPhi));
2367: PetscCall(DMRestoreGlobalVector(dm, &phiFull));
2368: PetscCall(PetscLogEventEnd(ctx->ESolveEvent, snes, sw, 0, 0));
2370: PetscCall(DMGetDS(dm, &ds));
2371: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
2372: PetscCall(DMSwarmSortGetAccess(sw));
2373: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2374: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
2376: PetscCall(PetscLogEventBegin(ctx->ETabEvent, snes, sw, 0, 0));
2377: PetscTabulation tab;
2378: PetscReal *pcoord, *refcoord;
2379: PetscFEGeom *chunkgeom = NULL;
2380: PetscInt maxNcp = 0;
2382: for (PetscInt c = cStart; c < cEnd; ++c) {
2383: PetscInt Ncp;
2385: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, c, &Ncp));
2386: maxNcp = PetscMax(maxNcp, Ncp);
2387: }
2388: PetscCall(DMGetWorkArray(dm, maxNcp * dim, MPIU_REAL, &refcoord));
2389: PetscCall(DMGetWorkArray(dm, maxNcp * dim, MPIU_REAL, &pcoord));
2390: PetscCall(PetscFECreateTabulation(fe, 1, maxNcp, refcoord, 1, &tab));
2391: for (PetscInt c = cStart; c < cEnd; ++c) {
2392: PetscScalar *clPhi = NULL;
2393: PetscInt *points;
2394: PetscInt Ncp;
2396: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
2397: for (PetscInt cp = 0; cp < Ncp; ++cp)
2398: for (PetscInt d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
2399: {
2400: PetscCall(PetscFEGeomGetChunk(ctx->fegeom, c - cStart, c - cStart + 1, &chunkgeom));
2401: for (PetscInt i = 0; i < Ncp; ++i) {
2402: const PetscReal x0[3] = {-1., -1., -1.};
2403: CoordinatesRealToRef(dim, dim, x0, chunkgeom->v, chunkgeom->invJ, &pcoord[dim * i], &refcoord[dim * i]);
2404: }
2405: }
2406: PetscCall(PetscFEComputeTabulation(fe, Ncp, refcoord, 1, tab));
2407: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
2408: for (PetscInt cp = 0; cp < Ncp; ++cp) {
2409: const PetscInt p = points[cp];
2411: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] = 0.;
2412: PetscCall(PetscFEInterpolateAtPoints_Static(fe, tab, clPhi, chunkgeom, cp, &E[p * dim]));
2413: PetscCall(PetscFEPushforward(fe, chunkgeom, 1, &E[p * dim]));
2414: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] *= -1.0;
2415: }
2416: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
2417: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
2418: }
2419: PetscCall(DMRestoreWorkArray(dm, maxNcp * dim, MPIU_REAL, &pcoord));
2420: PetscCall(DMRestoreWorkArray(dm, maxNcp * dim, MPIU_REAL, &refcoord));
2421: PetscCall(PetscTabulationDestroy(&tab));
2422: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
2423: PetscCall(DMSwarmSortRestoreAccess(sw));
2424: PetscCall(DMRestoreLocalVector(dm, &locPhi));
2425: PetscCall(PetscFEGeomRestoreChunk(ctx->fegeom, 0, 1, &chunkgeom));
2426: PetscCall(PetscLogEventEnd(ctx->ETabEvent, snes, sw, 0, 0));
2427: PetscFunctionReturn(PETSC_SUCCESS);
2428: }
2430: static PetscErrorCode ComputeFieldAtParticles(SNES snes, DM sw)
2431: {
2432: AppCtx *ctx;
2433: Mat M_p;
2434: PetscReal *E;
2435: PetscInt dim, Np;
2437: PetscFunctionBegin;
2440: PetscCall(DMGetDimension(sw, &dim));
2441: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2442: PetscCall(DMGetApplicationContext(sw, &ctx));
2444: PetscCall(DMSwarmSetCellDMActive(sw, "moments"));
2445: // TODO: Could share sort context with space cellDM
2446: PetscCall(DMSwarmMigrate(sw, PETSC_FALSE));
2447: PetscCall(DMCreateMassMatrix(sw, ctx->dmPot, &M_p));
2448: PetscCall(DMSwarmSetCellDMActive(sw, "space"));
2450: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
2451: PetscCall(PetscArrayzero(E, Np * dim));
2452: ctx->validE = PETSC_TRUE;
2454: switch (ctx->em) {
2455: case EM_COULOMB:
2456: PetscCall(ComputeFieldAtParticles_Coulomb(snes, sw, E));
2457: break;
2458: case EM_PRIMAL:
2459: PetscCall(ComputeFieldAtParticles_Primal(snes, sw, M_p, E));
2460: break;
2461: case EM_MIXED:
2462: PetscCall(ComputeFieldAtParticles_Mixed(snes, sw, M_p, E));
2463: break;
2464: case EM_NONE:
2465: break;
2466: default:
2467: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No solver for electrostatic model %s", EMTypes[ctx->em]);
2468: }
2469: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
2470: PetscCall(MatDestroy(&M_p));
2471: PetscFunctionReturn(PETSC_SUCCESS);
2472: }
2474: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, PetscCtx ctx)
2475: {
2476: DM sw;
2477: SNES snes = ((AppCtx *)ctx)->snes;
2478: const PetscScalar *u;
2479: PetscScalar *g;
2480: PetscReal *E, m_p = 1., q_p = -1.;
2481: PetscInt dim, d, Np, p;
2483: PetscFunctionBeginUser;
2484: PetscCall(TSGetDM(ts, &sw));
2485: PetscCall(ComputeFieldAtParticles(snes, sw));
2487: PetscCall(DMGetDimension(sw, &dim));
2488: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2489: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
2490: PetscCall(VecGetArrayRead(U, &u));
2491: PetscCall(VecGetArray(G, &g));
2492: Np /= 2 * dim;
2493: for (p = 0; p < Np; ++p) {
2494: for (d = 0; d < dim; ++d) {
2495: g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d];
2496: g[(p * 2 + 1) * dim + d] = q_p * E[p * dim + d] / m_p;
2497: }
2498: }
2499: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
2500: PetscCall(VecRestoreArrayRead(U, &u));
2501: PetscCall(VecRestoreArray(G, &g));
2502: PetscFunctionReturn(PETSC_SUCCESS);
2503: }
2505: /* J_{ij} = dF_i/dx_j
2506: J_p = ( 0 1)
2507: (-w^2 0)
2508: TODO Now there is another term with w^2 from the electric field. I think we will need to invert the operator.
2509: Perhaps we can approximate the Jacobian using only the cellwise P-P gradient from Coulomb
2510: */
2511: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, PetscCtx ctx)
2512: {
2513: DM sw;
2514: const PetscReal *coords, *vel;
2515: PetscInt dim, d, Np, p, rStart;
2517: PetscFunctionBeginUser;
2518: PetscCall(TSGetDM(ts, &sw));
2519: PetscCall(DMGetDimension(sw, &dim));
2520: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2521: PetscCall(MatGetOwnershipRange(J, &rStart, NULL));
2522: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
2523: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&vel));
2524: Np /= 2 * dim;
2525: for (p = 0; p < Np; ++p) {
2526: // TODO This is not right because dv/dx has the electric field in it
2527: PetscScalar vals[4] = {0., 1., -1., 0.};
2529: for (d = 0; d < dim; ++d) {
2530: const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
2531: PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES));
2532: }
2533: }
2534: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
2535: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&vel));
2536: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
2537: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
2538: PetscFunctionReturn(PETSC_SUCCESS);
2539: }
2541: static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, void *Ctx)
2542: {
2543: AppCtx *ctx = (AppCtx *)Ctx;
2544: DM sw;
2545: const PetscScalar *v;
2546: PetscScalar *xres;
2547: PetscInt Np, p, d, dim;
2549: PetscFunctionBeginUser;
2550: PetscCall(PetscLogEventBegin(ctx->RhsXEvent, ts, 0, 0, 0));
2551: PetscCall(TSGetDM(ts, &sw));
2552: PetscCall(DMGetDimension(sw, &dim));
2553: PetscCall(VecGetLocalSize(Xres, &Np));
2554: PetscCall(VecGetArrayRead(V, &v));
2555: PetscCall(VecGetArray(Xres, &xres));
2556: Np /= dim;
2557: for (p = 0; p < Np; ++p) {
2558: for (d = 0; d < dim; ++d) xres[p * dim + d] = v[p * dim + d];
2559: }
2560: PetscCall(VecRestoreArrayRead(V, &v));
2561: PetscCall(VecRestoreArray(Xres, &xres));
2562: PetscCall(PetscLogEventEnd(ctx->RhsXEvent, ts, 0, 0, 0));
2563: PetscFunctionReturn(PETSC_SUCCESS);
2564: }
2566: static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, void *Ctx)
2567: {
2568: DM sw;
2569: AppCtx *ctx = (AppCtx *)Ctx;
2570: SNES snes = ((AppCtx *)ctx)->snes;
2571: const PetscScalar *x;
2572: PetscScalar *vres;
2573: PetscReal *E, m_p, q_p;
2574: PetscInt Np, p, dim, d;
2575: Parameter *param;
2577: PetscFunctionBeginUser;
2578: PetscCall(PetscLogEventBegin(ctx->RhsVEvent, ts, 0, 0, 0));
2579: PetscCall(TSGetDM(ts, &sw));
2580: PetscCall(ComputeFieldAtParticles(snes, sw));
2582: PetscCall(DMGetDimension(sw, &dim));
2583: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
2584: PetscCall(PetscBagGetData(ctx->bag, ¶m));
2585: m_p = ctx->masses[0] * param->m0;
2586: q_p = ctx->charges[0] * param->q0;
2587: PetscCall(VecGetLocalSize(Vres, &Np));
2588: PetscCall(VecGetArrayRead(X, &x));
2589: PetscCall(VecGetArray(Vres, &vres));
2590: Np /= dim;
2591: for (p = 0; p < Np; ++p) {
2592: for (d = 0; d < dim; ++d) vres[p * dim + d] = q_p * E[p * dim + d] / m_p;
2593: }
2594: PetscCall(VecRestoreArrayRead(X, &x));
2595: /*
2596: Synchronized, ordered output for parallel/sequential test cases.
2597: In the 1D (on the 2D mesh) case, every y component should be zero.
2598: */
2599: if (ctx->checkVRes) {
2600: PetscBool pr = ctx->checkVRes > 1 ? PETSC_TRUE : PETSC_FALSE;
2601: PetscInt step;
2603: PetscCall(TSGetStepNumber(ts, &step));
2604: if (pr) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "step: %" PetscInt_FMT "\n", step));
2605: for (PetscInt p = 0; p < Np; ++p) {
2606: if (pr) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Residual: %.12g %.12g\n", (double)PetscRealPart(vres[p * dim + 0]), (double)PetscRealPart(vres[p * dim + 1])));
2607: 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]));
2608: }
2609: if (pr) PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
2610: }
2611: PetscCall(VecRestoreArray(Vres, &vres));
2612: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
2613: PetscCall(PetscLogEventEnd(ctx->RhsVEvent, ts, 0, 0, 0));
2614: PetscFunctionReturn(PETSC_SUCCESS);
2615: }
2617: /* Discrete Gradients Formulation: S, F, gradF (G) */
2618: PetscErrorCode RHSJacobianS(TS ts, PetscReal t, Vec U, Mat S, PetscCtx ctx)
2619: {
2620: PetscScalar vals[4] = {0., 1., -1., 0.};
2621: DM sw;
2622: PetscInt dim, d, Np, p, rStart;
2624: PetscFunctionBeginUser;
2625: PetscCall(TSGetDM(ts, &sw));
2626: PetscCall(DMGetDimension(sw, &dim));
2627: PetscCall(VecGetLocalSize(U, &Np));
2628: PetscCall(MatGetOwnershipRange(S, &rStart, NULL));
2629: Np /= 2 * dim;
2630: for (p = 0; p < Np; ++p) {
2631: for (d = 0; d < dim; ++d) {
2632: const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
2633: PetscCall(MatSetValues(S, 2, rows, 2, rows, vals, INSERT_VALUES));
2634: }
2635: }
2636: PetscCall(MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY));
2637: PetscCall(MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY));
2638: PetscFunctionReturn(PETSC_SUCCESS);
2639: }
2641: PetscErrorCode RHSObjectiveF(TS ts, PetscReal t, Vec U, PetscScalar *F, void *Ctx)
2642: {
2643: AppCtx *ctx = (AppCtx *)Ctx;
2644: DM sw;
2645: Vec phi;
2646: const PetscScalar *u;
2647: PetscInt dim, Np, cStart, cEnd;
2648: PetscReal *vel, *coords, m_p = 1.;
2650: PetscFunctionBeginUser;
2651: PetscCall(TSGetDM(ts, &sw));
2652: PetscCall(DMGetDimension(sw, &dim));
2653: PetscCall(DMPlexGetHeightStratum(ctx->dmPot, 0, &cStart, &cEnd));
2655: PetscCall(DMGetNamedGlobalVector(ctx->dmPot, "phi", &phi));
2656: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view_dg"));
2657: PetscCall(computeFieldEnergy(ctx->dmPot, phi, F));
2658: PetscCall(DMRestoreNamedGlobalVector(ctx->dmPot, "phi", &phi));
2660: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
2661: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&vel));
2662: PetscCall(DMSwarmSortGetAccess(sw));
2663: PetscCall(VecGetArrayRead(U, &u));
2664: PetscCall(VecGetLocalSize(U, &Np));
2665: Np /= 2 * dim;
2666: for (PetscInt c = cStart; c < cEnd; ++c) {
2667: PetscInt *points;
2668: PetscInt Ncp;
2670: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
2671: for (PetscInt cp = 0; cp < Ncp; ++cp) {
2672: const PetscInt p = points[cp];
2673: const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);
2675: *F += 0.5 * m_p * v2;
2676: }
2677: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
2678: }
2679: PetscCall(VecRestoreArrayRead(U, &u));
2680: PetscCall(DMSwarmSortRestoreAccess(sw));
2681: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
2682: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&vel));
2683: PetscFunctionReturn(PETSC_SUCCESS);
2684: }
2686: /* dF/dx = q E dF/dv = v */
2687: PetscErrorCode RHSFunctionG(TS ts, PetscReal t, Vec U, Vec G, PetscCtx ctx)
2688: {
2689: DM sw;
2690: SNES snes = ((AppCtx *)ctx)->snes;
2691: const PetscReal *coords, *vel, *E;
2692: const PetscScalar *u;
2693: PetscScalar *g;
2694: PetscReal m_p = 1., q_p = -1.;
2695: PetscInt dim, d, Np, p;
2697: PetscFunctionBeginUser;
2698: PetscCall(TSGetDM(ts, &sw));
2699: PetscCall(DMGetDimension(sw, &dim));
2700: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2701: PetscCall(VecGetArrayRead(U, &u));
2702: PetscCall(VecGetArray(G, &g));
2704: PetscLogEvent COMPUTEFIELD;
2705: PetscCall(PetscLogEventRegister("COMPFIELDATPART", TS_CLASSID, &COMPUTEFIELD));
2706: PetscCall(PetscLogEventBegin(COMPUTEFIELD, 0, 0, 0, 0));
2707: PetscCall(ComputeFieldAtParticles(snes, sw));
2708: PetscCall(PetscLogEventEnd(COMPUTEFIELD, 0, 0, 0, 0));
2709: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
2710: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&vel));
2711: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
2712: for (p = 0; p < Np; ++p) {
2713: for (d = 0; d < dim; ++d) {
2714: g[(p * 2 + 0) * dim + d] = -(q_p / m_p) * E[p * dim + d];
2715: g[(p * 2 + 1) * dim + d] = m_p * u[(p * 2 + 1) * dim + d];
2716: }
2717: }
2718: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
2719: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
2720: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&vel));
2721: PetscCall(VecRestoreArrayRead(U, &u));
2722: PetscCall(VecRestoreArray(G, &g));
2723: PetscFunctionReturn(PETSC_SUCCESS);
2724: }
2726: static PetscErrorCode CreateSolution(TS ts)
2727: {
2728: DM sw;
2729: Vec u;
2730: PetscInt dim, Np;
2732: PetscFunctionBegin;
2733: PetscCall(TSGetDM(ts, &sw));
2734: PetscCall(DMGetDimension(sw, &dim));
2735: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2736: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
2737: PetscCall(VecSetBlockSize(u, dim));
2738: PetscCall(VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE));
2739: PetscCall(VecSetUp(u));
2740: PetscCall(TSSetSolution(ts, u));
2741: PetscCall(VecDestroy(&u));
2742: PetscFunctionReturn(PETSC_SUCCESS);
2743: }
2745: static PetscErrorCode SetProblem(TS ts)
2746: {
2747: AppCtx *ctx;
2748: DM sw;
2750: PetscFunctionBegin;
2751: PetscCall(TSGetDM(ts, &sw));
2752: PetscCall(DMGetApplicationContext(sw, &ctx));
2753: // Define unified system for (X, V)
2754: {
2755: Mat J;
2756: PetscInt dim, Np;
2758: PetscCall(DMGetDimension(sw, &dim));
2759: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2760: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
2761: PetscCall(MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE));
2762: PetscCall(MatSetBlockSize(J, 2 * dim));
2763: PetscCall(MatSetFromOptions(J));
2764: PetscCall(MatSetUp(J));
2765: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, ctx));
2766: PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, ctx));
2767: PetscCall(MatDestroy(&J));
2768: }
2769: /* Define split system for X and V */
2770: {
2771: Vec u;
2772: IS isx, isv, istmp;
2773: const PetscInt *idx;
2774: PetscInt dim, Np, rstart;
2776: PetscCall(TSGetSolution(ts, &u));
2777: PetscCall(DMGetDimension(sw, &dim));
2778: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2779: PetscCall(VecGetOwnershipRange(u, &rstart, NULL));
2780: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 0, 2, &istmp));
2781: PetscCall(ISGetIndices(istmp, &idx));
2782: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isx));
2783: PetscCall(ISRestoreIndices(istmp, &idx));
2784: PetscCall(ISDestroy(&istmp));
2785: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 1, 2, &istmp));
2786: PetscCall(ISGetIndices(istmp, &idx));
2787: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isv));
2788: PetscCall(ISRestoreIndices(istmp, &idx));
2789: PetscCall(ISDestroy(&istmp));
2790: PetscCall(TSRHSSplitSetIS(ts, "position", isx));
2791: PetscCall(TSRHSSplitSetIS(ts, "momentum", isv));
2792: PetscCall(ISDestroy(&isx));
2793: PetscCall(ISDestroy(&isv));
2794: PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, ctx));
2795: PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, ctx));
2796: }
2797: // Define symplectic formulation U_t = S . G, where G = grad F
2798: {
2799: PetscCall(TSDiscGradSetFormulation(ts, RHSJacobianS, RHSObjectiveF, RHSFunctionG, ctx));
2800: }
2801: PetscFunctionReturn(PETSC_SUCCESS);
2802: }
2804: static PetscErrorCode DMSwarmTSRedistribute(TS ts)
2805: {
2806: DM sw;
2807: Vec u;
2808: PetscReal t, maxt, dt;
2809: PetscInt n, maxn;
2811: PetscFunctionBegin;
2812: PetscCall(TSGetDM(ts, &sw));
2813: PetscCall(TSGetTime(ts, &t));
2814: PetscCall(TSGetMaxTime(ts, &maxt));
2815: PetscCall(TSGetTimeStep(ts, &dt));
2816: PetscCall(TSGetStepNumber(ts, &n));
2817: PetscCall(TSGetMaxSteps(ts, &maxn));
2819: PetscCall(TSReset(ts));
2820: PetscCall(TSSetDM(ts, sw));
2821: PetscCall(TSSetFromOptions(ts));
2822: PetscCall(TSSetTime(ts, t));
2823: PetscCall(TSSetMaxTime(ts, maxt));
2824: PetscCall(TSSetTimeStep(ts, dt));
2825: PetscCall(TSSetStepNumber(ts, n));
2826: PetscCall(TSSetMaxSteps(ts, maxn));
2828: PetscCall(CreateSolution(ts));
2829: PetscCall(SetProblem(ts));
2830: PetscCall(TSGetSolution(ts, &u));
2831: PetscFunctionReturn(PETSC_SUCCESS);
2832: }
2834: PetscErrorCode line(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *Ctx)
2835: {
2836: DM sw, cdm;
2837: PetscInt Np;
2838: PetscReal low[2], high[2];
2839: AppCtx *ctx = (AppCtx *)Ctx;
2841: sw = ctx->swarm;
2842: PetscCall(DMSwarmGetCellDM(sw, &cdm));
2843: // Get the bounding box so we can equally space the particles
2844: PetscCall(DMGetLocalBoundingBox(cdm, low, high));
2845: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2846: // shift it by h/2 so nothing is initialized directly on a boundary
2847: x[0] = ((high[0] - low[0]) / Np) * (p + 0.5);
2848: x[1] = 0.;
2849: return PETSC_SUCCESS;
2850: }
2852: /*
2853: InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values.
2855: Input Parameters:
2856: + ts - The TS
2857: - useInitial - Flag to also set the initial conditions to the current coordinates and velocities and setup the problem
2859: Output Parameters:
2860: . u - The initialized solution vector
2862: Level: advanced
2864: .seealso: InitializeSolve()
2865: */
2866: static PetscErrorCode InitializeSolveAndSwarm(TS ts, PetscBool useInitial)
2867: {
2868: DM sw;
2869: Vec u, gc, gv;
2870: IS isx, isv;
2871: PetscInt dim;
2872: AppCtx *ctx;
2874: PetscFunctionBeginUser;
2875: PetscCall(TSGetDM(ts, &sw));
2876: PetscCall(DMGetApplicationContext(sw, &ctx));
2877: PetscCall(DMGetDimension(sw, &dim));
2878: if (useInitial) {
2879: PetscReal v0[2] = {1., 0.};
2880: if (ctx->perturbed_weights) {
2881: PetscCall(InitializeParticles_PerturbedWeights(sw, ctx));
2882: } else {
2883: PetscCall(DMSwarmComputeLocalSizeFromOptions(sw));
2884: PetscCall(DMSwarmInitializeCoordinates(sw));
2885: PetscCall(DMSwarmInitializeVelocitiesFromOptions(sw, v0));
2886: }
2887: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
2888: PetscCall(DMSwarmTSRedistribute(ts));
2889: }
2890: PetscCall(DMSetUp(sw));
2891: PetscCall(TSGetSolution(ts, &u));
2892: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
2893: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
2894: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
2895: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
2896: PetscCall(VecISCopy(u, isx, SCATTER_FORWARD, gc));
2897: PetscCall(VecISCopy(u, isv, SCATTER_FORWARD, gv));
2898: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
2899: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
2900: PetscFunctionReturn(PETSC_SUCCESS);
2901: }
2903: static PetscErrorCode InitializeSolve(TS ts, Vec u)
2904: {
2905: PetscFunctionBegin;
2906: PetscCall(TSSetSolution(ts, u));
2907: PetscCall(InitializeSolveAndSwarm(ts, PETSC_TRUE));
2908: PetscFunctionReturn(PETSC_SUCCESS);
2909: }
2911: static PetscErrorCode MigrateParticles(TS ts)
2912: {
2913: DM sw, cdm;
2914: const PetscReal *L;
2915: AppCtx *ctx;
2917: PetscFunctionBeginUser;
2918: PetscCall(TSGetDM(ts, &sw));
2919: PetscCall(DMGetApplicationContext(sw, &ctx));
2920: PetscCall(DMViewFromOptions(sw, NULL, "-migrate_view_pre"));
2921: {
2922: Vec u, gc, gv, position, momentum;
2923: IS isx, isv;
2924: PetscReal *pos, *mom;
2926: PetscCall(TSGetSolution(ts, &u));
2927: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
2928: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
2929: PetscCall(VecGetSubVector(u, isx, &position));
2930: PetscCall(VecGetSubVector(u, isv, &momentum));
2931: PetscCall(VecGetArray(position, &pos));
2932: PetscCall(VecGetArray(momentum, &mom));
2933: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
2934: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
2935: PetscCall(VecISCopy(u, isx, SCATTER_REVERSE, gc));
2936: PetscCall(VecISCopy(u, isv, SCATTER_REVERSE, gv));
2938: PetscCall(DMSwarmGetCellDM(sw, &cdm));
2939: PetscCall(DMGetPeriodicity(cdm, NULL, NULL, &L));
2940: PetscCheck(L, PetscObjectComm((PetscObject)cdm), PETSC_ERR_ARG_WRONG, "Mesh must be periodic");
2941: if ((L[0] || L[1]) >= 0.) {
2942: PetscReal *x, *v, upper[3], lower[3];
2943: PetscInt Np, dim;
2945: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2946: PetscCall(DMGetDimension(cdm, &dim));
2947: PetscCall(DMGetBoundingBox(cdm, lower, upper));
2948: PetscCall(VecGetArray(gc, &x));
2949: PetscCall(VecGetArray(gv, &v));
2950: for (PetscInt p = 0; p < Np; ++p) {
2951: for (PetscInt d = 0; d < dim; ++d) {
2952: if (pos[p * dim + d] < lower[d]) {
2953: x[p * dim + d] = pos[p * dim + d] + (upper[d] - lower[d]);
2954: } else if (pos[p * dim + d] > upper[d]) {
2955: x[p * dim + d] = pos[p * dim + d] - (upper[d] - lower[d]);
2956: } else {
2957: x[p * dim + d] = pos[p * dim + d];
2958: }
2959: 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]);
2960: v[p * dim + d] = mom[p * dim + d];
2961: }
2962: }
2963: PetscCall(VecRestoreArray(gc, &x));
2964: PetscCall(VecRestoreArray(gv, &v));
2965: }
2966: PetscCall(VecRestoreArray(position, &pos));
2967: PetscCall(VecRestoreArray(momentum, &mom));
2968: PetscCall(VecRestoreSubVector(u, isx, &position));
2969: PetscCall(VecRestoreSubVector(u, isv, &momentum));
2970: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
2971: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
2972: }
2973: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
2974: PetscInt step;
2976: PetscCall(TSGetStepNumber(ts, &step));
2977: if (!(step % ctx->remapFreq)) {
2978: // Monitor electric field before we destroy it
2979: PetscReal ptime;
2980: PetscInt step;
2982: PetscCall(TSGetStepNumber(ts, &step));
2983: PetscCall(TSGetTime(ts, &ptime));
2984: if (ctx->efield_monitor) PetscCall(MonitorEField(ts, step, ptime, NULL, ctx));
2985: if (ctx->poisson_monitor) PetscCall(MonitorPoisson(ts, step, ptime, NULL, ctx));
2986: PetscCall(DMSwarmRemap(sw));
2987: ctx->validE = PETSC_FALSE;
2988: }
2989: // This MUST come last, since it recreates the subswarms and they must DMClone() the new swarm
2990: PetscCall(DMSwarmTSRedistribute(ts));
2991: PetscCall(InitializeSolveAndSwarm(ts, PETSC_FALSE));
2992: PetscFunctionReturn(PETSC_SUCCESS);
2993: }
2995: int main(int argc, char **argv)
2996: {
2997: DM dm, sw;
2998: TS ts;
2999: Vec u;
3000: PetscReal dt;
3001: PetscInt maxn;
3002: AppCtx ctx;
3004: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
3005: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
3006: PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &ctx.bag));
3007: PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
3008: PetscCall(CreatePoisson(dm, &ctx));
3009: PetscCall(CreateMomentFields(dm, &ctx));
3010: PetscCall(CreateSwarm(dm, &ctx, &sw));
3011: PetscCall(SetupParameters(PETSC_COMM_WORLD, &ctx));
3012: PetscCall(InitializeConstants(sw, &ctx));
3013: PetscCall(DMSetApplicationContext(sw, &ctx));
3015: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
3016: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
3017: PetscCall(TSSetDM(ts, sw));
3018: PetscCall(TSSetMaxTime(ts, 0.1));
3019: PetscCall(TSSetTimeStep(ts, 0.00001));
3020: PetscCall(TSSetMaxSteps(ts, 100));
3021: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
3023: if (ctx.efield_monitor) PetscCall(TSMonitorSet(ts, MonitorEField, &ctx, NULL));
3024: if (ctx.moment_monitor) PetscCall(TSMonitorSet(ts, MonitorMoments, &ctx, NULL));
3025: if (ctx.moment_field_monitor) PetscCall(TSMonitorSet(ts, MonitorMomentFields, &ctx, NULL));
3026: if (ctx.initial_monitor) PetscCall(TSMonitorSet(ts, MonitorInitialConditions, &ctx, NULL));
3027: if (ctx.positions_monitor) PetscCall(TSMonitorSet(ts, MonitorPositions_2D, &ctx, NULL));
3028: if (ctx.poisson_monitor) PetscCall(TSMonitorSet(ts, MonitorPoisson, &ctx, NULL));
3029: if (ctx.velocity_monitor >= 0) PetscCall(TSMonitorSet(ts, MonitorVelocity, &ctx, NULL));
3031: PetscCall(TSSetFromOptions(ts));
3032: PetscCall(TSGetTimeStep(ts, &dt));
3033: PetscCall(TSGetMaxSteps(ts, &maxn));
3034: ctx.steps = maxn;
3035: ctx.stepSize = dt;
3036: PetscCall(SetupContext(dm, sw, &ctx));
3037: PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
3038: PetscCall(TSSetPostStep(ts, MigrateParticles));
3039: PetscCall(CreateSolution(ts));
3040: PetscCall(TSGetSolution(ts, &u));
3041: PetscCall(TSComputeInitialCondition(ts, u));
3042: PetscCall(CheckNonNegativeWeights(sw, &ctx));
3043: PetscCall(TSSolve(ts, NULL));
3045: if (ctx.checkLandau) {
3046: // We should get a lookup table based on charge density and \hat k
3047: const PetscReal gammaEx = -0.15336;
3048: const PetscReal omegaEx = 1.4156;
3049: const PetscReal tol = 1e-2;
3051: PetscCheck(PetscAbsReal((ctx.gamma - gammaEx) / gammaEx) < tol, PETSC_COMM_WORLD, PETSC_ERR_LIB, "Invalid Landau gamma %g != %g", ctx.gamma, gammaEx);
3052: PetscCheck(PetscAbsReal((ctx.omega - omegaEx) / omegaEx) < tol, PETSC_COMM_WORLD, PETSC_ERR_LIB, "Invalid Landau omega %g != %g", ctx.omega, omegaEx);
3053: }
3055: PetscCall(SNESDestroy(&ctx.snes));
3056: PetscCall(DMDestroy(&ctx.dmN));
3057: PetscCall(ISDestroy(&ctx.isN));
3058: PetscCall(MatDestroy(&ctx.MN));
3059: PetscCall(DMDestroy(&ctx.dmP));
3060: PetscCall(ISDestroy(&ctx.isP));
3061: PetscCall(MatDestroy(&ctx.MP));
3062: PetscCall(DMDestroy(&ctx.dmE));
3063: PetscCall(ISDestroy(&ctx.isE));
3064: PetscCall(MatDestroy(&ctx.ME));
3065: PetscCall(DMDestroy(&ctx.dmMom));
3066: PetscCall(DMDestroy(&ctx.dmPot));
3067: PetscCall(ISDestroy(&ctx.isPot));
3068: PetscCall(MatDestroy(&ctx.M));
3069: PetscCall(PetscFEGeomDestroy(&ctx.fegeom));
3070: PetscCall(TSDestroy(&ts));
3071: PetscCall(DMDestroy(&sw));
3072: PetscCall(DMDestroy(&dm));
3073: PetscCall(DestroyContext(&ctx));
3074: PetscCall(PetscFinalize());
3075: return 0;
3076: }
3078: /*TEST
3080: build:
3081: requires: !complex double
3083: # This tests that we can compute the correct decay rate and frequency
3084: # For gold runs, use -dm_plex_box_faces 160 -vdm_plex_box_faces 450 -remap_dm_plex_box_faces 80,150 -ts_max_steps 1000
3085: # -remap_freq 100 -emax_start_step 50 -emax_solve_step 100
3086: testset:
3087: args: -cosine_coefficients 0.01 -charges -1. -perturbed_weights -total_weight 1. \
3088: -dm_plex_dim 1 -dm_plex_box_faces 80 -dm_plex_box_lower 0. -dm_plex_box_upper 12.5664 \
3089: -dm_plex_box_bd periodic -dm_plex_hash_location \
3090: -vdm_plex_dim 1 -vdm_plex_box_faces 220 -vdm_plex_box_lower -6 -vdm_plex_box_upper 6 \
3091: -vpetscspace_degree 2 -vdm_plex_hash_location \
3092: -remap_freq 1 -dm_swarm_remap_type pfak -remap_dm_plex_dim 2 -remap_dm_plex_simplex 0 \
3093: -remap_dm_plex_box_faces 40,110 -remap_dm_plex_box_bd periodic,none \
3094: -remap_dm_plex_box_lower 0.,-6. -remap_dm_plex_box_upper 12.5664,6. \
3095: -remap_petscspace_degree 1 -remap_dm_plex_hash_location \
3096: -ftop_ksp_type lsqr -ftop_pc_type none -ftop_ksp_rtol 1.e-14 -ptof_pc_type lu \
3097: -em_type primal -petscspace_degree 1 -em_snes_atol 1.e-12 -em_snes_error_if_not_converged \
3098: -em_ksp_error_if_not_converged -em_pc_type svd -em_proj_pc_type lu \
3099: -ts_time_step 0.03 -ts_max_steps 2 -ts_max_time 100 \
3100: -emax_tao_type brgn -emax_tao_max_it 100 -emax_tao_brgn_regularization_type l2pure \
3101: -emax_tao_brgn_regularizer_weight 1e-5 -emax_tao_brgn_subsolver_tao_bnk_ksp_rtol 1e-12 \
3102: -emax_start_step 1 -emax_solve_step 1 \
3103: -output_step 1 -efield_monitor quiet
3105: test:
3106: suffix: landau_damping_1d_bs
3107: args: -ts_type basicsymplectic -ts_basicsymplectic_type 1
3109: test:
3110: suffix: landau_damping_1d_dg
3111: args: -ts_type discgrad -ts_discgrad_type average -snes_type qn
3113: TEST*/