Changes between Initial Version and Version 1 of cypress/Programming/StokesFlow/ex43.c


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

migrated content from ccs wiki

Legend:

Unmodified
Added
Removed
Modified
  • cypress/Programming/StokesFlow/ex43.c

    v1 v1  
     1= Ex43.c =
     2{{{#!C
     3/*
     4  Large Scale Computing
     5  Stokes Flow in a Cavity
     6  MPI version
     7  ex43.c
     8*/
     9#include<stdio.h>
     10#include<stdlib.h>
     11#include<math.h>
     12#include<time.h>
     13#include"mpi.h"
     14#include"stokeslet2d_dist.h"
     15 
     16#define EPSILON 0.005 /* blob size*/
     17#define INTGRID 51    /* # of x-grid lines for internal velocity */
     18 
     19/* min and max macros */
     20#define max(a,b)        (((a) > (b)) ? (a) : (b))
     21#define min(a,b)        (((a) < (b)) ? (a) : (b))
     22 
     23double t_cost(void );
     24int main(int , char **);
     25int main(int argc, char **argv){
     26  double dp;
     27  int numpdepth;
     28  int numpwidth;
     29  int numparticles;
     30  double *loc; /* particle location */
     31  double *vel; /* particle velocity */
     32  double *foc; /* particle force */
     33  double *mat; /* 2Nx2N dense matrix */
     34  double cloc[2]; /* internal points */
     35  double cvel[2]; /* internal velocity */
     36  int mystart,myend,mysize;
     37  int *ostart,*oend;
     38  int myid,numproc;
     39  int i,j,li;
     40  int nyg,cpu;
     41  FILE *fp;
     42  char fn[256];
     43 
     44  /* Initialize */
     45  MPI_Init(&argc,&argv);
     46 
     47  /* get myid and # of processors */
     48  MPI_Comm_size(MPI_COMM_WORLD,&numproc);
     49  MPI_Comm_rank(MPI_COMM_WORLD,&myid);
     50 
     51  if (argc < 2){
     52    if (myid == 0) printf("Usage:%s [Depth of Cavity]\n",argv[0]);
     53    MPI_Abort(MPI_COMM_WORLD, -1);
     54  }
     55 
     56  /* get inputed depth */
     57  dp = atof(argv[1]);
     58 
     59  /* # of particles in depth */
     60  numpdepth = (int)(dp / EPSILON + 0.5);
     61 
     62  /* # of particles in width */
     63  numpwidth = (int)(1.0 / EPSILON + 0.5);
     64 
     65  /* total # od particles */
     66  numparticles = numpdepth * 2 + numpwidth * 2;
     67  if (myid == 0) printf("Total # of Particles=%d\n",numparticles);
     68 
     69  /* divide blobs into # of cpus */
     70  mystart = (numparticles / numproc) * myid;
     71  if (numparticles % numproc > myid){
     72    mystart += myid;
     73    myend = mystart + (numparticles / numproc) + 1;
     74  }else{
     75    mystart += numparticles % numproc;
     76    myend = mystart + (numparticles / numproc);
     77  }
     78  mysize = myend - mystart;
     79  printf("CPU%d %d ~ %d : size=%d\n",myid,mystart,myend,mysize);
     80 
     81  /* collect start end info */
     82  ostart = (int *)calloc(numproc, sizeof(int));
     83  oend = (int *)calloc(numproc, sizeof(int));
     84  MPI_Allgather(&mysize, 1, MPI_INT, oend, 1, MPI_INT, MPI_COMM_WORLD); 
     85  ostart[0] = 0;
     86  for (cpu = 1 ; cpu < numproc ; cpu++){
     87    ostart[cpu] = ostart[cpu-1] + oend[cpu-1];
     88    oend[cpu-1] += ostart[cpu-1];
     89  }
     90  oend[numproc-1] += ostart[numproc-1];
     91 
     92  /* Allocate Space */
     93  /* DIM=2 defined by 'stokeslet2d.h' */
     94  loc = (double *)calloc(mysize * DIM,sizeof(double));
     95  vel = (double *)calloc(mysize * DIM,sizeof(double));
     96  foc = (double *)calloc(mysize * DIM,sizeof(double));
     97  mat = (double *)calloc(mysize * DIM * numparticles * DIM,sizeof(double));
     98  if ((loc == (double *)NULL) ||
     99      (vel == (double *)NULL) ||
     100      (foc == (double *)NULL) ||
     101      (mat == (double *)NULL) ){
     102    printf("Can't allocate memory!!!\n");
     103    exit(-1);
     104  }
     105 
     106  /* Set location & velocity of particles (blob)*/
     107  for (i = mystart ; i < myend ; i++){
     108    li = i - mystart;
     109    foc[li * DIM] = 0.0;
     110    foc[li * DIM + 1] = 0.0;
     111    if ((i >= 0) && (i <= numpwidth)){ /* top */
     112      loc[li * DIM] = -0.5 + EPSILON * (double)i; /* x */
     113      loc[li * DIM + 1] = 0.0; /* y */
     114      vel[li * DIM] = 1.0;
     115      vel[li * DIM + 1] = 0.0;
     116    }else{
     117      if (i <= (numpwidth + numpdepth)){ /* right wall */     
     118        loc[li * DIM] = 0.5; /* x */
     119        loc[li * DIM + 1] = -EPSILON * (double)(i - numpwidth); /* y */
     120        vel[li * DIM] = 0.0;
     121        vel[li * DIM + 1] = 0.0;
     122      }else{
     123        if (i <= (2 * numpwidth + numpdepth)){ /* bottom */
     124          loc[li * DIM] = 0.5 - EPSILON * (double)(i - (numpwidth + numpdepth)); /* x */
     125          loc[li * DIM + 1] = -EPSILON * numpdepth; /* y */
     126          vel[li * DIM] = 0.0;
     127          vel[li * DIM + 1] = 0.0;       
     128        }else{                 /* left wall */   
     129          loc[li * DIM] = -0.5; /* x */
     130          loc[li * DIM + 1] = -EPSILON * (double)((2 * numpwidth + 2 * numpdepth) - i); /* y */
     131          vel[li * DIM] = 0.0;
     132          vel[li * DIM + 1] = 0.0;       
     133        }
     134      }
     135    }
     136  }
     137 
     138 
     139  /* make 2Nx2N Matrix */
     140  t_cost();
     141  slet2d_mkMatrix_dist(numparticles,ostart,oend,loc,EPSILON,mat);
     142  if (myid == 0) printf("setting matrix %f sec\n",t_cost());
     143 
     144  /* Sovle linear ststem */
     145  t_cost();
     146  slet2d_solve_CG_dist(numparticles,ostart,oend,mat,vel,foc);
     147  if (myid == 0) printf("Sovle linear ststem %f sec\n",t_cost());
     148 
     149  /* deallocate big matrix */
     150  free(mat);
     151 
     152  /* out particle (blob) data */
     153  sprintf(fn,"particle%02d.dat",myid);
     154  fp = fopen(fn,"w");
     155  for (i = 0 ; i < mysize ; i++){
     156    fprintf(fp,"%e %e %e %e %e %e\n",
     157            loc[i * DIM],loc[i * DIM + 1],
     158            vel[i * DIM],vel[i * DIM + 1],
     159            foc[i * DIM],foc[i * DIM + 1]);
     160  }
     161  fclose(fp);
     162 
     163  /*
     164     compute internal velocity
     165  */
     166  t_cost();
     167  if (myid == 0){
     168    fp = fopen("res.dat","w");
     169  }
     170  nyg = (int)(EPSILON * (double)(numpdepth * (INTGRID - 1)));
     171  for (j = 0 ; j < nyg ; j++){
     172    for (i = 0 ; i < INTGRID ; i++){
     173      cloc[0] = -0.5 + (double)i / (double)(INTGRID - 1);
     174      cloc[1] = -(double)j / (double)(INTGRID - 1);
     175      slet2d_velocity_dist(numparticles, ostart,oend,loc, foc, EPSILON, cloc, cvel);       
     176      if (myid == 0){
     177        fprintf(fp,"%e %e %e %e\n",cloc[0], cloc[1],cvel[0], cvel[1]);     
     178      }
     179    }
     180  }
     181  if (myid == 0){
     182    fclose(fp); 
     183  }
     184  if (myid == 0)printf("Compute internal velocity %f sec\n",t_cost());
     185 
     186 
     187  /* free */
     188  free(loc);
     189  free(vel);
     190  free(foc);
     191 
     192  printf("CPU%d is Done\n",myid);
     193 
     194  MPI_Finalize();
     195}
     196 
     197/* for timing */
     198double t_cost(void ){
     199  static clock_t ptm;
     200  static int fl = 0;
     201  clock_t tm;
     202  double sec;
     203 
     204  tm = clock();
     205  if (fl == 0) ptm = tm;
     206  fl = 1;
     207 
     208  sec = (double)(tm - ptm) / (double)CLOCKS_PER_SEC;
     209 
     210  ptm = tm;
     211 
     212  return sec;
     213}
     214stokeslet2d_dist.c
     215
     216/*
     217  Large Scale Computing
     218  Stokeslet for 2D
     219  MPI Version
     220 
     221  Ricardo Cortez,"The Method of Regularized Stokeslets",
     222  (2001), SIAM J. Sci. Comput., Vol.23, No.4, pp.1204
     223 
     224  assuming mu = 1
     225*/
     226#include<stdio.h>
     227#include<stdlib.h>
     228#include<math.h>
     229#include"mpi.h"
     230#include"stokeslet2d_dist.h"
     231 
     232/* min and max macros */
     233#define max(a,b)        (((a) > (b)) ? (a) : (b))
     234#define min(a,b)        (((a) < (b)) ? (a) : (b))
     235 
     236double term1(double ,double );
     237double term2(double ,double );
     238/* Blas Function */
     239double ddot_(int *,double *,int *,double *,int *);
     240double dnrm2_(int *, double *, int *);
     241 
     242/* make Matrix */
     243void slet2d_mkMatrix_dist(int np,int *ostart,int *oend,double *loc,double ep,double *mat){
     244  int i,j;
     245  int li,lj;
     246  double r;
     247  double dx,dy;
     248  double tr1,tr2;
     249  int mysize,maxsize;
     250  int stend[2];
     251  int myid,numproc;
     252  double *oloc,*ploc;
     253  int cpu;
     254 
     255  /* get myid and # of processors */
     256  MPI_Comm_size(MPI_COMM_WORLD,&numproc);
     257  MPI_Comm_rank(MPI_COMM_WORLD,&myid);
     258 
     259  /* get max size */
     260  maxsize = 0;
     261  for (cpu = 0 ; cpu < numproc ; cpu++){
     262    maxsize = max(oend[cpu] - ostart[cpu],maxsize);
     263  }
     264  mysize = oend[myid] - ostart[myid];
     265 
     266  /* alloc */
     267  oloc = (double *)calloc(maxsize * DIM,sizeof(double));
     268 
     269  /*zeros mat*/
     270  for (i = 0 ; i < DIM * DIM * np * mysize ; i++) mat[i] = 0.0;
     271 
     272  /* make mat */
     273  for (cpu = 0 ; cpu < numproc ; cpu++){
     274    /* get location vector from others */
     275    ploc = oloc;
     276    if (myid == cpu) ploc = loc;
     277    MPI_Bcast(ploc, DIM * (oend[cpu] - ostart[cpu]), MPI_DOUBLE, cpu, MPI_COMM_WORLD);
     278 
     279    /* compute local matrix*/
     280    for (i = ostart[cpu] ; i < oend[cpu] ; i++){   
     281      for (j = ostart[myid] ; j < oend[myid] ; j++){
     282        li = i - ostart[cpu];
     283        lj = j - ostart[myid];
     284        dx = ploc[li * DIM    ] - loc[lj * DIM    ];
     285        dy = ploc[li * DIM + 1] - loc[lj * DIM + 1];
     286        r = sqrt(dx * dx + dy * dy);
     287 
     288        tr1 = term1(r,ep) / (4.0 * M_PI);
     289        tr2 = term2(r,ep) / (4.0 * M_PI);
     290 
     291        mat[i * DIM +      lj * DIM      * np * DIM] += -tr1 + tr2 * dx * dx;
     292        mat[i * DIM +     (lj * DIM + 1) * np * DIM] +=        tr2 * dx * dy;
     293        mat[i * DIM + 1 +  lj * DIM      * np * DIM] +=        tr2 * dx * dy;
     294        mat[i * DIM + 1 + (lj * DIM + 1) * np * DIM] += -tr1 + tr2 * dy * dy;     
     295      }
     296    }
     297  }
     298 
     299  free(oloc);
     300}
     301 
     302/* Sovle Liear ststem */
     303/* CG */
     304void slet2d_solve_CG_dist(int np, int *ostart, int *oend, double *mat, double *b, double *x){
     305  double tol = 1.0e-8;
     306  double *r;
     307  double *p;
     308  double *Ap;
     309  double *work;
     310  double rr,pAp,rr1,alpha,beta;
     311  int i,j,count;
     312  double w,normb;
     313  int n;
     314  int ms,me;
     315  int cpu;
     316  int myid,numproc;
     317  int li,lj;
     318 
     319  /* get myid and # of processors */
     320  MPI_Comm_size(MPI_COMM_WORLD,&numproc);
     321  MPI_Comm_rank(MPI_COMM_WORLD,&myid);
     322 
     323  n = DIM * np;
     324  ms = DIM * ostart[myid];
     325  me = DIM * oend[myid];
     326 
     327  /* Allocate Working Spaces */
     328  r = (double *)calloc(me - ms, sizeof(double));
     329  p = (double *)calloc(me - ms, sizeof(double));
     330  Ap = (double *)calloc(me - ms, sizeof(double));
     331  work = (double *)calloc(n, sizeof(double));
     332  if ((r == NULL) ||
     333      (p == NULL) ||
     334      (Ap == NULL) ||
     335      (work == NULL)){
     336    printf("memory allocation failed\n");
     337    return ;
     338  }
     339 
     340  /* initialize */
     341  for (i = 0 ; i < me - ms ; i++){
     342    r[i] = p[i] = Ap[i] = 0.0;
     343  }
     344 
     345  /* |b| */
     346  w = 0.0;
     347  for (i = 0 ; i < me - ms ; i++){
     348    w += b[i] * b[i]; /* w = b.b */
     349  }
     350  MPI_Allreduce(&w, &normb, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); 
     351  normb = sqrt(normb);
     352 
     353  /* compute r0 = b - Ax */
     354  for (i = 0 ; i < n ; i++) work[i] = 0.0;
     355  for (i = ms ; i < me ; i++){
     356    li = i - ms;
     357    work[i] = b[li];
     358  }
     359  for (j = ms ; j < me ; j++){
     360    lj = j - ms;
     361    for (i = 0 ; i < n ; i++){
     362      work[i] -= mat[i + lj * n] * x[lj];
     363    }
     364  }
     365  for (cpu = 0 ; cpu < numproc ; cpu++){
     366    MPI_Reduce(&work[DIM * ostart[cpu]], r, DIM * (oend[cpu] - ostart[cpu]),
     367               MPI_DOUBLE, MPI_SUM, cpu, MPI_COMM_WORLD);
     368  }
     369 
     370  w = 0.0;
     371  for (i = 0 ; i < me - ms ; i++){
     372    p[i] = r[i];  /* p = r */
     373    w += r[i] * r[i]; /* rr = r.r */
     374  }
     375  MPI_Allreduce(&w, &rr, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
     376 
     377  /* cg iteration */
     378  count = 0;
     379  while(rr  > tol * tol * normb * normb){   
     380    /* Ap = A*p */
     381    for (i = 0 ; i < n ; i++) work[i] = 0.0;
     382    for (j = ms ; j < me ; j++){
     383      lj = j - ms;
     384      for (i = 0 ; i < n ; i++){
     385        work[i] += mat[i + lj * n] * p[lj];
     386      }
     387    }
     388    for (cpu = 0 ; cpu < numproc ; cpu++){
     389      MPI_Reduce(&work[DIM * ostart[cpu]], Ap, DIM * (oend[cpu] - ostart[cpu]),
     390                 MPI_DOUBLE, MPI_SUM, cpu, MPI_COMM_WORLD);
     391    }
     392 
     393    /* pAp = p.Ap */
     394    w = 0.0;
     395    for (i = 0 ; i < me - ms ; i++){
     396      w += p[i] * Ap[i];
     397    }
     398    MPI_Allreduce(&w, &pAp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
     399 
     400    /* alpha = r.r / p.Ap */
     401    alpha = rr / pAp;
     402 
     403    /* Beta */
     404    for (i = 0 ; i < me - ms ; i++){
     405      x[i] += alpha * p[i];  /* x += alpha * p */
     406      r[i] -= alpha * Ap[i]; /* r -= alpha * Ap */
     407    }   
     408 
     409    /* rr1 = r.r */     
     410    w = 0.0;
     411    for (i = 0 ; i < me - ms ; i++){
     412      w += r[i] * r[i];
     413    }
     414    MPI_Allreduce(&w, &rr1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
     415 
     416    beta = rr1 / rr;
     417 
     418    /* p = r + beta * p */
     419    for (i = 0 ; i < me - ms ; i++){
     420      p[i] = r[i] + beta * p[i];
     421    }
     422 
     423    rr = rr1;
     424    count++;
     425  }
     426  if (myid == 0) printf("count=%d\n",count);
     427 
     428  /* Deallocate Working Spaces */
     429  free(r);
     430  free(p);
     431  free(Ap);
     432  free(work);
     433}
     434 
     435 
     436/* compute velocity */
     437void slet2d_velocity_dist(int np, int *ostart,int *oend,double *loc, double *foc, double ep,
     438                          double *cloc,double *cvel){
     439  int i,p;
     440  double r;
     441  double dx,dy;
     442  double tr1,tr2;
     443  double ov[2];
     444  int myid,numproc;
     445  int li;
     446 
     447  /* get myid and # of processors */
     448  MPI_Comm_size(MPI_COMM_WORLD,&numproc);
     449  MPI_Comm_rank(MPI_COMM_WORLD,&myid);
     450 
     451  /* zeros */
     452  ov[0] = 0.0;
     453  ov[1] = 0.0;
     454 
     455  /* loop for partilces (blob) */
     456  for (i = ostart[myid] ; i < oend[myid] ; i++){
     457    li = i - ostart[myid];
     458    dx = cloc[0] - loc[li * DIM    ];
     459    dy = cloc[1] - loc[li * DIM + 1];
     460    r = sqrt(dx * dx + dy * dy);
     461 
     462    tr1 = term1(r,ep) / (4.0 * M_PI);
     463    tr2 = term2(r,ep) / (4.0 * M_PI);       
     464 
     465    tr2 *= foc[li * DIM] * dx + foc[li * DIM + 1] * dy;
     466 
     467    ov[0] += -foc[li * DIM    ] * tr1 + tr2 * dx;
     468    ov[1] += -foc[li * DIM + 1] * tr1 + tr2 * dy;
     469  }
     470 
     471  /* sum up to CPU0 */
     472  MPI_Reduce(ov, cvel, 2, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); 
     473}
     474 
     475double term1(double r,double ep){
     476  double sq;
     477 
     478  sq = sqrt(r * r + ep * ep);
     479 
     480  return log(sq + ep) - ep * (sq + 2.0 * ep) / (sq + ep) / sq;
     481}
     482 
     483 
     484double term2(double r,double ep){
     485  double sq;
     486 
     487  sq = sqrt(r * r + ep * ep);
     488 
     489  return (sq + 2.0 * ep) / (sq + ep) / (sq + ep) / sq;
     490}
     491}}}