Merge branch 'master' into binaries

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

problems/binaries.in

@@ -0,0 +1,74 @@
+# problem name and parameters
+#
+problem           = "binaries"
+gamma             =  1.1d+00
+mstar             =  2.0d+00
+mcomp             =  1.0d+00
+rstar             =  5.0d-02
+rcomp             =  5.0d-02
+dstar             =  1.0d+07
+dratio            =  1.0d+02
+msstar            =  3.0d+00
+mscomp            =  1.5d+01
+tstar             =  2.0d+01
+tcomp             =  1.0d+01
+vstar             =  5.0d-01
+vcomp             =  2.5d+00
+distance          =  2.0d-01
+eccentricity      =  0.9d+00
+period            =  1.0d+02
+buni              =  1.0d-03
+
+# physics
+#
+equation_system   = "mhd"
+equation_of_state = "adi"
+
+# numerical methods
+#
+time_advance      = "ssprk(m,2)"
+stages            = 4
+riemann_solver    = "hlld"
+reconstruction    = "tvd"
+limiter           = "mc"
+fix_positivity    = "off"
+clip_extrema      = "off"
+
+# mesh parameters
+#
+xblocks           =  3
+yblocks           =  2
+zblocks           =  1
+xmin              = -4.5d+00
+xmax              =  4.5d+00
+ymin              = -3.0d+00
+ymax              =  3.0d+00
+zmin              = -1.5d+00
+zmax              =  1.5d+00
+
+# refinement control
+#
+ncells            = 16
+nghosts           = 2
+minlev            = 1
+maxlev            = 5
+
+# boundary conditions
+#
+xlbndry           = "open"
+xubndry           = "open"
+ylbndry           = "open"
+yubndry           = "open"
+zlbndry           = "open"
+zubndry           = "open"
+enable_shapes     = "on"
+
+# runtime control parameters
+#
+tmax              = 1.2d+02
+cfl               = 3.0d-01
+
+# data output control
+#
+precise_snapshots = "on"
+snapshot_interval = 1.0d-01

sources/user_problem.F90

@@ -31,6 +31,56 @@ module user_problem
 
   implicit none
 
+! flag indicating the stars are orbiting
+!
+  logical, save :: orbiting = .false.
+
+! default problem parameter values
+!
+  real(kind=8), save :: mstar    = 2.00d+00
+  real(kind=8), save :: mcomp    = 1.00d+00
+  real(kind=8), save :: rstar    = 1.00d-01
+  real(kind=8), save :: rcomp    = 1.00d-01
+  real(kind=8), save :: rbuffer  = 8.00d-02
+  real(kind=8), save :: dstar    = 1.00d+06
+  real(kind=8), save :: dratio   = 1.00d+02
+  real(kind=8), save :: msstar   = 3.00d+00
+  real(kind=8), save :: mscomp   = 1.50d+01
+  real(kind=8), save :: tstar    = 0.00d+00
+  real(kind=8), save :: tcomp    = 0.00d+00
+  real(kind=8), save :: omstar   = 0.00d+00
+  real(kind=8), save :: omcomp   = 0.00d+00
+  real(kind=8), save :: vstar    = 5.00d-01
+  real(kind=8), save :: vcomp    = 2.50d+00
+  real(kind=8), save :: xshock   = 0.00d+00
+  real(kind=8), save :: cshock   = 0.00d+00
+  real(kind=8), save :: xshift   = 0.00d+00
+  real(kind=8), save :: dist     = 2.00d-01
+  real(kind=8), save :: period   = 1.00d+02
+  real(kind=8), save :: ecc      = 9.00d-01
+  real(kind=8), save :: buni     = 1.00d-03
+
+  real(kind=8), save :: tol      = 1.00d-14
+  integer     , save :: maxit    = 20
+
+  real(kind=8), save :: dcomp, pstar, pcomp
+  real(kind=8), save :: r2stari, r2compi, r2staro, r2compo
+  real(kind=8), save :: acomp , bcomp
+  real(kind=8), save :: om
+
+  real(kind=8), save :: ms, mc
+  real(kind=8), save :: asemi, omega, ean
+  real(kind=8), save :: uxs, uxc
+  real(kind=8), save :: tprev
+
+  real(kind=8), save :: xps =  1.0d+00, yps = 0.0d+00
+  real(kind=8), save :: xpc = -1.0d+00, ypc = 0.0d+00
+  real(kind=8), save :: uys, uyc, xrs, yrs, xrc, yrc
+
+  real(kind=8)       :: ida = 1.0d+00, idm = 1.0d+00
+
+  integer, parameter :: idx = NDIMS - 1
+
   private
   public :: initialize_user_problem, finalize_user_problem
   public :: user_statistics
