Nov 27, 2014 at 12:08pm UTC
Hey guys,
I am having this error and I can not solve it. I am new to C++ and I need to solve this it's about a project.
Here is my code:
#include <iostream>
#include <fstream>
using namespace std;
int main(){
//lattice numbers in x, y direction
int Lx = 100;
int Ly = 100;
const int sMod = 8; //speed model
//declaration of distribution functions, density and domain
float f[sMod][Lx][Ly];
float rho[Lx][Ly];
float feq[9];
float x[Lx];
float y[Ly];
float csq, sum, alpha, omega, Tw, u, v, ck, density, T;
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;
}
//cout << x[99] << endl;
u = 0.1;
v = 0.4;
Tw = 1.0;
ck = dx/dt;
csq = ck * ck;
alpha = 1.0;
omega = 1.0/((3.*alpha/(dt*csq))+0.5);
cout << "csq: " << csq << ", omega: " << omega << endl;
int mstep = 400; //the total time steps number
density = 0.0;
for(j = 0; j < Ly; j++)
{
for(i = 0; i < Lx; i++)
{
for(k = 0; k < 9; k++)
{
f[k][i][j] = w[k] * density;
if(i == 0){
f[k][i][j] = w[k] * Tw;
}
}
}
}
cout << 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++)
{
feq[0] = w[0] * rho[i][j];
feq[1] = w[1] * rho[i][j] * (1.0 + 3.0 * (u/ck));
feq[2] = w[2] * rho[i][j] * (1.0 + 3.0 * (v/ck));
feq[3] = w[3] * rho[i][j] * (1.0 - 3.0 * (u/ck));
feq[4] = w[4] * rho[i][j] * (1.0 - 3.0 * (v/ck));
feq[5] = w[5] * rho[i][j] * (1.0 + 3.0 * ((u+v)/ck));
feq[6] = w[6] * rho[i][j] * (1.0 + 3.0 * ((-u+v)/ck));
feq[7] = w[7] * rho[i][j] * (1.0 - 3.0 * ((u+v)/ck));
feq[8] = w[8] * rho[i][j] * (1.0 + 3.0 * ((u-v)/ck));
for(k = 0; k < 9; k++)
{
f[k][i][j] = omega * feq[k] + (1. - omega) * f[k][i][j];
}
}
}
//streaming
for(j = Ly-1; j > 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-1; j > 0; j--)
{
for(i = Lx-1; i > 0; i--)
{
f[1][i][j] = f[1][i-1][j];
f[5][i][j] = f[5][i-1][j-1];
}
}
for(j = 0; j < Ly; j++)
{
for(i = Lx-1; 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, Twall is given
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];
}
for(j = 0; j < Ly; j++)
{
f[6][Lx][j] = -f[8][Lx][j];
f[3][Lx][j] = -f[1][Lx][j];
f[7][Lx][j] = -f[5][Lx][j];
f[2][Lx][j] = -f[4][Lx][j];
f[0][Lx][j] = 0.0;
}
for(i = 0; i < Lx; i++)
{
f[8][i][Ly] = -f[6][i][Ly];
f[7][i][Ly] = -f[5][i][Ly];
f[4][i][Ly] = -f[2][i][Ly];
f[1][i][Ly] = -f[3][i][Ly];
f[0][i][Ly] = 0.0;
}
for(i = 0; i < Lx; i++)
{
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];
}
T = 0.0;
for(i = 0; i < Lx; i++)
{
f[2][i][0] = -f[4][i][0];
f[6][i][0] = -f[8][i][0];
f[5][i][0] = -f[7][i][0];
f[1][i][0] = -f[3][i][0];
f[0][i][0] = 0.0;
}
}//end of main loop
for(j = 0; j < Ly; j++)
{
for(i = 0; i < Lx; i++)
{
sum = 0.0;
for(k = 0;k < 9; k++)
{
sum = sum + f[k][i][j];
}
rho[i][j] = sum;
}
}
//cout << rho[0][Ly/2] << endl;
mFile = fopen("d2q9lbmadvdiff.csv","w");
fprintf(mFile,"\t\t----TITLE= D2Q9 RESULTS----\n");
fprintf(mFile," VARIABLES = Temperature distribution\n");
//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;
}
Nov 27, 2014 at 12:23pm UTC
The problem is this:
for (k = 0; k < 9 ; k++) // Note: 9 is out of bounds
check where 9 occurs and use sMod
instead
f[8][i][j] = f[8][i-1][j+1];
8 as an index is out of bounds, change the value of sMod
-> 9