| 1 | = Ex62.c = |
| 2 | {{{#!C |
| 3 | /* |
| 4 | Large Scale Computing |
| 5 | Particle Strength Exchange (PSE) method |
| 6 | Sovle Convective-Diffusion Equation in Stagnation Point Flow |
| 7 | |
| 8 | ex62.c |
| 9 | */ |
| 10 | |
| 11 | #include<stdio.h> |
| 12 | #include<stdlib.h> |
| 13 | #include<math.h> |
| 14 | #include<time.h> |
| 15 | #include<string.h> |
| 16 | |
| 17 | #define NUMPX 80 /* number of particles in x-direction */ |
| 18 | #define BP 2 /* number of particle lines outside of the bottom */ |
| 19 | #define MAXTIME 2.0 /* Max Time */ |
| 20 | #define INTVTIME 0.1 /* Interval to export data */ |
| 21 | #define PE 10.0 |
| 22 | |
| 23 | #define DELTA (1.0 / (double)(NUMPX - 1)) /* Blob Size*/ |
| 24 | #define DT (DELTA * DELTA * PE * 0.1) /* a time step */ |
| 25 | #define CONST 1 |
| 26 | #define NORMAL 0 |
| 27 | |
| 28 | /* Diffusion Kernel */ |
| 29 | double DKernel(double r2){ |
| 30 | double Delta2; |
| 31 | Delta2 = DELTA * DELTA; |
| 32 | return 4.0 * exp(-(r2 / Delta2)) / (M_PI * Delta2); |
| 33 | } |
| 34 | |
| 35 | /* Compute Diffusion */ |
| 36 | void Diffusion(int tnump,int *pt_flag, |
| 37 | double *pt_x,double *pt_y, |
| 38 | double *pt_u,double *pt_v, |
| 39 | double *pt_u0,double *pt_v0, |
| 40 | double *pt_cnc,double *pt_lap){ |
| 41 | double r2,m,Delta2; |
| 42 | int i,j; |
| 43 | |
| 44 | Delta2 = DELTA * DELTA; |
| 45 | |
| 46 | for (i = 0 ; i < tnump ; i++){ |
| 47 | pt_lap[i] = 0.0; |
| 48 | } |
| 49 | |
| 50 | /* Compute Laplacian */ |
| 51 | for (i = 0 ; i < tnump ; i++){ |
| 52 | for (j = 0 ; j < tnump ; j++){ |
| 53 | r2 = (pt_x[i] - pt_x[j]) * (pt_x[i] - pt_x[j]) |
| 54 | + (pt_y[i] - pt_y[j]) * (pt_y[i] - pt_y[j]); |
| 55 | |
| 56 | m = DKernel(r2); |
| 57 | pt_lap[i] += (pt_cnc[j] - pt_cnc[i]) * m; |
| 58 | } |
| 59 | } |
| 60 | |
| 61 | /* update C */ |
| 62 | for (i = 0 ; i < tnump ; i++){ |
| 63 | if (pt_flag[i] == NORMAL){ |
| 64 | pt_cnc[i] += pt_lap[i] * DT / PE; |
| 65 | } |
| 66 | } |
| 67 | } |
| 68 | |
| 69 | void ComputeVelocity(int tnump,int *pt_flag, |
| 70 | double *pt_x,double *pt_y, |
| 71 | double *pt_u,double *pt_v, |
| 72 | double *pt_u0,double *pt_v0, |
| 73 | double *pt_cnc,double *pt_lap){ |
| 74 | int i; |
| 75 | |
| 76 | /* Stagnation Point Flow */ |
| 77 | for (i = 0 ; i < tnump ; i++){ |
| 78 | pt_u[i] = pt_x[i]; |
| 79 | pt_v[i] = -pt_y[i]; |
| 80 | } |
| 81 | } |
| 82 | |
| 83 | void Convection(int tnump,int *pt_flag, |
| 84 | double *pt_x,double *pt_y, |
| 85 | double *pt_u,double *pt_v, |
| 86 | double *pt_u0,double *pt_v0, |
| 87 | double *pt_cnc,double *pt_lap){ |
| 88 | int i; |
| 89 | |
| 90 | /* 2nd-Order Adams-Bashforth Scheme */ |
| 91 | for (i = 0 ; i < tnump ; i++){ |
| 92 | if (pt_flag[i] == NORMAL){ |
| 93 | pt_x[i] += (1.5 * pt_u[i] - 0.5 * pt_u0[i]) * DT; |
| 94 | pt_y[i] += (1.5 * pt_v[i] - 0.5 * pt_v0[i]) * DT; |
| 95 | pt_u0[i] = pt_u[i]; |
| 96 | pt_v0[i] = pt_v[i]; |
| 97 | } |
| 98 | } |
| 99 | } |
| 100 | |
| 101 | void CheckBounds(int tnump,int *pt_flag, |
| 102 | double *pt_x,double *pt_y, |
| 103 | double *pt_u,double *pt_v, |
| 104 | double *pt_u0,double *pt_v0, |
| 105 | double *pt_cnc,double *pt_lap){ |
| 106 | int i; |
| 107 | double ox,oy; |
| 108 | |
| 109 | /* Replace Particles out of right side onto top*/ |
| 110 | for (i = 0 ; i < tnump ; i++){ |
| 111 | if (pt_x[i] > 1.0){ |
| 112 | ox = pt_x[i]; |
| 113 | oy = pt_y[i]; |
| 114 | pt_x[i] = oy; |
| 115 | pt_y[i] = 2.0 - ox; |
| 116 | pt_u0[i] = pt_v[i]; |
| 117 | pt_v0[i] = -pt_u[i]; |
| 118 | pt_cnc[i] = 1.0; |
| 119 | } |
| 120 | } |
| 121 | } |
| 122 | |
| 123 | void DataOut(int n,int tnump,int *pt_flag, |
| 124 | double *pt_x,double *pt_y, |
| 125 | double *pt_u,double *pt_v, |
| 126 | double *pt_u0,double *pt_v0, |
| 127 | double *pt_cnc,double *pt_lap){ |
| 128 | char fn[256]; |
| 129 | FILE *fp; |
| 130 | int i; |
| 131 | |
| 132 | sprintf(fn,"res%05d.dat",n); |
| 133 | fp = fopen(fn,"w"); |
| 134 | for (i = 0 ; i < tnump ; i++){ |
| 135 | fprintf(fp,"%e %e %e\n",pt_x[i],pt_y[i],pt_cnc[i]); |
| 136 | } |
| 137 | fclose(fp); |
| 138 | } |
| 139 | |
| 140 | /* for timing */ |
| 141 | double t_cost(void ){ |
| 142 | static clock_t ptm; |
| 143 | static int fl = 0; |
| 144 | clock_t tm; |
| 145 | double sec; |
| 146 | |
| 147 | tm = clock(); |
| 148 | if (fl == 0) ptm = tm; |
| 149 | fl = 1; |
| 150 | |
| 151 | sec = (double)(tm - ptm) / (double)CLOCKS_PER_SEC; |
| 152 | |
| 153 | ptm = tm; |
| 154 | |
| 155 | return sec; |
| 156 | } |
| 157 | |
| 158 | int main(int arc,char *argv[]){ |
| 159 | double *pt_x,*pt_y; |
| 160 | double *pt_u,*pt_v; |
| 161 | double *pt_u0,*pt_v0; |
| 162 | double *pt_cnc; |
| 163 | double *pt_lap; |
| 164 | int *pt_flag; |
| 165 | int tnump; |
| 166 | int i,j,k; |
| 167 | int tcount; |
| 168 | double time; |
| 169 | double ctime; |
| 170 | |
| 171 | /* total number of particles */ |
| 172 | tnump = NUMPX * (NUMPX + BP); |
| 173 | printf("Number of particle=%d\n",tnump); |
| 174 | printf("Time Step=%e\n",DT); |
| 175 | printf("Max Steps=%e\n",MAXTIME / DT); |
| 176 | |
| 177 | /* allocate space */ |
| 178 | pt_x = (double *)calloc(tnump,sizeof(double)); |
| 179 | pt_y = (double *)calloc(tnump,sizeof(double)); |
| 180 | pt_u = (double *)calloc(tnump,sizeof(double)); |
| 181 | pt_v = (double *)calloc(tnump,sizeof(double)); |
| 182 | pt_u0 = (double *)calloc(tnump,sizeof(double)); |
| 183 | pt_v0 = (double *)calloc(tnump,sizeof(double)); |
| 184 | pt_cnc = (double *)calloc(tnump,sizeof(double)); |
| 185 | pt_lap = (double *)calloc(tnump,sizeof(double)); |
| 186 | pt_flag = (int *)calloc(tnump,sizeof(int)); |
| 187 | |
| 188 | /* Initialize */ |
| 189 | for (i = 0 ; i < tnump ; i++){ |
| 190 | pt_x[i] = (double)rand() / (double)RAND_MAX; |
| 191 | pt_y[i] = -2.0 * DELTA + (1.0 + 2.0 * DELTA) * (double)rand() / (double)RAND_MAX; |
| 192 | pt_u0[i] = pt_u[i] = pt_x[i]; |
| 193 | pt_v0[i] = pt_v[i] = -pt_y[i]; |
| 194 | pt_cnc[i] = 0.0; |
| 195 | pt_flag[i] = NORMAL; |
| 196 | if (pt_y[i] <= 0.0){ |
| 197 | pt_flag[i] = CONST; |
| 198 | } |
| 199 | } |
| 200 | |
| 201 | /* compute */ |
| 202 | t_cost(); |
| 203 | tcount = 0; |
| 204 | time = 0.0; |
| 205 | while(1){ |
| 206 | /* Compute Diffusion */ |
| 207 | Diffusion(tnump,pt_flag,pt_x,pt_y,pt_u,pt_v,pt_u0,pt_v0,pt_cnc,pt_lap); |
| 208 | |
| 209 | /* Compute Velocity */ |
| 210 | ComputeVelocity(tnump,pt_flag,pt_x,pt_y,pt_u,pt_v,pt_u0,pt_v0,pt_cnc,pt_lap); |
| 211 | |
| 212 | /* Update Particle Location */ |
| 213 | Convection(tnump,pt_flag,pt_x,pt_y,pt_u,pt_v,pt_u0,pt_v0,pt_cnc,pt_lap); |
| 214 | |
| 215 | /* Check Bounds */ |
| 216 | CheckBounds(tnump,pt_flag,pt_x,pt_y,pt_u,pt_v,pt_u0,pt_v0,pt_cnc,pt_lap); |
| 217 | |
| 218 | ctime = t_cost(); |
| 219 | |
| 220 | /* increment time */ |
| 221 | time += DT; |
| 222 | |
| 223 | /* Data Out */ |
| 224 | if (time >= INTVTIME * (double)tcount){ |
| 225 | DataOut(tcount,tnump,pt_flag,pt_x,pt_y,pt_u,pt_v,pt_u0,pt_v0,pt_cnc,pt_lap); |
| 226 | printf("Data Out at time=%e, %e (sec/step)\n",time,ctime); |
| 227 | tcount++; |
| 228 | } |
| 229 | /* |
| 230 | break; |
| 231 | */ |
| 232 | if (time >= MAXTIME) break; |
| 233 | } |
| 234 | |
| 235 | /* Cleanup */ |
| 236 | free(pt_x); |
| 237 | free(pt_y); |
| 238 | free(pt_u); |
| 239 | free(pt_v); |
| 240 | free(pt_u0); |
| 241 | free(pt_v0); |
| 242 | free(pt_cnc); |
| 243 | free(pt_lap); |
| 244 | free(pt_flag); |
| 245 | } |
| 246 | }}} |