Google Git
Sign in
gem5 / public / gem5-resources / a534a55c08210e1e116a542985d4fabf05b8750b / . / src / npb / disk-image / npb / npb-hooks / NPB3.3.1 / NPB3.3-MPI / SP / rhs.f
blob: b687f12f83614cfbbd75f1838e1abaf82cea9080 [file] [log] [blame]
c---------------------------------------------------------------------
c---------------------------------------------------------------------
subroutine compute_rhs
c---------------------------------------------------------------------
c---------------------------------------------------------------------
include 'header.h'
integer c, i, j, k, m
double precision aux, rho_inv, uijk, up1, um1, vijk, vp1, vm1,
> wijk, wp1, wm1
if (timeron) call timer_start(t_rhs)
c---------------------------------------------------------------------
c loop over all cells owned by this node
c---------------------------------------------------------------------
do c = 1, ncells
c---------------------------------------------------------------------
c compute the reciprocal of density, and the kinetic energy,
c and the speed of sound.
c---------------------------------------------------------------------
do k = -1, cell_size(3,c)
do j = -1, cell_size(2,c)
do i = -1, cell_size(1,c)
rho_inv = 1.0d0/u(i,j,k,1,c)
rho_i(i,j,k,c) = rho_inv
us(i,j,k,c) = u(i,j,k,2,c) * rho_inv
vs(i,j,k,c) = u(i,j,k,3,c) * rho_inv
ws(i,j,k,c) = u(i,j,k,4,c) * rho_inv
square(i,j,k,c) = 0.5d0* (
> u(i,j,k,2,c)*u(i,j,k,2,c) +
> u(i,j,k,3,c)*u(i,j,k,3,c) +
> u(i,j,k,4,c)*u(i,j,k,4,c) ) * rho_inv
qs(i,j,k,c) = square(i,j,k,c) * rho_inv
c---------------------------------------------------------------------
c (don't need speed and ainx until the lhs computation)
c---------------------------------------------------------------------
aux = c1c2*rho_inv* (u(i,j,k,5,c) - square(i,j,k,c))
aux = dsqrt(aux)
speed(i,j,k,c) = aux
ainv(i,j,k,c) = 1.0d0/aux
end do
end do
end do
c---------------------------------------------------------------------
c copy the exact forcing term to the right hand side; because
c this forcing term is known, we can store it on the whole of every
c cell, including the boundary
c---------------------------------------------------------------------
do m = 1, 5
do k = 0, cell_size(3,c)-1
do j = 0, cell_size(2,c)-1
do i = 0, cell_size(1,c)-1
rhs(i,j,k,m,c) = forcing(i,j,k,m,c)
end do
end do
end do
end do
c---------------------------------------------------------------------
c compute xi-direction fluxes
c---------------------------------------------------------------------
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do j = start(2,c), cell_size(2,c)-end(2,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
uijk = us(i,j,k,c)
up1 = us(i+1,j,k,c)
um1 = us(i-1,j,k,c)
rhs(i,j,k,1,c) = rhs(i,j,k,1,c) + dx1tx1 *
> (u(i+1,j,k,1,c) - 2.0d0*u(i,j,k,1,c) +
> u(i-1,j,k,1,c)) -
> tx2 * (u(i+1,j,k,2,c) - u(i-1,j,k,2,c))
rhs(i,j,k,2,c) = rhs(i,j,k,2,c) + dx2tx1 *
> (u(i+1,j,k,2,c) - 2.0d0*u(i,j,k,2,c) +
> u(i-1,j,k,2,c)) +
> xxcon2*con43 * (up1 - 2.0d0*uijk + um1) -
> tx2 * (u(i+1,j,k,2,c)*up1 -
> u(i-1,j,k,2,c)*um1 +
> (u(i+1,j,k,5,c)- square(i+1,j,k,c)-
> u(i-1,j,k,5,c)+ square(i-1,j,k,c))*
> c2)
rhs(i,j,k,3,c) = rhs(i,j,k,3,c) + dx3tx1 *
> (u(i+1,j,k,3,c) - 2.0d0*u(i,j,k,3,c) +
> u(i-1,j,k,3,c)) +
> xxcon2 * (vs(i+1,j,k,c) - 2.0d0*vs(i,j,k,c) +
> vs(i-1,j,k,c)) -
> tx2 * (u(i+1,j,k,3,c)*up1 -
> u(i-1,j,k,3,c)*um1)
rhs(i,j,k,4,c) = rhs(i,j,k,4,c) + dx4tx1 *
> (u(i+1,j,k,4,c) - 2.0d0*u(i,j,k,4,c) +
> u(i-1,j,k,4,c)) +
> xxcon2 * (ws(i+1,j,k,c) - 2.0d0*ws(i,j,k,c) +
> ws(i-1,j,k,c)) -
> tx2 * (u(i+1,j,k,4,c)*up1 -
> u(i-1,j,k,4,c)*um1)
rhs(i,j,k,5,c) = rhs(i,j,k,5,c) + dx5tx1 *
> (u(i+1,j,k,5,c) - 2.0d0*u(i,j,k,5,c) +
> u(i-1,j,k,5,c)) +
> xxcon3 * (qs(i+1,j,k,c) - 2.0d0*qs(i,j,k,c) +
> qs(i-1,j,k,c)) +
> xxcon4 * (up1*up1 - 2.0d0*uijk*uijk +
> um1*um1) +
> xxcon5 * (u(i+1,j,k,5,c)*rho_i(i+1,j,k,c) -
> 2.0d0*u(i,j,k,5,c)*rho_i(i,j,k,c) +
> u(i-1,j,k,5,c)*rho_i(i-1,j,k,c)) -
> tx2 * ( (c1*u(i+1,j,k,5,c) -
> c2*square(i+1,j,k,c))*up1 -
> (c1*u(i-1,j,k,5,c) -
> c2*square(i-1,j,k,c))*um1 )
end do
end do
end do
c---------------------------------------------------------------------
c add fourth order xi-direction dissipation
c---------------------------------------------------------------------
if (start(1,c) .gt. 0) then
i = 1
do m = 1, 5
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do j = start(2,c), cell_size(2,c)-end(2,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c)- dssp *
> ( 5.0d0*u(i,j,k,m,c) - 4.0d0*u(i+1,j,k,m,c) +
> u(i+2,j,k,m,c))
end do
end do
end do
i = 2
do m = 1, 5
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do j = start(2,c), cell_size(2,c)-end(2,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
> (-4.0d0*u(i-1,j,k,m,c) + 6.0d0*u(i,j,k,m,c) -
> 4.0d0*u(i+1,j,k,m,c) + u(i+2,j,k,m,c))
end do
end do
end do
endif
do m = 1, 5
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do j = start(2,c), cell_size(2,c)-end(2,c)-1
do i = 3*start(1,c),cell_size(1,c)-3*end(1,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
> ( u(i-2,j,k,m,c) - 4.0d0*u(i-1,j,k,m,c) +
> 6.0*u(i,j,k,m,c) - 4.0d0*u(i+1,j,k,m,c) +
> u(i+2,j,k,m,c) )
end do
end do
end do
end do
if (end(1,c) .gt. 0) then
i = cell_size(1,c)-3
do m = 1, 5
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do j = start(2,c), cell_size(2,c)-end(2,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
> ( u(i-2,j,k,m,c) - 4.0d0*u(i-1,j,k,m,c) +
> 6.0d0*u(i,j,k,m,c) - 4.0d0*u(i+1,j,k,m,c) )
end do
end do
end do
i = cell_size(1,c)-2
do m = 1, 5
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do j = start(2,c), cell_size(2,c)-end(2,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
> ( u(i-2,j,k,m,c) - 4.d0*u(i-1,j,k,m,c) +
> 5.d0*u(i,j,k,m,c) )
end do
end do
end do
endif
c---------------------------------------------------------------------
c compute eta-direction fluxes
c---------------------------------------------------------------------
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do j = start(2,c), cell_size(2,c)-end(2,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
vijk = vs(i,j,k,c)
vp1 = vs(i,j+1,k,c)
vm1 = vs(i,j-1,k,c)
rhs(i,j,k,1,c) = rhs(i,j,k,1,c) + dy1ty1 *
> (u(i,j+1,k,1,c) - 2.0d0*u(i,j,k,1,c) +
> u(i,j-1,k,1,c)) -
> ty2 * (u(i,j+1,k,3,c) - u(i,j-1,k,3,c))
rhs(i,j,k,2,c) = rhs(i,j,k,2,c) + dy2ty1 *
> (u(i,j+1,k,2,c) - 2.0d0*u(i,j,k,2,c) +
> u(i,j-1,k,2,c)) +
> yycon2 * (us(i,j+1,k,c) - 2.0d0*us(i,j,k,c) +
> us(i,j-1,k,c)) -
> ty2 * (u(i,j+1,k,2,c)*vp1 -
> u(i,j-1,k,2,c)*vm1)
rhs(i,j,k,3,c) = rhs(i,j,k,3,c) + dy3ty1 *
> (u(i,j+1,k,3,c) - 2.0d0*u(i,j,k,3,c) +
> u(i,j-1,k,3,c)) +
> yycon2*con43 * (vp1 - 2.0d0*vijk + vm1) -
> ty2 * (u(i,j+1,k,3,c)*vp1 -
> u(i,j-1,k,3,c)*vm1 +
> (u(i,j+1,k,5,c) - square(i,j+1,k,c) -
> u(i,j-1,k,5,c) + square(i,j-1,k,c))
> *c2)
rhs(i,j,k,4,c) = rhs(i,j,k,4,c) + dy4ty1 *
> (u(i,j+1,k,4,c) - 2.0d0*u(i,j,k,4,c) +
> u(i,j-1,k,4,c)) +
> yycon2 * (ws(i,j+1,k,c) - 2.0d0*ws(i,j,k,c) +
> ws(i,j-1,k,c)) -
> ty2 * (u(i,j+1,k,4,c)*vp1 -
> u(i,j-1,k,4,c)*vm1)
rhs(i,j,k,5,c) = rhs(i,j,k,5,c) + dy5ty1 *
> (u(i,j+1,k,5,c) - 2.0d0*u(i,j,k,5,c) +
> u(i,j-1,k,5,c)) +
> yycon3 * (qs(i,j+1,k,c) - 2.0d0*qs(i,j,k,c) +
> qs(i,j-1,k,c)) +
> yycon4 * (vp1*vp1 - 2.0d0*vijk*vijk +
> vm1*vm1) +
> yycon5 * (u(i,j+1,k,5,c)*rho_i(i,j+1,k,c) -
> 2.0d0*u(i,j,k,5,c)*rho_i(i,j,k,c) +
> u(i,j-1,k,5,c)*rho_i(i,j-1,k,c)) -
> ty2 * ((c1*u(i,j+1,k,5,c) -
> c2*square(i,j+1,k,c)) * vp1 -
> (c1*u(i,j-1,k,5,c) -
> c2*square(i,j-1,k,c)) * vm1)
end do
end do
end do
c---------------------------------------------------------------------
c add fourth order eta-direction dissipation
c---------------------------------------------------------------------
if (start(2,c) .gt. 0) then
j = 1
do m = 1, 5
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c)- dssp *
> ( 5.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j+1,k,m,c) +
> u(i,j+2,k,m,c))
end do
end do
end do
j = 2
do m = 1, 5
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
> (-4.0d0*u(i,j-1,k,m,c) + 6.0d0*u(i,j,k,m,c) -
> 4.0d0*u(i,j+1,k,m,c) + u(i,j+2,k,m,c))
end do
end do
end do
endif
do m = 1, 5
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do j = 3*start(2,c), cell_size(2,c)-3*end(2,c)-1
do i = start(1,c),cell_size(1,c)-end(1,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
> ( u(i,j-2,k,m,c) - 4.0d0*u(i,j-1,k,m,c) +
> 6.0*u(i,j,k,m,c) - 4.0d0*u(i,j+1,k,m,c) +
> u(i,j+2,k,m,c) )
end do
end do
end do
end do
if (end(2,c) .gt. 0) then
j = cell_size(2,c)-3
do m = 1, 5
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
> ( u(i,j-2,k,m,c) - 4.0d0*u(i,j-1,k,m,c) +
> 6.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j+1,k,m,c) )
end do
end do
end do
j = cell_size(2,c)-2
do m = 1, 5
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
> ( u(i,j-2,k,m,c) - 4.d0*u(i,j-1,k,m,c) +
> 5.d0*u(i,j,k,m,c) )
end do
end do
end do
endif
c---------------------------------------------------------------------
c compute zeta-direction fluxes
c---------------------------------------------------------------------
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do j = start(2,c), cell_size(2,c)-end(2,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
wijk = ws(i,j,k,c)
wp1 = ws(i,j,k+1,c)
wm1 = ws(i,j,k-1,c)
rhs(i,j,k,1,c) = rhs(i,j,k,1,c) + dz1tz1 *
> (u(i,j,k+1,1,c) - 2.0d0*u(i,j,k,1,c) +
> u(i,j,k-1,1,c)) -
> tz2 * (u(i,j,k+1,4,c) - u(i,j,k-1,4,c))
rhs(i,j,k,2,c) = rhs(i,j,k,2,c) + dz2tz1 *
> (u(i,j,k+1,2,c) - 2.0d0*u(i,j,k,2,c) +
> u(i,j,k-1,2,c)) +
> zzcon2 * (us(i,j,k+1,c) - 2.0d0*us(i,j,k,c) +
> us(i,j,k-1,c)) -
> tz2 * (u(i,j,k+1,2,c)*wp1 -
> u(i,j,k-1,2,c)*wm1)
rhs(i,j,k,3,c) = rhs(i,j,k,3,c) + dz3tz1 *
> (u(i,j,k+1,3,c) - 2.0d0*u(i,j,k,3,c) +
> u(i,j,k-1,3,c)) +
> zzcon2 * (vs(i,j,k+1,c) - 2.0d0*vs(i,j,k,c) +
> vs(i,j,k-1,c)) -
> tz2 * (u(i,j,k+1,3,c)*wp1 -
> u(i,j,k-1,3,c)*wm1)
rhs(i,j,k,4,c) = rhs(i,j,k,4,c) + dz4tz1 *
> (u(i,j,k+1,4,c) - 2.0d0*u(i,j,k,4,c) +
> u(i,j,k-1,4,c)) +
> zzcon2*con43 * (wp1 - 2.0d0*wijk + wm1) -
> tz2 * (u(i,j,k+1,4,c)*wp1 -
> u(i,j,k-1,4,c)*wm1 +
> (u(i,j,k+1,5,c) - square(i,j,k+1,c) -
> u(i,j,k-1,5,c) + square(i,j,k-1,c))
> *c2)
rhs(i,j,k,5,c) = rhs(i,j,k,5,c) + dz5tz1 *
> (u(i,j,k+1,5,c) - 2.0d0*u(i,j,k,5,c) +
> u(i,j,k-1,5,c)) +
> zzcon3 * (qs(i,j,k+1,c) - 2.0d0*qs(i,j,k,c) +
> qs(i,j,k-1,c)) +
> zzcon4 * (wp1*wp1 - 2.0d0*wijk*wijk +
> wm1*wm1) +
> zzcon5 * (u(i,j,k+1,5,c)*rho_i(i,j,k+1,c) -
> 2.0d0*u(i,j,k,5,c)*rho_i(i,j,k,c) +
> u(i,j,k-1,5,c)*rho_i(i,j,k-1,c)) -
> tz2 * ( (c1*u(i,j,k+1,5,c) -
> c2*square(i,j,k+1,c))*wp1 -
> (c1*u(i,j,k-1,5,c) -
> c2*square(i,j,k-1,c))*wm1)
end do
end do
end do
c---------------------------------------------------------------------
c add fourth order zeta-direction dissipation
c---------------------------------------------------------------------
if (start(3,c) .gt. 0) then
k = 1
do m = 1, 5
do j = start(2,c), cell_size(2,c)-end(2,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c)- dssp *
> ( 5.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j,k+1,m,c) +
> u(i,j,k+2,m,c))
end do
end do
end do
k = 2
do m = 1, 5
do j = start(2,c), cell_size(2,c)-end(2,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
> (-4.0d0*u(i,j,k-1,m,c) + 6.0d0*u(i,j,k,m,c) -
> 4.0d0*u(i,j,k+1,m,c) + u(i,j,k+2,m,c))
end do
end do
end do
endif
do m = 1, 5
do k = 3*start(3,c), cell_size(3,c)-3*end(3,c)-1
do j = start(2,c), cell_size(2,c)-end(2,c)-1
do i = start(1,c),cell_size(1,c)-end(1,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
> ( u(i,j,k-2,m,c) - 4.0d0*u(i,j,k-1,m,c) +
> 6.0*u(i,j,k,m,c) - 4.0d0*u(i,j,k+1,m,c) +
> u(i,j,k+2,m,c) )
end do
end do
end do
end do
if (end(3,c) .gt. 0) then
k = cell_size(3,c)-3
do m = 1, 5
do j = start(2,c), cell_size(2,c)-end(2,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
> ( u(i,j,k-2,m,c) - 4.0d0*u(i,j,k-1,m,c) +
> 6.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j,k+1,m,c) )
end do
end do
end do
k = cell_size(3,c)-2
do m = 1, 5
do j = start(2,c), cell_size(2,c)-end(2,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
> ( u(i,j,k-2,m,c) - 4.d0*u(i,j,k-1,m,c) +
> 5.d0*u(i,j,k,m,c) )
end do
end do
end do
endif
do m = 1, 5
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do j = start(2,c), cell_size(2,c)-end(2,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
rhs(i,j,k,m,c) = rhs(i,j,k,m,c) * dt
end do
end do
end do
end do
end do
if (timeron) call timer_stop(t_rhs)
return
end