Actual source code: ex35.c

  1: static char help[] = "Test of Colorized Scatter Plot.\n";

  3: #include "petscdraw.h"
  4: #include "petscvec.h"
  5: #include "petscis.h"

  7: typedef struct {
  8:   PetscInt Np; /* total number of particles */
  9:   PetscInt dim;
 10:   PetscInt dim_inp;
 11: } AppCtx;

 13: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 14: {
 16:   options->dim     = 2;
 17:   options->dim_inp = 2;
 18:   options->Np      = 100;
 19:   PetscOptionsBegin(comm, "", "Test of colorized scatter plot", "");
 20:   PetscOptionsInt("-Np", "Number of particles", "ex35.c", options->Np, &options->Np, PETSC_NULL);
 21:   PetscOptionsInt("-dim", "Number of dimensions", "ex35.c", options->dim_inp, &options->dim_inp, PETSC_NULL);
 22:   PetscOptionsEnd();
 23:   return 0;
 24: }

 26: /*
 27:   ref: http://www.mimirgames.com/articles/programming/approximations-of-the-inverse-error-function/
 28: */
 29: PetscReal erfinv(PetscReal x)
 30: {
 31:   PetscReal *ck, r = 0.;
 32:   PetscInt   maxIter = 100;

 34:   PetscCalloc1(maxIter, &ck);
 35:   ck[0] = 1;
 36:   r     = ck[0] * ((PetscSqrtReal(PETSC_PI) / 2.) * x);
 37:   for (PetscInt k = 1; k < maxIter; ++k) {
 38:     for (PetscInt m = 0; m <= k - 1; ++m) {
 39:       PetscReal denom = (m + 1.) * (2. * m + 1.);
 40:       ck[k] += (ck[m] * ck[k - 1 - m]) / denom;
 41:     }
 42:     PetscReal temp = 2. * k + 1.;
 43:     r += (ck[k] / temp) * PetscPowReal((PetscSqrtReal(PETSC_PI) / 2.) * x, 2. * k + 1.);
 44:   }
 45:   PETSC_COMM_SELF, PetscFree(ck);
 46:   return r;
 47: }

 49: int main(int argc, char **argv)
 50: {
 51:   PetscInt      p, dim, Np;
 52:   PetscScalar  *randVecNums;
 53:   PetscReal     speed, value, *x, *v;
 54:   PetscRandom   rngx, rng1, rng2;
 55:   Vec           randVec, subvecvx, subvecvy;
 56:   IS            isvx, isvy;
 57:   AppCtx        user;
 58:   PetscDrawAxis axis;
 59:   PetscDraw     positionDraw;
 60:   PetscDrawSP   positionDrawSP;
 61:   MPI_Comm      comm;

 64:   PetscInitialize(&argc, &argv, NULL, help);
 65:   comm = PETSC_COMM_WORLD;
 66:   ProcessOptions(comm, &user);

 68:   Np  = user.Np;
 69:   dim = user.dim;

 71:   PetscMalloc2(Np * dim, &x, Np * dim, &v);

 73:   PetscRandomCreate(comm, &rngx);
 74:   PetscRandomSetInterval(rngx, 0., 1.);
 75:   PetscRandomSetFromOptions(rngx);
 76:   PetscRandomSetSeed(rngx, 1034);
 77:   PetscRandomSeed(rngx);

 79:   PetscRandomCreate(comm, &rng1);
 80:   PetscRandomSetInterval(rng1, 0., 1.);
 81:   PetscRandomSetFromOptions(rng1);
 82:   PetscRandomSetSeed(rng1, 3084);
 83:   PetscRandomSeed(rng1);

 85:   PetscRandomCreate(comm, &rng2);
 86:   PetscRandomSetInterval(rng2, 0., 1.);
 87:   PetscRandomSetFromOptions(rng2);
 88:   PetscRandomSetSeed(rng2, 2397);
 89:   PetscRandomSeed(rng2);

 91:   /* Set particle positions and velocities */
 92:   if (user.dim_inp == 1) {
 93:     for (p = 0; p < Np; ++p) {
 94:       PetscReal temp;
 95:       PetscRandomGetValueReal(rngx, &value);
 96:       x[p * dim]     = value;
 97:       x[p * dim + 1] = 0.;
 98:       temp           = erfinv(2 * value - 1);
 99:       v[p * dim]     = temp;
100:       v[p * dim + 1] = 0.;
101:     }
102:   } else if (user.dim_inp == 2) {
103:     /*
104:      Use Box-Muller to sample a distribution of velocities for the maxwellian.
105:      https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
106:     */
107:     VecCreate(comm, &randVec);
108:     VecSetSizes(randVec, PETSC_DECIDE, Np * dim);
109:     VecSetFromOptions(randVec);

111:     ISCreateStride(comm, Np * dim / 2, 0, 2, &isvx);
112:     ISCreateStride(comm, Np * dim / 2, 1, 2, &isvy);
113:     VecGetSubVector(randVec, isvx, &subvecvx);
114:     VecGetSubVector(randVec, isvy, &subvecvy);
115:     VecSetRandom(subvecvx, rng1);
116:     VecSetRandom(subvecvy, rng2);
117:     VecRestoreSubVector(randVec, isvx, &subvecvx);
118:     VecRestoreSubVector(randVec, isvy, &subvecvy);
119:     VecGetArray(randVec, &randVecNums);

121:     for (p = 0; p < Np; ++p) {
122:       PetscReal u1, u2, mag, zx, zy;

124:       u1 = PetscRealPart(randVecNums[p * dim]);
125:       u2 = PetscRealPart(randVecNums[p * dim + 1]);

127:       x[p * dim]     = u1;
128:       x[p * dim + 1] = u2;

130:       mag = PetscSqrtReal(-2.0 * PetscLogReal(u1));

132:       zx = mag * PetscCosReal(2 * PETSC_PI * u2) + 0;
133:       zy = mag * PetscSinReal(2 * PETSC_PI * u2) + 0;

135:       v[p * dim]     = zx;
136:       v[p * dim + 1] = zy;
137:     }
138:     ISDestroy(&isvx);
139:     ISDestroy(&isvy);
140:     VecDestroy(&subvecvx);
141:     VecDestroy(&subvecvy);
142:     VecDestroy(&randVec);
143:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Do not support dimension %" PetscInt_FMT, dim);

145:   PetscDrawCreate(comm, NULL, "monitor_particle_positions", 0, 0, 400, 300, &positionDraw);
146:   PetscDrawSetFromOptions(positionDraw);
147:   PetscDrawSPCreate(positionDraw, 10, &positionDrawSP);
148:   PetscDrawSPSetDimension(positionDrawSP, 1);
149:   PetscDrawSPGetAxis(positionDrawSP, &axis);
150:   PetscDrawSPReset(positionDrawSP);
151:   PetscDrawAxisSetLabels(axis, "Particles", "x", "y");
152:   PetscDrawSetSave(positionDraw, "ex35_pos.ppm");
153:   PetscDrawSPReset(positionDrawSP);
154:   PetscDrawSPSetLimits(positionDrawSP, 0, 1, 0, 1);
155:   for (p = 0; p < Np; ++p) {
156:     speed = PetscSqrtReal(PetscSqr(v[p * dim]) + PetscSqr(v[p * dim + 1]));
157:     PetscDrawSPAddPointColorized(positionDrawSP, &x[p * dim], &x[p * dim + 1], &speed);
158:   }
159:   PetscDrawSPDraw(positionDrawSP, PETSC_TRUE);
160:   PetscDrawSave(positionDraw);

162:   PetscFree2(x, v);
163:   PetscRandomDestroy(&rngx);
164:   PetscRandomDestroy(&rng1);
165:   PetscRandomDestroy(&rng2);

167:   PetscDrawSPDestroy(&positionDrawSP);
168:   PetscDrawDestroy(&positionDraw);
169:   PetscFinalize();
170:   return 0;
171: }

173: /*TEST
174:    test:
175:      suffix: 1D
176:      args: -Np 50\
177:      -dim 1
178:    test:
179:      suffix: 2D
180:      args: -Np 50\
181:      -dim 2
182: TEST*/