@@ -58,6 +108,15 @@ module user_problem
 !
   subroutine initialize_user_problem(problem, rcount, verbose, status)
 
+    use algebra    , only : quadratic
+    use constants  , only : pi2
+    use equations  , only : ipr
+    use equations  , only : adiabatic_index
+    use helpers    , only : print_section, print_parameter
+    use mesh       , only : setup_problem
+    use parameters , only : get_parameter
+    use shapes     , only : update_shapes
+
     implicit none
 
     character(len=64), intent(in)  :: problem
@@ -65,10 +124,205 @@ module user_problem
     logical          , intent(in)  :: verbose
     integer          , intent(out) :: status
 
+    character(len=32) :: orb = "off"
+
+    integer                    :: n
+    real(kind=8)               :: as, ac
+    real(kind=8), dimension(3) :: a
+    real(kind=8), dimension(2) :: x
+!
 !-------------------------------------------------------------------------------
 !
     status = 0
 
+! star masses, radia, densities and sonic Mach numbers
+!
+    call get_parameter("mstar"       , mstar )
+    call get_parameter("mcomp"       , mcomp )
+    call get_parameter("rstar"       , rstar )
+    call get_parameter("rcomp"       , rcomp )
+    call get_parameter("rbuffer"     , rbuffer)
+    call get_parameter("dstar"       , dstar )
+    call get_parameter("dratio"      , dratio)
+    call get_parameter("msstar"      , msstar)
+    call get_parameter("mscomp"      , mscomp)
+
+! wind speeds
+!
+    call get_parameter("vstar"       , vstar )
+    call get_parameter("vcomp"       , vcomp )
+
+! orbit parameters
+!
+    call get_parameter("xshift"      , xshift)
+    call get_parameter("distance"    , dist  )
+    call get_parameter("period"      , period)
+    call get_parameter("eccentricity", ecc   )
+
+! the initial shock curvature
+!
+    call get_parameter("shock_curvature", cshock)
+
+! star rotation periods
+!
+    call get_parameter("tstar"       , tstar )
+    call get_parameter("tcomp"       , tcomp )
+
+! magnetic field profile parameters
+!
+    call get_parameter("buni"        , buni  )
+
+! Kepler's equation solver parameters
+!
+    call get_parameter("orbiting"    , orb   )
+    call get_parameter("tolerance"   , tol   )
+    call get_parameter("maxit"       , maxit )
+
+! check if stars should be orbiting
+!
+    select case(trim(orb))
+    case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
+      orbiting = .true.
+    case default
+      orbiting = .false.
+    end select
+
+! indices
+!
+    ida = idx * adiabatic_index
+    idm = 5.0d-01 * NDIMS
+
+! correct rbuffer
+!
+    rbuffer = max(0.0d+00, rbuffer)
+
+! calculate the square of radia
+!
+    r2stari =  rstar**2
+    r2compi =  rcomp**2
+    r2staro = (rstar + rbuffer)**2
+    r2compo = (rcomp + rbuffer)**2
+
+! calculate densities and pressures
+!
+    dcomp = dstar / dratio
+    if (ipr > 0) then
+      pstar = (vstar / msstar)**2 * dstar / adiabatic_index
+      pcomp = (vcomp / mscomp)**2 * dcomp / adiabatic_index
+    end if
+
+! calculate the rotation speeds
+!
+    if (tstar > 0.0d+00) omstar = 1.0d+00 / tstar
+    if (tcomp > 0.0d+00) omcomp = 1.0d+00 / tcomp
+
+! calculate orbit parameters
+!
+    acomp  = dist / (1.0d+00 - ecc)
+    bcomp  = acomp * sqrt(1.0d+00 - ecc * ecc)
+    om     = (pi2 / period) / (1.0d+00 - ecc)
+
+    if (.not. orbiting) then
+
+! fix positions if the start are not orbiting
+!
+      xps    = xshift + 0.5d+00 * acomp
+      yps    = 0.0d+00
+      xpc    = xshift - 0.5d+00 * acomp
+      ypc    = 0.0d+00
+      uxs    = 0.0d+00
+      uys    = 0.0d+00
+      uxc    = 0.0d+00
+      uyc    = 0.0d+00
+
+    else
+
+! calculate initial positions and velocities
+!
+      xps    =   mcomp / (mstar + mcomp) * dist + xshift
+      xpc    = - mstar / (mstar + mcomp) * dist + xshift
+      uys    =   xps * om
+      uyc    = - xpc * om
+
+! calculate mass fractions
+!
+      ms = mcomp / (mstar + mcomp)
+      mc = mstar / (mstar + mcomp)
+
+! calculate orbit parameters
+!
+      asemi  = dist / (1.0d+00 - ecc)
+      omega  = pi2 / period
+      ean    = 0.0d+00
+
+! set the initial positions
+!
+      xrs =   ms * dist + xshift
+      yrs =   0.0d+00
+      xrc = - mc * dist + xshift
+      yrc =   0.0d+00
+
+    end if
+
+    as = 1.0d+00 / (vstar * dstar * rstar**2)
+    ac = 1.0d+00 / (vcomp * dcomp * rcomp**2)
+
+    a(1) = as * xps**2 - ac * xpc**2
+    a(2) = 2.0d+00 * (ac * xpc - as * xps)
+    a(3) = as - ac
+
+    n = quadratic(a, x)
+    if (a(3) >= 0.0d+00) then
+      xshock = minval(x)
+    else
+      xshock = maxval(x)
+    end if
+
+! set the previous time
+!
+    tprev  = 0.0d+00
+
+! proceed if no errors
+!
+    if (status == 0) then
+
+! print information about the user problem setup
+!
+      call print_section(verbose, "Parameters (* - set, otherwise calculated)")
+      call print_parameter(verbose, '(*) adiabatic_index'     , adiabatic_index)
+      call print_parameter(verbose, '(*) dratio = dstar/dcomp', dratio)
+      call print_parameter(verbose, '(*) dstar'               , dstar)
+      call print_parameter(verbose, '( ) dcomp'               , dcomp)
+      call print_parameter(verbose, '(*) vstar'               , vstar)
+      call print_parameter(verbose, '(*) vcomp'               , vcomp)
+      call print_parameter(verbose, '(*) msstar'              , msstar)
+      call print_parameter(verbose, '(*) mscomp'              , mscomp)
+      call print_parameter(verbose, '( ) pstar'               , pstar)
+      call print_parameter(verbose, '( ) pcomp'               , pcomp)
+      call print_parameter(verbose, '(*) buni'                , buni)
+      call print_parameter(verbose, '(*) rstar'               , rstar)
+      call print_parameter(verbose, '(*) rcomp'               , rcomp)
+      call print_parameter(verbose, '(*) rbuffer'             , rbuffer)
+      call print_parameter(verbose, '(*) tstar'               , tstar)
+      call print_parameter(verbose, '(*) tcomp'               , tcomp)
+      if (orbiting) then
+        call print_parameter(verbose, '(*) distance'    , dist)
+        call print_parameter(verbose, '(*) eccentricity', ecc)
+        call print_parameter(verbose, '( ) acomp'       , acomp)
+        call print_parameter(verbose, '( ) bcomp'       , bcomp)
+      else
+        call print_parameter(verbose, '(*) distance'    , acomp)
+      end if
+      call print_parameter(verbose, '( ) shock position' , xshock)
+      call print_parameter(verbose, '(*) shock curvature', cshock)
+
+    end if ! status
+
+! associate subroutines to problem setup and shape update
+!
+    setup_problem => setup_user_problem
+    update_shapes => update_user_shapes
+
 !-------------------------------------------------------------------------------
 !
   end subroutine initialize_user_problem
