= Ex62.c = 
{{{#!C
/*
  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);
}
}}}