Seismic Model code runs in CodeBlocks but not in Qt.

[Flavio Mesquita]

Hy everyone,
I´m away from the forums for a time, working on a academic research project that I´m willing to share with the community. I´ll describe the idea, I got an old code written in fortran and i´m translating it to C++ in order to make it more resourcefull, like making an interface, being able to use GPU to accelerate the processing, create a plot of the data files generated, etc.
The code gets a underground earth velocity model, then u fire an acoustic pulse, it propagates in this earth model and reflects at the boundaries, so the receivers at the surface get this reflections back and creates a seismogram that later on can be used as a input for another model or simply processed in order to visualize this underground.
So i wrote a code in CodeBlocks, got a very simple 3 layer model, adn propagated one pulse. It works altough the results are still kind of weird, the result images have an inclination of 45 degress, it should be horizontal. Anyway I tried to put it in Qt, but it gives a segmentation fault , but in CodeBlocks it was not giving this error.
"Value conversion issue
implicit conversion changes signedness: int to unsigned int"

Also it gives lots of warnings like:
"-use of old style cast"
What is the new style then?

If anyone is interested the code is shared here:

////Propag-Bin.c
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>



double fon_d2g(double t, double fc)
//*** Seismic source function***
//Source: Ricker(Gaussian second derivative )
{
double amp;
double pi=3.141592653589793238460;
amp = - pi * (pi * fc * t) * (pi * fc * t);
amp = exp(amp);
amp *= 1.0 - 2.0 * pi * (pi * fc * t) * (pi * fc * t);
return (amp);
}

