Merge pull request #608 from danieljprice/gr_sink

.github/workflows/test.yml

@@ -31,7 +31,7 @@ jobs:
           - ['test', '']
           - ['testkd', '']
           - ['testdust', 'dust']
-          - ['testgr', 'gr']
+          - ['testgr', 'gr ptmass']
           - ['testgrav', 'gravity ptmass setstar']
           - ['testgrowth', 'dustgrowth']
           - ['testnimhd', 'nimhd']

AUTHORS

@@ -13,8 +13,8 @@ David Liptai <[email protected]>
 Lionel Siess <[email protected]>
 Fangyi (Fitz) Hu <[email protected]>
 Yann Bernard <[email protected]>
-Daniel Mentiplay <[email protected]>
 Megha Sharma <[email protected]>
+Daniel Mentiplay <[email protected]>
 Arnaud Vericel <[email protected]>
 Mark Hutchison <[email protected]>
 Mats Esseldeurs <[email protected]>
@@ -60,6 +60,7 @@ Benoit Commercon <[email protected]>
 Christopher Russell <[email protected]>
 Giulia Ballabio <[email protected]>
 Joe Fisher <[email protected]>
+Kateryna Andrych <[email protected]>
 Maxime Lombart <[email protected]>
 Orsola De Marco <[email protected]>
 Shunquan Huang <[email protected]>

src/main/checksetup.f90

@@ -410,6 +410,10 @@ subroutine check_setup(nerror,nwarn,restart)
 !
  if (gr) call check_gr(npart,nerror,xyzh,vxyzu)
 !
+!--check sink GR setup
+!
+ if (gr) call check_gr(nptmass,nerror,xyzmh_ptmass,vxyz_ptmass)
+!
 !--check radiation setup
 !
  if (do_radiation) call check_setup_radiation(npart,nerror,nwarn,radprop,rad)
@@ -544,12 +548,6 @@ subroutine check_setup_ptmass(nerror,nwarn,hmin)
 
  isoblate = .false.
 
- if (gr .and. nptmass > 0) then
-    print*,' ERROR: nptmass = ',nptmass, ' should be = 0 for GR'
-    nwarn = nwarn + 1
-    return
- endif
-
  if (nptmass < 0) then
     print*,' ERROR: nptmass = ',nptmass, ' should be >= 0 '
     nerror = nerror + 1

src/main/config.F90

@@ -262,8 +262,10 @@ module dim
  integer :: maxgr = 0
 #ifdef GR
  logical, parameter :: gr = .true.
+ integer, parameter :: maxptmassgr = maxptmass
 #else
  logical, parameter :: gr = .false.
+ integer, parameter :: maxptmassgr = 0
 #endif
 
 !---------------------

src/main/cons2prim.f90

@@ -20,7 +20,7 @@ module cons2prim
 !   Liptai & Price (2019), MNRAS 485, 819-842
 !   Ballabio et al. (2018), MNRAS 477, 2766-2771
 !
-! :Owner: Elisabeth Borchert
+! :Owner: Megha Sharma
 !
 ! :Runtime parameters: None
 !
@@ -29,7 +29,7 @@ module cons2prim
 !
  implicit none
 
- public :: cons2primall,cons2prim_everything
+ public :: cons2primall,cons2prim_everything,cons2primall_sink
  public :: prim2consall,prim2consi
 
  private
@@ -44,17 +44,17 @@ module cons2prim
 !  (density,velocity,internal energy), for ALL particles
 !+
 !----------------------------------------------------------------------
-subroutine prim2consall(npart,xyzh,metrics,vxyzu,dens,pxyzu,use_dens)
+subroutine prim2consall(npart,xyzh,metrics,vxyzu,pxyzu,use_dens,dens,use_sink)
  use part, only:isdead_or_accreted,ien_type,eos_vars,igasP,igamma,itemp
  use eos,  only:gamma,ieos
  integer, intent(in)  :: npart
  real,    intent(in)  :: xyzh(:,:),metrics(:,:,:,:),vxyzu(:,:)
- real,    intent(inout) :: dens(:)
  real,    intent(out) :: pxyzu(:,:)
- logical, intent(in), optional :: use_dens
+ real,    intent(inout), optional :: dens(:)
+ logical, intent(in),    optional :: use_dens, use_sink
  logical :: usedens
  integer :: i
