Changes between Initial Version and Version 1 of cypress/Programming/ParticleStrengthExchangeMethod/ex62.c

Timestamp:

05/14/15 14:46:44 (10 years ago)

Author:

cmaggio

Comment:

—

cypress/Programming/ParticleStrengthExchangeMethod/ex62.c

@@ +1 @@
+= 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);
+}
+}}}