simulation.c.bak

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;
extern int nprocs, proc;


extern double befcommusrtime, aftcommusrtime,srcommutime, befcommualtime ,aftcommualtime ,befgathtime , aftgathtime ,gathtime ,befbcasttime,aftbcasttime,bcasttime, befundicomtime,aftundicomtime.undicommutime,commualtime, totalcommutime, totalcomputtime;

/* 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;
    float add, beta_2, beta_mod;
    float p0 = 0.0;
    
    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;
	/* the value of proce is start from 0*/
	int right = imax/nprocs*(proc+1);
		if(imax%nprocs!=0)
		   if(proc==nprocs-1)
			right =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++) {


      /*      for (i = 1; i <= imax; i++) { */
	          for (i = 1+imax/nprocs*proc; i <= right; 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]
                            );
                    }
                } /* end of j */
            } /* end of i */

			befcommusrtime = MPI_Wtime();

			/*communication of data with each proc*/
			if(proc<nprocs-1){
		/*sent data to next proc. because the row has been divided, then the number of value +2 
			the first proc is 0, so here proc should +1*/
			

			MPI_Send(&p[imax/nprocs*(proc+1)][0],jmax+2,MPI_FLOAT,proc+1,111,MPI_COMM_WORLD);
			
					}
		/*recieve data to last proc. because the row has been divided, then the number of value +2  */	
			if(proc>0){
		
			MPI_Recv(&p[imax/nprocs*proc][0],jmax+2,MPI_FLOAT,proc-1,111,MPI_COMM_WORLD,&status);
		
				}
		/*at the same time this proc should sent the boundy data to back as well*/
		    if(proc>0) {  
			
		    MPI_Send(&p[1+imax/nprocs*proc][0],jmax+2,MPI_FLOAT,proc-1,111,MPI_COMM_WORLD);
		
				}
			/*the last one should recieve the data back as well*/
			if(proc<nprocs-1){
			
			MPI_Recv(&p[1+imax/nprocs*(proc+1)][0],jmax+2,MPI_FLOAT,proc+1,111,MPI_COMM_WORLD,&status);
		
					}

			aftcommusrtime= MPI_Wtime();
			srcommutime = srcommutime + (aftcommusrtime- befcommusrtime) ;
			




        } /* end of rb */


        
        /* Partial computation of residual */
        *res = 0.0;
        for (i = 1+imax/nprocs*proc; i <= imax/nprocs*(proc+1); 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;
                }
            }
        }
		
		befcommualtime = MPI_Wtime();
		MPI_Allreduce(res,res,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
	    aftcommualtime = MPI_Wtime();
		commualtime= commualtime + (aftcommualtime -befcommualtime);

        *res = sqrt((*res)/ifull)/p0;

        /* convergence? */
        if (*res<eps) break;
    } /* end of iter */

 	/*get all the data together and put them in to root 0*/
	
	befgathtime = MPI_Wtime();
	MPI_Gather(&p[1+imax/nprocs*proc][0],imax/nprocs*(jmax+2),MPI_FLOAT,&p[1+imax/nprocs*proc][0],imax/nprocs*(jmax+2),MPI_FLOAT,0,MPI_COMM_WORLD);
	aftgathtime = MPI_Wtime();
	gathtime = gathtime + (aftgathtime-befgathtime);

	
	befundicomtime= MPI_Wtime();
	if(imax%nprocs !=0)
	{
		if(proc==nprocs-1)
			MPI_Send(&p[imax-imax%nprocs+1][0],imax%nprocs*(jmax+2),MPI_FLOAT,0,333,MPI_COMM_WORLD);
		if(proc==0)
			MPI_Recv(&p[imax-imax%nprocs+1][0],imax%nprocs*(jmax+2),MPI_FLOAT,nprocs-1,333,MPI_COMM_WORLD,&status);
	}
	aftundicomtime = MPI_Wtime();
	undicommutime = undicommutime + (aftundicomtime- befundicomtime);


		/*broadcast the data in root 0 to all the process*/

			befbcasttime = MPI_Wtime();
			MPI_Bcast(&p[0][0],(imax+2)*(jmax+2),MPI_FLOAT,0,MPI_COMM_WORLD);
			aftbcasttime = MPI_Wtime();
			bcasttime=bcasttime + (aftbcasttime-befbcasttime);
			totalcommutime= totalcommutime + srcommutime + commualtime + gathtime + undicommutime + bcasttime;
			

 
    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 */
    }
}