orbitingcylinders.f90

self-gravitating orbiting cylinder test for spectral solver

self-gravitating orbiting cylinder test for spectral solver

Author

Manuel Jung

This part is adapted (and translated) from the PhD-thesis from Jung (2016, [25] ). It contains two different types of tests with cylinders. Both originally come from Chan (2006, [8] ):

1. A static simulation of the gravitational potential in a
   • razor-thin disk.
   • gaussian distribution in vertical direction.
2. A dynamic evolution of the cylinders.

Density-Potential Pairs

This part compares the solution of the potential by the gravitational spectral solver with an existing analytical result. Since the module only calculates the potential within the simulation area, the potential produced by the boundaries need to be neglectable. However, the potential from outside the potential can calculated with additional gravitational modules.

Assuming a density distributions around a central point

\[ \Sigma_{(r_k,\varphi_k)}(r,\varphi) = \frac{1}{2 \pi \sigma^2} \exp{\left( - \frac{R_k^2}{\sigma} \right)}, \]

with \( \sigma \) a measurement for the width of and \( R_k \) the distance of a coordinate point form the center of the distribution. For a razor thin disk, this yields the potential (see [8]):

\[ \Phi_{(r_k, \varphi_k)} = - \frac{G}{\sigma} \left(I_0(y_k)K_1(y_k) - I_1(y_k)K_0(y_k)\right), \]

where \( y_k = \frac{R_k}{2 \sigma} \). \( I_n \) and \( K_n \) are the modified Bessel functions of first and second kind with integral order (see functions ).

The density distribution in our example is

\[ \Sigma(r,\varphi) = 2 \Sigma_{(1,10^{-3})}(r,\varphi) + 0.5 \Sigma_{(1,\pi+10^{-3})} (r,\varphi) + \Sigma_{(0.9,\frac{3}{4}\pi)} (r,\varphi), \]

with the according potential

\[ \Phi(r,\varphi) = 2 \Phi_{(1,10^{-3})}(r,\varphi) + 0.5 \Phi_{(1,\pi+10^{-3})} (r,\varphi) + \Phi_{(0.9,\frac{3}{4}\pi)} (r,\varphi). \]

Additionally, the implementation is tested for a vertical gaussian distribution, with a scale height similar to \( \sigma =0.1 \). with the same setup as above. This gives us a solution for the potential

\[ \Phi(r_k, \varphi_k) = - \frac{1}{R_k} \mathrm{erf}\left( \frac{R_k}{\sqrt{2} \sigma} \right). \]

Below the relative error for the potential is shown, solved for a razor-thin disk (left) and a vertical gaussian distribution (right). For the first case we see a maximum error of \( 10^{-3} \) at steep gradients and \( 10^{-5} \) for the second case.

Self-Gravitating Rotating Cylinders

A time-dependent test from [8], which generates two cylinders with gaussian density distribution in polar grid placed opposite towards each other. A background field adds a rigid rotation. The thereby arising centrifugal force has to be balanced from the self-gravitation of the mass-distribution. This is a very sophisticated setup for the accuracy of the acceleration by self-gravitation and the general angular-momentum conservation. The setup is

\[ \varrho(r,\varphi) = \frac{10^{-2}}{\pi \left( r_{\mathrm{max}}^2 - r_{\mathrm{min}}^2 \right)} + 0.99 \left( \frac{1}{2 \pi \sigma^2} \exp{\left(- \frac{R_1}{2 \sigma^2} \right)} + \frac{1}{2 \pi \sigma} \exp{\left(- \frac{R_2}{2 \sigma^2} \right)}\right) \\ p(r,\varphi) = \frac{10^{-2}}{\pi \left( r_{\mathrm{max}}^2 - r_{\mathrm{min}}^2 \right)} + \frac{G}{2 \pi \sigma^2} \left(\mathrm{Ei}\left(- \frac{R_1}{\sigma^2} \right) - \mathrm{Ei}\left(- \frac{R_1}{2 \sigma^2} \right) + \mathrm{Ei}\left(- \frac{R_2}{\sigma^2} \right) - \mathrm{Ei}\left(- \frac{R_2}{2 \sigma^2} \right)\right), \]

