solutions.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: solutions.f90                                                     #
!#                                                                           #
!# Copyright (C) 2013 Manuel Jung <mjung@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.                 #
!#                                                                           #
!#############################################################################

!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!

MODULE solutions
  USE roots, ONLY: getroot
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  PRIVATE
    REAL :: rho_l,u_l,p_l,a_l,b_l,c_l,&
            rho_r,u_r,p_r,a_r,b_r,c_r,gamma
    REAL :: gam_p1, gam_m1, r0, rt, vmin, vmax, &
            vshock, n1, n2, n3, n4, n5
    REAL :: r, xi, vxi, gxi, zxi, zgxi, xacc
    REAL :: e0, rho0, p1
  !--------------------------------------------------------------------------!
  PUBLIC :: &
    riemann, &
    sedov
  !--------------------------------------------------------------------------!
CONTAINS
  SUBROUTINE riemann(x0,gamma_,rho_l_,u_l_,p_l_,rho_r_,u_r_,p_r_,t,x,pvar)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL                :: gamma_,t,x0,rho_l_,u_l_,p_l_,rho_r_,u_r_,p_r_
    REAL,DIMENSION(:)   :: x
    REAL,DIMENSION(:,:) :: pvar
    !------------------------------------------------------------------------!
    INTEGER :: i
    Real    :: rho_lstar, rho_rstar, p_star, u_star, &
               c_lstar, c_rstar,           &
               s_l, s_hl, s_tl, s_r, s_hr, s_tr
    INTEGER :: error
    !------------------------------------------------------------------------!
    gamma = gamma_ !1.4                  !heat capacity ratio

    rho_l = rho_l_
    rho_r = rho_r_

    u_l = u_l_
    u_r = u_r_

    p_l = p_l_
    p_r = p_r_

    ! sound speeds of far regions
    c_l = sqrt(gamma*p_l/rho_l)
    c_r = sqrt(gamma*p_r/rho_r)

    a_l = 2.0/((gamma+1.0)*rho_l)
    a_r = 2.0/((gamma+1.0)*rho_r)

    b_l = p_l*(gamma-1.0)/(gamma+1.0)
    b_r = p_r*(gamma-1.0)/(gamma+1.0)

    ! calculate p_star by findung the root of function f
    ! upper limit is just a somewhat big number. Is this really always
    ! sufficient?
    CALL getroot(f, 0., 2.e+3, p_star, error)
    !p_star = GetRoot_andbjo(f, 0., 1.E+4,1.E-6)

    u_star = 0.5*((u_l + u_r) + ( f_x(p_star, rho_r, u_r, p_r, a_r, b_r, c_r) &
                                  - f_x(p_star, rho_l, u_l, p_l, a_l, b_l, c_l)))

    IF (p_star > p_l) THEN
      !left shock
      rho_lstar = rho_l * ( ((p_star/p_l) + (gamma-1.0)/(gamma+1.0)) &
                            / (((gamma-1.0)/(gamma+1.0))*(p_star/p_l) + 1.0))
    ELSE
      !left rarefaction
      rho_lstar = rho_l * (p_star/p_l)**(1.0/gamma)
    END IF

    IF (p_star > p_r) THEN
      !right shock
      rho_rstar = rho_r * ( ((p_star/p_r) + (gamma-1.0)/(gamma+1.0)) &
                           /(((gamma-1.0)/(gamma+1.0))*(p_star/p_r) + 1.0))
    ELSE
      !right rarefaction
      rho_rstar = rho_r * (p_star/p_r)**(1.0/gamma)
    END IF

    ! sound speeds of star regions
    c_lstar = c_l * (p_star/p_l)**((gamma-1.0)/(2.0*gamma))
    c_rstar = c_r * (p_star/p_r)**((gamma-1.0)/(2.0*gamma))

    ! shock speeds
    s_l = u_l - c_l * sqrt( ((gamma+1.0)/(2.0*gamma))*(p_star/p_l) + (gamma-1.0)/(2.0*gamma))
    s_r = u_r + c_r * sqrt( ((gamma+1.0)/(2.0*gamma))*(p_star/p_r) + (gamma-1.0)/(2.0*gamma))

    ! head speed of rarefraction waves
    s_hl = u_l - c_l
    s_hr = u_r + c_r
    ! tail speed of rarefraction waves
    s_tl = u_star - c_lstar
    s_tr = u_star + c_rstar

    DO i=1,SIZE(x)
      CALL sample(x(i),pvar(i,1),pvar(i,2),pvar(i,3))
    END DO

    CONTAINS
      SUBROUTINE sample(x,rho,u,p)
        IMPLICIT NONE
        !----------------------------------------------------------------------!
        REAL        :: x, rho, u, p
        !----------------------------------------------------------------------!
        REAL        :: S
        !----------------------------------------------------------------------!
        INTENT(IN)  :: x
        INTENT(OUT) :: rho, u, p
        !----------------------------------------------------------------------!

        !Sampling the solution
        s = (x-x0)/t
        IF (s <= u_star) THEN
          !left side of contact
          IF (p_star > p_l) THEN
            !left shock
            IF (s < s_l) THEN
              !far left region
              rho = rho_l
              u   = u_l
              p   = p_l
            ELSE
              !left star region
              rho = rho_lstar
              u   = u_star
              p   = p_star
            END IF
          ELSE
            !left fan
            IF (s < s_hl) THEN
              !far left region
              rho = rho_l
              u   = u_l
              p   = p_l
            ELSE IF (s > s_tl) THEN
              !left star region
              rho = rho_lstar
              u   = u_star
              p   = p_star
            ELSE
              !left fan region
              rho = rho_l * ( 2.0/(gamma+1.0) + ((gamma-1.0)/((gamma+1.0)*c_l)) &
                *(u_l - (x-x0)/t) ) ** (2.0/(gamma-1.0))
              u   = 2.0/(gamma+1.0) * (c_l + (gamma-1.0)*u_l/2.0 + (x-x0)/t)
              p   = p_l * ( 2.0/(gamma+1.0) + ((gamma-1.0)/((gamma+1.0)*c_l)) &
                * (u_l - (x-x0)/t) ) ** (2.0*gamma/(gamma-1.0))
            END IF
          END IF
        ELSE
          !right side of contact
          IF (p_star > p_r) THEN
            !right shock
            IF (s > s_r) THEN
              rho = rho_r
              u   = u_r
              p   = p_r
            ELSE
              rho = rho_rstar
              u   = u_star
              p   = p_star
            END IF
          ELSE
            !right fan
            IF (s > s_hr) THEN
              rho = rho_r
              u   = u_r
              p   = p_r
            ELSE IF (s < s_tr) THEN
              rho = rho_rstar
              u   = u_star
              p   = p_star
            ELSE
              rho = rho_r * ( 2.0/(gamma+1.0) - ((gamma-1.0)/((gamma+1.0)*c_r)) &
                * (u_r - (x-x0)/t) ) ** (2.0/(gamma-1.0))
              u = 2.0/(gamma+1.0) * ( -c_r + (gamma-1.0)*u_r/2.0 + (x-x0)/t)
              p = p_r * ( 2.0/(gamma+1.0) - ((gamma-1.0)/((gamma+1.0)*c_r)) &
                * (u_r - (x-x0)/t) ) ** (2.0*gamma/(gamma-1.0))
            END IF
          END IF
        END IF
      END SUBROUTINE sample
  END SUBROUTINE riemann

  SUBROUTINE sedov(gamma_, E0_, rho0_, P1_, time, dim, x, pvar)
    IMPLICIT NONE
    !----------------------------------------------------------------------!
    REAL, INTENT(IN)                    :: gamma_, e0_, rho0_, p1_, time
    INTEGER, INTENT(IN)                 :: dim
    REAL, DIMENSION(:), INTENT(IN)      :: x
    REAL, DIMENSION(:,:), INTENT(INOUT) :: pvar
    !----------------------------------------------------------------------!
    REAL, PARAMETER                     :: eps = 1.0d-09
    INTEGER, PARAMETER                  :: maxit = 1.0d+08
    INTEGER                             :: i
    !----------------------------------------------------------------------!
    gamma = gamma_
    e0    = e0_
    rho0  = rho0_
    p1    = p1_
    gam_p1 = gamma + 1.0
    gam_m1 = gamma - 1.0

    SELECT CASE(dim)
    CASE(2)
      r0 = 1.0
      rt = r0*(e0*time**2/rho0)**0.25
      vmin = 1.0 / gamma
      vmax = 2.0 / gamma
      vshock = 0.5 *rt/time
      n1 = -2.0
      n2 = 2.0*gam_m1/gamma
      n3 = 1.0/gamma
      n4 = -n1 / (2.0-gamma)
      n5 = -2.0 / (2.0-gamma)
    CASE(3)
      r0 = 1.033 ! for gamma=7/5, see Padmanabhan: Theo. Astro., Vol 1,  p.409
      rt = r0*(e0*time**2/rho0)**0.2
      vmin = 1./gamma
      vmax = 5./(3.*gamma-1.)
      vshock = 2.*rt/(5.*time)
      n1 = -(13.*gamma**2.-7.*gamma+12.) / ((3.*gamma-1.)*(2.*gamma+1.))
      n2 = 5.*gam_m1 / (2.*gamma+1.)
      n3 = 3. / (2.*gamma+1.)
      n4 = -n1 / (2.-gamma)
      n5 = -2. / (2.-gamma)
    END SELECT
  ! accuracy
  xacc = eps
  
  ! main loop
  DO i=1, SIZE(x)
     r = x(i)
     xi = r/rt
     IF (xi.LE.1.0) THEN
        ! Newton-Raphson to solve the implicit equation
        IF (dim.EQ.2) THEN
           vxi = getroot_test(funcd2d,vmin,vmax*0.99,xacc)
           gxi = gam_p1/gam_m1*(gam_p1/gam_m1*(gamma*vxi-1.0))**n3 &
                *(gam_p1/2.0*(2.0-gamma*vxi))**n4 &
                *(gam_p1/gam_m1*(1.0-vxi))**n5
        ELSE
           vxi = getroot_test(funcd3d,vmin,vmax*0.99,xacc)
           gxi = gam_p1/gam_m1*(gam_p1/gam_m1*(gamma*vxi-1.))**n3 &
                *(gam_p1/(7.-gamma)*(5.-(3.*gamma-1.)*vxi))**n4 &
                *(gam_p1/gam_m1*(1.-vxi))**n5
        END IF
        zxi = gamma*gam_m1*(1.-vxi)*vxi*vxi/(2.*(gamma*vxi-1.))
