Merge pull request #404 from danieljprice/damping-bcs

Damping bcs

build/Makefile

@@ -639,7 +639,7 @@ SRCLIBSETUP=physcon.f90 geometry.f90 random.f90 utils_tables.f90 utils_vectors.f
             utils_binary.f90 set_binary.f90 set_flyby.f90 \
             set_hierarchical_utils.f90 set_hierarchical.f90 \
             set_unifdis.f90 set_sphere.f90 set_shock.f90 \
-            set_dust.f90 libsetup.f90 
+            set_dust.f90 libsetup.f90
 OBJLIBSETUP=${SRCLIBSETUP:.f90=.o}
 
 libsetup: $(OBJLIBSETUP)
@@ -692,9 +692,9 @@ SRCTESTS=utils_testsuite.f90 ${TEST_FASTMATH} test_kernel.f90 \
          test_dust.F90 test_growth.F90 test_smol.F90 \
          test_nonidealmhd.F90 directsum.f90 test_gravity.F90 \
          test_derivs.F90 test_cooling.f90 test_eos_stratified.f90 \
-	 test_eos.f90 test_externf.f90 test_rwdump.f90 \
+         test_eos.f90 test_externf.f90 test_rwdump.f90 \
          test_step.F90 test_indtstep.F90 set_disc.F90 test_setdisc.F90 \
-         test_hierarchical.f90 \
+         test_hierarchical.f90 test_damping.f90 \
          test_link.F90 test_kdtree.F90 test_part.f90 test_ptmass.f90 test_luminosity.F90\
          test_gnewton.f90 test_corotate.f90 test_geometry.f90 \
          ${SRCTESTMPI} test_sedov.F90 test_poly.f90 test_radiation.F90 \

src/main/damping.f90

@@ -29,7 +29,7 @@ module damping
 !
  implicit none
 
- public  :: calc_damp,apply_damp
+ public  :: calc_damp,apply_damp,get_damp_fac_disc
  public  :: write_options_damping,read_options_damping
 
  private
@@ -55,11 +55,12 @@ subroutine calc_damp(time, damp_fac)
  use physcon, only:pi
  real, intent(out)   :: damp_fac
  real, intent(in)    :: time
- real                :: tau1, tau2, tdyn_star
+ real                :: tau1, tau2, tdyn_star, orbital_period
 
  select case(idamp)
  case(3)
-    damp_fac = 1./(damp*2.*pi*sqrt(r1in**3)) ! fraction of orbital time at r=r1in with G=M=1
+    orbital_period = 2.*pi*sqrt(r1in**3)  ! G=M=1
+    damp_fac = damp/orbital_period ! fraction of orbital time at r=r1in with G=M=1
  case(2)
     tdyn_star = tdyn_s / utime
     tau1 = tdyn_star * 0.1

src/tests/test_damping.f90