and \( \mathbf{v} = r \mathbf{e}_{\varphi} \). The function \( \mathrm{Ei}(x) := \int_{-\infty}^x \frac{\exp{t}}{2} \mathrm{d}t \) is the exponential integral, \( R_1 \), \( R_2 \) are the radii from the center of the cylinders and \( \sigma = 0.1 \) a measure for spreading of the cylinders.

Additional Settings:
angular velocity: \( \Omega = 1.0 \)
heat capacity ratio: \( \gamma = \frac{5}{3} \)
resolution: \( N_r \times N_{\varphi} = 256 \times 1024 \)
simulation time: \( t_{\mathrm{sim}} = 100 \)
gravitational constant: \( G = 1.0 \)
\( r_{\mathrm{min}}, r_{\mathrm{max}} = 0.2, 1.8 \)

The images below show the results at \( t = 0, 100 \). The cylinders retain their shape very well. Within \( 100 \) orbits of the cylinders in the rotating frame there is a loss of \( 0.06 \% \) of total angular momentum. Compared with [8], we have a \( 200 \times \) better angular momentum conservation in one time step.

References:

• [8] Chan, Chi-kwan; Psaltis, Dimitrios; Özel, Feryal, 2006 "Spectral Methods for Time-dependent Studies of Accretion Flows. II. Two-dimensional Hydrodynamic Disks with Self-Gravity" http://adsabs.harvard.edu/abs/2006ApJ...645..506C
• [25] Jung, Manuel. “Multiskalensimulationen von Schwarzlochakkretion in Balkengalaxien.” PhD-thesis, Christian-Albrechts Universität Kiel, CAU Kiel, 2016. http://macau.uni-kiel.de/receive/dissertation_diss_00019153.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: orbitingcylinders.f90                                             #
!#                                                                           #
!# Copyright (C) 2013-2024                                                   #
!# Manuel Jung <mjung@astrophysik.uni-kiel.de>                               #
!# Jannes Klee <jklee@astrophysik.uni-kiel.de>                               #
!#                                                                           #
!# This program is free software; you can redistribute it and/or modify      #
!# it under the terms of the GNU General Public License as published by      #
!# the Free Software Foundation; either version 2 of the License, or (at     #
!# your option) any later version.                                           #
!#                                                                           #
!# This program is distributed in the hope that it will be useful, but       #
!# WITHOUT ANY WARRANTY; without even the implied warranty of                #
!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or        #
!# NON INFRINGEMENT.  See the GNU General Public License for more            #
!# details.                                                                  #
!#                                                                           #
!# You should have received a copy of the GNU General Public License         #
!# along with this program; if not, write to the Free Software               #
!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################
 
