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*/