MODULE Generic_GetRHS
! ---- Computes the RHS as the derivative of the flux vector
! ---- as the derivate of the x (1) and y (2) and z (3) directions and
! ---- saves them seperately in rhs(nvars,direction,:,:,:)
!
! ---- The locations where WENO is used are stored in the feild
! ---- 'dcflag'
INTERFACE GetRHS
MODULE PROCEDURE OneDGetRHS
MODULE PROCEDURE TwoDGetRHS
MODULE PROCEDURE ThreeDGetRHS
END INTERFACE
CONTAINS
SUBROUTINE OneDGetRHS(rk, ux,vx,rhs,fxi,dcflag,ifilter)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
use cles_interfaces
! ---- Shared Procedures
USE Generic_GetFlux
USE Generic_AcousticBC
USE Generic_CellWallFlux
USE Generic_Transport
USE Generic_ViscousAtCellWalls
USE Generic_FDFluxInterps
USE Generic_FDInterpFlux
USE Generic_Source
IMPLICIT NONE
include 'cles.i'
INTEGER, INTENT(IN) :: rk, ifilter
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi)
DOUBLE PRECISION, INTENT(OUT) :: vx(nvars,ixlo:ixhi)
DOUBLE PRECISION, INTENT(OUT) :: rhs(nvars,1:nx)
DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ixlo:ixhi,1)
INTEGER, INTENT(INOUT) :: dcflag(1:nx+1,1)
DOUBLE PRECISION, ALLOCATABLE :: fx(:,:)
INTEGER :: direction
INTEGER :: iret
INTEGER :: i
ALLOCATE ( fx(nvars,ixlo:ixhi) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'OneDGetRHS/fx', iret)
endif
! ---- initialize flux interpolates and mask
fxi = 0.0D0
rhs = 0.0D0
! ---- Convert conserved variables to Primative
vx(1,:) = ux(1,:)
do i=2,nvars
vx(i,:) = ux(i,:)/ux(1,:)
enddo
call cles_eqstate(ux, ncomps, vx, nvars, mx, CLES_FALSE)
! ---- compute the transport terms
IF (useViscous.eq.CLES_TRUE) CALL GetTransport(ux, vx)
! ---- Initialize dcflag
if ( rk .eq. 1 ) then
DO direction = 1,1,1
call cles_dcflag(ux,vx,dcflag,ncomps,nvars, &
ixlo,ixhi, nx, dx, direction, enoOrder)
enddo
endif
! ---- Add body force source term
IF(useSource.eq.CLES_TRUE) CALL AddSource(ux,vx,rhs)
! ---- compute the flux interpolates in
! ---- the x (1) direction
DO direction = 1,1,1
! ---- compute the flux vector for this direction
CALL GetFlux(ux,vx,fx,direction)
! ---- compute the interpolated fluxes on the cell walls
CALL CellWallFlux(ux,vx,fx,fxi,dcflag,direction,ifilter)
enddo
! ---- apply the characteristic boundary conditions
CALL AcousticBC(fxi, ux, vx, rhs, dcflag)
! the RHS must be computed outside of the previouse loop to
! give a chance to the characteristic boundary conditions to
! perform the decomposition and cancel incoming contributions
! from possible source terms
do direction = 1, 1, 1
IF (useViscous.eq.CLES_TRUE) THEN
! ---- Modify the flux at the cell wall to include
! ---- the viscous terms.
CALL ViscousFluxModification(fxi, ux, vx, direction)
END IF
! ---- Compute the flux derivatives from the interpolates
! ---- note this is only done for the interior (not Ghost)
! ---- cells.
CALL AccumulateFluxDiffs(rhs, fxi, direction)
! ---- finish the 'direction' loop
END DO
DEALLOCATE(fx, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'OneDGetRHS/fx', iret)
endif
END SUBROUTINE OneDGetRHS
SUBROUTINE TwoDGetRHS(rk, ux,vx,rhs,fxi,dcflag, ifilter)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
use cles_interfaces
! ---- Shared Procedures
USE Generic_GetFlux
USE Generic_AcousticBC
USE Generic_CellWallFlux
USE Generic_Transport
USE Generic_ViscousAtCellWalls
USE Generic_FDFluxInterps
USE Generic_FDInterpFlux
USE Generic_Source
IMPLICIT NONE
include 'cles.i'
INTEGER, INTENT(IN) :: rk, ifilter
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi)
DOUBLE PRECISION, INTENT(OUT) :: vx(nvars,ixlo:ixhi,iylo:iyhi)
DOUBLE PRECISION, INTENT(OUT) :: rhs(nvars,1:nx,1:ny)
DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,2)
INTEGER, INTENT(INOUT) :: dcflag(1:nx+1,1:ny+1,2)
DOUBLE PRECISION, ALLOCATABLE :: fx(:,:,:)
INTEGER :: iret
INTEGER :: direction
INTEGER :: i,j, d, n
ALLOCATE ( fx(nvars,ixlo:ixhi,iylo:iyhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'TwoDGetRHS/fx', iret)
endif
! ---- initialize flux interpolates and mask
fxi = 0.0D0
rhs = 0.0D0
! ---- Convert conserved variables to Primative
vx(1,:,:) = ux(1,:,:)
do i=2,nvars
vx(i,:,:) = ux(i,:,:)/ux(1,:,:)
enddo
n = mx*my
call cles_eqstate(ux, ncomps, vx, nvars, n, CLES_FALSE)
! --- compute the transport terms
IF (useViscous.eq.CLES_TRUE) CALL GetTransport(ux, vx)
! ---- Initialize dcflag
if ( rk .eq. 1 ) then
DO direction = 1,2,1
call cles_dcflag(ux,vx,dcflag,ncomps,nvars, &
ixlo,ixhi,iylo,iyhi, nx, ny, dx, dy, &
direction, enoOrder)
enddo
endif
! ---- Add body force source term
IF(useSource.eq.CLES_TRUE) CALL AddSource(ux,vx,rhs)
! ---- compute the flux interpolates in
! ---- the x (1) direction, and then in the y (2) direction
DO direction = 1,2,1
! ---- compute the inviscid flux vector for this direction
CALL GetFlux(ux,vx,fx,direction)
! ---- compute the interpolated fluxes on the cell walls
CALL CellWallFlux(ux,vx,fx,fxi,dcflag,direction,ifilter)
enddo
! ---- apply the characteristic boundary conditions
CALL AcousticBC(fxi, ux, vx, rhs, dcflag)
! the RHS must be computed outside of the previouse loop to
! give a chance to the characteristic boundary conditions to
! perform the decomposition and cancel incoming contributions
! from possible source terms
do direction = 1, 2, 1
IF (useViscous.eq.CLES_TRUE) THEN
! ---- Modify the flux at the cell wall to include
! ---- the viscous terms.
CALL ViscousFluxModification(fxi, ux, vx, direction)
END IF
! ---- Compute the flux derivatives from the interpolates
! ---- note this is only done for the interior (not Ghost)
! ---- cells.
CALL AccumulateFluxDiffs(rhs, fxi, direction)
END DO
DEALLOCATE ( fx , stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'TwoDGetRHS/fx', iret)
endif
END SUBROUTINE TwoDGetRHS
SUBROUTINE ThreeDGetRHS(rk, ux,vx,rhs,fxi,dcflag, ifilter)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
use cles_interfaces
! ---- Shared Procedures
USE Generic_GetFlux
USE Generic_AcousticBC
USE Generic_CellWallFlux
USE Generic_Transport
USE Generic_ViscousAtCellWalls
USE Generic_FDFluxInterps
USE Generic_FDInterpFlux
USE Generic_Source
IMPLICIT NONE
include 'cles.i'
INTEGER, INTENT(IN) :: rk, ifilter
DOUBLE PRECISION, INTENT(INOUT) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION, INTENT(OUT) :: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION, INTENT(OUT) :: rhs(nvars,1:nx,1:ny,1:nz)
DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi,3)
INTEGER, INTENT(INOUT) :: dcflag(1:nx+1,1:ny+1,1:nz+1,3)
DOUBLE PRECISION, ALLOCATABLE :: fx(:,:,:,:)
INTEGER :: iret
INTEGER :: direction
INTEGER :: i,j,k, d, n
call cleslog_log_enter('ThreeDGetRHS')
ALLOCATE ( fx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDGetRHS/fx', iret)
endif
! ---- initialize flux interpolates and mask
fxi(:,:,:,:,:) = 0.0D0
rhs(:,:,:,:) = 0.0D0
! ---- Convert conserved variables to primitive
vx(1,:,:,:) = ux(1,:,:,:)
do i=2,nvars
vx(i,:,:,:) = ux(i,:,:,:)/ux(1,:,:,:)
enddo
! ---- set up viscous terms if needed (only temperature is required)
! ---- we assume that the viscosity does not depend explicitly on the
! ---- pressure. Because it is not defined at this point. We determine
! ---- it after sgs terms are computed so the actual temperature and
! ---- pressure contain the effect of the sgs kinetic energy
IF (useViscous.eq.CLES_TRUE) CALL GetTransport(ux, vx)
! ---- Computes the model terms, and makes the
! ---- SGS fix to the pressure
! ---- also stores the scalar variance in
! ---- the extended vector of state
if (useLES.eq.CLES_TRUE) then
call InitializeLES(ux,vx)
! need to fix the values of sgskin on the ghost cells.
! this is required for consistency since the users have
! no control over what will live there. Only the LES module
! is capable of determining the state here.
call FixSgsBoundaries(ux, vx)
endif
! ---- now the primatives have the sgs correction
n = mx*my*mz
! compute the pressure and temperature from updated sgske
call cles_eqstate(ux, ncomps, vx, nvars, n, useLES)
! ---- Initialize dcflag
if ( rk .eq. 1 ) then
DO direction = 1,3,1
call cles_dcflag(ux,vx,dcflag,ncomps,nvars, &
ixlo,ixhi,iylo,iyhi,izlo,izhi,nx,ny,nz, &
dx,dy,dz,direction,enoOrder)
enddo
endif
! ---- Add body force source term
IF(useSource.eq.CLES_TRUE) CALL AddSource(ux,vx,rhs)
! ---- compute the flux interpolates in
! ---- the x (1) direction, the y (2) direction, the z (3) direction
DO direction = 1,3,1
! ---- compute the inviscid flux vector for this direction
CALL GetFlux(ux,vx,fx,direction)
! ---- compute the interpolated fluxes on the cell walls
CALL CellWallFlux(ux,vx,fx,fxi,dcflag,direction,ifilter)
enddo
! ---- apply the characteristic boundary conditions
CALL AcousticBC(fxi, ux, vx, rhs, dcflag)
! the RHS must be computed outside of the previouse loop to
! give a chance to the characteristic boundary conditions to
! perform the decomposition and cancel incoming contributions
! from possible source terms
DO direction = 1,3,1
if (useViscous.eq.CLES_TRUE) THEN
! ---- Modify the flux at the cell wall to include
! ---- the viscous terms.
CALL ViscousFluxModification(fxi, ux, vx, direction)
endif
! ---- modify the flux at the cell walls to account for viscousity
! ---- the subgrid controbution to the fluxes
IF (useLES.eq.CLES_TRUE) CALL AddSgsFlux(fxi,ux,vx,direction)
! ---- Compute the flux derivatives from the interpolates
! ---- note this is only done for the interior (not Ghost)
! ---- cells.
CALL AccumulateFluxDiffs(rhs, fxi, direction)
END DO
DEALLOCATE ( fx ,stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDGetRHS/fx', iret)
endif
call cleslog_log_exit('ThreeDGetRHS')
END SUBROUTINE ThreeDGetRHS
SUBROUTINE FixSgsBoundaries(ux, vx)
USE mesh
USE array_bounds
USE method_parms
implicit none
include 'cles.i'
DOUBLE PRECISION, INTENT(INOUT) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION, INTENT(INOUT) :: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
INTEGER :: i,j,k
if ( useLES .eq. CLES_FALSE ) return
if ( bc_ixlow .ne. CLES_PATCH_CORE ) then
do k=1,nz
do j=1,ny
do i=ixlo+1,1-bc_ixlow
! retreive original pressure
call cles_eqstate(ux(:,i,j,k), ncomps, vx(:,i,j,k), &
nvars, 1, CLES_FALSE)
! add sgs term to total energy
call cles_inveqst(ux(:,i,j,k),ncomps,vx(:,i,j,k), &
nvars,1,useLES)
enddo
enddo
enddo
endif
if ( bc_ixup .ne. CLES_PATCH_CORE ) then
do k=1,nz
do j=1,ny
do i=nx+bc_ixup,ixhi-1
! retreive original pressure
call cles_eqstate(ux(:,i,j,k), ncomps, vx(:,i,j,k), &
nvars, 1, CLES_FALSE)
! add sgs term to total energy
call cles_inveqst(ux(:,i,j,k),ncomps,vx(:,i,j,k), &
nvars,1,useLES)
enddo
enddo
enddo
endif
if ( bc_iylow .ne. CLES_PATCH_CORE ) then
do k=1,nz
do j=iylo+1,1-bc_iylow
do i=1,nx
! retreive original pressure
call cles_eqstate(ux(:,i,j,k), ncomps, vx(:,i,j,k), &
nvars, 1, CLES_FALSE)
! add sgs term to total energy
call cles_inveqst(ux(:,i,j,k),ncomps,vx(:,i,j,k), &
nvars,1,useLES)
enddo
enddo
enddo
endif
if ( bc_iyup .ne. CLES_PATCH_CORE ) then
do k=1,nz
do j=ny+bc_iyup,iyhi-1
do i=1,nx
! retreive original pressure
call cles_eqstate(ux(:,i,j,k), ncomps, vx(:,i,j,k), &
nvars, 1, CLES_FALSE)
! add sgs term to total energy
call cles_inveqst(ux(:,i,j,k),ncomps,vx(:,i,j,k), &
nvars,1,useLES)
enddo
enddo
enddo
endif
if ( bc_izlow .ne. CLES_PATCH_CORE ) then
do k=izlo+1,1-bc_izlow
do j=1,ny
do i=1,nx
! retreive original pressure
call cles_eqstate(ux(:,i,j,k), ncomps, vx(:,i,j,k), &
nvars, 1, CLES_FALSE)
! add sgs term to total energy
call cles_inveqst(ux(:,i,j,k),ncomps,vx(:,i,j,k), &
nvars,1,useLES)
enddo
enddo
enddo
endif
if ( bc_izup .ne. CLES_PATCH_CORE ) then
do k=nz+bc_izup,izhi-1
do j=1,ny
do i=1,nx
! retreive original pressure
call cles_eqstate(ux(:,i,j,k), ncomps, vx(:,i,j,k), &
nvars, 1, CLES_FALSE)
! add sgs term to total energy
call cles_inveqst(ux(:,i,j,k),ncomps,vx(:,i,j,k), &
nvars,1,useLES)
enddo
enddo
enddo
endif
END SUBROUTINE FixSgsBoundaries
END MODULE Generic_GetRHS