subroutine cleslog_error(type, msg, ival)
implicit none
include 'cles.i'
integer :: type, ival
character(len=*) :: msg
character(len=4) :: val
if ( type .eq. CLESLOG_ERROR_ALLOCATE ) then
call cleslog_trace_enter('Error: Memory allocation failed'//char(0))
else if ( type .eq. CLESLOG_ERROR_DEALLOCATE ) then
call cleslog_trace_enter('Error: Memory deallocation failed'//char(0))
endif
write(val, '(I4)') ival
call cleslog_trace_enter(msg//'|Code='//val//char(0))
stop
end subroutine cleslog_error
subroutine cleslog_log_enter(name)
implicit none
character(len=*) :: name
call cleslog_trace_enter(name//char(0))
end subroutine cleslog_log_enter
subroutine cleslog_log_exit(name)
implicit none
character(len=*) :: name
call cleslog_trace_exit(name//char(0))
end subroutine cleslog_log_exit
subroutine cles_test_pack_dir(ipos, jpos, kpos)
implicit none
integer :: ipos, jpos, kpos
double precision :: x, y, z
character(len=72) :: slog
if ( ipos .ne. 0 ) call cles_xLocation(ipos, x)
if ( jpos .ne. 0 ) call cles_yLocation(jpos, y)
if ( kpos .ne. 0 ) call cles_zLocation(kpos, z)
if ( kpos .eq. 0 ) then
if ( jpos .eq. 0 ) then
write(slog, 100) 'Coordinates (x=',x,')'
else
write(slog, 200) 'Coordinates (x=',x,' y=',y,')'
endif
else
write(slog, 300) 'Coordinates (x=',x,' y=',y,' z=',z,')'
endif
call cleslog_log(slog//char(0))
100 FORMAT(A15,E12.4,A1)
200 FORMAT(A15,E12.4,A3,E12.4,A1)
300 FORMAT(A15,E12.4,A3,E12.4,A3,E12.4,A1)
end subroutine cles_test_pack_dir
subroutine cles_test_getflux(ux, vx, flux)
use method_parms
use Generic_FluxF
implicit none
include 'cles.i'
double precision :: ux(ncomps), flux(nvars)
double precision :: vx(nvars)
integer :: i
vx(1) = ux(1)
do i=2,nvars
vx(i) = ux(i)/ux(1)
enddo
call cles_eqstate(ux, ncomps, vx, nvars, 1, CLES_FALSE)
call GetFluxF(ux, vx, flux)
return
end subroutine cles_test_getflux
subroutine cles_test_eigen(uxref, Rinv, Rdir, lamda)
use method_parms
use Generic_EigenSystem
use Generic_EvalGamma
use cles_interfaces
implicit none
! construct A
double precision :: uxref(ncomps), vxref(nvars), fluxref(nvars)
double precision :: ux(ncomps), vx(nvars), flux(nvars)
double precision :: A(nvars, nvars), Rinv(nvars, nvars)
double precision :: Rdir(nvars, nvars), lamda(nvars)
double precision :: tmp, du, dflux, dnrm, gamma, err
integer :: i, j, is
character(len=nvars*12+2) :: slog
call cles_test_getflux(uxref, vxref, fluxref)
call GetGamma(uxref, vxref, gamma)
write(slog,*) 'Evaluating eigensystem for variables'
call cleslog_log(slog//char(0))
write(slog,'(6(E12.4))') (uxref(i), i=1,nvars)
call cleslog_log(slog//char(0))
write(slog,'(6(E12.4))') (vxref(i), i=1,nvars)
call cleslog_log(slog//char(0))
write(slog,'(A6,E12.4)') ' Gamma', gamma
call cleslog_log(slog//char(0))
write(slog,*) 'Matrix A'
call cleslog_log(slog//char(0))
do j=1, nvars
ux = uxref
du = ux(j)*1.0d-3
if ( abs(du) .lt. 1.0d-4 ) du = 1.0d-4
if ( j .gt. 5 ) then ! scalars
if ( ux(j) + du .gt. ux(1) ) du = -du
endif
ux(j) = ux(j) + du
call cles_test_getflux(ux, vx, flux)
do i=1, nvars
dflux = flux(i)-fluxref(i)
A(i,j) = dflux/du
enddo
enddo
do i=1,nvars
write (slog,'(6(E12.4))') (A(i,j), j=1,nvars)
call cleslog_log(slog//char(0))
enddo
! This bit of code is useful when trying to figure out the scaling
! of the eigenvectors. Since this is arbitrary, in some conditions
! the wrong scale for the eigenvectors leads to accumulation of
! numerical errors because the eigenvectors are not perfectly orthogonal
! numerically. This happens easily with reactive mixtures because the
! enthalpy of formation of species is typically very large with respect
! to the flow kinetic energy and it ends up polluting the eigenvectors.
write(slog,*) 'Test '
call cleslog_log(slog//char(0))
do j=1, nvars
do i=1, nvars
tmp = SUM(Rinv(i,:)*Rdir(:,j))
if ( i .eq. j ) tmp = tmp - 1.0d0
if ( ABS(tmp) .gt. 1.0d-10 ) then
write(slog,'(I2,1X,I2,E12.4)') i, j, tmp
call cleslog_log(slog//char(0))
do is=1, nvars
write(slog,'(A2,I2,2(E12.4))') '=>', is, Rinv(i,is), Rdir(is,j)
call cleslog_log(slog//char(0))
enddo
endif
enddo
enddo
write(slog,*) 'Test '
call cleslog_log(slog//char(0))
do j=1, nvars
do i=1, nvars
tmp = SUM(Rdir(i,:)*Rinv(:,j))
if ( i .eq. j ) tmp = tmp - 1.0d0
if ( ABS(tmp) .gt. 1.0d-10 ) then
write(slog,'(I2,1X,I2,E12.4)') i, j, tmp
call cleslog_log(slog//char(0))
do is=1, nvars
write(slog,'(A2,I2,2(E12.4))') '=>', is, Rdir(i,is), Rinv(is,j)
call cleslog_log(slog//char(0))
enddo
endif
enddo
enddo
write(slog,*) 'Test '
call cleslog_log(slog//char(0))
do j=1, nvars
do i=1, nvars
Rdir(i,j) = Rdir(i,j)*lamda(j)
enddo
enddo
do i=1, nvars
do j=1, nvars
tmp = SUM(Rdir(i,:)*Rinv(:,j))
dnrm = abs(A(i,j))
err = abs(A(i,j)-tmp)
if ( err .gt. 1.0d-3*max(dnrm,1.0d-10) ) then
err = err/max(dnrm,abs(tmp))
write(slog,'(I2,1X,I2,3(E12.4))') i, j, A(i,j), tmp, err
call cleslog_log(slog//char(0))
endif
enddo
enddo
return
end subroutine cles_test_eigen
subroutine cles_test_tcdstencil()
use method_parms
use Interp_coeffs
implicit none
include 'cles.i'
integer :: nx, ilo, ihi, ib
double precision, allocatable :: a(:), flux(:), der(:), b(:)
character(len=72) :: slog
write(slog,'(A16,I2)') 'Stencil width : ', stencil
call cleslog_log(slog//char(0))
write(slog,'(A17,I2)') 'Boundary width : ', tcd_bndry_width
call cleslog_log(slog//char(0))
call cleslog_log('<----- Divergence flux test ----->'//char(0))
! test n > 2 x boundary width stencil
nx = stencil + 2 * tcd_bndry_width
ilo = 1-nghost
ihi = nx+nghost
write(slog,'(A16,I2)') 'Wide patch, n : ', nx
call cleslog_log(slog//char(0))
allocate (a(ilo:ihi))
allocate (der(nx))
allocate (flux(nx+1))
do ib=0, tcd_bndry_width
write(slog,*) 'Testing ', ib, '/CORE bndry'
call cleslog_log(slog//char(0))
call cles_test_tcdwide(nx, ilo, ihi, a, der, flux, &
ib, CLES_PATCH_CORE)
write(slog,*) 'Testing CORE/', ib, ' bndry'
call cleslog_log(slog//char(0))
call cles_test_tcdwide(nx, ilo, ihi, a, der, flux, &
CLES_PATCH_CORE, ib)
enddo
deallocate(a)
deallocate(der)
deallocate(flux)
do nx=1, tcd_bndry_width+1
ilo = 1-nghost
ihi = nx+nghost
write(slog,'(A17,I2)') 'Short patch, n : ', nx
call cleslog_log(slog//char(0))
allocate (a(ilo:ihi))
allocate (der(nx))
allocate (flux(nx+1))
do ib=0,tcd_bndry_width
write(slog,'(A14,I2)') 'Boundary mark:', ib
call cleslog_log(slog//char(0))
call cles_test_tcdwide(nx, ilo, ihi, a, der, flux, &
ib, CLES_PATCH_CORE)
call cles_test_tcdwide(nx, ilo, ihi, a, der, flux, &
CLES_PATCH_CORE, ib)
enddo
deallocate(a)
deallocate(der)
deallocate(flux)
enddo
!----------------------------------------------------
call cleslog_log('<----- Skew flux tests ----->'//char(0))
! test n > 2 x boundary width stencil
nx = stencil + 2 * tcd_bndry_width
ilo = 1-nghost
ihi = nx+nghost
write(slog,'(A16,I2)') 'Wide patch, n : ', nx
call cleslog_log(slog//char(0))
allocate (a(ilo:ihi))
allocate (b(ilo:ihi))
allocate (der(nx))
allocate (flux(nx+1))
do ib=0, tcd_bndry_width
write(slog,*) 'Testing ', ib, '/CORE bndry'
call cleslog_log(slog//char(0))
call cles_test_tcdwideskew(nx, ilo, ihi, a, b, der, flux, &
ib, CLES_PATCH_CORE)
write(slog,*) 'Testing CORE/', ib, ' bndry'
call cleslog_log(slog//char(0))
call cles_test_tcdwideskew(nx, ilo, ihi, a, b, der, flux, &
CLES_PATCH_CORE, ib)
enddo
deallocate(a)
deallocate(b)
deallocate(der)
deallocate(flux)
do nx=1, tcd_bndry_width+1
ilo = 1-nghost
ihi = nx+nghost
write(slog,'(A17,I2)') 'Short patch, n : ', nx
call cleslog_log(slog//char(0))
allocate (a(ilo:ihi))
allocate (b(ilo:ihi))
allocate (der(nx))
allocate (flux(nx+1))
do ib=0,tcd_bndry_width
write(slog,'(A14,I2)') 'Boundary mark:', ib
call cleslog_log(slog//char(0))
call cles_test_tcdwideskew(nx, ilo, ihi, a, b, der, flux, &
ib, CLES_PATCH_CORE)
call cles_test_tcdwideskew(nx, ilo, ihi, a, b, der, flux, &
CLES_PATCH_CORE, ib)
enddo
deallocate(a)
deallocate(b)
deallocate(der)
deallocate(flux)
enddo
return
end subroutine cles_test_tcdstencil
subroutine cles_test_tcdwide(nx, ilo, ihi, a, der, flux, ixlo, ixup)
use method_parms
use Generic_FDInterp
implicit none
include 'cles.i'
integer :: nx, ilo, ihi, ixlo, ixup
double precision :: a(ilo:ihi), flux(nx+1), der(nx)
integer :: ord, i
double precision :: center
double precision :: error(nx,order)
character(len=72) :: slog
if ( mod(nx,2) .eq. 0 ) then
center = 0.5d0*(nx+1)
else
center = 0.5d0*nx
endif
do ord=1, order
do i=ilo, ihi
a(i) = (i-center)**ord
enddo
call FDInterpolate(1, a, flux, ilo, ihi, nx, ixlo, ixup)
der(:) = flux(2:nx+1)-flux(1:nx)
if ( ord .eq. 1 ) then
error(:,ord) = der(:)-1.0d0 ! exact result
else
do i=1, nx
error(i,ord) = der(i) - ord*(i-center)**(ord-1)
enddo
endif
enddo
! print results
write(slog,'(A5,I9,9(I12))') 'i/ord', (ord, ord=1, order)
call cleslog_log(slog//char(0))
do i=1, nx
write(slog,'(I2,10(E12.4))') i, (error(i,ord), ord=1,order)
call cleslog_log(slog//char(0))
enddo
return
end subroutine cles_test_tcdwide
subroutine cles_test_tcdwideskew(nx, ilo, ihi, a, b, der, flux, ixlo, ixup)
use method_parms
use Generic_FDSkewInterp
implicit none
include 'cles.i'
integer :: nx, ilo, ihi, ixlo, ixup
double precision :: a(ilo:ihi), b(ilo:ihi), flux(nx+1), der(nx)
integer :: orda, ordb, i
double precision :: center
double precision :: error(nx,order*order), exa(nx), exb(nx)
character(len=72) :: slog
character(len=7) :: id(order*order)
if ( mod(nx,2) .eq. 0 ) then
center = 0.5d0*(nx+1)
else
center = 0.5d0*nx
endif
do ordb=1, order
do i=ilo, ihi
b(i) = (i-center)**ordb
enddo
if ( ordb .eq. 1 ) then
exb(:) = 1.0d0 ! exact result
else
do i=1, nx
exb(i) = ordb*(i-center)**(ordb-1)
enddo
endif
do orda=1, order
do i=ilo, ihi
a(i) = (i-center)**orda
enddo
call FDSkewInterpolate(a, b, flux, ilo, ihi, nx, ixlo, ixup)
der(:) = flux(2:nx+1)-flux(1:nx)
if ( orda .eq. 1 ) then
exa(:) = 1.0d0 ! exact result
else
do i=1, nx
exa(i) = orda*(i-center)**(orda-1)
enddo
endif
error(:,orda+(ordb-1)*order) = der(:)-a(1:nx)*exb(:)-b(1:nx)*exa(:)
enddo
enddo
! print results
do ordb=1, order
do orda=1, order
write(id(orda+(ordb-1)*order),'(A1,I2,A1,I2,A1)') &
'(', orda, ',', ordb, ')'
enddo
enddo
write(slog,'(A5,2X,A7,9(5X,A7))') 'i/ord', (id(orda), orda=1, order*order)
call cleslog_log(slog//char(0))
do i=1, nx
write(slog,'(I2,10(E12.4))') i, (error(i,orda), orda=1,order*order)
call cleslog_log(slog//char(0))
enddo
return
end subroutine cles_test_tcdwideskew