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


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

migrated content from ccs wiki

Legend:

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

    v1 v1  
     1= Ex32.c =
     2{{{#!C
     3/*
     4  Large Scale Computing
     5  Stokes Flow in a Cavity
     6  ex32.c
     7 
     8  Need Lapack and Blas Libraries
     9*/
     10#include<stdio.h>
     11#include<stdlib.h>
     12#include<math.h>
     13#include<time.h>
     14#include"stokeslet2d.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; /* array for internal points */
     35  double *cvel; /* array for internal velocity */
     36  int i,j;
     37  int nyg;
     38  FILE *fp;
     39 
     40  if (argc < 2){
     41    printf("Usage:%s [Depth of Cavity]\n",argv[0]);
     42    exit(-1);
     43  }
     44 
     45  /* get inputed depth */
     46  dp = atof(argv[1]);
     47 
     48  /* # of particles in depth */
     49  numpdepth = (int)(dp / EPSILON + 0.5);
     50 
     51  /* # of particles in width */
     52  numpwidth = (int)(1.0 / EPSILON + 0.5);
     53 
     54  /* total # od particles */
     55  numparticles = numpdepth * 2 + numpwidth * 2;
     56  printf("Total # of Particles=%d\n",numparticles);
     57 
     58  /* Allocate Space */
     59  /* DIM=2 defined by 'stokeslet2d.h' */
     60  loc = (double *)calloc(numparticles * DIM,sizeof(double));
     61  vel = (double *)calloc(numparticles * DIM,sizeof(double));
     62  foc = (double *)calloc(numparticles * DIM,sizeof(double));
     63  mat = (double *)calloc(numparticles * DIM * numparticles * DIM,sizeof(double));
     64  if ((loc == (double *)NULL) ||
     65      (vel == (double *)NULL) ||
     66      (foc == (double *)NULL) ||
     67      (mat == (double *)NULL) ){
     68    printf("Can't allocate memory!!!\n");
     69    exit(-1);
     70  }
     71 
     72  /* set location & velocity of particles (blob)*/
     73  for (i = 0 ; i < numparticles ; i++){
     74    foc[i * DIM] = 0.0;
     75    foc[i * DIM + 1] = 0.0;
     76    if ((i >= 0) && (i <= numpwidth)){ /* top */
     77      loc[i * DIM] = -0.5 + EPSILON * (double)i; /* x */
     78      loc[i * DIM + 1] = 0.0; /* y */
     79      vel[i * DIM] = 1.0;
     80      vel[i * DIM + 1] = 0.0;
     81    }else{
     82      if (i <= (numpwidth + numpdepth)){ /* right wall */     
     83        loc[i * DIM] = 0.5; /* x */
     84        loc[i * DIM + 1] = -EPSILON * (double)(i - numpwidth); /* y */ 
     85        vel[i * DIM] = 0.0;
     86        vel[i * DIM + 1] = 0.0;
     87      }else{
     88        if (i <= (2 * numpwidth + numpdepth)){ /* bottom */
     89          loc[i * DIM] = 0.5 - EPSILON * (double)(i - (numpwidth + numpdepth)); /* x */
     90          loc[i * DIM + 1] = -EPSILON * numpdepth; /* y */
     91          vel[i * DIM] = 0.0;
     92          vel[i * DIM + 1] = 0.0;         
     93        }else{                 /* left wall */   
     94          loc[i * DIM] = -0.5; /* x */
     95          loc[i * DIM + 1] = -EPSILON * (double)((2 * numpwidth + 2 * numpdepth) - i); /* y */ 
     96          vel[i * DIM] = 0.0;
     97          vel[i * DIM + 1] = 0.0;         
     98        }
     99      }
     100    }
     101  }
     102 
     103  /* make 2Nx2N Matrix */
     104  t_cost();
     105  slet2d_mkMatrix(numparticles,loc,EPSILON,mat);
     106  printf("setting matrix %f sec\n",t_cost());
     107 
     108  /* Sovle linear ststem */
     109  slet2d_solve(numparticles,mat,vel,foc);
     110  //slet2d_solve_CG(numparticles,mat,vel,foc);
     111  //slet2d_solve_GMRES(numparticles,mat,vel,foc);
     112 
     113  printf("Sovle linear ststem %f sec\n",t_cost());
     114 
     115  /* check the solution */
     116  slet2d_mkMatrix(numparticles,loc,EPSILON,mat);
     117  check_solution(numparticles,mat,vel,foc);
     118 
     119  /* deallocate big matrix */
     120  free(mat);
     121 
     122  /* out particle (blob) data */
     123  fp = fopen("particle.dat","w");
     124  for (i = 0 ; i < numparticles ; i++){
     125    fprintf(fp,"%e %e %e %e %e %e\n",
     126            loc[i * DIM],loc[i * DIM + 1],
     127            vel[i * DIM],vel[i * DIM + 1],
     128            foc[i * DIM],foc[i * DIM + 1]);
     129  }
     130  fclose(fp);
     131 
     132  /*
     133     compute internal velocity
     134  */
     135  nyg = (int)(EPSILON * (double)(numpdepth * (INTGRID - 1)));
     136  cvel = (double *)calloc(INTGRID * nyg * DIM, sizeof(double)); 
     137  cloc = (double *)calloc(INTGRID * nyg * DIM, sizeof(double));
     138 
     139  if ((cloc == (double *)NULL) ||
     140      (cvel == (double *)NULL)){
     141    printf("Can't allocate memory!!!\n");
     142    exit(-1);
     143  }
     144  /* setting location */
     145  for (j = 0 ; j < nyg ; j++){
     146    for (i = 0 ; i < INTGRID ; i++){
     147      cloc[DIM * (i + INTGRID * j)    ] = -0.5 + (double)i / (double)(INTGRID - 1);
     148      cloc[DIM * (i + INTGRID * j) + 1] = -(double)j / (double)(INTGRID - 1);
     149    }
     150  }
     151 
     152  /* compute velocities */
     153  t_cost();
     154  slet2d_velocity(numparticles, loc, foc, EPSILON, INTGRID * nyg, cloc, cvel); 
     155  printf("Compute internal velocity %f sec\n",t_cost());
     156 
     157  /* out velocities */
     158  fp = fopen("res.dat","w");
     159  for (j = 0 ; j < nyg ; j++){
     160    for (i = 0 ; i < INTGRID ; i++){
     161      fprintf(fp,"%e %e %e %e\n",
     162              cloc[DIM * (i + INTGRID * j)    ], cloc[DIM * (i + INTGRID * j) + 1],
     163              cvel[DIM * (i + INTGRID * j)    ], cvel[DIM * (i + INTGRID * j) + 1]);
     164    }
     165  }
     166  fclose(fp);
     167 
     168  /* free */
     169  free(cloc);
     170  free(cvel);
     171 
     172  free(loc);
     173  free(vel);
     174  free(foc);
     175}
     176 
     177/* for timing */
     178double t_cost(void ){
     179  static clock_t ptm;
     180  static int fl = 0;
     181  clock_t tm;
     182  double sec;
     183 
     184  tm = clock();
     185  if (fl == 0) ptm = tm;
     186  fl = 1;
     187 
     188  sec = (double)(tm - ptm) / (double)CLOCKS_PER_SEC;
     189 
     190  ptm = tm;
     191 
     192  return sec;
     193}
     194stokeslet2d.c
     195
     196/*
     197  Large Scale Computing
     198  Stokeslet for 2D
     199 
     200  Ricardo Cortez,"The Method of Regularized Stokeslets",
     201  (2001), SIAM J. Sci. Comput., Vol.23, No.4, pp.1204
     202 
     203  assuming mu = 1
     204*/
     205#include<stdio.h>
     206#include<stdlib.h>
     207#include<math.h>
     208#include"stokeslet2d.h"
     209#include"gmres.h"
     210 
     211double term1(double ,double );
     212double term2(double ,double );
     213/* Blas Function */
     214double ddot_(int *,double *,int *,double *,int *);
     215double dnrm2_(int *, double *, int *);
     216 
     217/* make Matrix */
     218void slet2d_mkMatrix(int np,double *loc,double ep,double *mat){
     219  int i,j;
     220  double r;
     221  double dx,dy;
     222  double tr1,tr2;
     223 
     224  /*zeros mat*/
     225  for (i = 0 ; i < DIM * DIM * np * np ; i++) mat[i] = 0.0;
     226 
     227  /* make mat */
     228  /* Because of symmetric, we can loop over half */
     229  for (i = 0 ; i < np ; i++){
     230    for (j = i ; j < np ; j++){
     231      dx = loc[i * DIM    ] - loc[j * DIM    ];
     232      dy = loc[i * DIM + 1] - loc[j * DIM + 1];
     233      r = sqrt(dx * dx + dy * dy);
     234 
     235      tr1 = term1(r,ep) / (4.0 * M_PI);
     236      tr2 = term2(r,ep) / (4.0 * M_PI);
     237 
     238      /* diagoanl elements and upper half */
     239      mat[i * DIM +      j * DIM      * np * DIM] += -tr1 + tr2 * dx * dx;
     240      mat[i * DIM +     (j * DIM + 1) * np * DIM] +=        tr2 * dx * dy;
     241      mat[i * DIM + 1 +  j * DIM      * np * DIM] +=        tr2 * dx * dy;
     242      mat[i * DIM + 1 + (j * DIM + 1) * np * DIM] += -tr1 + tr2 * dy * dy;     
     243      if (i != j){
     244        /* lower half */
     245        mat[j * DIM +      i * DIM      * np * DIM] += -tr1 + tr2 * dx * dx;
     246        mat[j * DIM +     (i * DIM + 1) * np * DIM] +=        tr2 * dx * dy;
     247        mat[j * DIM + 1 +  i * DIM      * np * DIM] +=        tr2 * dx * dy;
     248        mat[j * DIM + 1 + (i * DIM + 1) * np * DIM] += -tr1 + tr2 * dy * dy;
     249      }
     250    }
     251  }
     252}
     253 
     254/* Sovle Liear ststem */
     255/* Use Lapack */
     256void slet2d_solve(int np,double *mat,double *vel,double *foc){
     257  int nrhs;
     258  int n;
     259  int lda;
     260  int ldb;
     261  int *ipiv;
     262  int i;
     263  int info;
     264 
     265  /*solve Ax=b */
     266  /* A (amat) is n x n matrix */
     267  nrhs = 1; /* # of RHS is 1 */
     268  n = DIM * np;
     269  lda = DIM * np;
     270  ldb = DIM * np;
     271  ipiv = calloc(DIM * np,sizeof(int)); /* The pivot indices that define the permutation matrix P;
     272                                          row i of the matrix was interchanged with row IPIV(i). */
     273  if (ipiv == NULL) exit(-1);
     274 
     275  /* set RHS */
     276  for (i = 0 ; i < n ; i++) foc[i] = vel[i];
     277 
     278  /* solve with LAPACK */
     279  dgesv_(&n, &nrhs, mat, &lda, ipiv, foc, &ldb, &info);
     280 
     281  /* check */
     282  /* if (info == 0) printf("successfully done\n"); */
     283  if (info < 0){
     284    printf("the %d-th argument had an illegal value\n",-info);
     285    exit(-1);
     286  }
     287  if (info > 0){
     288    printf("U(%d,%d) is exactly zero.\n",info,info);
     289    exit(-1);
     290  }
     291 
     292  /* free */
     293  free(ipiv);
     294}
     295 
     296/* CG use blas */
     297void slet2d_solve_CG(int np,double *mat,double *b,double *x){
     298  double tol = 1.0e-8;
     299  double *r;
     300  double *p;
     301  double *Ap;
     302  double rr,pAp,rr1,alpha,beta;
     303  int i,j,count;
     304  double al,bt,normb;
     305  char tr[1] = "N";
     306  int inc,n;
     307 
     308  n = DIM * np;
     309  inc = 1;
     310 
     311  /* Check if the solution is trivial. */
     312  normb = dnrm2_(&n, b, &inc);
     313  if (normb <= 0.0){
     314    for (i = 0 ; i < n ; i++) x[i] = 0.0;
     315    return;
     316  } 
     317 
     318  /* Allocate Working Spaces */
     319  r = (double *)calloc(n, sizeof(double));
     320  p = (double *)calloc(n, sizeof(double));
     321  Ap = (double *)calloc(n, sizeof(double));
     322  if ((r == NULL) || (p == NULL) || (Ap == NULL)){
     323    printf("memory allocation failed\n");
     324    return ;
     325  }
     326 
     327  /* compute r0 = b - Ax */
     328  al = -1.0;
     329  bt = 1.0;
     330  dcopy_(&n,b,&inc,r,&inc); /* r = b */
     331  dgemv_(tr,&n,&n,&al,mat,&n,x,&inc,&bt,r,&inc); 
     332 
     333  rr = 0.0;
     334  dcopy_(&n,r,&inc,p,&inc); /* p = r */
     335  rr = ddot_(&n,r,&inc,r,&inc); /* rr = r.r */
     336  printf("rr=%e\n",rr);
     337 
     338  /* cg iteration */
     339  count = 0;
     340  while(rr  > tol * tol * normb * normb){   
     341    // Ap = A*p
     342    al = 1.0;
     343    bt = 0.0;
     344    dgemv_(tr,&n,&n,&al,mat,&n,p,&inc,&bt,Ap,&inc); 
     345 
     346    /* alpha = r.r / p.Ap */
     347    pAp = 0.0;
     348    pAp = ddot_(&n,p,&inc,Ap,&inc); /* pAp = p.Ap */
     349    alpha = rr / pAp;
     350 
     351    /* Beta */
     352    rr1 = 0.0;
     353    daxpy_(&n,&alpha,p,&inc,x,&inc); /* x += alpha * p */
     354    al = -alpha;
     355    daxpy_(&n,&al,Ap,&inc,r,&inc); /* r -= alpha * Ap */
     356    rr1 = ddot_(&n,r,&inc,r,&inc); /* rr1 = r.r */     
     357 
     358    beta = rr1 / rr;
     359    /* p = r + beta * p  no blas routine :( */
     360    for (i = 0 ; i < n ; i++){
     361      p[i] = r[i] + beta * p[i];
     362    }
     363 
     364    rr = rr1;
     365    count++;
     366  }
     367 
     368  /* Deallocate Working Spaces */
     369  free(r);
     370  free(p);
     371  free(Ap);
     372}
     373 
     374void check_solution(int np,double *mat,double *b,double *x){
     375  double *r;
     376  double al,bt,rr;
     377  char tr[1] = "N";
     378  int inc,n;
     379 
     380  n = DIM * np;
     381  inc = 1;
     382 
     383  /* Allocate Working Spaces */
     384  r = (double *)calloc(n, sizeof(double)); 
     385 
     386  /* compute r0 = b - Ax */
     387  al = -1.0;
     388  bt = 1.0;
     389  dcopy_(&n,b,&inc,r,&inc); /* r = b */
     390  dgemv_(tr,&n,&n,&al,mat,&n,x,&inc,&bt,r,&inc); 
     391 
     392  /* rr = |r| */
     393  rr = dnrm2_(&n, r, &inc); 
     394 
     395  printf("|b - Ax|=%e\n",rr);
     396  free(r);
     397}
     398 
     399/* Use Lapack */
     400void slet2d_solve_GMRES(int np,double *mat,double *vel,double *foc){
     401  int n;
     402  int restart;
     403  int max_iterations;
     404  double tol;
     405 
     406  /* set GMRES parameters */
     407  restart = 100;
     408  max_iterations = 100;
     409  tol = 1.0e-8;
     410 
     411  /*solve Ax=b */
     412  /* A (amat) is n x n matrix */
     413  n = DIM * np;
     414 
     415  /* solve with GMRES */
     416  gmres_solve(mat, foc, vel, n, restart, &tol, &max_iterations);
     417}
     418 
     419/* compute velocity */
     420void slet2d_velocity(int np, double *loc, double *foc, double ep, int nump, double *cloc,double *cvel){
     421  int i,p;
     422  double r;
     423  double dx,dy;
     424  double tr1,tr2;
     425 
     426  /* loop for grid */
     427  for (p = 0 ; p < nump ; p++){
     428    /* zeros */
     429    cvel[p * DIM    ] = 0.0;
     430    cvel[p * DIM + 1] = 0.0;
     431 
     432    /* loop for partilces (blob) */
     433    for (i = 0 ; i < np ; i++){
     434      dx = cloc[p * DIM    ] - loc[i * DIM    ];
     435      dy = cloc[p * DIM + 1] - loc[i * DIM + 1];
     436      r = sqrt(dx * dx + dy * dy);
     437 
     438      tr1 = term1(r,ep) / (4.0 * M_PI);
     439      tr2 = term2(r,ep) / (4.0 * M_PI);       
     440 
     441      tr2 *= foc[i * DIM] * dx + foc[i * DIM + 1] * dy;
     442 
     443      cvel[p * DIM    ] += -foc[i * DIM    ] * tr1 + tr2 * dx;
     444      cvel[p * DIM + 1] += -foc[i * DIM + 1] * tr1 + tr2 * dy;
     445    }
     446  }
     447}
     448 
     449double term1(double r,double ep){
     450  double sq;
     451 
     452  sq = sqrt(r * r + ep * ep);
     453 
     454  return log(sq + ep) - ep * (sq + 2.0 * ep) / (sq + ep) / sq;
     455}
     456 
     457 
     458double term2(double r,double ep){
     459  double sq;
     460 
     461  sq = sqrt(r * r + ep * ep);
     462 
     463  return (sq + 2.0 * ep) / (sq + ep) / (sq + ep) / sq;
     464}}}}