(CFD) Problem with C code for 1D linear convection equation

[mcgrane5]

TL;DR Summary

Hi im writing a c program on linear convection in 1D and im having problems with my final for loop. It returns bizarre results and i dont know why

I have the same problem coded in python and it works perfectly. Could anyone offer a fix to my code. Its only the last for loop where im having problems. I have attached my C code below as well as the working python script to help anyone debug

This is my c program below. Again everything works as expected up until the highlighted lines or the second for loop. When i run my program to
output the arrays contents i get bizarre numbers that are not correct.

# include <math.h>
# include <stdlib.h>
# include <stdio.h>
# include <time.h>int main ( )
{
    //eqyation to be solved 1D linear advection -- du/dt + c(du/dx) = 0
    //initial conditions u(x,0) = u0(x)

    //after discretisation using forward difference time and backward differemce space
    //update equation becomes u[i] = u0[i] - c * dt/dx * (u0[i] - u[i - 1]);
    FILE *fpointer = NULL;
    fpointer = fopen("bbb.txt", "w");

    int nx = 41; //num of grid points

    double dx = 2.0 / (nx - 1.0); //magnitude of the spacing between grid points

    int nt = 25;//nt is the number of timesteps

    double dt = 0.25; //the amount of time each timestep covers

    int c = 1;  //assume wavespeed
    //double u0_array;

    //set up our initial conditions. The initial velocity u_0
    //is 2 across the interval of 0.5 <x < 1 and u_0 = 1 everywhere else.
    //we will define an array of ones
    double* u0_array = (double*)calloc(nx, sizeof(double));

    for (int i = 0; i < nx-1; i++)
    {
        if (i * dx >= 0.5 && i * dx <= 1)
        {
            u0_array[i] = 2;
        }else
        {
             u0_array[i] = 1;
        }

    }

   //test code to make sure array is populated correctly
    for (int i = 0; i < nx; i++)
    {
        //u0_array[i] = 2;
        //printf("%f, ", u0_array[i]);

    }

    //make a temporary array that allows us to store the solution
    double* usol_array = (double*)calloc(nx, sizeof(double));

    //apply numerical scheme forward difference in
    //time an backward difference in space
    for (int i = 0; i < nt; i++)
    {
        usol_array[i] = u0_array[i];

        for (int j = 0; j < nx-1; j++)
        {
            //usol_array[i] = u0_array[i];
            u0_array[j] = usol_array[j] - c * dt/dx * (usol_array[j] - usol_array[j - 1]);
            printf("%f\n", u0_array[i]);
        }
        //first loop iterates through each time step
    }
        

    return EXIT_SUCCESS;
}

BELOW IS THE WORKING PYTHON SCRIPT TO HELP ANYONE DEBUG MY C CODE

import numpy                     
import matplotlib.pyplot as plt     
import time, sys   

nx = 41  # try changing this number from 41 to 81 and Run All ... what happens?
dx = 2 / (nx-1)
nt = 25    #nt is the number of timesteps we want to calculate
dt = .025  #dt is the amount of time each timestep covers (delta t)
c = 1

u = numpy.ones(nx)      #numpy function ones()
u[int(.5 / dx):int(1 / dx + 1)] = 2  #setting u = 2 between 0.5 and 1 as per our I.C.s
print(u)

un = numpy.ones(nx) #initialize a temporary array

for n in range(nt):  #loop for values of n from 0 to nt, so it will run nt times
    un = u.copy()
    print(un)##copy the existing values of u into un
    for i in range(1, nx): ## you can try commenting this line and...
    #for i in range(nx): ## ... uncommenting this line and see what happens!
        u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])
        #print(u[i])

plt.plot(numpy.linspace(0, 2, nx), u);
plt.show()

 

[bigfooted]

In your c file, manually compute what the value is that you get for u0_array when j=0 and you should see the problem.