1 #include "petscdm.h"
2 static char help[] = "Tests for particle initialization using the KS test\n\n";
3
4 #include <petscdmswarm.h>
5 #include <petscdmplex.h>
6 #include <petsc/private/petscfeimpl.h>
7
8 /*
9 View the KS test using
10
11 -ks_monitor draw -draw_size 500,500 -draw_pause 3
12
13 Set total number to n0 / Mp = 3e14 / 1e12 = 300 using -dm_swarm_num_particles 300
14
15 */
16
17 #define BOLTZMANN_K 1.380649e-23 /* J/K */
18
19 typedef struct {
20 PetscReal mass[2]; /* Electron, Sr+ Mass [kg] */
21 PetscReal T[2]; /* Electron, Ion Temperature [K] */
22 PetscReal v0[2]; /* Species mean velocity in 1D */
23 } AppCtx;
24
ProcessOptions(MPI_Comm comm,AppCtx * options)25 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
26 {
27 PetscFunctionBegin;
28 options->mass[0] = 9.10938356e-31; /* Electron Mass [kg] */
29 options->mass[1] = 87.62 * 1.66054e-27; /* Sr+ Mass [kg] */
30 options->T[0] = 1.; /* Electron Temperature [K] */
31 options->T[1] = 25.; /* Sr+ Temperature [K] */
32 options->v0[0] = PetscSqrtReal(BOLTZMANN_K * options->T[0] / options->mass[0]); /* electron mean velocity in 1D */
33 options->v0[1] = PetscSqrtReal(BOLTZMANN_K * options->T[1] / options->mass[1]); /* ion mean velocity in 1D */
34
35 PetscOptionsBegin(comm, "", "KS Test Options", "DMPLEX");
36 PetscOptionsEnd();
37 PetscFunctionReturn(PETSC_SUCCESS);
38 }
39
CreateMesh(MPI_Comm comm,DM * dm)40 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
41 {
42 PetscFunctionBeginUser;
43 PetscCall(DMCreate(comm, dm));
44 PetscCall(DMSetType(*dm, DMPLEX));
45 PetscCall(DMSetFromOptions(*dm));
46 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
CreateSwarm(DM dm,AppCtx * user,DM * sw)50 static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
51 {
52 PetscInt dim;
53
54 PetscFunctionBeginUser;
55 PetscCall(DMGetDimension(dm, &dim));
56 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
57 PetscCall(DMSetType(*sw, DMSWARM));
58 PetscCall(DMSetDimension(*sw, dim));
59 PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
60 PetscCall(DMSwarmSetCellDM(*sw, dm));
61 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
62 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
63 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT));
64 PetscCall(DMSwarmFinalizeFieldRegister(*sw));
65 PetscCall(DMSwarmComputeLocalSizeFromOptions(*sw));
66 PetscCall(DMSwarmInitializeCoordinates(*sw));
67 PetscCall(DMSwarmInitializeVelocitiesFromOptions(*sw, user->v0));
68 PetscCall(DMSetFromOptions(*sw));
69 PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
70 PetscCall(DMViewFromOptions(*sw, NULL, "-swarm_view"));
71 PetscFunctionReturn(PETSC_SUCCESS);
72 }
73
TestDistribution(DM sw,PetscReal confidenceLevel,AppCtx * user)74 static PetscErrorCode TestDistribution(DM sw, PetscReal confidenceLevel, AppCtx *user)
75 {
76 DM dm;
77 Vec locx, locv, locw;
78 PetscProbFn *cdf;
79 PetscReal alpha, gmin, gmax;
80 PetscInt dim;
81 MPI_Comm comm;
82
83 PetscFunctionBeginUser;
84 PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
85 PetscCall(DMSwarmGetCellDM(sw, &dm));
86 PetscCall(DMGetBoundingBox(dm, &gmin, &gmax));
87 PetscCall(DMGetDimension(sw, &dim));
88 switch (dim) {
89 case 1:
90 cdf = PetscCDFMaxwellBoltzmann1D;
91 break;
92 case 2:
93 cdf = PetscCDFMaxwellBoltzmann2D;
94 break;
95 case 3:
96 cdf = PetscCDFMaxwellBoltzmann3D;
97 break;
98 default:
99 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported", dim);
100 }
101 PetscCall(DMSwarmCreateLocalVectorFromField(sw, DMSwarmPICField_coor, &locx));
102 PetscCall(DMSwarmCreateLocalVectorFromField(sw, "velocity", &locv));
103 PetscCall(DMSwarmCreateLocalVectorFromField(sw, "w_q", &locw));
104 PetscCall(VecViewFromOptions(locx, NULL, "-x_view"));
105 PetscCall(VecViewFromOptions(locv, NULL, "-v_view"));
106 PetscCall(VecViewFromOptions(locw, NULL, "-w_view"));
107 // Need to rescale coordinates to compare to distribution
108 PetscCall(VecScale(locx, 1. / gmax));
109 PetscCall(PetscProbComputeKSStatisticWeighted(locx, locw, PetscCDFConstant1D, &alpha));
110 if (alpha < confidenceLevel) PetscCall(PetscPrintf(comm, "The KS test for X rejects the null hypothesis at level %.2g (%.2g)\n", (double)confidenceLevel, (double)alpha));
111 else PetscCall(PetscPrintf(comm, "The KS test for X accepts the null hypothesis at level %.2g (%.2g)\n", (double)confidenceLevel, (double)alpha));
112 PetscCall(PetscProbComputeKSStatisticWeighted(locv, locw, PetscCDFGaussian1D, &alpha));
113 if (alpha < confidenceLevel) PetscCall(PetscPrintf(comm, "The KS test for V rejects the null hypothesis at level %.2g (%.2g)\n", (double)confidenceLevel, (double)alpha));
114 else PetscCall(PetscPrintf(comm, "The KS test for V accepts the null hypothesis at level %.2g (%.2g)\n", (double)confidenceLevel, (double)alpha));
115 PetscCall(PetscProbComputeKSStatisticMagnitude(locv, cdf, &alpha));
116 if (alpha < confidenceLevel) PetscCall(PetscPrintf(comm, "The KS test for |V| rejects the null hypothesis at level %.2g (%.2g)\n", (double)confidenceLevel, (double)alpha));
117 else PetscCall(PetscPrintf(comm, "The KS test for |V| accepts the null hypothesis at level %.2g (%.2g)\n", (double)confidenceLevel, (double)alpha));
118 PetscCall(DMSwarmDestroyLocalVectorFromField(sw, DMSwarmPICField_coor, &locx));
119 PetscCall(DMSwarmDestroyLocalVectorFromField(sw, "velocity", &locv));
120 PetscCall(DMSwarmDestroyLocalVectorFromField(sw, "w_q", &locw));
121 PetscFunctionReturn(PETSC_SUCCESS);
122 }
123
main(int argc,char ** argv)124 int main(int argc, char **argv)
125 {
126 DM dm, sw;
127 AppCtx user;
128
129 PetscFunctionBeginUser;
130 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
131 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
132 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm));
133 PetscCall(CreateSwarm(dm, &user, &sw));
134 PetscCall(TestDistribution(sw, 0.01, &user));
135 PetscCall(DMDestroy(&sw));
136 PetscCall(DMDestroy(&dm));
137 PetscCall(PetscFinalize());
138 return 0;
139 }
140
141 /*TEST
142
143 build:
144 requires: !complex double
145
146 testset:
147 requires: ks !complex
148 args: -dm_plex_dim 1 -dm_plex_box_lower -4 -dm_plex_box_upper 4 -dm_plex_box_faces 40 -dm_swarm_num_particles 10000
149
150 test:
151 suffix: 0_constant
152 args: -dm_swarm_coordinate_density constant
153
154 test:
155 suffix: 0_gaussian
156 args: -dm_swarm_coordinate_density gaussian
157
158 TEST*/
159