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: {
15: PetscFunctionBeginUser;
16: options->dim = 2;
17: options->dim_inp = 2;
18: options->Np = 100;
19: PetscOptionsBegin(comm, "", "Test of colorized scatter plot", "");
20: PetscCall(PetscOptionsInt("-Np", "Number of particles", "ex35.c", options->Np, &options->Np, NULL));
21: PetscCall(PetscOptionsInt("-dim", "Number of dimensions", "ex35.c", options->dim_inp, &options->dim_inp, NULL));
22: PetscOptionsEnd();
23: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscCall(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: PetscCallAbort(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;
63: PetscFunctionBeginUser;
64: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
65: comm = PETSC_COMM_WORLD;
66: PetscCall(ProcessOptions(comm, &user));
68: Np = user.Np;
69: dim = user.dim;
71: PetscCall(PetscMalloc2(Np * dim, &x, Np * dim, &v));
73: PetscCall(PetscRandomCreate(comm, &rngx));
74: PetscCall(PetscRandomSetInterval(rngx, 0., 1.));
75: PetscCall(PetscRandomSetFromOptions(rngx));
76: PetscCall(PetscRandomSetSeed(rngx, 1034));
77: PetscCall(PetscRandomSeed(rngx));
79: PetscCall(PetscRandomCreate(comm, &rng1));
80: PetscCall(PetscRandomSetInterval(rng1, 0., 1.));
81: PetscCall(PetscRandomSetFromOptions(rng1));
82: PetscCall(PetscRandomSetSeed(rng1, 3084));
83: PetscCall(PetscRandomSeed(rng1));
85: PetscCall(PetscRandomCreate(comm, &rng2));
86: PetscCall(PetscRandomSetInterval(rng2, 0., 1.));
87: PetscCall(PetscRandomSetFromOptions(rng2));
88: PetscCall(PetscRandomSetSeed(rng2, 2397));
89: PetscCall(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: PetscCall(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: PetscCall(VecCreate(comm, &randVec));
108: PetscCall(VecSetSizes(randVec, PETSC_DECIDE, Np * dim));
109: PetscCall(VecSetFromOptions(randVec));
111: PetscCall(ISCreateStride(comm, Np * dim / 2, 0, 2, &isvx));
112: PetscCall(ISCreateStride(comm, Np * dim / 2, 1, 2, &isvy));
113: PetscCall(VecGetSubVector(randVec, isvx, &subvecvx));
114: PetscCall(VecGetSubVector(randVec, isvy, &subvecvy));
115: PetscCall(VecSetRandom(subvecvx, rng1));
116: PetscCall(VecSetRandom(subvecvy, rng2));
117: PetscCall(VecRestoreSubVector(randVec, isvx, &subvecvx));
118: PetscCall(VecRestoreSubVector(randVec, isvy, &subvecvy));
119: PetscCall(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: PetscCall(ISDestroy(&isvx));
139: PetscCall(ISDestroy(&isvy));
140: PetscCall(VecDestroy(&subvecvx));
141: PetscCall(VecDestroy(&subvecvy));
142: PetscCall(VecDestroy(&randVec));
143: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Do not support dimension %" PetscInt_FMT, dim);
145: PetscCall(PetscDrawCreate(comm, NULL, "monitor_particle_positions", 0, 0, 400, 300, &positionDraw));
146: PetscCall(PetscDrawSetFromOptions(positionDraw));
147: PetscCall(PetscDrawSPCreate(positionDraw, 10, &positionDrawSP));
148: PetscCall(PetscDrawSPSetDimension(positionDrawSP, 1));
149: PetscCall(PetscDrawSPGetAxis(positionDrawSP, &axis));
150: PetscCall(PetscDrawSPReset(positionDrawSP));
151: PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "x", "y"));
152: PetscCall(PetscDrawSetSave(positionDraw, "ex35_pos.ppm"));
153: PetscCall(PetscDrawSPReset(positionDrawSP));
154: PetscCall(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: PetscCall(PetscDrawSPAddPointColorized(positionDrawSP, &x[p * dim], &x[p * dim + 1], &speed));
158: }
159: PetscCall(PetscDrawSPDraw(positionDrawSP, PETSC_TRUE));
160: PetscCall(PetscDrawSave(positionDraw));
162: PetscCall(PetscFree2(x, v));
163: PetscCall(PetscRandomDestroy(&rngx));
164: PetscCall(PetscRandomDestroy(&rng1));
165: PetscCall(PetscRandomDestroy(&rng2));
167: PetscCall(PetscDrawSPDestroy(&positionDrawSP));
168: PetscCall(PetscDrawDestroy(&positionDraw));
169: PetscCall(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*/