为什么我的 c++ 代码会进入无限循环

我为 LBM 编写了一个C++代码来求解扩散方程，并在 Ubuntu 中用 g++ 编译了它。它编译得很好，但是当我尝试执行程序时，它并没有结束。我以为该程序以某种方式有一个无限循环，但我无法解决在哪里？提前谢谢你们。

这是我的代码：

#include <iostream>
#include <fstream>
using namespace std;
int main(){
    const int Lx = 100;
    const int Ly = 100;
    int sMod = 8; //speed model 
    float f[sMod][Lx+1][Ly+1];
    float rho[Lx+1][Ly+1];
    float feq, sum;
    float x[Lx+1];
    float y[Ly+1];
    float csq, alpha, omega,Tw;
    float w[9] = {4./9., 1./9., 1./9., 1./9., 1./9., 1./36., 1./36., 1./36., 1./36.};
    int i, j, dt, dx, k, dy; 
    FILE * mFile;
    dt = 1.0;
    dx = 1.0;
    dy = dx;
    x[0] = 0.0;
    for(i = 1; i < Lx; i++){
        x[i] = x[i-1] + dx;
        }
    y[0] = 0.0;
    for(j = 1; j < Ly; j++){
        y[j] = y[j-1] + dy;
        }
    Tw = 1.0;
    csq = (dx * dx)/(dt * dt);
    alpha = 0.25;
    omega = 1.0/((3.*alpha/(dt*csq))+0.5);
    cout << "csq: " << csq << ", omega: " << omega << endl;
    int mstep = 2; 
    for(j = 0; j <= Ly; j++)
    {
        for(i = 0; i <= Lx; i++)
        {
        rho[i][j] = 0.0;
        //cout << "1: " << rho[0][Ly/2] << endl;
        }
    }
        //cout << "5: " << rho[0][Ly/2] << endl;
    for(j = 0; j <= Ly; j++)
    {   
            //cout << "2: " << rho[0][Ly/2] << endl;
        for(i = 0; i <= Lx; i++)
        {
            //cout << "3: " << rho[0][Ly/2] << endl;
            for(k = 0; k <= sMod; k++)
            {
                            //cout << "4: " << rho[0][Ly/2] << endl;
                f[k][i][j] = w[k] * rho[i][j];
                if(i == 0){
                    f[k][i][j] = w[k] * Tw;
                    }
            }
        }
    }
            cout << "2: " << rho[0][Ly/2] << endl;
    for(k = 1; k <= mstep; k++)
    {
        for(j = 0; j <= Ly; j++)
        {
            for(i = 0; i <= Lx; i++)
            {   
                sum = 0.0;
                for(k = 0;k <= sMod; k++)
                {
                    sum = sum + f[k][i][j];
                }
                rho[i][j] = sum;
            }
        }
        cout << rho[0][Ly/2] << endl;
    for(j = 0; j <= Ly; j++)
    {   
        for(i = 0; i <= Lx; i++)
        {
            for(k = 0;k <= sMod; k++)
            {
                feq = w[k] * rho[i][j];
                f[k][i][j] = omega * feq + (1. - omega) * f[k][i][j];
            }
        }
    }
    /*-------
    streaming
    --------*/
    for(j = Ly; i >= 0; j--)
    {
        for(i = 0; i <= Lx; i++)
        {       
            f[2][i][j] = f[2][i][j-1];
            f[6][i][j] = f[6][i+1][j-1];
        }
    }
    for(j = Ly; i >= 0; j--)
    {
        for(i = Lx; i >= 0; i--)
        {                       
            f[1][i][j] = f[1][i-1][j]; //right to left
            f[5][i][j] = f[5][i-1][j-1];
        }
    }

    for(j = 0; i <= Ly; j--)
    {
        for(i = Lx; i >= 0; i--)
        {   
            f[4][i][j] = f[4][i][j+1]; 
            f[8][i][j] = f[8][i-1][j+1];
        }
    }
    for(j = 0; j <= Ly; j++)
    {   
        for(i = 0; i <= Lx; i++)
        {   
            f[3][i][j] = f[3][i+1][j]; 
            f[7][i][j] = f[7][i+1][j+1];
        }
    }
    //boundary conditions
    for(j = 0; j <= Ly; j++)
    {
        f[1][0][j] = w[1]*Tw + w[3]*Tw - f[3][0][j];
        f[5][0][j] = w[5]*Tw + w[7]*Tw - f[7][0][j];
        f[8][0][j] = w[8]*Tw + w[6]*Tw - f[6][0][j];
        f[3][Lx][j] = -f[1][Lx][j];
        f[6][Lx][j] = -f[8][Lx][j];
        f[7][Lx][j] = -f[5][Lx][j];
    }
    for(i = 0; i <= Lx; i++)
    {
        f[4][i][Ly] = -f[2][i][Ly];
        f[7][i][Ly] = -f[5][i][Ly];
        f[8][i][Ly] = -f[6][i][Ly];
        f[1][i][0] = f[1][i][1];
        f[2][i][0] = f[2][i][1];
        f[3][i][0] = f[3][i][1];
        f[4][i][0] = f[4][i][1];
        f[5][i][0] = f[5][i][1];
        f[6][i][0] = f[6][i][1];
        f[7][i][0] = f[7][i][1];
        f[8][i][0] = f[8][i][1];
    }
    }
    for(j = 0; j <= Ly; j++)
    {   
        for(i = 0; i <= Lx; i++)
        {   
            sum = 0.0;
            for(k = 0;k <= sMod; k++)
            {
                sum = sum + f[k][i][j];
            }
        rho[i][j] = sum;
        }
    }
    cout << rho[0][Ly/2] << endl;
    mFile = fopen("lbmdiffusiond2q9.csv","w");
    fprintf(mFile,"tt----TITLE= D2Q9 RESULTS----n");
    fprintf(mFile," VARIABLES = RHOn");
    //myfile << "ZONE " << "I=" << Lx+1 << " J=" << Ly+1 << " F=POINT" << endl;
    for(i = 0; i <= Lx; i++)
    {
        for(j = 0; j <= Ly; j++)
        {   
    //myfile << (dx/Lx) * i << " " << (dy/Ly) * j << " " << rho[i][j] << endl;
            fprintf(mFile,"%.6f ",rho[i][j]);   
        }
        fprintf(mFile,"n");
    }   
    fclose(mFile);      
return 0;
}