mpi_d2q9.c

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

    u /= c;
    v /= c;
    u2 = 1-1.5*(u*u + v*v);

    f0[0] = 4.0/9*rho*u2;
    for(k=1; k<Q; k++)
    {
        eu = e[k][0]*u + e[k][1]*v;
        eu2 = eu*eu;
        f0[k] = rho/9*(3*eu + 4.5*eu2 + u2);
    }
    for(k=5; k<Q; k++)
        f0[k] *= 0.25;
*/
    double temp_u2, temp_rho, uav, vmu;

    u /= c;
    v /= c;
    temp_u2 = 1-1.5*(u*u + v*v);
    temp_rho = 4.0*rho/9.0;
    uav = u+v;
    vmu = v-u;

    f0[0] = temp_rho * temp_u2;

    temp_rho *= 0.25;
    f0[1] = temp_rho * (temp_u2 + 3*u + 4.5*u*u);
    f0[2] = temp_rho * (temp_u2 + 3*v + 4.5*v*v);
    f0[3] = f0[1] - temp_rho * 6*u;
    f0[4] = f0[2] - temp_rho * 6*v;

    temp_rho *= 0.25;
    f0[5] = temp_rho * (temp_u2 + 3*uav + 4.5*uav*uav);
    f0[6] = temp_rho * (temp_u2 + 3*vmu + 4.5*vmu*vmu);
    f0[7] = f0[5] - temp_rho * 6*uav;
    f0[8] = f0[6] - temp_rho * 6*vmu;
}

void collision(int numprocs, int myid)
{
    int i, j, k, nx, x0, x1;

    nx = Nx/numprocs;
    if(myid==0)
    {
        x0 = 0; x1 = nx-1;
    }
    else
    {
        x0 = 1; x1 = nx;
    }
    for(i=x0; i<=x1; i++)
    {
        for(j=0; j<Ny; j++)
        {
            if(boundary[i][j]==1)
            {
                u[i][j] = U0;
                v[i][j] = V0;
                rho[i][j] = Rho0;
                continue;
            }
            else if(boundary[i][j]==3)
            {
                u[i][j] = 0;
                v[i][j] = 0;
                continue;
            }

            rho[i][j] = f[i][j][0];
            u[i][j] = v[i][j] = 0;
            for(k=1; k<Q; k++)
            {
                u[i][j] += e[k][0] * f[i][j][k];
                v[i][j] += e[k][1] * f[i][j][k];
                rho[i][j] += f[i][j][k];
            }
            u[i][j] *= c;
            v[i][j] *= c;
            if(rho[i][j]>rho_max) rho_max = rho[i][j];
            if(rho[i][j]>1e-6)
            {
                u[i][j] /= rho[i][j];
                v[i][j] /= rho[i][j];
            }
            else
            {
                u[i][j] = 0;
                v[i][j] = 0;
            }
        }
    }                        // 计算宏观量

    for(i=x0; i<=x1; i++)
    {
        for(j=0; j<Ny; j++)
        {
            feq(u[i][j], v[i][j], rho[i][j], f0);
            for(k=0; k<Q; k++)
                f1[i][j][k] = (1-omega)*f[i][j][k] + omega*f0[k];
        }
    }
}

