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 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 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 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 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; 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 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