!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
PROGRAM orbitingcylinders
  USE fosite_mod
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! problem setup
  REAL, PARAMETER :: sigma       = 0.1
  REAL, DIMENSION(2), PARAMETER :: x1 = (/ 1.0, 0. /)
  REAL, DIMENSION(2), PARAMETER :: x2 = (/ 1.0, pi /)
  REAL, DIMENSION(2), PARAMETER :: x3 = (/ 0.9, 3./4.*pi /)
  REAL, DIMENSION(2), PARAMETER :: x4 = (/ 1.0, -pi/2. /)
 
  ! general constants
  REAL, PARAMETER :: GN          = 1.
  REAL, PARAMETER :: P0          = 1.
  REAL, PARAMETER :: P           = 10.
 
  REAL, PARAMETER :: flaring     = 0.05
  REAL, PARAMETER :: RG          = 8.31
  REAL, PARAMETER :: MU          = 6.02e-04
  REAL, PARAMETER :: GAMMA       = 5./3.
 
  REAL, PARAMETER :: RMIN        = 0.4
  REAL, PARAMETER :: RMAX        = 2.0
  INTEGER, PARAMETER :: YRES     = 128
  INTEGER, PARAMETER :: XRES     = yres/4
  INTEGER, PARAMETER :: ONUM     = p * 1
  REAL, PARAMETER    :: omega    = 1.0
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE :: Sim
  INTEGER :: green
  CHARACTER(LEN=1)  :: str
  REAL, DIMENSION(2) :: pot_err
  REAL, DIMENSION(:,:,:), POINTER :: numpot => null(), anapot => null()
  !--------------------------------------------------------------------------!
 
  DO green=1,3
    ALLOCATE(sim)
    CALL sim%InitFosite()
    ! basic configuration for all tests
    CALL makeconfig(sim,sim%config)
 
    ! set data file name and Green's function type
    WRITE(str,'(I1)') green
    CALL setattr(sim%config, "/datafile/filename",("orbitingcylinders-test" // str))
    CALL setattr(sim%config, "/sources/grav/self/green",green)
 
    SELECT CASE(green)
    CASE(1,2)
      CALL setattr(sim%config, "/timedisc/stoptime",1e-5*p0)
      CALL setattr(sim%config, "/datafile/count",1)
    CASE(3)
      CALL setattr(sim%config, "/mesh/xmin", 0.2)
      CALL setattr(sim%config, "/mesh/xmax", 1.8)
      CALL setattr(sim%config, "/mesh/inum", 64)
      CALL setattr(sim%config, "/mesh/jnum", 192)
      CALL setattr(sim%config, "/timedisc/stoptime",p*p0)
      CALL setattr(sim%config, "/datafile/count",onum)
    END SELECT
 
    CALL sim%Setup()
    CALL initdata(sim%Timedisc,sim%Mesh,sim%Physics,sim%Fluxes,sim%Sources)
    CALL sim%FirstStep()
 
    SELECT CASE(green)
    CASE(1,2)
!       CALL Sim%Datafile%WriteDataset(Sim%Mesh,Sim%Physics,Sim%Fluxes,&
!                           Sim%Timedisc,Sim%config,Sim%IO)
      ! compute max of local/global relative error
      pot_err(1) = maxval(abs((numpot(:,:,:)-anapot(:,:,:))/anapot(:,:,:)), &
         mask=sim%Mesh%without_ghost_zones%mask3d(:,:,:))
      pot_err(2) = sum(abs(numpot-anapot),mask=sim%Mesh%without_ghost_zones%mask3d(:,:,:)) &
             /sum(abs(anapot),mask=sim%Mesh%without_ghost_zones%mask3d(:,:,:))
      sim%aborted=.false.
      CALL sim%Finalize()
      print *,"max local relative error: ", pot_err(1)
      print *,"global relative error: ", pot_err(2)
    CASE(3)
      CALL sim%Run()
      CALL sim%Finalize()
    END SELECT
 
    DEALLOCATE(sim)
  END DO
 
  IF (ASSOCIATED(anapot)) DEALLOCATE(anapot)
 
 
CONTAINS
 
  SUBROUTINE makeconfig(Sim,config)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fosite), INTENT(INOUT) :: Sim
    TYPE(Dict_TYP),POINTER :: config
    !------------------------------------------------------------------------!
    ! Local variable declaration
    TYPE(Dict_TYP),POINTER :: mesh,boundary,timedisc,datafile,&
                              fluxes,physics,grav,rotframe,sources
    !------------------------------------------------------------------------!
    !------------------------------------------------------------------------!
    physics => dict( &
        "problem" / euler,&
        "mu" / mu, &
        "gamma" / gamma, &
        "units" / geometrical)
 
    fluxes => dict( &
        "fluxtype"  / kt, &!HLLC, &
        "order" / linear, &
        "variables" / primitive, &
        "limiter" / monocent, &
        "theta" / 1.2)
 
    mesh => dict( &
        "meshtype" / midpoint, &
        "geometry" / cylindrical, &
        "xmin" / rmin, &
        "xmax" / rmax, &
!        "geometry" / LOGPOLAR, &
!        "xmin" / LOG(RMIN/1.), &
!        "xmax" / LOG(RMAX/1.), &
        "gparam" / 1., &
        "omega" / omega, &
        "inum" / xres, &
        "jnum" / yres, &
        "knum" / 1, &
        "ymin" / (0.), &
        "ymax" / (2.*pi), &
        "zmin" / (0.), &
        "zmax" / (0.), &
        "decomposition" / (/ -1, 1/))!, &
!        "output/dl" / 1, &
!        "output/radius" / 1, &
!        "output/volume" / 1 )
 
    boundary => dict( &
        "western" / noslip, &
        "eastern" / noslip, &
        "southern" / periodic, &
        "northern" / periodic, &
        "bottomer" / reflecting, &
        "topper" / reflecting)
 
    grav => dict( &
        "stype" / gravity, &
        "output/accel" / 1, &
        "self/gtype" / spectral, &
        "self/green" / 1, &
        "self/sigma" / sigma)
