/*
Large Scale Computing
Particle Strength Exchange (PSE) method
Sovle Convective-Diffusion Equation in Stagnation Point Flow
ex62.c
*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include<string.h>
#define NUMPX 80 /* number of particles in x-direction */
#define BP 2 /* number of particle lines outside of the bottom */
#define MAXTIME 2.0 /* Max Time */
#define INTVTIME 0.1 /* Interval to export data */
#define PE 10.0
#define DELTA (1.0 / (double)(NUMPX - 1)) /* Blob Size*/
#define DT (DELTA * DELTA * PE * 0.1) /* a time step */
#define CONST 1
#define NORMAL 0
/* Diffusion Kernel */
double DKernel(double r2){
double Delta2;
Delta2 = DELTA * DELTA;
return 4.0 * exp(-(r2 / Delta2)) / (M_PI * Delta2);
}
/* Compute Diffusion */
void Diffusion(int tnump,int *pt_flag,
double *pt_x,double *pt_y,
double *pt_u,double *pt_v,
double *pt_u0,double *pt_v0,
double *pt_cnc,double *pt_lap){
double r2,m,Delta2;
int i,j;
Delta2 = DELTA * DELTA;
for (i = 0 ; i < tnump ; i++){
pt_lap[i] = 0.0;
}
/* Compute Laplacian */
for (i = 0 ; i < tnump ; i++){
for (j = 0 ; j < tnump ; j++){
r2 = (pt_x[i] - pt_x[j]) * (pt_x[i] - pt_x[j])
+ (pt_y[i] - pt_y[j]) * (pt_y[i] - pt_y[j]);
m = DKernel(r2);
pt_lap[i] += (pt_cnc[j] - pt_cnc[i]) * m;
}
}
/* update C */
for (i = 0 ; i < tnump ; i++){
if (pt_flag[i] == NORMAL){
pt_cnc[i] += pt_lap[i] * DT / PE;
}
}
}
void ComputeVelocity(int tnump,int *pt_flag,
double *pt_x,double *pt_y,
double *pt_u,double *pt_v,
double *pt_u0,double *pt_v0,
double *pt_cnc,double *pt_lap){
int i;
/* Stagnation Point Flow */
for (i = 0 ; i < tnump ; i++){
pt_u[i] = pt_x[i];
pt_v[i] = -pt_y[i];
}
}
void Convection(int tnump,int *pt_flag,
double *pt_x,double *pt_y,
double *pt_u,double *pt_v,
double *pt_u0,double *pt_v0,
double *pt_cnc,double *pt_lap){
int i;
/* 2nd-Order Adams-Bashforth Scheme */
for (i = 0 ; i < tnump ; i++){
if (pt_flag[i] == NORMAL){
pt_x[i] += (1.5 * pt_u[i] - 0.5 * pt_u0[i]) * DT;
pt_y[i] += (1.5 * pt_v[i] - 0.5 * pt_v0[i]) * DT;
pt_u0[i] = pt_u[i];
pt_v0[i] = pt_v[i];
}
}
}
void CheckBounds(int tnump,int *pt_flag,
double *pt_x,double *pt_y,
double *pt_u,double *pt_v,
double *pt_u0,double *pt_v0,
double *pt_cnc,double *pt_lap){
int i;
double ox,oy;
/* Replace Particles out of right side onto top*/
for (i = 0 ; i < tnump ; i++){
if (pt_x[i] > 1.0){
ox = pt_x[i];
oy = pt_y[i];
pt_x[i] = oy;
pt_y[i] = 2.0 - ox;
pt_u0[i] = pt_v[i];
pt_v0[i] = -pt_u[i];
pt_cnc[i] = 1.0;
}
}
}
void DataOut(int n,int tnump,int *pt_flag,
double *pt_x,double *pt_y,
double *pt_u,double *pt_v,
double *pt_u0,double *pt_v0,
double *pt_cnc,double *pt_lap){
char fn[256];
FILE *fp;
int i;
sprintf(fn,"res%05d.dat",n);
fp = fopen(fn,"w");
for (i = 0 ; i < tnump ; i++){
fprintf(fp,"%e %e %e\n",pt_x[i],pt_y[i],pt_cnc[i]);
}
fclose(fp);
}
/* for timing */
double t_cost(void ){
static clock_t ptm;
static int fl = 0;
clock_t tm;
double sec;
tm = clock();
if (fl == 0) ptm = tm;
fl = 1;
sec = (double)(tm - ptm) / (double)CLOCKS_PER_SEC;
ptm = tm;
return sec;
}
int main(int arc,char *argv[]){
double *pt_x,*pt_y;
double *pt_u,*pt_v;
double *pt_u0,*pt_v0;
double *pt_cnc;
double *pt_lap;
int *pt_flag;
int tnump;
int i,j,k;
int tcount;
double time;
double ctime;
/* total number of particles */
tnump = NUMPX * (NUMPX + BP);
printf("Number of particle=%d\n",tnump);
printf("Time Step=%e\n",DT);
printf("Max Steps=%e\n",MAXTIME / DT);
/* allocate space */
pt_x = (double *)calloc(tnump,sizeof(double));
pt_y = (double *)calloc(tnump,sizeof(double));
pt_u = (double *)calloc(tnump,sizeof(double));
pt_v = (double *)calloc(tnump,sizeof(double));
pt_u0 = (double *)calloc(tnump,sizeof(double));
pt_v0 = (double *)calloc(tnump,sizeof(double));
pt_cnc = (double *)calloc(tnump,sizeof(double));
pt_lap = (double *)calloc(tnump,sizeof(double));
pt_flag = (int *)calloc(tnump,sizeof(int));
/* Initialize */
for (i = 0 ; i < tnump ; i++){
pt_x[i] = (double)rand() / (double)RAND_MAX;
pt_y[i] = -2.0 * DELTA + (1.0 + 2.0 * DELTA) * (double)rand() / (double)RAND_MAX;
pt_u0[i] = pt_u[i] = pt_x[i];
pt_v0[i] = pt_v[i] = -pt_y[i];
pt_cnc[i] = 0.0;
pt_flag[i] = NORMAL;
if (pt_y[i] <= 0.0){
pt_flag[i] = CONST;
}
}
/* compute */
t_cost();
tcount = 0;
time = 0.0;
while(1){
/* Compute Diffusion */
Diffusion(tnump,pt_flag,pt_x,pt_y,pt_u,pt_v,pt_u0,pt_v0,pt_cnc,pt_lap);
/* Compute Velocity */
ComputeVelocity(tnump,pt_flag,pt_x,pt_y,pt_u,pt_v,pt_u0,pt_v0,pt_cnc,pt_lap);
/* Update Particle Location */
Convection(tnump,pt_flag,pt_x,pt_y,pt_u,pt_v,pt_u0,pt_v0,pt_cnc,pt_lap);
/* Check Bounds */
CheckBounds(tnump,pt_flag,pt_x,pt_y,pt_u,pt_v,pt_u0,pt_v0,pt_cnc,pt_lap);
ctime = t_cost();
/* increment time */
time += DT;
/* Data Out */
if (time >= INTVTIME * (double)tcount){
DataOut(tcount,tnump,pt_flag,pt_x,pt_y,pt_u,pt_v,pt_u0,pt_v0,pt_cnc,pt_lap);
printf("Data Out at time=%e, %e (sec/step)\n",time,ctime);
tcount++;
}
/*
break;
*/
if (time >= MAXTIME) break;
}
/* Cleanup */
free(pt_x);
free(pt_y);
free(pt_u);
free(pt_v);
free(pt_u0);
free(pt_v0);
free(pt_cnc);
free(pt_lap);
free(pt_flag);
}