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