!        "self/output/potential" / 1)
 
    rotframe => dict( &
               "stype" / rotating_frame, &
               "x"     / 0.0, &
               "y"     / 0.0, &
               "z"     / 0.0)
 
    sources => dict( &
!              "rotframe"        / rotframe, &
              "grav"            / grav)
 
    timedisc => dict( &
        "method" / ssprk, &
        "rhstype" / 1, &
!        "dtmax" / 8., &
!        "fargo" / 1, &
!        "method" / RK_FEHLBERG, &
!        "order" / 5, &
        "tol_rel" / 1.0e-3, &
        "cfl" / 0.3, &
        "stoptime" / p0, &
        "dtlimit" / 1.0e-10, &
        "maxiter" / 2000000000)
 
    datafile => dict( &
        "fileformat" / xdmf, &
        "filename" / "orbitingcylinders", &
        "count" / onum)
 
    config => dict( &
        "physics" / physics, &
        "fluxes" / fluxes, &
        "mesh" / mesh, &
        "boundary" / boundary, &
        "sources" / sources, &
        "timedisc" / timedisc, &
        "datafile" / datafile)
 
  END SUBROUTINE makeconfig
 
 
  SUBROUTINE initdata(Timedisc,Mesh,Physics,Fluxes,Sources)
#ifdef HAVE_FGSL
    USE fgsl
#endif
    USE functions, ONLY : ei,bessel_i0,bessel_i1,bessel_k0,bessel_k1,bessel_k0e
    USE sources_base_mod, ONLY : sources_base
    USE sources_gravity_mod, ONLY : sources_gravity
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
    CLASS(mesh_base),     INTENT(IN)    :: Mesh
    CLASS(physics_base),  INTENT(INOUT) :: Physics
    CLASS(fluxes_base),   INTENT(IN)    :: Fluxes
    CLASS(sources_list), ALLOCATABLE, INTENT(IN) :: Sources
    !------------------------------------------------------------------------!
    ! Local variable declaration
    INTEGER           :: i,j,k,dir,ig
    REAL,DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: r1,r2,r3
    REAL, DIMENSION(:,:,:), POINTER :: r
    CLASS(sources_base), POINTER :: sp => null()
    CLASS(gravity_base), POINTER :: gp => null()
    !------------------------------------------------------------------------!
    r => mesh%radius%bcenter
    ! get gravity spectral solver
    IF (ALLOCATED(sources)) THEN
      sp => sources%GetSourcesPointer(gravity)
    ELSE
      CALL physics%Error("orbitingcylinders::InitData","no source terms initialized")
    END IF
    IF (ASSOCIATED(sp)) THEN
      SELECT TYPE (sp)
      CLASS IS(sources_gravity)
        numpot => sp%pot%data4d(:,:,:,1)
        gp => sp%glist%GetGravityPointer(spectral)
      END SELECT
    ELSE
      CALL physics%Error("orbitingcylinders::InitData","gravity module not initialized")
    END IF
    IF (.NOT.ASSOCIATED(numpot)) &
      CALL physics%Error("orbitingcylinders::InitData","no pointer to numerical gravitational potential found")
    IF (.NOT.ASSOCIATED(gp)) &
      CALL physics%Error("orbitingcylinders::InitData","no spectral gravity solver initialized")
 
    IF (.NOT.ASSOCIATED(anapot)) ALLOCATE(anapot(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX))
 
    SELECT CASE(green)
    CASE(1)
      SELECT TYPE (pvar => timedisc%pvar)
      CLASS IS(statevector_euler)
        r1(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x1(1),x1(2))
        r2(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x3(1),x3(2))
        r3(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x4(1),x4(2))
        pvar%density%data3d(:,:,:) &
        = 2.0 / (2.*pi*sigma**2) * exp(-sqrt(r1(:,:,:))/sigma) &
        + 0.5 / (2.*pi*sigma**2) * exp(-sqrt(r2(:,:,:))/sigma) &
        + 1.0 / (2.*pi*sigma**2) * exp(-sqrt(r3(:,:,:))/sigma)
        pvar%pressure%data1d(:) = 1.0
        pvar%velocity%data1d(:) = 0.0
 
        IF (ASSOCIATED(anapot)) THEN
          r1(:,:,:) = sqrt(rsq1(mesh%radius%faces(:,:,:,1),mesh%curv%faces(:,:,:,1,2),x1(1),x1(2))) / (2.*sigma)
          r2(:,:,:) = sqrt(rsq1(mesh%radius%faces(:,:,:,1),mesh%curv%faces(:,:,:,1,2),x3(1),x3(2))) / (2.*sigma)
          r3(:,:,:) = sqrt(rsq1(mesh%radius%faces(:,:,:,1),mesh%curv%faces(:,:,:,1,2),x4(1),x4(2))) / (2.*sigma)
          DO k=mesh%KGMIN,mesh%KGMAX
            DO j=mesh%JGMIN,mesh%JGMAX
              DO i=mesh%IGMIN,mesh%IGMAX
                anapot(i,j,k) &
