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

Timestamp:

05/14/15 14:52:05 (9 years ago)

Author:

cmaggio

Comment:

migrated content from ccs wiki

cypress/Programming/StokesFlow/ex32.c

@@ +1 @@
+= Ex32.c =
+{{{#!C
+/*
+  Large Scale Computing
+  Stokes Flow in a Cavity
+  ex32.c
+ 
+  Need Lapack and Blas Libraries
+*/
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+#include<time.h>
+#include"stokeslet2d.h"
+ 
+#define EPSILON 0.005 /* blob size*/
+#define INTGRID 51    /* # of x-grid lines for internal velocity */
+ 
+/* min and max macros */
+#define max(a,b)        (((a) > (b)) ? (a) : (b))
+#define min(a,b)        (((a) < (b)) ? (a) : (b))
+ 
+double t_cost(void );
+int main(int , char **);
+int main(int argc, char **argv){
+  double dp;
+  int numpdepth;
+  int numpwidth;
+  int numparticles; 
+  double *loc; /* particle location */
+  double *vel; /* particle velocity */
+  double *foc; /* particle force */
+  double *mat; /* 2Nx2N dense matrix */
+  double *cloc; /* array for internal points */
+  double *cvel; /* array for internal velocity */
+  int i,j;
+  int nyg;
+  FILE *fp;
+ 
+  if (argc < 2){
+    printf("Usage:%s [Depth of Cavity]\n",argv[0]);
+    exit(-1);
+  }
+ 
+  /* get inputed depth */
+  dp = atof(argv[1]);
+ 
+  /* # of particles in depth */
+  numpdepth = (int)(dp / EPSILON + 0.5);
+ 
+  /* # of particles in width */
+  numpwidth = (int)(1.0 / EPSILON + 0.5);
+ 
+  /* total # od particles */
+  numparticles = numpdepth * 2 + numpwidth * 2;
+  printf("Total # of Particles=%d\n",numparticles);
+ 
+  /* Allocate Space */
+  /* DIM=2 defined by 'stokeslet2d.h' */
+  loc = (double *)calloc(numparticles * DIM,sizeof(double));
+  vel = (double *)calloc(numparticles * DIM,sizeof(double));
+  foc = (double *)calloc(numparticles * DIM,sizeof(double));
+  mat = (double *)calloc(numparticles * DIM * numparticles * DIM,sizeof(double));
+  if ((loc == (double *)NULL) ||
+      (vel == (double *)NULL) || 
+      (foc == (double *)NULL) || 
+      (mat == (double *)NULL) ){
+    printf("Can't allocate memory!!!\n");
+    exit(-1);
+  }
+ 
+  /* set location & velocity of particles (blob)*/
+  for (i = 0 ; i < numparticles ; i++){
+    foc[i * DIM] = 0.0;
+    foc[i * DIM + 1] = 0.0;
+    if ((i >= 0) && (i <= numpwidth)){ /* top */
+      loc[i * DIM] = -0.5 + EPSILON * (double)i; /* x */
+      loc[i * DIM + 1] = 0.0; /* y */
+      vel[i * DIM] = 1.0;
+      vel[i * DIM + 1] = 0.0;
+    }else{
+      if (i <= (numpwidth + numpdepth)){ /* right wall */      
+        loc[i * DIM] = 0.5; /* x */
+        loc[i * DIM + 1] = -EPSILON * (double)(i - numpwidth); /* y */  
+        vel[i * DIM] = 0.0;
+        vel[i * DIM + 1] = 0.0;
+      }else{
+        if (i <= (2 * numpwidth + numpdepth)){ /* bottom */
+          loc[i * DIM] = 0.5 - EPSILON * (double)(i - (numpwidth + numpdepth)); /* x */
+          loc[i * DIM + 1] = -EPSILON * numpdepth; /* y */
+          vel[i * DIM] = 0.0;
+          vel[i * DIM + 1] = 0.0;         
+        }else{                 /* left wall */    
+          loc[i * DIM] = -0.5; /* x */
+          loc[i * DIM + 1] = -EPSILON * (double)((2 * numpwidth + 2 * numpdepth) - i); /* y */  
+          vel[i * DIM] = 0.0;
+          vel[i * DIM + 1] = 0.0;         
+        }
+      }
+    }
+  }
+ 
+  /* make 2Nx2N Matrix */
+  t_cost();
+  slet2d_mkMatrix(numparticles,loc,EPSILON,mat);
+  printf("setting matrix %f sec\n",t_cost());
+ 
+  /* Sovle linear ststem */
+  slet2d_solve(numparticles,mat,vel,foc);
+  //slet2d_solve_CG(numparticles,mat,vel,foc);
+  //slet2d_solve_GMRES(numparticles,mat,vel,foc);
+ 
+  printf("Sovle linear ststem %f sec\n",t_cost());
+ 
+  /* check the solution */
+  slet2d_mkMatrix(numparticles,loc,EPSILON,mat);
+  check_solution(numparticles,mat,vel,foc);
+ 
+  /* deallocate big matrix */
+  free(mat);
+ 
+  /* out particle (blob) data */
+  fp = fopen("particle.dat","w");
+  for (i = 0 ; i < numparticles ; i++){
+    fprintf(fp,"%e %e %e %e %e %e\n",
+            loc[i * DIM],loc[i * DIM + 1],
+            vel[i * DIM],vel[i * DIM + 1],
+            foc[i * DIM],foc[i * DIM + 1]);
+  }
+  fclose(fp);
+ 
+  /* 
+     compute internal velocity 
+  */
+  nyg = (int)(EPSILON * (double)(numpdepth * (INTGRID - 1)));
+  cvel = (double *)calloc(INTGRID * nyg * DIM, sizeof(double));  
+  cloc = (double *)calloc(INTGRID * nyg * DIM, sizeof(double));
+ 
+  if ((cloc == (double *)NULL) ||
+      (cvel == (double *)NULL)){
+    printf("Can't allocate memory!!!\n");
+    exit(-1);
+  }
+  /* setting location */
+  for (j = 0 ; j < nyg ; j++){
+    for (i = 0 ; i < INTGRID ; i++){
+      cloc[DIM * (i + INTGRID * j)    ] = -0.5 + (double)i / (double)(INTGRID - 1);
+      cloc[DIM * (i + INTGRID * j) + 1] = -(double)j / (double)(INTGRID - 1);
+    }
+  }
+ 
+  /* compute velocities */
+  t_cost();
+  slet2d_velocity(numparticles, loc, foc, EPSILON, INTGRID * nyg, cloc, cvel);  
+  printf("Compute internal velocity %f sec\n",t_cost());
+ 
+  /* out velocities */
+  fp = fopen("res.dat","w");
+  for (j = 0 ; j < nyg ; j++){
+    for (i = 0 ; i < INTGRID ; i++){
+      fprintf(fp,"%e %e %e %e\n",
+              cloc[DIM * (i + INTGRID * j)    ], cloc[DIM * (i + INTGRID * j) + 1],
+              cvel[DIM * (i + INTGRID * j)    ], cvel[DIM * (i + INTGRID * j) + 1]);
+    }
+  }
+  fclose(fp);
+ 
+  /* free */
+  free(cloc);
+  free(cvel);
+ 
+  free(loc);
+  free(vel);
+  free(foc);
+}
+ 
+/* 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;
+}
+stokeslet2d.c
+
+/*
+  Large Scale Computing
+  Stokeslet for 2D
+ 
+  Ricardo Cortez,"The Method of Regularized Stokeslets", 
+  (2001), SIAM J. Sci. Comput., Vol.23, No.4, pp.1204
+ 
+  assuming mu = 1
+*/
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+#include"stokeslet2d.h"
+#include"gmres.h"
+ 
+double term1(double ,double );
+double term2(double ,double );
+/* Blas Function */
+double ddot_(int *,double *,int *,double *,int *);
+double dnrm2_(int *, double *, int *);
+ 
+/* make Matrix */
+void slet2d_mkMatrix(int np,double *loc,double ep,double *mat){
+  int i,j;
+  double r;
+  double dx,dy;
+  double tr1,tr2;
+ 
+  /*zeros mat*/
+  for (i = 0 ; i < DIM * DIM * np * np ; i++) mat[i] = 0.0;
+ 
+  /* make mat */
+  /* Because of symmetric, we can loop over half */
+  for (i = 0 ; i < np ; i++){
+    for (j = i ; j < np ; j++){
+      dx = loc[i * DIM    ] - loc[j * DIM    ];
+      dy = loc[i * DIM + 1] - loc[j * DIM + 1];
+      r = sqrt(dx * dx + dy * dy);
+ 
+      tr1 = term1(r,ep) / (4.0 * M_PI);
+      tr2 = term2(r,ep) / (4.0 * M_PI); 
+ 
+      /* diagoanl elements and upper half */
+      mat[i * DIM +      j * DIM      * np * DIM] += -tr1 + tr2 * dx * dx;
+      mat[i * DIM +     (j * DIM + 1) * np * DIM] +=        tr2 * dx * dy;
+      mat[i * DIM + 1 +  j * DIM      * np * DIM] +=        tr2 * dx * dy;
+      mat[i * DIM + 1 + (j * DIM + 1) * np * DIM] += -tr1 + tr2 * dy * dy;      
+      if (i != j){
+        /* lower half */
+        mat[j * DIM +      i * DIM      * np * DIM] += -tr1 + tr2 * dx * dx;
+        mat[j * DIM +     (i * DIM + 1) * np * DIM] +=        tr2 * dx * dy;
+        mat[j * DIM + 1 +  i * DIM      * np * DIM] +=        tr2 * dx * dy;
+        mat[j * DIM + 1 + (i * DIM + 1) * np * DIM] += -tr1 + tr2 * dy * dy;
+      }
+    }
+  }
+}
+ 
+/* Sovle Liear ststem */
+/* Use Lapack */
+void slet2d_solve(int np,double *mat,double *vel,double *foc){
+  int nrhs;
+  int n;
+  int lda;
+  int ldb;
+  int *ipiv;
+  int i;
+  int info;
+ 
+  /*solve Ax=b */
+  /* A (amat) is n x n matrix */
+  nrhs = 1; /* # of RHS is 1 */
+  n = DIM * np;
+  lda = DIM * np;
+  ldb = DIM * np;
+  ipiv = calloc(DIM * np,sizeof(int)); /* The pivot indices that define the permutation matrix P;
+                                          row i of the matrix was interchanged with row IPIV(i). */
+  if (ipiv == NULL) exit(-1);
+ 
+  /* set RHS */
+  for (i = 0 ; i < n ; i++) foc[i] = vel[i];
+ 
+  /* solve with LAPACK */
+  dgesv_(&n, &nrhs, mat, &lda, ipiv, foc, &ldb, &info);
+ 
+  /* check */
+  /* if (info == 0) printf("successfully done\n"); */
+  if (info < 0){
+    printf("the %d-th argument had an illegal value\n",-info);
+    exit(-1);
+  }
+  if (info > 0){
+    printf("U(%d,%d) is exactly zero.\n",info,info);
+    exit(-1);
+  }
+ 
+  /* free */
+  free(ipiv);
+}
+ 
+/* CG use blas */
+void slet2d_solve_CG(int np,double *mat,double *b,double *x){
+  double tol = 1.0e-8;
+  double *r;
+  double *p;
+  double *Ap;
+  double rr,pAp,rr1,alpha,beta;
+  int i,j,count;
+  double al,bt,normb;
+  char tr[1] = "N";
+  int inc,n;
+ 
+  n = DIM * np;
+  inc = 1;
+ 
+  /* Check if the solution is trivial. */
+  normb = dnrm2_(&n, b, &inc);
+  if (normb <= 0.0){
+    for (i = 0 ; i < n ; i++) x[i] = 0.0;
+    return;
+  }  
+ 
+  /* Allocate Working Spaces */
+  r = (double *)calloc(n, sizeof(double));
+  p = (double *)calloc(n, sizeof(double));
+  Ap = (double *)calloc(n, sizeof(double));
+  if ((r == NULL) || (p == NULL) || (Ap == NULL)){
+    printf("memory allocation failed\n");
+    return ;
+  }
+ 
+  /* compute r0 = b - Ax */
+  al = -1.0;
+  bt = 1.0;
+  dcopy_(&n,b,&inc,r,&inc); /* r = b */
+  dgemv_(tr,&n,&n,&al,mat,&n,x,&inc,&bt,r,&inc);  
+ 
+  rr = 0.0;
+  dcopy_(&n,r,&inc,p,&inc); /* p = r */
+  rr = ddot_(&n,r,&inc,r,&inc); /* rr = r.r */ 
+  printf("rr=%e\n",rr);
+ 
+  /* cg iteration */
+  count = 0;
+  while(rr  > tol * tol * normb * normb){    
+    // Ap = A*p
+    al = 1.0;
+    bt = 0.0;
+    dgemv_(tr,&n,&n,&al,mat,&n,p,&inc,&bt,Ap,&inc);  
+ 
+    /* alpha = r.r / p.Ap */
+    pAp = 0.0;
+    pAp = ddot_(&n,p,&inc,Ap,&inc); /* pAp = p.Ap */ 
+    alpha = rr / pAp;
+ 
+    /* Beta */
+    rr1 = 0.0;
+    daxpy_(&n,&alpha,p,&inc,x,&inc); /* x += alpha * p */
+    al = -alpha;
+    daxpy_(&n,&al,Ap,&inc,r,&inc); /* r -= alpha * Ap */
+    rr1 = ddot_(&n,r,&inc,r,&inc); /* rr1 = r.r */     
+ 
+    beta = rr1 / rr;
+    /* p = r + beta * p  no blas routine :( */
+    for (i = 0 ; i < n ; i++){
+      p[i] = r[i] + beta * p[i];
+    }
+ 
+    rr = rr1;
+    count++;
+  }
+ 
+  /* Deallocate Working Spaces */
+  free(r); 
+  free(p); 
+  free(Ap);
+}
+ 
+void check_solution(int np,double *mat,double *b,double *x){
+  double *r;
+  double al,bt,rr;
+  char tr[1] = "N";
+  int inc,n;
+ 
+  n = DIM * np;
+  inc = 1;
+ 
+  /* Allocate Working Spaces */
+  r = (double *)calloc(n, sizeof(double));  
+ 
+  /* compute r0 = b - Ax */
+  al = -1.0;
+  bt = 1.0;
+  dcopy_(&n,b,&inc,r,&inc); /* r = b */
+  dgemv_(tr,&n,&n,&al,mat,&n,x,&inc,&bt,r,&inc);  
+ 
+  /* rr = |r| */
+  rr = dnrm2_(&n, r, &inc);  
+ 
+  printf("|b - Ax|=%e\n",rr);
+  free(r);
+}
+ 
+/* Use Lapack */
+void slet2d_solve_GMRES(int np,double *mat,double *vel,double *foc){
+  int n;
+  int restart;
+  int max_iterations;
+  double tol;
+ 
+  /* set GMRES parameters */
+  restart = 100;
+  max_iterations = 100;
+  tol = 1.0e-8;
+ 
+  /*solve Ax=b */
+  /* A (amat) is n x n matrix */
+  n = DIM * np;
+ 
+  /* solve with GMRES */
+  gmres_solve(mat, foc, vel, n, restart, &tol, &max_iterations);
+}
+ 
+/* compute velocity */
+void slet2d_velocity(int np, double *loc, double *foc, double ep, int nump, double *cloc,double *cvel){
+  int i,p;
+  double r;
+  double dx,dy;
+  double tr1,tr2;
+ 
+  /* loop for grid */
+  for (p = 0 ; p < nump ; p++){
+    /* zeros */
+    cvel[p * DIM    ] = 0.0;
+    cvel[p * DIM + 1] = 0.0;
+ 
+    /* loop for partilces (blob) */
+    for (i = 0 ; i < np ; i++){
+      dx = cloc[p * DIM    ] - loc[i * DIM    ];
+      dy = cloc[p * DIM + 1] - loc[i * DIM + 1];
+      r = sqrt(dx * dx + dy * dy);
+ 
+      tr1 = term1(r,ep) / (4.0 * M_PI);
+      tr2 = term2(r,ep) / (4.0 * M_PI);       
+ 
+      tr2 *= foc[i * DIM] * dx + foc[i * DIM + 1] * dy; 
+ 
+      cvel[p * DIM    ] += -foc[i * DIM    ] * tr1 + tr2 * dx;
+      cvel[p * DIM + 1] += -foc[i * DIM + 1] * tr1 + tr2 * dy;
+    }
+  }
+}
+ 
+double term1(double r,double ep){
+  double sq;
+ 
+  sq = sqrt(r * r + ep * ep);
+ 
+  return log(sq + ep) - ep * (sq + 2.0 * ep) / (sq + ep) / sq;
+}
+ 
+ 
+double term2(double r,double ep){
+  double sq;
+ 
+  sq = sqrt(r * r + ep * ep);
+ 
+  return (sq + 2.0 * ep) / (sq + ep) / (sq + ep) / sq;
+}}}}