@@ -91,7 +345,7 @@ module user_problem
     implicit none
 
     integer, intent(out) :: status
-
+!
 !-------------------------------------------------------------------------------
 !
     status = 0
@@ -105,7 +359,8 @@ module user_problem
 ! subroutine SETUP_USER_PROBLEM:
 ! -----------------------------
 !
-!   Subroutine sets the initial conditions for the user specific problem.
+!   Subroutine sets the initial conditions for the wind interaction in
+!   the binary star problem.
 !
 !   Arguments:
 !
@@ -116,14 +371,180 @@ module user_problem
 !
   subroutine setup_user_problem(pdata)
 
-    use blocks, only : block_data
+    use blocks     , only : block_data
+    use coordinates, only : nn => bcells
+    use coordinates, only : ax, ay, az
+    use equations  , only : prim2cons, magnetized
+    use equations  , only : nv
+    use equations  , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
 
     implicit none
 
     type(block_data), pointer, intent(inout) :: pdata
 
+    integer       :: i, j, k
+    real(kind=8)  :: xs, ys, zs = 0.0d+00
+    real(kind=8)  :: xc, yc, zc = 0.0d+00
+    real(kind=8)  :: rs2, rc2, rs, rc, rd, ra
+    real(kind=8)  :: dns, prs, vxs, vys, vzs, bzs
+    real(kind=8)  :: dnc, prc, vxc, vyc, vzc, bzc
+
+    real(kind=8), dimension(nv,nn) :: q, u
+    real(kind=8), dimension(nn)    :: x, y
+#if NDIMS == 3
+    real(kind=8), dimension(nn)    :: z
+#endif /* NDIMS == 3 */
+!
 !-------------------------------------------------------------------------------
 !