- real    :: pri,tempi
+ real    :: pri,tempi,xyzhi(4),vxyzui(4),densi
 
 !  By default, use the smoothing length to compute primitive density, and then compute the conserved variables.
 !  (Alternatively, use the provided primitive density to compute conserved variables.
@@ -65,25 +65,37 @@ subroutine prim2consall(npart,xyzh,metrics,vxyzu,dens,pxyzu,use_dens)
     usedens = .false.
  endif
 
-!$omp parallel do default (none) &
-!$omp shared(xyzh,metrics,vxyzu,dens,pxyzu,npart,usedens,ien_type,eos_vars,gamma,ieos) &
-!$omp private(i,pri,tempi)
+ !$omp parallel do default (none) &
+ !$omp shared(xyzh,metrics,vxyzu,dens,pxyzu,npart,usedens,ien_type,eos_vars,gamma,ieos,use_sink,use_dens) &
+ !$omp private(i,pri,tempi,xyzhi,vxyzui,densi)
  do i=1,npart
-    if (.not.isdead_or_accreted(xyzh(4,i))) then
-       call prim2consi(xyzh(:,i),metrics(:,:,:,i),vxyzu(:,i),dens(i),pri,tempi,pxyzu(:,i),usedens,ien_type)
 
-       ! save eos vars for later use
-       eos_vars(igasP,i)  = pri
-       eos_vars(itemp,i)  = tempi
-       if (vxyzu(4,i) > 0. .and. ieos == 12) then
-          eos_vars(igamma,i) = 1. + pri/(dens(i)*vxyzu(4,i))
-       else
-          ! prevent getting NaN or Infinity when u = 0
-          eos_vars(igamma,i) = gamma
+    if (present(use_sink)) then
+       xyzhi(1:3) = xyzh(1:3,i) ! save positions
+       xyzhi(4) = xyzh(5,i) ! save smoothing length, h
+       vxyzui(1:3) = vxyzu(1:3,i)
+       vxyzui(4) = 0. ! assume energy as 0. for sink
+       densi = 1.
+       call prim2consi(xyzhi,metrics(:,:,:,i),vxyzui,pri,tempi,pxyzu(:,i),ien_type,&
+                   use_sink=use_sink,dens_i=densi) ! this returns temperature and pressure as 0.
+    else
+       if (.not.isdead_or_accreted(xyzh(4,i))) then
+          call prim2consi(xyzh(:,i),metrics(:,:,:,i),vxyzu(:,i),pri,tempi,pxyzu(:,i),ien_type,&
+                   use_dens=usedens,dens_i=dens(i))
+
+          ! save eos vars for later use
+          eos_vars(igasP,i)  = pri
+          eos_vars(itemp,i)  = tempi
+          if (vxyzu(4,i) > 0. .and. ieos == 12) then
+             eos_vars(igamma,i) = 1. + pri/(dens(i)*vxyzu(4,i))
+          else
+             ! prevent getting NaN or Infinity when u = 0
+             eos_vars(igamma,i) = gamma
+          endif
        endif
     endif
  enddo
-!$omp end parallel do
+ !$omp end parallel do
 
 end subroutine prim2consall
 
@@ -93,19 +105,21 @@ end subroutine prim2consall
 !  for a single SPH particle
 !+
 !----------------------------------------------------------------------
-subroutine prim2consi(xyzhi,metrici,vxyzui,dens_i,pri,tempi,pxyzui,use_dens,ien_type)
+subroutine prim2consi(xyzhi,metrici,vxyzui,pri,tempi,pxyzui,ien_type,use_dens,use_sink,dens_i)
  use cons2primsolver, only:primitive2conservative
  use utils_gr,        only:h2dens
  use eos,             only:equationofstate,ieos
  real, dimension(4), intent(in)  :: xyzhi, vxyzui
  real,               intent(in)  :: metrici(:,:,:)
- real, intent(inout)             :: dens_i,pri,tempi
+ real, intent(inout)             :: pri,tempi
  integer,            intent(in)  :: ien_type
  real, dimension(4), intent(out) :: pxyzui
- logical, intent(in), optional   :: use_dens
+ logical, intent(in), optional   :: use_dens,use_sink
+ real, intent(inout), optional   :: dens_i
  logical :: usedens
  real    :: rhoi,ui,xyzi(1:3),vi(1:3),pondensi,spsoundi,densi
 
+ pondensi = 0.
  !  By default, use the smoothing length to compute primitive density, and then compute the conserved variables.
  !  (Alternatively, use the provided primitive density to compute conserved variables.
  !   Depends whether you have prim dens prior or not.)
@@ -118,13 +132,20 @@ subroutine prim2consi(xyzhi,metrici,vxyzui,dens_i,pri,tempi,pxyzui,use_dens,ien_
  xyzi = xyzhi(1:3)
  vi   = vxyzui(1:3)
  ui   = vxyzui(4)
+
  if (usedens) then
     densi = dens_i
  else
-    call h2dens(densi,xyzhi,metrici,vi) ! Compute dens from h
-    dens_i = densi                      ! Feed the newly computed dens back out of the routine
+    if (present(use_sink)) then
+       densi    = 1.    ! using a value of 0. results in NaN values for the pxyzui array.
+       pondensi = 0.
+    else
+       call h2dens(densi,xyzhi,metrici,vi) ! Compute dens from h
+       dens_i = densi ! Feed the newly computed dens back out of the routine
+       call equationofstate(ieos,pondensi,spsoundi,densi,xyzi(1),xyzi(2),xyzi(3),tempi,ui)
+    endif
  endif
- call equationofstate(ieos,pondensi,spsoundi,densi,xyzi(1),xyzi(2),xyzi(3),tempi,ui)
+
  pri = pondensi*densi
  call primitive2conservative(xyzi,metrici,vi,densi,ui,pri,rhoi,pxyzui(1:3),pxyzui(4),ien_type)
 
@@ -173,7 +194,6 @@ subroutine cons2primall(npart,xyzh,metrics,pxyzu,vxyzu,dens,eos_vars)
 
        call conservative2primitive(xyzh(1:3,i),metrics(:,:,:,i),vxyzu(1:3,i),dens(i),vxyzu(4,i), &
                                   p_guess,tempi,gammai,rhoi,pxyzu(1:3,i),pxyzu(4,i),ierr,ien_type)
-
        ! store results
        eos_vars(igasP,i)  = p_guess
        eos_vars(ics,i)    = get_spsound(ieos,xyzh(1:3,i),dens(i),vxyzu(:,i),gammai)
@@ -191,6 +211,49 @@ subroutine cons2primall(npart,xyzh,metrics,pxyzu,vxyzu,dens,eos_vars)
 
 end subroutine cons2primall
 
+!----------------------------------------------------------------------
+!+
+!  Conservative to primitive routines (for GR sink particles):
+!  Solve for primitive variables (density,velocity,internal energy)
+!  from the evolved/conservative variables (rho*,momentum,entropy)
+!+
+!----------------------------------------------------------------------
+subroutine cons2primall_sink(nptmass,xyzmh_ptmass,metrics_ptmass,pxyzu_ptmass,vxyz_ptmass,eos_vars)
+ use cons2primsolver, only:conservative2primitive
+ use io,              only:fatal
+ use part,            only:ien_type
+ integer, intent(in)    :: nptmass
+ real,    intent(in)    :: pxyzu_ptmass(:,:),xyzmh_ptmass(:,:),metrics_ptmass(:,:,:,:)
+ real,    intent(inout) :: vxyz_ptmass(:,:)
+ real,    intent(out), optional   :: eos_vars(:,:)
+ integer :: i, ierr
+ real    :: p_guess,rhoi,tempi,gammai,eni,densi
+
+!$omp parallel do default (none) &
+!$omp shared(xyzmh_ptmass,metrics_ptmass,vxyz_ptmass,pxyzu_ptmass,nptmass,ien_type) &
+!$omp private(i,ierr,p_guess,rhoi,tempi,gammai,eni,densi)
+ do i=1,nptmass
+    p_guess = 0.
+    tempi   = 0.
+    gammai  = 0.
+    rhoi    = 1.
+    densi   = 1.
+    ! conservative 2 primitive
+    call conservative2primitive(xyzmh_ptmass(1:3,i),metrics_ptmass(:,:,:,i),vxyz_ptmass(1:3,i),densi,eni, &
+                              p_guess,tempi,gammai,rhoi,pxyzu_ptmass(1:3,i),pxyzu_ptmass(4,i),ierr,ien_type)
+
+    if (ierr > 0) then
+       print*,' pmom =',pxyzu_ptmass(1:3,i)
+       print*,' rho* =',rhoi
+       print*,' en   =',eni
+       call fatal('cons2prim','could not solve rootfinding',i)
+    endif
+
+ enddo
+!$omp end parallel do
+
+end subroutine cons2primall_sink
+
 !-----------------------------------------------------------------------------
 !+
 !  Solve for primitive variables (v,u,P,B,dustfrac) from evolved variables

src/main/deriv.F90

@@ -221,6 +221,7 @@ end subroutine derivs
 !  this should NOT be called during timestepping, it is useful
 !  for when one requires just a single call to evaluate derivatives
 !  and store them in the global shared arrays
+!  does not work for sink GR yet
 !+
 !--------------------------------------
 subroutine get_derivs_global(tused,dt_new,dt)
@@ -246,7 +247,7 @@ subroutine get_derivs_global(tused,dt_new,dt)
  ! update conserved quantities in the GR code
  if (gr) then
     call init_metric(npart,xyzh,metrics)
-    call prim2consall(npart,xyzh,metrics,vxyzu,dens,pxyzu,use_dens=.false.)
+    call prim2consall(npart,xyzh,metrics,vxyzu,pxyzu,use_dens=.false.,dens=dens)
  endif
 
  ! evaluate derivatives