mpi_d2q9.c

格子Boltzmann方法 格子Boltzmann方法是为了保留格子气自动机方法的优点

// 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)        // 临时使用的处理字符串函数
{    // 删除str中的空格并将小写英文字母转换成大写字母，并将等号转换成'\0'，返回其后面的下标值
    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++)        // 注意循环的顺序
        {
            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)    // 流场流体入口处
            {
                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;