Actual source code: ex8.c
1: #include "petscdm.h"
2: static char help[] = "Tests for particle initialization using the KS test\n\n";
4: #include <petscdmswarm.h>
5: #include <petscdmplex.h>
6: #include <petsc/private/petscfeimpl.h>
8: /*
9: View the KS test using
11: -ks_monitor draw -draw_size 500,500 -draw_pause 3
13: Set total number to n0 / Mp = 3e14 / 1e12 = 300 using -dm_swarm_num_particles 300
15: */
17: #define BOLTZMANN_K 1.380649e-23 /* J/K */
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;
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 */
35: PetscOptionsBegin(comm, "", "KS Test Options", "DMPLEX");
36: PetscOptionsEnd();
37: PetscFunctionReturn(PETSC_SUCCESS);
38: }
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: }
50: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
51: {
52: PetscInt dim;
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: }
74: static PetscErrorCode TestDistribution(DM sw, PetscReal confidenceLevel, AppCtx *user)
75: {
76: DM dm;
77: Vec locx, locv, locw;
78: PetscProbFunc cdf;
79: PetscReal alpha, gmin, gmax;
80: PetscInt dim;
81: MPI_Comm comm;
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: }
124: int main(int argc, char **argv)
125: {
126: DM dm, sw;
127: AppCtx user;
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: }
141: /*TEST
143: build:
144: requires: !complex double
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
150: test:
151: suffix: 0_constant
152: args: -dm_swarm_coordinate_density constant
154: test:
155: suffix: 0_gaussian
156: args: -dm_swarm_coordinate_density gaussian
158: TEST*/