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 PetscErrorCode ierr; 16 17 PetscFunctionBeginUser; 18 options->dim = 2; 19 options->dim_inp = 2; 20 options->Np = 100; 21 22 ierr = PetscOptionsBegin(comm, "", "Test of colorized scatter plot", "");PetscCall(ierr); 23 PetscCall(PetscOptionsInt("-Np", "Number of particles", "ex35.c", options->Np, &options->Np, PETSC_NULL)); 24 PetscCall(PetscOptionsInt("-dim", "Number of dimensions", "ex35.c", options->dim_inp, &options->dim_inp, PETSC_NULL)); 25 ierr = PetscOptionsEnd();PetscCall(ierr); 26 PetscFunctionReturn(0); 27 } 28 29 /* 30 ref: http://www.mimirgames.com/articles/programming/approximations-of-the-inverse-error-function/ 31 */ 32 PetscReal erfinv(PetscReal x) 33 { 34 PetscReal *ck, r = 0.; 35 PetscInt maxIter = 100; 36 37 PetscCall(PetscCalloc1(maxIter,&ck)); 38 ck[0] = 1; 39 r = ck[0]*((PetscSqrtReal(PETSC_PI)/2.)*x); 40 for (PetscInt k = 1; k < maxIter; ++k){ 41 for (PetscInt m = 0; m <= k-1; ++m){ 42 PetscReal denom = (m+1.)*(2.*m+1.); 43 ck[k] += (ck[m]*ck[k-1-m])/denom; 44 } 45 PetscReal temp = 2.*k+1.; 46 r += (ck[k]/temp)*PetscPowReal((PetscSqrtReal(PETSC_PI)/2.)*x,2.*k+1.); 47 } 48 PetscCallAbort(PETSC_COMM_SELF,PetscFree(ck)); 49 return r; 50 } 51 52 int main(int argc, char **argv) 53 { 54 PetscInt p, dim, Np; 55 PetscScalar *randVecNums; 56 PetscReal speed, value, *x, *v; 57 PetscRandom rngx, rng1, rng2; 58 Vec randVec, subvecvx, subvecvy; 59 IS isvx, isvy; 60 AppCtx user; 61 PetscDrawAxis axis; 62 PetscDraw positionDraw; 63 PetscDrawSP positionDrawSP; 64 MPI_Comm comm; 65 66 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 67 comm = PETSC_COMM_WORLD; 68 PetscCall(ProcessOptions(comm, &user)); 69 70 Np = user.Np; 71 dim = user.dim; 72 73 PetscCall(PetscMalloc2(Np*dim, &x, Np*dim, &v)); 74 75 PetscCall(PetscRandomCreate(comm, &rngx)); 76 PetscCall(PetscRandomSetInterval(rngx, 0., 1.)); 77 PetscCall(PetscRandomSetFromOptions(rngx)); 78 PetscCall(PetscRandomSetSeed(rngx, 1034)); 79 PetscCall(PetscRandomSeed(rngx)); 80 81 PetscCall(PetscRandomCreate(comm, &rng1)); 82 PetscCall(PetscRandomSetInterval(rng1, 0., 1.)); 83 PetscCall(PetscRandomSetFromOptions(rng1)); 84 PetscCall(PetscRandomSetSeed(rng1, 3084)); 85 PetscCall(PetscRandomSeed(rng1)); 86 87 PetscCall(PetscRandomCreate(comm, &rng2)); 88 PetscCall(PetscRandomSetInterval(rng2, 0., 1.)); 89 PetscCall(PetscRandomSetFromOptions(rng2)); 90 PetscCall(PetscRandomSetSeed(rng2, 2397)); 91 PetscCall(PetscRandomSeed(rng2)); 92 93 /* Set particle positions and velocities */ 94 if (user.dim_inp == 1) { 95 for (p = 0; p < Np; ++p) { 96 PetscReal temp; 97 PetscCall(PetscRandomGetValueReal(rngx, &value)); 98 x[p*dim] = value; 99 x[p*dim+1] = 0.; 100 temp = erfinv(2*value-1); 101 v[p*dim] = temp; 102 v[p*dim+1] = 0.; 103 } 104 } else if (user.dim_inp == 2) { 105 /* 106 Use Box-Muller to sample a distribution of velocities for the maxwellian. 107 https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform 108 */ 109 PetscCall(VecCreate(comm,&randVec)); 110 PetscCall(VecSetSizes(randVec,PETSC_DECIDE, Np*dim)); 111 PetscCall(VecSetFromOptions(randVec)); 112 113 PetscCall(ISCreateStride(comm, Np*dim/2, 0, 2, &isvx)); 114 PetscCall(ISCreateStride(comm, Np*dim/2, 1, 2, &isvy)); 115 PetscCall(VecGetSubVector(randVec, isvx, &subvecvx)); 116 PetscCall(VecGetSubVector(randVec, isvy, &subvecvy)); 117 PetscCall(VecSetRandom(subvecvx, rng1)); 118 PetscCall(VecSetRandom(subvecvy, rng2)); 119 PetscCall(VecRestoreSubVector(randVec, isvx, &subvecvx)); 120 PetscCall(VecRestoreSubVector(randVec, isvy, &subvecvy)); 121 PetscCall(VecGetArray(randVec, &randVecNums)); 122 123 for (p = 0; p < Np; ++p) { 124 PetscReal u1, u2, mag, zx, zy; 125 126 u1 = PetscRealPart(randVecNums[p*dim]); 127 u2 = PetscRealPart(randVecNums[p*dim+1]); 128 129 x[p*dim] = u1; 130 x[p*dim+1] = u2; 131 132 mag = PetscSqrtReal(-2.0 * PetscLogReal(u1)); 133 134 zx = mag * PetscCosReal(2*PETSC_PI*u2) + 0; 135 zy = mag * PetscSinReal(2*PETSC_PI*u2) + 0; 136 137 v[p*dim] = zx; 138 v[p*dim+1] = zy; 139 } 140 PetscCall(ISDestroy(&isvx)); 141 PetscCall(ISDestroy(&isvy)); 142 PetscCall(VecDestroy(&subvecvx)); 143 PetscCall(VecDestroy(&subvecvy)); 144 PetscCall(VecDestroy(&randVec)); 145 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Do not support dimension %D", dim); 146 147 PetscCall(PetscDrawCreate(comm, NULL, "monitor_particle_positions", 0,0,400,300, &positionDraw)); 148 PetscCall(PetscDrawSetFromOptions(positionDraw)); 149 PetscCall(PetscDrawSPCreate(positionDraw, 10, &positionDrawSP)); 150 PetscCall(PetscDrawSPSetDimension(positionDrawSP,1)); 151 PetscCall(PetscDrawSPGetAxis(positionDrawSP, &axis)); 152 PetscCall(PetscDrawSPReset(positionDrawSP)); 153 PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "x", "y")); 154 PetscCall(PetscDrawSetSave(positionDraw, "ex35_pos.ppm")); 155 PetscCall(PetscDrawSPReset(positionDrawSP)); 156 PetscCall(PetscDrawSPSetLimits(positionDrawSP, 0, 1, 0, 1)); 157 for (p = 0; p < Np; ++p) { 158 speed = PetscSqrtReal(PetscSqr(v[p*dim]) + PetscSqr(v[p*dim+1])); 159 PetscCall(PetscDrawSPAddPointColorized(positionDrawSP, &x[p*dim], &x[p*dim+1], &speed)); 160 } 161 PetscCall(PetscDrawSPDraw(positionDrawSP, PETSC_TRUE)); 162 PetscCall(PetscDrawSave(positionDraw)); 163 164 PetscCall(PetscFree2(x, v)); 165 PetscCall(PetscRandomDestroy(&rngx)); 166 PetscCall(PetscRandomDestroy(&rng1)); 167 PetscCall(PetscRandomDestroy(&rng2)); 168 169 PetscCall(PetscDrawSPDestroy(&positionDrawSP)); 170 PetscCall(PetscDrawDestroy(&positionDraw)); 171 PetscCall(PetscFinalize()); 172 return 0; 173 } 174 175 /*TEST 176 test: 177 suffix: 1D 178 args: -Np 50\ 179 -dim 1 180 test: 181 suffix: 2D 182 args: -Np 50\ 183 -dim 2 184 TEST*/ 185