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 **)&param));
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 **)&param));
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 **)&param));
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*/