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

Have you tried running your code line by line in a debugger while monitoring the values of all variables, in order to determine at which point your program stops behaving as intended? If you did not try this, then you may want to read this: [What is a debugger and how can it help me diagnose problems?](https://stackoverflow.com/q/25385173/12149471) You may also want to read this: [How to debug small programs?](https://ericlippert.com/2014/03/05/how-to-debug-small-programs/). — Andreas Wenzel, Jan 17 '22 at 13:20
One of your array indexes probably runs off the array. I note, for example, that `T1` and `T2` have `n` and `m` as their first dimension, respectively, but further down you use loops running indexes up to to `m` for both of them. — 500 - Internal Server Error, Jan 17 '22 at 13:21
@500-InternalServerError, in this case m=n, I am note sure if this is the cause of the problem. — Cairo Martins, Jan 17 '22 at 14:01
Your code is not fully ISO C compliant, as you are using the macro constant `M_PI`, which seems to be [a GCC extension](https://www.gnu.org/software/libc/manual/html_node/Mathematical-Constants.html). Other compilers may possibly support this extension too. You may want to mention in your question that you are using platform-specific extensions, and mention which compiler you are using. This may make it easier for other people to reproduce your problem. — Andreas Wenzel, Jan 17 '22 at 14:31

chux - Reinstate Monica · Accepted Answer · 2022-01-17T13:55:59.633

1

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.

edited Jan 17 '22 at 13:55

answered Jan 17 '22 at 13:49

chux - Reinstate Monica

143,097
13
135
256

Thanks for the appointments, which were added to the main code. However, for a larger number of "nx", the code still crashes. Is there a better way to declare the arrays? – Cairo Martins Jan 17 '22 at 13:59
@CairoMartins no, there is no better way, but there is a way to not use out of range indexes. See answer below. – Jabberwocky Jan 17 '22 at 14:00
@Jabberwocky Yes there are better ways, yet OP's code is rife with indexing issues and those issues need prompt attention. – chux - Reinstate Monica Jan 17 '22 at 14:02
@CairoMartins the remarks in this answer are all correct, yet the out of range problem is still not resolved. If your codes appears to work now, it might be just anorhte manifestation of undefind behaviour. – Jabberwocky Jan 17 '22 at 14:05

score 0 · Answer 2 · answered Jan 17 '22 at 13:59

0

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.

answered Jan 17 '22 at 13:59

Jabberwocky

48,281
17
65
115

Thanks for the notes, it actually was a typo when writing my question. I have made the changes, but my code still won't work for nx>70. – Cairo Martins Jan 17 '22 at 14:36
Then it's probably time to start learning how to use your debugger. I've used my debugger to spot this problem and it took me 5 minutes. – Jabberwocky Jan 17 '22 at 14:55

Loop returns -1073741819 (0xC0000005) in C

2 Answers2