@@ -0,0 +1,97 @@
+!--------------------------------------------------------------------------!
+! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. !
+! Copyright (c) 2007-2023 The Authors (see AUTHORS)                        !
+! See LICENCE file for usage and distribution conditions                   !
+! http://phantomsph.bitbucket.io/                                          !
+!--------------------------------------------------------------------------!
+module testdamping
+!
+! Unit tests of damping module
+!
+! :References: exoALMA comparison project
+!
+! :Owner: Daniel Price
+!
+! :Runtime parameters: None
+!
+! :Dependencies: io
+!
+ implicit none
+ public :: test_damping
+
+ private
+
+contains
+!-----------------------------------------------------------------------
+!+
+!   Unit tests of the routines in the damping module
+!+
+!-----------------------------------------------------------------------
+subroutine test_damping(ntests,npass)
+ use io, only:id,master
+ integer, intent(inout) :: ntests,npass
+
+ if (id==master) write(*,"(/,a,/)") '--> TESTING DAMPING TERMS'
+
+ call test_damping_disc(ntests,npass)
+
+ if (id==master) write(*,"(/,a)") '<-- DAMPING TESTS COMPLETE'
+
+end subroutine test_damping
+
+!-----------------------------------------------------------------------
+!+
+!   Unit tests of the disc damping boundary conditions
+!+
+!-----------------------------------------------------------------------
+subroutine test_damping_disc(ntests,npass)
+ use damping,   only:idamp,damp,get_damp_fac_disc,r1in,r2in,r1out,r2out,calc_damp
+ use physcon,   only:pi
+ use testutils, only:checkval,update_test_scores
+ integer, intent(inout) :: ntests,npass
+ integer :: nfail(1),ipos
+ real :: xyz(3),v0(3),fac,r,damp_fac
+ real :: omega,t_orb,t_damp,time
+ real, parameter :: tol = 3.e-16
+ real :: rpos(4)
+ character(len=5) :: label(4)
+
+ ! select disc damping
+ idamp = 3
+ time = 0.
+ damp = 0.01
+
+ ! set positions of damping zones
+ r1in = 0.3; r2in = 2.
+ r1out = 2.52; r2out = 3.0
+ call calc_damp(time,damp_fac)
+
+ ! set position at inner edge of first damping zone
+ rpos  = [r1in,r2out,r1out,r2in]
+ label = ['r1in ','r2out','r1out','r2in ']
+
+ do ipos=1,size(rpos)
+    r = rpos(ipos)
+    xyz = [r,0.,0.]
+    fac = damp_fac*get_damp_fac_disc(xyz,v0)
+
+    ! check the damping time is the correct multiple
+    ! of the orbital time at this radius
+    omega = sqrt(1./r**3)
+    t_orb = 2.*pi/omega
+
+    if (ipos==3 .or. ipos==4) then
+       call checkval(fac,0.,tol,nfail(1),'damping = 0 at r='//trim(label(ipos)))
+    else
+       t_damp = 1./fac
+       call checkval(t_damp,t_orb/damp,tol,nfail(1),'t_damp = f*t_orb at r='//trim(label(ipos)))
+    endif
+    call update_test_scores(ntests,nfail,npass)
+
+    call checkval(v0(2),r*omega,tol,nfail(1),'v = v_kep at r='//trim(label(ipos)))
+    call update_test_scores(ntests,nfail,npass)
+ enddo
+
+end subroutine test_damping_disc
+
+end module testdamping

src/tests/test_externf.f90

@@ -43,6 +43,7 @@ subroutine test_externf(ntests,npass)
  use units,    only:set_units
  use physcon,  only:pc,solarm
  use mpidomain,only:i_belong
+ use kernel,   only:hfact_default
  integer, intent(inout) :: ntests,npass
  integer                :: i,iextf,nfail1,ierr
  logical                :: dotest1,dotest2,dotest3,accreted
@@ -65,7 +66,7 @@ subroutine test_externf(ntests,npass)
 !
  xmini(:) = -100.
  xmaxi(:) = 100.
- hfact      = 1.2
+ hfact      = hfact_default
  accradius1 = 100.  ! should be >6 for Lense-Thirring to pass
  psep  = (xmaxi(1) - xmini(1))/10.
  npart = 0

src/tests/testsuite.F90

@@ -62,21 +62,22 @@ subroutine testsuite(string,first,last,ntests,npass,nfail)
  use testeos,      only:test_eos
  use testcooling,  only:test_cooling
  use testgeometry, only:test_geometry
+ use testpoly,     only:test_poly
+ use testdamping,  only:test_damping
+ use testradiation,only:test_radiation
 #ifdef MPI
  use testmpi,      only:test_mpi
 #endif
  use timing,       only:get_timings,print_time
- use mpiutils,      only:barrier_mpi
- use testradiation, only:test_radiation
- use dim,           only:do_radiation
- use testpoly,      only:test_poly
+ use mpiutils,     only:barrier_mpi
+ use dim,          only:do_radiation
  character(len=*), intent(in)    :: string
  logical,          intent(in)    :: first,last
  integer,          intent(inout) :: ntests,npass,nfail
  logical :: testall,dolink,dokdtree,doderivs,dokernel,dostep,dorwdump,dosmol
  logical :: doptmass,dognewton,dosedov,doexternf,doindtstep,dogravity,dogeom
  logical :: dosetdisc,doeos,docooling,dodust,donimhd,docorotate,doany,dogrowth
- logical :: dogr,doradiation,dopart,dopoly,dompi,dohier
+ logical :: dogr,doradiation,dopart,dopoly,dompi,dohier,dodamp
 #ifdef FINVSQRT
  logical :: usefsqrt,usefinvsqrt
 #endif
@@ -127,6 +128,7 @@ subroutine testsuite(string,first,last,ntests,npass,nfail)
  dopoly     = .false.
  dompi      = .false.
  dohier     = .false.
+ dodamp     = .false.
 
  if (index(string,'deriv')     /= 0) doderivs  = .true.
  if (index(string,'grav')      /= 0) dogravity = .true.
@@ -146,10 +148,11 @@ subroutine testsuite(string,first,last,ntests,npass,nfail)
  if (index(string,'poly')      /= 0) dopoly    = .true.
  if (index(string,'mpi')       /= 0) dompi     = .true.
  if (index(string,'hier')      /= 0) dohier    = .true.
+ if (index(string,'damp')      /= 0) dodamp    = .true.
 
  doany = any((/doderivs,dogravity,dodust,dogrowth,donimhd,dorwdump,&
                doptmass,docooling,dogeom,dogr,dosmol,doradiation,&
-               dopart,dopoly,dohier/))
+               dopart,dopoly,dohier,dodamp/))
 
  select case(trim(string))
  case('kernel','kern')
@@ -367,13 +370,20 @@ subroutine testsuite(string,first,last,ntests,npass,nfail)
     call set_default_options_testsuite(iverbose) ! restore defaults
  endif
 !
-!--test of geometry module
+!--test of polynomial solvers
 !
  if (dopoly.or.testall) then
     call test_poly(ntests,npass)
     call set_default_options_testsuite(iverbose) ! restore defaults
  endif
 !
+!--test of damping module
+!
+ if (dodamp.or.testall) then
+    call test_damping(ntests,npass)
+    call set_default_options_testsuite(iverbose) ! restore defaults
+ endif
+!
 !--test of radiation module
 !
  if (doradiation.or.testall) then