!
! Correct fluxes at boundaries for subsonic flows to avoid
! ill posedness and allow large temporal integrations.
! Creator: Carlos Pantano
! Copyright: California Institute of Technology (2004-2006)
! Notes: This module is part of the CLES patch solver.
MODULE Generic_AcousticBC
SAVE
INTEGER CLESLOG_CBC
INTEGER CLESLOG_CBC_EIGENSYSTEM
INTEGER CLESLOG_CBC_REFPRESSURE
PARAMETER (CLESLOG_CBC=1)
PARAMETER (CLESLOG_CBC_EIGENSYSTEM=1)
PARAMETER (CLESLOG_CBC_REFPRESSURE=2)
! Determine if acoustic BC should be use
! Use characteristic acoustic BC and incomming as suggested by
! Thompson; Rudy & Strikwerda; Poinsot & Lele; Stanley & Sarkar;
! Giles; Rowley & Colonius. (many theoretical contributions) but
! using the SAT procesure of Carpenter, Gotlieb and Abarbanel.
INTEGER :: cbc_active = 0
INTEGER :: cbc_mode = 2 ! default is 2d mode
DOUBLE PRECISION, ALLOCATABLE :: cbc_force_matrix(:,:)
DOUBLE PRECISION :: cbc_sigma = 0.25d0
DOUBLE PRECISION :: cbc_Mach = 0.0d0
DOUBLE PRECISION :: cbc_factor = 0.125d0
! ---- Adds acoustic forcing to the state flux to mimmic the
! ---- physical effects that a real infinite domain would
! ---- have on a subsonic outflow.
INTERFACE AcousticBC
MODULE PROCEDURE OneDAcousticBC
MODULE PROCEDURE TwoDAcousticBC
MODULE PROCEDURE ThreeDAcousticBC
END INTERFACE
INTERFACE AcousticViscousBC
MODULE PROCEDURE OneDAcousticViscousBC
MODULE PROCEDURE TwoDAcousticViscousBC
MODULE PROCEDURE ThreeDAcousticViscousBC
END INTERFACE
CONTAINS
SUBROUTINE SingleCSAT(ux, vx, uxref, source, sat, flux, p, dl, tau, &
cbc_type, ipos, jpos, kpos)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
use cles_interfaces
USE Generic_Source
USE Generic_EigenSystem
USE Generic_OuterBoundary
USE Interp_coeffs
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps), dl, p
DOUBLE PRECISION, INTENT(IN) :: vx(nvars), tau(nvars)
DOUBLE PRECISION, INTENT(IN) :: uxref(ncomps)
DOUBLE PRECISION, INTENT(INOUT) :: flux(nvars,3)
DOUBLE PRECISION, INTENT(INOUT) :: source(nvars)
DOUBLE PRECISION, INTENT(OUT) :: sat(nvars)
INTEGER, INTENT(IN) :: ipos, jpos, kpos, cbc_type
DOUBLE PRECISION :: a, c, fx, fy, fz
DOUBLE PRECISION :: lamda(nvars), Rdir(nvars,nvars), Rinv(nvars,nvars)
DOUBLE PRECISION :: vec(nvars), dU(nvars), DeltaG(nvars)
INTEGER :: i, k
CHARACTER*56 slog
call cleslog_log_enter('SingleCSAT')
Call EigenSystem(ux, ux, vx, vx, Rinv, Rdir, lamda)
if ( cleslog_active(CLESLOG_CBC, CLESLOG_CBC_EIGENSYSTEM) &
.eq. CLES_TRUE ) then
call cles_test_pack_dir(ipos, jpos, kpos)
call cles_test_eigen(ux, Rinv, Rdir, lamda)
call cleslog_log_flush()
endif
c = (lamda(nvars)-lamda(1))*0.5d0 ! for example
DeltaG(:) = 0.0
do k=1,3
do i=1,nvars
vec(i) = SUM(Rinv(i,:)*flux(:,k))
enddo
flux(:,k) = vec(:)
enddo
if ( useSource .eq. CLES_TRUE ) then
! do i=1, nvars
! dU(i) = SUM(Rinv(i,:)*source(:))
! enddo
! source = dU
endif
! save characteristic incomming strength
fx = flux(1,1)
fy = flux(1,2)
fz = flux(1,3)
flux(:,:) = 0.0d0
! Point #1
! The only values that should be non-zero are those values that
! correspond to characteristic waves trying to get inside the domain
! for which we must know their value.
if ( lamda(nvars) .lt. 0.0 ) then
DeltaG(:) = lamda(:)
! if ( useSource .eq. CLES_TRUE ) source(:) = 0.0d0
else if ( lamda(2) .lt. 0.0 ) then
DeltaG(1:nvars-1) = lamda(1:nvars-1)
! if ( useSource .eq. CLES_TRUE ) source(1:nvars-1) = 0.0d0
else if ( lamda(1) .lt. 0.0 ) then
DeltaG(1) = lamda(1)
! if ( useSource .eq. CLES_TRUE ) source(1) = 0.0d0
else
! Nothing to do
endif
a = tcd_bndry(1,1)/dl
! Add penalty terms for those directions that are incomming
dU(1:nvars) = ux(1:nvars)-uxref(1:nvars)
do i=1, nvars
vec(i) = DeltaG(i)*a*tau(i)*SUM(Rinv(i,:)*dU(:))
enddo
if ( cbc_type .eq. CLES_CBC_OUTFLOW .and. lamda(1) .lt. 0.0d0 &
.and. cbc_Mach .lt. 1.0d0 ) then
if ( cleslog_active(CLESLOG_CBC, CLESLOG_CBC_REFPRESSURE) &
.eq. CLES_TRUE ) then
call cles_test_pack_dir(ipos, jpos, kpos)
write(slog, '(A10,E12.4)') 'Pressure =', p
call cleslog_log(slog//char(0))
call cleslog_log_flush()
endif
vec(1) = 0.0d0
! if ( useSource .eq. CLES_TRUE ) source(1) = 0.0d0
! cancel outgoing differential part
flux(1,1) = cbc_factor*(vx(5)-p)/((xrg-xlg)*c*Rdir(1,1))-fx
if ( cbc_mode .eq. CLES_CBC_MODE2D ) then
flux(1,2) = - fy
flux(1,3) = - fz
endif
endif
if ( useSource .eq. CLES_TRUE ) then
! do i=1, nvars
! dU(i) = SUM(Rdir(i,:)*source(:))
! enddo
! source = dU
endif
! Multiply by R (make conservative)
do i=1, nvars
sat(i) = SUM(Rdir(i,:)*vec(:))
enddo
! rotate back flux part
do k=1,3
do i=1,nvars
vec(i) = SUM(Rdir(i,:)*flux(:,k))
enddo
flux(:,k) = vec(:)
enddo
call cleslog_log_exit('SingleCSAT')
END SUBROUTINE SingleCSAT
subroutine DirectRotateXYZ(ux, vx, uxref, src, flux, P, bc)
use method_parms
use Generic_Source
implicit none
include 'cles.i'
double precision, intent(inout) :: ux(ncomps), vx(nvars)
double precision, intent(inout) :: uxref(ncomps), src(nvars)
double precision, intent(inout) :: flux(nvars,3)
double precision, intent(in) :: P(3,3)
integer, intent(in) :: bc
integer i, k
double precision vec(3)
! vector rotation
do i=1,3
vec(i) = SUM(ux(2:4)*P(:,i))
enddo
ux(2:4) = vec(:)
do i=1,3
vec(i) = SUM(vx(2:4)*P(:,i))
enddo
vx(2:4) = vec(:)
do i=1,3
vec(i) = SUM(uxref(2:4)*P(:,i))
enddo
uxref(2:4) = vec(:)
if ( useSource .eq. CLES_TRUE ) then
do i=1,3
vec(i) = SUM(src(2:4)*P(:,i))
enddo
src(2:4) = vec(:)
endif
! change of independent variables
do i=1,nvars
do k=1,3
vec(k) = SUM(P(:,k)*flux(i,:))
enddo
flux(i,:) = vec(:)
enddo
! rotation
do k=1, 3
do i=1,3
vec(i) = SUM(flux(2:4,k)*P(:,i))
enddo
flux(2:4,k) = vec(:)
enddo
end subroutine DirectRotateXYZ
subroutine InverseRotateXYZ(src, sat, flux, P, bc)
use method_parms
use Generic_Source
implicit none
include 'cles.i'
double precision, intent(inout) :: src(nvars), sat(nvars)
double precision, intent(inout) :: flux(nvars,3)
double precision, intent(in) :: P(3,3)
integer :: bc
double precision vec(3)
integer i, k
do i=1,3
vec(i) = SUM(P(i,:)*sat(2:4))
enddo
sat(2:4) = vec(:)
if ( useSource .eq. CLES_TRUE ) then
do i=1,3
vec(i) = SUM(P(i,:)*src(2:4))
enddo
src(2:4) = vec(:)
endif
! rotation first
do k=1, 3
do i=1,3
vec(i) = SUM(P(i,:)*flux(2:4,k))
enddo
flux(2:4,k) = vec(:)
enddo
! change of independent variables
do i=1,nvars
do k=1,3
vec(k) = SUM(flux(i,:)*P(k,:))
enddo
flux(i,:) = vec(:)
enddo
end subroutine InverseRotateXYZ
subroutine setproj(theta, phi, proj)
implicit none
double precision theta, phi, proj(3,3)
double precision ctheta, stheta, cphi, sphi
! projection matrix
ctheta = cos(theta)
stheta = sin(theta)
cphi = cos(phi)
sphi = sin(phi)
proj(1,1) = ctheta*cphi
proj(1,2) = -stheta
proj(1,3) = -sphi*ctheta
proj(2,1) = cphi*stheta
proj(2,2) = ctheta
proj(2,3) = -sphi*stheta
proj(3,1) = sphi
proj(3,2) = 0.0d0
proj(3,3) = cphi
return
end subroutine setproj
SUBROUTINE OneDAcousticBC(fxi, ux, vx, source, dcflag)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
USE Generic_Source
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi)
DOUBLE PRECISION,INTENT(INOUT) :: fxi(ncomps,ixlo:ixhi,1)
! the cell centered flux
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi)
DOUBLE PRECISION, INTENT(INOUT) :: source(nvars,1:nx)
INTEGER, INTENT(IN) :: dcflag(1:nx+1,1)
INTEGER :: ip
DOUBLE PRECISION :: uxloc(ncomps), vxloc(nvars), uxref(ncomps)
DOUBLE PRECISION :: sat(nvars), flux(nvars,3), Proj(3,3)
if ( cbc_active .ne. CLES_TRUE ) return
!-----------------------------------------------------------
!
! LEFT BOUNDARY
!
!-----------------------------------------------------------
if ( cbc_ixlow .ne. CLES_CBC_NONE ) then
Proj = 0.0d0
Proj(1,1) = -1.0d0
Proj(2,2) = -1.0d0
Proj(3,3) = 1.0d0
flux(:,1) = (fxi(1:nvars,2,1)-fxi(1:nvars,1,1))/dx
flux(:,2:3) = 0.0d0
uxloc = ux(:,1)
vxloc = vx(:,1)
uxref = ux(:,0)
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(1,1), flux, Proj, &
cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(1,1), sat, flux, vx(5,0), &
dx, cbc_force_matrix(:,1), cbc_ixlow, 1, 0, 0)
! rotation back
call InverseRotateXYZ(source(1,1), sat, flux, Proj, cbc_ixlow)
! Correct flux to the last point
source(:,1) = source(:,1) - sat(:) - flux(:,1)
endif
!-----------------------------------------------------------
!
! RIGHT BOUNDARY
!
!-----------------------------------------------------------
if ( cbc_ixup .ne. CLES_CBC_NONE ) then
flux(:,1) = (fxi(1:nvars,nx+1,1)-fxi(1:nvars,nx,1))/dx
flux(:,2:3) = 0.0d0
call SingleCSAT(ux(1,nx), vx(1,nx), ux(1,nx+1), source(1,nx), sat, &
flux, vx(5,nx+1), dx, cbc_force_matrix(:,2), cbc_ixup, nx, 0, 0)
! Correct flux to the last point
source(:,nx) = source(:,nx) - sat(:) - flux(:,1)
endif
RETURN
END SUBROUTINE OneDAcousticBC
SUBROUTINE TwoDAcousticBC(fxi, ux, vx, source, dcflag)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
USE Generic_Source
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi)
DOUBLE PRECISION, INTENT(INOUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,2)
! the cell centered flux
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi)
DOUBLE PRECISION, INTENT(INOUT) :: source(nvars,1:nx,1:ny)
INTEGER, INTENT(IN) :: dcflag(1:nx+1,1:ny+1,2)
INTEGER :: i, j, ip
DOUBLE PRECISION :: uxloc(ncomps), vxloc(nvars), uxref(ncomps)
DOUBLE PRECISION :: Proj(3,3), flux(nvars,3), sat(nvars)
INTEGER :: ilow, iup, jlow, jup
DOUBLE PRECISION :: p, dl, isqr2
if ( cbc_active .ne. CLES_TRUE ) return
isqr2 = 1.0d0/sqrt(2.0d0)
!-----------------------------------------------------------
!
! LEFT BOUNDARY
!
!-----------------------------------------------------------
if ( cbc_ixlow .ne. CLES_CBC_NONE ) then
jlow = 0; jup = 0
! lower corner
if ( cbc_iylow .eq. cbc_ixlow ) jlow = 1
! upper corner
if ( cbc_iyup .eq. cbc_ixlow ) jup = 1
! low priority corner (just avoid)
if ( cbc_ixlow .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1
if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1
endif
Proj = 0.0d0
Proj(1,1) = -1.0d0
Proj(2,2) = -1.0d0
Proj(3,3) = 1.0d0
do j=1+jlow, ny-jup
flux(:,1) = (fxi(1:nvars,2,j,1) -fxi(1:nvars,1,j,1))/dx
flux(:,2) = (fxi(1:nvars,1,j+1,2)-fxi(1:nvars,1,j,2))/dy
flux(:,3) = 0.0d0
uxloc = ux(:,1,j)
vxloc = vx(:,1,j)
uxref = ux(:,0,j)
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,j), &
flux, Proj, cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,1,j), sat, &
flux, vx(5,0,j), dx, cbc_force_matrix(:,1), cbc_ixlow, &
1, j, 0)
! rotation back
call InverseRotateXYZ(source(:,1,j), sat, flux, Proj, cbc_ixlow)
! Correct flux to the last point
source(:,1,j) = source(:,1,j) - sat(:) - flux(:,1) - flux(:,2)
enddo
endif
!-----------------------------------------------------------
!
! RIGHT BOUNDARY
!
!-----------------------------------------------------------
if ( cbc_ixup .ne. CLES_CBC_NONE ) then
jlow = 0; jup = 0
! lower corner
if ( cbc_iylow .eq. cbc_ixup ) jlow = 1
! upper corner
if ( cbc_iyup .eq. cbc_ixup ) jup = 1
! low priority corner (just avoid)
if ( cbc_ixup .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1
if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1
endif
do j=1+jlow, ny-jup
flux(:,1) = (fxi(1:nvars,nx+1,j,1)-fxi(1:nvars,nx,j,1))/dx
flux(:,2) = (fxi(1:nvars,nx,j+1,2)-fxi(1:nvars,nx,j,2))/dy
flux(:,3) = 0.0d0
call SingleCSAT(ux(:,nx,j), vx(:,nx,j), ux(:,nx+1,j), &
source(:,nx,j), sat, flux, vx(5,nx+1,j), dx, &
cbc_force_matrix(:,2), cbc_ixup, nx, j, 0)
! Correct flux to the last point
source(:,nx,j) = source(:,nx,j) - sat(:) -flux(:,1) -flux(:,2)
enddo
endif
!-----------------------------------------------------------
!
! BOTTOM BOUNDARY
!
!-----------------------------------------------------------
if ( cbc_iylow .ne. CLES_CBC_NONE ) then
ilow = 0; iup = 0
! lower corner
if ( cbc_ixlow .eq. cbc_iylow ) ilow = 1
! upper corner
if ( cbc_ixup .eq. cbc_iylow ) iup = 1
! low priority corner (just avoid)
if ( cbc_iylow .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1
if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1
endif
Proj = 0.0d0
Proj(1,2) = 1.0d0
Proj(2,1) = -1.0d0
Proj(3,3) = 1.0d0
do i=1+ilow, nx-iup
flux(:,1) = (fxi(1:nvars,i+1,1,1)-fxi(1:nvars,i,1,1))/dx
flux(:,2) = (fxi(1:nvars,i,2,2) -fxi(1:nvars,i,1,2))/dy
flux(:,3) = 0.0d0
uxloc = ux(:,i,1)
vxloc = vx(:,i,1)
uxref = ux(:,i,0)
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,1), flux, &
Proj, cbc_iylow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,i,1), sat, &
flux, vx(5,i,0), dy, cbc_force_matrix(:,3), cbc_iylow, &
i, 1, 0)
! rotation back
call InverseRotateXYZ(source(:,i,1), sat, flux, Proj, cbc_iylow)
! Correct flux to the last point
source(:,i,1) = source(:,i,1) - sat(:)- flux(:,1)- flux(:,2)
enddo
endif
!-----------------------------------------------------------
!
! TOP BOUNDARY
!
!-----------------------------------------------------------
if ( cbc_iyup .ne. CLES_CBC_NONE ) then
ilow = 0; iup = 0
! lower corner
if ( cbc_ixlow .eq. cbc_iyup ) ilow = 1
! upper corner
if ( cbc_ixup .eq. cbc_iyup ) iup = 1
! low priority corner (just avoid)
if ( cbc_iyup .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1
if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1
endif
Proj = 0.0d0
Proj(1,2) = -1.0d0
Proj(2,1) = 1.0d0
Proj(3,3) = 1.0d0
do i=1+ilow, nx-iup
flux(:,1) = (fxi(1:nvars,i+1,ny,1)-fxi(1:nvars,i,ny,1))/dx
flux(:,2) = (fxi(1:nvars,i,ny+1,2)-fxi(1:nvars,i,ny,2))/dy
flux(:,3) = 0.0d0
uxloc = ux(:,i,ny)
vxloc = vx(:,i,ny)
uxref = ux(:,i,ny+1)
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,ny), &
flux, Proj, cbc_iyup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,i,ny), sat, &
flux, vx(5,i,ny+1), dy, cbc_force_matrix(:,4), &
cbc_iyup, i, ny, 0)
! rotation back
call InverseRotateXYZ(source(:,i,ny), sat, flux, Proj, cbc_iyup)
! Correct flux to the last point
source(:,i,ny) = source(:,i,ny) - sat(:)- flux(:,1)- flux(:,2)
enddo
endif
!-----------------------------------------------------------
!
! LOWER LEFT CORNER
!
!-----------------------------------------------------------
if ( cbc_ixlow .eq. cbc_iylow .and. cbc_ixlow .ne. CLES_CBC_NONE ) then
Proj = 0.0d0
Proj(1,1) = -isqr2
Proj(1,2) = isqr2
Proj(2,1) = -isqr2
Proj(2,2) = -isqr2
Proj(3,3) = 1.0d0
flux(:,1) = (fxi(1:nvars,2,1,1)-fxi(1:nvars,1,1,1))/dx
flux(:,2) = (fxi(1:nvars,1,2,2)-fxi(1:nvars,1,1,2))/dy
flux(:,3) = 0.0d0
uxloc = ux(:,1,1)
vxloc = vx(:,1,1)
uxref = (ux(:,1,0)+ux(:,0,1))*0.5d0
dl = sqrt(dx*dy)
p = (vx(5,1,0)+vx(5,0,1))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,1), flux, &
Proj, cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,1,1), sat, flux, &
p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, 1, 0)
! rotation back
call InverseRotateXYZ(source(:,1,1), sat, flux, Proj, cbc_ixlow)
! Correct flux to the last point
source(:,1,1) = source(:,1,1) - sat(:)- flux(:,1) -flux(:,2)
endif
!-----------------------------------------------------------
!
! UPPER LEFT CORNER
!
!-----------------------------------------------------------
if ( cbc_ixlow .eq. cbc_iyup .and. cbc_ixlow .ne. CLES_CBC_NONE ) then
Proj = 0.0d0
Proj(1,1) = -isqr2
Proj(1,2) = -isqr2
Proj(2,1) = isqr2
Proj(2,2) = -isqr2
Proj(3,3) = 1.0d0
flux(:,1) = (fxi(1:nvars,2,ny,1) -fxi(1:nvars,1,ny,1))/dx
flux(:,2) = (fxi(1:nvars,1,ny+1,2)-fxi(1:nvars,1,ny,2))/dy
flux(:,3) = 0.0d0
uxloc = ux(:,1,ny)
vxloc = vx(:,1,ny)
uxref = (ux(:,1,ny+1)+ux(:,0,ny))*0.5d0
dl = sqrt(dx*dy)
p = (vx(5,1,ny+1)+vx(5,0,ny))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,ny), flux, &
Proj, cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,1,ny), sat, flux, &
p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, ny, 0)
! rotation back
call InverseRotateXYZ(source(:,1,ny), sat, flux, Proj, cbc_ixlow)
! Correct flux to the last point
source(:,1,ny) = source(:,1,ny) - sat(:)- flux(:,1) -flux(:,2)
endif
!-----------------------------------------------------------
!
! LOWER RIGHT CORNER
!
!-----------------------------------------------------------
if ( cbc_ixup .eq. cbc_iylow .and. cbc_ixup .ne. CLES_CBC_NONE ) then
Proj = 0.0d0
Proj(1,1) = isqr2
Proj(1,2) = isqr2
Proj(2,1) = -isqr2
Proj(2,2) = isqr2
Proj(3,3) = 1.0d0
flux(:,1) = (fxi(1:nvars,nx+1,1,1)-fxi(1:nvars,nx,1,1))/dx
flux(:,2) = (fxi(1:nvars,nx,2,2) -fxi(1:nvars,nx,1,2))/dy
flux(:,3) = 0.0d0
uxloc = ux(:,nx,1)
vxloc = vx(:,nx,1)
uxref = (ux(:,nx,0)+ux(:,nx+1,1))*0.5d0
dl = sqrt(dx*dy)
p = (vx(5,nx,0)+vx(5,nx+1,1))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,1), flux, &
Proj, cbc_ixup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,1), sat, flux, &
p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, 1, 0)
! rotation back
call InverseRotateXYZ(source(:,nx,1), sat, flux, Proj, cbc_ixup)
! Correct flux to the last point
source(:,nx,1) = source(:,nx,1) - sat(:)- flux(:,1) -flux(:,2)
endif
!-----------------------------------------------------------
!
! UPPER RIGHT CORNER
!
!-----------------------------------------------------------
if ( cbc_ixup .eq. cbc_iyup .and. cbc_ixup .ne. CLES_CBC_NONE ) then
Proj = 0.0d0
Proj(1,1) = isqr2
Proj(1,2) = -isqr2
Proj(2,1) = isqr2
Proj(2,2) = isqr2
Proj(3,3) = 1.0d0
flux(:,1) = (fxi(1:nvars,nx+1,ny,1)-fxi(1:nvars,nx,ny,1))/dx
flux(:,2) = (fxi(1:nvars,nx,ny+1,2)-fxi(1:nvars,nx,ny,2))/dy
flux(:,3) = 0.0d0
uxloc = ux(:,nx,ny)
vxloc = vx(:,nx,ny)
uxref = (ux(:,nx,ny+1)+ux(:,nx+1,ny))*0.5d0
dl = sqrt(dx*dy)
p = (vx(5,nx,ny+1)+vx(5,nx+1,ny))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,ny), flux, &
Proj, cbc_ixup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,ny), sat, flux, &
p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, ny, 0)
! rotation back
call InverseRotateXYZ(source(:,nx,ny), sat, flux, Proj, cbc_ixup)
! Correct flux to the last point
source(:,nx,ny) = source(:,nx,ny) - sat(:)- flux(:,1) - flux(:,2)
endif
RETURN
END SUBROUTINE TwoDAcousticBC
SUBROUTINE ThreeDAcousticBC(fxi, ux, vx, source, dcflag)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
USE Generic_Source
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION, INTENT(INOUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi,3)
! the cell centered flux
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION, INTENT(INOUT) :: source(nvars,1:nx,1:ny,1:nz)
INTEGER, INTENT(IN) :: dcflag(1:nx+1,1:ny+1,1:nz+1,3)
INTEGER :: i, j, k, ip
INTEGER :: ilow, iup, jlow, jup, klow, kup
DOUBLE PRECISION :: uxloc(ncomps), vxloc(nvars), uxref(ncomps)
DOUBLE PRECISION :: sat(nvars), flux(nvars,3), Proj(3,3), dl, p
DOUBLE PRECISION :: isqr2, pi, phi
if ( cbc_active .ne. CLES_TRUE ) return
isqr2 = 1.0d0/sqrt(2.0d0)
pi = 4.0d0*atan(1.0d0)
phi = asin(1.0d0/sqrt(3.0d0))
!-----------------------------------------------------------
!
! LEFT BOUNDARY
!
!-----------------------------------------------------------
if ( cbc_ixlow .ne. CLES_CBC_NONE ) then
jlow = 0; jup = 0
! lower corner
if ( cbc_iylow .eq. cbc_ixlow ) jlow = 1
! upper corner
if ( cbc_iyup .eq. cbc_ixlow ) jup = 1
klow = 0; kup = 0
! lower corner
if ( cbc_izlow .eq. cbc_ixlow ) klow = 1
! upper corner
if ( cbc_izup .eq. cbc_ixlow ) kup = 1
! low priority corner (just avoid)
if ( cbc_ixlow .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1
if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1
if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1
if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1
endif
Proj = 0.0d0
Proj(1,1) = -1.0d0
Proj(2,2) = -1.0d0
Proj(3,3) = 1.0d0
do k=1+klow, nz-kup
do j=1+jlow, ny-jup
flux(:,1) = (fxi(1:nvars,2,j,k,1) -fxi(1:nvars,1,j,k,1))/dx
flux(:,2) = (fxi(1:nvars,1,j+1,k,2)-fxi(1:nvars,1,j,k,2))/dy
flux(:,3) = (fxi(1:nvars,1,j,k+1,3)-fxi(1:nvars,1,j,k,3))/dz
uxloc = ux(:,1,j,k)
vxloc = vx(:,1,j,k)
uxref = ux(:,0,j,k)
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,j,k), &
flux, Proj, cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,1,j,k), sat, &
flux, vx(5,0,j,k), dx, cbc_force_matrix(:,1), cbc_ixlow, &
1, j, k)
! rotation back
call InverseRotateXYZ(source(:,1,j,k), sat, flux, Proj, &
cbc_ixlow)
! Correct flux to the last point
source(:,1,j,k) = source(:,1,j,k) - sat(:) &
-flux(:,1) -flux(:,2) -flux(:,3)
enddo
enddo
endif
!-----------------------------------------------------------
!
! RIGHT BOUNDARY
!
!-----------------------------------------------------------
if ( cbc_ixup .ne. CLES_CBC_NONE ) then
jlow = 0; jup = 0
! lower corner
if ( cbc_iylow .eq. cbc_ixup ) jlow = 1
! upper corner
if ( cbc_iyup .eq. cbc_ixup ) jup = 1
klow = 0; kup = 0
! lower corner
if ( cbc_izlow .eq. cbc_ixup ) klow = 1
! upper corner
if ( cbc_izup .eq. cbc_ixup ) kup = 1
! low priority corner (just avoid)
if ( cbc_ixup .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1
if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1
if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1
if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1
endif
do k=1+klow, nz-kup
do j=1+jlow, ny-jup
flux(:,1) = (fxi(1:nvars,nx+1,j,k,1)-fxi(1:nvars,nx,j,k,1))/dx
flux(:,2) = (fxi(1:nvars,nx,j+1,k,2)-fxi(1:nvars,nx,j,k,2))/dy
flux(:,3) = (fxi(1:nvars,nx,j,k+1,3)-fxi(1:nvars,nx,j,k,3))/dz
call SingleCSAT(ux(:,nx,j,k), vx(:,nx,j,k), ux(:,nx+1,j,k), &
source(:,nx,j,k), sat, flux, vx(5,nx+1,j,k), dx, &
cbc_force_matrix(:,2), cbc_ixup, nx, j, k)
! Correct flux to the last point
source(:,nx,j,k) = source(:,nx,j,k) - sat(:) &
- flux(:,1)- flux(:,2)- flux(:,3)
enddo
enddo
endif
!-----------------------------------------------------------
!
! BOTTOM BOUNDARY
!
!-----------------------------------------------------------
if ( cbc_iylow .ne. CLES_CBC_NONE ) then
ilow = 0; iup = 0
! lower corner
if ( cbc_ixlow .eq. cbc_iylow ) ilow = 1
! upper corner
if ( cbc_ixup .eq. cbc_iylow ) iup = 1
klow = 0; kup = 0
! lower corner
if ( cbc_izlow .eq. cbc_iylow ) klow = 1
! upper corner
if ( cbc_izup .eq. cbc_iylow ) kup = 1
! low priority corner (just avoid)
if ( cbc_iylow .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1
if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1
if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1
if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1
endif
Proj = 0.0d0
Proj(1,2) = 1.0d0
Proj(2,1) = -1.0d0
Proj(3,3) = 1.0d0
do k=1+klow, nz-kup
do i=1+ilow, nx-iup
flux(:,1) = (fxi(1:nvars,i+1,1,k,1)-fxi(1:nvars,i,1,k,1))/dx
flux(:,2) = (fxi(1:nvars,i,2,k,2) -fxi(1:nvars,i,1,k,2))/dy
flux(:,3) = (fxi(1:nvars,i,1,k+1,3)-fxi(1:nvars,i,1,k,3))/dz
uxloc = ux(:,i,1,k)
vxloc = vx(:,i,1,k)
uxref = ux(:,i,0,k)
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,1,k), &
flux, Proj, cbc_iylow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,i,1,k), sat, &
flux, vx(5,i,0,k), dy, cbc_force_matrix(:,3), cbc_iylow, &
i, 1, k)
! rotation back
call InverseRotateXYZ(source(:,i,1,k), sat, flux, Proj, &
cbc_iylow)
! Correct flux to the last point
source(:,i,1,k) = source(:,i,1,k) - sat(:)- flux(:,1)&
- flux(:,2)- flux(:,3)
enddo
enddo
endif
!-----------------------------------------------------------
!
! TOP BOUNDARY
!
!-----------------------------------------------------------
if ( cbc_iyup .ne. CLES_CBC_NONE ) then
ilow = 0; iup = 0
! lower corner
if ( cbc_ixlow .eq. cbc_iyup ) ilow = 1
! upper corner
if ( cbc_ixup .eq. cbc_iyup ) iup = 1
klow = 0; kup = 0
! lower corner
if ( cbc_izlow .eq. cbc_iyup ) klow = 1
! upper corner
if ( cbc_izup .eq. cbc_iyup ) kup = 1
! low priority corner (just avoid)
if ( cbc_iyup .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1
if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1
if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1
if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1
endif
Proj = 0.0d0
Proj(1,2) = -1.0d0
Proj(2,1) = 1.0d0
Proj(3,3) = 1.0d0
do k=1+klow, nz-kup
do i=1+ilow, nx-iup
flux(:,1) = (fxi(1:nvars,i+1,ny,k,1)-fxi(1:nvars,i,ny,k,1))/dx
flux(:,2) = (fxi(1:nvars,i,ny+1,k,2)-fxi(1:nvars,i,ny,k,2))/dy
flux(:,3) = (fxi(1:nvars,i,ny,k+1,3)-fxi(1:nvars,i,ny,k,3))/dz
uxloc = ux(:,i,ny,k)
vxloc = vx(:,i,ny,k)
uxref = ux(:,i,ny+1,k)
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,ny,k), &
flux, Proj, cbc_iyup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,i,ny,k), sat, &
flux, vx(5,i,ny+1,k), dy, cbc_force_matrix(:,4), &
cbc_iyup, i, ny, k)
! rotation back
call InverseRotateXYZ(source(:,i,ny,k), sat, flux, Proj, &
cbc_iyup)
! Correct flux to the last point
source(:,i,ny,k) = source(:,i,ny,k) - sat(:)- flux(:,1)&
- flux(:,2)- flux(:,3)
enddo
enddo
endif
!-----------------------------------------------------------
!
! BACK BOUNDARY
!
!-----------------------------------------------------------
if ( cbc_izlow .ne. CLES_CBC_NONE ) then
ilow = 0; iup = 0
! lower corner
if ( cbc_ixlow .eq. cbc_izlow ) ilow = 1
! upper corner
if ( cbc_ixup .eq. cbc_izlow ) iup = 1
jlow = 0; jup = 0
! lower corner
if ( cbc_iylow .eq. cbc_izlow ) jlow = 1
! upper corner
if ( cbc_iyup .eq. cbc_izlow ) jup = 1
! low priority corner (just avoid)
if ( cbc_izlow .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1
if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1
if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1
if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1
endif
Proj = 0.0d0
Proj(1,3) = 1.0d0
Proj(2,2) = 1.0d0
Proj(3,1) = -1.0d0
do j=1+jlow, ny-jup
do i=1+ilow, nx-iup
flux(:,1) = (fxi(1:nvars,i+1,j,1,1)-fxi(1:nvars,i,j,1,1))/dx
flux(:,2) = (fxi(1:nvars,i,j+1,1,2)-fxi(1:nvars,i,j,1,2))/dy
flux(:,3) = (fxi(1:nvars,i,j,2,3) -fxi(1:nvars,i,j,1,3))/dz
uxloc = ux(:,i,j,1)
vxloc = vx(:,i,j,1)
uxref = ux(:,i,j,0)
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,j,1), &
flux, Proj, cbc_izlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,i,j,1), sat, &
flux, vx(5,i,j,0), dz, cbc_force_matrix(:,5), cbc_izlow, &
i, j, 1)
! rotation back
call InverseRotateXYZ(source(:,i,j,1), sat, flux, Proj, &
cbc_izlow)
! Correct flux to the last point
source(:,i,j,1) = source(:,i,j,1) - sat(:)- flux(:,1)&
- flux(:,2)- flux(:,3)
enddo
enddo
endif
!-----------------------------------------------------------
!
! FRONT BOUNDARY
!
!-----------------------------------------------------------
if ( cbc_izup .ne. CLES_CBC_NONE ) then
ilow = 0; iup = 0
! lower corner
if ( cbc_ixlow .eq. cbc_izup ) ilow = 1
! upper corner
if ( cbc_ixup .eq. cbc_izup ) iup = 1
jlow = 0; jup = 0
! lower corner
if ( cbc_iylow .eq. cbc_izup ) jlow = 1
! upper corner
if ( cbc_iyup .eq. cbc_izup ) jup = 1
! low priority corner (just avoid)
if ( cbc_izup .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1
if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1
if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1
if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1
endif
Proj = 0.0d0
Proj(1,3) = -1.0d0
Proj(2,2) = 1.0d0
Proj(3,1) = 1.0d0
do j=1+jlow, ny-jup
do i=1+ilow, nx-iup
flux(:,1) = (fxi(1:nvars,i+1,j,nz,1)-fxi(1:nvars,i,j,nz,1))/dx
flux(:,2) = (fxi(1:nvars,i,j+1,nz,2)-fxi(1:nvars,i,j,nz,2))/dy
flux(:,3) = (fxi(1:nvars,i,j,nz+1,3)-fxi(1:nvars,i,j,nz,3))/dz
uxloc = ux(:,i,j,nz)
vxloc = vx(:,i,j,nz)
uxref = ux(:,i,j,nz+1)
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,j,nz), &
flux, Proj, cbc_izup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,i,j,nz), sat, &
flux, vx(5,i,j,nz+1), dz, cbc_force_matrix(:,6), &
cbc_izup, i, j, nz)
! rotation back
call InverseRotateXYZ(source(:,i,j,nz), sat, flux, Proj, &
cbc_izup)
! Correct flux to the last point
source(:,i,j,nz) = source(:,i,j,nz) - sat(:)- flux(:,1)&
- flux(:,2)- flux(:,3)
enddo
enddo
endif
!-----------------------------------------------------------
!
! LOWER LEFT EDGE
!
!-----------------------------------------------------------
if ( cbc_ixlow .eq. cbc_iylow .and. cbc_ixlow .ne. CLES_CBC_NONE ) then
klow = 0; kup = 0
! lower corner
if ( cbc_izlow .eq. cbc_ixlow ) klow = 1
! upper corner
if ( cbc_izup .eq. cbc_ixlow ) kup = 1
! low priority corner (just avoid)
if ( cbc_ixlow .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1
if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1
endif
Proj = 0.0d0
Proj(1,1) = -isqr2
Proj(1,2) = isqr2
Proj(2,1) = -isqr2
Proj(2,2) = -isqr2
Proj(3,3) = 1.0d0
dl = sqrt(dx*dy)
do k=1+klow, nz-kup
flux(:,1) = (fxi(1:nvars,2,1,k,1)-fxi(1:nvars,1,1,k,1))/dx
flux(:,2) = (fxi(1:nvars,1,2,k,2)-fxi(1:nvars,1,1,k,2))/dy
flux(:,3) = (fxi(1:nvars,1,1,k+1,3)-fxi(1:nvars,1,1,k,3))/dz
uxloc = ux(:,1,1,k)
vxloc = vx(:,1,1,k)
uxref = (ux(:,1,0,k)+ux(:,0,1,k))*0.5d0
p = (vx(5,1,0,k)+vx(5,0,1,k))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,1,k), flux, &
Proj, cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,1,1,k), sat, &
flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, 1, k)
! rotation back
call InverseRotateXYZ(source(:,1,1,k), sat, flux, Proj, cbc_ixlow)
! Correct flux to the last point
source(:,1,1,k) = source(:,1,1,k) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
enddo
endif
!-----------------------------------------------------------
!
! UPPER LEFT EDGE
!
!-----------------------------------------------------------
if ( cbc_ixlow .eq. cbc_iyup .and. cbc_ixlow .ne. CLES_CBC_NONE ) then
klow = 0; kup = 0
! lower corner
if ( cbc_izlow .eq. cbc_ixlow ) klow = 1
! upper corner
if ( cbc_izup .eq. cbc_ixlow ) kup = 1
! low priority corner (just avoid)
if ( cbc_ixlow .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1
if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1
endif
Proj = 0.0d0
Proj(1,1) = -isqr2
Proj(1,2) = -isqr2
Proj(2,1) = isqr2
Proj(2,2) = -isqr2
Proj(3,3) = 1.0d0
dl = sqrt(dx*dy)
do k=1+klow, nz-kup
flux(:,1) = (fxi(1:nvars,2,ny,k,1) -fxi(1:nvars,1,ny,k,1))/dx
flux(:,2) = (fxi(1:nvars,1,ny+1,k,2)-fxi(1:nvars,1,ny,k,2))/dy
flux(:,3) = (fxi(1:nvars,1,ny,k+1,3)-fxi(1:nvars,1,ny,k,3))/dz
uxloc = ux(:,1,ny,k)
vxloc = vx(:,1,ny,k)
uxref = (ux(:,1,ny+1,k)+ux(:,0,ny,k))*0.5d0
p = (vx(5,1,ny+1,k)+vx(5,0,ny,k))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,ny,k), flux,&
Proj, cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,1,ny,k), sat, &
flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, ny, k)
! rotation back
call InverseRotateXYZ(source(:,1,ny,k), sat, flux, Proj, &
cbc_ixlow)
! Correct flux to the last point
source(:,1,ny,k) = source(:,1,ny,k) - sat(:) &
- flux(:,1) -flux(:,2) -flux(:,3)
enddo
endif
!-----------------------------------------------------------
!
! LOWER RIGHT EDGE
!
!-----------------------------------------------------------
if ( cbc_ixup .eq. cbc_iylow .and. cbc_ixup .ne. CLES_CBC_NONE ) then
klow = 0; kup = 0
! lower corner
if ( cbc_izlow .eq. cbc_ixup ) klow = 1
! upper corner
if ( cbc_izup .eq. cbc_ixup ) kup = 1
! low priority corner (just avoid)
if ( cbc_ixup .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1
if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1
endif
Proj = 0.0d0
Proj(1,1) = isqr2
Proj(1,2) = isqr2
Proj(2,1) = -isqr2
Proj(2,2) = isqr2
Proj(3,3) = 1.0d0
dl = sqrt(dx*dy)
do k=1+klow, nz-kup
flux(:,1) = (fxi(1:nvars,nx+1,1,k,1)-fxi(1:nvars,nx,1,k,1))/dx
flux(:,2) = (fxi(1:nvars,nx,2,k,2) -fxi(1:nvars,nx,1,k,2))/dy
flux(:,3) = (fxi(1:nvars,nx,1,k+1,3)-fxi(1:nvars,nx,1,k,3))/dz
uxloc = ux(:,nx,1,k)
vxloc = vx(:,nx,1,k)
uxref = (ux(:,nx,0,k)+ux(:,nx+1,1,k))*0.5d0
p = (vx(5,nx,0,k)+vx(5,nx+1,1,k))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,1,k), flux,&
Proj, cbc_ixup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,1,k), sat, &
flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, 1, k)
! rotation back
call InverseRotateXYZ(source(:,nx,1,k), sat, flux, Proj, cbc_ixup)
! Correct flux to the last point
source(:,nx,1,k) = source(:,nx,1,k) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
enddo
endif
!-----------------------------------------------------------
!
! UPPER RIGHT EDGE
!
!-----------------------------------------------------------
if ( cbc_ixup .eq. cbc_iyup .and. cbc_ixup .ne. CLES_CBC_NONE ) then
klow = 0; kup = 0
! lower corner
if ( cbc_izlow .eq. cbc_ixup ) klow = 1
! upper corner
if ( cbc_izup .eq. cbc_ixup ) kup = 1
! low priority corner (just avoid)
if ( cbc_ixup .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_izlow .eq. CLES_CBC_INFLOW ) klow = 1
if ( cbc_izup .eq. CLES_CBC_INFLOW ) kup = 1
endif
Proj = 0.0d0
Proj(1,1) = isqr2
Proj(1,2) = -isqr2
Proj(2,1) = isqr2
Proj(2,2) = isqr2
Proj(3,3) = 1.0d0
dl = sqrt(dx*dy)
do k=1+klow, nz-kup
flux(:,1) = (fxi(1:nvars,nx+1,ny,k,1)-fxi(1:nvars,nx,ny,k,1))/dx
flux(:,2) = (fxi(1:nvars,nx,ny+1,k,2)-fxi(1:nvars,nx,ny,k,2))/dy
flux(:,3) = (fxi(1:nvars,nx,ny,k+1,3)-fxi(1:nvars,nx,ny,k,3))/dz
uxloc = ux(:,nx,ny,k)
vxloc = vx(:,nx,ny,k)
uxref = (ux(:,nx,ny+1,k)+ux(:,nx+1,ny,k))*0.5d0
p = (vx(5,nx,ny+1,k)+vx(5,nx+1,ny,k))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,ny,k), &
flux, Proj, cbc_ixup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,ny,k), sat, &
flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, ny, k)
! rotation back
call InverseRotateXYZ(source(:,nx,ny,k), sat, flux, Proj, &
cbc_ixup)
! Correct flux to the last point
source(:,nx,ny,k) = source(:,nx,ny,k) - sat(:) &
- flux(:,1) - flux(:,2) -flux(:,3)
enddo
endif
!-----------------------------------------------------------
!
! BACK LEFT EDGE
!
!-----------------------------------------------------------
if ( cbc_ixlow .eq. cbc_izlow .and. cbc_ixlow .ne. CLES_CBC_NONE ) then
jlow = 0; jup = 0
! lower corner
if ( cbc_iylow .eq. cbc_ixlow ) jlow = 1
! upper corner
if ( cbc_iyup .eq. cbc_ixlow ) jup = 1
! low priority corner (just avoid)
if ( cbc_ixlow .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1
if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1
endif
Proj = 0.0d0
Proj(1,1) = -isqr2
Proj(1,3) = isqr2
Proj(3,1) = -isqr2
Proj(3,3) = -isqr2
Proj(2,2) = 1.0d0
dl = sqrt(dx*dz)
do j=1+jlow, ny-jup
flux(:,1) = (fxi(1:nvars,2,j,1,1) -fxi(1:nvars,1,j,1,1))/dx
flux(:,2) = (fxi(1:nvars,1,j+1,1,2)-fxi(1:nvars,1,j,1,2))/dy
flux(:,3) = (fxi(1:nvars,1,j,2,3) -fxi(1:nvars,1,j,1,3))/dz
uxloc = ux(:,1,j,1)
vxloc = vx(:,1,j,1)
uxref = (ux(:,1,j,0)+ux(:,0,j,1))*0.5d0
p = (vx(5,1,j,0)+vx(5,0,j,1))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,j,1), flux, &
Proj, cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,1,j,1), sat, &
flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, j, 1)
! rotation back
call InverseRotateXYZ(source(:,1,j,1), sat, flux, Proj, cbc_ixlow)
! Correct flux to the last point
source(:,1,j,1) = source(:,1,j,1) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
enddo
endif
!-----------------------------------------------------------
!
! FRONT LEFT EDGE
!
!-----------------------------------------------------------
if ( cbc_ixlow .eq. cbc_izup .and. cbc_ixlow .ne. CLES_CBC_NONE ) then
jlow = 0; jup = 0
! lower corner
if ( cbc_iylow .eq. cbc_ixlow ) jlow = 1
! upper corner
if ( cbc_iyup .eq. cbc_ixlow ) jup = 1
! low priority corner (just avoid)
if ( cbc_ixlow .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1
if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1
endif
Proj = 0.0d0
Proj(1,1) = -isqr2
Proj(1,3) = -isqr2
Proj(3,1) = isqr2
Proj(3,3) = -isqr2
Proj(2,2) = 1.0d0
dl = sqrt(dx*dz)
do j=1+jlow, ny-jup
flux(:,1) = (fxi(1:nvars,2,j,nz,1) -fxi(1:nvars,1,j,nz,1))/dx
flux(:,2) = (fxi(1:nvars,1,j+1,nz,2)-fxi(1:nvars,1,j,nz,2))/dy
flux(:,3) = (fxi(1:nvars,1,j,nz+1,3)-fxi(1:nvars,1,j,nz,3))/dz
uxloc = ux(:,1,j,nz)
vxloc = vx(:,1,j,nz)
uxref = (ux(:,1,j,nz+1)+ux(:,0,j,nz))*0.5d0
p = (vx(5,1,j,nz+1)+vx(5,0,j,nz))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,j,nz), flux,&
Proj, cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,1,j,nz), sat, &
flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, j, nz)
! rotation back
call InverseRotateXYZ(source(:,1,j,nz), sat, flux, Proj, &
cbc_ixlow)
! Correct flux to the last point
source(:,1,j,nz) = source(:,1,j,nz) - sat(:) &
- flux(:,1) -flux(:,2) -flux(:,3)
enddo
endif
!-----------------------------------------------------------
!
! BACK RIGHT EDGE
!
!-----------------------------------------------------------
if ( cbc_ixup .eq. cbc_izlow .and. cbc_ixup .ne. CLES_CBC_NONE ) then
jlow = 0; jup = 0
! lower corner
if ( cbc_iylow .eq. cbc_ixup ) jlow = 1
! upper corner
if ( cbc_iyup .eq. cbc_ixup ) jup = 1
! low priority corner (just avoid)
if ( cbc_ixup .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1
if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1
endif
Proj = 0.0d0
Proj(1,1) = isqr2
Proj(1,3) = isqr2
Proj(3,1) = -isqr2
Proj(3,3) = isqr2
Proj(2,2) = 1.0d0
dl = sqrt(dx*dz)
do j=1+jlow, ny-jup
flux(:,1) = (fxi(1:nvars,nx+1,j,1,1)-fxi(1:nvars,nx,j,1,1))/dx
flux(:,2) = (fxi(1:nvars,nx,j+1,1,2)-fxi(1:nvars,nx,j,1,2))/dy
flux(:,3) = (fxi(1:nvars,nx,j,2,3) -fxi(1:nvars,nx,j,1,3))/dz
uxloc = ux(:,nx,j,1)
vxloc = vx(:,nx,j,1)
uxref = (ux(:,nx,j,0)+ux(:,nx+1,j,1))*0.5d0
p = (vx(5,nx,j,0)+vx(5,nx+1,j,1))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,j,1), flux,&
Proj, cbc_ixup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,j,1), sat, &
flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, j, 1)
! rotation back
call InverseRotateXYZ(source(:,nx,j,1), sat, flux, Proj, cbc_ixup)
! Correct flux to the last point
source(:,nx,j,1) = source(:,nx,j,1) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
enddo
endif
!-----------------------------------------------------------
!
! FRONT RIGHT EDGE
!
!-----------------------------------------------------------
if ( cbc_ixup .eq. cbc_izup .and. cbc_ixup .ne. CLES_CBC_NONE ) then
jlow = 0; jup = 0
! lower corner
if ( cbc_iylow .eq. cbc_ixup ) jlow = 1
! upper corner
if ( cbc_iyup .eq. cbc_ixup ) jup = 1
! low priority corner (just avoid)
if ( cbc_ixup .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_iylow .eq. CLES_CBC_INFLOW ) jlow = 1
if ( cbc_iyup .eq. CLES_CBC_INFLOW ) jup = 1
endif
Proj = 0.0d0
Proj(1,1) = isqr2
Proj(1,3) = -isqr2
Proj(3,1) = isqr2
Proj(3,3) = isqr2
Proj(2,2) = 1.0d0
dl = sqrt(dx*dz)
do j=1+jlow, ny-jup
flux(:,1) = (fxi(1:nvars,nx+1,j,nz,1)-fxi(1:nvars,nx,j,nz,1))/dx
flux(:,2) = (fxi(1:nvars,nx,j+1,nz,2)-fxi(1:nvars,nx,j,nz,2))/dy
flux(:,3) = (fxi(1:nvars,nx,j,nz+1,3)-fxi(1:nvars,nx,j,nz,3))/dz
uxloc = ux(:,nx,j,nz)
vxloc = vx(:,nx,j,nz)
uxref = (ux(:,nx,j,nz+1)+ux(:,nx+1,j,nz))*0.5d0
p = (vx(5,nx,j,nz+1)+vx(5,nx+1,j,nz))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,j,nz), &
flux, Proj, cbc_ixup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,j,nz), sat, &
flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, j, nz)
! rotation back
call InverseRotateXYZ(source(:,nx,j,nz), sat, flux, Proj, &
cbc_ixup)
! Correct flux to the last point
source(:,nx,j,nz) = source(:,nx,j,nz) - sat(:) &
- flux(:,1) - flux(:,2) -flux(:,3)
enddo
endif
!-----------------------------------------------------------
!
! BACK BOTTOM EDGE
!
!-----------------------------------------------------------
if ( cbc_iylow .eq. cbc_izlow .and. cbc_iylow .ne. CLES_CBC_NONE ) then
ilow = 0; iup = 0
! lower corner
if ( cbc_ixlow .eq. cbc_iylow ) ilow = 1
! upper corner
if ( cbc_ixup .eq. cbc_iylow ) iup = 1
! low priority corner (just avoid)
if ( cbc_iylow .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1
if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1
endif
Proj = 0.0d0
Proj(1,2) = 1.0d0
Proj(2,1) = -isqr2
Proj(2,3) = -isqr2
Proj(3,1) = -isqr2
Proj(3,3) = isqr2
dl = sqrt(dy*dz)
do i=1+ilow, nx-iup
flux(:,1) = (fxi(1:nvars,i+1,1,1,1)-fxi(1:nvars,i,1,1,1))/dx
flux(:,2) = (fxi(1:nvars,i,2,1,2) -fxi(1:nvars,i,1,1,2))/dy
flux(:,3) = (fxi(1:nvars,i,1,2,3) -fxi(1:nvars,i,1,1,3))/dz
uxloc = ux(:,i,1,1)
vxloc = vx(:,i,1,1)
uxref = (ux(:,i,1,0)+ux(:,i,0,1))*0.5d0
p = (vx(5,i,1,0)+vx(5,i,0,1))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,1,1), flux, &
Proj, cbc_iylow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,i,1,1), sat, &
flux, p, dl, cbc_force_matrix(:,3), cbc_iylow, i, 1, 1)
! rotation back
call InverseRotateXYZ(source(:,i,1,1), sat, flux, Proj, cbc_iylow)
! Correct flux to the last point
source(:,i,1,1) = source(:,i,1,1) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
enddo
endif
!-----------------------------------------------------------
!
! FRONT BOTTOM EDGE
!
!-----------------------------------------------------------
if ( cbc_iylow .eq. cbc_izup .and. cbc_iylow .ne. CLES_CBC_NONE ) then
ilow = 0; iup = 0
! lower corner
if ( cbc_ixlow .eq. cbc_iylow ) ilow = 1
! upper corner
if ( cbc_ixup .eq. cbc_iylow ) iup = 1
! low priority corner (just avoid)
if ( cbc_iylow .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1
if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1
endif
Proj = 0.0d0
Proj(1,2) = 1.0d0
Proj(2,1) = -isqr2
Proj(2,3) = isqr2
Proj(3,1) = isqr2
Proj(3,3) = isqr2
dl = sqrt(dy*dz)
do i=1+ilow, nx-iup
flux(:,1) = (fxi(1:nvars,i+1,1,nz,1)-fxi(1:nvars,i,1,nz,1))/dx
flux(:,2) = (fxi(1:nvars,i,2,nz,2) -fxi(1:nvars,i,1,nz,2))/dy
flux(:,3) = (fxi(1:nvars,i,1,nz+1,3)-fxi(1:nvars,i,1,nz,3))/dz
uxloc = ux(:,i,1,nz)
vxloc = vx(:,i,1,nz)
uxref = (ux(:,i,1,nz+1)+ux(:,i,0,nz))*0.5d0
p = (vx(5,i,1,nz+1)+vx(5,i,0,nz))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,1,nz), flux,&
Proj, cbc_iylow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,i,1,nz), sat, &
flux, p, dl, cbc_force_matrix(:,3), cbc_iylow, i, 1, nz)
! rotation back
call InverseRotateXYZ(source(:,i,1,nz), sat, flux, Proj, &
cbc_iylow)
! Correct flux to the last point
source(:,i,1,nz) = source(:,i,1,nz) - sat(:) &
- flux(:,1) -flux(:,2) -flux(:,3)
enddo
endif
!-----------------------------------------------------------
!
! BACK TOP EDGE
!
!-----------------------------------------------------------
if ( cbc_iyup .eq. cbc_izlow .and. cbc_iyup .ne. CLES_CBC_NONE ) then
ilow = 0; iup = 0
! lower corner
if ( cbc_ixlow .eq. cbc_iyup ) ilow = 1
! upper corner
if ( cbc_ixup .eq. cbc_iyup ) iup = 1
! low priority corner (just avoid)
if ( cbc_iyup .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1
if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1
endif
Proj = 0.0d0
Proj(1,2) = -1.0d0
Proj(2,1) = isqr2
Proj(2,3) = isqr2
Proj(3,1) = -isqr2
Proj(3,3) = isqr2
dl = sqrt(dy*dz)
do i=1+ilow, nx-iup
flux(:,1) = (fxi(1:nvars,i+1,ny,1,1)-fxi(1:nvars,i,ny,1,1))/dx
flux(:,2) = (fxi(1:nvars,i,ny+1,1,2)-fxi(1:nvars,i,ny,1,2))/dy
flux(:,3) = (fxi(1:nvars,i,ny,2,3) -fxi(1:nvars,i,ny,1,3))/dz
uxloc = ux(:,i,ny,1)
vxloc = vx(:,i,ny,1)
uxref = (ux(:,i,ny,0)+ux(:,i,ny+1,1))*0.5d0
p = (vx(5,i,ny,0)+vx(5,i,ny+1,1))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,ny,1), flux,&
Proj, cbc_iyup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,i,ny,1), sat, &
flux, p, dl, cbc_force_matrix(:,4), cbc_iyup, i, ny, 1)
! rotation back
call InverseRotateXYZ(source(:,i,ny,1), sat, flux, Proj, cbc_iyup)
! Correct flux to the last point
source(:,i,ny,1) = source(:,i,ny,1) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
enddo
endif
!-----------------------------------------------------------
!
! FRONT TOP EDGE
!
!-----------------------------------------------------------
if ( cbc_iyup .eq. cbc_izup .and. cbc_iyup .ne. CLES_CBC_NONE ) then
ilow = 0; iup = 0
! lower corner
if ( cbc_ixlow .eq. cbc_iyup ) ilow = 1
! upper corner
if ( cbc_ixup .eq. cbc_iyup ) iup = 1
! low priority corner (just avoid)
if ( cbc_iyup .eq. CLES_CBC_OUTFLOW ) then
if ( cbc_ixlow .eq. CLES_CBC_INFLOW ) ilow = 1
if ( cbc_ixup .eq. CLES_CBC_INFLOW ) iup = 1
endif
Proj = 0.0d0
Proj(1,2) = -1.0d0
Proj(2,1) = isqr2
Proj(2,3) = -isqr2
Proj(3,1) = isqr2
Proj(3,3) = isqr2
dl = sqrt(dy*dz)
do i=1+ilow, nx-iup
flux(:,1) = (fxi(1:nvars,i+1,ny,nz,1)-fxi(1:nvars,i,ny,nz,1))/dx
flux(:,2) = (fxi(1:nvars,i,ny+1,nz,2)-fxi(1:nvars,i,ny,nz,2))/dy
flux(:,3) = (fxi(1:nvars,i,ny,nz+1,3)-fxi(1:nvars,i,ny,nz,3))/dz
uxloc = ux(:,i,ny,nz)
vxloc = vx(:,i,ny,nz)
uxref = (ux(:,i,ny,nz+1)+ux(:,i,ny+1,nz))*0.5d0
p = (vx(5,i,ny,nz+1)+vx(5,i,ny+1,nz))*0.5d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,i,ny,nz), &
flux, Proj, cbc_iyup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,i,ny,nz), sat, &
flux, p, dl, cbc_force_matrix(:,4), cbc_iyup, i, ny, nz)
! rotation back
call InverseRotateXYZ(source(:,i,ny,nz), sat, flux, Proj, &
cbc_iyup)
! Correct flux to the last point
source(:,i,ny,nz) = source(:,i,ny,nz) - sat(:) &
- flux(:,1) - flux(:,2) -flux(:,3)
enddo
endif
!-----------------------------------------------------------
!
! LEFT-BOTTOM-BACK CORNER
!
!-----------------------------------------------------------
if ( cbc_ixlow .eq. cbc_iylow .and. cbc_ixlow .ne. CLES_CBC_NONE &
.and. cbc_iylow .eq. cbc_izlow ) then
call setproj(-3.0d0*pi/4.0d0,-phi, Proj)
flux(:,1) = (fxi(1:nvars,2,1,1,1)-fxi(1:nvars,1,1,1,1))/dx
flux(:,2) = (fxi(1:nvars,1,2,1,2)-fxi(1:nvars,1,1,1,2))/dy
flux(:,3) = (fxi(1:nvars,1,1,2,3)-fxi(1:nvars,1,1,1,3))/dz
uxloc = ux(:,1,1,1)
vxloc = vx(:,1,1,1)
uxref = (ux(:,0,1,1)+ux(:,1,0,1)+ux(:,1,1,0))/3.0d0
dl = (dx*dy*dz)**(1.0d0/3.0d0)
p = (vx(5,0,1,1)+vx(5,1,0,1)+vx(5,1,1,0))/3.0d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,1,1), flux, &
Proj, cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,1,1,1), sat, &
flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, 1, 1)
! rotation back
call InverseRotateXYZ(source(:,1,1,1), sat, flux, Proj, cbc_ixlow)
! Correct flux to the last point
source(:,1,1,1) = source(:,1,1,1) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
endif
!-----------------------------------------------------------
!
! RIGHT-BOTTOM-BACK CORNER
!
!-----------------------------------------------------------
if ( cbc_ixup .eq. cbc_iylow .and. cbc_ixup .ne. CLES_CBC_NONE &
.and. cbc_iylow .eq. cbc_izlow ) then
call setproj(-pi/4.0d0,-phi, Proj)
flux(:,1) = (fxi(1:nvars,nx+1,1,1,1)-fxi(1:nvars,nx,1,1,1))/dx
flux(:,2) = (fxi(1:nvars,nx,2,1,2) -fxi(1:nvars,nx,1,1,2))/dy
flux(:,3) = (fxi(1:nvars,nx,1,2,3) -fxi(1:nvars,nx,1,1,3))/dz
uxloc = ux(:,nx,1,1)
vxloc = vx(:,nx,1,1)
uxref = (ux(:,nx+1,1,1)+ux(:,nx,0,1)+ux(:,nx,1,0))/3.0d0
dl = (dx*dy*dz)**(1.0d0/3.0d0)
p = (vx(5,nx+1,1,1)+vx(5,nx,0,1)+vx(5,nx,1,0))/3.0d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,1,1), flux, &
Proj, cbc_ixup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,1,1), sat, &
flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, 1, 1)
! rotation back
call InverseRotateXYZ(source(:,nx,1,1), sat, flux, Proj, cbc_ixup)
! Correct flux to the last point
source(:,nx,1,1) = source(:,nx,1,1) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
endif
!-----------------------------------------------------------
!
! LEFT-TOP-BACK CORNER
!
!-----------------------------------------------------------
if ( cbc_ixlow .eq. cbc_iyup .and. cbc_ixlow .ne. CLES_CBC_NONE &
.and. cbc_iyup .eq. cbc_izlow ) then
call setproj(3.0d0*pi/4.0d0,-phi, Proj)
flux(:,1) = (fxi(1:nvars,2,ny,1,1) -fxi(1:nvars,1,ny,1,1))/dx
flux(:,2) = (fxi(1:nvars,1,ny+1,1,2)-fxi(1:nvars,1,ny,1,2))/dy
flux(:,3) = (fxi(1:nvars,1,ny,2,3) -fxi(1:nvars,1,ny,1,3))/dz
uxloc = ux(:,1,ny,1)
vxloc = vx(:,1,ny,1)
uxref = (ux(:,0,ny,1)+ux(:,1,ny+1,1)+ux(:,1,ny,0))/3.0d0
dl = (dx*dy*dz)**(1.0d0/3.0d0)
p = (vx(5,0,ny,1)+vx(5,1,ny+1,1)+vx(5,1,ny,0))/3.0d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,ny,1), flux, &
Proj, cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,1,ny,1), sat, &
flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, ny, 1)
! rotation back
call InverseRotateXYZ(source(:,1,ny,1), sat, flux, Proj, cbc_ixlow)
! Correct flux to the last point
source(:,1,ny,1) = source(:,1,ny,1) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
endif
!-----------------------------------------------------------
!
! LEFT-BOTTOM-FRONT CORNER
!
!-----------------------------------------------------------
if ( cbc_ixlow .eq. cbc_iylow .and. cbc_ixlow .ne. CLES_CBC_NONE &
.and. cbc_iylow .eq. cbc_izup ) then
call setproj(-3.0d0*pi/4.0d0,phi, Proj)
flux(:,1) = (fxi(1:nvars,2,1,nz,1) -fxi(1:nvars,1,1,nz,1))/dx
flux(:,2) = (fxi(1:nvars,1,2,nz,2) -fxi(1:nvars,1,1,nz,2))/dy
flux(:,3) = (fxi(1:nvars,1,1,nz+1,3)-fxi(1:nvars,1,1,nz,3))/dz
uxloc = ux(:,1,1,nz)
vxloc = vx(:,1,1,nz)
uxref = (ux(:,0,1,nz)+ux(:,1,0,nz)+ux(:,1,1,nz+1))/3.0d0
dl = (dx*dy*dz)**(1.0d0/3.0d0)
p = (vx(5,0,1,nz)+vx(5,1,0,nz)+vx(5,1,1,nz+1))/3.0d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,1,nz), flux, &
Proj, cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,1,1,nz), sat, &
flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, 1, nz)
! rotation back
call InverseRotateXYZ(source(:,1,1,nz), sat, flux, Proj, cbc_ixlow)
! Correct flux to the last point
source(:,1,1,nz) = source(:,1,1,nz) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
endif
!-----------------------------------------------------------
!
! RIGHT-TOP-BACK CORNER
!
!-----------------------------------------------------------
if ( cbc_ixup .eq. cbc_iyup .and. cbc_ixup .ne. CLES_CBC_NONE &
.and. cbc_iyup .eq. cbc_izlow ) then
call setproj(pi/4.0d0,-phi, Proj)
flux(:,1) = (fxi(1:nvars,nx+1,ny,1,1)-fxi(1:nvars,nx,ny,1,1))/dx
flux(:,2) = (fxi(1:nvars,nx,ny+1,1,2)-fxi(1:nvars,nx,ny,1,2))/dy
flux(:,3) = (fxi(1:nvars,nx,ny,2,3) -fxi(1:nvars,nx,ny,1,3))/dz
uxloc = ux(:,nx,ny,1)
vxloc = vx(:,nx,ny,1)
uxref = (ux(:,nx+1,ny,1)+ux(:,nx,ny+1,1)+ux(:,nx,ny,0))/3.0d0
dl = (dx*dy*dz)**(1.0d0/3.0d0)
p = (vx(5,nx+1,ny,1)+vx(5,nx,ny+1,1)+vx(5,nx,ny,0))/3.0d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,ny,1), flux, &
Proj, cbc_ixup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,ny,1), sat, &
flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, ny, 1)
! rotation back
call InverseRotateXYZ(source(:,nx,ny,1), sat, flux, Proj, cbc_ixup)
! Correct flux to the last point
source(:,nx,ny,1) = source(:,nx,ny,1) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
endif
!-----------------------------------------------------------
!
! LEFT-TOP-FRONT CORNER
!
!-----------------------------------------------------------
if ( cbc_ixlow .eq. cbc_iyup .and. cbc_ixlow .ne. CLES_CBC_NONE &
.and. cbc_iyup .eq. cbc_izup ) then
call setproj(3.0d0*pi/4.0d0,phi, Proj)
flux(:,1) = (fxi(1:nvars,2,ny,nz,1) -fxi(1:nvars,1,ny,nz,1))/dx
flux(:,2) = (fxi(1:nvars,1,ny+1,nz,2)-fxi(1:nvars,1,ny,nz,2))/dy
flux(:,3) = (fxi(1:nvars,1,ny,nz+1,3)-fxi(1:nvars,1,ny,nz,3))/dz
uxloc = ux(:,1,ny,nz)
vxloc = vx(:,1,ny,nz)
uxref = (ux(:,0,ny,nz)+ux(:,1,ny+1,nz)+ux(:,1,ny,nz+1))/3.0d0
dl = (dx*dy*dz)**(1.0d0/3.0d0)
p = (vx(5,0,ny,nz)+vx(5,1,ny+1,nz)+vx(5,1,ny,nz+1))/3.0d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,1,ny,nz), flux, &
Proj, cbc_ixlow)
call SingleCSAT(uxloc, vxloc, uxref, source(:,1,ny,nz), sat, &
flux, p, dl, cbc_force_matrix(:,1), cbc_ixlow, 1, ny, nz)
! rotation back
call InverseRotateXYZ(source(:,1,ny,nz), sat, flux, Proj, cbc_ixlow)
! Correct flux to the last point
source(:,1,ny,nz) = source(:,1,ny,nz) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
endif
!-----------------------------------------------------------
!
! RIGHT-BOTTOM-FRONT CORNER
!
!-----------------------------------------------------------
if ( cbc_ixup .eq. cbc_iylow .and. cbc_ixup .ne. CLES_CBC_NONE &
.and. cbc_iylow .eq. cbc_izup ) then
call setproj(-pi/4.0d0,phi, Proj)
flux(:,1) = (fxi(1:nvars,nx+1,1,nz,1)-fxi(1:nvars,nx,1,nz,1))/dx
flux(:,2) = (fxi(1:nvars,nx,2,nz,2) -fxi(1:nvars,nx,1,nz,2))/dy
flux(:,3) = (fxi(1:nvars,nx,1,nz+1,3)-fxi(1:nvars,nx,1,nz,3))/dz
uxloc = ux(:,nx,1,nz)
vxloc = vx(:,nx,1,nz)
uxref = (ux(:,nx+1,1,nz)+ux(:,nx,0,nz)+ux(:,nx,1,nz+1))/3.0d0
dl = (dx*dy*dz)**(1.0d0/3.0d0)
p = (vx(5,nx+1,1,nz)+vx(5,nx,0,nz)+vx(5,nx,1,nz+1))/3.0d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,1,nz), flux, &
Proj, cbc_ixup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,1,nz), sat, &
flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, 1, nz)
! rotation back
call InverseRotateXYZ(source(:,nx,1,nz), sat, flux, Proj, cbc_ixup)
! Correct flux to the last point
source(:,nx,1,nz) = source(:,nx,1,nz) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
endif
!-----------------------------------------------------------
!
! RIGHT-TOP-FRONT CORNER
!
!-----------------------------------------------------------
if ( cbc_ixup .eq. cbc_iyup .and. cbc_ixup .ne. CLES_CBC_NONE &
.and. cbc_iyup .eq. cbc_izup ) then
call setproj(pi/4.0d0,phi, Proj)
flux(:,1) = (fxi(1:nvars,nx+1,ny,nz,1)-fxi(1:nvars,nx,ny,nz,1))/dx
flux(:,2) = (fxi(1:nvars,nx,ny+1,nz,2)-fxi(1:nvars,nx,ny,nz,2))/dy
flux(:,3) = (fxi(1:nvars,nx,ny,nz+1,3)-fxi(1:nvars,nx,ny,nz,3))/dz
uxloc = ux(:,nx,ny,nz)
vxloc = vx(:,nx,ny,nz)
uxref = (ux(:,nx+1,ny,nz)+ux(:,nx,ny+1,nz)+ux(:,nx,ny,nz+1))/3.0d0
dl = (dx*dy*dz)**(1.0d0/3.0d0)
p = (vx(5,nx+1,ny,nz)+vx(5,nx,ny+1,nz)+vx(5,nx,ny,nz+1))/3.0d0
! rotate
call DirectRotateXYZ(uxloc, vxloc, uxref, source(:,nx,ny,nz), flux, &
Proj, cbc_ixup)
call SingleCSAT(uxloc, vxloc, uxref, source(:,nx,ny,nz), sat, &
flux, p, dl, cbc_force_matrix(:,2), cbc_ixup, nx, ny, nz)
! rotation back
call InverseRotateXYZ(source(:,nx,ny,nz), sat, flux, Proj, cbc_ixup)
! Correct flux to the last point
source(:,nx,ny,nz) = source(:,nx,ny,nz) - sat(:) &
- flux(:,1) -flux(:,2) - flux(:,3)
endif
RETURN
END SUBROUTINE ThreeDAcousticBC
SUBROUTINE OneDAcousticViscousBC(is, fivx, direction)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
USE Interp_coeffs
USE Generic_EvalGamma
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(INOUT) :: fivx(ixlo:ixhi)
INTEGER, INTENT(IN) :: direction, is
logical :: c1, c2
IF ( cbc_active .NE. CLES_TRUE ) RETURN
IF ( direction .EQ. CLES_CBC_XDIR ) THEN
IF ( cbc_ixlow .NE. CLES_CBC_NONE ) THEN
! Inflow
! set the flux on left-most cell wall
! so that in this boundary cell
! df/dx = ( fivx(1) - fivx(2) )/dx =0
c1 = (cbc_ixlow .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction)
c2 = (cbc_ixlow .NE. CLES_CBC_INFLOW) .and. (is .NE. direction)
if ( c1 .or. c2 ) fivx(1) = fivx(2)
ENDIF
IF ( cbc_ixup .NE. CLES_CBC_NONE ) THEN
! Inflow
! set the flux on right-most cell wall
! so that in this boundary cell
! df/dx = ( fivx(nx+1) - fivx(nx) )/dx =0
c1 = (cbc_ixup .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction)
c2 = (cbc_ixup .NE. CLES_CBC_INFLOW) .and. (is .NE. direction)
if ( c1 .or. c2 ) fivx(nx+1) = fivx(nx)
ENDIF
ENDIF
RETURN
END SUBROUTINE OneDAcousticViscousBC
SUBROUTINE TwoDAcousticViscousBC(is, fivx, direction)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
USE Interp_coeffs
USE Generic_EvalGamma
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(INOUT) :: fivx(ixlo:ixhi,iylo:iyhi)
INTEGER, INTENT(IN) :: direction, is
logical :: c1, c2
IF ( cbc_active .NE. CLES_TRUE ) RETURN
IF ( direction .EQ. CLES_CBC_XDIR ) THEN
IF ( cbc_ixlow .NE. CLES_CBC_NONE ) THEN
! Inflow
! set the flux on left-most cell wall
! so that in this boundary cell
! df/dx = ( fivx(1,:) - fivx(2,:) )/dx =0
c1 = (cbc_ixlow .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction)
c2 = (cbc_ixlow .NE. CLES_CBC_INFLOW) .and. (is .NE. direction)
if ( c1 .or. c2 ) fivx(1,:) = fivx(2,:)
ENDIF
IF ( cbc_ixup .NE. CLES_CBC_NONE ) THEN
! Inflow
! set the flux on right-most cell wall
! so that in this boundary cell
! df/dx = ( fivx(nx+1,:) - fivx(nx,:) )/dx =0
c1 = (cbc_ixup .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction)
c2 = (cbc_ixup .NE. CLES_CBC_INFLOW) .and. (is .NE. direction)
if ( c1 .or. c2 ) fivx(nx+1,:) = fivx(nx,:)
ENDIF
ENDIF
IF ( direction .EQ. CLES_CBC_YDIR ) THEN
IF ( cbc_iylow .NE. CLES_CBC_NONE ) THEN
! Inflow
! set the flux on left-most cell wall
! so that in this boundary cell
! df/dy = ( fivx(:,1) - fivx(:,2) )/dy =0
c1 = (cbc_iylow .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction)
c2 = (cbc_iylow .NE. CLES_CBC_INFLOW) .and. (is .NE. direction)
if ( c1 .or. c2 ) fivx(:,1) = fivx(:,2)
ENDIF
IF ( cbc_iyup .NE. CLES_CBC_NONE ) THEN
! Inflow
! set the flux on right-most cell wall
! so that in this boundary cell
! df/dy = ( fivx(:,ny+1) - fivx(:,ny) )/dy =0
c1 = (cbc_iyup .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction)
c2 = (cbc_iyup .NE. CLES_CBC_INFLOW) .and. (is .NE. direction)
if ( c1 .or. c2 ) fivx(:,ny+1) = fivx(:,ny)
ENDIF
ENDIF
RETURN
END SUBROUTINE TwoDAcousticViscousBC
SUBROUTINE ThreeDAcousticViscousBC(is, fivx, direction)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
USE Interp_coeffs
USE Generic_EvalGamma
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(INOUT) :: fivx(ixlo:ixhi,iylo:iyhi,izlo:izhi)
INTEGER, INTENT(IN) :: direction, is
logical :: c1, c2
IF ( cbc_active .NE. CLES_TRUE ) RETURN
IF ( direction .EQ. CLES_CBC_XDIR ) THEN
IF ( cbc_ixlow .NE. CLES_CBC_NONE ) THEN
! Inflow
! set the flux on left-most cell wall
! so that in this boundary cell
! df/dx = ( fivx(1,:,:) - fivx(2,:,:) )/dx =0
c1 = (cbc_ixlow .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction)
c2 = (cbc_ixlow .NE. CLES_CBC_INFLOW) .and. (is .NE. direction)
if ( c1 .or. c2 ) fivx(1,:,:) = fivx(2,:,:)
ENDIF
IF ( cbc_ixup .NE. CLES_CBC_NONE ) THEN
! Inflow
! set the flux on right-most cell wall
! so that in this boundary cell
! df/dx = ( fivx(nx+1,:,:) - fivx(nx,:,:) )/dx =0
c1 = (cbc_ixup .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction)
c2 = (cbc_ixup .NE. CLES_CBC_INFLOW) .and. (is .NE. direction)
if ( c1 .or. c2 ) fivx(nx+1,:,:) = fivx(nx,:,:)
ENDIF
ENDIF
IF ( direction .EQ. CLES_CBC_YDIR ) THEN
IF ( cbc_iylow .NE. CLES_CBC_NONE ) THEN
! Inflow
! set the flux on left-most cell wall
! so that in this boundary cell
! df/dy = ( fivx(:,1,:) - fivx(:,2,:) )/dy =0
c1 = (cbc_iylow .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction)
c2 = (cbc_iylow .NE. CLES_CBC_INFLOW) .and. (is .NE. direction)
if ( c1 .or. c2 ) fivx(:,1,:) = fivx(:,2,:)
ENDIF
IF ( cbc_iyup .NE. CLES_CBC_NONE ) THEN
! Inflow
! set the flux on right-most cell wall
! so that in this boundary cell
! df/dy = ( fivx(:,ny+1,:) - fivx(:,ny,:) )/dy =0
c1 = (cbc_iyup .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction)
c2 = (cbc_iyup .NE. CLES_CBC_INFLOW) .and. (is .NE. direction)
if ( c1 .or. c2 ) fivx(:,ny+1,:) = fivx(:,ny,:)
ENDIF
ENDIF
IF ( direction .EQ. CLES_CBC_ZDIR ) THEN
IF ( cbc_izlow .NE. CLES_CBC_NONE ) THEN
! Inflow
! set the flux on left-most cell wall
! so that in this boundary cell
! df/dz = ( fivx(:,:,1) - fivx(:,:,2) )/dz =0
c1 = (cbc_izlow .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction)
c2 = (cbc_izlow .NE. CLES_CBC_INFLOW) .and. (is .NE. direction)
if ( c1 .or. c2 ) fivx(:,:,1) = fivx(:,:,2)
ENDIF
IF ( cbc_izup .NE. CLES_CBC_NONE ) THEN
! Inflow
! set the flux on right-most cell wall
! so that in this boundary cell
! df/dz = ( fivx(:,:,nz+1) - fivx(:,:,nz) )/dz =0
c1 = (cbc_izup .EQ. CLES_CBC_INFLOW) .and. (is .EQ. direction)
c2 = (cbc_izup .NE. CLES_CBC_INFLOW) .and. (is .NE. direction)
if ( c1 .or. c2 ) fivx(:,:,nz+1) = fivx(:,:,nz)
ENDIF
ENDIF
RETURN
END SUBROUTINE ThreeDAcousticViscousBC
END MODULE Generic_AcousticBC
SUBROUTINE DontUseAcousticBC(ival)
use mesh
USE Generic_AcousticBC
USE method_parms
IMPLICIT NONE
include 'cles.i'
INTEGER, INTENT(IN) :: ival
integer iret
cbc_active = CLES_FALSE
if ( ALLOCATED(cbc_force_matrix) ) then
DEALLOCATE(cbc_force_matrix, stat=iret)
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'DontUseAcousticBC/cbc_force_matrix', iret)
endif
endif
RETURN
END SUBROUTINE DontUseAcousticBC
SUBROUTINE UseAcousticBC(ival)
use mesh
USE Generic_AcousticBC
USE method_parms
IMPLICIT NONE
include 'cles.i'
INTEGER, INTENT(IN) :: ival
integer iret
cbc_active = CLES_TRUE
cbc_mode = ival
ALLOCATE(cbc_force_matrix(nvars,6), stat=iret)
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'UseAcousticBC/cbc_force_matrix', iret)
endif
cbc_force_matrix(:,:) = 0.0d0
RETURN
END SUBROUTINE UseAcousticBC
SUBROUTINE SetAcousticBoundary(idir, iside, itype)
use method_parms
IMPLICIT NONE
include 'cles.i'
INTEGER, INTENT(IN) :: idir, iside, itype
if ( idir .lt. 1 .or. idir .gt. 3 ) then
print *, 'AcousticBC Error: boundary direction incorrect'
stop
endif
if ( iside .lt. CLES_CBC_LEFT .or. iside .gt. CLES_CBC_RIGHT ) then
print *, 'AcousticBC Error: boundary side incorrect'
stop
endif
cbc_direction(iside + (idir-1)*2) = itype
END SUBROUTINE SetAcousticBoundary
SUBROUTINE SetAcousticForceMatrix(idir, iside, vec)
USE Generic_AcousticBC
USE method_parms
IMPLICIT NONE
include 'cles.i'
INTEGER, INTENT(IN) :: idir, iside
DOUBLE PRECISION, INTENT(IN) :: vec(nvars)
if ( idir .lt. 1 .or. idir .gt. 3 ) then
print *, 'AcousticBC Error: boundary direction incorrect'
stop
endif
if ( iside .lt. CLES_CBC_LEFT .or. iside .gt. CLES_CBC_RIGHT ) then
print *, 'AcousticBC Error: boundary side incorrect'
stop
endif
cbc_force_matrix(:,iside + (idir-1)*2) = vec(:)
END SUBROUTINE SetAcousticForceMatrix
FUNCTION CBCIsActive()
USE Generic_AcousticBC
IMPLICIT NONE
include 'cles.i'
LOGICAL :: CBCIsActive
if ( cbc_active .eq. CLES_TRUE ) then
CBCIsActive = .TRUE.
else
CBCIsActive = .FALSE.
endif
RETURN
END FUNCTION CBCIsActive
SUBROUTINE SetAcousticImpedance(sigma)
USE Generic_AcousticBC
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: sigma
cbc_sigma = sigma
cbc_factor = cbc_sigma*(1.0d0-cbc_Mach*cbc_Mach)*0.5d0
END SUBROUTINE SetAcousticImpedance
SUBROUTINE SetAcousticMach(Ma1)
USE Generic_AcousticBC
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: Ma1
cbc_Mach = Ma1
cbc_factor = cbc_sigma*(1.0d0-cbc_Mach*cbc_Mach)*0.5d0
END SUBROUTINE SetAcousticMach