Nozzle 1D Euler Code using C
 

Can anyone check this code, it's not converging.
Beyond 2000 iterations, it's 'nan'.

Thanks

#include "math.h"
#include "stdio.h"
#include "iostream"
using namespace std;


#include "math.h"
#include "stdio.h"
#include "iostream"
using namespace std;

double S[1001], U1[1001], U2[1001], U3[1001], u[1001], T[1001], gama=1.4, a[1001], M[1001], B[1001], fp_1[1001], fp_2[1001], fp_3[1001],
fn_1[1001], fn_2[1001], fn_3[1001], delta_t, CFL = 0.1, delta_x = 0.1, U01[1001], U02[1001], U03[1001], delt[1001], minimum, resc[1001],
maxc, res_limit = 1e-010;

int i, R=287, n, nmax=1000;

void initialization()
{
for (i=1; i<11; i++)
{
S[i] = 0.130014 + 0.032271*tanh(0.8*i*0.1 - 4);
U1[i] = 1.219;
U2[i] = 428.826;
U3[i] = 195065;
u[i] = U2[i]/U1[i];
T[i] = (U3[i]/U1[i])*((gama-1)/R);
a[i] = sqrt(gama*R*T[i]);
M[i] = u[i]/a[i];
B[i] = U1[i]*R*T[i];

fp_1[i] = 0.25*U1[i]*a[i]*pow(M[i]+1 , 2);
fp_2[i] = 0.25*U1[i]*a[i]*pow(M[i]+1, 2) * ((gama-1)*u[i] + 2*a[i])/gama;
fp_3[i] = 0.25*U1[i]*a[i]*pow(M[i]+1, 2) * pow(((gama-1)*u[i] + 2*a[i]),2) / (2*(1.4*1.4 - 1));
fn_1[i] = -0.25*U1[i]*a[i]*pow(M[i]-1 , 2);
fn_2[i] = -0.25*U1[i]*a[i]*pow(M[i]-1 , 2) * ((gama-1)*u[i] - 2*a[i])/gama;
fn_3[i] = -0.25*U1[i]*a[i]*pow(M[i]-1, 2) * pow(((gama-1)*u[i] - 2*a[i]),2) / (2*(1.4*1.4 - 1));

delta_t = delta_x * CFL / (u[1] + a[1]);

}
}

void update()
{
for (i=2; i<11; i++)
{
U01[i] = U1[i];
U02[i] = U2[i];
U03[i] = U3[i];
U1[i] = U01[i] + (delta_t/delta_x * S[i]) * (S[i]*fp_1[i] - S[i-1]*fp_1[i-1] + S[i+1]*fn_1[i+1] - S[i]*fn_1[i]);
U2[i] = U02[i] + (delta_t/delta_x * S[i]) * (S[i]*fp_2[i] - S[i-1]*fp_2[i-1] + S[i+1]*fn_2[i+1] - S[i]*fn_2[i]+ B[i+1]*S[i+1] - B[i]*S[i]);
U3[i] = U03[i] + (delta_t/delta_x * S[i]) * (S[i]*fp_3[i] - S[i-1]*fp_3[i-1] + S[i+1]*fn_3[i+1] - S[i]*fn_3[i]);

u[i] = U2[i]/U1[i];
T[i] = (U3[i]/U1[i])*((gama-1)/R);
a[i] = sqrt(gama*R*T[i]);
M[i] = u[i]/a[i];
B[i] = U1[i]*R*T[i];

fp_1[i] = 0.25*U1[i]*a[i]*pow(M[i]+1 , 2);
fp_2[i] = 0.25*U1[i]*a[i]*pow(M[i]+1, 2) * ((gama-1)*u[i] + 2*a[i])/gama;
fp_3[i] = 0.25*U1[i]*a[i]*pow(M[i]+1, 2) * pow(((gama-1)*u[i] + 2*a[i]),2) / (2*(1.4*1.4 - 1));
fn_1[i] = -0.25*U1[i]*a[i]*pow(M[i]-1 , 2);
fn_2[i] = -0.25*U1[i]*a[i]*pow(M[i]-1 , 2) * ((gama-1)*u[i] - 2*a[i])/gama;
fn_3[i] = -0.25*U1[i]*a[i]*pow(M[i]-1, 2) * pow(((gama-1)*u[i] - 2*a[i]),2) / (2*(1.4*1.4 - 1));

resc[i] = U1[i] - U01[i];
if (resc[i]<0)
resc[i] = -resc[i];
if(i==2)
maxc=resc[i];
else
{
if (maxc<resc[i])
maxc=resc[i];
}
}
}

void time_step()
{
for (int i=2; i<11; i++)
{
delt[i] = delta_x/fabs((u[i])+a[i]);
if(i==2)
minimum = delt[i];
else
{
if(delt[i]<minimum)
minimum = delt[i];
}
}
delta_t = CFL*minimum;
}


int main()
{
initialization();
for (n=1; n<nmax; n++)
{
update();
time_step();

/*if (maxc >= res_limit )
{
n==n;
}
else
{
nmax=n;
}*/

}

for (i=2; i<10; i++)
{
//cout<<resc[i]<<" "<<U1[i]<<" "<<U01[i]<<endl;
cout<<i<<" "<<n<<" "<<M[i]<<" "<<T[i]<<" "<<maxc<<endl;
}


return 0;
}

Just skimming through it; can see that the index [i+1] in update() function at i = 10, would simply evaluate to zero... that seems problematic.. seems like the differencing scheme needs to be corrected given how you're coding the problem...

Hi,
I've added the following terms in update, but still it gives me the same problem.

fn_1[11]=fn_1[10];
fn_2[11]=fn_2[10];
fn_3[11]=fn_3[10];

Not sure if that would completely solve the problem. Have you had a look at these lines? What are S[i+1] and B[i+1] at i = 10?

I added this term to update()

fn_1[11]=fn_1[10];
fn_2[11]=fn_2[10];
fn_3[11]=fn_3[10];
B[11] = B[10];

S[i+1] would compute from initialization(), as it would be constant for all the points.