mod_hd_ppm.t

Go to the documentation of this file.

!> Hydrodynamics PPM module
module mod_hd_ppm
  use mod_hd_phys
 
  implicit none
  private
 
  public :: hd_ppm_init
 
contains
 
  subroutine hd_ppm_init()
    use mod_physics_ppm
 
    phys_ppm_flatcd => hd_ppm_flatcd
    phys_ppm_flatsh => hd_ppm_flatsh
  end subroutine hd_ppm_init
 
  subroutine hd_ppm_flatcd(ixI^L,ixO^L,ixL^L,ixR^L,w,d2w,drho,dp)
    use mod_global_parameters
 
    integer, intent(in)             :: ixI^L, ixO^L, ixL^L, ixR^L
    double precision, intent(in)    :: w(ixI^S, nw), d2w(ixG^T, 1:nwflux)
    double precision, intent(inout) :: drho(ixG^T), dp(ixG^T)
 
    if(hd_energy) then
      drho(ixo^s) = hd_gamma*abs(d2w(ixo^s, rho_))&
           /min(w(ixl^s, rho_), w(ixr^s, rho_))
      dp(ixo^s) = abs(d2w(ixo^s, e_))/min(w(ixl^s, e_), w(ixr^s, e_))
    else
      call mpistop("PPM with flatcd=.true. can not be used with -eos = iso !")
    end if
  end subroutine hd_ppm_flatcd
 
  subroutine hd_ppm_flatsh(ixI^L,ixO^L,ixLL^L,ixL^L,ixR^L,ixRR^L,idims,w,drho,dp,dv)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)             :: ixI^L, ixO^L, ixLL^L, ixL^L, ixR^L, ixRR^L
    integer, intent(in)             :: idims
    double precision, intent(in)    :: w(ixI^S, nw)
    double precision, intent(inout) :: drho(ixG^T), dp(ixG^T), dv(ixG^T)
    double precision                :: v(ixG^T)
 
    if(hd_energy) then
      ! eq. B15, page 218, Mignone and Bodo 2005, ApJS (beta1)
      where (abs(w(ixrr^s, e_)-w(ixll^s, e_))>smalldouble)
         drho(ixo^s) = abs((w(ixr^s, e_)-w(ixl^s, e_))&
              /(w(ixrr^s, e_)-w(ixll^s, e_)))
      else where
         drho(ixo^s) = zero
      end where
 
      !  eq. B76, page 48, Miller and Collela 2002, JCP 183, 26
      !  use "dp" to save squared sound speed, assuming primitives
      dp(ixo^s) =(hd_gamma*w(ixo^s, e_)/w(ixo^s, rho_))
 
      dp(ixo^s) = abs(w(ixr^s, e_)-w(ixl^s, e_))&
           /(w(ixo^s, rho_)*dp(ixo^s))
      v(ixi^s)  = w(ixi^s, mom(idims))
      call gradient(v, ixi^l, ixo^l, idims, dv)
    else
      call mpistop("PPM with flatsh=.true. can not be used with -eos = iso !")
    end if
  end subroutine hd_ppm_flatsh
 
end module mod_hd_ppm