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: options->mass[0] = 9.10938356e-31; /* Electron Mass [kg] */
27: options->mass[1] = 87.62 * 1.66054e-27; /* Sr+ Mass [kg] */
28: options->T[0] = 1.; /* Electron Temperature [K] */
29: options->T[1] = 25.; /* Sr+ Temperature [K] */
30: options->v0[0] = PetscSqrtReal(BOLTZMANN_K * options->T[0] / options->mass[0]); /* electron mean velocity in 1D */
31: options->v0[1] = PetscSqrtReal(BOLTZMANN_K * options->T[1] / options->mass[1]); /* ion mean velocity in 1D */
33: PetscOptionsBegin(comm, "", "KS Test Options", "DMPLEX");
34: PetscOptionsEnd();
35: return 0;
36: }
38: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
39: {
41: DMCreate(comm, dm);
42: DMSetType(*dm, DMPLEX);
43: DMSetFromOptions(*dm);
44: DMViewFromOptions(*dm, NULL, "-dm_view");
45: return 0;
46: }
48: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
49: {
50: PetscInt dim;
53: DMGetDimension(dm, &dim);
54: DMCreate(PetscObjectComm((PetscObject)dm), sw);
55: DMSetType(*sw, DMSWARM);
56: DMSetDimension(*sw, dim);
57: DMSwarmSetType(*sw, DMSWARM_PIC);
58: DMSwarmSetCellDM(*sw, dm);
59: DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR);
60: DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL);
61: DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT);
62: DMSwarmFinalizeFieldRegister(*sw);
63: DMSwarmComputeLocalSizeFromOptions(*sw);
64: DMSwarmInitializeCoordinates(*sw);
65: DMSwarmInitializeVelocitiesFromOptions(*sw, user->v0);
66: DMSetFromOptions(*sw);
67: PetscObjectSetName((PetscObject)*sw, "Particles");
68: DMViewFromOptions(*sw, NULL, "-swarm_view");
69: return 0;
70: }
72: static PetscErrorCode TestDistribution(DM sw, PetscReal confidenceLevel, AppCtx *user)
73: {
74: Vec locv;
75: PetscProbFunc cdf;
76: PetscReal alpha;
77: PetscInt dim;
78: MPI_Comm comm;
81: PetscObjectGetComm((PetscObject)sw, &comm);
82: DMGetDimension(sw, &dim);
83: switch (dim) {
84: case 1:
85: cdf = PetscCDFMaxwellBoltzmann1D;
86: break;
87: case 2:
88: cdf = PetscCDFMaxwellBoltzmann2D;
89: break;
90: case 3:
91: cdf = PetscCDFMaxwellBoltzmann3D;
92: break;
93: default:
94: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported", dim);
95: }
96: DMSwarmCreateLocalVectorFromField(sw, "velocity", &locv);
97: PetscProbComputeKSStatistic(locv, cdf, &alpha);
98: DMSwarmDestroyLocalVectorFromField(sw, "velocity", &locv);
99: if (alpha < confidenceLevel) PetscPrintf(comm, "The KS test accepts the null hypothesis at level %.2g\n", (double)confidenceLevel);
100: else PetscPrintf(comm, "The KS test rejects the null hypothesis at level %.2g (%.2g)\n", (double)confidenceLevel, (double)alpha);
101: return 0;
102: }
104: int main(int argc, char **argv)
105: {
106: DM dm, sw;
107: AppCtx user;
110: PetscInitialize(&argc, &argv, NULL, help);
111: ProcessOptions(PETSC_COMM_WORLD, &user);
112: CreateMesh(PETSC_COMM_WORLD, &dm);
113: CreateSwarm(dm, &user, &sw);
114: TestDistribution(sw, 0.05, &user);
115: DMDestroy(&sw);
116: DMDestroy(&dm);
117: PetscFinalize();
118: return 0;
119: }
121: /*TEST
123: test:
124: suffix: 0
125: requires: ks !complex
126: 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}}
128: TEST*/