USER_PROBLEM: Implement absorption/diffusion layer near Y boundaries.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>

sources/user_problem.F90

@@ -77,14 +77,16 @@ module user_problem
   real(kind=8), save :: ylo  = -9.00d+99
   real(kind=8), save :: yup  =  9.00d+99
 
-  real(kind=8), save :: aamp = 1.00d-02
+  real(kind=8), save :: tabs = 1.00d+00
-  real(kind=8), save :: acut = 5.00d-01
+  real(kind=8), save :: adif = 5.00d-01
-  real(kind=8), save :: adec = 5.00d-01
+  real(kind=8), save :: acut = 1.00d+00
+  real(kind=8), save :: adec = 1.00d+00
   real(kind=8), save :: yabs = 9.00d+99
 
   integer(kind=4), save :: runit = 33
 
-  logical, save :: resistive = .false.
+  logical, save :: absorption = .false.
+  logical, save :: resistive  = .false.
 
 ! flag indicating if the gravitational source term is enabled
 !
@@ -152,6 +154,7 @@ module user_problem
 ! local variables
 !
     character(len=64) :: perturbation = "noise", append = "off", fname
+    character(len=64) :: enable_absorption = "off"
     logical           :: flag
     integer           :: n, kd
     real(kind=8)      :: thh, fc, ka, ydis = 9.00d+99
@@ -426,7 +429,15 @@ module user_problem
 
 ! absorption parameters
 !
-      call get_parameter("aamp", aamp)
+      call get_parameter("enable_shapes", enable_absorption)
+      select case(trim(enable_absorption))
+      case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
+        absorption = .true.
+      case default
+        absorption = .false.
+      end select
+      call get_parameter("tabs", tabs)
+      call get_parameter("adif", adif)
       call get_parameter("acut", acut)
       call get_parameter("adec", adec)
       yabs = ymax - abs(acut)
@@ -468,10 +479,14 @@ 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, '(*) aamp'        , aamp)
+      call print_parameter(verbose, '(*) ydistance (Vrec)', ydis)
-      call print_parameter(verbose, '(*) acut'        , acut)
+      if (absorption) then
-      call print_parameter(verbose, '(*) adec'        , adec)
+        call print_parameter(verbose, '(*) tabs (absorption)', tabs)
-      call print_parameter(verbose, '(*) ydistance'   , ydis)
+        call print_parameter(verbose, '(*) adif (diffusion)' , adif)
+        call print_parameter(verbose, '(*) acut'             , acut)
+        call print_parameter(verbose, '(*) adec'             , adec)
+        call print_parameter(verbose, '( ) yabs'             , yabs)
+      end if
     end if ! status
 
 #ifdef PROFILE
@@ -1049,10 +1064,11 @@ module user_problem
     use blocks     , only : block_data
     use constants  , only : pi
     use coordinates, only : nn => bcells
-    use coordinates, only : ymax, ay
+    use coordinates, only : ymax, ay, adx, ady, adz
     use equations  , only : nv
     use equations  , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
     use equations  , only : prim2cons
+    use operators  , only : laplace
 
     implicit none
 
@@ -1060,10 +1076,16 @@ module user_problem
     real(kind=8)             , intent(in)    :: time, dt
 
     integer       :: j, k = 1
-    real(kind=8)  :: fl, fr
+    real(kind=8)  :: fl, fr, fd, fa, fb
 
-    real(kind=8), dimension(nv,nn) :: q, u
+    real(kind=8), dimension(3)        :: dh = 1.0d+00
-    real(kind=8), dimension(nn)    :: yc, fc
+    real(kind=8), dimension(nn)       :: yc, fc
+    real(kind=8), dimension(nv,nn)    :: q, u
+#if NDIMS == 3
+    real(kind=8), dimension(nn,nn,nn) :: q2
+#else /* NDIMS == 3 */
+    real(kind=8), dimension(nn,nn, 1) :: q2
+#endif /* NDIMS == 3 */
 
 !-------------------------------------------------------------------------------
 !
@@ -1075,19 +1097,35 @@ module user_problem
 
     if (yc(1) < -yabs .or. yc(nn) > yabs) then
 
+      dh(1) = adx(pdata%meta%level)
+      dh(2) = ady(pdata%meta%level)
+#if NDIMS == 3
+      dh(3) = adz(pdata%meta%level)
+#endif /* NDIMS == 3 */
+
+      fa = dt / tabs
+      fb = 5.0d-01 * adif * minval(dh(1:NDIMS))**2
+
+      call laplace(dh, pdata%q(ivy,:,:,:), q2)
+
       fc(:) = (abs(yc(:)) - yabs) / adec
 
       do j = 1, nn
-        if (fc(j) > 0.0d+00 .and. fc(j) < 1.0d+00) then
+        if (fc(j) >= 0.0d+00) then
-          fr = aamp * sin(5.0d-01 * pi * yc(j))**2
+          if (fc(j) <= 1.0d+00) then
+            fr = fa * sin(5.0d-01 * pi * yc(j))**2
+          else
+            fr = fa
+          end if
           fl = 1.0d+00 - fr
+          fd = fr * fb
 
 #if NDIMS == 3
           do k = 1, nn
 #endif /* NDIMS == 3 */
             q(idn,:) = fl * pdata%q(idn,:,j,k) + fr * dens
             q(ivx,:) = fl * pdata%q(ivx,:,j,k)
-            q(ivy,:) = pdata%q(ivy,:,j,k)
+            q(ivy,:) =      pdata%q(ivy,:,j,k) + fd * q2(:,j,k)
             q(ivz,:) = fl * pdata%q(ivz,:,j,k)
             q(ibx,:) = fl * pdata%q(ibx,:,j,k) + fr * sign(bamp, yc(j))
             q(iby,:) = fl * pdata%q(iby,:,j,k)
@@ -1101,27 +1139,6 @@ module user_problem
             pdata%u(:,:,j,k) = u(:,:)
 #if NDIMS == 3
           end do
-#endif /* NDIMS == 3 */
-        else if (fc(j) >= 1.0d+00) then
-#if NDIMS == 3
-          do k = 1, nn
-#endif /* NDIMS == 3 */
-            q(idn,:) = dens
-            q(ivx,:) = 0.0d+00
-            q(ivy,:) = pdata%q(ivy,:,j,k)
-            q(ivz,:) = 0.0d+00
-            q(ibx,:) = sign(bamp, yc(j))
-            q(iby,:) = 0.0d+00
-            q(ibz,:) = bgui
-            q(ibp,:) = 0.0d+00
-            if (ipr > 0) q(ipr,:) = pres
-
-            call prim2cons(q(:,:), u(:,:))
-
-            pdata%q(:,:,j,k) = q(:,:)
-            pdata%u(:,:,j,k) = u(:,:)
-#if NDIMS == 3
-          end do
 #endif /* NDIMS == 3 */
         end if
       end do