mod_rho_phys.t

Go to the documentation of this file.

!> Module containing the physics routines for scalar advection
module mod_rho_phys
 
  implicit none
  private
 
  integer, protected, public          :: rho_       = 1
 
  !> Whether particles module is added
  logical, public, protected              :: rho_particles = .false.
 
  double precision, protected, public :: rho_v(^nd) = 1.0d0
 
  ! Public methods
  public :: rho_phys_init
  public :: rho_get_v
 
contains
 
  subroutine rho_params_read(files)
    use mod_global_parameters, only: unitpar
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /rho_list/ rho_v, rho_particles
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status='old')
       read(unitpar, rho_list, end=111)
111    close(unitpar)
    end do
 
  end subroutine rho_params_read
 
  !> Write this module's parameters to a snapsoht
  subroutine rho_write_info(fh)
    use mod_global_parameters
    integer, intent(in)                 :: fh
    integer, parameter                  :: n_par = ^nd
    double precision                    :: values(n_par)
    character(len=name_len)             :: names(n_par)
    integer, dimension(MPI_STATUS_SIZE) :: st
    integer                             :: er
    integer                             :: idim
 
    call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
 
    do idim=1,ndim
      write(names(idim),'(a,i1)') "v",idim
      values(idim) = rho_v(idim)
    end do
    call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
    call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
  end subroutine rho_write_info
 
  subroutine rho_phys_init()
    use mod_global_parameters
    use mod_physics
    use mod_particles, only: particles_init
 
    call rho_params_read(par_files)
 
    physics_type = "rho"
    phys_energy  = .false.
    phys_req_diagonal = .false.
    use_particles = rho_particles
 
    allocate(start_indices(number_species),stop_indices(number_species))
    ! set the index of the first flux variable for species 1
    start_indices(1)=1
    rho_ = var_set_rho()
 
    ! set number of variables which need update ghostcells
    nwgc=nwflux
 
    ! set the index of the last flux variable for species 1
    stop_indices(1)=nwflux
 
    ! Check whether custom flux types have been defined
    if (.not. allocated(flux_type)) then
       allocate(flux_type(ndir, nw))
       flux_type = flux_default
    else if (any(shape(flux_type) /= [ndir, nw])) then
       call mpistop("phys_check error: flux_type has wrong shape")
    end if
 
    phys_get_cmax        => rho_get_cmax
    phys_get_cbounds     => rho_get_cbounds
    phys_get_flux        => rho_get_flux
    phys_add_source_geom => rho_add_source_geom
    phys_to_conserved    => rho_to_conserved
    phys_to_primitive    => rho_to_primitive
    phys_get_dt          => rho_get_dt
    phys_write_info      => rho_write_info
 
    ! Initialize particles module
    if (rho_particles) then
       call particles_init()
       phys_req_diagonal = .true.
    end if
 
  end subroutine rho_phys_init
 
  subroutine rho_to_conserved(ixI^L, ixO^L, w, x)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(inout) :: w(ixI^S, nw)
    double precision, intent(in)    :: x(ixI^S, 1:^ND)
 
    ! Do nothing (primitive and conservative are equal for rho module)
  end subroutine rho_to_conserved
 
  subroutine rho_to_primitive(ixI^L, ixO^L, w, x)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(inout) :: w(ixI^S, nw)
    double precision, intent(in)    :: x(ixI^S, 1:^ND)
 
    ! Do nothing (primitive and conservative are equal for rho module)
  end subroutine rho_to_primitive
 
  subroutine rho_get_v(w, x, ixI^L, ixO^L, idim, v, centered)
    use mod_global_parameters
    use mod_geometry
    logical, intent(in)           :: centered
    integer, intent(in)           :: ixi^l, ixo^l, idim
    double precision, intent(in)  :: w(ixi^s, nw), x(ixi^s, 1:^nd)
    double precision, intent(out) :: v(ixi^s)
 
    double precision :: dtheta, dphi, halfdtheta, halfdphi, invdtheta, invdphi
    {^ifthreed
    double precision :: appcosphi(ixi^s), appsinphi(ixi^s), &
                        appcosthe(ixi^s), appsinthe(ixi^s)
    }
 
    select case (coordinate)
    case (cylindrical)
       {^ifoned 
       call mpistop("advection in 1D cylindrical not available")
       }
       {^iftwod 
       ! advection in 2D cylindrical: polar grid: to v_R, v_varphi
       if(centered)then
          select case (idim)
             case (1) ! radial velocity
                v(ixo^s) = rho_v(1)*dcos(x(ixo^s,2))+rho_v(2)*dsin(x(ixo^s,2))
             case (2) ! v_varphi
                v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,2))+rho_v(2)*dcos(x(ixo^s,2))
          end select
       else
          ! assumed uniform in varphi=theta
          dtheta=x(ixomin1,ixomin2+1,2)-x(ixomin1,ixomin2,2)
          halfdtheta=0.5d0*dtheta
          invdtheta=1.0d0/dtheta
          select case (idim)
             case (1) ! radial velocity
                v(ixo^s) =( rho_v(1)*( dsin(x(ixo^s,2)+halfdtheta) &
                                      -dsin(x(ixo^s,2)-halfdtheta)) &
                           +rho_v(2)*(-dcos(x(ixo^s,2)+halfdtheta) &
                                      +dcos(x(ixo^s,2)-halfdtheta)))*invdtheta
             case (2) ! v_varphi
                v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,2)+halfdtheta) &
                          +rho_v(2)*dcos(x(ixo^s,2)+halfdtheta)
          end select
       endif
       }
       {^ifthreed
       ! advection in 3D cylindrical: convert to v_R, v_Z, v_varphi
       if(centered)then
          select case (idim)
             case (1) ! v_R velocity
                v(ixo^s) = rho_v(1)*dcos(x(ixo^s,3))+rho_v(2)*dsin(x(ixo^s,3))
             case (2) ! v_Z velocity
                v(ixo^s) = rho_v(3)
             case (3) ! v_varphi velocity
                v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,3))+rho_v(2)*dcos(x(ixo^s,3))
          end select
       else
          ! assumed uniform in varphi=theta
          dtheta=x(ixomin1,ixomin2,ixomin3+1,3)-x(ixomin1,ixomin2,ixomin3,3)
          halfdtheta=0.5d0*dtheta
          invdtheta=1.0d0/dtheta
          select case (idim)
             case (1) ! radial velocity
                v(ixo^s) =( rho_v(1)*( dsin(x(ixo^s,3)+halfdtheta) &
                                      -dsin(x(ixo^s,3)-halfdtheta)) &
                           +rho_v(2)*(-dcos(x(ixo^s,3)+halfdtheta) &
                                      +dcos(x(ixo^s,3)-halfdtheta)))*invdtheta
             case (2) ! v_Z velocity
                v(ixo^s) = rho_v(3)
             case (3) ! v_varphi
                v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,3)+halfdtheta) &
                          +rho_v(2)*dcos(x(ixo^s,3)+halfdtheta)
          end select
       endif
       }
    case (spherical)
       {^ifoned 
       call mpistop("advection in 1D spherical not available")
       }
       {^iftwod 
       call mpistop("advection in 2D spherical not available")
       }
       {^ifthreed
       ! advection in 3D spherical: convert to v_r, v_theta, v_phi
       if(centered)then
          select case (idim)
             case (1) ! radial velocity
                v(ixo^s) = rho_v(1)*dsin(x(ixo^s,2))*dcos(x(ixo^s,3)) &
                          +rho_v(2)*dsin(x(ixo^s,2))*dsin(x(ixo^s,3)) &
                          +rho_v(3)*dcos(x(ixo^s,2)) 
             case (2) ! theta velocity
                v(ixo^s) = rho_v(1)*dcos(x(ixo^s,2))*dcos(x(ixo^s,3)) &
                          +rho_v(2)*dcos(x(ixo^s,2))*dsin(x(ixo^s,3)) &
                          -rho_v(3)*dsin(x(ixo^s,2)) 
             case (3) ! phi velocity
                v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,3)) &
                          +rho_v(2)*dcos(x(ixo^s,3)) 
          end select
       else
          ! assumed uniform in theta and phi
          dtheta=x(ixomin1,ixomin2+1,ixomin3,2)-x(ixomin1,ixomin2,ixomin3,2)
          dphi=x(ixomin1,ixomin2,ixomin3+1,3)-x(ixomin1,ixomin2,ixomin3,3)
          halfdtheta=0.5d0*dtheta
          halfdphi=0.5d0*dphi
          invdtheta=1.0d0/dtheta
          invdphi=1.0d0/dphi
          select case (idim)
             case (1) ! radial velocity
                appcosphi(ixo^s)=( dsin(x(ixo^s,3)+halfdphi) &
                                  -dsin(x(ixo^s,3)-halfdphi))*invdphi 
                appsinphi(ixo^s)=(-dcos(x(ixo^s,3)+halfdphi) &
                                  +dcos(x(ixo^s,3)-halfdphi))*invdphi 
                appcosthe(ixo^s)=(dsin(x(ixo^s,2)+halfdtheta)**2 &
                                 -dsin(x(ixo^s,2)-halfdtheta)**2) &
                          /(4.0d0*dabs(dsin(x(ixo^s,2)))*dsin(halfdtheta))
                appsinthe(ixo^s)= &
                           (-dsin(x(ixo^s,2)+halfdtheta)*dcos(x(ixo^s,2)+halfdtheta)  & 
                            +dsin(x(ixo^s,2)-halfdtheta)*dcos(x(ixo^s,2)-halfdtheta)  & 
                            +dtheta)/(4.0d0*dabs(dsin(x(ixo^s,2)))*dsin(halfdtheta))
                v(ixo^s) = rho_v(1)*appsinthe(ixo^s)*appcosphi(ixo^s) &
                          +rho_v(2)*appsinthe(ixo^s)*appsinphi(ixo^s) &
                          +rho_v(3)*appcosthe(ixo^s) 
             case (2) ! theta velocity
                appcosphi(ixo^s)=( dsin(x(ixo^s,3)+halfdphi) &
                                  -dsin(x(ixo^s,3)-halfdphi))*invdphi 
                appsinphi(ixo^s)=(-dcos(x(ixo^s,3)+halfdphi) &
                                  +dcos(x(ixo^s,3)-halfdphi))*invdphi 
                v(ixo^s) = rho_v(1)*dcos(x(ixo^s,2)+halfdtheta)*appcosphi(ixo^s) &
                          +rho_v(2)*dcos(x(ixo^s,2)+halfdtheta)*appsinphi(ixo^s) &
                          -rho_v(3)*dsin(x(ixo^s,2)+halfdtheta) 
             case (3) ! phi velocity
                v(ixo^s) =-rho_v(1)*dsin(x(ixo^s,3)+halfdphi) &
                          +rho_v(2)*dcos(x(ixo^s,3)+halfdphi) 
          end select
       endif
       }
    case default
       v(ixo^s) = rho_v(idim)
    end select
  end subroutine rho_get_v
 
  !> Calculate simple v component
  subroutine rho_get_v_idim(w,x,ixI^L,ixO^L,idim,v)
    use mod_global_parameters
 
    integer, intent(in)           :: ixI^L, ixO^L, idim
    double precision, intent(in)  :: w(ixI^S,nw), x(ixI^S,1:ndim)
    double precision, intent(out) :: v(ixI^S)
 
    v(ixo^s) = rho_v(idim)
 
  end subroutine rho_get_v_idim
 
  subroutine rho_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
    use mod_global_parameters
    integer, intent(in)                       :: ixI^L, ixO^L, idim
    double precision, intent(in)              :: w(ixI^S, nw), x(ixI^S, 1:^ND)
    double precision, intent(inout)           :: cmax(ixI^S)
 
    call rho_get_v(w, x, ixi^l, ixo^l, idim, cmax, .true.)
 
    cmax(ixo^s) = abs(cmax(ixo^s))
 
  end subroutine rho_get_cmax
 
  subroutine rho_get_cbounds(wLC, wRC, wLp, wRp, x, ixI^L, ixO^L, idim,Hspeed, cmax, cmin)
    use mod_global_parameters
    use mod_variables
    integer, intent(in)             :: ixI^L, ixO^L, idim
    double precision, intent(in)    :: wLC(ixI^S, nw), wRC(ixI^S,nw)
    double precision, intent(in)    :: wLp(ixI^S, nw), wRp(ixI^S,nw)
    double precision, intent(in)    :: x(ixI^S, 1:^ND)
    double precision, intent(inout) :: cmax(ixI^S,1:number_species)
    double precision, intent(inout), optional :: cmin(ixI^S,1:number_species)
    double precision, intent(in)    :: Hspeed(ixI^S,1:number_species)
 
    ! If get_v depends on w, the first argument should be some average over the
    ! left and right state
    call rho_get_v(wlc, x, ixi^l, ixo^l, idim, cmax(ixi^s,1), .false.)
 
    if (present(cmin)) then
       cmin(ixo^s,1) = min(cmax(ixo^s,1), zero)
       cmax(ixo^s,1) = max(cmax(ixo^s,1), zero)
    else
       cmax(ixo^s,1) = maxval(abs(cmax(ixo^s,1)))
    end if
 
  end subroutine rho_get_cbounds
 
  subroutine rho_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: dx^D, x(ixI^S, 1:^ND)
    double precision, intent(in)    :: w(ixI^S, 1:nw)
    double precision, intent(inout) :: dtnew
 
    dtnew = bigdouble
  end subroutine rho_get_dt
 
  ! There is nothing to add to the transport flux in the transport equation
  subroutine rho_get_flux(wC, w, x, ixI^L, ixO^L, idim, f)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L, idim
    double precision, intent(in)    :: wC(ixI^S, 1:nw)
    double precision, intent(in)    :: w(ixI^S, 1:nw)
    double precision, intent(in)    :: x(ixI^S, 1:^ND)
    double precision, intent(out)   :: f(ixI^S, nwflux)
    double precision                :: v(ixI^S)
 
    call rho_get_v(wc, x, ixi^l, ixo^l, idim, v, .false.)
 
    f(ixo^s, rho_) = w(ixo^s, rho_) * v(ixo^s)
  end subroutine rho_get_flux
 
  subroutine rho_add_source_geom(qdt, dtfactor, ixI^L, ixO^L, wCT,wprim, w, x)
 
    ! Add geometrical source terms to w
    ! There are no geometrical source terms in the transport equation
 
    use mod_global_parameters
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: qdt, dtfactor, x(ixI^S, 1:^ND)
    double precision, intent(inout) :: wCT(ixI^S, 1:nw),wprim(ixI^S,1:nw), w(ixI^S, 1:nw)
 
  end subroutine rho_add_source_geom
 
end module mod_rho_phys