mod_gravity.t

Go to the documentation of this file.

!> Module for including gravity in (magneto)hydrodynamics simulations
module mod_gravity
  implicit none
 
  !> source split or not
  logical :: grav_split= .false.
 
contains
  !> Read this module's parameters from a file
  subroutine grav_params_read(files)
    use mod_global_parameters, only: unitpar
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /grav_list/ grav_split
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, grav_list, end=111)
111    close(unitpar)
    end do
 
  end subroutine grav_params_read
 
  !> Initialize the module
  subroutine gravity_init()
    use mod_global_parameters
    integer :: nwx,idir
 
    call grav_params_read(par_files)
 
  end subroutine gravity_init
 
  !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
  subroutine gravity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
       energy,qsourcesplit,active)
    use mod_global_parameters
    use mod_usr_methods
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt, x(ixI^S,1:ndim)
    double precision, intent(in)    :: wCT(ixI^S,1:nw)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    logical, intent(in) :: energy,qsourcesplit
    logical, intent(inout) :: active
    integer                         :: idim
 
    double precision :: gravity_field(ixI^S,ndim)
 
    if(qsourcesplit .eqv. grav_split) then
      active = .true.
 
      if (.not. associated(usr_gravity)) then
        write(*,*) "mod_usr.t: please point usr_gravity to a subroutine"
        write(*,*) "like the phys_gravity in mod_usr_methods.t"
        call mpistop("gravity_add_source: usr_gravity not defined")
      else
        call usr_gravity(ixi^l,ixo^l,wct,x,gravity_field)
      end if
  
      do idim = 1, ndim
        w(ixo^s,iw_mom(idim)) = w(ixo^s,iw_mom(idim)) &
              + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,iw_rho)
        if(energy) then
          if(iw_equi_rho>0) then
            w(ixo^s,iw_e)=w(ixo^s,iw_e) &
              + qdt * gravity_field(ixo^s,idim) *&
               wct(ixo^s,iw_mom(idim))/(wct(ixo^s,iw_rho)+block%equi_vars(ixo^s,iw_equi_rho,0))*&
              wct(ixo^s,iw_rho)
 
          else
            w(ixo^s,iw_e)=w(ixo^s,iw_e) &
              + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,iw_mom(idim))
          endif
        end if
      end do
    end if
 
  end subroutine gravity_add_source
 
  subroutine gravity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
    use mod_global_parameters
    use mod_usr_methods
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: dx^D, x(ixI^S,1:ndim), w(ixI^S,1:nw)
    double precision, intent(inout) :: dtnew
 
    double precision                :: dxinv(1:ndim), max_grav
    integer                         :: idim
 
    double precision :: gravity_field(ixI^S,ndim)
 
    ^d&dxinv(^d)=one/dx^d;
 
    if(.not. associated(usr_gravity)) then
      write(*,*) "mod_usr.t: please point usr_gravity to a subroutine"
      write(*,*) "like the phys_gravity in mod_usr_methods.t"
      call mpistop("gravity_get_dt: usr_gravity not defined")
    else
      call usr_gravity(ixi^l,ixo^l,w,x,gravity_field)
    end if
 
    do idim = 1, ndim
      max_grav = maxval(abs(gravity_field(ixo^s,idim)))
      max_grav = max(max_grav, epsilon(1.0d0))
      dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
    end do
 
  end subroutine gravity_get_dt
 
end module mod_gravity