1 static char help[] = "Tests for particle initialization using the KS test\n\n"; 2 3 #include <petscdmswarm.h> 4 #include <petscdmplex.h> 5 #include <petsc/private/petscfeimpl.h> 6 7 /* 8 View the KS test using 9 10 -ks_monitor draw -draw_size 500,500 -draw_pause 3 11 12 Set total number to n0 / Mp = 3e14 / 1e12 = 300 using -dm_swarm_num_particles 300 13 14 */ 15 16 #define BOLTZMANN_K 1.380649e-23 /* J/K */ 17 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; 23 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 */ 33 34 PetscOptionsBegin(comm, "", "KS Test Options", "DMPLEX"); 35 PetscOptionsEnd(); 36 PetscFunctionReturn(PETSC_SUCCESS); 37 } 38 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 } 48 49 static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw) 50 { 51 PetscInt dim; 52 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 } 72 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; 80 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 } 104 105 int main(int argc, char **argv) 106 { 107 DM dm, sw; 108 AppCtx user; 109 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 } 121 122 /*TEST 123 124 build: 125 requires: !complex double 126 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}} 131 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}} 136 137 TEST*/ 138