mod_source_srhd.f90 Source File

!> Module for handling split source terms (split from the fluxes)


module mod_source_srhd

  use mod_variables
  use mod_physics_vars
  use mod_global_parameters, only: ndim
  
  implicit none
  public

  !> How to apply dimensional splitting to the source terms, see
  !> @ref discretization.md
  !> defaulting to sfs
  integer :: sourcesplit = 0
  
  integer, parameter :: sourcesplit_sfs    = 0
  integer, parameter :: sourcesplit_sf     = 1
  integer, parameter :: sourcesplit_ssf    = 2
  integer, parameter :: sourcesplit_ssfss  = 3

  public :: add_split_source

contains

! instantiate the templated functions here for inlining:
subroutine addsource_local(qdt, dtfactor, qtC, wCT, wCTprim, qt, wnew, x, dr, &
    qsourcesplit)
  use mod_usr, only: addsource_usr
  

  real(dp), intent(in)     :: qdt, dtfactor, qtC, qt
  real(dp), intent(in)     :: wCT(nw_phys), wCTprim(nw_phys)
  real(dp), intent(in)     :: x(1:ndim), dr(ndim)
  real(dp), intent(inout)  :: wnew(nw_phys)
  logical, intent(in)      :: qsourcesplit

  if (.not. qsourcesplit) then 
     !---------------------------------
     ! unsplit sources
     !---------------------------------

     call addsource_usr(qdt, qt, wCT, wCTprim, wnew, x, .false.)

  !!!else
     !---------------------------------
     ! split sources     
     !---------------------------------
        
  end if

end subroutine addsource_local
subroutine addsource_nonlocal(qdt, dtfactor, qtC, wCTprim, qt, wnew, x, dx, idir, &
     qsourcesplit)
  

  real(dp), intent(in)     :: qdt, dtfactor, qtC, qt
  real(dp), intent(in)     :: wCTprim(nw_phys,5)
  real(dp), intent(in)     :: x(1:ndim), dx(1:ndim)
  real(dp), intent(inout)  :: wnew(nw_phys)
  integer, intent(in)      :: idir
  logical, intent(in)      :: qsourcesplit


end subroutine addsource_nonlocal
subroutine addsource_compact(qdt, dtfactor, qtC, wCTprim1, wCTprim2, wCTprim3, qt, wnew, x, dx, &
     qsourcesplit)
  

  real(dp), intent(in)     :: qdt, dtfactor, qtC, qt
  real(dp), intent(in)     :: wCTprim1(nw_phys,3),wCTprim2(nw_phys,3),wCTprim3(nw_phys,3)
  real(dp), intent(in)     :: x(1:ndim), dx(1:ndim)
  real(dp), intent(inout)  :: wnew(nw_phys)
  logical, intent(in)      :: qsourcesplit


end subroutine addsource_compact
  pure subroutine to_primitive(u)
    use mod_con2prim
    
    real(dp), intent(inout) :: u(nw_phys)

    real(dp) :: rho,rhoh,pth,E
    real(dp) :: ssqr

    !! begin: call srhd_get_auxiliary(u)
    ssqr=(u(iw_mom(1))**2+u(iw_mom(2))**2+u(iw_mom(3))**2)
    if(srhd_eos)then
       call con2prim_eos(u(lfac_),u(xi_),u(iw_rho),ssqr,u(iw_e))
    else
       call con2prim(u(lfac_),u(xi_),u(iw_rho),ssqr,u(iw_e))
    endif
    !! end: call srhd_get_auxiliary(u)

    rho=u(iw_rho)/u(lfac_)
    rhoh=u(xi_)/u(lfac_)**2
    !! begin: call srhd_get_pressure_eos(rho,rhoh,pth,E)
    if(srhd_eos)then
      E = (rhoh+dsqrt(rhoh**2+(srhd_gamma**2-1.0_dp)*rho**2)) &
         /(srhd_gamma+1.0_dp)
      pth = 0.5d0*gamma_1* (E-rho*(rho/E))
    else
      pth =(rhoh-rho)/gamma_to_gamma_1
    endif 
    !! end: call srhd_get_pressure_eos(rho,rhoh,pth,E)

    u(iw_rho)=rho
    u(iw_mom(1))=u(lfac_)*u(iw_mom(1))/u(xi_)
    u(iw_mom(2))=u(lfac_)*u(iw_mom(2))/u(xi_)
    u(iw_mom(3))=u(lfac_)*u(iw_mom(3))/u(xi_)
    u(iw_e)=pth

  end subroutine to_primitive


  subroutine add_split_source(prior)
    use mod_global_parameters
    use mod_ghostcells_update
    use mod_physics, only: phys_req_diagonal, phys_global_source_after
    use mod_comm_lib, only: mpistop

    logical, intent(in)    :: prior
    real(dp)               :: qt
    real(dp)               :: dr(ndim), xloc(ndim)
    ! Physical cell size handed to the source term. dr is the *logical*
    ! spacing, which is not a length in a curvilinear system and, under
    ! LOG_RADIUS, is not even proportional to one.
    real(dp)               :: dloc(ndim)
    integer                :: iigrid, n, ix1, ix2, ix3
    integer                :: ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3
    real(dp)               :: wprim(nw_phys), wnew(nw_phys), wCT(nw_phys)
    logical                :: src_active
    
    src_active = .false.
    
    if ((.not.prior).and.(sourcesplit==sourcesplit_sf .or. &
       sourcesplit==sourcesplit_ssf)) return

    if (prior) then
       qt=global_time
    else
       qt=global_time+dt
    end if

    if(any_source_split) then
       ! add split source terms
       
       ixOmin1=ixGlo1+nghostcells;ixOmin2=ixGlo2+nghostcells
       ixOmin3=ixGlo3+nghostcells;ixOmax1=ixGhi1-nghostcells
       ixOmax2=ixGhi2-nghostcells;ixOmax3=ixGhi3-nghostcells;

      select case (sourcesplit)
      case (sourcesplit_sfs)
          
         do iigrid=1,igridstail_active
            n = igrids_active(iigrid)            
            dr  = rnode(rpdx1_:rnodehi, n)
            
            
            do ix3=ixOmin3,ixOmax3 
               do ix2=ixOmin2,ixOmax2 
                  do ix1=ixOmin1,ixOmax1
                     
                     xloc(1:ndim) = bgeo%x(ix1, ix2, ix3, 1:ndim, n)
                     dloc = dr
                     wCT   = bg(1)%w(ix1,ix2,ix3, 1:nw_phys, n)
                     wnew  = wCT
                     wprim = wCT
                     call to_primitive( wprim )

                     call addsource_local(0.5d0*dt,&
                          0.5d0, qt, wCT,&
                          wprim, qt, wnew, xloc, dloc, .true. )
                     bg(1)%w(ix1,ix2,ix3, 1:nw_phys, n) = wnew(1:nw_phys)

                     
                end do
             end do
          end do
        end do
        src_active = .true.
      case default
         write(unitterm,*)'No such type/not yet implemented sourcesplit=',sourcesplit
         call mpistop("Error: Unknown type of sourcesplit!")
      end select

      if (.not. prior .and. associated(phys_global_source_after)) then
         call phys_global_source_after(dt, qt, src_active)
      end if

      if (src_active) then
         call getbc(qt,0.d0,ps,iwstart,nwgc,phys_req_diagonal)
      end if
      
   end if 

  end subroutine add_split_source

end module mod_source_srhd