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

 34:   PetscOptionsBegin(comm, "", "KS Test Options", "DMPLEX");
 35:   PetscOptionsEnd();
 36:   PetscFunctionReturn(PETSC_SUCCESS);
 37: }

 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: }

 49: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
 50: {
 51:   PetscInt dim;

 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: }

 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;

 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: }

105: int main(int argc, char **argv)
106: {
107:   DM     dm, sw;
108:   AppCtx user;

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: }

122: /*TEST

124:   build:
125:     requires: !complex double

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

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

137: TEST*/