[Vishish]

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;
}

[fresty]

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...

[Vishish]

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];

[fresty]

Originally Posted by Vishish

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]);

}

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?

[Vishish]

Originally Posted by fresty

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.