#ifdef HAVE_FGSL
                = -2.0 * r1(i,j,k)/sigma*(fgsl_sf_bessel_ic0(r1(i,j,k))*fgsl_sf_bessel_kc1(r1(i,j,k))&
                    -fgsl_sf_bessel_ic1(r1(i,j,k))*fgsl_sf_bessel_kc0(r1(i,j,k))) &
                -0.5 * r2(i,j,k)/sigma*(fgsl_sf_bessel_ic0(r2(i,j,k))*fgsl_sf_bessel_kc1(r2(i,j,k))&
                    -fgsl_sf_bessel_ic1(r2(i,j,k))*fgsl_sf_bessel_kc0(r2(i,j,k))) &
                -1.0 * r3(i,j,k)/sigma*(fgsl_sf_bessel_ic0(r3(i,j,k))*fgsl_sf_bessel_kc1(r3(i,j,k))&
                    -fgsl_sf_bessel_ic1(r3(i,j,k))*fgsl_sf_bessel_kc0(r3(i,j,k)))
#else
                = -2.0 * r1(i,j,k)/sigma*(bessel_i0(r1(i,j,k))*bessel_k1(r1(i,j,k))&
                    -bessel_i1(r1(i,j,k))*bessel_k0(r1(i,j,k))) &
                -0.5 * r2(i,j,k)/sigma*(bessel_i0(r2(i,j,k))*bessel_k1(r2(i,j,k))&
                    -bessel_i1(r2(i,j,k))*bessel_k0(r2(i,j,k))) &
                -1.0 * r3(i,j,k)/sigma*(bessel_i0(r3(i,j,k))*bessel_k1(r3(i,j,k))&
                    -bessel_i1(r3(i,j,k))*bessel_k0(r3(i,j,k)))
