Loop returns -1073741819 (0xC0000005) in C

Question

I have to solve a Finite Difference problem (2d Heat Transfer, non-transient) in C, however when I try to solve this for a (nx x m) arrays, for nx >70, the program crashes and returns -1073741819 (0xC0000005). What should I do? By the end, the script should calculate temperature fields in a composite medium. This code tries to put on evidence the temperature at surface T1[m-1][:] (its last element).

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



int main(){

    int i, j,k;
    int nx = 200;   // número de pontos em x
    int n = 50;
    int m=n;      // número de pontos em b1 e b2


    srand((unsigned)time(NULL));

    // Parâmetros de simulação

    float b1 = 0.01;         // altura do meio 1, m
    float b2 = 0.01;         // altura do meio 2, m
    float a  = 0.04;         // largura das placas, m
    float q  = -100000.0;      // fonte de calor na face superior [W/m2]
    float k1e = 54.0;          // condutividade térmica do material 1 [W/mC]
    float k2 = 54.0;           // condutividade térmica do material 1 [W/mC]
    float Tb = 0;            //temperatura de referência .C
    float hmax = 1000;

    float x[nx], y1[m], y2[n];
    float dx = a/(nx-1);
    float dy = b1/(n-1);
    y1[0] = b1; x[0] = 0.0; y2[0] = 0.0;

    for (i=1; i<nx; i++){
        x[i] = x[i-1]+dx;
    }
    for (j=1; j<m; j++){
        y1[j]= y1[j-1] + dy;
        y2[j]= y2[j-1] + dy;
    }

    float T1[n][nx], T2[m][nx];
    float T1old[n][nx], T2old[m][nx];

    //===============  Cálculo de hc =================
    // Caso 1
    //char case = 'CP1h1 '

    float hc[nx];

    for (i=0; i<nx;i++){
        if (x[i]<a/4)
            hc[i] = hmax;
        else if (x[i] > 3*a/4)
            hc[i] = hmax;
        else
            hc[i] = 0.;

    }
//==============================================


    int nexp = 1;
    float diff = 0.0, maxdiff = 0.0;
    float err1 = 10.0;
    float err2 = 10.0;
    float tol = pow(10.0,-12);
    int contador = 0;



    float k1=k1e;
    for (int k=0; k<nexp; k++){
    //=======================================================================================
        // Valor sintético de k1
        float sigmak1 = 0.00;
        float sigma = sigmak1*k1e;
        float u = (float)rand()/ RAND_MAX;
        float v = (float)rand()/ RAND_MAX;
        float eps = (2*M_PI*v)*sqrt(-2*log(u));
        float k1 = k1e + eps*sigma;
        //======================================================================================


        float k1dy=k1/dy;
        float k1k2=k1/k2;
        float twoqdyk1=2*q*dy/k1;
        for (i=0; i<nx; i++){
            for (j=0; j<m; j++){
                T1[j][i] = 0.;
                T2[j][i] = 0.;
                T1old[j][i] = 0.;
                T2old[j][i] = 0.;
            }
        }

        diff = abs(T1old[0][0] - T1[0][0]);
        maxdiff = diff;
        contador = 0;

        while (contador<14000){
            contador = contador +1;
            for (i=0; i<nx; i++){
                for (j=0; j<nx;j++){
                    diff = abs(T1old[j][i] - T1[j][i]);
                    T1old[j][i] = T1[j][i];
                    T2old[j][i] = T2[j][i];
                    printf("%f \n",diff);
                    if (diff>=maxdiff){
                        maxdiff = diff;
                    }
                }
            }
            err1 = maxdiff;
            //printf("%i \n", contador);

            for (i=0; i<nx;i++){
                for(j=0; j<m;j++){
                    if (i==0 && j!=0 && j!=m-1){ //condição de contorno direita -> exclui j=0 e j=n-1
                        T1[j][i] = 0.25*(T1old[j-i][i] + T1old[j+i][i] + 2*T1old[i+1][j]);
                        T2[j][i] = 0.25*(T2old[j-i][i] + T2old[j+i][i] + 2*T2old[i+1][j]);
                    }
                    else if (i==nx-1 & j!=0 && j!=m-1){ //condição de contorno direita -> exclui j=0 e j=m-1
                        T1[j][i] = 0.25*(T1old[j-i][i] + T1old[j+i][i] + 2*T1old[i-1][j]);
                        T2[j][i] = 0.25*(T2old[j-i][i] + T2old[j+i][i] + 2*T2old[i-1][j]);
                    }
                    else if(j==0){ //condição de contorno inferior -> inclui i=0 e i=nx-1
                        T1[j][i] = (1/(k1dy + hc[i]))*(hc[i]*T2old[n-1][i] + (k1dy)*(T1old[j+1][i]));
                        T2[j][i] = Tb;

                    }
                    else if(j==m-1){//condição de contorno superior -> inclui i=0 e i=nx-1
                        T1[j][i] = 0.25*(2*T1old[j-1][i] -(2*twoqdyk1)+ T1old[j][i-1] + T1old[j][i-1]);
                        T2[j][i] = (-k1k2)*(T1old[0][i]-T1[1][i]) + T2old[j-1][i];

                    }
                    else{ //pontos internos
                        T1[j][i] = 0.25*(T1old[j-1][i]+ T1old[j+1][i]+ T1old[j][i-1] + T1old[j][i-1]);
                        T2[j][i] = 0.25*(T2old[j-1][i]+ T2old[j+1][i]+ T2old[j][i-1] + T2old[j][i-1]);
                    }

                }

            }
        }

    }




    for (i=1; i<nx;i++){
        printf("%f \n     ", T1[m-1][i]);
    }
    return 0;
}

Answer 1

Look closely at this piece of code:

  for (i = 0; i < nx; i++) {
    for (j = 0; j < nx; j++) {
      diff = abs(T1old[j][i] - T1[j][i]);

The inner loop varies j from 0 to nx but it should be n.

  for (i = 0; i < nx; i++) {
    for (j = 0; j < n; j++) {   // << use n instead of nx
      diff = abs(T1old[j][i] - T1[j][i]);

Definition of T1:

float T1[n][nx]

Sou what you have here is accessing an array with an out of range index which result in undefined behaviour.

Answer 2

At least these issues

Not all warnings enabled

Save time, enable them all.

int absolute

abs(T1old[0][0] - T1[0][0]); determines the int absolute value. This leads to incorrect functionality and undefined behavior with values much outside the int range.

Use a floating point function: fabsf(T1old[0][0] - T1[0][0]);

log(0.0} possible which is -INF

Below may generate u == 0.0

float u = (float)rand()/ RAND_MAX;
float eps = (2*M_PI*v)*sqrt(-2*log(u));  // Bad

Perhaps ?

float u = (rand() + 0.5f)/ (RAND_MAX + 1u);

Why float?

Code mixes double function calls and constants with float variables. Recommend to just use double throughout.

Debug tips

Comment out srand((unsigned)time(NULL)); until code fully runs at least once correctly.

Print floating point values with %g or %e rather than %f - at least during debug. It is more informative and less noisy.