!        IF (xi.GT.0.01) THEN
!           zgxi = 0.5*gam*gam_p1*(gam_p1/gam_m1)**(n3+n5) * vxi**2 &
!                *(gam*vxi-1.)**(n3-1.) * (1.-vxi)**(n5+1.) &
!                *(gam_p1/(7.-gam)*(5.-(3.*gam-1.)*vxi))**n4
!        ELSE
!           zgxi = 0.5*gam*gam_m1*(1.-vxi)*vxi**2 * xi**(5/n2*(n3-1.))
!        END IF
!        CALL funcd2D(vxi,fv,dfv)
!        WRITE (*,"(5(ES15.7))") xi,vxi,gxi,zxi,fv
!        STOP
     ELSE
        vxi = 0.
        gxi = 1.
        zxi = ((2.+dim)*time/(2.*r))**2 * gamma*p1/rho0
        zgxi= zxi
     END IF
     ! set primitive variables
     pvar(i,1) = rho0*gxi
     pvar(i,2) = vshock*xi*vxi
     pvar(i,3) = pvar(i,1)/gamma*zxi*(2.*r/((2.+dim)*time))**2
!     pvar(i,3) = Control%rho0/gam*(2.*r/(5.*Control%time))**2 * zgxi
  END DO

  CONTAINS

  FUNCTION getroot_test(funcd,x1,x2,xacc) RESULT(root)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN) :: x1,x2,xacc
    REAL :: root
    !------------------------------------------------------------------------!
    INTERFACE
       SUBROUTINE funcd(x,fx,dfx)
         IMPLICIT NONE
         REAL, INTENT(IN)  :: x
         REAL, INTENT(OUT) :: fx,dfx
       END SUBROUTINE funcd
    END INTERFACE
    !------------------------------------------------------------------------!
    REAL    :: fm,dfm,fl,dfl,fr,dfr
    REAL    :: xm,xl,xr,dx
    INTEGER :: i
    !------------------------------------------------------------------------!
    ! compute left and right function values
    xl = min(x1,x2)
    xr = max(x1,x2)    
    CALL funcd(xl,fl,dfl)
    CALL funcd(xr,fr,dfr)
    ! check if root is within the interval [x1,x2]
    IF (fl*fr.GT.0.0) THEN
       WRITE (*,*) xl, xr,fl, fr,gamma, "Error: f(x1)*f(x2) should be < 0, aborting!"
       stop
    END IF
    ! main loop
    DO i=1,maxit
  !     WRITE (*,"(4(ES20.12))") xl,fl,xr,fr
       ! regular falsi
       dx = fl*(xl-xr)/(fl-fr)
       xm = xl - dx
       root = xm
       CALL funcd(xm,fm,dfm)
       ! check abort criteron
       IF (abs(fm).LT.xacc) THEN
          EXIT
       END IF
       IF (fm*fl.GT.0.0) THEN
          xl=xm
          fl=fm
       ELSE
          xr=xm
          fr=fm
       END IF
    END DO
    IF (i.GT.maxit) THEN
       WRITE (*,*) "WARNING: limit of iterations exceeded!"
    END IF
  END FUNCTION getroot_test

  END SUBROUTINE sedov


  ! the root of this function gives p_star
  !REAL FUNCTION f(p)
  PURE SUBROUTINE f(p,fx,plist)
    IMPLICIT NONE
    !----------------------------------------------------------------------!
    REAL, INTENT(IN)                         :: p
    REAL, INTENT(OUT)                        :: fx
    REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
    !----------------------------------------------------------------------!
    fx = f_x(p, rho_l, u_l, p_l, a_l, b_l, c_l) &
       + f_x(p, rho_r, u_r, p_r, a_r, b_r, c_r) &
       + (u_r - u_l)
  END SUBROUTINE f

  !left side or right side function
  REAL PURE FUNCTION f_x(p, rho_x, u_x, p_x, A_x, B_x, c_x)
    IMPLICIT NONE
    !----------------------------------------------------------------------!
    REAL, INTENT(IN) :: p, rho_x, u_x, p_x, a_x, b_x, c_x
    !----------------------------------------------------------------------!
    IF (p > p_x) THEN
      f_x = (p-p_x) * sqrt(a_x/(p+b_x))
    ELSE
      f_x = (2.0*c_x/(gamma-1.0)) * ((p/p_x)**((gamma-1.0)/(2.0*gamma)) - 1.0)
    END IF
  END FUNCTION f_x