#endif
              END DO
            END DO
          END DO
        END IF
      END SELECT
      SELECT TYPE(sp)
      CLASS IS(sources_gravity)
        CALL sp%UpdateGravity(mesh,physics,fluxes,timedisc%pvar,0.0,0.0)
      END SELECT
    CASE(2)
      SELECT TYPE (pvar => timedisc%pvar)
      CLASS IS(statevector_euler)
        r1(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x1(1),x1(2))
        r2(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x3(1),x3(2))
        r3(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x4(1),x4(2))
        pvar%density%data3d(:,:,:) &
        = 2.0 / (2.*pi*sigma**2) * exp(-r1(:,:,:)/(2.*sigma**2)) &
        + 0.5 / (2.*pi*sigma**2) * exp(-r2(:,:,:)/(2.*sigma**2)) &
        + 1.0 / (2.*pi*sigma**2) * exp(-r3(:,:,:)/(2.*sigma**2))
        pvar%pressure%data1d(:) = 1.0
        pvar%velocity%data1d(:) = 0.0
 
        IF (ASSOCIATED(anapot)) THEN
          r1(:,:,:) = rsq1(mesh%radius%faces(:,:,:,1),mesh%curv%faces(:,:,:,1,2),x1(1),x1(2))
          r2(:,:,:) = rsq1(mesh%radius%faces(:,:,:,1),mesh%curv%faces(:,:,:,1,2),x3(1),x3(2))
          r3(:,:,:) = rsq1(mesh%radius%faces(:,:,:,1),mesh%curv%faces(:,:,:,1,2),x4(1),x4(2))
          anapot(:,:,:) &
          = - 2.0 / sqrt(r1(:,:,:)) * erf(sqrt(r1(:,:,:) / 2.) / sigma) &
            - 0.5 / sqrt(r2(:,:,:)) * erf(sqrt(r2(:,:,:) / 2.) / sigma) &
            - 1.0 / sqrt(r3(:,:,:)) * erf(sqrt(r3(:,:,:) / 2.) / sigma)
        END IF
      END SELECT
      SELECT TYPE(sp)
      CLASS IS(sources_gravity)
        CALL sp%UpdateGravity(mesh,physics,fluxes,timedisc%pvar,0.0,0.0)
      END SELECT
    CASE(3)
      SELECT TYPE (pvar => timedisc%pvar)
      CLASS IS(statevector_euler)
        r1(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x1(1),x1(2))
        r2(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x2(1),x2(2))
        pvar%density%data3d(:,:,:)  = 0.02/(3.2*pi) &
          + 0.99 * (  exp(-0.5*(r1(:,:,:)/sigma**2))/(2.*pi*sigma**2) &
                  + exp(-0.5*(r2(:,:,:)/sigma**2))/(2.*pi*sigma**2))
 
        pvar%pressure%data3d(:,:,:) = 0.02/(3.2*pi) &
          + gn / (2.*pi*sigma**2) &
             * (  ei(-(r1(:,:,:)/sigma**2)) - ei(-0.5*(r1(:,:,:)/sigma**2)) &
              + ei(-(r2(:,:,:)/sigma**2)) - ei(-0.5*(r2(:,:,:)/sigma**2)))
 
        pvar%velocity%data4d(:,:,:,1) = 0.0
        pvar%velocity%data4d(:,:,:,2) = r(:,:,:)*(1.0 - omega)
      END SELECT
    END SELECT
 
    DO dir=west,east
      SELECT TYPE(bound => timedisc%Boundary%boundary(dir)%p)
      CLASS IS (boundary_noslip)
        DO k=mesh%KMIN,mesh%KMAX
          DO j=mesh%JMIN,mesh%JMAX
            DO ig=1,mesh%GNUM
              SELECT CASE (dir)
              CASE(west)
                  i = mesh%IMIN-ig
              CASE(east)
                  i = mesh%IMAX+ig
              END SELECT
              bound%data(ig,j,k,physics%YVELOCITY) &
                = timedisc%pvar%data4d(i,j,k,physics%YVELOCITY)
            END DO
          END DO
        END DO
      END SELECT
    END DO
 
    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
  END SUBROUTINE initdata
 
  FUNCTION rsq(Mesh,r,x) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(mesh_base) :: Mesh
    REAL, DIMENSION(2) :: x
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: r
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: res
    !------------------------------------------------------------------------!
    INTENT(IN) :: x
    !------------------------------------------------------------------------!
    res = r**2 + x(1)**2 - 2.*r*x(1) * cos(mesh%curv%bcenter(:,:,:,2)-x(2))
  END FUNCTION rsq
 
  ELEMENTAL FUNCTION rsq1(r,phi,r0,phi0) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN) :: r, phi, r0, phi0
    REAL :: res
    !------------------------------------------------------------------------!
    res = r**2 + r0**2 - 2.*r*r0 * cos(phi-phi0)
  END FUNCTION rsq1
 
END PROGRAM orbitingcylinders