+#if NDIMS == 2
+    k = 1
+#endif /* NDIMS == 2 */
+
+    x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
+    y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
+#if NDIMS == 3
+    z(:) = pdata%meta%zmin + az(pdata%meta%level,:)
+#endif /* NDIMS == 3 */
+
+! set magnetic field components
+!
+    if (magnetized) then
+
+      q(ibx,:) = 0.0d+00
+      q(iby,:) = 0.0d+00
+      q(ibz,:) = buni
+      q(ibp,:) = 0.0d+00
+
+    end if
+
+#if NDIMS == 3
+! iterate over all positions in the YZ plane
+!
+    do k = 1, nn
+
+! calculate the Z coordinates of the central and companion stars
+!
+      zs = z(k)
+      zc = z(k)
+#endif /* NDIMS == 3 */
+
+      do j = 1, nn
+
+! calculate the Y coordinates of the central and companion stars
+!
+        ys = y(j) - yps
+        yc = y(j) - ypc
+
+! distance from the line connecting both stars
+!
+#if NDIMS == 3
+        ra = y(j)**2 + z(k)**2
+#else /* NDIMS == 3 */
+        ra = y(j)**2
+#endif /* NDIMS == 3 */
+
+! sweep along the X coordinate
+!
+        do i = 1, nn
+
+! calculate the X coordinates of the central and companion stars
+!
+          xs = x(i) - xps
+          xc = x(i) - xpc
+
+! calculate the distances from the centers of the central and companion stars
+!
+#if NDIMS == 3
+          rs2 = max(xs * xs + ys * ys + zs * zs, 1.0d-16)
+          rc2 = max(xc * xc + yc * yc + zc * zc, 1.0d-16)
+#else /* NDIMS == 3 */
+          rs2 = max(xs * xs + ys * ys, 1.0d-16)
+          rc2 = max(xc * xc + yc * yc, 1.0d-16)
+#endif /* NDIMS == 3 */
+          rs  = sqrt(rs2)
+          rc  = sqrt(rc2)
+
+! calculate profiles from the central star
+!
+          rd  = max(1.0d+00, rs / rstar)
+
+          dns = dstar / rd**idx
+          if (ipr > 0) prs = pstar / rd**ida
+          if (ibz > 0) bzs = buni  / rd**idm
+
+! set the central star wind velocity
+!
+          vxs = vstar * xs / rs
+          vys = vstar * ys / rs
+          vzs = vstar * zs / rs
+
+! calculate profiles from the companion star
+!
+          rd = max(1.0d+00, rc / rcomp)
+
+          dnc = dcomp / rd**idx
+          if (ipr > 0) prc = pcomp / rd**ida
+          if (ibz > 0) bzc = buni  / rd**idm
+
+! set the companion star wind velocity
+!
+          vxc = vcomp * xc / rc
+          vyc = vcomp * yc / rc
+          vzc = vcomp * zc / rc
+
+! set the variables
+!
+          if (rs2 <= r2stari) then
+            q(idn,i) = dns
+            q(ivx,i) = vxs * rs / rstar - omstar * ys
+            q(ivy,i) = vys * rs / rstar + omstar * xs + uys
+            q(ivz,i) = vzs * rs / rstar
+            if (ipr > 0) q(ipr,i) = prs
+            if (ibz > 0) q(ibz,i) = bzs
+          else if (rc2 <= r2compi) then
+            q(idn,i) = dnc
+            q(ivx,i) = vxc * rs / rcomp - omcomp * yc
+            q(ivy,i) = vyc * rs / rcomp + omcomp * xc + uyc
+            q(ivz,i) = vzc * rs / rcomp
+            if (ipr > 0) q(ipr,i) = prc
+            if (ibz > 0) q(ibz,i) = bzc
+          else
+            if (x(i) > xshock - cshock * ra) then
+              q(idn,i) = dns
+              q(ivx,i) = vxs
+              q(ivy,i) = vys
+              q(ivz,i) = vzs
+              if (ipr > 0) q(ipr,i) = prs
+              if (ibz > 0) q(ibz,i) = bzs
+            else
+              q(idn,i) = dnc
+              q(ivx,i) = vxc
+              q(ivy,i) = vyc
+              q(ivz,i) = vzc
+              if (ipr > 0) q(ipr,i) = prc
+              if (ibz > 0) q(ibz,i) = bzc
+            end if
+          end if
+
+        end do ! i = 1, nn
+
+! convert the primitive variables to conservative ones
+!
+        call prim2cons(q(:,:), u(:,:))
+
+! copy the primitive variables to the current block
+!
+        pdata%q(:,:,j,k) = q(:,:)
+
+! copy the conserved variables to the current block
+!
+        pdata%u(:,:,j,k) = u(:,:)
+
+      end do ! j = 1, nn
+#if NDIMS == 3
+    end do ! k = 1, nn
+#endif /* NDIMS == 3 */
 
 !-------------------------------------------------------------------------------
 !
