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 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 */ 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 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 test: 179 suffix: 2D 180 args: -Np 50\ 181 -dim 2 182 TEST*/ 183