I have a program for the propagation of seismic waves, when I compile it is fine, but when I run it I get this error:
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0 0x7F65346DEE08
#1 0x7F65346DDF90
#2 0x7F65340264AF
#3 0x4032F6 in MAIN__ at ondas.f90:?
The problem happens when I add more variables to the program:
program vix
implicit none
integer::i,j,k,m,n,aux,nt,time,mpmlx,npmlz,ntz,mtx
double precision::dx,dt,dz,pi,alfa,t0,f0,xi,xf,ti,tf,zi,zf,LENGTHX,LENGTHZ,VPMAX,CON_ESTA
double precision,allocatable,dimension(:)::f,t
double precision,allocatable,dimension(:,:)::p,mu,lambda,VS,VP,S_Vx,S_Vz
double precision,allocatable,dimension(:,:,:)::v_x,v_z,o_xx,o_xz,o_zz
character(20) tt, title
open(unit=1,file='viex.dat',status='replace',action='write')
time=1
!Dimensiones Espaciales
LENGTHX=2000.d0
LENGTHZ=2000.d0
dx=10.d0
dz=10.d0
m=int(LENGTHX/dx)
n=int(LENGTHX/dz)
!Dimensiones Temporales
ti=0.d0
tf=2.0d0
dt=0.003d0
nt=int(tf/dt)+1
!Dimensiones de la frontera absorbente
mpmlx=20
npmlz=20
mtx=m+mpmlx*2
ntz=n+npmlz
!----------------Velocities Model-------------------!
allocate(VP(n,m),VS(n,m),p(n,m))
call lectura('VelocitiesP.dat',VP,n, m)
call lectura('VelocitiesS.dat',VS,n, m)
call lectura('Densities.dat',p,n,m)
!---------------Constantes elasticas-----------------!
allocate(mu(n,m),lambda(n,m))
mu=p*VS**2
lambda=p*(VP**2)-2.0d0*mu
allocate (t(nt),f(nt),o_xx(n,m,2),o_zz(n,m,2),o_xz(n,m,2))
allocate(v_x(n,m,2),v_z(n,m,2))
v_x=0.0
v_z=0.0
o_xx=0.0
o_zz=0.0
o_xz=0.0
!---------------Condicion de estabilidad-------------!
VPMAX=maxval(VP)
CON_ESTA=VPMAX*dt*sqrt((1.d0/dx**2)+(1.d0/dz**2))
if (CON_ESTA.GE.1) then
print*, 'Error en codicion de estabilidad'
print*, 'La condicion de estabilidad debe ser menor a 1'
print*, CON_ESTA
else !Realiza todo el proceso de propagación
!!!allocate ,(S_Vx(m,nt),S_Vz(m,nt))
S_Vx=0.d0
S_Vz=0.d0
!-----------------Unidades para sismogramas----------------!
open(unit=98,file='sismograma_x.dat',status='replace')
open(unit=99,file='sismograma_z.dat',status='replace')
!------------------Fuente Gaussiana----------------!
pi=4.0*atan(1.0) !pi
t0=0.25 !periodo dominante
f0=1.d0/t0 !frecuancia dominante
alfa=2.d0*((pi*f0)**2) !alfa
do k=1,nt !inicio del ciclo de tiempo
print*,'muestra de tiempo',k
t(k)=ti+(k-1)*dt
f(k)=-2.0*alfa*(t(k)-t0)*exp(-alfa*(t(k)-t0)**2) !fuente inicial
o_xx(100,100,1)=o_xx(100,100,1)+f(k)
o_zz(100,100,1)=o_zz(100,100,1)+f(k)
!------------------VZ------------!!-----------------VX-----------!
do j=1,m
do i=1,n
v_z(i,j,2)=v_z(i,j,1)+(1.0/p(i,j))*(dt/dz*(o_zz(i+1,j,1)-o_zz(i,j,1))+dt/dx*(o_xz(i,j+1,1)-o_xz(i,j,1)))
end do
end do
do j=2,m
do i=1,n
v_x(i,j,2)=v_x(i,j,1)+(1.0/p(i,j))*(dt/dx*(o_xx(i,j,1)-o_xx(i,j-1,1))+dt/dz*(o_xz(i,j,1)-o_xz(i-1,j,1)))
end do
end do
!----------------OXX------OZZ--------OXZ---------------------!
do j=1,m
do i=2,n
o_xx(i,j,2)=o_xx(i,j,1)+(lambda(i,j)+2.0d0*mu(i,j))*(dt/dx)*(v_x(i,j+1,2)-v_x(i,j,2))&
+lambda(i,j)*(dt/dz)*(v_z(i,j,2)-v_z(i-1,j,2))
o_zz(i,j,2)=o_zz(i,j,1)+(lambda(i,j)+2.0d0*mu(i,j))*(dt/dx)*(v_z(i,j,2)-v_z(i-1,j,2))&
+lambda(i,j)*(dt/dz)*(v_x(i,j+1,2)-v_x(i,j,2))
enddo
enddo
do i=1,m
do j=1,n
o_xz(i,j,2)=o_xz(i,j,1)+mu(i,j)*((dt/dz)*(v_x(i+1,j,2)-v_x(i,j,2))+(dt/dx)*(v_z(i,j,2)&
-v_z(i,j-1,2)))
enddo
enddo
write(1,*),t(k),f(k),o_xx(5,5,1),o_zz(5,5,1),v_x(i,1,1)
!-----------Sismogramas-------------!
write(98,*)(v_x(i,2,2),i=1,m)
write(99,*)(v_z(i,2,2),i=1,m)
!------------Snapshots--------------!
if (k==time) then
write(tt,'(i4)')time
title='t'//trim(adjustl(tt))//'.txt'
open(unit=10,file=title,status='replace')
do i=1,n
write(10,*)(o_xx(i,j,2),j=1,m)
end do
time=time+40
close(10)
end if
v_x(:,:,1)=v_x(:,:,2)
v_z(:,:,1)=v_z(:,:,2)
o_xx(:,:,1)=o_xx(:,:,2)
o_zz(:,:,1)=o_zz(:,:,2)
o_xz(:,:,1)=o_xz(:,:,2)
end do !final del ciclo de tiempo
end if !final de la estabilidad
call system ("gnuplot Script_vix.gp")
close(1)
end program
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!.........SUBRUTINA PARA LEER..........................
subroutine lectura(Archivo,Arreglo,z,x)
implicit none
integer::i,j
integer, intent (in)::z,x
double precision, dimension (1:z,1:x), intent(out)::Arreglo
character (len=*), intent (in)::Archivo
open (unit=100,file=Archivo,status='old')
do i=1, z
read(100,*)(Arreglo(i,j),j=1,x)
end do
close(100)
end subroutine lectura
program vix
implicit none
integer::i,j,k,m,n,aux,nt,time,mpmlx,npmlz,ntz,mtx
double precision::dx,dt,dz,pi,alfa,t0,f0,xi,xf,ti,tf,zi,zf,LENGTHX,LENGTHZ,VPMAX,CON_ESTA
double precision,allocatable,dimension(:)::f,t
double precision,allocatable,dimension(:,:)::p,mu,lambda,VS,VP,S_Vx,S_Vz
double precision,allocatable,dimension(:,:,:)::v_x,v_z,o_xx,o_xz,o_zz
character(20) tt, title
open(unit=1,file='viex.dat',status='replace',action='write')
time=1
!Dimensiones Espaciales
LENGTHX=2000.d0
LENGTHZ=2000.d0
dx=10.d0
dz=10.d0
m=int(LENGTHX/dx)
n=int(LENGTHX/dz)
!Dimensiones Temporales
ti=0.d0
tf=2.0d0
dt=0.003d0
nt=int(tf/dt)+1
!Dimensiones de la frontera absorbente
mpmlx=20
npmlz=20
mtx=m+mpmlx*2
ntz=n+npmlz
!----------------Velocities Model-------------------!
allocate(VP(n,m),VS(n,m),p(n,m))
call lectura('VelocitiesP.dat',VP,n, m)
call lectura('VelocitiesS.dat',VS,n, m)
call lectura('Densities.dat',p,n,m)
!---------------Constantes elasticas-----------------!
allocate(mu(n,m),lambda(n,m))
mu=p*VS**2
lambda=p*(VP**2)-2.0d0*mu
allocate (t(nt),f(nt),o_xx(n,m,2),o_zz(n,m,2),o_xz(n,m,2))
allocate(v_x(n,m,2),v_z(n,m,2))
v_x=0.0
v_z=0.0
o_xx=0.0
o_zz=0.0
o_xz=0.0
!---------------Condicion de estabilidad-------------!
VPMAX=maxval(VP)
CON_ESTA=VPMAX*dt*sqrt((1.d0/dx**2)+(1.d0/dz**2))
if (CON_ESTA.GE.1) then
print*, 'Error en codicion de estabilidad'
print*, 'La condicion de estabilidad debe ser menor a 1'
print*, CON_ESTA
else !Realiza todo el proceso de propagación
!!!allocate ,(S_Vx(m,nt),S_Vz(m,nt))
S_Vx=0.d0
S_Vz=0.d0
!-----------------Unidades para sismogramas----------------!
open(unit=98,file='sismograma_x.dat',status='replace')
open(unit=99,file='sismograma_z.dat',status='replace')
!------------------Fuente Gaussiana----------------!
pi=4.0*atan(1.0) !pi
t0=0.25 !periodo dominante
f0=1.d0/t0 !frecuancia dominante
alfa=2.d0*((pi*f0)**2) !alfa
do k=1,nt !inicio del ciclo de tiempo
print*,'muestra de tiempo',k
t(k)=ti+(k-1)*dt
f(k)=-2.0*alfa*(t(k)-t0)*exp(-alfa*(t(k)-t0)**2) !fuente inicial
o_xx(100,100,1)=o_xx(100,100,1)+f(k)
o_zz(100,100,1)=o_zz(100,100,1)+f(k)
!------------------VZ------------!!-----------------VX-----------!
do j=1,m
do i=1,n
v_z(i,j,2)=v_z(i,j,1)+(1.0/p(i,j))*(dt/dz*(o_zz(i+1,j,1)-o_zz(i,j,1))+dt/dx*(o_xz(i,j+1,1)-o_xz(i,j,1)))
end do
end do
do j=2,m
do i=1,n
v_x(i,j,2)=v_x(i,j,1)+(1.0/p(i,j))*(dt/dx*(o_xx(i,j,1)-o_xx(i,j-1,1))+dt/dz*(o_xz(i,j,1)-o_xz(i-1,j,1)))
end do
end do
!----------------OXX------OZZ--------OXZ---------------------!
do j=1,m
do i=2,n
o_xx(i,j,2)=o_xx(i,j,1)+(lambda(i,j)+2.0d0*mu(i,j))*(dt/dx)*(v_x(i,j+1,2)-v_x(i,j,2))&
+lambda(i,j)*(dt/dz)*(v_z(i,j,2)-v_z(i-1,j,2))
o_zz(i,j,2)=o_zz(i,j,1)+(lambda(i,j)+2.0d0*mu(i,j))*(dt/dx)*(v_z(i,j,2)-v_z(i-1,j,2))&
+lambda(i,j)*(dt/dz)*(v_x(i,j+1,2)-v_x(i,j,2))
enddo
enddo
do i=1,m
do j=1,n
o_xz(i,j,2)=o_xz(i,j,1)+mu(i,j)*((dt/dz)*(v_x(i+1,j,2)-v_x(i,j,2))+(dt/dx)*(v_z(i,j,2)&
-v_z(i,j-1,2)))
enddo
enddo
write(1,*),t(k),f(k),o_xx(5,5,1),o_zz(5,5,1),v_x(i,1,1)
!-----------Sismogramas-------------!
write(98,*)(v_x(i,2,2),i=1,m)
write(99,*)(v_z(i,2,2),i=1,m)
!------------Snapshots--------------!
if (k==time) then
write(tt,'(i4)')time
title='t'//trim(adjustl(tt))//'.txt'
open(unit=10,file=title,status='replace')
do i=1,n
write(10,*)(o_xx(i,j,2),j=1,m)
end do
time=time+40
close(10)
end if
v_x(:,:,1)=v_x(:,:,2)
v_z(:,:,1)=v_z(:,:,2)
o_xx(:,:,1)=o_xx(:,:,2)
o_zz(:,:,1)=o_zz(:,:,2)
o_xz(:,:,1)=o_xz(:,:,2)
end do !final del ciclo de tiempo
end if !final de la estabilidad
call system ("gnuplot Script_vix.gp")
close(1)
end program
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!.........SUBRUTINA PARA LEER..........................
subroutine lectura(Archivo,Arreglo,z,x)
implicit none
integer::i,j
integer, intent (in)::z,x
double precision, dimension (1:z,1:x), intent(out)::Arreglo
character (len=*), intent (in)::Archivo
open (unit=100,file=Archivo,status='old')
do i=1, z
read(100,*)(Arreglo(i,j),j=1,x)
end do
close(100)
end subroutine lectura