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