> Hello,
>
> I've run into a segmentation fault runtime error with a fortran code
> that I've been writing, which I no doubt think has to do with some
> out-of-bounds arrays I've declared and used, but I'm not certain of
> where my mistake lies. I was hoping someone here more proficient at
> programming might be able to spot and point out what my mistake has
> been and proffer a short explanation if they have some time. I've
> included my program below. I'd appreciate any help any of you may have
> to offer.
>
> Thanks.
>
>
>
> program prandtl meyer
>
> implicit none
>
> integer :: i, j
> integer, parameter :: imax = 100
> integer, parameter :: jmax = 40
>
> real :: pi
> real, parameter :: cfl = 0.5
> real, parameter :: cy = 0.6
> real, parameter :: gam = 1.4
> real, parameter :: R = 287.
> real, parameter :: H1 = 40.
> real, parameter :: L = 65.
> real, parameter :: e = 10.
> real, parameter :: theta = 5.352
> real, dimension(0:jmax,0:imax) :: x, h, u, v, rho, p, temp, mach
> real, dimension(0:jmax,0:imax) :: f1, f2, f3, f4, g1, g2, g3, g4
> real, dimension(0:jmax,0:imax) :: y, ys, eta, deta_dx, a, b, c
> real, dimension(0:jmax,0:imax) :: rho_bar, f1_bar, f2_bar, f3_bar
> real, dimension(0:jmax,0:imax) :: f4_bar, g1_bar, g2_bar, g3_bar
> real, dimension(0:jmax,0:imax) :: g4_bar, mu, xi_step, df1_dxi
> real, dimension(0:jmax,0:imax) :: df2_dxi, df3_dxi, sf1, sf2, sf3
> real, dimension(0:jmax,0:imax) :: sf4, a_bar, b_bar, c_bar
> real, dimension(0:jmax,0:imax) :: u_bar, v_bar, p_bar, temp_bar
> real, dimension(0:jmax,0:imax) :: df4_dxi, df1_bar_dxi
> real, dimension(0:jmax,0:imax) :: df2_bar_dxi, df3_bar_dxi
> real, dimension(0:jmax,0:imax) :: df4_bar_dxi, df1_dxi_av
> real, dimension(0:jmax,0:imax) :: df2_dxi_av, df3_dxi_av
> real, dimension(0:jmax,0:imax) :: df4_dxi_av, sf1_bar
> real, dimension(0:jmax,0:imax) :: sf2_bar, sf3_bar, sf4_bar
> real, dimension(0:jmax,0:imax) :: thetam, thetap, den
>
> ! Specify initial geometry and flow conditions.
>
> ! write(*,*)'=======================================================
> ! &====================================='
> ! write(*,*)' j u, m/s v, m/s rho, kg/s p,
> ! & N/m**2 T, K M'
> ! write(*,*)'=======================================================
> ! &====================================='
>
> pi = 4.*atan( 1.0 )
>
> do i = 0, imax
> x(j,i) = L/imax*i

