From 7ac4f2f8a705b4e0eb7bfe4253fe2d6ae5b67dc9 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 8 Mar 2017 13:20:59 -0300
Subject: [PATCH] USER_PROBLEM: Implement reconnection boundaries along X.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
src/user_problem.F90 | 135 +++++++++++++++++++++++++++++++++++++++++++
1 file changed, 135 insertions(+)
diff --git a/src/user_problem.F90 b/src/user_problem.F90
index 084a551..6e0dc35 100644
--- a/src/user_problem.F90
+++ b/src/user_problem.F90
@@ -569,7 +569,9 @@ module user_problem
! import external procedures and variables
!
use coordinates , only : im, jm, km
+ use coordinates , only : ib, ibl, ie, ieu
use equations , only : nv
+ use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
! local variables are not implicit by default
!
@@ -585,6 +587,13 @@ module user_problem
real(kind=8), dimension(1:jm) , intent(in) :: y
real(kind=8), dimension(1:km) , intent(in) :: z
real(kind=8), dimension(1:nv,1:im,1:jm,1:km), intent(inout) :: qn
+
+! local variables
+!
+ integer :: im2, im1, i, ip1, ip2
+ integer :: jm2, jm1, j, jp1, jp2
+ integer :: km2, km1, k, kp1, kp2
+ real(kind=8) :: dx, dy, dz, dxy, dxz
!
!-------------------------------------------------------------------------------
!
@@ -594,6 +603,132 @@ module user_problem
call start_timer(imb)
#endif /* PROFILE */
+! process case with magnetic field, otherwise revert to standard outflow
+!
+ if (ibx > 0) then
+
+! get the cell sizes and their ratios
+!
+ dx = x(2) - x(1)
+ dy = y(2) - y(1)
+#if NDIMS == 3
+ dz = z(2) - z(1)
+#endif /* NDIMS == 3 */
+ dxy = dx / dy
+#if NDIMS == 3
+ dxz = dx / dz
+#endif /* NDIMS == 3 */
+
+! process left and right side boundary separatelly
+!
+ if (ic == 1) then
+
+! iterate over left-side ghost layers
+!
+ do i = ibl, 1, -1
+
+! calculate neighbor cell indices
+!
+ ip1 = min(im, i + 1)
+ ip2 = min(im, i + 2)
+
+! iterate over boundary layer
+!
+ do k = kl, ku
+#if NDIMS == 3
+ km2 = max( 1, k - 2)
+ km1 = max( 1, k - 1)
+ kp1 = min(km, k + 1)
+ kp2 = min(km, k + 2)
+#endif /* NDIMS == 3 */
+ do j = jl, ju
+ jm2 = max( 1, j - 2)
+ jm1 = max( 1, j - 1)
+ jp1 = min(jm, j + 1)
+ jp2 = min(jm, j + 2)
+
+! make the normal derivative zero
+!
+ qn(1:nv,i,j,k) = qn(1:nv,ib,j,k)
+
+! prevent the inflow
+!
+ qn(ivx,i,j,k) = min(0.0d+00, qn(ivx,ib,j,k))
+
+! update the normal component of magnetic field from divergence-free condition
+!
+ qn(ibx,i,j,k) = qn(ibx,ip2,j,k) &
+ + (qn(iby,ip1,jp1,k) - qn(iby,ip1,jm1,k)) * dxy
+#if NDIMS == 3
+ qn(ibx,i,j,k) = qn(ibx,i ,j,k) &
+ + (qn(ibz,ip1,j,kp1) - qn(ibz,ip1,j,km1)) * dxz
+#endif /* NDIMS == 3 */
+ qn(ibp,i,j,k) = 0.0d+00
+ end do ! j = jl, ju
+ end do ! k = kl, ku
+ end do ! i = ibl, 1, -1
+
+ else ! ic == 1
+
+! iterate over right-side ghost layers
+!
+ do i = ieu, im
+
+! calculate neighbor cell indices
+!
+ im1 = max( 1, i - 1)
+ im2 = max( 1, i - 2)
+
+! iterate over boundary layer
+!
+ do k = kl, ku
+#if NDIMS == 3
+ km1 = max( 1, k - 1)
+ kp1 = min(km, k + 1)
+ km2 = max( 1, k - 2)
+ kp2 = min(km, k + 2)
+#endif /* NDIMS == 3 */
+ do j = jl, ju
+ jm1 = max( 1, j - 1)
+ jp1 = min(jm, j + 1)
+ jm2 = max( 1, j - 2)
+ jp2 = min(jm, j + 2)
+
+! make the normal derivative zero
+!
+ qn(1:nv,i,j,k) = qn(1:nv,ie,j,k)
+
+! prevent the inflow
+!
+ qn(ivx,i,j,k) = max(0.0d+00, qn(ivx,ie,j,k))
+
+! update the normal component of magnetic field from divergence-free condition
+!
+ qn(ibx,i,j,k) = qn(ibx,im2,j,k) &
+ + (qn(iby,im1,jm1,k) - qn(iby,im1,jp1,k)) * dxy
+#if NDIMS == 3
+ qn(ibx,i,j,k) = qn(ibx,i ,j,k) &
+ + (qn(ibz,im1,j,km1) - qn(ibz,im1,j,kp1)) * dxz
+#endif /* NDIMS == 3 */
+ qn(ibp,i,j,k) = 0.0d+00
+ end do ! j = jl, ju
+ end do ! k = kl, ku
+ end do ! i = ieu, im
+ end if ! ic == 1
+ else ! ibx > 0
+ if (ic == 1) then
+ do i = ibl, 1, -1
+ qn(1:nv,i,jl:ju,kl:ku) = qn(1:nv,ib,jl:ju,kl:ku)
+ qn(ivx ,i,jl:ju,kl:ku) = min(0.0d+00, qn(ivx,ib,jl:ju,kl:ku))
+ end do ! i = ibl, 1, -1
+ else
+ do i = ieu, im
+ qn(1:nv,i,jl:ju,kl:ku) = qn(1:nv,ie,jl:ju,kl:ku)
+ qn(ivx ,i,jl:ju,kl:ku) = max(0.0d+00, qn(ivx,ie,jl:ju,kl:ku))
+ end do ! i = ieu, im
+ end if
+ end if ! ibx > 0
+
#ifdef PROFILE
! stop accounting time for the boundary update
!