SUBROUTINE funcd2d(y,fy,dfy)
  IMPLICIT NONE
  !------------------------------------------------------------------------!
  REAL, INTENT(IN)  :: y
  REAL, INTENT(OUT) :: fy,dfy
  !------------------------------------------------------------------------!
  REAL :: Ay,dAy,By,dBy,Cy,dCy
  !------------------------------------------------------------------------!
  
  ay = 2./(gam_p1*y)
  day= -ay/y
  by = gam_p1/2. * (2.-gamma*y)
  dby= -gam_p1/2. * gamma
  cy = gam_p1/gam_m1*(gamma*y-1.0)
  dcy= gamma*gam_p1/gam_m1
  fy = ay**2 * by**n1 * cy**n2 - xi**4
  dfy= day*by*cy + ay*dby*cy + ay*by*dcy
!  PRINT "(A,4(ES14.6))", "xi,n1,n2       = ", xi,n1,n2
!  PRINT "(A,3(ES14.6))", "A(y),B(y),C(y) = ", Ay,By,Cy
!  PRINT "(A,3(ES14.6))", "y,f(y),df(y)   = ", y,fy,dfy
!  PRINT "(A,3(ES14.6))", "gamma,gam_p1,gam_m1   = ", gamma,gam_p1,gam_m1
END SUBROUTINE funcd2d


