diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index 6bfbbe6..e27b46c 100644 --- a/sources/user_problem.F90 +++ b/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 blocks , only : block_data + use constants , only : pi + 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 */