// mpi_D2Q9.c#include "mpi.h"#include "mpi_D2Q9.h"#include "MyArrays.h"#include <stdio.h>#include <math.h>#include <sys/times.h>#include <sys/time.h>int D=2, Q=9, Nx, Ny, n;char e[9][2], e1[9], **boundary, bounceback, filename[80]="管道流.lbm";double ***f, ***f1, **rho, **u, **v, *f0;double mu, Re, omega, dx, dt, t, c, c_s;double Lx, Ly, L0, U0, V0, Rho0, rho_max;double CPU_time(){    #ifndef CLK_TCK        #define CLK_TCK 100.0    #endif    struct tms t;    times(&t);    return (1.0*t.tms_utime + t.tms_stime)/CLK_TCK;}void ShowBdy(int nx, int ny){    int i, j;    for(j=0; j<ny; j++)    {        for(i=0; i<nx; i++)            printf("%d", boundary[i][j]);        printf("\n");    }}int MyUpper(char *str)        // 臨時(shí)使用的處理字符串函數(shù){    // 刪除str中的空格并將小寫英文字母轉(zhuǎn)換成大寫字母，并將等號(hào)轉(zhuǎn)換成'\0'，返回其后面的下標(biāo)值    int i=0, j=0, k=0;    for(i=j=0; str[i]; i++)    {        if(str[i]!=' ')        {            if(str[i]>='a' && str[i]<='z')                str[j] = str[i] - 'a'+'A';            else if(str[i]=='=')            {                str[j] = '\0';                k = j + 1;            }            else                str[j] = str[i];            j++;        }    }    str[j] = '\0';    return k;}void ReadParaFile(char *filename){    int i, j;    char str[80];    FILE *fp;    if((fp=fopen(filename, "rt"))==NULL) exit(1);    mu = 1.0;    L0 = -1.0;    bounceback = 0;    while(1)    {        fgets(str, 80, fp);        i = MyUpper(str);//        if(strcmp(str, "D")==0)//            sscanf(str+i, "%d", &D);//        else if(strcmp(str, "Q")==0)//            sscanf(str+i, "%d", &Q);//        else         if(strcmp(str, "DOMAINSIZE")==0)            sscanf(str+i, "%lf,%lf", &Lx, &Ly);        else if(strcmp(str, "L0")==0)            sscanf(str+i, "%lf", &L0);        else if(strcmp(str, "BOUNCEBACK")==0)            sscanf(str+i, "%d", &bounceback);        else if(strcmp(str, "MESHSIZE")==0)            sscanf(str+i, "%d,%d", &Nx, &Ny);        else if(strcmp(str, "MU")==0)            sscanf(str+i, "%lf", &mu);        else if(strcmp(str, "RE")==0)            sscanf(str+i, "%lf", &Re);        else if(strcmp(str, "OMEGA")==0)            sscanf(str+i, "%lf", &omega);        else if(strcmp(str, "BOUNDARYDATA")==0)            break;    }    if(L0<0) L0 = Ly;    boundary=(char**)Calloc2DArray(Nx, Ny, sizeof(char));    for(j=0; j<Ny; j++)    {        for(i=0; i<Nx; i++)        // 注意循環(huán)的順序        {            fscanf(fp, "%c", str);            boundary[i][j] = str[0] - '0';            }        fgets(str, 80, fp);    }    fclose(fp);    ShowBdy(Nx, Ny);}void Finish(){    if(boundary) Free2DArray((void**)boundary);    if(f0)  Free1DArray((void*)f0);    if(u)   Free2DArray((void**)v);    if(v)   Free2DArray((void**)u);    if(rho) Free2DArray((void**)rho);    if(f1)  Free3DArray((void***)f1);    if(f)   Free3DArray((void***)f);}void InitData(int numprocs, int myid){    int i, nx, packsize, position, extra;    char packbuf[100];    MPI_Status status;    e[0][0] = 0;  e[0][1] = 0;    e[1][0] = 1;  e[1][1] = 0;    e[2][0] = 0;  e[2][1] = 1;    e[3][0] =-1;  e[3][1] = 0;    e[4][0] = 0;  e[4][1] =-1;    e[5][0] = 1;  e[5][1] = 1;    e[6][0] =-1;  e[6][1] = 1;    e[7][0] =-1;  e[7][1] =-1;    e[8][0] = 1;  e[8][1] =-1;    e1[0] = 0;    e1[1] = 3;    e1[2] = 4;    e1[3] = 1;    e1[4] = 2;    e1[5] = 7;    e1[6] = 8;    e1[7] = 5;    e1[8] = 6;    if(myid==0)    {        packsize = 0;        MPI_Pack(&Nx, 1, MPI_INT, packbuf, 100, &packsize, MPI_COMM_WORLD);        MPI_Pack(&Ny, 1, MPI_INT, packbuf, 100, &packsize, MPI_COMM_WORLD);        MPI_Pack(&bounceback, 1, MPI_CHAR, packbuf, 100, &packsize, MPI_COMM_WORLD);        MPI_Pack(&Lx, 1, MPI_DOUBLE, packbuf, 100, &packsize, MPI_COMM_WORLD);        MPI_Pack(&Ly, 1, MPI_DOUBLE, packbuf, 100, &packsize, MPI_COMM_WORLD);        MPI_Pack(&mu, 1, MPI_DOUBLE, packbuf, 100, &packsize, MPI_COMM_WORLD);        MPI_Pack(&Re, 1, MPI_DOUBLE, packbuf, 100, &packsize, MPI_COMM_WORLD);        MPI_Pack(&L0, 1, MPI_DOUBLE, packbuf, 100, &packsize, MPI_COMM_WORLD);        MPI_Pack(&omega, 1, MPI_DOUBLE, packbuf, 100, &packsize, MPI_COMM_WORLD);    }    MPI_Bcast(&packsize, 1, MPI_INT, 0, MPI_COMM_WORLD);    MPI_Bcast(packbuf, packsize, MPI_PACKED, 0, MPI_COMM_WORLD);    if(myid)    {        position = 0;        MPI_Unpack(packbuf, packsize, &position, &Nx, 1, MPI_INT, MPI_COMM_WORLD);        MPI_Unpack(packbuf, packsize, &position, &Ny, 1, MPI_INT, MPI_COMM_WORLD);        MPI_Unpack(packbuf, packsize, &position, &bounceback, 1, MPI_CHAR, MPI_COMM_WORLD);        MPI_Unpack(packbuf, packsize, &position, &Lx, 1, MPI_DOUBLE, MPI_COMM_WORLD);        MPI_Unpack(packbuf, packsize, &position, &Ly, 1, MPI_DOUBLE, MPI_COMM_WORLD);        MPI_Unpack(packbuf, packsize, &position, &mu, 1, MPI_DOUBLE, MPI_COMM_WORLD);        MPI_Unpack(packbuf, packsize, &position, &Re, 1, MPI_DOUBLE, MPI_COMM_WORLD);        MPI_Unpack(packbuf, packsize, &position, &L0, 1, MPI_DOUBLE, MPI_COMM_WORLD);        MPI_Unpack(packbuf, packsize, &position, &omega, 1, MPI_DOUBLE, MPI_COMM_WORLD);    }    nx = Nx/numprocs;    extra = (myid==0 || myid==numprocs-1) ? 1 : 2;    D = 2;    Q = 9;    if(myid==0)    {        for(i=1; i<numprocs-1; i++)            MPI_Send(boundary[i*nx-1], (nx+2)*Ny, MPI_CHAR, i, 100*i, MPI_COMM_WORLD);        MPI_Send(boundary[i*nx-1], (nx+1)*Ny, MPI_CHAR, i, 100*i, MPI_COMM_WORLD);        rho = (double **)Calloc2DArray(Nx, Ny, sizeof(double));        u   = (double **)Calloc2DArray(Nx, Ny, sizeof(double));        v   = (double **)Calloc2DArray(Nx, Ny, sizeof(double));    }    else    {        boundary=(char**)Calloc2DArray(nx+extra, Ny, sizeof(char));        MPI_Recv(boundary[0], (nx+extra)*Ny, MPI_CHAR, 0, 100*myid, MPI_COMM_WORLD, &status);        rho = (double **)Calloc2DArray(nx+extra, Ny, sizeof(double));        u   = (double **)Calloc2DArray(nx+extra, Ny, sizeof(double));        v   = (double **)Calloc2DArray(nx+extra, Ny, sizeof(double));    }    f   = (double***)Calloc3DArray(nx+extra, Ny, Q, sizeof(double));    f1  = (double***)Calloc3DArray(nx+extra, Ny, Q, sizeof(double));    f0  = (double  *)Calloc1DArray(Q, sizeof(double));    dx = Lx/Nx;    dt = (2/omega-1)*dx*dx/(6*mu);    c = dx/dt;    c_s = c/sqrt(3);    U0 = mu*Re/L0;    V0 = 0;    Rho0 = 1;    t = 0;    n = 0;    init_macro(numprocs, myid);    init_micro(numprocs, myid);}void init_macro(int numprocs, int myid){    int i, j, nx;    if(myid==0 || myid==numprocs-1)        nx = Nx/numprocs + 1;    else        nx = Nx/numprocs + 2;    for(i=0; i<nx; i++)    {        for(j=0; j<Ny; j++)        {            if(boundary[i][j]==1)    // 流場(chǎng)流體入口處            {                u[i][j] = U0;                v[i][j] = V0;            }            else                    // 其它                u[i][j] = v[i][j] = 0;            rho[i][j] = Rho0;        }    }}void init_micro(int numprocs, int myid){    int i, j, nx;    if(myid==0 || myid==numprocs-1)        nx = Nx/numprocs + 1;    else        nx = Nx/numprocs + 2;    for(i=0; i<nx; i++)        for(j=0; j<Ny; j++)            feq(u[i][j], v[i][j], rho[i][j], f[i][j]);}void CollectData(int numprocs, int myid){    int i, nx;    MPI_Status status;    nx = Nx/numprocs;    if(myid==0)    {        for(i=1; i<numprocs; i++)        {            MPI_Recv(u[i*nx], nx*Ny, MPI_DOUBLE, i, 1000*i, MPI_COMM_WORLD, &status);            MPI_Recv(v[i*nx], nx*Ny, MPI_DOUBLE, i, 1000*i, MPI_COMM_WORLD, &status);            MPI_Recv(rho[i*nx], nx*Ny, MPI_DOUBLE, i, 1000*i, MPI_COMM_WORLD, &status);        }    }    else    {        MPI_Send(u[1], nx*Ny, MPI_DOUBLE, 0, 1000*myid, MPI_COMM_WORLD);        MPI_Send(v[1], nx*Ny, MPI_DOUBLE, 0, 1000*myid, MPI_COMM_WORLD);        MPI_Send(rho[1], nx*Ny, MPI_DOUBLE, 0, 1000*myid, MPI_COMM_WORLD);    }}void feq(double u, double v, double rho, double *f0){/*    double u2, eu, eu2;    int k;