void Jacobi(int numprocs, int myid)
{
    int nx=Nx/numprocs;
    MPI_Status status;

    if(myid==0)
    {
        MPI_Send(f1[nx-1][0], Ny*Q, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
        MPI_Recv(f1[nx][0],   Ny*Q, MPI_DOUBLE, 1, 500, MPI_COMM_WORLD, &status);
    }
    else if(myid==numprocs-1)
    {
        MPI_Send(f1[1][0], Ny*Q, MPI_DOUBLE, myid-1, 500*myid, MPI_COMM_WORLD);
        MPI_Recv(f1[0][0], Ny*Q, MPI_DOUBLE, myid-1, 500*(myid-1), MPI_COMM_WORLD, &status);
    }
    else
    {
        MPI_Send(f1[1][0],    Ny*Q, MPI_DOUBLE, myid-1, 500*myid, MPI_COMM_WORLD);
        MPI_Send(f1[nx][0],   Ny*Q, MPI_DOUBLE, myid+1, 500*myid, MPI_COMM_WORLD);
        MPI_Recv(f1[0][0],    Ny*Q, MPI_DOUBLE, myid-1, 500*(myid-1), MPI_COMM_WORLD, &status);
        MPI_Recv(f1[nx+1][0], Ny*Q, MPI_DOUBLE, myid+1, 500*(myid+1), MPI_COMM_WORLD, &status);
    }
}

void flow(int numprocs, int myid)
{
    int i, j, k, i1, j1, nx, x0, x1;
    double f00[9];

    nx = Nx/numprocs;
    if(myid==0)
    {
        x0 = 0; x1 = nx-1;
    }
    else
    {
        x0 = 1; x1 = nx;
    }

    if(bounceback)    // 流动步，反弹边界处理
    {
        for(i=x0; i<=x1; i++)
        {
            for(j=0; j<Ny; j++)
            {
                switch(boundary[i][j])
                {
                case 1:
                    feq(u[i][j], v[i][j], rho[i][j], f1[i][j]); break;
                case 2:
                    for(k=0; k<Q; k++)
                        f1[i][j][k] = f1[i-1][j][k];            break;
                case 3:
                    for(k=0; k<Q; k++)
                    {
                        i1 = i - e[k][0];
                        j1 = j - e[k][1];
                        if(i1>=0 && i1<=x1 && j1>=0 && j1<Ny && boundary[i1][j1]==0)
                            f1[i][j][k] = f1[i1][j1][e1[k]];
                    }                                            break;
                }
            }
        }

        for(i=x0; i<=x1; i++)
        {
            for(j=0; j<Ny; j++)
            {
                if(boundary[i][j]==0)
                {
                    for(k=0; k<Q; k++)
                        f[i][j][k] = f1[i-e[k][0]][j-e[k][1]][k];
                }
            }
        }
    }
    else        // 流动步，外推边界处理
    {
        for(i=x0; i<=x1; i++)
        {
            for(j=0; j<Ny; j++)
            {
                switch(boundary[i][j])
                {
                case 0:
                    for(k=0; k<Q; k++)
                    {
                        i1 = i - e[k][0];
                        j1 = j - e[k][1];
                        if(i1<0 || i1>x1 || j1<0 || j1>=Ny)
                        {
                            i1 = i + e[k][0];
                            j1 = j + e[k][1];
                            if(i1<0 || i1>x1 || j1<0 || j1>=Ny)    // 角点
                            {
                                f[i][j][k] = f1[i][j][k];
                            }
                            else    // 边点（外推）
                            {
                                feq(u[i][j], v[i][j], rho[i][j], f0);
                                feq(u[i1][j1], v[i1][j1], rho[i1][j1], f00);
                                f[i][j][k] = f1[i1][j1][k] + f0[k] - f00[k];
                            }
                        }
                        else
                            f[i][j][k] = f1[i1][j1][k];
                    }
                    break;
                case 1:
                    feq(u[i][j], v[i][j], rho[i][j], f[i][j]);
                    break;
                case 2:
                    for(k=0; k<Q; k++)
                        f[i][j][k] = f1[i-1][j][k];
                    break;                        // 出口处理:无穷远边界
                default:
                    feq(u[i][j], v[i][j], rho[i][j], f[i][j]);
                    break;
                }
            }
        }
    }
}

void evolution(int numprocs, int myid)
{
    collision(numprocs, myid);
    Jacobi(numprocs, myid);
    flow(numprocs, myid);
    t += dt;
    n++;
}

void SaveFile()    // 将速度(u,v)、密度 rho 写入相应的文件
{
    FILE *fpu, *fpv, *fprho;
    int nx=Nx, ny=Ny;
    int i, j;

    fpu = fopen("u.dat", "wb");
    fpv = fopen("v.dat", "wb");
    fprho = fopen("rho.dat", "wb");

    fwrite(&nx, 1, sizeof(int), fpu);    // 按二进制格式存放（前2个整数皆为网格规模）
    fwrite(&ny, 1, sizeof(int), fpu);
    fwrite(u[0], nx*ny, sizeof(u[0][0]), fpu);

    fwrite(&nx, 1, sizeof(int), fpv);
    fwrite(&ny, 1, sizeof(int), fpv);
    fwrite(v[0], nx*ny, sizeof(v[0][0]), fpv);

    fwrite(&nx, 1, sizeof(int), fprho);
    fwrite(&ny, 1, sizeof(int), fprho);
    fwrite(rho[0], nx*ny, sizeof(rho[0][0]), fprho);

    fclose(fpu); fclose(fpv); fclose(fprho);

    fpu = fopen("u.txt", "wb");
    fpv = fopen("v.txt", "wb");
    fprho = fopen("rho.txt", "wb");

    for(j=0; j<Ny; j++)        // 写成文本文件格式
    {
        for(i=0; i<Nx; i++)
        {
            fprintf(fpu, "%f ", u[i][j]);
            fprintf(fpv, "%f ", v[i][j]);
            fprintf(fprho, "%f ", rho[i][j]);
        }
        fprintf(fpu, "\n");
        fprintf(fpv, "\n");
        fprintf(fprho, "\n");
    }
    fclose(fpu); fclose(fpv); fclose(fprho);
}