Merge pull request #568 from rebeccagmartin/master

added sink surface force from Ayliffe & Bate (2010)

src/main/ptmass.F90

@@ -59,6 +59,7 @@ module ptmass
 
  ! settings affecting routines in module (read from/written to input file)
  integer, public :: icreate_sinks = 0
+ integer, public :: isink_potential = 0
  real,    public :: rho_crit_cgs  = 1.e-10
  real,    public :: r_crit = 5.e-3
  real,    public :: h_acc  = 1.e-3
@@ -69,7 +70,6 @@ module ptmass
  real,    public :: r_merge_cond    = 0.0     ! sinks will merge if bound within this radius
  real,    public :: f_crit_override = 0.0     ! 1000.
 
-
  logical, public :: use_regnbody    = .false. ! subsystems switch
  logical, public :: use_fourthorder = .true.
  integer, public :: n_force_order   = 3
@@ -160,7 +160,7 @@ subroutine get_accel_sink_gas(nptmass,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi,phi,
  real                             :: ftmpxi,ftmpyi,ftmpzi
  real                             :: dx,dy,dz,rr2,ddr,dr3,f1,f2,pmassj,J2,shat(3),Rsink
  real                             :: hsoft,hsoft1,hsoft21,q2i,qi,psoft,fsoft
- real                             :: fxj,fyj,fzj,dsx,dsy,dsz
+ real                             :: fxj,fyj,fzj,dsx,dsy,dsz,fac,r
  integer                          :: j
  logical                          :: tofrom,extrap
  !
@@ -241,14 +241,30 @@ subroutine get_accel_sink_gas(nptmass,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi,phi,
 
        ! acceleration of gas due to point mass particle
        f1     = pmassj*dr3
+
+       ! acceleration of sink from gas
+       if (tofrom) f2 = pmassi*dr3
+
+       ! modified potential
+       select case (isink_potential)
+       case(1)
+          ! Ayliffe & Bate (2010) equation 2 (prevent accretion on to sink)
+          Rsink = xyzmh_ptmass(iReff,j)
+          r=1./ddr
+          if (Rsink > 0. .and. r < 2*Rsink) then
+             fac = (1. - (2. - r/Rsink)**4)
+             f1 = f1*fac
+             f2 = f2*fac
+             phi = phi - pmassj*(r**3/3.-4.*r**2*Rsink+24.*r*Rsink**2 &
+                  -16.*Rsink**4/r-32.*Rsink**3*log(r))/Rsink**4
+          endif
+       end select
+
        ftmpxi = ftmpxi - dx*f1
        ftmpyi = ftmpyi - dy*f1
        ftmpzi = ftmpzi - dz*f1
        phi    = phi    - pmassj*ddr      ! potential (GM/r)
 
-       ! acceleration of sink from gas
-       if (tofrom) f2 = pmassi*dr3
-
        ! additional accelerations due to oblateness
        if (abs(J2) > 0.) then
           shat = unitvec(xyzmh_ptmass(ispinx:ispinz,j))
@@ -1934,6 +1950,7 @@ subroutine write_options_ptmass(iunit)
  integer, intent(in) :: iunit
 
  write(iunit,"(/,a)") '# options controlling sink particles'
+ call write_inopt(isink_potential,'isink_potential','sink potential(0=1/r,1=surf)',iunit)
  if (gravity) then
     call write_inopt(icreate_sinks,'icreate_sinks','allow automatic sink particle creation',iunit)
     if (icreate_sinks > 0) then
@@ -1980,6 +1997,9 @@ subroutine read_options_ptmass(name,valstring,imatch,igotall,ierr)
     read(valstring,*,iostat=ierr) icreate_sinks
     ngot = ngot + 1
     if (icreate_sinks < 0) call fatal(label,'sink creation option out of range')
+ case('isink_potential')
+    read(valstring,*,iostat=ierr) isink_potential
+    ngot = ngot + 1
  case('rho_crit_cgs')
     read(valstring,*,iostat=ierr) rho_crit_cgs
     if (rho_crit_cgs < 0.) call fatal(label,'rho_crit < 0')

src/tests/test_ptmass.f90

@@ -40,7 +40,7 @@ subroutine test_ptmass(ntests,npass,string)
  integer, intent(inout) :: ntests,npass
  integer :: itmp,ierr,itest,istart
  logical :: do_test_binary,do_test_accretion,do_test_createsink,do_test_softening
- logical :: do_test_chinese_coin,do_test_merger
+ logical :: do_test_chinese_coin,do_test_merger,do_test_potential
  logical :: testall
 
  if (id==master) write(*,"(/,a,/)") '--> TESTING PTMASS MODULE'
@@ -50,6 +50,7 @@ subroutine test_ptmass(ntests,npass,string)
  do_test_createsink = .false.
  do_test_softening = .false.
  do_test_merger = .false.
+ do_test_potential = .false.
  do_test_chinese_coin = .false.
  testall = .false.
  istart = 1
@@ -64,6 +65,8 @@ subroutine test_ptmass(ntests,npass,string)
     do_test_softening = .true.
  case('ptmassmerger')
     do_test_merger = .true.
+ case('ptmasspotential')
+    do_test_potential = .true.
  case('ptmasschinchen','ptmasscoin','chinchen','coin','chinesecoin')
     do_test_chinese_coin = .true.
  case('ptmassfsi','fsi')
@@ -111,6 +114,10 @@ subroutine test_ptmass(ntests,npass,string)
     if (do_test_merger .or. testall) call test_merger(ntests,npass)
  enddo
  !
+ !  Test of sink particle potentials
+ !
+ if (do_test_potential .or. testall) call test_sink_potential(ntests,npass)
+ !
  !  Tests of accrete_particle routine
  !
  if (do_test_accretion .or. testall) then
@@ -1111,6 +1118,61 @@ subroutine test_merger(ntests,npass)
 
 end subroutine test_merger
 
+!-----------------------------------------------------------------------
+!+
+!  Test sink particle surface force, simply that the acceleration
+!  is the gradient of the potential
+!+
+!-----------------------------------------------------------------------
+subroutine test_sink_potential(ntests,npass)
+ use io,         only:id,master
+ use testutils,  only:checkval,update_test_scores
+ use ptmass,     only:get_accel_sink_gas,isink_potential
+ use part,       only:npart,npartoftype,nptmass,xyzmh_ptmass,ihacc,iReff
+ use units,      only:set_units
+ integer, intent(inout) :: ntests,npass
+ integer :: nfailed(1)
+ real :: phi1,phi,eps,x0(3)
+ real :: dphidx,hi,xi,yi,zi,dumxi,dumyi,dumzi,fxi,fyi,fzi,rp
+
+ if (id==master) write(*,"(/,a)") '--> testing sink particle surface force'
+ nptmass = 1
+ npart = 0
+ npartoftype = 0
+ hi = 0.
+ x0 = [100.,100.,100.]
+ rp = 2.
+ isink_potential = 1
+ ! place a single point mass at a random location
+ xyzmh_ptmass(:,:)     = 0.
+ xyzmh_ptmass(1:3,1)   = x0
+ xyzmh_ptmass(4,1)     = 3.14159
+ xyzmh_ptmass(ihacc,1) = 0.
+ xyzmh_ptmass(iReff,1) = rp  ! surface radius = 2
+
+ call set_units(mass=1.d0,dist=1.d0,G=1.d0)
+
+ ! evaluate sink-gas acceleration at some position
+ xi = x0(1) + 1.00001*rp
+ yi = x0(2) + 1.*rp
+ zi = x0(3) + 1.*rp
+ fxi = 0.; fyi = 0.; fzi = 0.; phi = 0.
+ call get_accel_sink_gas(nptmass,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi,phi)
+ ! evaluate sink-gas acceleration at some position + epsilon
+ eps = 1.e-6
+ dumxi = 0.; dumyi = 0.; dumzi = 0.; phi1 = 0.
+ call get_accel_sink_gas(nptmass,xi+eps,yi,zi,hi,xyzmh_ptmass,dumxi,dumyi,dumzi,phi1)
+ ! get the derivative of phi and check it equals the acceleration
+ dphidx = -(phi1 - phi)/eps
+
+ call checkval(dphidx,fxi,3.3e-8,nfailed(1),'dphi/dx = acceleration')
+ call update_test_scores(ntests,nfailed(1:1),npass)
+
+ ! reset options
+ isink_potential = 0
+
+end subroutine test_sink_potential
+
 !-----------------------------------------------------------------------
 !+
 !  Helper function used in sink particle creation test

src/tests/testsuite.F90

@@ -447,7 +447,7 @@ subroutine testsuite(string,first,last,ntests,npass,nfail)
           "|_| /_/   \_\____/____/ "
 
        write(*,"(a)") 'TEST SUITE PASSED'
-       call system("say OK")
+       call system("say fantastic!")
     else
        write(*,"(5(a,/))") &
           " _____ _    ___ _     ", &