@@ -177,15 +598,194 @@ module user_problem
 !
   subroutine update_user_shapes(pdata, time, dt)
 
-    use blocks, only : block_data
+    use blocks     , only : block_data
+    use coordinates, only : nn => bcells
+    use coordinates, only : ax, ay, az
+    use equations  , only : nv
+    use equations  , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
+    use equations  , only : prim2cons
 
     implicit none
 
     type(block_data), pointer, intent(inout) :: pdata
     real(kind=8)             , intent(in)    :: time, dt
 
+    integer       :: i, j, k
+    real(kind=8)  :: xs, ys, zs
+    real(kind=8)  :: xc, yc, zc
+    real(kind=8)  :: dv, ds
+    real(kind=8)  :: sn, cs, man, res
+    real(kind=8)  :: rs2, rc2, rs, rc, rd
+
+    real(kind=8), dimension(nv) :: qs, qc
+    real(kind=8), dimension(nn) :: x, y
+#if NDIMS == 3
+    real(kind=8), dimension(nn) :: z
+#endif /* NDIMS == 3 */
+
 !-------------------------------------------------------------------------------
 !
+! if time changes, we have to recalculate the star positions
+!
+    if (time > tprev .and. orbiting) then
+
+! solve the Kepler's equation to obtain the new value of eccentric anomaly
+!
+      man = omega * time
+      res = 1.0d+00
+      i   = 1
+      do while(abs(res) >= tol .and. i <= maxit)
+        res = (ean - ecc * sin(ean) - man) / (1.0d+00 - ecc * cos(ean))
+        ean = ean - res
+        i   = i + 1
+      end do
+      if (abs(res) >= tol) then
+        print *, "Kepler's equations could not be solved!"
+        print *, "it: ", i, ", E: ", ean, ", dE: ", res
+        print *, ""
+      end if
+
+! calculate trigonometric functions of the true anomaly
+!
+      dv = 1.0d+00 - ecc * cos(ean)
+      sn = sqrt(1.0d+00 - ecc**2) * sin(ean) / dv
+      cs = (cos(ean) - ecc) / dv
+
+! calculate the position and velocity of the companion star
+!
+      ds  =   asemi * dv
+      rs  =   ms * ds
+      rc  =   mc * ds
+      xps =   rs * cs + xshift
+      yps =   rs * sn
+      xpc = - rc * cs + xshift
+      ypc = - rc * sn
+      ds  = time - tprev
+      uxs = (xps - xrs) / ds
+      uys = (yps - yrs) / ds
+      uxc = (xpc - xrc) / ds
+      uyc = (ypc - yrc) / ds
+
+! update tprev, previous positions
+!
+      tprev = time
+      xrs   = xps
+      yrs   = yps
+      xrc   = xpc
+      yrc   = ypc
+
+    end if ! time /= tprev
+
+    x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
+    y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
+#if NDIMS == 3
+    z(:) = pdata%meta%zmin + az(pdata%meta%level,:)
+#endif /* NDIMS == 3 */
+
+! set density and pressure in the stars' interiors
+!
+    qs(idn) = dstar
+    qc(idn) = dcomp
+    if (ipr > 0) then
+      qs(ipr) = pstar
+      qc(ipr) = pcomp
+    end if
+
+! set magnetic field in the stars' interiors
+!
+    if (ibx > 0) then
+      qs(ibx) = 0.0d+00
+      qs(iby) = 0.0d+00
+      qs(ibz) = buni
+      qs(ibp) = 0.0d+00
+      qc(ibx) = 0.0d+00
+      qc(iby) = 0.0d+00
+      qc(ibz) = buni
+      qc(ibp) = 0.0d+00
+    end if
+
+#if NDIMS == 3
+    do k = 1, nn
+      zs = z(k)
+      zc = z(k)
+#else /* NDIMS == 3 */
+    do k = 1, 1
+      zs = 0.0d+00
+      zc = 0.0d+00
+#endif /* NDIMS == 3 */
+
+      do j = 1, nn
+        ys = y(j) - yps
+        yc = y(j) - ypc
+
+        do i = 1, nn
+          xs = x(i) - xps
+          xc = x(i) - xpc
+
+! calculate distances from the centers of the central and companion stars
+!
+          rs2 = xs * xs + ys * ys + zs * zs
+
+          if (rs2 <= r2staro) then
+
+            rs = max(sqrt(rs2), 1.0d-16)
+            rd = max(1.0d+00, rs / rstar)
+
+            qs(idn) = dstar / rd**idx
+            if (ipr > 0) qs(ipr) = pstar / rd**ida
+            if (ibz > 0) qs(ibz) = buni  / rd**idm
+
+            if (rs2 <= r2stari) then
+              qs(ivx) = vstar * xs / rstar
+              qs(ivy) = vstar * ys / rstar
+              qs(ivz) = vstar * zs / rstar
+            else
+              qs(ivx) = vstar * xs / rs
+              qs(ivy) = vstar * ys / rs
+              qs(ivz) = vstar * zs / rs
+            end if
+
+            qs(ivx) = qs(ivx) - omstar * ys + uxs
+            qs(ivy) = qs(ivy) + omstar * xs + uys
+
+            pdata%q(:,i,j,k) = qs(:)
+
+            call prim2cons(qs(:), pdata%u(:,i,j,k))
+
+          end if
+
+          rc2 = xc * xc + yc * yc + zc * zc
+
+          if (rc2 <= r2compo) then
+
+            rc = max(sqrt(rc2), 1.0d-16)
+            rd = max(1.0d+00, rc / rcomp)
+
+            qc(idn) = dcomp / rd**idx
+            if (ipr > 0) qc(ipr) = pcomp / rd**ida
+            if (ibz > 0) qc(ibz) = buni  / rd**idm
+
+            if (rc2 <= r2compi) then
+              qc(ivx) = vcomp * xc / rcomp
+              qc(ivy) = vcomp * yc / rcomp
+              qc(ivz) = vcomp * zc / rcomp
+            else
+              qc(ivx) = vcomp * xc / rc
+              qc(ivy) = vcomp * yc / rc
+              qc(ivz) = vcomp * zc / rc
+            end if
+
+            qc(ivx) = qc(ivx) - omcomp * yc + uxc
+            qc(ivy) = qc(ivy) + omcomp * xc + uyc
+
+            pdata%q(:,i,j,k) = qc(:)
+
+            call prim2cons(qc(:), pdata%u(:,i,j,k))
+          end if
+        end do
+
+      end do
+    end do
 
 !-------------------------------------------------------------------------------
 !