Hi,
I’ve been working on a project for quite some time now and finally have my code compiling. The problem, is that since I work over ssh I don’t have access to a debugger (or at least I don’t know how to use one through ssh). I recently translated my code from straight fortran 90 to cudafortran. It compiles but my simulation doesn’t work properly now (it just exits when the routine I coded in cudafortran gets called). This leads me to believe that my kernel isn’t working. I’m not sure how to debug the kernel, but I was hoping that someone could take a look at it and possibly tell me where I’m going wrong or if there are any obvious issues. My code is below
!!****if* source/physics/Hydro/HydroMain/split/PPM/PPMKernel/rieman
!!
!! NAME
!!
!! rieman
!!
!! SYNOPSIS
!!
!! call rieman ( integer(IN) :: numIntCells,
!! integer(IN) :: numCells,
!! real(OUT) :: rhoav(numCells),
!! real(OUT) :: uav(numCells),
!! real(OUT) :: utav(numCells),
!! real(OUT) :: uttav(numCells),
!! real(OUT) :: pav(numCells),
!! real(OUT) :: urell(numCells),
!! real(IN) :: ugrdl(numCells),
!! real(IN) :: game(numCells),
!! real(OUT) :: gameav(numCells),
!! real(OUT) :: xnav(numCells, hy_numXn),
!! real(IN) :: x(numCells))
!!
!!
!!
!! DESCRIPTION
!!
!! Solve riemann shock tube problem for a general equation of state using
!! the method of Colella and Glaz. Use a two shock approximation, and
!! linearly interpolation between the head and tail of a rarefaction to
!! treat rarefactions.
!!
!! Take as input the effective left and right states, obtained by
!! integrating the parabolic reconstructions of the data over the domain
!! of dependence of each of the characteristics on the respective side
!! of the interface. This is accomplished by states(). Return the
!! solution to Riemann's problem on axis -- this is used in computing
!! the fluxes.
!!
!! The Riemann problem for the Euler's equation produces 4 states,
!! separated by the three characteristics (u - cs, u, u + cs):
!!
!!
!! l_1 t l_2 l_3
!! \ ^ . /
!! \ *L | . *R /
!! \ | . /
!! \ | . /
!! L \ | . / R
!! \ | . /
!! \ |. /
!! \|./
!! ----------+----------------> x
!!
!! l_1 = u - cs eigenvalue
!! l_2 = u eigenvalue (contact)
!! l_3 = u + cs eigenvalue
!!
!! only density jumps across l_2
!!
!! References:
!!
!! CG: Colella & Glaz 1985, JCP, 59, 264.
!!
!! CW: Colella & Woodward 1984, JCP, 54, 174.
!!
!! Fry: Fryxell et al. 2000, ApJS, 131, 273.
!!
!! Toro: Toro 1999, ``Riemann Solvers and Numerical Methods for Fluid
!! Dynamcs: A Practical Introduction, 2nd Ed.'', Springer-Verlag
!!
!! ARGUMENTS
!!
!! numIntCells :
!! numCells :
!! rhoav :
!! uav :
!! utav :
!! uttav :
!! pav :
!! urell :
!! ugrdl :
!! game :
!! gameav :
!! xnav :
!! x :
!!
!!
!!
!!***
module riemanCalc
use cudafor
implicit none
#include "constants.h"
#include "Flash.h"
real, device, allocatable, DIMENSION(:)::urght
real, device, allocatable, DIMENSION(:)::uavDev
real, device, allocatable, DIMENSION(:)::rhoavDev
real, device, allocatable, DIMENSION(:)::utavDev
real, device, allocatable, DIMENSION(:)::uttavDev
real, device, allocatable, DIMENSION(:)::pavDev
real, device, allocatable, DIMENSION(:)::urellDev
real, device, allocatable, DIMENSION(:)::gameavDev
real, device, allocatable, DIMENSION(:)::val1
real, device, allocatable, DIMENSION(:)::val2
real, device, allocatable, DIMENSION(:)::ugrdlDev
real, device, allocatable, DIMENSION(:)::gameDev
real, device, allocatable, DIMENSION(:)::numXn
real, device, allocatable, DIMENSION(:)::gmelft
real, device, allocatable, DIMENSION(:)::gmergt
real, device, allocatable, DIMENSION(:)::plft
real, device, allocatable, DIMENSION(:)::prght
real, device, allocatable, DIMENSION(:)::clft
real, device, allocatable, DIMENSION(:)::crght
real, device, allocatable, DIMENSION(:)::ulft
real, device, allocatable, DIMENSION(:)::vlft
real, device, allocatable, DIMENSION(:)::vrght
real, device, allocatable, DIMENSION(:)::utlft
real, device, allocatable, DIMENSION(:)::uttlft
real, device, allocatable, DIMENSION(:)::utrght
real, device, allocatable, DIMENSION(:)::uttrgt
real, device, allocatable, DIMENSION(:,:)::xnlft
real, device, allocatable, DIMENSION(:,:)::xnrght
real, device, allocatable, DIMENSION(:)::smallp
real, device, allocatable, DIMENSION(:)::smallu
real, device, allocatable, DIMENSION(:)::smlrho
real, device, allocatable, DIMENSION(:)::nriem
real, device, allocatable, DIMENSION(:)::gmclft
real, device, allocatable, DIMENSION(:)::gmcrgt
real, device, allocatable, DIMENSION(:)::pstor
real, device, allocatable, DIMENSION(:)::riemanTol
real, device, allocatable, DIMENSION(:)::errorFlag
real, device, allocatable, DIMENSION(:,:)::xnavDev
real, device, allocatable, DIMENSION(:) :: wlft, wrght, pstar, ustar, vstar, cestar, &
rhostr, westar, ps, us, uts, utts, vs, rhos, ces, ws, wes, &
gmstar, games, gamcs
real, device, allocatable, DIMENSION(:) :: pstar1, pstar2, gmstrl, gmstrr, &
& wlft1, wrght1, gmin, gmax, &
& gamfac, aux
real, device, allocatable, DIMENSION(:) :: scrch11, scrch2, scrch3, scrch4
real, device, allocatable :: ge, gc, ustrl1, ustrr1, ustrl2, ustrr2, &
& delu1, delu2, pressErr,mxVal
real, parameter :: small_dp = 1.e2 * epsilon(1.e0)
contains
subroutine rieman (numIntCells, numCells, &
rhoav, uav, utav, uttav, pav, &
urell, ugrdl, game, gameav, xnav, x)
use Hydro_data, ONLY: hy_numXn, &
hy_gmelft, hy_gmergt, &
hy_plft, hy_prght, &
hy_clft, hy_crght, &
hy_ulft, hy_urght, &
hy_vlft, hy_vrght, &
hy_utlft, hy_utrght, &
hy_uttlft, hy_uttrgt, &
hy_xnlft, hy_xnrght, &
hy_smallp, hy_smallu, &
hy_smlrho, hy_nriem, &
hy_gmclft, hy_gmcrgt,hy_pstor, &
hy_riemanTol
use Driver_interface, ONLY : Driver_abortFlash
implicit none
!! Arguments ----------------------------------
integer, intent (IN) :: numIntCells,numCells
real, intent(IN), DIMENSION(numCells) :: x
real, intent(IN), DIMENSION(numCells) :: ugrdl, game
real, intent(OUT), DIMENSION(numCells) :: uav, rhoav, utav, uttav, pav, &
urell, gameav
real, intent(OUT), DIMENSION(numCells,hy_numXn) :: xnav
integer :: i, j, k, n,ierr,counter
character(len=1), save :: dirs(3) = (/ 'x', 'y', 'z' /)
!! Local variable ---------------------------------
!--------------------------------------------------------------------------
! We carry around two adiabatic indices (gammas). gamc is usually
! referred to as gamma_1, and is d(log P)/d(log rho) (CG Eq. 7). game
! is CG Eq. 8.
! calculate limits on gamma based on the values in the neighboring zones
! gamfac is the max gamma factor in CG Eq. 31, and is used over and over,
! so store it for efficiency.
integer::numIntCells5
real, DIMENSION(numIntCells5)::DummyArray
type(dim3)::dimGrid,dimBlock
numIntCells5=numIntCells+5
val2=numIntCells5
do i=0,numCells-1
DummyArray(i)=0
enddo
allocate(ge)
allocate(gc)
allocate(ustrl1)
allocate(ustrr1)
allocate(ustrl2)
allocate(ustrr2)
allocate(delu1)
allocate(delu2)
allocate(pressErr)
allocate(mxVal)
allocate(pstar1(numCells))
allocate(pstar2(numCells))
allocate(gmstrl(numCells))
allocate(gmstrr(numCells))
allocate(wlft1(numCells))
allocate(wrght1(numCells))
allocate(gmin(numCells))
allocate(gmax(numCells))
allocate(gamfac(numCells))
allocate(aux(numCells))
allocate(scrch11(numCells))
allocate(scrch2(numCells))
allocate(scrch3(numCells))
allocate(scrch4(numCells))
allocate(wlft(numCells))
allocate(wrght(numCells))
allocate(pstar(numCells))
allocate(ustar(numCells))
allocate(vstar(numCells))
allocate(cestar(numCells))
allocate(rhostr(numCells))
allocate(westar(numCells))
allocate(ps(numCells))
allocate(us(numCells))
allocate(uts(numCells))
allocate(utts(numCells))
allocate(vs(numCells))
allocate(rhos(numCells))
allocate(ces(numCells))
allocate(ws(numCells))
allocate(wes(numCells))
allocate(gmstar(numCells))
allocate(games(numCells))
allocate(gamcs(numCells))
allocate(xnavDev(numCells,hy_numXn))
allocate(uavDev(numCells))
allocate(rhoavDev(numCells))
allocate(utavDev(numCells))
allocate(uttavDev(numCells))
allocate(pavDev(numCells))
allocate(urellDev(numCells))
allocate(gameavDev(numCells))
allocate(val1(numCells))
allocate(val2(numCells))
allocate(ugrdlDev(numCells))
allocate(gameDev(numCells))
allocate(numXn(numCells))
allocate(gmelft(numCells))
allocate(gmergt(numCells))
allocate(plft(numCells))
allocate(prght(numCells))
allocate(clft(numCells))
allocate(crght(numCells))
allocate(ulft(numCells))
allocate(vlft(numCells))
allocate(vrght(numCells))
allocate(utlft(numCells))
allocate(uttlft(numCells))
allocate(utrght(numCells))
allocate(uttrgt(numCells))
allocate(xnlft(MAXCELLS,numCells))
allocate(xnrght(MAXCELLS,numCells))
allocate(smallp(numCells))
allocate(smallu(numCells))
allocate(smlrho(numCells))
allocate(nriem(numCells))
allocate(gmclft(numCells))
allocate(gmcrgt(numCells))
allocate(pstor(numCells))
allocate(riemanTol(numCells))
!allocate(errorFlag(numCells))
allocate(urght(numCells))
pstar1=DummyArray(0:numCells - 1)
pstar2=DummyArray(0:numCells - 1)
gmstrl=DummyArray(0:numCells - 1)
gmstrr=DummyArray(0:numCells - 1)
wlft1=DummyArray(0:numCells - 1)
wrght1=DummyArray(0:numCells - 1)
gmin=DummyArray(0:numCells - 1)
gmax=DummyArray(0:numCells - 1)
gamfac=DummyArray(0:numCells - 1)
aux=DummyArray(0:numCells - 1)
scrch11=DummyArray(0:numCells - 1)
scrch2=DummyArray(0:numCells - 1)
scrch3=DummyArray(0:numCells - 1)
scrch4=DummyArray(0:numCells - 1)
wlft=DummyArray(0:numCells - 1)
wrght=DummyArray(0:numCells - 1)
pstar=DummyArray(0:numCells - 1)
ustar=DummyArray(0:numCells - 1)
vstar=DummyArray(0:numCells - 1)
cestar=DummyArray(0:numCells - 1)
rhostr=DummyArray(0:numCells - 1)
westar=DummyArray(0:numCells - 1)
ps=DummyArray(0:numCells - 1)
us=DummyArray(0:numCells - 1)
uts=DummyArray(0:numCells - 1)
utts=DummyArray(0:numCells - 1)
vs=DummyArray(0:numCells - 1)
rhos=DummyArray(0:numCells - 1)
ces=DummyArray(0:numCells - 1)
ws=DummyArray(0:numCells - 1)
wes=DummyArray(0:numCells - 1)
gmstar=DummyArray(0:numCells - 1)
games=DummyArray(0:numCells - 1)
gamcs=DummyArray(0:numCells - 1)
xnavDev=xnav(0:numCells - 1,0:hy_numXn)
utavDev=utav(0:numCells - 1)
uttavDev=uttav(0:numCells - 1)
pavDev=pav(0:numCells - 1)
urellDev=urell(0:numCells - 1)
gameavDev=gameav(0:numCells - 1)
rhoavDev=rhoav(0:numCells - 1)
uavDev=uav(0:numCells - 1)
ugrdlDev=ugrdl(0:numCells - 1)
gameDev=game(0:numCells - 1)
numXn=hy_numXn
gmelft=hy_gmelft(0:numCells - 1)
gmergt=hy_gmergt(0:numCells - 1)
plft=hy_plft(0:numCells - 1)
prght=hy_prght(0:numCells - 1)
clft=hy_clft(0:numCells - 1)
crght=hy_crght(0:numCells - 1)
ulft=hy_ulft(0:numCells - 1)
vlft=hy_vlft(0:numCells - 1)
vrght=hy_vrght(0:numCells - 1)
utlft=hy_utlft(0:numCells - 1)
utrght=hy_utrght(0:numCells - 1)
uttlft=hy_uttlft(0:numCells - 1)
uttrgt=hy_uttrgt(0:numCells - 1)
xnlft=hy_xnlft(0:MAXCELLS,0:numCells - 1)
xnrght=hy_xnrght(0:MAXCELLS,0:numCells - 1)
smallp=hy_smallp
smallu=hy_smallu
smlrho=hy_smlrho
nriem=hy_nriem
gmclft=hy_gmclft(0:numCells - 1)
gmcrgt=hy_gmcrgt(0:numCells - 1)
pstor=hy_pstor(0:numCells - 1)
riemanTol=hy_riemanTol
urght=hy_urght(0:numCells - 1)
dimGrid=dim3(numIntCells5/16,numIntCells5/16,1)
dimBlock=dim3(16,16,1)
call rieman_kernel<<<dimGrid,dimBlock>>>()
uav=uavDev(0:numCells - 1)
rhoav=rhoavDev(0:numCells - 1)
utav=utavDev(0:numCells - 1)
uttav=uttavDev(0:numCells - 1)
pav=pavDev(0:numCells - 1)
urell=urellDev(0:numCells - 1)
gameav=gameavDev(0:numCells - 1)
xnav=xnavDev(0:numCells - 1-1,0:hy_numXn-1)
deallocate(ge)
deallocate(gc)
deallocate(ustrl1)
deallocate(ustrr1)
deallocate(ustrl2)
deallocate(ustrr2)
deallocate(delu1)
deallocate(delu2)
deallocate(pressErr)
deallocate(mxVal)
deallocate(pstar1)
deallocate(pstar2)
deallocate(gmstrl)
deallocate(gmstrr)
deallocate(wlft1)
deallocate(wrght1)
deallocate(gmin)
deallocate(gmax)
deallocate(gamfac)
deallocate(aux)
deallocate(scrch11)
deallocate(scrch2)
deallocate(scrch3)
deallocate(scrch4)
deallocate(wlft)
deallocate(wrght)
deallocate(pstar)
deallocate(ustar)
deallocate(vstar)
deallocate(cestar)
deallocate(rhostr)
deallocate(westar)
deallocate(ps)
deallocate(us)
deallocate(uts)
deallocate(utts)
deallocate(vs)
deallocate(rhos)
deallocate(ces)
deallocate(ws)
deallocate(wes)
deallocate(gmstar)
deallocate(games)
deallocate(gamcs)
deallocate(xnavDev)
deallocate(uavDev)
deallocate(rhoavDev)
deallocate(utavDev)
deallocate(uttavDev)
deallocate(pavDev)
deallocate(urellDev)
deallocate(gameavDev)
deallocate(val1)
deallocate(val2)
deallocate(ugrdlDev)
deallocate(gameDev)
deallocate(numXn)
deallocate(gmelft)
deallocate(gmergt)
deallocate(plft)
deallocate(prght)
deallocate(clft)
deallocate(crght)
deallocate(ulft)
deallocate(vlft)
deallocate(vrght)
deallocate(utlft)
deallocate(uttlft)
deallocate(utrght)
deallocate(uttrgt)
deallocate(xnlft)
deallocate(xnrght)
deallocate(smallp)
deallocate(smallu)
deallocate(smlrho)
deallocate(nriem)
deallocate(gmclft)
deallocate(gmcrgt)
deallocate(pstor)
deallocate(riemanTol)
!allocate(errorFlag)
deallocate(urght)
return
end subroutine rieman
attributes (global) subroutine rieman_kernel()
implicit none
integer :: i, n
do i = 5, val2(1)
aux(i) = sqrt (0.5e0 * (gameDev(i) - 1.0e0) / gameDev(i))
ge = 0.5e0 * (gmelft(i) + gmergt(i))
gc = 0.5e0 * (gmclft(i) + gmcrgt(i))
gamfac(i) = (1.e0 - ge / gc) * (ge - 1.e0)
gmin(i) = min (gameDev(i-1), gameDev(i), gameDev(i+1))
gmax(i) = max (gameDev(i-1), gameDev(i), gameDev(i+1))
enddo
! construct first guess for secant iteration by assuming that the nonlinear
! wave speed is equal to the sound speed -- the resulting expression is the
! same as Toro, Eq. 9.28 in the Primitive Variable Riemann Solver (PVRS).
! See also Fry Eq. 72.
do i = 5, val2(1)
pstar1(i) = prght(i) - plft(i) -crght(i) * (urght(i) - ulft(i))
pstar1(i) = plft(i) + pstar1(i) * (clft(i) / (clft(i) + crght(i)))
pstar1(i) = max (smallp(1), pstar1(i))
enddo
! calculate approximation jump in gamma acrosss the interface based on the
! first guess for the pressure jump. There is a left and right 'star' region,
! so we need gamma add both places. Use CG Eq. 31 and 32, with definitions
! as in CG Eq. 33.
do i = 5, val2(1)
gmstrl(i) = gamfac(i) * (pstar1(i) - plft(i))
gmstrl(i) = gmelft(i) + 2.e0 * gmstrl(i) / (pstar1(i) + plft(i))
gmstrr(i) = gamfac(i) * (pstar1(i) - prght(i))
gmstrr(i) = gmergt(i) + 2.e0 * gmstrr(i) / (pstar1(i) + prght(i))
gmstrl(i) = max (gmin(i), min( gmstrl(i), gmax(i)))
gmstrr(i) = max (gmin(i), min( gmstrr(i), gmax(i)))
enddo
! calculate nonlinear wave speeds for the left and right moving waves based
! on the first guess for the pressure jump. Again, there is a left and a
! right wave speed. Compute this using CG Eq. 34.
do i = 5, val2(1)
scrch11(i) = pstar1(i) - (gmstrl(i) - 1.e0) * plft(i) &
& / (gmelft(i) - 1.e0)
if (scrch11(i) .EQ. 0.e0) scrch11(i) = smallp(1)
wlft1(i) = pstar1(i) + 0.5e0 * (gmstrl(i) - 1.e0) &
& * (pstar1(i) + plft(i))
wlft1(i) = (pstar1(i) - plft(i)) * wlft1(i) / (vlft(i) * scrch11(i))
wlft1(i) = sqrt(abs(wlft1(i)))
scrch2(i) = pstar1(i) - (gmstrr(i) - 1.e0) * prght(i) /(gmergt(i) - 1.e0)
if (scrch2(i) .EQ. 0.e0) scrch2(i) = smallp(1)
wrght1(i) = pstar1(i) + 0.5e0 * (gmstrr(i) - 1.e0) &
& * (pstar1(i) + prght(i))
wrght1(i) = (pstar1(i) - prght(i)) * wrght1(i) / (vrght(i) * scrch2(i))
wrght1(i) = sqrt(abs(wrght1(i)))
! if the pressure jump is small, the wave speed is just the sound speed
if (abs (pstar1(i) - plft(i)) < small_dp*(pstar1(i) + plft(i))) wlft1(i) = clft(i)
wlft1(i) = max (wlft1(i), aux(i) * clft(i))
if (abs (pstar1(i) - prght(i)) < small_dp*((pstar1(i) + prght(i)))) wrght1(i) = crght(i)
wrght1(i) = max (wrght1(i), aux(i) * crght(i))
enddo
! construct second guess for the pressure using the nonlinear wave speeds
! from the first guess. This is basically the same thing we did to get
! pstar1, except now we are using the better wave speeds instead of the
! sound speed.
do i = 5, val2(1)
pstar2(i) = prght(i) - plft(i) - wrght1(i) * (urght(i) - ulft(i))
pstar2(i) = plft(i) + pstar2(i) * wlft1(i) / (wlft1(i) + wrght1(i))
pstar2(i) = max (smallp(1), pstar2(i))
enddo
! begin the secant iteration -- see CG Eq. 17 for details. We will continue to
! interate for convergence until the error falls below tol (in which case,
! things are good), or we hit hy_nriem iterations (in which case we have a
! problem, and we spit out an error).
do i = 5, val2(1)
pstor(1) = pstar1(i)
pstor(2) = pstar2(i)
do n = 1, nriem(1)
! new values for the gamma at the "star" state -- again, using CG Eq. 31
gmstrl(i) = gamfac(i) * (pstar2(i) - plft(i))
gmstrl(i) = gmelft(i) + 2.e0 * gmstrl(i) / (pstar2(i) + plft(i))
gmstrr(i) = gamfac(i) * (pstar2(i) - prght(i))
gmstrr(i) = gmergt(i) + 2.e0 * gmstrr(i) / (pstar2(i) + prght(i))
gmstrl(i) = max (gmin(i), min (gmax(i), gmstrl(i)))
gmstrr(i) = max (gmin(i), min (gmax(i), gmstrr(i)))
! new nonlinear wave speeds, using CG Eq. 34 and the updated gammas
scrch11(i) = pstar2(i) - (gmstrl(i) - 1.e0) * plft(i) &
& / (gmelft(i) - 1.e0)
if (scrch11(i) .EQ. 0.e0) scrch11(i) = smallp(1)
wlft(i) = pstar2(i) + 0.5e0 * (gmstrl(i) - 1.e0) &
& * (pstar2(i) + plft(i))
wlft(i) = (pstar2(i) - plft(i)) * wlft(i) / (vlft(i) * scrch11(i))
wlft(i) = sqrt(abs(wlft(i)))
scrch2(i) = pstar2(i) - (gmstrr(i) - 1.e0) * prght(i) /(gmergt(i) - 1.e0)
if (scrch2(i) .EQ. 0.e0) scrch2(i) = smallp(1)
wrght(i) = pstar2(i) + 0.5e0 * (gmstrr(i) - 1.e0) &
& * (pstar2(i) + prght(i))
wrght(i) = (pstar2(i) - prght(i)) * wrght(i) / (vrght(i) * scrch2(i))
wrght(i) = sqrt(abs(wrght(i)))
! if the pressure jump is small, the wave speed is just the sound speed
if (abs (pstar2(i) - plft(i)) < small_dp*(pstar2(i) + plft(i))) wlft(i) = clft(i)
wlft(i) = max (wlft(i), aux(i) * clft(i))
if (abs (pstar2(i) - prght(i)) < small_dp*(pstar2(i) + prght(i))) wrght(i) =crght(i)
wrght(i) = max (wrght(i), aux(i) * crght(i))
! compute the velocities in the "star" state -- using CG Eq. 18 -- ustrl2 and
! ustrr2 are the velocities they define there. ustrl1 and ustrl2 seem to be
! the velocities at the last time, since pstar1 is the old 'star' pressure, and
! wlft1 is the old wave speed.
ustrl1 = ulft(i) - (pstar1(i) - plft(i)) / wlft1(i)
ustrr1 = urght(i) + (pstar1(i) - prght(i)) / wrght1(i)
ustrl2 = ulft(i) - (pstar2(i) - plft(i)) / wlft(i)
ustrr2 = urght(i) + (pstar2(i) - prght(i)) / wrght(i)
delu1 = ustrl1 - ustrr1
delu2 = ustrl2 - ustrr2
scrch11(i) = delu2 - delu1
if (abs(pstar2(i)-pstar1(i)) .le. smallp(1)) scrch11(i) = 0.e0
if (abs(scrch11(i)) .lt. smallu(1)) then
delu2 = 0.e0
scrch11(i) = 1.e0
endif
! pressure at the "star" state -- using CG Eq. 18
pstar(i) = pstar2(i) - delu2 * (pstar2(i) - pstar1(i)) / scrch11(i)
pstar(i) = max (smallp(1), pstar(i))
! check for convergence of iteration, hy_riemanTol is a run-time parameter
pressErr = abs(pstar(i)-pstar2(i)) / pstar(i)
if (pressErr .lt. riemanTol(1)) goto 10
! reset variables for next iteration
pstar1(i) = pstar2(i)
pstar2(i) = pstar(i)
pstor(n+2) = pstar(i)
wlft1(i) = wlft(i)
wrght1(i) = wrght(i)
enddo
n = n - 1
! print error message and stop code if iteration fails to converge
! print *, ' '
! print *, 'Nonconvergence in subroutine rieman'
! print *, ' '
! print *, 'Zone index = ', i
! print *, 'Zone center = ', x(i)
!print *, 'Iterations tried = ', n+2
! print *, 'Pressure error = ', pres_err
! print *, 'rieman_tol = ', hy_riemanTol
! print *, ' '
! print *, 'pL = ', hy_plft(i), ' pR =', hy_prght(i)
! print *, 'uL = ', hy_ulft(i), ' uR =', hy_urght(i)
! print *, 'cL = ', hy_clft(i), ' cR =', hy_crght(i)
! print *, 'gamma_eL = ', hy_gmelft(i), ' gamma_eR =', hy_gmergt(i)
! print *, 'gamma_cL = ', hy_gmclft(i), ' gamma_cR =', hy_gmcrgt(i)
! print *, ' '
! print *, 'Iteration history:'
! print *, ' '
! print '(A4, 2X, A20)', 'n', 'p*'
! do j = 1, n+2
! print '(I4, 2X, E20.12)', j, hy_pstor(j)
! enddo
! print *, ' '
! print *, 'Terminating execution.'
! land here if the iterations have converged
return
10 continue
enddo
! end of secant iteration
! calculate fluid velocity for the "star" state -- this comes from the shock
! jump equations, Fry Eq. 68 and 69. The ustar velocity can be computed
! using either the jump eq. for a left moving or right moving shock -- we use
! the average of the two.
do i = 5, val2(1)
scrch3(i) = ulft (i) - (pstar(i) - plft(i)) / wlft(i)
scrch4(i) = urght(i) + (pstar(i) - prght(i)) / wrght(i)
ustar(i) = 0.5e0 * (scrch3(i) + scrch4(i))
enddo
! account for grid velocity
do i = 5, val2(1)
urellDev(i) = ustar(i) - ugrdlDev(i)
scrch11(i) = sign (1.e0, urellDev(i))
enddo
! decide which state is located at the zone iterface based on the values
! of the wave speeds. This is just saying that if ustar > 0, then the state
! is U_L. if ustar < 0, then the state on the axis is U_R.
do i = 5, val2(1)
scrch2(i) = 0.5e0 * ( 1.e0 + scrch11(i))
scrch3(i) = 0.5e0 * ( 1.e0 - scrch11(i))
ps(i) = plft(i) * scrch2(i) + prght(i) * scrch3(i)
us(i) = ulft(i) * scrch2(i) + urght(i) * scrch3(i)
uts(i) =utlft(i) * scrch2(i) + utrght(i) * scrch3(i)
utts(i) =uttlft(i) * scrch2(i) + uttrgt(i) * scrch3(i)
vs(i) = vlft(i) * scrch2(i) + vrght(i) * scrch3(i)
games(i) = gmelft(i) * scrch2(i) + gmergt(i) * scrch3(i)
gamcs(i) = gmclft(i) * scrch2(i) + gmcrgt(i) * scrch3(i)
rhos(i) = 1.e0 / vs(i)
rhos(i) = max (smlrho(1), rhos(i))
vs(i) = 1.e0 / rhos(i)
ws(i) = wlft(i) * scrch2(i) + wrght(i) * scrch3(i)
ces(i) = sqrt (gamcs(i) * ps(i) * vs(i))
! compute rhostar, using the shock jump condition (Fry Eq. 80)
vstar(i) = vs(i) - (pstar(i) - ps(i)) / ws(i) / ws(i)
rhostr(i) = 1.e0 / vstar(i)
cestar(i) = sqrt (gamcs(i) * pstar(i) * vstar(i))
! compute some factors, Fry Eq. 81 and 82
wes(i) = ces(i) - scrch11(i) * us(i)
westar(i) = cestar(i) - scrch11(i) * ustar(i)
scrch4(i) = ws(i) * vs(i) - scrch11(i) * us(i)
if (pstar(i) - ps(i) .ge. 0.e0) then
wes(i) = scrch4(i)
westar(i) = scrch4(i)
endif
wes(i) = wes(i) + scrch11(i) * ugrdlDev(i)
westar(i) = westar(i) + scrch11(i) * ugrdlDev(i)
enddo
do i = 5, val2(1)
gamfac(i) = (1.e0 - games(i) / gamcs(i)) * (games(i) - 1.e0)
gmstar(i) = gamfac(i) * (pstar(i) - ps(i))
gmstar(i) = games(i) + 2.e0 * gmstar(i) / (pstar(i) + ps(i))
gmstar(i) = max (gmin(i), min (gmax(i), gmstar(i)))
enddo
do n = 1, numXn(1)
do i = 5, val2(1)
xnavDev(i,n) =xnlft(i,n) * scrch2(i) + xnrght(i,n) * scrch3(i)
enddo
enddo
! compute correct state for rarefaction fan by linear interpolation
do i = 5, val2(1)
scrch11(i) = max (wes(i) - westar(i), wes(i) + westar(i), smallu(i))
scrch11(i) = (wes(i) + westar(i)) / scrch11(i)
scrch11(i) = 0.5e0 * (1.e0 + scrch11(i))
scrch2(i) = 1.e0 - scrch11(i)
rhoavDev(i) = scrch11(i) * rhostr(i) + scrch2(i) * rhos (i)
uavDev (i) = scrch11(i) * ustar(i) + scrch2(i) * us(i)
utavDev (i) = uts(i)
uttavDev(i) = utts(i)
pavDev(i) = scrch11(i) * pstar(i) + scrch2(i) * ps(i)
gameavDev(i) = scrch11(i) * gmstar(i) + scrch2(i) * games(i)
enddo
do i = 5, val2(1)
if (westar(i) .ge. 0.e0) then
rhoavDev(i) = rhostr(i)
uavDev(i) = ustar(i)
pavDev(i) = pstar(i)
gameavDev(i) = gmstar(i)
endif
if (wes(i) .lt. 0.e0) then
rhoavDev(i) = rhos(i)
uavDev(i) = us(i)
pavDev(i) = ps(i)
gameavDev(i) = games(i)
endif
urellDev(i) = uavDev(i) - ugrdlDev(i)
enddo
return
end subroutine rieman_kernel
end module riemanCalc
[\code]