> do j = 0, jmax
> y(j,i) = H1/jmax*j
> if ( x(j,i) >= 0. .and. x(j,i) <= e ) then
> h(j,i) = 40.
> ys(j,i) = 0.
> eta(j,i) = ( y(j,i) - ys(j,i) )/h(j,i)
> deta_dx(j,i) = 0.
> else
> h(j,i) = H1 + ( x(j,i) - e )*tan( theta*pi/180. )
> ys(j,i) = -( x(j,i) - e )*tan( theta*pi/180. )
> eta(j,i) = ( y(j,i) - ys(j,i) )/h(j,i)
> deta_dx(j,i) = ( 1 - eta(j,i) )*tan( theta*pi/180. )/
> & h(j,i)
> end if
> p(j,i) = 1.01e+05
> rho(j,i) = 1.23
> temp(j,i) = 286.1
> mach(j,i) = 2.
> u(j,i) = mach(j,i)*sqrt( gam*R*temp(j,i) )
> v(j,i) = 0.
> f1(j,i) = rho(j,i)*u(j,i)
> f2(j,i) = p(j,i) + rho(j,i)*u(j,i)**2
> f3(j,i) = rho(j,i)*u(j,i)*v(j,i)
> f4(j,i) = gam/( gam - 1. )*p(j,i)*u(j,i) + rho(j,i)*
> & u(j,i)*( u(j,i) + v(j,i) )/2.
> g1(j,i) = rho(j,i)*f3(j,i)/f1(j,i)
> g2(j,i) = f3(j,i)
> g3(j,i) = rho(j,i)*( f3(j,i)/f1(j,i) )**2 + f2(j,i) -
> & ( f1(j,i)**2 )/rho(j,i)
> g4(j,i) = gam/( gam - 1. )*( f2(j,i) - ( f1(j,i)**2 )/
> & rho(j,i) )*f3(j,i)/f1(j,i) + rho(j,i)/2.*f3(j,i)/f1(j,i)
> & *( ( f1(j,i)/rho(j,i) )**2 + ( f3(j,i)/f1(j,i) )**2 )
> end do
>
> end do
>
> ! do j = 0, jmax
> ! write(*,5) j, u(j,0), v(j,0), rho(j,0), p(j,0),
> ! & temp(j,0), mach(j,0)
> ! 5 format(1x,i2,5x,6(e10.3,5x))
> ! end do
>
> write(*,10) x(1,15), y(1,15), h(1,15), ys(1,15), eta(1,15),
> & deta_dx(1,15)
> 10 format(1x,6(f6.3,3x))
>
> do i = 0, imax-1
>
> ! Begin predictor step for Maccormack's predictor-corrector
> ! technique for 2-D Prandtl-Meyer expansion wave flow.
>
> do j = 1, jmax-1
> df1_dxi(j,i) = deta_dx(j,i)*( f1(j,i) - f1(j+1,i) )/
> & ( eta(j,i) - eta(j+1,i) ) + 1./h(j,i)*( g1(j,i) -
> & g1(j+1,i) )/( eta(j,i) - eta(j+1,i) )
>
> df2_dxi(j,i) = deta_dx(j,i)*( f2(j,i) - f2(j+1,i) )/
> & ( eta(j,i) - eta(j+1,i) ) + 1./h(j,i)*( g2(j,i) -
> & g2(j+1,i) )/( eta(j,i) - eta(j+1,i) )
>
> df3_dxi(j,i) = deta_dx(j,i)*( f3(j,i) - f3(j+1,i) )/
> & ( eta(j,i) - eta(j+1,i) ) + 1./h(j,i)*( g3(j,i) -
> & g3(j+1,i) )/( eta(j,i) - eta(j+1,i) )

