USER_PROBLEM: Implement an absoption layer near Y-boundaries.

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

sources/user_problem.F90

@@ -77,6 +77,11 @@ 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 :: acut = 5.00d-01
+  real(kind=8), save :: adec = 5.00d-01
+  real(kind=8), save :: yabs = 9.00d+99
+
   integer(kind=4), save :: runit = 33
 
   logical, save :: resistive = .false.
@@ -419,6 +424,13 @@ module user_problem
       ylo  = - ydis
       yup  =   ydis
 
+! absorption parameters
+!
+      call get_parameter("aamp", aamp)
+      call get_parameter("acut", acut)
+      call get_parameter("adec", adec)
+      yabs = ymax - abs(acut)
+
     end if ! status
 
 ! proceed if no errors
@@ -456,6 +468,9 @@ 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, '(*) acut'        , acut)
+      call print_parameter(verbose, '(*) adec'        , adec)
       call print_parameter(verbose, '(*) ydistance'   , ydis)
     end if ! status
 
@@ -1031,30 +1046,88 @@ module user_problem
 !
   subroutine update_user_shapes(pdata, time, dt)
 
-! include external procedures and variables
+    use blocks     , only : block_data
-!
+    use constants  , only : pi
-    use blocks, only : block_data
+    use coordinates, only : nn => bcells
+    use coordinates, only : ymax, ay
+    use equations  , only : nv
+    use equations  , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
+    use equations  , only : prim2cons
 
-! local variables are not implicit by default
-!
     implicit none
 
-! subroutine arguments
-!
     type(block_data), pointer, intent(inout) :: pdata
     real(kind=8)             , intent(in)    :: time, dt
 
+    integer       :: j, k = 1
+    real(kind=8)  :: fl, fr
+
+    real(kind=8), dimension(nv,nn) :: q, u
+    real(kind=8), dimension(nn)    :: yc, fc
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
-! start accounting time for the shape update
-!
     call start_timer(ims)
 #endif /* PROFILE */
 
+    yc(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
+
+    if (yc(1) < -yabs .or. yc(nn) > yabs) then
+
+      fc(:) = (abs(yc(:)) - yabs) / adec
+
+      do j = 1, nn
+        if (fc(j) > 0.0d+00 .and. fc(j) < 1.0d+00) then
+          fr = aamp * sin(5.0d-01 * pi * yc(j))**2
+          fl = 1.0d+00 - fr
+
+#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(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)
+            q(ibz,:) = fl * pdata%q(ibz,:,j,k) + fr * bgui
+            q(ibp,:) = fl * pdata%q(ibp,:,j,k)
+            if (ipr > 0) q(ipr,:) = fl * pdata%q(ipr,:,j,k) + fr * pres
+
+            call prim2cons(q(:,:), u(:,:))
+
+            pdata%q(:,:,j,k) = q(:,:)
+            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
+    end if
+
 #ifdef PROFILE
-! stop accounting time for the shape update
-!
     call stop_timer(ims)
 #endif /* PROFILE */