SUBROUTINE funcd3d(y,fy,dfy)
  IMPLICIT NONE
  !------------------------------------------------------------------------!
  REAL, INTENT(IN)  :: y
  REAL, INTENT(OUT) :: fy,dfy
  !------------------------------------------------------------------------!
  REAL :: Ay,dAy,By,dBy,Cy,dCy
  !------------------------------------------------------------------------!

    ay = 2./(gam_p1*y)
    day= -ay/y
    by = abs(gam_p1/(7.-gamma)*(5.-(3*gamma-1.)*y))
    dby= -gam_p1/(7.-gamma)*(3*gamma-1.)
    cy = gam_p1/gam_m1*(gamma*y-1.)
    dcy= gamma*gam_p1/gam_m1
    IF (abs(gamma*y-1.).LT.1.0d-20) THEN
       fy = -xi**5
    ELSE
       fy = ay**2 * by**n1 * cy**n2 - xi**5
    END IF
    fy = ay**2 * by**n1 * cy**n2 - xi**5
    dfy= day*by*cy + ay*dby*cy + ay*by*dcy
!    PRINT "(A,ES14.6)","xi = ", xi
!    PRINT "(A,3(ES14.6))", "A(y),B(y),C(y) = ", Ay,By,Cy
!    PRINT "(A,3(ES14.6))", "y,f(y),df(y)   = ", y,fy,dfy
END SUBROUTINE funcd3d


END MODULE