simulation.c

C语言编写的并行计算柏松矩阵

#include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include "datadef.h"
    #include "init.h"
    #include <mpi.h>
    
    #define max(x,y) ((x)>(y)?(x):(y))
    #define min(x,y) ((x)<(y)?(x):(y))

    
    extern int *ileft, *iright,cell;
    extern int nprocs, proc;

    float inbuf, outbuf;
    int count;
    extern double TotalComuTime,AllreduceTime,GatherTime,BcastTime,RecComuTime,SendComuTime;

    extern double SStartTime,RStartTime,ComuStartTime,GatherStart,BcastStart,AllreduceStart;
    /* Computation of tentative velocity field (f, g) */
    void compute_tentative_velocity(float **u, float **v, float **f, float **g,
    char **flag, int imax, int jmax, float del_t, float delx, float dely,
    float gamma, float Re)
    {
    int  i, j;
    float du2dx, duvdy, duvdx, dv2dy, laplu, laplv;
    
    for (i=1; i<=imax-1; i++) {
        for (j=1; j<=jmax; j++) {
        /* only if both adjacent cells are fluid cells */
        if ((flag[i][j] & C_F) && (flag[i+1][j] & C_F)) {
            du2dx = ((u[i][j]+u[i+1][j])*(u[i][j]+u[i+1][j])+
            gamma*fabs(u[i][j]+u[i+1][j])*(u[i][j]-u[i+1][j])-
            (u[i-1][j]+u[i][j])*(u[i-1][j]+u[i][j])-
            gamma*fabs(u[i-1][j]+u[i][j])*(u[i-1][j]-u[i][j]))
            /(4.0*delx);
            duvdy = ((v[i][j]+v[i+1][j])*(u[i][j]+u[i][j+1])+
            gamma*fabs(v[i][j]+v[i+1][j])*(u[i][j]-u[i][j+1])-
            (v[i][j-1]+v[i+1][j-1])*(u[i][j-1]+u[i][j])-
            gamma*fabs(v[i][j-1]+v[i+1][j-1])*(u[i][j-1]-u[i][j]))
            /(4.0*dely);
            laplu = (u[i+1][j]-2.0*u[i][j]+u[i-1][j])/delx/delx+
            (u[i][j+1]-2.0*u[i][j]+u[i][j-1])/dely/dely;
    
            f[i][j] = u[i][j]+del_t*(laplu/Re-du2dx-duvdy);
        } else {
            f[i][j] = u[i][j];
        }
        }
    }
    
    for (i=1; i<=imax; i++) {
        for (j=1; j<=jmax-1; j++) {
        /* only if both adjacent cells are fluid cells */
        if ((flag[i][j] & C_F) && (flag[i][j+1] & C_F)) {
            duvdx = ((u[i][j]+u[i][j+1])*(v[i][j]+v[i+1][j])+
            gamma*fabs(u[i][j]+u[i][j+1])*(v[i][j]-v[i+1][j])-
            (u[i-1][j]+u[i-1][j+1])*(v[i-1][j]+v[i][j])-
            gamma*fabs(u[i-1][j]+u[i-1][j+1])*(v[i-1][j]-v[i][j]))
            /(4.0*delx);
            dv2dy = ((v[i][j]+v[i][j+1])*(v[i][j]+v[i][j+1])+
            gamma*fabs(v[i][j]+v[i][j+1])*(v[i][j]-v[i][j+1])-
            (v[i][j-1]+v[i][j])*(v[i][j-1]+v[i][j])-
            gamma*fabs(v[i][j-1]+v[i][j])*(v[i][j-1]-v[i][j]))
            /(4.0*dely);
    
            laplv = (v[i+1][j]-2.0*v[i][j]+v[i-1][j])/delx/delx+
            (v[i][j+1]-2.0*v[i][j]+v[i][j-1])/dely/dely;
    
            g[i][j] = v[i][j]+del_t*(laplv/Re-duvdx-dv2dy);
        } else {
            g[i][j] = v[i][j];
        }
        }
    }
    
    /* f & g at external boundaries */
    for (j=1; j<=jmax; j++) {
        f[0][j]    = u[0][j];
        f[imax][j] = u[imax][j];
    }
    for (i=1; i<=imax; i++) {
        g[i][0]    = v[i][0];
        g[i][jmax] = v[i][jmax];
    }
    }
    
    
    /* Calculate the right hand side of the pressure equation */
    void compute_rhs(float **f, float **g, float **rhs, char **flag, int imax,
    int jmax, float del_t, float delx, float dely)
    {
    int i, j;
    
    for (i=1;i<=imax;i++) {
        for (j=1;j<=jmax;j++) {
        if (flag[i][j] & C_F) {
            /* only for fluid and non-surface cells */
            rhs[i][j] = (
                (f[i][j]-f[i-1][j])/delx +
                (g[i][j]-g[i][j-1])/dely
                ) / del_t;
        }
        }
    }
    }
    
    
    /* Red/Black SOR to solve the poisson equation */
    int poisson(float **p, float **rhs, char **flag, int imax, int jmax,
    float delx, float dely, float eps, int itermax, float omega,
    float *res, int ifull)
    {
   int  i, j,iter;
     int pleft,lpleft,npleft, pright,lpright,npright;

    float add, beta_2, beta_mod;
    float p0 = 0.0;
    int pnext,pprev;
    pnext=proc+1;
    pprev=proc-1;
    int rb; /* Red-black value. */
    
    float rdx2 = 1.0/(delx*delx);
    float rdy2 = 1.0/(dely*dely);
    beta_2 = -omega/(2.0*(rdx2+rdy2));
    MPI_Status status;
    pleft=proc*(imax/nprocs)+1; /*the first column of the matrix of the current proc */
    npleft=(proc+1)*(imax/nprocs)+1;
    lpleft=(proc-1)*(imax/nprocs)+1;
    pright=(proc+1)*(imax/nprocs);/*the last column of the matrix of the current proc*/
    npright=(proc+2)*(imax/nprocs);
    lpright=(proc)*(imax/nprocs);
    /*if imax can not be divided evenly by the N processors, then the rest of the matrix belongs to the last processor */
     if(imax%nprocs!=0)
     {
        if(proc==nprocs-1)
        {pright =imax;
        }
     }

    /* Calculate sum of squares */

    for (i = 1; i <=imax; i++) {
        for (j=1; j<=jmax; j++) {
        if (flag[i][j] & C_F) { p0 += p[i][j]*p[i][j]; }
        }
    }
    
    p0 = sqrt(p0/ifull);
    if (p0 < 0.0001) { p0 = 1.0; }
    
    /* Red/Black SOR-iteration */

    for (iter = 0; iter < itermax; iter++) {
        for (rb = 0; rb <= 1; rb++) {
        /*Divide the matrix by column*/
        for (i = pleft; i <= pright; i++) {
            for (j = 1; j <= jmax; j++) {
            if ((i+j) % 2 != rb) { continue; }
            if (flag[i][j] == (C_F | B_NSEW)) {
                /* five point star for interior fluid cells */
                p[i][j] = (1.-omega)*p[i][j] - 
                beta_2*(
                    (p[i+1][j]+p[i-1][j])*rdx2
                    + (p[i][j+1]+p[i][j-1])*rdy2
                    -  rhs[i][j]
                );
            } else if (flag[i][j] & C_F) { 
                /* modified star near boundary */
                beta_mod = -omega/((eps_E+eps_W)*rdx2+(eps_N+eps_S)*rdy2);
                p[i][j] = (1.-omega)*p[i][j] -
                beta_mod*(
                    (eps_E*p[i+1][j]+eps_W*p[i-1][j])*rdx2
                    + (eps_N*p[i][j+1]+eps_S*p[i][j-1])*rdy2
                    - rhs[i][j]
                );
            }
            cell++;
            } /* end of j */
        } /* end of i */


       SStartTime=MPI_Wtime();
        /*if the current proccesor is not the last proccesor
           then send the last column of the current matrix to the
           first column of the next matrix of the proccesor */

     if(proc<nprocs-1){MPI_Send(&p[pright][0],jmax+2,MPI_FLOAT,pnext,99,MPI_COMM_WORLD); }

          /*if the current proccesor is not the first proccesor
        then receive the last column of the previous proccesor.*/

     if(proc>0){ MPI_Recv(&p[pleft-1][0],jmax+2,MPI_FLOAT,pprev,99,MPI_COMM_WORLD,&status);}
        /*at the same time the current processor will also send the first column of the matrix
         back to the previous processor  */
     if(proc>0) { MPI_Send(&p[pleft][0],jmax+2,MPI_FLOAT,pprev,99,MPI_COMM_WORLD);}

/*if the current proccesor is not the last proccesor, then receive the last column of the matrix of the next processor*/
    if(proc<nprocs-1){ MPI_Recv(&p[pright+1][0],jmax+2,MPI_FLOAT,pnext,99,MPI_COMM_WORLD,&status);}

      } /* end of rb */
      SendComuTime=SendComuTime+(MPI_Wtime()-SStartTime);
    /*
       }  end of rb
  / if (proc==0)
    {

    MPI_Send(p[pright],jmax+2,MPI_FLOAT,des,proc,MPI_COMM_WORLD);
    MPI_Recv(p[npleft],jmax+2,MPI_FLOAT,des,des,MPI_COMM_WORLD,&status);
    }
    if(proc!=0 &&proc!=nprocs-1){
    MPI_Recv(p[lpright],jmax+2,MPI_FLOAT,src,src,MPI_COMM_WORLD,&status);
    MPI_Send(p[pright],jmax+2,MPI_FLOAT,src,proc,MPI_COMM_WORLD);
    MPI_Send(p[pleft],jmax+2,MPI_FLOAT,des,proc,MPI_COMM_WORLD);
    MPI_Recv(p[npleft],jmax+2,MPI_FLOAT,des,des,MPI_COMM_WORLD,&status);


    }
    if (proc==nprocs-1)
    {
    MPI_Recv(p[lpright],jmax+2,MPI_FLOAT,src,src,MPI_COMM_WORLD,&status);
    MPI_Send(p[pleft],jmax+2,MPI_FLOAT,src,proc,MPI_COMM_WORLD);
    }   */




    /*if(imax%nprocs!=0 && proc==nprocs-1&&proc!=0)
    {
    pright+=imax%nprocs;
     MPI_Status status;
    if( proc==nprocs-1){
     MPI_Recv(p[pright],jmax+2,MPI_FLOAT,(proc-1+nprocs)/nprocs,100,MPI_COMM_WORLD,&status);
     MPI_Send(p[pleft],jmax+2,MPI_FLOAT,(proc+1)/nprocs,100,MPI_COMM_WORLD);}
    else
     {

     MPI_Send(p[pright],jmax+2,MPI_FLOAT,(proc+1)/nprocs,100,MPI_COMM_WORLD);
     MPI_Recv(p[pleft],jmax+2,MPI_FLOAT,(proc-1+nprocs)/nprocs,100,MPI_COMM_WORLD,&status);
     }
    }
   if(imax%nprocs==0)
    {
    if(proc!=0){
    MPI_Status status;
    MPI_Send(p[pright],jmax+2,MPI_FLOAT,(proc+1)/nprocs,99,MPI_COMM_WORLD);
    MPI_Recv(p[pleft],jmax+2,MPI_FLOAT,(proc-1+nprocs)/nprocs,99,MPI_COMM_WORLD,&status);
     }
     }


      */

        /* Partial computation of residual */
        *res = 0.0;


        for (i = pleft; i <= pright; i++) {
        for (j = 1; j <= jmax; j++) {
            if (flag[i][j] & C_F) {
            /* only fluid cells */
            add = (eps_E*(p[i+1][j]-p[i][j]) - 
                eps_W*(p[i][j]-p[i-1][j])) * rdx2  +
                (eps_N*(p[i][j+1]-p[i][j]) -
                eps_S*(p[i][j]-p[i][j-1])) * rdy2  -  rhs[i][j];
    
                  *res += add*add;
            }
        }
        }
          AllreduceStart=MPI_Wtime();

          MPI_Allreduce(res,res,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
          AllreduceTime= AllreduceTime+(MPI_Wtime()-AllreduceStart);


        *res = sqrt((*res)/ifull)/p0;
    
        /* convergence? */
        if (*res<eps) break;
    } /* end of iter */
    /*Gather all the data and put it into proccesor 0*/
    GatherStart=MPI_Wtime();
    MPI_Gather(&p[pleft][0],imax/nprocs*(jmax+2),MPI_FLOAT,&p[pleft][0],imax/nprocs*(jmax+2),MPI_FLOAT,0,MPI_COMM_WORLD);
    GatherTime= GatherTime+(MPI_Wtime()-GatherStart);
   /*MPI_Bcast(p[0],(pright-pleft+1)*(jmax+2),MPI_FLOAT,0, MPI_COMM_WORLD);  */

   /* int lproc,ir,il;
    if(proc==0){
    for (lproc = 1; lproc <= nprocs;lproc++){
    il=lproc*(imax/nprocs)+1;
    ir=(proc+1)*(imax/nprocs);


    MPI_Status status;
    MPI_Recv(p[ir],(ir-il+1)*(jmax+2),MPI_FLOAT,lproc,lproc+10,MPI_COMM_WORLD, &status);
    */


    /*if(lproc!=0)
    {MPI_Send(p[ir],(ir-il+1)*(jmax+2),MPI_FLOAT,0,lproc,MPI_COMM_WORLD);
    }
    }
    }*/

    if(imax%nprocs !=0)
    {   if(proc==nprocs-1)  { MPI_Send(&p[imax-imax%nprocs+1][0],imax%nprocs*(jmax+2),MPI_FLOAT,0,999,MPI_COMM_WORLD);}
        if(proc==0)  {MPI_Recv(&p[imax-imax%nprocs+1][0],imax%nprocs*(jmax+2),MPI_FLOAT,nprocs-1,999,MPI_COMM_WORLD,&status);}
    }
    /*the *u and *v will control by the matrix size of the *p so we need to broadcast all the data to every process*/
    BcastStart=MPI_Wtime();
    MPI_Bcast(&p[0][0],(imax+2)*(jmax+2), MPI_FLOAT, 0,MPI_COMM_WORLD);
    BcastTime= BcastTime+(MPI_Wtime()-BcastStart);

    return iter;


    }
    
    /* Update the velocity values based on the tentative
    * velocity values and the new pressure matrix
    */
    void update_velocity(float **u, float **v, float **f, float **g, float **p,
    char **flag, int imax, int jmax, float del_t, float delx, float dely)
    {
    int i, j;
    
    for (i=1; i<=imax-1; i++) {
        for (j=1; j<=jmax; j++) {
        /* only if both adjacent cells are fluid cells */
        if ((flag[i][j] & C_F) && (flag[i+1][j] & C_F)) {
            u[i][j] = f[i][j]-(p[i+1][j]-p[i][j])*del_t/delx;
        }
        }
    }
    for (i=1; i<=imax; i++) {
        for (j=1; j<=jmax-1; j++) {
        /* only if both adjacent cells are fluid cells */
        if ((flag[i][j] & C_F) && (flag[i][j+1] & C_F)) {
            v[i][j] = g[i][j]-(p[i][j+1]-p[i][j])*del_t/dely;
        }
        }
    }
    }
    
    
    /* Set the timestep size so that we satisfy the Courant-Friedrichs-Lewy
    * conditions (ie no particle moves more than one cell width in one
    * timestep). Otherwise the simulation becomes unstable.
    */
    void set_timestep_interval(float *del_t, int imax, int jmax, float delx,
    float dely, float **u, float **v, float Re, float tau)
    {
    int i, j;
    float umax, vmax, deltu, deltv, deltRe; 
    
    /* del_t satisfying CFL conditions */
    if (tau >= 1.0e-10) { /* else no time stepsize control */
        umax = 1.0e-10;
        vmax = 1.0e-10; 
        for (i=0; i<=imax+1; i++) {
        for (j=1; j<=jmax+1; j++) {
            umax = max(fabs(u[i][j]), umax);
        }
        }
        for (i=1; i<=imax+1; i++) {
        for (j=0; j<=jmax+1; j++) {
            vmax = max(fabs(v[i][j]), vmax);
        }
        }
    
        deltu = delx/umax;
        deltv = dely/vmax; 
        deltRe = 1/(1/(delx*delx)+1/(dely*dely))*Re/2.0;
    
        if (deltu<deltv) {
        *del_t = min(deltu, deltRe);
        } else {
        *del_t = min(deltv, deltRe);
        }
        *del_t = tau * (*del_t); /* multiply by safety factor */
    }
    }