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