//Dhruv Balwada(MATS 379)
//Program for solving Shallow water Equations
//Using the Method of Characterstics
//For constant ocean depth(still working on variable depth)
#include
#include
int main()
{
float dt, dx;
float g=9.81;
float w,a,m,r,l;
float u[1000][2],c[1000][2];
float ita, n[1000], h[1000];
int i,j,k;
FILE *fp;
fp=fopen("wave_eqn_data.dat","w");
// Characterstics
dt=1;
dx=100;
a=0.5;
w=4.398e-3;
m=0;
l=6.28/100000;
//Check CFL condtion
r=dt/dx;
if(r>0.5)
exit(0);
//Initialisation of the array
for(i=0; i<2; i++)
{
for(j=0; j<1000; j++)
{
u[j][i]=0;
c[j][i]=0;
}
}
//The values of h
for(i=0; i<1000; i++)
{
h[i]=m*(1000-i)*dx+200;
}
//The values of c
for(i=0; i<1000; i++)
{
c[i][0]=sqrt(g*h[i]);
}
// Big Time Loop
for(k=0; k<4000; k++)
{
//Boundary Condition Left
ita=a*sin(w*dt*k);
u[0][1]=ita*sqrt(g/h[0]);
c[0][1]=c[0][0]+0.5*(u[0][1]-u[0][0])-m*dt/2+(u[0][0]-c[0][0])*(dt/dx)*(0.5*(u[1][0]-u[0][0])-(c[1][0]-c[0][0]));
//Boundary Condition Right
u[999][1]=u[999][0]-(dt/(2*dx))*(3*u[999][0]-m*dt*k)*(u[999][0]-u[998][0])+m*dt;
c[999][1]=c[999][0]-(dt/dx)*(3*c[999][0]+m*k*dt)*(c[999][0]-c[998][0]);
//Chunk Calculation
for(i=1; i<999; i++)
{
u[i][1]=u[i][0]+(dt/dx)*((u[i][0]+c[i][0])*(0.5*u[i-1][0]-0.5*u[i][0]+c[i-1][0]-c[i][0])-(u[i][0]-c[i][0])*(0.5*u[i+1][0]-0.5*u[i][0]-c[i+1][0]+c[i][0])+m*dx);
c[i][1]=c[i][0]+(dt/(2*dx))*((u[i][0]+c[i][0])*(0.5*u[i-1][0]-0.5*u[i][0]+c[i-1][0]-c[i][0])+(u[i][0]-c[i][0])*(0.5*u[i+1][0]-0.5*u[i][0]-c[i+1][0]+c[i][0]));
}
//Find ita at all the points on the grid
for(j=0; j<1000; j++)
{
n[j]=u[j][1];
}
// Write to data file
if(k%100==0)
{
for(i=0; i<1000; i=i+20)
{
fprintf(fp,"%f",n[i]);
fprintf(fp,"\n");
}
fprintf(fp,"\n");
}
//Moving the Data in time
for(j=0; j<1000; j++)
{
u[j][0]=u[j][1];
c[j][0]=c[j][1];
}
}
fclose(fp);
}
