(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)

chux · Accepted Answer · 2021-02-09 16:01:04Z

2

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?

edited Feb 9, 2021 at 16:01

answered Feb 9, 2021 at 14:49

chux

158k17 gold badges160 silver badges311 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

mcgrane5 Over a year ago

hiya. Thank you this got rid of the -1.#IND00 problem. Now the output is a bit bazzar with really long double numbers but its a step in the right direction nonetheless. Thank you very much. I f you have any other suggestions that would be great.

chux Over a year ago

@mcgrane5 Tip: 1) It is more informative to print FP values with "%e" or "%g" than "%f" - especially when debugging. 2) enable all warnings. 3) When allocating use ptr = calloc(n, sizeof *ptr); idiom (Notice no types in that line)

Lutz Lehmann · Accepted Answer · 2021-02-09 15:18:22Z

0

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

answered Feb 9, 2021 at 15:18

Lutz Lehmann

26.1k2 gold badges26 silver badges57 bronze badges

2 Comments

mcgrane5 Over a year ago

Thank you. I rearranged my loops as by your code above but Its still not working as i expect. I have the same problem coded out in python and it works as i want. Ill continue to try debug it. thank you for your contribution though

mcgrane5 Over a year ago

hi there, just to let you know and perhaps anyone else who has a similar problem with a program of the same type i figured out the bug. I replaced the line usol_array[ i ] = u0_array[ i ] with memcpy( usol_array, u0_array, nx * sizeof(double) ) to copy the arrays contents that way and the program runs as expected. Thank you for all your help also

Collectives™ on Stack Overflow

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

2 Answers 2

2 Comments

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

2 Comments

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related