> sf1(j,i) = cy*( abs( p(j+1,i) - 2.*p(j,i) + p(j-1,i) ) )/
> & ( p(j+1,i) + 2.*p(j,i) + p(j-1,i) )*( f1(j+1,i) - 2.*
> & f1(j,i) + f1(j-1,i) )
>
> sf2(j,i) = cy*( abs( p(j+1,i) - 2.*p(j,i) + p(j-1,i) ) )/
> & ( p(j+1,i) + 2.*p(j,i) + p(j-1,i) )*( f2(j+1,i) - 2.*
> & f2(j,i) + f2(j-1,i) )
>
> sf3(j,i) = cy*( abs( p(j+1,i) - 2.*p(j,i) + p(j-1,i) ) )/
> & ( p(j+1,i) + 2.*p(j,i) + p(j-1,i) )*( f3(j+1,i) - 2.*
> & f3(j,i) + f3(j-1,i) )
>
> sf4(j,i) = cy*( abs( p(j+1,i) - 2.*p(j,i) + p(j-1,i) ) )/
> & ( p(j+1,i) + 2.*p(j,i) + p(j-1,i) )*( f4(j+1,i) - 2.*
> & f4(j,i) + f4(j-1,i) )
>
> mu(j,i) = asin( 1./mach(j,i) )
>
> thetap(j,i) = tan( ( theta + mu(j,i) )*pi/180. )
>
> thetam(j,i) = tan( ( theta - mu(j,i) )*pi/180. )
>
> den(j,i) = max( thetap(j,i), thetam(j,i) )
>
> xi_step(j,i) = cfl*( y(j,i) - y(j+1,i) )/den(j,i)
>
> f1_bar(j,i+1) = f1(j,i) + df1_dxi(j,i)*xi_step(j,i) +
> & sf1(j,i)
>
> f2_bar(j,i+1) = f2(j,i) + df2_dxi(j,i)*xi_step(j,i) +
> & sf2(j,i)
>
> f3_bar(j,i+1) = f3(j,i) + df3_dxi(j,i)*xi_step(j,i) +
> & sf3(j,i)
>
> f4_bar(j,i+1) = f4(j,i) + df4_dxi(j,i)*xi_step(j,i) +
> & sf4(j,i)
>
> a_bar(j,i+1) = ( f3_bar(j,i+1)**2 )/( 2.*f1_bar(j,i) ) -
> & f4_bar(j,i+1)
>
> b_bar(j,i+1) = gam/( gam - 1. )*f1_bar(j,i+1)*f2_bar(j,i)
>
> c_bar(j,i+1) = -( gam + 1. )/( 2.*( gam - 1. ) )*
> & f1_bar(j,i+1)**3
>
> rho_bar(j,i+1) = ( -b(j,i+1) + sqrt( b(j,i+1)**2 - 4*
> & a(j,i+1)*c(j,i+1) ) )/( 2.*a(j,i+1) )
>
> u_bar(j,i+1) = f1_bar(j,i+1)/rho_bar(j,i+1)
>
> v_bar(j,i+1) = f3_bar(j,i+1)/f1_bar(j,i+1)
>
> p_bar(j,i+1) = f2_bar(j,i+1) - f1_bar(j,i+1)*u_bar(j,i+1)
>
> temp_bar(j,i+1) = p_bar(j,i+1)/( rho_bar(j,i+1)*R )
>
> g1_bar(j,i+1) = rho_bar(j,i+1)*f3_bar(j,i+1)/f1_bar(j,i+1)
>
> g2_bar(j,i+1) = f3_bar(j,i+1)
>
> g3_bar(j,i+1) = rho_bar(j,i+1)*( f3_bar(j,i+1)/f1_bar(j,i+1)
> & )**2 + f2_bar(j,i+1) - ( f1_bar(j,i+1)**2 )/rho_bar(j,i+1)
>
> g4_bar(j,i+1) = gam/( gam - 1. )*( f2_bar(j,i+1) -
> & ( f1_bar(j,i+1)**2 )/rho_bar(j,i+1) )*f3_bar(j,i+1)/
> & f1_bar(j,i+1) + rho_bar(j,i+1)/2.*( f3_bar(j,i+1)/
> & f1_bar(j,i+1) )*( ( f1_bar(j,i+1)/rho_bar(j,i+1) )**2 +
> & ( f3_bar(j,i+1)/f1_bar(j,i+1) )**2 )
> end do
>
> ! Begin corrector step for Maccormack's predictor-corrector
> ! technique for 2-D Prandtl-Meyer expansion wave flow.
>
> do j = 1, jmax-1
> df1_bar_dxi(j,i+1) = deta_dx(j,i+1)*( f1_bar(j-1,i+1) -
> & f1_bar(j,i+1) )/( eta(j-1,i+1) - eta(j,i+1) ) + 1./h(j,i+1)
> & *( g1_bar(j-1,i+1) - g1_bar(j,i+1) )/( eta(j-1,i+1) -
> & eta(j,i+1) )
>
> df2_bar_dxi(j,i+1) = deta_dx(j,i+1)*( f2_bar(j-1,i+1) -
> & f2_bar(j,i+1) )/( eta(j-1,i+1) - eta(j,i+1) ) + 1./h(j,i+1)
> & *( g2_bar(j-1,i+1) - g2_bar(j,i+1) )/( eta(j-1,i+1) -
> & eta(j,i+1) )
>
> df3_bar_dxi(j,i+1) = deta_dx(j,i+1)*( f3_bar(j-1,i+1) -
> & f3_bar(j,i+1) )/( eta(j-1,i+1) - eta(j,i+1) ) + 1./h(j,i+1)
> & *( g3_bar(j-1,i+1) - g3_bar(j,i+1) )/( eta(j-1,i+1) -
> & eta(j,i+1) )
>
> df4_bar_dxi(j,i+1) = deta_dx(j,i+1)*( f4_bar(j-1,i+1) -
> & f4_bar(j,i+1) )/( eta(j-1,i+1) - eta(j,i+1) ) + 1./h(j,i+1)
> & *( g4_bar(j-1,i+1) - g4_bar(j,i+1) )/( eta(j-1,i+1) -
> & eta(j,i+1) )
>
> df1_dxi_av(j,i) = 0.5*( df1_dxi(j,i) + df1_bar_dxi(j,i+1) )
>
> df2_dxi_av(j,i) = 0.5*( df2_dxi(j,i) + df2_bar_dxi(j,i+1) )
>
> df3_dxi_av(j,i) = 0.5*( df3_dxi(j,i) + df3_bar_dxi(j,i+1) )
>
> df4_dxi_av(j,i) = 0.5*( df4_dxi(j,i) + df4_bar_dxi(j,i+1) )
>
> sf1_bar(j,i+1) = cy*( abs( p_bar(j+1,i+1) - 2.*p_bar(j,i+1)
> & + p_bar(j-1,i+1) ) )/( p_bar(j+1,i+1) + 2.*p_bar(j,i+1) +
> & p_bar(j-1,i+1) )*( f1_bar(j+1,i+1) - 2.*f1_bar(j,i+1) +
> & f1_bar(j-1,i+1) )
>
> sf2_bar(j,i+1) = cy*( abs( p_bar(j+1,i+1) - 2.*p_bar(j,i+1)
> & + p_bar(j-1,i+1) ) )/( p_bar(j+1,i+1) + 2.*p_bar(j,i+1) +
> & p_bar(j-1,i+1) )*( f2_bar(j+1,i+1) - 2.*f2_bar(j,i+1) +
> & f2_bar(j-1,i+1) )
>
> sf3_bar(j,i+1) = cy*( abs( p_bar(j+1,i+1) - 2.*p_bar(j,i+1)
> & + p_bar(j-1,i+1) ) )/( p_bar(j+1,i+1) + 2.*p_bar(j,i+1) +
> & p_bar(j-1,i+1) )*( f3_bar(j+1,i+1) - 2.*f3_bar(j,i+1) +
> & f3_bar(j-1,i+1) )
>
> sf4_bar(j,i+1) = cy*( abs( p_bar(j+1,i+1) - 2.*p_bar(j,i+1)
> & + p_bar(j-1,i+1) ) )/( p_bar(j+1,i+1) + 2.*p_bar(j,i+1) +
> & p_bar(j-1,i+1) )*( f4_bar(j+1,i+1) - 2.*f4_bar(j,i+1) +
> & f4_bar(j-1,i+1) )
>
> f1(j,i+1) = f1(j,i) + df1_dxi_av(j,i)*xi_step(j,i) +
> & sf1_bar(j,i+1)
>
> f2(j,i+1) = f2(j,i) + df2_dxi_av(j,i)*xi_step(j,i) +
> & sf2_bar(j,i+1)
>
> f3(j,i+1) = f3(j,i) + df3_dxi_av(j,i)*xi_step(j,i) +
> & sf3_bar(j,i+1)
>
> f4(j,i+1) = f4(j,i) + df4_dxi_av(j,i)*xi_step(j,i) +
> & sf4_bar(j,i+1)
>
> a(j,i+1) = ( f3(j,i+1)**2 )/( 2.*f1(j,i+1) ) - f4(j,i+1)
>
> b(j,i+1) = gam/( gam - 1. )*f1(j,i+1)*f2(j,i+1)
>
> c(j,i+1) = -( gam + 1. )/( 2.*( gam - 1. ) )*f1(j,i+1)**3
>
> rho(j,i+1) = ( -b(j,i+1) + sqrt( b(j,i+1)**2 - 4.*a(j,i+1)*
> & c(j,i+1) ) )/( 2.*a(j,i+1) )
>
> u(j,i+1) = f1(j,i+1)/rho(j,i+1)
>
> v(j,i+1) = f3(j,i+1)/f1(j,i+1)
>
> p(j,i+1) = f2(j,i+1) - f1(j,i+1)*u(j,i+1)
>
> temp(j,i+1) = p(j,i+1)/( rho(j,i+1)*R )
> end do
>
> end do
>
> stop
>
> end program prandtl meyer