Actual source code: ex8.c
1: static char help[] = "Tests for particle initialization using the KS test\n\n";
3: #include <petscdmswarm.h>
4: #include <petscdmplex.h>
5: #include <petsc/private/petscfeimpl.h>
7: /*
8: View the KS test using
10: -ks_monitor draw -draw_size 500,500 -draw_pause 3
12: Set total number to n0 / Mp = 3e14 / 1e12 = 300 using -dm_swarm_num_particles 300
14: */
16: #define BOLTZMANN_K 1.380649e-23 /* J/K */
18: typedef struct {
19: PetscReal mass[2]; /* Electron, Sr+ Mass [kg] */
20: PetscReal T[2]; /* Electron, Ion Temperature [K] */
21: PetscReal v0[2]; /* Species mean velocity in 1D */
22: } AppCtx;
24: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
25: {
26: PetscFunctionBegin;
27: options->mass[0] = 9.10938356e-31; /* Electron Mass [kg] */
28: options->mass[1] = 87.62 * 1.66054e-27; /* Sr+ Mass [kg] */
29: options->T[0] = 1.; /* Electron Temperature [K] */
30: options->T[1] = 25.; /* Sr+ Temperature [K] */
31: options->v0[0] = PetscSqrtReal(BOLTZMANN_K * options->T[0] / options->mass[0]); /* electron mean velocity in 1D */
32: options->v0[1] = PetscSqrtReal(BOLTZMANN_K * options->T[1] / options->mass[1]); /* ion mean velocity in 1D */
34: PetscOptionsBegin(comm, "", "KS Test Options", "DMPLEX");
35: PetscOptionsEnd();
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
40: {
41: PetscFunctionBeginUser;
42: PetscCall(DMCreate(comm, dm));
43: PetscCall(DMSetType(*dm, DMPLEX));
44: PetscCall(DMSetFromOptions(*dm));
45: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
50: {
51: PetscInt dim;
53: PetscFunctionBeginUser;
54: PetscCall(DMGetDimension(dm, &dim));
55: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
56: PetscCall(DMSetType(*sw, DMSWARM));
57: PetscCall(DMSetDimension(*sw, dim));
58: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
59: PetscCall(DMSwarmSetCellDM(*sw, dm));
60: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
61: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
62: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT));
63: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
64: PetscCall(DMSwarmComputeLocalSizeFromOptions(*sw));
65: PetscCall(DMSwarmInitializeCoordinates(*sw));
66: PetscCall(DMSwarmInitializeVelocitiesFromOptions(*sw, user->v0));
67: PetscCall(DMSetFromOptions(*sw));
68: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
69: PetscCall(DMViewFromOptions(*sw, NULL, "-swarm_view"));
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: static PetscErrorCode TestDistribution(DM sw, PetscReal confidenceLevel, AppCtx *user)
74: {
75: Vec locv;
76: PetscProbFunc cdf;
77: PetscReal alpha;
78: PetscInt dim;
79: MPI_Comm comm;
81: PetscFunctionBeginUser;
82: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
83: PetscCall(DMGetDimension(sw, &dim));
84: switch (dim) {
85: case 1:
86: cdf = PetscCDFMaxwellBoltzmann1D;
87: break;
88: case 2:
89: cdf = PetscCDFMaxwellBoltzmann2D;
90: break;
91: case 3:
92: cdf = PetscCDFMaxwellBoltzmann3D;
93: break;
94: default:
95: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported", dim);
96: }
97: PetscCall(DMSwarmCreateLocalVectorFromField(sw, "velocity", &locv));
98: PetscCall(PetscProbComputeKSStatistic(locv, cdf, &alpha));
99: PetscCall(DMSwarmDestroyLocalVectorFromField(sw, "velocity", &locv));
100: if (alpha < confidenceLevel) PetscCall(PetscPrintf(comm, "The KS test accepts the null hypothesis at level %.2g\n", (double)confidenceLevel));
101: else PetscCall(PetscPrintf(comm, "The KS test rejects the null hypothesis at level %.2g (%.2g)\n", (double)confidenceLevel, (double)alpha));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: int main(int argc, char **argv)
106: {
107: DM dm, sw;
108: AppCtx user;
110: PetscFunctionBeginUser;
111: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
112: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
113: PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm));
114: PetscCall(CreateSwarm(dm, &user, &sw));
115: PetscCall(TestDistribution(sw, 0.05, &user));
116: PetscCall(DMDestroy(&sw));
117: PetscCall(DMDestroy(&dm));
118: PetscCall(PetscFinalize());
119: return 0;
120: }
122: /*TEST
124: build:
125: requires: !complex double
127: test:
128: suffix: 0
129: requires: ks !complex
130: args: -dm_plex_dim 1 -dm_plex_box_lower -1 -dm_plex_box_upper 1 -dm_swarm_num_particles 375 -dm_swarm_coordinate_density {{constant gaussian}}
132: test:
133: suffix: 1
134: requires: ks !complex
135: args: -dm_plex_dim 1 -dm_plex_box_lower -10 -dm_plex_box_upper 10 -dm_plex_box_faces 20 -dm_swarm_num_particles 375 -dm_swarm_coordinate_density {{constant gaussian}}
137: TEST*/