Actual source code: ex2.c
1: static char help[] = "Landau Damping/Two Stream instability test using Vlasov-Poisson equations\n";
3: /*
4: TODO:
5: - Cache mesh geometry
6: - Move electrostatic solver to MG+SVD
8: 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"
9: According to Lukas, good damping results come at ~16k particles per cell
11: To visualize the maximum electric field use
13: -efield_monitor
15: To monitor velocity moments of the distribution use
17: -ptof_pc_type lu -moments_monitor
19: To monitor the particle positions in phase space use
21: -positions_monitor
23: To monitor the charge density, E field, and potential use
25: -poisson_monitor
27: To monitor the remapping field use
29: -remap_uf_view draw
31: To visualize the swarm distribution use
33: -ts_monitor_hg_swarm
35: To visualize the particles, we can use
37: -ts_monitor_sp_swarm -ts_monitor_sp_swarm_retain 0 -ts_monitor_sp_swarm_phase 1 -draw_size 500,500
39: For a Landau Damping verification run, we use
41: # Physics
42: -cosine_coefficients 0.01 -dm_swarm_num_species 1 -charges -1. -perturbed_weights -total_weight 1.
43: # Spatial Mesh
44: -dm_plex_dim 1 -dm_plex_box_faces 40 -dm_plex_box_lower 0. -dm_plex_box_upper 12.5664 -dm_plex_box_bd periodic -dm_plex_hash_location
45: # Velocity Mesh
46: -vdm_plex_dim 1 -vdm_plex_box_faces 40 -vdm_plex_box_lower -6 -vdm_plex_box_upper 6 -vpetscspace_degree 2 -vdm_plex_hash_location
47: # Remap Space
48: -dm_swarm_remap_type pfak -remap_freq 100
49: -remap_dm_plex_dim 2 -remap_dm_plex_simplex 0 -remap_dm_plex_box_faces 20,20 -remap_dm_plex_box_bd periodic,none -remap_dm_plex_box_lower 0.,-6.
50: -remap_dm_plex_box_upper 12.5664,6. -remap_petscspace_degree 1 -remap_dm_plex_hash_location
51: # Remap Solve
52: -ftop_ksp_type lsqr -ftop_pc_type none -ftop_ksp_rtol 1.e-14 -ptof_pc_type lu
53: # EM Solve
54: -em_type primal -petscspace_degree 1 -em_snes_atol 1.e-12 -em_snes_error_if_not_converged -em_ksp_error_if_not_converged -em_pc_type svd -em_proj_pc_type lu
55: # Timestepping
56: -ts_type basicsymplectic -ts_basicsymplectic_type 1 -ts_dt 0.03 -ts_max_steps 1500 -ts_max_time 100
57: # Monitoring
58: -output_step 1 -check_vel_res -efield_monitor -poisson_monitor -positions_monitor -dm_swarm_print_coords 0 -remap_uf_view draw
59: -ftop_ksp_lsqr_monitor -ftop_ksp_converged_reason
61: */
62: #include <petscts.h>
63: #include <petscdmplex.h>
64: #include <petscdmswarm.h>
65: #include <petscfe.h>
66: #include <petscds.h>
67: #include <petscbag.h>
68: #include <petscdraw.h>
69: #include <petsc/private/dmpleximpl.h>
70: #include <petsc/private/petscfeimpl.h>
71: #include <petsc/private/dmswarmimpl.h>
72: #include "petscdm.h"
73: #include "petscdmlabel.h"
75: PETSC_EXTERN PetscErrorCode stream(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
76: PETSC_EXTERN PetscErrorCode line(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
78: const char *EMTypes[] = {"primal", "mixed", "coulomb", "none", "EMType", "EM_", NULL};
79: typedef enum {
80: EM_PRIMAL,
81: EM_MIXED,
82: EM_COULOMB,
83: EM_NONE
84: } EMType;
86: typedef enum {
87: V0,
88: X0,
89: T0,
90: M0,
91: Q0,
92: PHI0,
93: POISSON,
94: VLASOV,
95: SIGMA,
96: NUM_CONSTANTS
97: } ConstantType;
98: typedef struct {
99: PetscScalar v0; /* Velocity scale, often the thermal velocity */
100: PetscScalar t0; /* Time scale */
101: PetscScalar x0; /* Space scale */
102: PetscScalar m0; /* Mass scale */
103: PetscScalar q0; /* Charge scale */
104: PetscScalar kb;
105: PetscScalar epsi0;
106: PetscScalar phi0; /* Potential scale */
107: PetscScalar poissonNumber; /* Non-Dimensional Poisson Number */
108: PetscScalar vlasovNumber; /* Non-Dimensional Vlasov Number */
109: PetscReal sigma; /* Nondimensional charge per length in x */
110: } Parameter;
112: typedef struct {
113: PetscBag bag; // Problem parameters
114: PetscBool error; // Flag for printing the error
115: PetscInt remapFreq; // Number of timesteps between remapping
116: PetscBool efield_monitor; // Flag to show electric field monitor
117: PetscBool moment_monitor; // Flag to show distribution moment monitor
118: PetscBool positions_monitor; // Flag to show particle positins at each time step
119: PetscBool poisson_monitor; // Flag to display charge, E field, and potential at each solve
120: PetscBool initial_monitor; // Flag to monitor the initial conditions
121: PetscInt velocity_monitor; // Cell to monitor the velocity distribution for
122: PetscBool perturbed_weights; // Uniformly sample x,v space with gaussian weights
123: PetscInt ostep; // Print the energy at each ostep time steps
124: PetscInt numParticles;
125: PetscReal timeScale; /* Nondimensionalizing time scale */
126: PetscReal charges[2]; /* The charges of each species */
127: PetscReal masses[2]; /* The masses of each species */
128: PetscReal thermal_energy[2]; /* Thermal Energy (used to get other constants)*/
129: PetscReal cosine_coefficients[2]; /*(alpha, k)*/
130: PetscReal totalWeight;
131: PetscReal stepSize;
132: PetscInt steps;
133: PetscReal initVel;
134: EMType em; // Type of electrostatic model
135: SNES snes; // EM solver
136: DM dmPot; // The DM for potential
137: Mat fftPot; // Fourier Transform operator for the potential
138: Vec fftX, fftY; // FFT vectors with phases added (complex parts)
139: IS fftReal; // The indices for real parts
140: IS isPot; // The IS for potential, or NULL in primal
141: Mat M; // The finite element mass matrix for potential
142: PetscFEGeom *fegeom; // Geometric information for the DM cells
143: PetscDrawHG drawhgic_x; // Histogram of the particle weight in each X cell
144: PetscDrawHG drawhgic_v; // Histogram of the particle weight in each X cell
145: PetscDrawHG drawhgcell_v; // Histogram of the particle weight in a given cell
146: PetscBool validE; // Flag to indicate E-field in swarm is valid
147: PetscReal drawlgEmin; // The minimum lg(E) to plot
148: PetscDrawLG drawlgE; // Logarithm of maximum electric field
149: PetscDrawSP drawspE; // Electric field at particle positions
150: PetscDrawSP drawspX; // Particle positions
151: PetscViewer viewerRho; // Charge density viewer
152: PetscViewer viewerRhoHat; // Charge density Fourier Transform viewer
153: PetscViewer viewerPhi; // Potential viewer
154: DM swarm;
155: PetscRandom random;
156: PetscBool twostream;
157: PetscBool checkweights;
158: PetscInt checkVRes; // Flag to check/output velocity residuals for nightly tests
160: PetscLogEvent RhsXEvent, RhsVEvent, ESolveEvent, ETabEvent;
161: } AppCtx;
163: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
164: {
165: PetscFunctionBeginUser;
166: PetscInt d = 2;
167: PetscInt maxSpecies = 2;
168: options->error = PETSC_FALSE;
169: options->remapFreq = 1;
170: options->efield_monitor = PETSC_FALSE;
171: options->moment_monitor = PETSC_FALSE;
172: options->initial_monitor = PETSC_FALSE;
173: options->perturbed_weights = PETSC_FALSE;
174: options->poisson_monitor = PETSC_FALSE;
175: options->positions_monitor = PETSC_FALSE;
176: options->velocity_monitor = -1;
177: options->ostep = 100;
178: options->timeScale = 2.0e-14;
179: options->charges[0] = -1.0;
180: options->charges[1] = 1.0;
181: options->masses[0] = 1.0;
182: options->masses[1] = 1000.0;
183: options->thermal_energy[0] = 1.0;
184: options->thermal_energy[1] = 1.0;
185: options->cosine_coefficients[0] = 0.01;
186: options->cosine_coefficients[1] = 0.5;
187: options->initVel = 1;
188: options->totalWeight = 1.0;
189: options->drawhgic_x = NULL;
190: options->drawhgic_v = NULL;
191: options->drawhgcell_v = NULL;
192: options->drawlgEmin = -6;
193: options->drawlgE = NULL;
194: options->drawspE = NULL;
195: options->drawspX = NULL;
196: options->viewerRho = NULL;
197: options->viewerRhoHat = NULL;
198: options->viewerPhi = NULL;
199: options->em = EM_COULOMB;
200: options->snes = NULL;
201: options->dmPot = NULL;
202: options->fftPot = NULL;
203: options->fftX = NULL;
204: options->fftY = NULL;
205: options->fftReal = NULL;
206: options->isPot = NULL;
207: options->M = NULL;
208: options->numParticles = 32768;
209: options->twostream = PETSC_FALSE;
210: options->checkweights = PETSC_FALSE;
211: options->checkVRes = 0;
213: PetscOptionsBegin(comm, "", "Landau Damping and Two Stream options", "DMSWARM");
214: PetscCall(PetscOptionsBool("-error", "Flag to print the error", "ex2.c", options->error, &options->error, NULL));
215: PetscCall(PetscOptionsInt("-remap_freq", "Number", "ex2.c", options->remapFreq, &options->remapFreq, NULL));
216: PetscCall(PetscOptionsBool("-efield_monitor", "Flag to plot log(max E) over time", "ex2.c", options->efield_monitor, &options->efield_monitor, NULL));
217: PetscCall(PetscOptionsReal("-efield_min_monitor", "Minimum E field to plot", "ex2.c", options->drawlgEmin, &options->drawlgEmin, NULL));
218: PetscCall(PetscOptionsBool("-moments_monitor", "Flag to show moments table", "ex2.c", options->moment_monitor, &options->moment_monitor, NULL));
219: PetscCall(PetscOptionsBool("-ics_monitor", "Flag to show initial condition histograms", "ex2.c", options->initial_monitor, &options->initial_monitor, NULL));
220: PetscCall(PetscOptionsBool("-positions_monitor", "The flag to show particle positions", "ex2.c", options->positions_monitor, &options->positions_monitor, NULL));
221: PetscCall(PetscOptionsBool("-poisson_monitor", "The flag to show charges, Efield and potential solve", "ex2.c", options->poisson_monitor, &options->poisson_monitor, NULL));
222: PetscCall(PetscOptionsInt("-velocity_monitor", "Cell to show velocity histograms", "ex2.c", options->velocity_monitor, &options->velocity_monitor, NULL));
223: PetscCall(PetscOptionsBool("-twostream", "Run two stream instability", "ex2.c", options->twostream, &options->twostream, NULL));
224: PetscCall(PetscOptionsBool("-perturbed_weights", "Flag to run uniform sampling with perturbed weights", "ex2.c", options->perturbed_weights, &options->perturbed_weights, NULL));
225: PetscCall(PetscOptionsBool("-check_weights", "Ensure all particle weights are positive", "ex2.c", options->checkweights, &options->checkweights, NULL));
226: PetscCall(PetscOptionsInt("-check_vel_res", "Check particle velocity residuals for nightly tests", "ex2.c", options->checkVRes, &options->checkVRes, NULL));
227: PetscCall(PetscOptionsInt("-output_step", "Number of time steps between output", "ex2.c", options->ostep, &options->ostep, NULL));
228: PetscCall(PetscOptionsReal("-timeScale", "Nondimensionalizing time scale", "ex2.c", options->timeScale, &options->timeScale, NULL));
229: PetscCall(PetscOptionsReal("-initial_velocity", "Initial velocity of perturbed particle", "ex2.c", options->initVel, &options->initVel, NULL));
230: PetscCall(PetscOptionsReal("-total_weight", "Total weight of all particles", "ex2.c", options->totalWeight, &options->totalWeight, NULL));
231: PetscCall(PetscOptionsRealArray("-cosine_coefficients", "Amplitude and frequency of cosine equation used in initialization", "ex2.c", options->cosine_coefficients, &d, NULL));
232: PetscCall(PetscOptionsRealArray("-charges", "Species charges", "ex2.c", options->charges, &maxSpecies, NULL));
233: PetscCall(PetscOptionsEnum("-em_type", "Type of electrostatic solver", "ex2.c", EMTypes, (PetscEnum)options->em, (PetscEnum *)&options->em, NULL));
234: PetscOptionsEnd();
236: PetscCall(PetscLogEventRegister("RhsX", TS_CLASSID, &options->RhsXEvent));
237: PetscCall(PetscLogEventRegister("RhsV", TS_CLASSID, &options->RhsVEvent));
238: PetscCall(PetscLogEventRegister("ESolve", TS_CLASSID, &options->ESolveEvent));
239: PetscCall(PetscLogEventRegister("ETab", TS_CLASSID, &options->ETabEvent));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: static PetscErrorCode SetupContext(DM dm, DM sw, AppCtx *user)
244: {
245: MPI_Comm comm;
247: PetscFunctionBeginUser;
248: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
249: if (user->efield_monitor) {
250: PetscDraw draw;
251: PetscDrawAxis axis;
253: PetscCall(PetscDrawCreate(comm, NULL, "Max Electric Field", 0, 300, 400, 300, &draw));
254: PetscCall(PetscDrawSetSave(draw, "ex2_Efield"));
255: PetscCall(PetscDrawSetFromOptions(draw));
256: PetscCall(PetscDrawLGCreate(draw, 1, &user->drawlgE));
257: PetscCall(PetscDrawDestroy(&draw));
258: PetscCall(PetscDrawLGGetAxis(user->drawlgE, &axis));
259: PetscCall(PetscDrawAxisSetLabels(axis, "Electron Electric Field", "time", "E_max"));
260: PetscCall(PetscDrawLGSetLimits(user->drawlgE, 0., user->steps * user->stepSize, user->drawlgEmin, 0.));
261: }
263: if (user->initial_monitor) {
264: PetscDraw drawic_x, drawic_v;
265: PetscDrawAxis axis1, axis2;
266: PetscReal dmboxlower[2], dmboxupper[2];
267: PetscInt dim, cStart, cEnd;
269: PetscCall(DMGetDimension(sw, &dim));
270: PetscCall(DMGetBoundingBox(dm, dmboxlower, dmboxupper));
271: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
273: PetscCall(PetscDrawCreate(comm, NULL, "monitor_initial_conditions_x", 0, 300, 400, 300, &drawic_x));
274: PetscCall(PetscDrawSetSave(drawic_x, "ex2_ic_x"));
275: PetscCall(PetscDrawSetFromOptions(drawic_x));
276: PetscCall(PetscDrawHGCreate(drawic_x, (int)dim, &user->drawhgic_x));
277: PetscCall(PetscDrawHGCalcStats(user->drawhgic_x, PETSC_TRUE));
278: PetscCall(PetscDrawHGGetAxis(user->drawhgic_x, &axis1));
279: PetscCall(PetscDrawHGSetNumberBins(user->drawhgic_x, (int)(cEnd - cStart)));
280: PetscCall(PetscDrawAxisSetLabels(axis1, "Initial X Distribution", "X", "weight"));
281: PetscCall(PetscDrawAxisSetLimits(axis1, dmboxlower[0], dmboxupper[0], 0, 0));
282: PetscCall(PetscDrawDestroy(&drawic_x));
284: PetscCall(PetscDrawCreate(comm, NULL, "monitor_initial_conditions_v", 400, 300, 400, 300, &drawic_v));
285: PetscCall(PetscDrawSetSave(drawic_v, "ex9_ic_v"));
286: PetscCall(PetscDrawSetFromOptions(drawic_v));
287: PetscCall(PetscDrawHGCreate(drawic_v, (int)dim, &user->drawhgic_v));
288: PetscCall(PetscDrawHGCalcStats(user->drawhgic_v, PETSC_TRUE));
289: PetscCall(PetscDrawHGGetAxis(user->drawhgic_v, &axis2));
290: PetscCall(PetscDrawHGSetNumberBins(user->drawhgic_v, 21));
291: PetscCall(PetscDrawAxisSetLabels(axis2, "Initial V_x Distribution", "V", "weight"));
292: PetscCall(PetscDrawAxisSetLimits(axis2, -6, 6, 0, 0));
293: PetscCall(PetscDrawDestroy(&drawic_v));
294: }
296: if (user->velocity_monitor >= 0) {
297: DM vdm;
298: DMSwarmCellDM celldm;
299: PetscDraw drawcell_v;
300: PetscDrawAxis axis;
301: PetscReal dmboxlower[2], dmboxupper[2];
302: PetscInt dim;
303: char title[PETSC_MAX_PATH_LEN];
305: PetscCall(DMSwarmGetCellDMByName(sw, "velocity", &celldm));
306: PetscCall(DMSwarmCellDMGetDM(celldm, &vdm));
307: PetscCall(DMGetDimension(vdm, &dim));
308: PetscCall(DMGetBoundingBox(vdm, dmboxlower, dmboxupper));
310: PetscCall(PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Cell %" PetscInt_FMT ": Velocity Distribution", user->velocity_monitor));
311: PetscCall(PetscDrawCreate(comm, NULL, title, 400, 300, 400, 300, &drawcell_v));
312: PetscCall(PetscDrawSetSave(drawcell_v, "ex2_cell_v"));
313: PetscCall(PetscDrawSetFromOptions(drawcell_v));
314: PetscCall(PetscDrawHGCreate(drawcell_v, (int)dim, &user->drawhgcell_v));
315: PetscCall(PetscDrawHGCalcStats(user->drawhgcell_v, PETSC_TRUE));
316: PetscCall(PetscDrawHGGetAxis(user->drawhgcell_v, &axis));
317: PetscCall(PetscDrawHGSetNumberBins(user->drawhgcell_v, 21));
318: PetscCall(PetscDrawAxisSetLabels(axis, "V_x Distribution", "V", "weight"));
319: PetscCall(PetscDrawAxisSetLimits(axis, dmboxlower[0], dmboxupper[0], 0, 0));
320: PetscCall(PetscDrawDestroy(&drawcell_v));
321: }
323: if (user->positions_monitor) {
324: PetscDraw draw;
325: PetscDrawAxis axis;
327: PetscCall(PetscDrawCreate(comm, NULL, "Particle Position", 0, 0, 400, 300, &draw));
328: PetscCall(PetscDrawSetSave(draw, "ex9_pos"));
329: PetscCall(PetscDrawSetFromOptions(draw));
330: PetscCall(PetscDrawSPCreate(draw, 10, &user->drawspX));
331: PetscCall(PetscDrawDestroy(&draw));
332: PetscCall(PetscDrawSPSetDimension(user->drawspX, 1));
333: PetscCall(PetscDrawSPGetAxis(user->drawspX, &axis));
334: PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "x", "v"));
335: PetscCall(PetscDrawSPReset(user->drawspX));
336: }
337: if (user->poisson_monitor) {
338: Vec rho, rhohat, phi;
339: PetscDraw draw;
340: PetscDrawAxis axis;
342: PetscCall(PetscDrawCreate(comm, NULL, "Electric_Field", 0, 0, 400, 300, &draw));
343: PetscCall(PetscDrawSetFromOptions(draw));
344: PetscCall(PetscDrawSetSave(draw, "ex9_E_spatial"));
345: PetscCall(PetscDrawSPCreate(draw, 10, &user->drawspE));
346: PetscCall(PetscDrawDestroy(&draw));
347: PetscCall(PetscDrawSPSetDimension(user->drawspE, 1));
348: PetscCall(PetscDrawSPGetAxis(user->drawspE, &axis));
349: PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "x", "E"));
350: PetscCall(PetscDrawSPReset(user->drawspE));
352: PetscCall(PetscViewerDrawOpen(comm, NULL, "Charge Density", 0, 0, 400, 300, &user->viewerRho));
353: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)user->viewerRho, "rho_"));
354: PetscCall(PetscViewerDrawGetDraw(user->viewerRho, 0, &draw));
355: PetscCall(PetscDrawSetSave(draw, "ex9_rho_spatial"));
356: PetscCall(PetscViewerSetFromOptions(user->viewerRho));
357: PetscCall(DMGetNamedGlobalVector(user->dmPot, "rho", &rho));
358: PetscCall(PetscObjectSetName((PetscObject)rho, "charge_density"));
359: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "rho", &rho));
361: PetscInt dim, N;
363: PetscCall(DMGetDimension(user->dmPot, &dim));
364: if (dim == 1) {
365: PetscCall(DMGetNamedGlobalVector(user->dmPot, "rhohat", &rhohat));
366: PetscCall(VecGetSize(rhohat, &N));
367: PetscCall(MatCreateFFT(comm, dim, &N, MATFFTW, &user->fftPot));
368: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "rhohat", &rhohat));
369: PetscCall(MatCreateVecs(user->fftPot, &user->fftX, &user->fftY));
370: PetscCall(ISCreateStride(PETSC_COMM_SELF, N, 0, 1, &user->fftReal));
371: }
373: PetscCall(PetscViewerDrawOpen(comm, NULL, "rhohat: Charge Density FT", 0, 0, 400, 300, &user->viewerRhoHat));
374: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)user->viewerRhoHat, "rhohat_"));
375: PetscCall(PetscViewerDrawGetDraw(user->viewerRhoHat, 0, &draw));
376: PetscCall(PetscDrawSetSave(draw, "ex9_rho_ft"));
377: PetscCall(PetscViewerSetFromOptions(user->viewerRhoHat));
378: PetscCall(DMGetNamedGlobalVector(user->dmPot, "rhohat", &rhohat));
379: PetscCall(PetscObjectSetName((PetscObject)rhohat, "charge_density_ft"));
380: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "rhohat", &rhohat));
382: PetscCall(PetscViewerDrawOpen(comm, NULL, "Potential", 400, 0, 400, 300, &user->viewerPhi));
383: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)user->viewerPhi, "phi_"));
384: PetscCall(PetscViewerDrawGetDraw(user->viewerPhi, 0, &draw));
385: PetscCall(PetscDrawSetSave(draw, "ex9_phi_spatial"));
386: PetscCall(PetscViewerSetFromOptions(user->viewerPhi));
387: PetscCall(DMGetNamedGlobalVector(user->dmPot, "phi", &phi));
388: PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
389: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "phi", &phi));
390: }
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: static PetscErrorCode DestroyContext(AppCtx *user)
395: {
396: PetscFunctionBeginUser;
397: PetscCall(PetscDrawHGDestroy(&user->drawhgic_x));
398: PetscCall(PetscDrawHGDestroy(&user->drawhgic_v));
399: PetscCall(PetscDrawHGDestroy(&user->drawhgcell_v));
401: PetscCall(PetscDrawLGDestroy(&user->drawlgE));
402: PetscCall(PetscDrawSPDestroy(&user->drawspE));
403: PetscCall(PetscDrawSPDestroy(&user->drawspX));
404: PetscCall(PetscViewerDestroy(&user->viewerRho));
405: PetscCall(PetscViewerDestroy(&user->viewerRhoHat));
406: PetscCall(MatDestroy(&user->fftPot));
407: PetscCall(VecDestroy(&user->fftX));
408: PetscCall(VecDestroy(&user->fftY));
409: PetscCall(ISDestroy(&user->fftReal));
410: PetscCall(PetscViewerDestroy(&user->viewerPhi));
412: PetscCall(PetscBagDestroy(&user->bag));
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: static PetscErrorCode CheckNonNegativeWeights(DM sw, AppCtx *user)
417: {
418: const PetscScalar *w;
419: PetscInt Np;
421: PetscFunctionBeginUser;
422: if (!user->checkweights) PetscFunctionReturn(PETSC_SUCCESS);
423: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&w));
424: PetscCall(DMSwarmGetLocalSize(sw, &Np));
425: 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]);
426: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
427: PetscFunctionReturn(PETSC_SUCCESS);
428: }
430: static PetscErrorCode computeVelocityFEMMoments(DM sw, PetscReal moments[], AppCtx *user)
431: {
432: DMSwarmCellDM celldm;
433: DM vdm;
434: Vec u[1];
435: const char *fields[1] = {"w_q"};
437: PetscFunctionBegin;
438: PetscCall(DMSwarmSetCellDMActive(sw, "velocity"));
439: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
440: PetscCall(DMSwarmCellDMGetDM(celldm, &vdm));
441: PetscCall(DMGetGlobalVector(vdm, &u[0]));
442: PetscCall(DMSwarmProjectFields(sw, vdm, 1, fields, u, SCATTER_FORWARD));
443: PetscCall(DMPlexComputeMoments(vdm, u[0], moments));
444: PetscCall(DMRestoreGlobalVector(vdm, &u[0]));
445: PetscCall(DMSwarmSetCellDMActive(sw, "space"));
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: 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[])
450: {
451: f0[0] = 0.;
452: for (PetscInt d = 0; d < dim; ++d) f0[0] += PetscSqr(u_x[uOff_x[0] + d * dim + d]);
453: }
455: static PetscErrorCode MonitorEField(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
456: {
457: AppCtx *user = (AppCtx *)ctx;
458: DM sw;
459: PetscScalar intESq;
460: PetscReal *E, *x, *weight;
461: PetscReal Enorm = 0., lgEnorm, lgEmax, sum = 0., Emax = 0., chargesum = 0.;
462: PetscReal pmoments[4]; /* \int f, \int v f, \int v^2 f */
463: PetscInt *species, dim, Np, gNp;
464: MPI_Comm comm;
466: PetscFunctionBeginUser;
467: if (step < 0 || !user->validE) PetscFunctionReturn(PETSC_SUCCESS);
468: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
469: PetscCall(TSGetDM(ts, &sw));
470: PetscCall(DMGetDimension(sw, &dim));
471: PetscCall(DMSwarmGetLocalSize(sw, &Np));
472: PetscCall(DMSwarmGetSize(sw, &gNp));
473: PetscCall(DMSwarmSortGetAccess(sw));
474: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
475: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
476: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
477: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
479: for (PetscInt p = 0; p < Np; ++p) {
480: for (PetscInt d = 0; d < 1; ++d) {
481: PetscReal temp = PetscAbsReal(E[p * dim + d]);
482: if (temp > Emax) Emax = temp;
483: }
484: Enorm += PetscSqrtReal(E[p * dim] * E[p * dim]);
485: sum += E[p * dim];
486: chargesum += user->charges[0] * weight[p];
487: }
488: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &Emax, 1, MPIU_REAL, MPIU_MAX, comm));
489: lgEnorm = Enorm != 0 ? PetscLog10Real(Enorm) : -16.;
490: lgEmax = Emax != 0 ? PetscLog10Real(Emax) : user->drawlgEmin;
492: PetscDS ds;
493: Vec phi;
495: PetscCall(DMGetNamedGlobalVector(user->dmPot, "phi", &phi));
496: PetscCall(DMGetDS(user->dmPot, &ds));
497: PetscCall(PetscDSSetObjective(ds, 0, &f0_grad_phi2));
498: PetscCall(DMPlexComputeIntegralFEM(user->dmPot, phi, &intESq, user));
499: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "phi", &phi));
501: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
502: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
503: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
504: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
505: PetscCall(PetscDrawLGAddPoint(user->drawlgE, &t, &lgEmax));
506: PetscCall(PetscDrawLGDraw(user->drawlgE));
507: PetscDraw draw;
508: PetscCall(PetscDrawLGGetDraw(user->drawlgE, &draw));
509: PetscCall(PetscDrawSave(draw));
511: PetscCall(DMSwarmComputeMoments(sw, "velocity", "w_q", pmoments));
512: 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));
513: PetscCall(DMViewFromOptions(sw, NULL, "-sw_efield_view"));
514: PetscFunctionReturn(PETSC_SUCCESS);
515: }
517: static PetscErrorCode MonitorMoments(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
518: {
519: AppCtx *user = (AppCtx *)ctx;
520: DM sw;
521: PetscReal pmoments[4], fmoments[4]; /* \int f, \int v f, \int v^2 f */
523: PetscFunctionBeginUser;
524: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS);
525: PetscCall(TSGetDM(ts, &sw));
527: PetscCall(DMSwarmComputeMoments(sw, "velocity", "w_q", pmoments));
528: PetscCall(computeVelocityFEMMoments(sw, fmoments, user));
530: 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]));
531: PetscFunctionReturn(PETSC_SUCCESS);
532: }
534: PetscErrorCode MonitorInitialConditions(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
535: {
536: AppCtx *user = (AppCtx *)ctx;
537: DM sw;
538: PetscDraw drawic_x, drawic_v;
539: PetscReal *weight, *pos, *vel;
540: PetscInt dim, Np;
542: PetscFunctionBegin;
543: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* -1 indicates interpolated solution */
544: if (step == 0) {
545: PetscCall(TSGetDM(ts, &sw));
546: PetscCall(DMGetDimension(sw, &dim));
547: PetscCall(DMSwarmGetLocalSize(sw, &Np));
549: PetscCall(PetscDrawHGReset(user->drawhgic_x));
550: PetscCall(PetscDrawHGGetDraw(user->drawhgic_x, &drawic_x));
551: PetscCall(PetscDrawClear(drawic_x));
552: PetscCall(PetscDrawFlush(drawic_x));
554: PetscCall(PetscDrawHGReset(user->drawhgic_v));
555: PetscCall(PetscDrawHGGetDraw(user->drawhgic_v, &drawic_v));
556: PetscCall(PetscDrawClear(drawic_v));
557: PetscCall(PetscDrawFlush(drawic_v));
559: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&pos));
560: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&vel));
561: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
562: for (PetscInt p = 0; p < Np; ++p) {
563: PetscCall(PetscDrawHGAddWeightedValue(user->drawhgic_x, pos[p * dim], weight[p]));
564: PetscCall(PetscDrawHGAddWeightedValue(user->drawhgic_v, vel[p * dim], weight[p]));
565: }
566: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&pos));
567: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&vel));
568: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
570: PetscCall(PetscDrawHGDraw(user->drawhgic_x));
571: PetscCall(PetscDrawHGSave(user->drawhgic_x));
572: PetscCall(PetscDrawHGDraw(user->drawhgic_v));
573: PetscCall(PetscDrawHGSave(user->drawhgic_v));
574: }
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: // Right now, make the complete velocity histogram
579: PetscErrorCode MonitorVelocity(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
580: {
581: AppCtx *user = (AppCtx *)ctx;
582: DM sw, dm;
583: Vec ks;
584: PetscProbFunc cdf;
585: PetscDraw drawcell_v;
586: PetscScalar *ksa;
587: PetscReal *weight, *vel;
588: PetscInt *pidx;
589: PetscInt dim, Npc, cStart, cEnd, cell = user->velocity_monitor;
591: PetscFunctionBegin;
592: PetscCall(TSGetDM(ts, &sw));
593: PetscCall(DMGetDimension(sw, &dim));
595: PetscCall(DMSwarmGetCellDM(sw, &dm));
596: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
597: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &ks));
598: PetscCall(PetscObjectSetName((PetscObject)ks, "KS Statistic by Cell"));
599: PetscCall(VecSetSizes(ks, cEnd - cStart, PETSC_DETERMINE));
600: PetscCall(VecSetFromOptions(ks));
601: switch (dim) {
602: case 1:
603: //cdf = PetscCDFMaxwellBoltzmann1D;
604: cdf = PetscCDFGaussian1D;
605: break;
606: case 2:
607: cdf = PetscCDFMaxwellBoltzmann2D;
608: break;
609: case 3:
610: cdf = PetscCDFMaxwellBoltzmann3D;
611: break;
612: default:
613: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported", dim);
614: }
616: PetscCall(PetscDrawHGReset(user->drawhgcell_v));
617: PetscCall(PetscDrawHGGetDraw(user->drawhgcell_v, &drawcell_v));
618: PetscCall(PetscDrawClear(drawcell_v));
619: PetscCall(PetscDrawFlush(drawcell_v));
621: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&vel));
622: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
623: PetscCall(DMSwarmSortGetAccess(sw));
624: PetscCall(VecGetArrayWrite(ks, &ksa));
625: for (PetscInt c = cStart; c < cEnd; ++c) {
626: Vec cellv, cellw;
627: PetscScalar *cella, *cellaw;
628: PetscReal totWgt = 0.;
630: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
631: PetscCall(VecCreate(PETSC_COMM_SELF, &cellv));
632: PetscCall(VecSetBlockSize(cellv, dim));
633: PetscCall(VecSetSizes(cellv, Npc * dim, Npc));
634: PetscCall(VecSetFromOptions(cellv));
635: PetscCall(VecCreate(PETSC_COMM_SELF, &cellw));
636: PetscCall(VecSetSizes(cellw, Npc, Npc));
637: PetscCall(VecSetFromOptions(cellw));
638: PetscCall(VecGetArrayWrite(cellv, &cella));
639: PetscCall(VecGetArrayWrite(cellw, &cellaw));
640: for (PetscInt q = 0; q < Npc; ++q) {
641: const PetscInt p = pidx[q];
642: if (c == cell) PetscCall(PetscDrawHGAddWeightedValue(user->drawhgcell_v, vel[p * dim], weight[p]));
643: for (PetscInt d = 0; d < dim; ++d) cella[q * dim + d] = vel[p * dim + d];
644: cellaw[q] = weight[p];
645: totWgt += weight[p];
646: }
647: PetscCall(VecRestoreArrayWrite(cellv, &cella));
648: PetscCall(VecRestoreArrayWrite(cellw, &cellaw));
649: PetscCall(VecScale(cellw, 1. / totWgt));
650: PetscCall(PetscProbComputeKSStatisticWeighted(cellv, cellw, cdf, &ksa[c - cStart]));
651: PetscCall(VecDestroy(&cellv));
652: PetscCall(VecDestroy(&cellw));
653: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
654: }
655: PetscCall(VecRestoreArrayWrite(ks, &ksa));
656: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&vel));
657: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
658: PetscCall(DMSwarmSortRestoreAccess(sw));
660: PetscReal minalpha, maxalpha;
661: PetscInt mincell, maxcell;
663: PetscCall(VecFilter(ks, PETSC_SMALL));
664: PetscCall(VecMin(ks, &mincell, &minalpha));
665: PetscCall(VecMax(ks, &maxcell, &maxalpha));
666: 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));
667: PetscCall(VecViewFromOptions(ks, NULL, "-ks_view"));
668: PetscCall(VecDestroy(&ks));
670: PetscCall(PetscDrawHGDraw(user->drawhgcell_v));
671: PetscCall(PetscDrawHGSave(user->drawhgcell_v));
672: PetscFunctionReturn(PETSC_SUCCESS);
673: }
675: static PetscErrorCode MonitorPositions_2D(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
676: {
677: AppCtx *user = (AppCtx *)ctx;
678: DM dm, sw;
679: PetscScalar *x, *v, *weight;
680: PetscReal lower[3], upper[3], speed;
681: const PetscInt *s;
682: PetscInt dim, cStart, cEnd, c;
684: PetscFunctionBeginUser;
685: if (step > 0 && step % user->ostep == 0) {
686: PetscCall(TSGetDM(ts, &sw));
687: PetscCall(DMSwarmGetCellDM(sw, &dm));
688: PetscCall(DMGetDimension(dm, &dim));
689: PetscCall(DMGetBoundingBox(dm, lower, upper));
690: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
691: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
692: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
693: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
694: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&s));
695: PetscCall(DMSwarmSortGetAccess(sw));
696: PetscCall(PetscDrawSPReset(user->drawspX));
697: PetscCall(PetscDrawSPSetLimits(user->drawspX, lower[0], upper[0], lower[1], upper[1]));
698: PetscCall(PetscDrawSPSetLimits(user->drawspX, lower[0], upper[0], -12, 12));
699: for (c = 0; c < cEnd - cStart; ++c) {
700: PetscInt *pidx, Npc, q;
701: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
702: for (q = 0; q < Npc; ++q) {
703: const PetscInt p = pidx[q];
704: if (s[p] == 0) {
705: speed = 0.;
706: for (PetscInt d = 0; d < dim; ++d) speed += PetscSqr(v[p * dim + d]);
707: speed = PetscSqrtReal(speed);
708: if (dim == 1) {
709: PetscCall(PetscDrawSPAddPointColorized(user->drawspX, &x[p * dim], &v[p * dim], &speed));
710: } else {
711: PetscCall(PetscDrawSPAddPointColorized(user->drawspX, &x[p * dim], &x[p * dim + 1], &speed));
712: }
713: } else if (s[p] == 1) {
714: PetscCall(PetscDrawSPAddPoint(user->drawspX, &x[p * dim], &v[p * dim]));
715: }
716: }
717: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
718: }
719: PetscCall(PetscDrawSPDraw(user->drawspX, PETSC_TRUE));
720: PetscDraw draw;
721: PetscCall(PetscDrawSPGetDraw(user->drawspX, &draw));
722: PetscCall(PetscDrawSave(draw));
723: PetscCall(DMSwarmSortRestoreAccess(sw));
724: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
725: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
726: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
727: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&s));
728: }
729: PetscFunctionReturn(PETSC_SUCCESS);
730: }
732: static PetscErrorCode MonitorPoisson(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
733: {
734: AppCtx *user = (AppCtx *)ctx;
735: DM dm, sw;
737: PetscFunctionBeginUser;
738: if (step > 0 && step % user->ostep == 0) {
739: PetscCall(TSGetDM(ts, &sw));
740: PetscCall(DMSwarmGetCellDM(sw, &dm));
742: if (user->validE) {
743: PetscScalar *x, *E, *weight;
744: PetscReal lower[3], upper[3], xval;
745: PetscDraw draw;
746: PetscInt dim, cStart, cEnd;
748: PetscCall(DMGetDimension(dm, &dim));
749: PetscCall(DMGetBoundingBox(dm, lower, upper));
750: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
752: PetscCall(PetscDrawSPReset(user->drawspE));
753: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
754: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
755: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
757: PetscCall(DMSwarmSortGetAccess(sw));
758: for (PetscInt c = 0; c < cEnd - cStart; ++c) {
759: PetscReal Eavg = 0.0;
760: PetscInt *pidx, Npc;
762: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
763: for (PetscInt q = 0; q < Npc; ++q) {
764: const PetscInt p = pidx[q];
765: Eavg += E[p * dim];
766: }
767: Eavg /= Npc;
768: xval = (c + 0.5) * ((upper[0] - lower[0]) / (cEnd - cStart));
769: PetscCall(PetscDrawSPAddPoint(user->drawspE, &xval, &Eavg));
770: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
771: }
772: PetscCall(PetscDrawSPDraw(user->drawspE, PETSC_TRUE));
773: PetscCall(PetscDrawSPGetDraw(user->drawspE, &draw));
774: PetscCall(PetscDrawSave(draw));
775: PetscCall(DMSwarmSortRestoreAccess(sw));
776: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
777: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
778: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
779: }
781: Vec rho, rhohat, phi;
783: PetscCall(DMGetNamedGlobalVector(user->dmPot, "rho", &rho));
784: PetscCall(DMGetNamedGlobalVector(user->dmPot, "rhohat", &rhohat));
785: PetscCall(VecView(rho, user->viewerRho));
786: PetscCall(VecISCopy(user->fftX, user->fftReal, SCATTER_FORWARD, rho));
787: PetscCall(MatMult(user->fftPot, user->fftX, user->fftY));
788: PetscCall(VecFilter(user->fftY, PETSC_SMALL));
789: PetscCall(VecViewFromOptions(user->fftX, NULL, "-real_view"));
790: PetscCall(VecViewFromOptions(user->fftY, NULL, "-fft_view"));
791: PetscCall(VecISCopy(user->fftY, user->fftReal, SCATTER_REVERSE, rhohat));
792: PetscCall(VecSetValue(rhohat, 0, 0., INSERT_VALUES)); // Remove large DC component
793: PetscCall(VecView(rhohat, user->viewerRhoHat));
794: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "rho", &rho));
795: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "rhohat", &rhohat));
797: PetscCall(DMGetNamedGlobalVector(user->dmPot, "phi", &phi));
798: PetscCall(VecView(phi, user->viewerPhi));
799: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "phi", &phi));
800: }
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
805: {
806: PetscBag bag;
807: Parameter *p;
809: PetscFunctionBeginUser;
810: /* setup PETSc parameter bag */
811: PetscCall(PetscBagGetData(ctx->bag, (void **)&p));
812: PetscCall(PetscBagSetName(ctx->bag, "par", "Vlasov-Poisson Parameters"));
813: bag = ctx->bag;
814: PetscCall(PetscBagRegisterScalar(bag, &p->v0, 1.0, "v0", "Velocity scale, m/s"));
815: PetscCall(PetscBagRegisterScalar(bag, &p->t0, 1.0, "t0", "Time scale, s"));
816: PetscCall(PetscBagRegisterScalar(bag, &p->x0, 1.0, "x0", "Space scale, m"));
817: PetscCall(PetscBagRegisterScalar(bag, &p->v0, 1.0, "phi0", "Potential scale, kg*m^2/A*s^3"));
818: PetscCall(PetscBagRegisterScalar(bag, &p->q0, 1.0, "q0", "Charge Scale, A*s"));
819: PetscCall(PetscBagRegisterScalar(bag, &p->m0, 1.0, "m0", "Mass Scale, kg"));
820: PetscCall(PetscBagRegisterScalar(bag, &p->epsi0, 1.0, "epsi0", "Permittivity of Free Space, kg"));
821: PetscCall(PetscBagRegisterScalar(bag, &p->kb, 1.0, "kb", "Boltzmann Constant, m^2 kg/s^2 K^1"));
823: PetscCall(PetscBagRegisterScalar(bag, &p->sigma, 1.0, "sigma", "Charge per unit area, C/m^3"));
824: PetscCall(PetscBagRegisterScalar(bag, &p->poissonNumber, 1.0, "poissonNumber", "Non-Dimensional Poisson Number"));
825: PetscCall(PetscBagRegisterScalar(bag, &p->vlasovNumber, 1.0, "vlasovNumber", "Non-Dimensional Vlasov Number"));
826: PetscCall(PetscBagSetFromOptions(bag));
827: {
828: PetscViewer viewer;
829: PetscViewerFormat format;
830: PetscBool flg;
832: PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg));
833: if (flg) {
834: PetscCall(PetscViewerPushFormat(viewer, format));
835: PetscCall(PetscBagView(bag, viewer));
836: PetscCall(PetscViewerFlush(viewer));
837: PetscCall(PetscViewerPopFormat(viewer));
838: PetscCall(PetscViewerDestroy(&viewer));
839: }
840: }
841: PetscFunctionReturn(PETSC_SUCCESS);
842: }
844: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
845: {
846: PetscFunctionBeginUser;
847: PetscCall(DMCreate(comm, dm));
848: PetscCall(DMSetType(*dm, DMPLEX));
849: PetscCall(DMSetFromOptions(*dm));
850: PetscCall(PetscObjectSetName((PetscObject)*dm, "space"));
851: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
853: // Cache the mesh geometry
854: DMField coordField;
855: IS cellIS;
856: PetscQuadrature quad;
857: PetscReal *wt, *pt;
858: PetscInt cdim, cStart, cEnd;
860: PetscCall(DMGetCoordinateField(*dm, &coordField));
861: PetscCheck(coordField, comm, PETSC_ERR_USER, "DM must have a coordinate field");
862: PetscCall(DMGetCoordinateDim(*dm, &cdim));
863: PetscCall(DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd));
864: PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cStart, cStart, 1, &cellIS));
865: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
866: PetscCall(PetscMalloc1(1, &wt));
867: PetscCall(PetscMalloc1(2, &pt));
868: wt[0] = 1.;
869: pt[0] = -1.;
870: pt[1] = -1.;
871: PetscCall(PetscQuadratureSetData(quad, cdim, 1, 1, pt, wt));
872: PetscCall(DMFieldCreateFEGeom(coordField, cellIS, quad, PETSC_FEGEOM_BASIC, &user->fegeom));
873: PetscCall(PetscQuadratureDestroy(&quad));
874: PetscCall(ISDestroy(&cellIS));
875: PetscFunctionReturn(PETSC_SUCCESS);
876: }
878: 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[])
879: {
880: f0[0] = -constants[SIGMA];
881: }
883: 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[])
884: {
885: PetscInt d;
886: for (d = 0; d < dim; ++d) f1[d] = u_x[d];
887: }
889: 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[])
890: {
891: PetscInt d;
892: for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
893: }
895: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
896: {
897: *u = 0.0;
898: return PETSC_SUCCESS;
899: }
901: /*
902: / I -grad\ / q \ = /0\
903: \-div 0 / \phi/ \f/
904: */
905: 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[])
906: {
907: for (PetscInt d = 0; d < dim; ++d) f0[d] += u[uOff[0] + d];
908: }
910: 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[])
911: {
912: for (PetscInt d = 0; d < dim; ++d) f1[d * dim + d] = u[uOff[1]];
913: }
915: 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[])
916: {
917: f0[0] += constants[SIGMA];
918: for (PetscInt d = 0; d < dim; ++d) f0[0] += u_x[uOff_x[0] + d * dim + d];
919: }
921: /* Boundary residual. Dirichlet boundary for u means u_bdy=p*n */
922: 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[])
923: {
924: for (PetscInt d = 0; d < dim; ++d) g0[d * dim + d] = 1.0;
925: }
927: 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[])
928: {
929: for (PetscInt d = 0; d < dim; ++d) g2[d * dim + d] = 1.0;
930: }
932: 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[])
933: {
934: for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
935: }
937: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
938: {
939: PetscFE fephi, feq;
940: PetscDS ds;
941: PetscBool simplex;
942: PetscInt dim;
944: PetscFunctionBeginUser;
945: PetscCall(DMGetDimension(dm, &dim));
946: PetscCall(DMPlexIsSimplex(dm, &simplex));
947: if (user->em == EM_MIXED) {
948: DMLabel label;
949: const PetscInt id = 1;
951: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "field_", PETSC_DETERMINE, &feq));
952: PetscCall(PetscObjectSetName((PetscObject)feq, "field"));
953: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "potential_", PETSC_DETERMINE, &fephi));
954: PetscCall(PetscObjectSetName((PetscObject)fephi, "potential"));
955: PetscCall(PetscFECopyQuadrature(feq, fephi));
956: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feq));
957: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fephi));
958: PetscCall(DMCreateDS(dm));
959: PetscCall(PetscFEDestroy(&fephi));
960: PetscCall(PetscFEDestroy(&feq));
962: PetscCall(DMGetLabel(dm, "marker", &label));
963: PetscCall(DMGetDS(dm, &ds));
965: PetscCall(PetscDSSetResidual(ds, 0, f0_q, f1_q));
966: PetscCall(PetscDSSetResidual(ds, 1, f0_phi_backgroundCharge, NULL));
967: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_qq, NULL, NULL, NULL));
968: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_qphi, NULL));
969: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_phiq, NULL, NULL));
971: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))zero, NULL, NULL, NULL));
973: } else {
974: MatNullSpace nullsp;
975: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &fephi));
976: PetscCall(PetscObjectSetName((PetscObject)fephi, "potential"));
977: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fephi));
978: PetscCall(DMCreateDS(dm));
979: PetscCall(DMGetDS(dm, &ds));
980: PetscCall(PetscDSSetResidual(ds, 0, ion_f0, laplacian_f1));
981: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, laplacian_g3));
982: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_TRUE, 0, NULL, &nullsp));
983: PetscCall(PetscObjectCompose((PetscObject)fephi, "nullspace", (PetscObject)nullsp));
984: PetscCall(MatNullSpaceDestroy(&nullsp));
985: PetscCall(PetscFEDestroy(&fephi));
986: }
987: PetscFunctionReturn(PETSC_SUCCESS);
988: }
990: static PetscErrorCode CreatePoisson(DM dm, AppCtx *user)
991: {
992: SNES snes;
993: Mat J;
994: MatNullSpace nullSpace;
996: PetscFunctionBeginUser;
997: PetscCall(CreateFEM(dm, user));
998: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
999: PetscCall(SNESSetOptionsPrefix(snes, "em_"));
1000: PetscCall(SNESSetDM(snes, dm));
1001: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, user));
1002: PetscCall(SNESSetFromOptions(snes));
1004: PetscCall(DMCreateMatrix(dm, &J));
1005: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_TRUE, 0, NULL, &nullSpace));
1006: PetscCall(MatSetNullSpace(J, nullSpace));
1007: PetscCall(MatNullSpaceDestroy(&nullSpace));
1008: PetscCall(SNESSetJacobian(snes, J, J, NULL, NULL));
1009: PetscCall(MatDestroy(&J));
1010: if (user->em == EM_MIXED) {
1011: const PetscInt potential = 1;
1013: PetscCall(DMCreateSubDM(dm, 1, &potential, &user->isPot, &user->dmPot));
1014: } else {
1015: user->dmPot = dm;
1016: PetscCall(PetscObjectReference((PetscObject)user->dmPot));
1017: }
1018: PetscCall(DMCreateMassMatrix(user->dmPot, user->dmPot, &user->M));
1019: PetscCall(DMPlexCreateClosureIndex(dm, NULL));
1020: user->snes = snes;
1021: PetscFunctionReturn(PETSC_SUCCESS);
1022: }
1024: PetscErrorCode PetscPDFPertubedConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
1025: {
1026: p[0] = (1 + 0.01 * PetscCosReal(0.5 * x[0])) / (2 * PETSC_PI);
1027: p[1] = (1 + 0.01 * PetscCosReal(0.5 * x[1])) / (2 * PETSC_PI);
1028: return PETSC_SUCCESS;
1029: }
1030: PetscErrorCode PetscPDFPertubedConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
1031: {
1032: p[0] = (1. + 0.01 * PetscCosReal(0.5 * x[0])) / (2 * PETSC_PI);
1033: return PETSC_SUCCESS;
1034: }
1036: PetscErrorCode PetscPDFCosine1D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
1037: {
1038: const PetscReal alpha = scale ? scale[0] : 0.0;
1039: const PetscReal k = scale ? scale[1] : 1.;
1040: p[0] = (1 + alpha * PetscCosReal(k * x[0]));
1041: return PETSC_SUCCESS;
1042: }
1044: PetscErrorCode PetscPDFCosine2D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
1045: {
1046: const PetscReal alpha = scale ? scale[0] : 0.;
1047: const PetscReal k = scale ? scale[0] : 1.;
1048: p[0] = (1 + alpha * PetscCosReal(k * (x[0] + x[1])));
1049: return PETSC_SUCCESS;
1050: }
1052: static PetscErrorCode CreateVelocityDM(DM sw, DM *vdm)
1053: {
1054: PetscFE fe;
1055: DMPolytopeType ct;
1056: PetscInt dim, cStart;
1057: const char *prefix = "v";
1059: PetscFunctionBegin;
1060: PetscCall(DMCreate(PETSC_COMM_SELF, vdm));
1061: PetscCall(DMSetType(*vdm, DMPLEX));
1062: PetscCall(DMPlexSetOptionsPrefix(*vdm, prefix));
1063: PetscCall(DMSetFromOptions(*vdm));
1064: PetscCall(PetscObjectSetName((PetscObject)*vdm, "velocity"));
1065: PetscCall(DMViewFromOptions(*vdm, NULL, "-dm_view"));
1067: PetscCall(DMGetDimension(*vdm, &dim));
1068: PetscCall(DMPlexGetHeightStratum(*vdm, 0, &cStart, NULL));
1069: PetscCall(DMPlexGetCellType(*vdm, cStart, &ct));
1070: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, prefix, PETSC_DETERMINE, &fe));
1071: PetscCall(PetscObjectSetName((PetscObject)fe, "distribution"));
1072: PetscCall(DMSetField(*vdm, 0, NULL, (PetscObject)fe));
1073: PetscCall(DMCreateDS(*vdm));
1074: PetscCall(PetscFEDestroy(&fe));
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: /*
1079: InitializeParticles_Centroid - Initialize a regular grid of particles.
1081: Input Parameters:
1082: + sw - The `DMSWARM`
1083: - force1D - Treat the spatial domain as 1D
1085: Notes:
1086: This functions sets the species, cellid, spatial coordinate, and velocity fields for all particles.
1088: It places one particle in the centroid of each cell in the implicit tensor product of the spatial
1089: and velocity meshes.
1090: */
1091: static PetscErrorCode InitializeParticles_Centroid(DM sw)
1092: {
1093: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1094: DMSwarmCellDM celldm;
1095: DM xdm, vdm;
1096: PetscReal vmin[3], vmax[3];
1097: PetscReal *x, *v;
1098: PetscInt *species, *cellid;
1099: PetscInt dim, xcStart, xcEnd, vcStart, vcEnd, Ns, Np, Npc, debug;
1100: PetscBool flg;
1101: MPI_Comm comm;
1102: const char *cellidname;
1104: PetscFunctionBegin;
1105: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1107: PetscOptionsBegin(comm, "", "DMSwarm Options", "DMSWARM");
1108: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1109: PetscCall(PetscOptionsInt("-dm_swarm_num_species", "The number of species", "DMSwarmSetNumSpecies", Ns, &Ns, &flg));
1110: if (flg) PetscCall(DMSwarmSetNumSpecies(sw, Ns));
1111: PetscCall(PetscOptionsBoundedInt("-dm_swarm_print_coords", "Debug output level for particle coordinate computations", "InitializeParticles", 0, &swarm->printCoords, NULL, 0));
1112: PetscCall(PetscOptionsBoundedInt("-dm_swarm_print_weights", "Debug output level for particle weight computations", "InitializeWeights", 0, &swarm->printWeights, NULL, 0));
1113: PetscOptionsEnd();
1114: debug = swarm->printCoords;
1116: PetscCall(DMGetDimension(sw, &dim));
1117: PetscCall(DMSwarmGetCellDM(sw, &xdm));
1118: PetscCall(DMPlexGetHeightStratum(xdm, 0, &xcStart, &xcEnd));
1120: PetscCall(DMSwarmGetCellDMByName(sw, "velocity", &celldm));
1121: PetscCall(DMSwarmCellDMGetDM(celldm, &vdm));
1122: PetscCall(DMPlexGetHeightStratum(vdm, 0, &vcStart, &vcEnd));
1124: // One particle per centroid on the tensor product grid
1125: Npc = (vcEnd - vcStart) * Ns;
1126: Np = (xcEnd - xcStart) * Npc;
1127: PetscCall(DMSwarmSetLocalSizes(sw, Np, 0));
1128: if (debug) {
1129: PetscInt gNp, gNc, Nc = xcEnd - xcStart;
1130: PetscCallMPI(MPIU_Allreduce(&Np, &gNp, 1, MPIU_INT, MPIU_SUM, comm));
1131: PetscCall(PetscPrintf(comm, "Global Np = %" PetscInt_FMT "\n", gNp));
1132: PetscCallMPI(MPIU_Allreduce(&Nc, &gNc, 1, MPIU_INT, MPIU_SUM, comm));
1133: PetscCall(PetscPrintf(comm, "Global X-cells = %" PetscInt_FMT "\n", gNc));
1134: PetscCall(PetscPrintf(comm, "Global V-cells = %" PetscInt_FMT "\n", vcEnd - vcStart));
1135: }
1137: // Set species and cellid
1138: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1139: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellidname));
1140: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1141: PetscCall(DMSwarmGetField(sw, cellidname, NULL, NULL, (void **)&cellid));
1142: for (PetscInt c = 0, p = 0; c < xcEnd - xcStart; ++c) {
1143: for (PetscInt s = 0; s < Ns; ++s) {
1144: for (PetscInt q = 0; q < Npc / Ns; ++q, ++p) {
1145: species[p] = s;
1146: cellid[p] = c;
1147: }
1148: }
1149: }
1150: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1151: PetscCall(DMSwarmRestoreField(sw, cellidname, NULL, NULL, (void **)&cellid));
1153: // Set particle coordinates
1154: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
1155: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
1156: PetscCall(DMSwarmSortGetAccess(sw));
1157: PetscCall(DMGetBoundingBox(vdm, vmin, vmax));
1158: PetscCall(DMGetCoordinatesLocalSetUp(xdm));
1159: for (PetscInt c = 0; c < xcEnd - xcStart; ++c) {
1160: const PetscInt cell = c + xcStart;
1161: PetscInt *pidx, Npc;
1162: PetscReal centroid[3], volume;
1164: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
1165: PetscCall(DMPlexComputeCellGeometryFVM(xdm, cell, &volume, centroid, NULL));
1166: for (PetscInt s = 0; s < Ns; ++s) {
1167: for (PetscInt q = 0; q < Npc / Ns; ++q) {
1168: const PetscInt p = pidx[q * Ns + s];
1170: for (PetscInt d = 0; d < dim; ++d) {
1171: x[p * dim + d] = centroid[d];
1172: v[p * dim + d] = vmin[0] + (q + 0.5) * ((vmax[0] - vmin[0]) / (Npc / Ns));
1173: }
1174: if (debug > 1) {
1175: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Particle %4" PetscInt_FMT " ", p));
1176: PetscCall(PetscPrintf(PETSC_COMM_SELF, " x: ("));
1177: for (PetscInt d = 0; d < dim; ++d) {
1178: if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
1179: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", x[p * dim + d]));
1180: }
1181: PetscCall(PetscPrintf(PETSC_COMM_SELF, ") v:("));
1182: for (PetscInt d = 0; d < dim; ++d) {
1183: if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
1184: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", v[p * dim + d]));
1185: }
1186: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
1187: }
1188: }
1189: }
1190: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
1191: }
1192: PetscCall(DMSwarmSortRestoreAccess(sw));
1193: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&x));
1194: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: /*
1199: InitializeWeights - Compute weight for each local particle
1201: Input Parameters:
1202: + sw - The `DMSwarm`
1203: . totalWeight - The sum of all particle weights
1204: . func - The PDF for the particle spatial distribution
1205: - param - The PDF parameters
1207: Notes:
1208: The PDF for velocity is assumed to be a Gaussian
1210: The particle weights are returned in the `w_q` field of `sw`.
1211: */
1212: static PetscErrorCode InitializeWeights(DM sw, PetscReal totalWeight, PetscProbFunc func, const PetscReal param[])
1213: {
1214: DM xdm, vdm;
1215: DMSwarmCellDM celldm;
1216: PetscScalar *weight;
1217: PetscQuadrature xquad;
1218: const PetscReal *xq, *xwq;
1219: const PetscInt order = 5;
1220: PetscReal xi0[3];
1221: PetscReal xwtot = 0., pwtot = 0.;
1222: PetscInt xNq;
1223: PetscInt dim, Ns, xcStart, xcEnd, vcStart, vcEnd, debug = ((DM_Swarm *)sw->data)->printWeights;
1224: MPI_Comm comm;
1225: PetscMPIInt rank;
1227: PetscFunctionBegin;
1228: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1229: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1230: PetscCall(DMGetDimension(sw, &dim));
1231: PetscCall(DMSwarmGetCellDM(sw, &xdm));
1232: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1233: PetscCall(DMPlexGetHeightStratum(xdm, 0, &xcStart, &xcEnd));
1234: PetscCall(DMSwarmGetCellDMByName(sw, "velocity", &celldm));
1235: PetscCall(DMSwarmCellDMGetDM(celldm, &vdm));
1236: PetscCall(DMPlexGetHeightStratum(vdm, 0, &vcStart, &vcEnd));
1238: // Setup Quadrature for spatial and velocity weight calculations
1239: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, order, -1.0, 1.0, &xquad));
1240: PetscCall(PetscQuadratureGetData(xquad, NULL, NULL, &xNq, &xq, &xwq));
1241: for (PetscInt d = 0; d < dim; ++d) xi0[d] = -1.0;
1243: // Integrate the density function to get the weights of particles in each cell
1244: PetscCall(DMGetCoordinatesLocalSetUp(vdm));
1245: PetscCall(DMSwarmSortGetAccess(sw));
1246: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
1247: for (PetscInt c = xcStart; c < xcEnd; ++c) {
1248: PetscReal xv0[3], xJ[9], xinvJ[9], xdetJ, xqr[3], xden, xw = 0.;
1249: PetscInt *pidx, Npc;
1250: PetscInt xNc;
1251: const PetscScalar *xarray;
1252: PetscScalar *xcoords = NULL;
1253: PetscBool xisDG;
1255: PetscCall(DMPlexGetCellCoordinates(xdm, c, &xisDG, &xNc, &xarray, &xcoords));
1256: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
1257: 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);
1258: PetscCall(DMPlexComputeCellGeometryFEM(xdm, c, NULL, xv0, xJ, xinvJ, &xdetJ));
1259: for (PetscInt q = 0; q < xNq; ++q) {
1260: // Transform quadrature points from ref space to real space
1261: CoordinatesRefToReal(dim, dim, xi0, xv0, xJ, &xq[q * dim], xqr);
1262: // Get probability density at quad point
1263: // No need to scale xqr since PDF will be periodic
1264: PetscCall((*func)(xqr, param, &xden));
1265: xw += xden * (xwq[q] * xdetJ);
1266: }
1267: xwtot += xw;
1268: if (debug) {
1269: IS globalOrdering;
1270: const PetscInt *ordering;
1272: PetscCall(DMPlexGetCellNumbering(xdm, &globalOrdering));
1273: PetscCall(ISGetIndices(globalOrdering, &ordering));
1274: 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));
1275: PetscCall(ISRestoreIndices(globalOrdering, &ordering));
1276: }
1277: // Set weights to be Gaussian in velocity cells
1278: for (PetscInt vc = vcStart; vc < vcEnd; ++vc) {
1279: const PetscInt p = pidx[vc * Ns + 0];
1280: PetscReal vw = 0.;
1281: PetscInt vNc;
1282: const PetscScalar *varray;
1283: PetscScalar *vcoords = NULL;
1284: PetscBool visDG;
1286: PetscCall(DMPlexGetCellCoordinates(vdm, vc, &visDG, &vNc, &varray, &vcoords));
1287: // TODO: Fix 2 stream Ask Joe
1288: // Two stream function from 1/2pi v^2 e^(-v^2/2)
1289: // 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.)));
1290: vw = 0.5 * (PetscErfReal(vcoords[1] / PetscSqrtReal(2.)) - PetscErfReal(vcoords[0] / PetscSqrtReal(2.)));
1292: weight[p] = totalWeight * vw * xw;
1293: pwtot += weight[p];
1294: PetscCheck(weight[p] <= 10., PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Particle %" PetscInt_FMT " weight exceeded 1: %g, %g, %g", p, xw, vw, totalWeight);
1295: PetscCall(DMPlexRestoreCellCoordinates(vdm, vc, &visDG, &vNc, &varray, &vcoords));
1296: if (debug > 1) PetscCall(PetscPrintf(comm, "particle %" PetscInt_FMT ": %g, vw: %g xw: %g\n", p, weight[p], vw, xw));
1297: }
1298: PetscCall(DMPlexRestoreCellCoordinates(xdm, c, &xisDG, &xNc, &xarray, &xcoords));
1299: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
1300: }
1301: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
1302: PetscCall(DMSwarmSortRestoreAccess(sw));
1303: PetscCall(PetscQuadratureDestroy(&xquad));
1305: if (debug) {
1306: PetscReal wtot[2] = {pwtot, xwtot}, gwtot[2];
1308: PetscCall(PetscSynchronizedFlush(comm, NULL));
1309: PetscCallMPI(MPIU_Allreduce(wtot, gwtot, 2, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
1310: PetscCall(PetscPrintf(comm, "particle weight sum = %1.10f cell weight sum = %1.10f\n", (double)gwtot[0], (double)gwtot[1]));
1311: }
1312: PetscFunctionReturn(PETSC_SUCCESS);
1313: }
1315: static PetscErrorCode InitializeParticles_PerturbedWeights(DM sw, AppCtx *user)
1316: {
1317: PetscReal scale[2] = {user->cosine_coefficients[0], user->cosine_coefficients[1]};
1318: PetscInt dim;
1320: PetscFunctionBegin;
1321: PetscCall(DMGetDimension(sw, &dim));
1322: PetscCall(InitializeParticles_Centroid(sw));
1323: PetscCall(InitializeWeights(sw, user->totalWeight, dim == 1 ? PetscPDFCosine1D : PetscPDFCosine2D, scale));
1324: PetscFunctionReturn(PETSC_SUCCESS);
1325: }
1327: static PetscErrorCode InitializeConstants(DM sw, AppCtx *user)
1328: {
1329: DM dm;
1330: PetscInt *species;
1331: PetscReal *weight, totalCharge = 0., totalWeight = 0., gmin[3], gmax[3], global_charge, global_weight;
1332: PetscInt Np, dim;
1334: PetscFunctionBegin;
1335: PetscCall(DMSwarmGetCellDM(sw, &dm));
1336: PetscCall(DMGetDimension(sw, &dim));
1337: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1338: PetscCall(DMGetBoundingBox(dm, gmin, gmax));
1339: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
1340: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1341: for (PetscInt p = 0; p < Np; ++p) {
1342: totalWeight += weight[p];
1343: totalCharge += user->charges[species[p]] * weight[p];
1344: }
1345: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
1346: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1347: {
1348: Parameter *param;
1349: PetscReal Area;
1351: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1352: switch (dim) {
1353: case 1:
1354: Area = (gmax[0] - gmin[0]);
1355: break;
1356: case 2:
1357: Area = (gmax[0] - gmin[0]) * (gmax[1] - gmin[1]);
1358: break;
1359: case 3:
1360: Area = (gmax[0] - gmin[0]) * (gmax[1] - gmin[1]) * (gmax[2] - gmin[2]);
1361: break;
1362: default:
1363: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", dim);
1364: }
1365: PetscCallMPI(MPIU_Allreduce(&totalWeight, &global_weight, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
1366: PetscCallMPI(MPIU_Allreduce(&totalCharge, &global_charge, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD));
1367: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "dim = %" PetscInt_FMT "\ttotalWeight = %f, user->charges[species[0]] = %f\ttotalCharge = %f, Total Area = %f\n", dim, (double)global_weight, (double)user->charges[0], (double)global_charge, (double)Area));
1368: param->sigma = PetscAbsReal(global_charge / (Area));
1370: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sigma: %g\n", (double)param->sigma));
1371: 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,
1372: (double)param->vlasovNumber));
1373: }
1374: /* Setup Constants */
1375: {
1376: PetscDS ds;
1377: Parameter *param;
1378: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1379: PetscScalar constants[NUM_CONSTANTS];
1380: constants[SIGMA] = param->sigma;
1381: constants[V0] = param->v0;
1382: constants[T0] = param->t0;
1383: constants[X0] = param->x0;
1384: constants[M0] = param->m0;
1385: constants[Q0] = param->q0;
1386: constants[PHI0] = param->phi0;
1387: constants[POISSON] = param->poissonNumber;
1388: constants[VLASOV] = param->vlasovNumber;
1389: PetscCall(DMGetDS(dm, &ds));
1390: PetscCall(PetscDSSetConstants(ds, NUM_CONSTANTS, constants));
1391: }
1392: PetscFunctionReturn(PETSC_SUCCESS);
1393: }
1395: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
1396: {
1397: DMSwarmCellDM celldm;
1398: DM vdm;
1399: PetscReal v0[2] = {1., 0.};
1400: PetscInt dim;
1402: PetscFunctionBeginUser;
1403: PetscCall(DMGetDimension(dm, &dim));
1404: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
1405: PetscCall(DMSetType(*sw, DMSWARM));
1406: PetscCall(DMSetDimension(*sw, dim));
1407: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
1408: PetscCall(DMSetApplicationContext(*sw, user));
1410: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
1411: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
1412: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT));
1413: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "E_field", dim, PETSC_REAL));
1415: const char *fieldnames[2] = {DMSwarmPICField_coor, "velocity"};
1417: PetscCall(DMSwarmCellDMCreate(dm, 2, fieldnames, 1, fieldnames, &celldm));
1418: PetscCall(DMSwarmAddCellDM(*sw, celldm));
1419: PetscCall(DMSwarmCellDMDestroy(&celldm));
1421: const char *vfieldnames[1] = {"w_q"};
1423: PetscCall(CreateVelocityDM(*sw, &vdm));
1424: PetscCall(DMSwarmCellDMCreate(vdm, 1, vfieldnames, 1, &fieldnames[1], &celldm));
1425: PetscCall(DMSwarmAddCellDM(*sw, celldm));
1426: PetscCall(DMSwarmCellDMDestroy(&celldm));
1427: PetscCall(DMDestroy(&vdm));
1429: DM mdm;
1431: PetscCall(DMClone(dm, &mdm));
1432: PetscCall(PetscObjectSetName((PetscObject)mdm, "moments"));
1433: PetscCall(DMCopyDisc(dm, mdm));
1434: PetscCall(DMSwarmCellDMCreate(mdm, 1, vfieldnames, 1, fieldnames, &celldm));
1435: PetscCall(DMDestroy(&mdm));
1436: PetscCall(DMSwarmAddCellDM(*sw, celldm));
1437: PetscCall(DMSwarmCellDMDestroy(&celldm));
1439: PetscCall(DMSetFromOptions(*sw));
1440: PetscCall(DMSetUp(*sw));
1442: PetscCall(DMSwarmSetCellDMActive(*sw, "space"));
1443: user->swarm = *sw;
1444: // TODO: This is redundant init since it is done in InitializeSolveAndSwarm, however DMSetUp() requires the local size be set
1445: if (user->perturbed_weights) {
1446: PetscCall(InitializeParticles_PerturbedWeights(*sw, user));
1447: } else {
1448: PetscCall(DMSwarmComputeLocalSizeFromOptions(*sw));
1449: PetscCall(DMSwarmInitializeCoordinates(*sw));
1450: PetscCall(DMSwarmInitializeVelocitiesFromOptions(*sw, v0));
1451: }
1452: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
1453: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
1454: PetscFunctionReturn(PETSC_SUCCESS);
1455: }
1457: static PetscErrorCode ComputeFieldAtParticles_Coulomb(SNES snes, DM sw, PetscReal E[])
1458: {
1459: AppCtx *user;
1460: PetscReal *coords;
1461: PetscInt *species, dim, Np, Ns;
1462: PetscMPIInt size;
1464: PetscFunctionBegin;
1465: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)snes), &size));
1466: PetscCheck(size == 1, PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Coulomb code only works in serial");
1467: PetscCall(DMGetDimension(sw, &dim));
1468: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1469: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1470: PetscCall(DMGetApplicationContext(sw, (void *)&user));
1472: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1473: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1474: for (PetscInt p = 0; p < Np; ++p) {
1475: PetscReal *pcoord = &coords[p * dim];
1476: PetscReal pE[3] = {0., 0., 0.};
1478: /* Calculate field at particle p due to particle q */
1479: for (PetscInt q = 0; q < Np; ++q) {
1480: PetscReal *qcoord = &coords[q * dim];
1481: PetscReal rpq[3], r, r3, q_q;
1483: if (p == q) continue;
1484: q_q = user->charges[species[q]] * 1.;
1485: for (PetscInt d = 0; d < dim; ++d) rpq[d] = pcoord[d] - qcoord[d];
1486: r = DMPlex_NormD_Internal(dim, rpq);
1487: if (r < PETSC_SQRT_MACHINE_EPSILON) continue;
1488: r3 = PetscPowRealInt(r, 3);
1489: for (PetscInt d = 0; d < dim; ++d) pE[d] += q_q * rpq[d] / r3;
1490: }
1491: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] = pE[d];
1492: }
1493: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1494: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1495: PetscFunctionReturn(PETSC_SUCCESS);
1496: }
1498: static PetscErrorCode ComputeFieldAtParticles_Primal(SNES snes, DM sw, Mat M_p, PetscReal E[])
1499: {
1500: DM dm;
1501: AppCtx *user;
1502: PetscDS ds;
1503: PetscFE fe;
1504: KSP ksp;
1505: Vec rhoRhs; // Weak charge density, \int phi_i rho
1506: Vec rho; // Charge density, M^{-1} rhoRhs
1507: Vec phi, locPhi; // Potential
1508: Vec f; // Particle weights
1509: PetscReal *coords;
1510: PetscInt dim, cStart, cEnd, Np;
1512: PetscFunctionBegin;
1513: PetscCall(DMGetApplicationContext(sw, (void *)&user));
1514: PetscCall(PetscLogEventBegin(user->ESolveEvent, snes, sw, 0, 0));
1515: PetscCall(DMGetDimension(sw, &dim));
1516: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1518: PetscCall(SNESGetDM(snes, &dm));
1519: PetscCall(DMGetGlobalVector(dm, &rhoRhs));
1520: PetscCall(PetscObjectSetName((PetscObject)rhoRhs, "Weak charge density"));
1521: PetscCall(DMGetNamedGlobalVector(user->dmPot, "rho", &rho));
1522: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
1523: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
1525: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
1526: PetscCall(MatViewFromOptions(user->M, NULL, "-m_view"));
1527: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
1529: PetscCall(MatMultTranspose(M_p, f, rhoRhs));
1530: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
1532: // Low-pass filter rhoRhs
1533: PetscInt window = 0;
1534: PetscCall(PetscOptionsGetInt(NULL, NULL, "-rho_average", &window, NULL));
1535: if (window) {
1536: PetscScalar *a;
1537: PetscInt n;
1538: PetscReal width = 2. * window + 1.;
1540: // This will only work for P_1
1541: // This works because integration against a test function is linear
1542: // This only works in serial since I need the periodic values (maybe use FFT)
1543: // Do a real integral against weight function for higher order
1544: PetscCall(VecGetLocalSize(rhoRhs, &n));
1545: PetscCall(VecGetArrayWrite(rhoRhs, &a));
1546: for (PetscInt i = 0; i < n; ++i) {
1547: PetscScalar avg = a[i];
1548: for (PetscInt j = 1; j <= window; ++j) avg += a[(i - j + n) % n] + a[(i + j) % n];
1549: a[i] = avg / width;
1550: //a[i] = (a[(i - 1 + n) % n] + a[i] + a[(i + 1) % n]) / 3.;
1551: }
1552: PetscCall(VecRestoreArrayWrite(rhoRhs, &a));
1553: }
1555: PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &ksp));
1556: PetscCall(KSPSetOptionsPrefix(ksp, "em_proj_"));
1557: PetscCall(KSPSetOperators(ksp, user->M, user->M));
1558: PetscCall(KSPSetFromOptions(ksp));
1559: PetscCall(KSPSolve(ksp, rhoRhs, rho));
1561: PetscCall(VecScale(rhoRhs, -1.0));
1563: PetscCall(VecViewFromOptions(rhoRhs, NULL, "-rho_view"));
1564: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "rho", &rho));
1565: PetscCall(KSPDestroy(&ksp));
1567: PetscCall(DMGetNamedGlobalVector(user->dmPot, "phi", &phi));
1568: PetscCall(VecSet(phi, 0.0));
1569: PetscCall(SNESSolve(snes, rhoRhs, phi));
1570: PetscCall(DMRestoreGlobalVector(dm, &rhoRhs));
1571: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
1573: PetscCall(DMGetLocalVector(dm, &locPhi));
1574: PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
1575: PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
1576: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "phi", &phi));
1577: PetscCall(PetscLogEventEnd(user->ESolveEvent, snes, sw, 0, 0));
1579: PetscCall(DMGetDS(dm, &ds));
1580: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
1581: PetscCall(DMSwarmSortGetAccess(sw));
1582: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1583: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1585: PetscCall(PetscLogEventBegin(user->ETabEvent, snes, sw, 0, 0));
1586: PetscTabulation tab;
1587: PetscReal *pcoord, *refcoord;
1588: PetscFEGeom *chunkgeom = NULL;
1589: PetscInt maxNcp = 0;
1591: for (PetscInt c = cStart; c < cEnd; ++c) {
1592: PetscInt Ncp;
1594: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, c, &Ncp));
1595: maxNcp = PetscMax(maxNcp, Ncp);
1596: }
1597: PetscCall(DMGetWorkArray(dm, maxNcp * dim, MPIU_REAL, &refcoord));
1598: PetscCall(DMGetWorkArray(dm, maxNcp * dim, MPIU_REAL, &pcoord));
1599: PetscCall(PetscFECreateTabulation(fe, 1, maxNcp, refcoord, 1, &tab));
1600: for (PetscInt c = cStart; c < cEnd; ++c) {
1601: PetscScalar *clPhi = NULL;
1602: PetscInt *points;
1603: PetscInt Ncp;
1605: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
1606: for (PetscInt cp = 0; cp < Ncp; ++cp)
1607: for (PetscInt d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
1608: {
1609: PetscCall(PetscFEGeomGetChunk(user->fegeom, c - cStart, c - cStart + 1, &chunkgeom));
1610: for (PetscInt i = 0; i < Ncp; ++i) {
1611: const PetscReal x0[3] = {-1., -1., -1.};
1612: CoordinatesRealToRef(dim, dim, x0, chunkgeom->v, chunkgeom->invJ, &pcoord[dim * i], &refcoord[dim * i]);
1613: }
1614: }
1615: PetscCall(PetscFEComputeTabulation(fe, Ncp, refcoord, 1, tab));
1616: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
1617: for (PetscInt cp = 0; cp < Ncp; ++cp) {
1618: const PetscReal *basisDer = tab->T[1];
1619: const PetscInt p = points[cp];
1621: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] = 0.;
1622: PetscCall(PetscFEFreeInterpolateGradient_Static(fe, basisDer, clPhi, dim, chunkgeom->invJ, NULL, cp, &E[p * dim]));
1623: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] *= -1.0;
1624: }
1625: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
1626: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
1627: }
1628: PetscCall(DMRestoreWorkArray(dm, maxNcp * dim, MPIU_REAL, &pcoord));
1629: PetscCall(DMRestoreWorkArray(dm, maxNcp * dim, MPIU_REAL, &refcoord));
1630: PetscCall(PetscTabulationDestroy(&tab));
1631: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1632: PetscCall(DMSwarmSortRestoreAccess(sw));
1633: PetscCall(DMRestoreLocalVector(dm, &locPhi));
1634: PetscCall(PetscFEGeomRestoreChunk(user->fegeom, 0, 1, &chunkgeom));
1635: PetscCall(PetscLogEventEnd(user->ETabEvent, snes, sw, 0, 0));
1636: PetscFunctionReturn(PETSC_SUCCESS);
1637: }
1639: static PetscErrorCode ComputeFieldAtParticles_Mixed(SNES snes, DM sw, Mat M_p, PetscReal E[])
1640: {
1641: DM dm;
1642: AppCtx *user;
1643: PetscDS ds;
1644: PetscFE fe;
1645: KSP ksp;
1646: Vec rhoRhs, rhoRhsFull; // Weak charge density, \int phi_i rho, and embedding in mixed problem
1647: Vec rho; // Charge density, M^{-1} rhoRhs
1648: Vec phi, locPhi, phiFull; // Potential and embedding in mixed problem
1649: Vec f; // Particle weights
1650: PetscReal *coords;
1651: PetscInt dim, cStart, cEnd, Np;
1653: PetscFunctionBegin;
1654: PetscCall(DMGetApplicationContext(sw, &user));
1655: PetscCall(PetscLogEventBegin(user->ESolveEvent, snes, sw, 0, 0));
1656: PetscCall(DMGetDimension(sw, &dim));
1657: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1659: PetscCall(SNESGetDM(snes, &dm));
1660: PetscCall(DMGetGlobalVector(user->dmPot, &rhoRhs));
1661: PetscCall(PetscObjectSetName((PetscObject)rhoRhs, "Weak charge density"));
1662: PetscCall(DMGetGlobalVector(dm, &rhoRhsFull));
1663: PetscCall(PetscObjectSetName((PetscObject)rhoRhsFull, "Weak charge density"));
1664: PetscCall(DMGetNamedGlobalVector(user->dmPot, "rho", &rho));
1665: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
1666: PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
1668: PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
1669: PetscCall(MatViewFromOptions(user->M, NULL, "-m_view"));
1670: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
1672: PetscCall(MatMultTranspose(M_p, f, rhoRhs));
1673: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
1675: PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &ksp));
1676: PetscCall(KSPSetOptionsPrefix(ksp, "em_proj"));
1677: PetscCall(KSPSetOperators(ksp, user->M, user->M));
1678: PetscCall(KSPSetFromOptions(ksp));
1679: PetscCall(KSPSolve(ksp, rhoRhs, rho));
1681: PetscCall(VecISCopy(rhoRhsFull, user->isPot, SCATTER_FORWARD, rhoRhs));
1682: //PetscCall(VecScale(rhoRhsFull, -1.0));
1684: PetscCall(VecViewFromOptions(rhoRhs, NULL, "-rho_view"));
1685: PetscCall(VecViewFromOptions(rhoRhsFull, NULL, "-rho_full_view"));
1686: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "rho", &rho));
1687: PetscCall(DMRestoreGlobalVector(user->dmPot, &rhoRhs));
1688: PetscCall(KSPDestroy(&ksp));
1690: PetscCall(DMGetGlobalVector(dm, &phiFull));
1691: PetscCall(DMGetNamedGlobalVector(user->dmPot, "phi", &phi));
1692: PetscCall(VecSet(phiFull, 0.0));
1693: PetscCall(SNESSolve(snes, rhoRhsFull, phiFull));
1694: PetscCall(DMRestoreGlobalVector(dm, &rhoRhsFull));
1695: PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
1697: PetscCall(VecISCopy(phiFull, user->isPot, SCATTER_REVERSE, phi));
1698: PetscCall(DMRestoreNamedGlobalVector(user->dmPot, "phi", &phi));
1700: PetscCall(DMGetLocalVector(dm, &locPhi));
1701: PetscCall(DMGlobalToLocalBegin(dm, phiFull, INSERT_VALUES, locPhi));
1702: PetscCall(DMGlobalToLocalEnd(dm, phiFull, INSERT_VALUES, locPhi));
1703: PetscCall(DMRestoreGlobalVector(dm, &phiFull));
1704: PetscCall(PetscLogEventEnd(user->ESolveEvent, snes, sw, 0, 0));
1706: PetscCall(DMGetDS(dm, &ds));
1707: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
1708: PetscCall(DMSwarmSortGetAccess(sw));
1709: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1710: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1712: PetscCall(PetscLogEventBegin(user->ETabEvent, snes, sw, 0, 0));
1713: PetscTabulation tab;
1714: PetscReal *pcoord, *refcoord;
1715: PetscFEGeom *chunkgeom = NULL;
1716: PetscInt maxNcp = 0;
1718: for (PetscInt c = cStart; c < cEnd; ++c) {
1719: PetscInt Ncp;
1721: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, c, &Ncp));
1722: maxNcp = PetscMax(maxNcp, Ncp);
1723: }
1724: PetscCall(DMGetWorkArray(dm, maxNcp * dim, MPIU_REAL, &refcoord));
1725: PetscCall(DMGetWorkArray(dm, maxNcp * dim, MPIU_REAL, &pcoord));
1726: PetscCall(PetscFECreateTabulation(fe, 1, maxNcp, refcoord, 1, &tab));
1727: for (PetscInt c = cStart; c < cEnd; ++c) {
1728: PetscScalar *clPhi = NULL;
1729: PetscInt *points;
1730: PetscInt Ncp;
1732: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
1733: for (PetscInt cp = 0; cp < Ncp; ++cp)
1734: for (PetscInt d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
1735: {
1736: PetscCall(PetscFEGeomGetChunk(user->fegeom, c - cStart, c - cStart + 1, &chunkgeom));
1737: for (PetscInt i = 0; i < Ncp; ++i) {
1738: const PetscReal x0[3] = {-1., -1., -1.};
1739: CoordinatesRealToRef(dim, dim, x0, chunkgeom->v, chunkgeom->invJ, &pcoord[dim * i], &refcoord[dim * i]);
1740: }
1741: }
1742: PetscCall(PetscFEComputeTabulation(fe, Ncp, refcoord, 1, tab));
1743: PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
1744: for (PetscInt cp = 0; cp < Ncp; ++cp) {
1745: const PetscInt p = points[cp];
1747: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] = 0.;
1748: PetscCall(PetscFEInterpolateAtPoints_Static(fe, tab, clPhi, chunkgeom, cp, &E[p * dim]));
1749: PetscCall(PetscFEPushforward(fe, chunkgeom, 1, &E[p * dim]));
1750: for (PetscInt d = 0; d < dim; ++d) E[p * dim + d] *= -1.0;
1751: }
1752: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
1753: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
1754: }
1755: PetscCall(DMRestoreWorkArray(dm, maxNcp * dim, MPIU_REAL, &pcoord));
1756: PetscCall(DMRestoreWorkArray(dm, maxNcp * dim, MPIU_REAL, &refcoord));
1757: PetscCall(PetscTabulationDestroy(&tab));
1758: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1759: PetscCall(DMSwarmSortRestoreAccess(sw));
1760: PetscCall(DMRestoreLocalVector(dm, &locPhi));
1761: PetscCall(PetscFEGeomRestoreChunk(user->fegeom, 0, 1, &chunkgeom));
1762: PetscCall(PetscLogEventEnd(user->ETabEvent, snes, sw, 0, 0));
1763: PetscFunctionReturn(PETSC_SUCCESS);
1764: }
1766: static PetscErrorCode ComputeFieldAtParticles(SNES snes, DM sw)
1767: {
1768: AppCtx *user;
1769: Mat M_p;
1770: PetscReal *E;
1771: PetscInt dim, Np;
1773: PetscFunctionBegin;
1776: PetscCall(DMGetDimension(sw, &dim));
1777: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1778: PetscCall(DMGetApplicationContext(sw, &user));
1780: PetscCall(DMSwarmSetCellDMActive(sw, "moments"));
1781: // TODO: Could share sort context with space cellDM
1782: PetscCall(DMSwarmMigrate(sw, PETSC_FALSE));
1783: PetscCall(DMCreateMassMatrix(sw, user->dmPot, &M_p));
1784: PetscCall(DMSwarmSetCellDMActive(sw, "space"));
1786: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
1787: PetscCall(PetscArrayzero(E, Np * dim));
1788: user->validE = PETSC_TRUE;
1790: switch (user->em) {
1791: case EM_COULOMB:
1792: PetscCall(ComputeFieldAtParticles_Coulomb(snes, sw, E));
1793: break;
1794: case EM_PRIMAL:
1795: PetscCall(ComputeFieldAtParticles_Primal(snes, sw, M_p, E));
1796: break;
1797: case EM_MIXED:
1798: PetscCall(ComputeFieldAtParticles_Mixed(snes, sw, M_p, E));
1799: break;
1800: case EM_NONE:
1801: break;
1802: default:
1803: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No solver for electrostatic model %s", EMTypes[user->em]);
1804: }
1805: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
1806: PetscCall(MatDestroy(&M_p));
1807: PetscFunctionReturn(PETSC_SUCCESS);
1808: }
1810: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx)
1811: {
1812: DM sw;
1813: SNES snes = ((AppCtx *)ctx)->snes;
1814: const PetscScalar *u;
1815: PetscScalar *g;
1816: PetscReal *E, m_p = 1., q_p = -1.;
1817: PetscInt dim, d, Np, p;
1819: PetscFunctionBeginUser;
1820: PetscCall(TSGetDM(ts, &sw));
1821: PetscCall(ComputeFieldAtParticles(snes, sw));
1823: PetscCall(DMGetDimension(sw, &dim));
1824: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1825: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
1826: PetscCall(VecGetArrayRead(U, &u));
1827: PetscCall(VecGetArray(G, &g));
1828: Np /= 2 * dim;
1829: for (p = 0; p < Np; ++p) {
1830: for (d = 0; d < dim; ++d) {
1831: g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d];
1832: g[(p * 2 + 1) * dim + d] = q_p * E[p * dim + d] / m_p;
1833: }
1834: }
1835: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
1836: PetscCall(VecRestoreArrayRead(U, &u));
1837: PetscCall(VecRestoreArray(G, &g));
1838: PetscFunctionReturn(PETSC_SUCCESS);
1839: }
1841: /* J_{ij} = dF_i/dx_j
1842: J_p = ( 0 1)
1843: (-w^2 0)
1844: TODO Now there is another term with w^2 from the electric field. I think we will need to invert the operator.
1845: Perhaps we can approximate the Jacobian using only the cellwise P-P gradient from Coulomb
1846: */
1847: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, void *ctx)
1848: {
1849: DM sw;
1850: const PetscReal *coords, *vel;
1851: PetscInt dim, d, Np, p, rStart;
1853: PetscFunctionBeginUser;
1854: PetscCall(TSGetDM(ts, &sw));
1855: PetscCall(DMGetDimension(sw, &dim));
1856: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1857: PetscCall(MatGetOwnershipRange(J, &rStart, NULL));
1858: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1859: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&vel));
1860: Np /= 2 * dim;
1861: for (p = 0; p < Np; ++p) {
1862: const PetscReal x0 = coords[p * dim + 0];
1863: const PetscReal vy0 = vel[p * dim + 1];
1864: const PetscReal omega = vy0 / x0;
1865: PetscScalar vals[4] = {0., 1., -PetscSqr(omega), 0.};
1867: for (d = 0; d < dim; ++d) {
1868: const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
1869: PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES));
1870: }
1871: }
1872: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1873: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&vel));
1874: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1875: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1876: PetscFunctionReturn(PETSC_SUCCESS);
1877: }
1879: static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
1880: {
1881: AppCtx *user = (AppCtx *)ctx;
1882: DM sw;
1883: const PetscScalar *v;
1884: PetscScalar *xres;
1885: PetscInt Np, p, d, dim;
1887: PetscFunctionBeginUser;
1888: PetscCall(PetscLogEventBegin(user->RhsXEvent, ts, 0, 0, 0));
1889: PetscCall(TSGetDM(ts, &sw));
1890: PetscCall(DMGetDimension(sw, &dim));
1891: PetscCall(VecGetLocalSize(Xres, &Np));
1892: PetscCall(VecGetArrayRead(V, &v));
1893: PetscCall(VecGetArray(Xres, &xres));
1894: Np /= dim;
1895: for (p = 0; p < Np; ++p) {
1896: for (d = 0; d < dim; ++d) xres[p * dim + d] = v[p * dim + d];
1897: }
1898: PetscCall(VecRestoreArrayRead(V, &v));
1899: PetscCall(VecRestoreArray(Xres, &xres));
1900: PetscCall(PetscLogEventEnd(user->RhsXEvent, ts, 0, 0, 0));
1901: PetscFunctionReturn(PETSC_SUCCESS);
1902: }
1904: static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
1905: {
1906: DM sw;
1907: AppCtx *user = (AppCtx *)ctx;
1908: SNES snes = ((AppCtx *)ctx)->snes;
1909: const PetscScalar *x;
1910: PetscScalar *vres;
1911: PetscReal *E, m_p, q_p;
1912: PetscInt Np, p, dim, d;
1913: Parameter *param;
1915: PetscFunctionBeginUser;
1916: PetscCall(PetscLogEventBegin(user->RhsVEvent, ts, 0, 0, 0));
1917: PetscCall(TSGetDM(ts, &sw));
1918: PetscCall(ComputeFieldAtParticles(snes, sw));
1920: PetscCall(DMGetDimension(sw, &dim));
1921: PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
1922: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1923: m_p = user->masses[0] * param->m0;
1924: q_p = user->charges[0] * param->q0;
1925: PetscCall(VecGetLocalSize(Vres, &Np));
1926: PetscCall(VecGetArrayRead(X, &x));
1927: PetscCall(VecGetArray(Vres, &vres));
1928: Np /= dim;
1929: for (p = 0; p < Np; ++p) {
1930: for (d = 0; d < dim; ++d) vres[p * dim + d] = q_p * E[p * dim + d] / m_p;
1931: }
1932: PetscCall(VecRestoreArrayRead(X, &x));
1933: /*
1934: Synchronized, ordered output for parallel/sequential test cases.
1935: In the 1D (on the 2D mesh) case, every y component should be zero.
1936: */
1937: if (user->checkVRes) {
1938: PetscBool pr = user->checkVRes > 1 ? PETSC_TRUE : PETSC_FALSE;
1939: PetscInt step;
1941: PetscCall(TSGetStepNumber(ts, &step));
1942: if (pr) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "step: %" PetscInt_FMT "\n", step));
1943: for (PetscInt p = 0; p < Np; ++p) {
1944: if (pr) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Residual: %.12g %.12g\n", (double)PetscRealPart(vres[p * dim + 0]), (double)PetscRealPart(vres[p * dim + 1])));
1945: 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]));
1946: }
1947: if (pr) PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
1948: }
1949: PetscCall(VecRestoreArray(Vres, &vres));
1950: PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
1951: PetscCall(PetscLogEventEnd(user->RhsVEvent, ts, 0, 0, 0));
1952: PetscFunctionReturn(PETSC_SUCCESS);
1953: }
1955: static PetscErrorCode CreateSolution(TS ts)
1956: {
1957: DM sw;
1958: Vec u;
1959: PetscInt dim, Np;
1961: PetscFunctionBegin;
1962: PetscCall(TSGetDM(ts, &sw));
1963: PetscCall(DMGetDimension(sw, &dim));
1964: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1965: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
1966: PetscCall(VecSetBlockSize(u, dim));
1967: PetscCall(VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE));
1968: PetscCall(VecSetUp(u));
1969: PetscCall(TSSetSolution(ts, u));
1970: PetscCall(VecDestroy(&u));
1971: PetscFunctionReturn(PETSC_SUCCESS);
1972: }
1974: static PetscErrorCode SetProblem(TS ts)
1975: {
1976: AppCtx *user;
1977: DM sw;
1979: PetscFunctionBegin;
1980: PetscCall(TSGetDM(ts, &sw));
1981: PetscCall(DMGetApplicationContext(sw, (void **)&user));
1982: // Define unified system for (X, V)
1983: {
1984: Mat J;
1985: PetscInt dim, Np;
1987: PetscCall(DMGetDimension(sw, &dim));
1988: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1989: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
1990: PetscCall(MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE));
1991: PetscCall(MatSetBlockSize(J, 2 * dim));
1992: PetscCall(MatSetFromOptions(J));
1993: PetscCall(MatSetUp(J));
1994: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, user));
1995: PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, user));
1996: PetscCall(MatDestroy(&J));
1997: }
1998: /* Define split system for X and V */
1999: {
2000: Vec u;
2001: IS isx, isv, istmp;
2002: const PetscInt *idx;
2003: PetscInt dim, Np, rstart;
2005: PetscCall(TSGetSolution(ts, &u));
2006: PetscCall(DMGetDimension(sw, &dim));
2007: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2008: PetscCall(VecGetOwnershipRange(u, &rstart, NULL));
2009: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 0, 2, &istmp));
2010: PetscCall(ISGetIndices(istmp, &idx));
2011: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isx));
2012: PetscCall(ISRestoreIndices(istmp, &idx));
2013: PetscCall(ISDestroy(&istmp));
2014: PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 1, 2, &istmp));
2015: PetscCall(ISGetIndices(istmp, &idx));
2016: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isv));
2017: PetscCall(ISRestoreIndices(istmp, &idx));
2018: PetscCall(ISDestroy(&istmp));
2019: PetscCall(TSRHSSplitSetIS(ts, "position", isx));
2020: PetscCall(TSRHSSplitSetIS(ts, "momentum", isv));
2021: PetscCall(ISDestroy(&isx));
2022: PetscCall(ISDestroy(&isv));
2023: PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, user));
2024: PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, user));
2025: }
2026: PetscFunctionReturn(PETSC_SUCCESS);
2027: }
2029: static PetscErrorCode DMSwarmTSRedistribute(TS ts)
2030: {
2031: DM sw;
2032: Vec u;
2033: PetscReal t, maxt, dt;
2034: PetscInt n, maxn;
2036: PetscFunctionBegin;
2037: PetscCall(TSGetDM(ts, &sw));
2038: PetscCall(TSGetTime(ts, &t));
2039: PetscCall(TSGetMaxTime(ts, &maxt));
2040: PetscCall(TSGetTimeStep(ts, &dt));
2041: PetscCall(TSGetStepNumber(ts, &n));
2042: PetscCall(TSGetMaxSteps(ts, &maxn));
2044: PetscCall(TSReset(ts));
2045: PetscCall(TSSetDM(ts, sw));
2046: PetscCall(TSSetFromOptions(ts));
2047: PetscCall(TSSetTime(ts, t));
2048: PetscCall(TSSetMaxTime(ts, maxt));
2049: PetscCall(TSSetTimeStep(ts, dt));
2050: PetscCall(TSSetStepNumber(ts, n));
2051: PetscCall(TSSetMaxSteps(ts, maxn));
2053: PetscCall(CreateSolution(ts));
2054: PetscCall(SetProblem(ts));
2055: PetscCall(TSGetSolution(ts, &u));
2056: PetscFunctionReturn(PETSC_SUCCESS);
2057: }
2059: PetscErrorCode line(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
2060: {
2061: DM sw, cdm;
2062: PetscInt Np;
2063: PetscReal low[2], high[2];
2064: AppCtx *user = (AppCtx *)ctx;
2066: sw = user->swarm;
2067: PetscCall(DMSwarmGetCellDM(sw, &cdm));
2068: // Get the bounding box so we can equally space the particles
2069: PetscCall(DMGetLocalBoundingBox(cdm, low, high));
2070: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2071: // shift it by h/2 so nothing is initialized directly on a boundary
2072: x[0] = ((high[0] - low[0]) / Np) * (p + 0.5);
2073: x[1] = 0.;
2074: return PETSC_SUCCESS;
2075: }
2077: /*
2078: InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values.
2080: Input Parameters:
2081: + ts - The TS
2082: - useInitial - Flag to also set the initial conditions to the current coordinates and velocities and setup the problem
2084: Output Parameters:
2085: . u - The initialized solution vector
2087: Level: advanced
2089: .seealso: InitializeSolve()
2090: */
2091: static PetscErrorCode InitializeSolveAndSwarm(TS ts, PetscBool useInitial)
2092: {
2093: DM sw;
2094: Vec u, gc, gv;
2095: IS isx, isv;
2096: PetscInt dim;
2097: AppCtx *user;
2099: PetscFunctionBeginUser;
2100: PetscCall(TSGetDM(ts, &sw));
2101: PetscCall(DMGetApplicationContext(sw, &user));
2102: PetscCall(DMGetDimension(sw, &dim));
2103: if (useInitial) {
2104: PetscReal v0[2] = {1., 0.};
2105: if (user->perturbed_weights) {
2106: PetscCall(InitializeParticles_PerturbedWeights(sw, user));
2107: } else {
2108: PetscCall(DMSwarmComputeLocalSizeFromOptions(sw));
2109: PetscCall(DMSwarmInitializeCoordinates(sw));
2110: PetscCall(DMSwarmInitializeVelocitiesFromOptions(sw, v0));
2111: }
2112: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
2113: PetscCall(DMSwarmTSRedistribute(ts));
2114: }
2115: PetscCall(DMSetUp(sw));
2116: PetscCall(TSGetSolution(ts, &u));
2117: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
2118: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
2119: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
2120: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
2121: PetscCall(VecISCopy(u, isx, SCATTER_FORWARD, gc));
2122: PetscCall(VecISCopy(u, isv, SCATTER_FORWARD, gv));
2123: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
2124: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
2125: PetscFunctionReturn(PETSC_SUCCESS);
2126: }
2128: static PetscErrorCode InitializeSolve(TS ts, Vec u)
2129: {
2130: PetscFunctionBegin;
2131: PetscCall(TSSetSolution(ts, u));
2132: PetscCall(InitializeSolveAndSwarm(ts, PETSC_TRUE));
2133: PetscFunctionReturn(PETSC_SUCCESS);
2134: }
2136: static PetscErrorCode MigrateParticles(TS ts)
2137: {
2138: DM sw, cdm;
2139: const PetscReal *L;
2140: AppCtx *ctx;
2142: PetscFunctionBeginUser;
2143: PetscCall(TSGetDM(ts, &sw));
2144: PetscCall(DMGetApplicationContext(sw, &ctx));
2145: PetscCall(DMViewFromOptions(sw, NULL, "-migrate_view_pre"));
2146: {
2147: Vec u, gc, gv, position, momentum;
2148: IS isx, isv;
2149: PetscReal *pos, *mom;
2151: PetscCall(TSGetSolution(ts, &u));
2152: PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
2153: PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
2154: PetscCall(VecGetSubVector(u, isx, &position));
2155: PetscCall(VecGetSubVector(u, isv, &momentum));
2156: PetscCall(VecGetArray(position, &pos));
2157: PetscCall(VecGetArray(momentum, &mom));
2158: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
2159: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
2160: PetscCall(VecISCopy(u, isx, SCATTER_REVERSE, gc));
2161: PetscCall(VecISCopy(u, isv, SCATTER_REVERSE, gv));
2163: PetscCall(DMSwarmGetCellDM(sw, &cdm));
2164: PetscCall(DMGetPeriodicity(cdm, NULL, NULL, &L));
2165: if ((L[0] || L[1]) >= 0.) {
2166: PetscReal *x, *v, upper[3], lower[3];
2167: PetscInt Np, dim;
2169: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2170: PetscCall(DMGetDimension(cdm, &dim));
2171: PetscCall(DMGetBoundingBox(cdm, lower, upper));
2172: PetscCall(VecGetArray(gc, &x));
2173: PetscCall(VecGetArray(gv, &v));
2174: for (PetscInt p = 0; p < Np; ++p) {
2175: for (PetscInt d = 0; d < dim; ++d) {
2176: if (pos[p * dim + d] < lower[d]) {
2177: x[p * dim + d] = pos[p * dim + d] + (upper[d] - lower[d]);
2178: } else if (pos[p * dim + d] > upper[d]) {
2179: x[p * dim + d] = pos[p * dim + d] - (upper[d] - lower[d]);
2180: } else {
2181: x[p * dim + d] = pos[p * dim + d];
2182: }
2183: 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]);
2184: v[p * dim + d] = mom[p * dim + d];
2185: }
2186: }
2187: PetscCall(VecRestoreArray(gc, &x));
2188: PetscCall(VecRestoreArray(gv, &v));
2189: }
2190: PetscCall(VecRestoreArray(position, &pos));
2191: PetscCall(VecRestoreArray(momentum, &mom));
2192: PetscCall(VecRestoreSubVector(u, isx, &position));
2193: PetscCall(VecRestoreSubVector(u, isv, &momentum));
2194: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
2195: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
2196: }
2197: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
2198: PetscInt step;
2200: PetscCall(TSGetStepNumber(ts, &step));
2201: if (!(step % ctx->remapFreq)) {
2202: // Monitor electric field before we destroy it
2203: PetscReal ptime;
2204: PetscInt step;
2206: PetscCall(TSGetStepNumber(ts, &step));
2207: PetscCall(TSGetTime(ts, &ptime));
2208: if (ctx->efield_monitor) PetscCall(MonitorEField(ts, step, ptime, NULL, ctx));
2209: if (ctx->poisson_monitor) PetscCall(MonitorPoisson(ts, step, ptime, NULL, ctx));
2210: PetscCall(DMSwarmRemap(sw));
2211: ctx->validE = PETSC_FALSE;
2212: }
2213: // This MUST come last, since it recreates the subswarms and they must DMClone() the new swarm
2214: PetscCall(DMSwarmTSRedistribute(ts));
2215: PetscCall(InitializeSolveAndSwarm(ts, PETSC_FALSE));
2216: PetscFunctionReturn(PETSC_SUCCESS);
2217: }
2219: int main(int argc, char **argv)
2220: {
2221: DM dm, sw;
2222: TS ts;
2223: Vec u;
2224: PetscReal dt;
2225: PetscInt maxn;
2226: AppCtx user;
2228: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
2229: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
2230: PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &user.bag));
2231: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
2232: PetscCall(CreatePoisson(dm, &user));
2233: PetscCall(CreateSwarm(dm, &user, &sw));
2234: PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
2235: PetscCall(InitializeConstants(sw, &user));
2236: PetscCall(DMSetApplicationContext(sw, &user));
2238: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
2239: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
2240: PetscCall(TSSetDM(ts, sw));
2241: PetscCall(TSSetMaxTime(ts, 0.1));
2242: PetscCall(TSSetTimeStep(ts, 0.00001));
2243: PetscCall(TSSetMaxSteps(ts, 100));
2244: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
2246: if (user.efield_monitor) PetscCall(TSMonitorSet(ts, MonitorEField, &user, NULL));
2247: if (user.moment_monitor) PetscCall(TSMonitorSet(ts, MonitorMoments, &user, NULL));
2248: if (user.initial_monitor) PetscCall(TSMonitorSet(ts, MonitorInitialConditions, &user, NULL));
2249: if (user.positions_monitor) PetscCall(TSMonitorSet(ts, MonitorPositions_2D, &user, NULL));
2250: if (user.poisson_monitor) PetscCall(TSMonitorSet(ts, MonitorPoisson, &user, NULL));
2251: if (user.velocity_monitor >= 0) PetscCall(TSMonitorSet(ts, MonitorVelocity, &user, NULL));
2253: PetscCall(TSSetFromOptions(ts));
2254: PetscCall(TSGetTimeStep(ts, &dt));
2255: PetscCall(TSGetMaxSteps(ts, &maxn));
2256: user.steps = maxn;
2257: user.stepSize = dt;
2258: PetscCall(SetupContext(dm, sw, &user));
2259: PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
2260: PetscCall(TSSetPostStep(ts, MigrateParticles));
2261: PetscCall(CreateSolution(ts));
2262: PetscCall(TSGetSolution(ts, &u));
2263: PetscCall(TSComputeInitialCondition(ts, u));
2264: PetscCall(CheckNonNegativeWeights(sw, &user));
2265: PetscCall(TSSolve(ts, NULL));
2267: PetscCall(SNESDestroy(&user.snes));
2268: PetscCall(DMDestroy(&user.dmPot));
2269: PetscCall(ISDestroy(&user.isPot));
2270: PetscCall(MatDestroy(&user.M));
2271: PetscCall(PetscFEGeomDestroy(&user.fegeom));
2272: PetscCall(TSDestroy(&ts));
2273: PetscCall(DMDestroy(&sw));
2274: PetscCall(DMDestroy(&dm));
2275: PetscCall(DestroyContext(&user));
2276: PetscCall(PetscFinalize());
2277: return 0;
2278: }
2280: /*TEST
2282: build:
2283: requires: !complex double
2285: # This tests that we can put particles in a box and compute the Coulomb force
2286: # Recommend -draw_size 500,500
2287: testset:
2288: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
2289: args: -dm_plex_dim 1 -dm_plex_simplex 0 -dm_plex_box_faces 20 \
2290: -dm_plex_box_lower 0 -dm_plex_box_upper 12.5664 -dm_plex_box_bd periodic \
2291: -vdm_plex_simplex 0 \
2292: -dm_swarm_coordinate_density constant -dm_swarm_num_particles 100 \
2293: -sigma 1.0e-8 -timeScale 2.0e-14 \
2294: -ts_type basicsymplectic -ts_basicsymplectic_type 1 \
2295: -output_step 50 -check_vel_res
2296: test:
2297: suffix: none_1d
2298: requires:
2299: args: -em_type none -error
2300: test:
2301: suffix: coulomb_1d
2302: args: -em_type coulomb
2304: # for viewers
2305: #-ts_monitor_sp_swarm_phase -ts_monitor_sp_swarm -em_snes_monitor -ts_monitor_sp_swarm_multi_species 0 -ts_monitor_sp_swarm_retain 0
2306: testset:
2307: nsize: {{1 2}}
2308: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
2309: args: -dm_plex_dim 1 -dm_plex_simplex 0 -dm_plex_box_faces 36 \
2310: -dm_plex_box_lower 0. -dm_plex_box_upper 12.5664 -dm_plex_box_bd periodic \
2311: -vdm_plex_dim 1 -vdm_plex_simplex 0 -vdm_plex_box_faces 10 \
2312: -vdm_plex_box_lower -3 -vdm_plex_box_upper 3 \
2313: -dm_swarm_num_species 1 -twostream -charges -1.,1. -sigma 1.0e-8 \
2314: -cosine_coefficients 0.01,0.5 -perturbed_weights -total_weight 1. \
2315: -ts_type basicsymplectic -ts_basicsymplectic_type 2 \
2316: -ts_dt 0.01 -ts_max_time 5 -ts_max_steps 10 \
2317: -em_snes_atol 1.e-15 -em_snes_error_if_not_converged -em_ksp_error_if_not_converged \
2318: -output_step 1 -check_vel_res -dm_swarm_print_coords 1 -dm_swarm_print_weights 1
2319: test:
2320: suffix: two_stream_c0
2321: args: -em_type primal -petscfe_default_quadrature_order 2 -petscspace_degree 2 -em_pc_type svd
2322: test:
2323: suffix: two_stream_rt
2324: requires: superlu_dist
2325: args: -em_type mixed \
2326: -potential_petscspace_degree 0 \
2327: -potential_petscdualspace_lagrange_use_moments \
2328: -potential_petscdualspace_lagrange_moment_order 2 \
2329: -field_petscspace_degree 1 -field_petscfe_default_quadrature_order 1 \
2330: -em_snes_error_if_not_converged \
2331: -em_ksp_type preonly -em_ksp_error_if_not_converged \
2332: -em_pc_type fieldsplit -em_pc_fieldsplit_type schur \
2333: -em_pc_fieldsplit_schur_fact_type full -em_pc_fieldsplit_schur_precondition full \
2334: -em_fieldsplit_field_pc_type lu \
2335: -em_fieldsplit_field_pc_factor_mat_solver_type superlu_dist \
2336: -em_fieldsplit_potential_pc_type svd
2338: # For an eyeball check, we use
2339: # -ts_max_steps 1000 -dm_plex_box_faces 10,1 -vdm_plex_box_faces 2000 -efield_monitor
2340: # For verification, we use
2341: # -ts_max_steps 1000 -dm_plex_box_faces 100,1 -vdm_plex_box_faces 8000 -dm_swarm_num_particles 800000 -monitor_efield
2342: # -ts_monitor_sp_swarm_multi_species 0 -ts_monitor_sp_swarm_retain 0 -ts_monitor_sp_swarm_phase 1 -draw_size 500,500
2343: testset:
2344: nsize: {{1 2}}
2345: requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
2346: args: -dm_plex_dim 1 -dm_plex_simplex 0 -dm_plex_box_faces 10 \
2347: -dm_plex_box_lower 0. -dm_plex_box_upper 12.5664 -dm_plex_box_bd periodic \
2348: -vdm_plex_dim 1 -vdm_plex_simplex 0 -vdm_plex_box_faces 10 \
2349: -vdm_plex_box_lower -10 -vdm_plex_box_upper 10 \
2350: -vpetscspace_degree 2 -vdm_plex_hash_location \
2351: -dm_swarm_num_species 1 -charges -1.,1. \
2352: -cosine_coefficients 0.01,0.5 -perturbed_weights -total_weight 1. \
2353: -ts_type basicsymplectic -ts_basicsymplectic_type 1 \
2354: -ts_dt 0.03 -ts_max_time 500 -ts_max_steps 1 \
2355: -em_snes_atol 1.e-12 -em_snes_error_if_not_converged -em_ksp_error_if_not_converged \
2356: -output_step 1 -check_vel_res -dm_swarm_print_coords 1 -dm_swarm_print_weights 1
2358: test:
2359: suffix: uniform_equilibrium_1d
2360: args: -cosine_coefficients 0.0,0.5 -em_type primal -petscspace_degree 1 -em_pc_type svd
2361: test:
2362: suffix: uniform_equilibrium_1d_real
2363: args: -dm_plex_dim 1 -dm_plex_simplex 1 -dm_plex_box_faces 10 \
2364: -dm_plex_box_lower 0. -dm_plex_box_upper 12.5664 -dm_plex_box_bd periodic \
2365: -cosine_coefficients 0.0 -em_type primal -petscspace_degree 1 -em_pc_type svd
2366: test:
2367: suffix: uniform_primal_1d
2368: args: -em_type primal -petscspace_degree 1 -em_pc_type svd
2369: test:
2370: suffix: uniform_primal_1d_real
2371: args: -dm_plex_dim 1 -dm_plex_simplex 1 -dm_plex_box_faces 10 \
2372: -dm_plex_box_lower 0. -dm_plex_box_upper 12.5664 -dm_plex_box_bd periodic \
2373: -cosine_coefficients 0.01 -em_type primal -petscspace_degree 1 -em_pc_type svd
2374: # NOT WORKING -ftop_pc_type bjacobi -ftop_sub_pc_type lu -ftop_sub_pc_factor_shift_type nonzero
2375: test:
2376: suffix: uniform_primal_1d_real_pfak
2377: nsize: 1
2378: args: -dm_plex_dim 1 -dm_plex_simplex 1 -dm_plex_box_faces 10 \
2379: -dm_plex_box_lower 0. -dm_plex_box_upper 12.5664 -dm_plex_box_bd periodic \
2380: -remap_dm_plex_dim 2 -remap_dm_plex_simplex 0 -remap_dm_plex_box_faces 10,10 -remap_dm_plex_box_bd periodic,none \
2381: -remap_dm_plex_box_lower 0.,-10. -remap_dm_plex_box_upper 12.5664,10. \
2382: -remap_petscspace_degree 2 -remap_dm_plex_hash_location \
2383: -remap_freq 1 -dm_swarm_remap_type pfak \
2384: -ftop_ksp_type lsqr -ftop_pc_type none -ftop_ksp_rtol 1.e-14 \
2385: -ptof_pc_type lu \
2386: -cosine_coefficients 0.01 -em_type primal -petscspace_degree 1 -em_pc_type svd -em_proj_pc_type lu
2387: test:
2388: requires: superlu_dist
2389: suffix: uniform_mixed_1d
2390: args: -em_type mixed \
2391: -potential_petscspace_degree 0 \
2392: -potential_petscdualspace_lagrange_use_moments \
2393: -potential_petscdualspace_lagrange_moment_order 2 \
2394: -field_petscspace_degree 1 \
2395: -em_snes_error_if_not_converged \
2396: -em_ksp_type preonly -em_ksp_error_if_not_converged \
2397: -em_pc_type fieldsplit -em_pc_fieldsplit_type schur \
2398: -em_pc_fieldsplit_schur_fact_type full -em_pc_fieldsplit_schur_precondition full \
2399: -em_fieldsplit_field_pc_type lu \
2400: -em_fieldsplit_field_pc_factor_mat_solver_type superlu_dist \
2401: -em_fieldsplit_potential_pc_type svd
2403: TEST*/