1 static char help[] = "Test of Colorized Scatter Plot.\n";
2
3 #include <petscdraw.h>
4 #include <petscvec.h>
5 #include <petscis.h>
6
7 typedef struct {
8 PetscInt Np; /* total number of particles */
9 PetscInt dim;
10 PetscInt dim_inp;
11 } AppCtx;
12
ProcessOptions(MPI_Comm comm,AppCtx * options)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 }
25
26 /*
27 ref: http://www.mimirgames.com/articles/programming/approximations-of-the-inverse-error-function/
28 */
erfinv(PetscReal x)29 PetscReal erfinv(PetscReal x)
30 {
31 PetscReal *ck, r = 0.;
32 PetscInt maxIter = 100;
33
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 }
48
main(int argc,char ** argv)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;
62
63 PetscFunctionBeginUser;
64 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
65 comm = PETSC_COMM_WORLD;
66 PetscCall(ProcessOptions(comm, &user));
67
68 Np = user.Np;
69 dim = user.dim;
70
71 PetscCall(PetscMalloc2(Np * dim, &x, Np * dim, &v));
72
73 PetscCall(PetscRandomCreate(comm, &rngx));
74 PetscCall(PetscRandomSetInterval(rngx, 0., 1.));
75 PetscCall(PetscRandomSetFromOptions(rngx));
76 PetscCall(PetscRandomSetSeed(rngx, 1034));
77 PetscCall(PetscRandomSeed(rngx));
78
79 PetscCall(PetscRandomCreate(comm, &rng1));
80 PetscCall(PetscRandomSetInterval(rng1, 0., 1.));
81 PetscCall(PetscRandomSetFromOptions(rng1));
82 PetscCall(PetscRandomSetSeed(rng1, 3084));
83 PetscCall(PetscRandomSeed(rng1));
84
85 PetscCall(PetscRandomCreate(comm, &rng2));
86 PetscCall(PetscRandomSetInterval(rng2, 0., 1.));
87 PetscCall(PetscRandomSetFromOptions(rng2));
88 PetscCall(PetscRandomSetSeed(rng2, 2397));
89 PetscCall(PetscRandomSeed(rng2));
90
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));
110
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));
120
121 for (p = 0; p < Np; ++p) {
122 PetscReal u1, u2, mag, zx, zy;
123
124 u1 = PetscRealPart(randVecNums[p * dim]);
125 u2 = PetscRealPart(randVecNums[p * dim + 1]);
126
127 x[p * dim] = u1;
128 x[p * dim + 1] = u2;
129
130 mag = PetscSqrtReal(-2.0 * PetscLogReal(u1));
131
132 zx = mag * PetscCosReal(2 * PETSC_PI * u2) + 0;
133 zy = mag * PetscSinReal(2 * PETSC_PI * u2) + 0;
134
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);
144
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));
161
162 PetscCall(PetscFree2(x, v));
163 PetscCall(PetscRandomDestroy(&rngx));
164 PetscCall(PetscRandomDestroy(&rng1));
165 PetscCall(PetscRandomDestroy(&rng2));
166
167 PetscCall(PetscDrawSPDestroy(&positionDrawSP));
168 PetscCall(PetscDrawDestroy(&positionDraw));
169 PetscCall(PetscFinalize());
170 return 0;
171 }
172
173 /*TEST
174 test:
175 suffix: 1D
176 args: -Np 50\
177 -dim 1
178 output_file: output/empty.out
179 test:
180 suffix: 2D
181 args: -Np 50\
182 -dim 2
183 output_file: output/empty.out
184 TEST*/
185