int main()
{
//Variables declaration
int Nx, Nz, Nt, Ns, ixf, izf, iprof, k_snap, k_sismo, dt_sismo, dt_snap, cerjan_n, i, j, k;
float v;
double Dt, h, fcorte, pi, nf, tf, fc, cput, cerjan_fat, *g, **vel, **sismo, **p1, **p2, **p3;
char *modelo, *sismo_arq;
FILE *parametros; //data input file name
//FILE *locais; //file that indicates where to save the files created
FILE *vel_arq; //velocity files name
FILE *varsnap; //snapshots file name
FILE *seismogram; //seismogram file name
clock_t start, end;
//Data input
parametros = fopen ("Parametros.txt", "r");
fscanf (parametros, "%d %d %lf %d %d %d %lf %d", &Nx, &Nz, &h, &ixf, &izf, &iprof, &Dt, &Nt);
//Nx x Nz = grid size
//h = receivers spacing
//ixf x izf = seismic source position
//iprof = seismogram reading depth
//Dt = delta_t
//Nt = time increments
//modelo = velocity file name
//sismo_arq = seismogram file name
pi=3.141592653589793238460; //Pi
fcorte = 40.0; //cutting frequency
cerjan_n = 100;
cerjan_fat = 0.00105; //AB Cerjan
//Number of time steps for**********************:
dt_snap = 400; //saves snapshots every 400 iterations
dt_sismo = 10; //Writes one seismogram sample every 10 iterations
Ns = Nt / dt_sismo; //Number of Seismogram samples
g = (double *) malloc (cerjan_n * sizeof(double));
vel = (double **) malloc (Nz * sizeof(double *));
for (i = 0; i < Nz; i++)
vel[i] = (double *) malloc (Nx * sizeof(double));
sismo = (double **) malloc (Ns * sizeof(double *));
for (i = 0; i < Ns; i++)
sismo[i] = (double *) malloc (Nx * sizeof(double));
p1 = (double **) malloc (Nz * sizeof(double *));
for (i = 0; i < Nz; i++)
p1[i] = (double *) malloc (Nx * sizeof(double));
p2 = (double **) malloc (Nz * sizeof(double *));
for (i = 0; i < Nz; i++)
    p2[i] = (double *) malloc (Nx * sizeof(double));
p3 = (double **) malloc (Nz * sizeof(double *));
for (i = 0; i < Nz; i++)
p3[i] = (double *) malloc (Nx * sizeof(double));
modelo = (char *) malloc (256 * sizeof(char));
sismo_arq = (char *) malloc (256 * sizeof(char));
modelo = "vp_model_400x300.bin";
sismo_arq = "Seismogram.bin";
fclose(parametros);
//p1(Nx,Nz) corresponds to time k-1
//p2(Nx,Nz) corresponds to time k
//p3(Nx,Nz) corresponds to time k+1

//Data Processing

//Processing time
start = clock();
//Source term calculation
nf = 4 * sqrt(pi) / (fcorte * Dt);
tf = 2 * sqrt(pi) / fcorte;
fc = fcorte / (3.0 * sqrt(pi));
//Dumping factors calculation
for (i = 0; i < cerjan_n; i++)
g[i] = exp( - pow(cerjan_fat * (cerjan_n - i), 2) );
printf("Source application length: Nf = %.4f\n", nf);
printf("Seismogram temporal samples number: Ns = %d\n", Ns);
//reading velocity grid
printf("\nSaving velocity model...\n");
vel_arq = fopen(modelo, "rb");
for (j = 0; j < Nx; j++)
for (i = 0; i < Nz; i++)
{
fread (&v, sizeof(float), 1, vel_arq);
vel[i][j] = (double) v;
}
fclose(vel_arq);
//Wavefield initialization (initial condition)
for (i = 0; i < Nz; i++)
for (j = 0; j < Nx; j++)
{
p1[i][j] = 0;
p2[i][j] = 0;
}
k_snap = 0; //Snapshots counter
k_sismo = 0; //Seismogram sample counter
printf("\nStart\n");
for (k = 0; k < Nt; k++)//Start time loop
{
//Printing the time loop on the screen
if (k % 50 == 0)
printf("n = %d\n", k);
//Applying the source
if ((double) k <= nf)
p1[izf][ixf] -= fon_d2g((float) (k - 1) * Dt - tf, fc);
//FDM operator (involving the whole grid, except boundaries)
//Fourth order spatial operator
for (i = 2; i < Nz - 2; i++)
    for (j = 2; j < Nx - 2; j++)
p3[i][j] = pow(vel[i][j] * (Dt/12*h), 2)  * (-(p2[i-2][j] + p2[i][j-2] + p2[i+2][j] + p2[i][j+2]) + 16 * (p2[i-1][j] + p2[i][j-1] + p2[i+1][j] + p2[i][j+1]) - 60 * p2[i][j]) + 2 * p2[i][j] - p1[i][j];
//Second order spatial operator
//p3[i][j] = pow(vel[i][j] * (Dt/h), 2)  * ( p2[i+1][j] - 2*p2[i][j] + p2[i-1][j] + p2[i][j+1] - 2*p2[i][j] + p2[i][j-1]) + 2*p2[i][j] - p1[i][j];

//Superior and inferior
for (j = 2; j < Nx - 2; j++)
{
p3[1][j] = (pow(vel[1][j] * Dt / h, 2)) * (p2[2][j] + p2[1][j+1] + p2[0][j] + p2[1][j-1] - 4 * p2[1][j]) + 2 * p2[1][j] - p1[1][j];
p3[Nz-2][j] = (pow(vel[Nz-2][j] * Dt / h, 2)) * (p2[Nz-1][j] + p2[Nz-2][j+1] + p2[Nz-3][j] + p2[Nz-2][j-1] - 4 * p2[Nz-2][j]) + 2 * p2[Nz-2][j] - p1[Nz-2][j];
}
//Left and Right
for (i = 1; i < Nz - 1; i++)
{
p3[i][1] = (pow(vel[i][1] * Dt / h, 2)) * (p2[i+1][1] + p2[i][2] + p2[i-1][1] + p2[i][0] - 4 * p2[i][1]) + 2 * p2[i][1] - p1[i][1];
p3[i][Nx-2] = (pow(vel[i][Nx-2] * Dt / h, 2)) * (p2[i+1][Nx-2] + p2[i][Nx-1] + p2[i-1][Nx-2] + p2[i][Nx-3] - 4 * p2[i][Nx-2]) + 2 * p2[i][Nx-2] - p1[i][Nx-2];
}
//Applying non-reflexive conditions on the boundaries
//Superior
for (j = 1; j < Nx - 1; j++)
p3[0][j] = 0; //Free surface condition
//Left
for (i = 0; i < Nz - 1; i++)
p3[i][0] = p2[i][0] + (Dt * vel[i][0] / h) * (p2[i][1] - p2[i][0]);
//Right
for (i = 0; i < Nz - 1; i++)
p3[i][Nx-1] = p2[i][Nx-1] + (Dt * vel[i][Nx-1] / h) * (p2[i][Nx-2] - p2[i][Nx-1]);
//Inferior
for (j = 0; j < Nx; j++)
p3[Nz-1][j] = p2[Nz-1][j] + (Dt * vel[Nz-1][j] / h) * (p2[Nz-2][j] - p2[Nz-1][j]);
//Applying dumping zone close to the boundaries
//Left
for (i = 0; i < Nz; i++)
for (j = 0; j < cerjan_n; j++)
{
p3[i][j] *= g[j];
p2[i][j] *= g[j];
}
//Right
for (i = 0; i < Nz; i++)
    for (j = 0; j < cerjan_n; j++)
{
p3[i][Nx - j - 1] *= g[j];
p2[i][Nx - j - 1] *= g[j];
}
//Inferior
for (i = 0; i < cerjan_n; i++)
for (j = 0; j < Nx; j++)
{
p3[Nz - i - 1][j] *= g[i];
p2[Nz - i - 1][j] *= g[i];
}
//Saves a seismogram sample every "dt_sismo" time steps

if (k % dt_sismo == 0)
{
for (j = 0; j < Nx; j++)
sismo[k_sismo][j] = p3[iprof][j];
k_sismo++;
}

//Snapshots print every "dt_snap" time intervals
if (k % dt_snap == 0 && k > 0)
{
char *nome;
nome = (char *) malloc (strlen(sismo_arq) * sizeof(char));
for (i = 0; i < strlen(sismo_arq)-14; i++)
nome[i] = sismo_arq[i];
nome[i] = 'S';
nome[i+1] = 'n';
nome[i+2] = 'a';
nome[i+3] = 'p';
nome[i+4] = '_';
nome[i+5] = '\0';
char snap[2];
itoa(k / dt_snap, snap, 10);
strcat(nome, snap);
strcat(nome, ".bin");
varsnap = fopen(nome, "wb");
printf("Snapshot = %s\n", nome);
for (j = 0; j < Nx; j++)
for (i = 0; i < Nz; i++)
{
v = (float) p3[i][j];
fwrite(&v, sizeof(float), 1, varsnap);
}
fclose(varsnap);
free (nome);
k_snap++;
}

//Wavefield update for the next time loop
for (i = 0; i < Nz; i++)
for (j = 0; j < Nx; j++)
p1[i][j] = p2[i][j];
for (i = 0; i < Nz; i++)
for (j = 0; j < Nx; j++)
p2[i][j] = p3[i][j];
}
//Data output
//Calculate processing time and prints it
end = clock();
cput = ((double) (end - start)) / CLOCKS_PER_SEC;
printf("\nExecution time: %.4f s.", cput);
seismogram = fopen(sismo_arq, "wb");
for(j = 0; j < Nx; j++)
for(i = 0; i < Ns; i++)
{
v = (float) sismo[i][j];
fwrite (&v, sizeof(float), 1, seismogram);
}
fclose(seismogram);
free(g);
free(vel);
free(sismo);
free(p1);
free(p2);
free(p3);
free(modelo);
free(sismo_arq);
return 0;
}

The input parameters are as follows:
400
300
3.5
150
6
6
0.00015
4000

they can be changed according to the model grid size used, receivers spacing, position of the source, depth of the receiver, dt, number of iterations.

[Flavio Mesquita]

I also put a bin viewer program for windows plus the simple model I´m using to test the code, plus the input parameters file on a github repository.
(https://github.com/demeflac/direct-wave-propagation)

[SGaist]

Hi,

int value = (int)someVariable; <= this is a C-style cast.
In C++ you should rather do something like int value = static_cast<int>(someVariable);

Where exactly do you get your crash ?

[Flavio Mesquita]

@SGaist Iḿ getting the crash when allocating memory for the velocity vector:
vel = (double **) malloc (Nz * sizeof(double *));

Also, do I literally have to write 'static_cast' in front of type?

[SGaist]

Are you implementing technique number 4 described here ?