(Computational fluid dynamics) Problem with arrays. Why is my C program outputting -1.#IND00

Question

I have a problem with a program i'm writing in C that solves the 1D linear convection equation. Basically i have initialized an two arrays. The first array (u0_array) is an array of ones with the arrays elements set equal to two over the interval or 0.5 < x < 1. The second array (usol_array) serves as a temporary array that the result or solution will be stored in.

The problem i am running into is the nested for loop at the end of my code. This loop applies the update equation required to calculate the next point, to each element in the array. When i run my script and try to print the results the output i get is just -1.IND00 for each iteration of the loop. (I am following Matlab style pseudo code which i also have attached below) I'm very new with C so my inexperience shows. I don't know why this is happening. If anyone could suggest a possible fix to this i would be very grateful. I have attached my code so far below with a few comments so you can follow my thought process. I have also attached the Matlab style pseudo code i'm followed.

# include <math.h>
# include <stdlib.h>
# include <stdio.h>
# include <time.h>


int main ( )
{
    //eqyation to be solved 1D linear convection -- 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]);


    int nx = 41; //num of grid points

    int dx = 2 / (nx - 1); //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


    //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; i++)
    {
        u0_array[i] = 1;
    }

    // set u = 2 between 0.5 and 1 as per initial conditions
    //note 0.5/dx = 10, 1/dx+1 = 21
    for (int i = 10; i < 21; 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++)
    {
        //first loop iterates through each time step
        usol_array[i] = u0_array[i];
        //printf("%f", usol_array[i]);
        

        //MY CODE WORKS FINE AS I WANT UP TO THIS LOOP
        //second array iterates through each grid point for that time step and applies
        //the update equation
        for (int j = 1; j < nx - 1; j++)
        {
            u0_array[j] = usol_array[j] - c * dt/dx * (usol_array[j] - usol_array[j - 1]);
            printf("%f, ", u0_array[j]);
        }
    }

    return EXIT_SUCCESS;
}

For reference also the pseudo code i am following is attached below 1D linear convection pseudo Code (Matlab Style)

Answer 1

Instead of integer division, use FP math.
This avoid the later division by 0 and -1.#IND00

// int dx = 2 / (nx - 1);   quotient is 0.
double dx = 2.0 / (nx - 1);

---

OP's code does not match comment

// u[i] = u0[i] - c * dt/dx * (u0[i] - u[i - 1]
u0_array[j] = usol_array[j] - c * dt/dx * (usol_array[j] - usol_array[j - 1]);

Easier to see if the redundant _array removed.

//                                     v correct?
// u[i] = u0[i] - c * dt/dx * (u0[i] - u[i - 1]
u0[j] = usol[j] - c * dt/dx * (usol[j] - usol[j - 1]);
//                                       ^^^^ Correct?

Answer 2

What you probably want is in matlab terms

for i = 1:nt
    usol = u0
    u0(2:nx) = usol(2:nx) - c*dt/dx*(usol(2:nx)-usol(1:nx-1))
end%for

This means that you have 2 inner loops over the space dimension, one for each of the two vector operations. These two separate loops have to be explicit in the C code

    //apply numerical scheme forward difference in
    //time an backward difference in space
    for (int i = 0; i < nt; i++) {
        //first loop iterates through each time step
        for (int j = 0; j < nx; j++) {
            usol_array[j] = u0_array[j];
            //printf("%f", usol_array[i]);
        } 
        

        //second loop iterates through each grid point for that time step 
        //and applies the update equation
        for (int j = 1; j < nx; j++)
        {
            u0_array[j] = usol_array[j] - c * dt/dx * (usol_array[j] - usol_array[j - 1]);
            printf("%f, ", u0_array[j]);
        }
    }