USER_PROBLEM: Rewrite y user boundaries.
Use boundary decay time to control the time scale of the boundary inflow values to decay to the upstream values. Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
771282e54e
commit
ff63044d1a
@ -45,7 +45,7 @@ module user_problem
|
||||
real(kind=8), save :: lund = 1.00d+00
|
||||
|
||||
real(kind=8), save :: dlta = 1.00d-16
|
||||
real(kind=8), save :: blim = 1.00d+00
|
||||
real(kind=8), save :: tdec = 1.00d+00
|
||||
|
||||
integer , save :: pert = 0
|
||||
integer , save :: nper = 10
|
||||
@ -107,7 +107,8 @@ module user_problem
|
||||
#endif /* NDIMS == 3 */
|
||||
use constants , only : pi2
|
||||
use coordinates, only : ng => nghosts, ady, xlen, zlen, ymax
|
||||
use equations , only : magnetized, csnd, csnd2
|
||||
use equations , only : magnetized, csnd, csnd2, qpbnd
|
||||
use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr
|
||||
use helpers , only : print_section, print_parameter, print_message
|
||||
use mesh , only : setup_problem
|
||||
use parameters , only : get_parameter
|
||||
@ -182,7 +183,13 @@ module user_problem
|
||||
call get_parameter("bamp", bamp)
|
||||
call get_parameter("bgui", bgui)
|
||||
|
||||
call get_parameter("blimit", blim)
|
||||
call get_parameter("boundary_decay_time", tdec)
|
||||
if (tdec <= 0.0d+00) then
|
||||
if (verbose) &
|
||||
call print_message(loc, "Boundary decay time must be larger than zero!")
|
||||
status = 1
|
||||
return
|
||||
end if
|
||||
|
||||
! calculate the maximum magnetic pressure, thermal pressure from the plasma-β
|
||||
! parameters, and the sound speed in the case of isothermal equations of state
|
||||
@ -195,7 +202,18 @@ module user_problem
|
||||
valf = sqrt(2.0d+00 * pmag / dens)
|
||||
lund = valf / max(tiny(eta), eta)
|
||||
dlta = lund**(- 1.0d+00 / 3.0d+00)
|
||||
blim = max(blim, ng * ady(1))
|
||||
|
||||
! fill up the boundary values
|
||||
!
|
||||
qpbnd(idn,2,:) = dens
|
||||
qpbnd(ivx,2,:) = 0.0d+00
|
||||
qpbnd(ivy,2,:) = 0.0d+00
|
||||
qpbnd(ivz,2,:) = 0.0d+00
|
||||
qpbnd(ibx,2,:) = [ - bamp, bamp ]
|
||||
qpbnd(iby,2,:) = 0.0d+00
|
||||
qpbnd(ibz,2,:) = bgui
|
||||
qpbnd(ibp,2,:) = 0.0d+00
|
||||
if (ipr > 0) qpbnd(ipr,2,:) = pres
|
||||
|
||||
! get the geometry parameters
|
||||
!
|
||||
@ -382,7 +400,6 @@ module user_problem
|
||||
call print_parameter(verbose, '( ) S (Lundquist number)', lund)
|
||||
end if
|
||||
call print_parameter(verbose, '(*) delta (thickness)', dlta)
|
||||
call print_parameter(verbose, '(*) blim' , blim)
|
||||
call print_parameter(verbose, '(*) perturbation', perturbation)
|
||||
call print_parameter(verbose, '(*) bper' , bper)
|
||||
call print_parameter(verbose, '(*) vper' , vper)
|
||||
@ -396,6 +413,7 @@ module user_problem
|
||||
call print_parameter(verbose, '(*) ycut' , ycut)
|
||||
call print_parameter(verbose, '(*) xdec' , xdec)
|
||||
call print_parameter(verbose, '(*) ydec' , ydec)
|
||||
call print_parameter(verbose, '(*) boundary decay time', tdec)
|
||||
call print_parameter(verbose, '(*) ydistance (Vrec)', ydis)
|
||||
if (absorption) then
|
||||
call print_parameter(verbose, '(*) tabs (absorption)', tabs)
|
||||
@ -1296,8 +1314,7 @@ module user_problem
|
||||
subroutine user_boundary_y(js, il, iu, kl, ku, t, dt, x, y, z, qn)
|
||||
|
||||
use coordinates, only : nn => bcells, nb, ne, nbl, neu
|
||||
use equations , only : magnetized, nv
|
||||
use equations , only : idn, ivy, ipr, ibx, iby, ibz, ibp
|
||||
use equations , only : magnetized, nv, ivy, qpbnd
|
||||
|
||||
implicit none
|
||||
|
||||
@ -1306,16 +1323,7 @@ module user_problem
|
||||
real(kind=8), dimension(:) , intent(in) :: x, y, z
|
||||
real(kind=8), dimension(:,:,:,:), intent(inout) :: qn
|
||||
|
||||
integer :: i, im1, ip1
|
||||
integer :: j, jm1, jp1, jm2, jp2
|
||||
integer :: k
|
||||
#if NDIMS == 3
|
||||
integer :: km1, kp1
|
||||
#endif /* NDIMS == 3 */
|
||||
real(kind=8) :: dx, dy, dyx
|
||||
#if NDIMS == 3
|
||||
real(kind=8) :: dz, dyz
|
||||
#endif /* NDIMS == 3 */
|
||||
integer :: i, j, k, jm, jp
|
||||
real(kind=8) :: fl, fr
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
@ -1326,130 +1334,44 @@ module user_problem
|
||||
k = 1
|
||||
#endif /* NDIMS == 2 */
|
||||
|
||||
dx = x(2) - x(1)
|
||||
dy = y(2) - y(1)
|
||||
#if NDIMS == 3
|
||||
dz = z(2) - z(1)
|
||||
#endif /* NDIMS == 3 */
|
||||
dyx = dy / dx
|
||||
#if NDIMS == 3
|
||||
dyz = dy / dz
|
||||
#endif /* NDIMS == 3 */
|
||||
fl = exp(- dt / tdec)
|
||||
fr = 1.0d+00 - fl
|
||||
|
||||
! process left and right side boundary separatelly
|
||||
!
|
||||
if (js == 1) then
|
||||
|
||||
! iterate over left-side ghost layers
|
||||
!
|
||||
do j = nbl, 1, -1
|
||||
|
||||
! calculate neighbor cell indices
|
||||
!
|
||||
jp1 = min(nn, j + 1)
|
||||
jp2 = min(nn, j + 2)
|
||||
|
||||
! calculate variable decay coefficients
|
||||
!
|
||||
fr = (dy * (nb - j - 5.0d-01)) / blim
|
||||
fl = 1.0d+00 - fr
|
||||
|
||||
! iterate over boundary layer
|
||||
!
|
||||
jp = j + 1
|
||||
#if NDIMS == 3
|
||||
do k = kl, ku
|
||||
km1 = max( 1, k - 1)
|
||||
kp1 = min(nn, k + 1)
|
||||
#endif /* NDIMS == 3 */
|
||||
do i = il, iu
|
||||
im1 = max( 1, i - 1)
|
||||
ip1 = min(nn, i + 1)
|
||||
|
||||
! make normal derivatives zero
|
||||
!
|
||||
qn(1:nv,i,j,k) = qn(1:nv,i,nb,k)
|
||||
|
||||
! decay density and pressure to their limits
|
||||
!
|
||||
qn(idn,i,j,k) = fl * qn(idn,i,nb,k) + fr * dens
|
||||
if (ipr > 0) qn(ipr,i,j,k) = fl * qn(ipr,i,nb,k) + fr * pres
|
||||
|
||||
! decay magnetic field to its limit
|
||||
!
|
||||
qn(ibx,i,j,k) = fl * qn(ibx,i,nb,k) - fr * bamp
|
||||
qn(ibz,i,j,k) = fl * qn(ibz,i,nb,k) + fr * bgui
|
||||
|
||||
! update By from div(B)=0
|
||||
!
|
||||
qn(iby,i,j,k) = qn(iby,i,jp2,k) &
|
||||
+ (qn(ibx,ip1,jp1,k) - qn(ibx,im1,jp1,k)) * dyx
|
||||
#if NDIMS == 3
|
||||
qn(iby,i,j,k) = qn(iby,i,j ,k) &
|
||||
+ (qn(ibz,i,jp1,kp1) - qn(ibz,i,jp1,km1)) * dyz
|
||||
#endif /* NDIMS == 3 */
|
||||
qn(ibp,i,j,k) = 0.0d+00
|
||||
qn( : ,i,j,k) = fl * qn( : ,i,jp,k) + fr * qpbnd(:,2,1)
|
||||
qn(ivy,i,j,k) = qn(ivy,i,jp,k)
|
||||
end do ! i = il, iu
|
||||
#if NDIMS == 3
|
||||
end do ! k = kl, ku
|
||||
#endif /* NDIMS == 3 */
|
||||
end do ! j = nbl, 1, -1
|
||||
|
||||
else ! js = 1
|
||||
|
||||
! iterate over right-side ghost layers
|
||||
!
|
||||
do j = neu, nn
|
||||
|
||||
! calculate neighbor cell indices
|
||||
!
|
||||
jm1 = max( 1, j - 1)
|
||||
jm2 = max( 1, j - 2)
|
||||
|
||||
! calculate variable decay coefficients
|
||||
!
|
||||
fr = (dy * (j - ne - 5.0d-01)) / blim
|
||||
fl = 1.0d+00 - fr
|
||||
|
||||
! iterate over boundary layer
|
||||
!
|
||||
jm = j - 1
|
||||
#if NDIMS == 3
|
||||
do k = kl, ku
|
||||
km1 = max( 1, k - 1)
|
||||
kp1 = min(nn, k + 1)
|
||||
#endif /* NDIMS == 3 */
|
||||
do i = il, iu
|
||||
im1 = max( 1, i - 1)
|
||||
ip1 = min(nn, i + 1)
|
||||
|
||||
! make normal derivatives zero
|
||||
!
|
||||
qn(1:nv,i,j,k) = qn(1:nv,i,ne,k)
|
||||
|
||||
! decay density and pressure to their limits
|
||||
!
|
||||
qn(idn,i,j,k) = fl * qn(idn,i,ne,k) + fr * dens
|
||||
if (ipr > 0) qn(ipr,i,j,k) = fl * qn(ipr,i,ne,k) + fr * pres
|
||||
|
||||
! decay magnetic field to its limit
|
||||
!
|
||||
qn(ibx,i,j,k) = fl * qn(ibx,i,ne,k) + fr * bamp
|
||||
qn(ibz,i,j,k) = fl * qn(ibz,i,ne,k) + fr * bgui
|
||||
|
||||
! update By from div(B)=0
|
||||
!
|
||||
qn(iby,i,j,k) = qn(iby,i,jm2,k) &
|
||||
+ (qn(ibx,im1,jm1,k) - qn(ibx,ip1,jm1,k)) * dyx
|
||||
#if NDIMS == 3
|
||||
qn(iby,i,j,k) = qn(iby,i,j ,k) &
|
||||
+ (qn(ibz,i,jm1,km1) - qn(ibz,i,jm1,kp1)) * dyz
|
||||
#endif /* NDIMS == 3 */
|
||||
qn(ibp,i,j,k) = 0.0d+00
|
||||
qn( : ,i,j,k) = fl * qn( : ,i,jm,k) + fr * qpbnd(:,2,2)
|
||||
qn(ivy,i,j,k) = qn(ivy,i,jm,k)
|
||||
end do ! i = il, iu
|
||||
#if NDIMS == 3
|
||||
end do ! k = kl, ku
|
||||
#endif /* NDIMS == 3 */
|
||||
end do ! j = neu, nn
|
||||
end if ! js = 1
|
||||
else ! ibx > 0
|
||||
|
||||
else ! magnetized
|
||||
|
||||
if (js == 1) then
|
||||
do j = nbl, 1, -1
|
||||
#if NDIMS == 3
|
||||
@ -1471,7 +1393,8 @@ module user_problem
|
||||
#endif /* NDIMS == 3 */
|
||||
end do ! j = neu, nn
|
||||
end if
|
||||
end if ! ibx > 0
|
||||
|
||||
end if ! magnetized
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user