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