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