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