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


Ignore:
Timestamp:
May 14, 2015 2:46:44 PM (7 years ago)
Author:
cmaggio
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • cypress/Programming/ParticleStrengthExchangeMethod/ex62.c

    v1 v1  
     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 */
     29double 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 */
     36void 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 
     69void 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 
     83void 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 
     101void 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 
     123void 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 */
     141double 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 
     158int 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}}}