网格问题代码

#include <iostream>                                      //声明头文件
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#include <ctime>

using namespace std;

//D3Q19 3维19个离散速度
const int Q = 19;
int e[Q][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { -1, 0, 0 }, { 0, 1, 0 }, { 0, -1, 0 }, { 0, 0, 1 }, { 0, 0, -1 }, { 1, 1, 0 }, { -1, -1, 0 }, { 1, -1, 0 },
{ -1, 1, 0 }, { 1, 0, 1 }, { -1, 0, -1 }, { 1, 0, -1 }, { -1, 0, 1 }, { 0, 1, 1 }, { 0, -1, -1 }, { 0, 1, -1 }, { 0, -1, 1 } };
double w[Q] = { 1. / 3, 1. / 18, 1. / 18, 1. / 18, 1. / 18, 1. / 18, 1. / 18, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36 };

//定义物理尺寸
const int NX = 100;
const int NY = 100;
const int NZ = 100;
const double M = NX*NY*NZ;                                    //10*10*10的立方体为网格区域

//lattice boltzmann method
double rho[NX + 1][NY + 1][NZ + 1], u[NX + 1][NY + 1][NZ + 1][3], u0[NX + 1][NX + 1][NZ + 1][3];     //声明密度，演化后速度，演化前速度
double velocity;                                                                                     //定义速度
double f[NX + 1][NY + 1][NZ + 1][Q], F[NX + 1][NY + 1][NZ + 1][Q], feqq[NX + 1][NY + 1][NZ + 1][Q], ftemp[NX + 1][NY + 1][NZ + 1][Q];//声明速度分布函数演化前与演化后的，平衡态分布函数，临时分布函数
int i, j, k, q, ip, jp, kp, n;                                                                       //
double c, Re, dx, dy, dz, LX, LY, LZ, dt, rho0, tau_f, niu, omega, error, temp;//声明格子速度、雷诺数、空间步长、方向长度、时间步长、初始密度、弛豫时间、运动粘度、离散速度比、临时变量
//QuartetStructureGenerationSet
double arrgrid[NX + 1][NY + 1][NZ + 1], solidBoundaryNode[NX + 1][NY + 1][NZ + 1];               //定义格点  以及边界点（三维）

double rho_in = 1.0;                                                             //内部密度
double rho_out = 0.5;                                                            //外部密度

void initilization();                                                          //声明初始化子函数
void QuartetStructureGenerationSet();                                          //声明QSGS
double feq(int q, double rho, double u[3]);                                   //
void Evolution();                                                              //
void outputdata(int m);                                                        //声明输出
void Error();
void boundary();                                                               //声明边界条件
void Macroscopic();
void output3d();

int main()
{
 initilization();
 QuartetStructureGenerationSet();
 output3d();
 for (n = 0;; n++)
 {
  Evolution();
  Macroscopic();
  boundary();

  if (n % 100 == 0)
  {
   Error();
   cout << "The " << n << " th computation result: The max relative error of uvw is: "
    << setiosflags(ios::scientific) << error << endl;
   if (n >= 500)
   {
    if (n % 500 == 0)
     outputdata(n);
    if (error<1.0e-8)
     break;
   }
  }
 }
 return 0;
}

void initilization()                                                //定义初始化子函数
{
 dx = 1.0;                                                   //空间步长、时间步长均为1
 dy = 1.0;
 dz = 1.0;
 dt = dx;
 c = dx / dt;
 LX = dx*NX;                                                 //空间方向长度
 LY = dy*NY;
 LZ = dz*NZ;
 tau_f = 1.35;                                               //弛豫时间定为1.35？（tau_f = 3*niu+0.5，niu未知）
 omega = 1.0 / tau_f;                                        //书P48 式子3.103
 rho0 = 0.5;                                                 //初始密度为0.5

 for (k = 0; k <= NZ; k++)
 {
  for (j = 0; j <= NY; j++)
  {
   for (i = 0; i <= NX; i++)                   //初始化三维速度、密度
   {
    u[i][j][k][0] = 0.0;
    u[i][j][k][1] = 0.0;
    u[i][j][k][2] = 0.0;
    rho[i][j][k] = rho0;

    for (q = 0; q < Q; q++)            //格点个数以及分布函数
    {
     f[i][j][k][q] = feq(q, rho[i][j][k], u[i][j][k]);
    }
   }
  }
 }
}

void QuartetStructureGenerationSet()
{
 double d1_6 = 0.02;                                                      //定义六个主方向（接触）生长概率为0.02
 double d7_18 = d1_6 / 2.;                                                //定义十二个次方向（水平）生长概率为0.01
 double d19_26 = d1_6 / 8.;                                               //定义六个第三方向（对角线）生长概率为0.0025
 double n = 0.8;                                                         //定义孔隙率为0.8
 double cdd = 0.001;                                                     //定义生长核分布概率为0.001
 double numtotal_need = (1 - n)*NX*NY*NZ;                                //定义生长相的数量
 int numsoild = 0;                                                       //定义初始生长核数量为0
 int Tnumsoild;                                                          //定义整型生长核总量
 double soild[M+1][3];
 srand((unsigned)time(NULL));

 //第一步，遍历所有网格，生成固相内核
 for (i = 0; i <= NX; i++)                                               //生成格点[0][0][0],,,,[0][0][10];[0][1][0],,,[0][1][10];;;;[10][10][10]
 {
  for (j = 0; j <= NY; j++)
  {
   for (k = 0; k <= NZ; k++)
   {
    arrgrid[i][j][k] = 0;                           //初始化每一个格点，初始化后，格点的坐标变为[0][0][0]，并赋值为0
    soild[numsoild][0] = 0;
    soild[numsoild][1] = 0;
    soild[numsoild][2] = 0;

    if (((rand() % 1000) / 1000.0) < cd