Add the option to integrate ptmass with 4th order scheme (FSI)

.github/workflows/build.yml

@@ -35,7 +35,7 @@ jobs:
       nbatch: ${{ steps.set-sequence.outputs.nbatch }}
     steps:
     - name: Check out repo
-      uses: actions/checkout@v3
+      uses: actions/checkout@v4
     - name: Generate sequence of batch numbers for normal tests, or run sequentially for scheduled tests
       id: set-sequence
       run: |

.github/workflows/krome.yml

@@ -32,7 +32,7 @@ jobs:
           compiler: ${{ matrix.toolchain.compiler }}
 
       - name: "Clone phantom"
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
 
       - name: "Clone krome"
         run: git clone https://bitbucket.org/tgrassi/krome.git krome

.github/workflows/mpi.yml

@@ -53,7 +53,7 @@ jobs:
         sudo apt-get --yes install gfortran openmpi-bin openmpi-common libopenmpi-dev
 
     # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
     - name: Compile with MPI
       run: make SETUP=${{ matrix.input[0] }} MPI=yes DEBUG=${{ matrix.debug }} OPENMP=${{ matrix.openmp }}  phantomtest

.github/workflows/test.yml

@@ -87,7 +87,7 @@ jobs:
         printenv >> $GITHUB_ENV
 
     - name: "Clone phantom"
-      uses: actions/checkout@v3
+      uses: actions/checkout@v4
 
     - name: "Compile phantom"
       run: make SETUP=${{ matrix.input[0] }} DEBUG=${{ matrix.debug }} phantomtest

build/Makefile

@@ -537,7 +537,7 @@ SOURCES= physcon.f90 ${CONFIG} ${SRCKERNEL} io.F90 units.f90 \
          ${SRCKROME} memory.f90 ${SRCREADWRITE_DUMPS}  \
          quitdump.f90 ptmass.F90 \
          readwrite_infile.F90 dens.F90 force.F90 deriv.F90 energies.F90 sort_particles.f90 \
-         utils_shuffleparticles.F90 evwrite.f90 step_leapfrog.F90 writeheader.F90 ${SRCAN} step_supertimestep.F90 \
+         utils_shuffleparticles.F90 evwrite.f90 substepping.F90 step_leapfrog.F90 writeheader.F90 ${SRCAN} step_supertimestep.F90 \
          mf_write.f90 evolve.F90 utils_orbits.f90 utils_linalg.f90 \
          checksetup.f90 initial.F90
 
@@ -679,6 +679,7 @@ else
     SRCTESTMPI =
 endif
 
+# 22/4/24: added setup_params to avoid weird build failure with ifort on Mac OS
 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 \
@@ -689,7 +690,7 @@ SRCTESTS=utils_testsuite.f90 ${TEST_FASTMATH} test_kernel.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 \
-         testsuite.F90 phantomtest.f90
+         testsuite.F90 setup_params.f90 phantomtest.f90
 
 ifeq (X$(SRCTEST), X)
 SRCTEST=${SRCTESTS}

build/Makefile_setups

@@ -342,6 +342,12 @@ ifeq ($(SETUP), galcen)
     KNOWN_SETUP=yes
 endif
 
+ifeq ($(SETUP), starcluster)
+#   Cluster of stars (ptmass) 
+    SETUPFILE= setup_starcluster.f90
+    KNOWN_SETUP=yes
+endif
+
 #--- Bondi accretion/wind ---------------------------
 ifeq ($(SETUP), bondi)
 #   Bondi accretion flow

data/starcluster/README

@@ -0,0 +1 @@
+file with mass position and velocity... (should have 7 column per ptmass)

docs/developer-guide/vscode.rst

@@ -2,14 +2,16 @@ Coding Phantom in VSCode or Cursor AI
 =====================================
 
 In VSCode or Cursor there are several helpful settings when editing Phantom source files. After installing a modern Fortran extension, to enforce the indentation conventions in Phantom you should use `findent <https://github.com/wvermin/findent>`_ as in the indentation engine:
+and pass it the same options as used in `the bots script <https://github.com/danieljprice/phantom/blob/master/scripts/bots.sh#L288>`_:
 
-.. image:: ../images/vscode-findent.png
+.. image:: ../images/vscode-findent-flags.png
   :width: 800
-  :alt: findent option in VSCode
+  :alt: findent flags in VSCode
 
-and pass it the same options as used in `the bots script <https://github.com/danieljprice/phantom/blob/master/scripts/bots.sh#L288>`_:
+and yes, you do have to type each flag in a separate box. Then it is useful to select the "format on save" option  in Settings->Text Editor->Formatting:
 
-.. image:: ../images/vscode-findent-flags.png
+.. image:: ../images/vscode-format-on-save.png
   :width: 800
   :alt: findent flags in VSCode
 
+Thanks to Yann Bernard for getting this working!

scripts/kernels.py

@@ -84,6 +84,8 @@ def getkernelfuncs(w,R):
     #--work out the integration constant for the potential
     #
     parg = list(pot.args)
+    lastarg = len(pot.args) - 1
+    parg[lastarg] = (sympify(-1/(q)),pot.args[lastarg].cond)
     if isinstance(pot, Piecewise):
        for i, (e, c) in reversed(list(enumerate(pot.args))):
            if i < len(pot.args) - 1:
@@ -98,15 +100,28 @@ def getkernelfuncs(w,R):
     #--derivative of potential with respect to h
     #
     dpotdh = pot
-    parg = list(pot.args)
+    pharg = list(pot.args)
     if isinstance(pot, Piecewise):
        for i, (e, c) in enumerate(pot.args):
            ep = simplify(-e - q*diff(e,q))
-           parg[i] = (ep, c)
-       tuple(parg)
-       dpotdh = Piecewise(*parg)
+           pharg[i] = (ep, c)
+       tuple(pharg)
+       dpotdh = Piecewise(*pharg)
 
-    return (dw, d2w, c1D, c2D, c3D, fsoft, pot, dpotdh)
+    #
+    #--kernel function needed in gradient acceleration 
+    #  for 4th order Forward Symplectic Integrator
+    #
+    farg = list(fsoft.args)
+    if isinstance(fsoft, Piecewise):
+       for i, (e, c) in enumerate(fsoft.args):
+           ep = simplify(q*diff(e,q) - e)
+           farg[i] = (ep, c)
+       tuple(farg)
+       gsoft = Piecewise(*farg)
+
+    #gsoft = piecewise_fold(simplify(diff(q*fsoft,q) - fsoft))
+    return (dw, d2w, c1D, c2D, c3D, fsoft, pot, dpotdh, gsoft)
 
 #---------------------------------------------
 # function to get the variance of the kernel
@@ -225,11 +240,12 @@ def printvariances(w,R):
 # function to print basic kernel information to the screen
 #-----------------------------------------------------------
 def printkernel(w,R):
-    dw, d2w, c1D, c2D, c3D, fsoft, pot, dpotdh = getkernelfuncs(w,R)
+    dw, d2w, c1D, c2D, c3D, fsoft, pot, dpotdh, gsoft = getkernelfuncs(w,R)
     print ("\n%s W:" %name)
     print (w)
     print ("\nFirst derivative:")
     print (dw)
+    #print (fmt(dw))
     print ("\n2nd derivative:")
     print (d2w)
     print ("\nnormalisation:")
@@ -241,6 +257,8 @@ def printkernel(w,R):
     avnorm = -pi/8*c2D*integrate(q*q*dw,(q,0,R))
     print (avnorm)
     printvariances(w,R)
+    print ("\n gradient acceleration term:")
+    print (gsoft)
     return
 
 #-------------------------------------------------------------
@@ -290,6 +308,10 @@ def fmt(e):
     # replace 15*x with 15.*x as long as it is not **15*x
     s = re.sub("(?!\*\d+)(\D\d+)\*","\g<1>.*", s)
 
+    # replace " 2)" with " 2.)"
+    # Use re.sub to replace " digit)" with " digit.)"
+    s = re.sub(r" (\d)\)", r" \1.)", s)
+
     f = sympify(s)
     #
     # expand if it makes it shorter
@@ -302,6 +324,9 @@ def fmt(e):
 
     # replace 1.4000000 with 1.4
     g = re.sub("(\.[1-9]*)(0+)(\D|$)","\g<1>\g<3>", g)
+    # replace " 2)" with " 2.)"
+    # Use re.sub to replace " digit)" with " digit.)"
+    g = re.sub(r" (\d)\)", r" \1.)", g)
 
     # only return simplify-ed strings if no fully expanded floats 0.345242545..
     if re.search("(\.\d\d\d\d\d+)",g):
@@ -614,7 +639,7 @@ def print_decl(w):
 #---------------------------------
 def printkernel_phantom(w,R,name):
     import datetime
-    dw, d2w, c1D, c2D, c3D, fsoft, pot, dpotdh = getkernelfuncs(w,R)
+    dw, d2w, c1D, c2D, c3D, fsoft, pot, dpotdh, gsoft = getkernelfuncs(w,R)
     w0 = w.subs(q,0)
     dpotdh0 = dpotdh.subs(q,0)
     #print("GOT dpotdh0",simplify(dpotdh0))
@@ -627,31 +652,26 @@ def printkernel_phantom(w,R,name):
     lb = "!"+"-"*62
     print ("!--------------------------------------------------------------------------!")
     print ("! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. !")
-    print ("! Copyright (c) 2007-2019 The Authors (see AUTHORS)                        !")
+    print ("! Copyright (c) 2007-2024 The Authors (see AUTHORS)                        !")
     print ("! See LICENCE file for usage and distribution conditions                   !")
-    print ("! http://phantomsph.bitbucket.io/                                          !")
+    print ("! http://phantomsph.github.io/                                             !")
     print ("!--------------------------------------------------------------------------!")
-    print ("!+")
-    print ("!  MODULE: kernel")
+    print ("module kernel")
     print ("!")
-    print ("!  DESCRIPTION:")
-    print ("!   This module implements the %s kernel" %name)
+    print ("! This module implements the %s kernel" %name)
     print ("!   DO NOT EDIT - auto-generated by kernels.py")
     print ("!")
-    print ("!  REFERENCES: None")
+    print ("! :References: None")
     print ("!")
-    print ("!  OWNER: Daniel Price")
+    print ("! :Owner: Daniel Price")
     print ("!")
-    print ("!  $Id:$")
+    print ("! :Runtime parameters: None")
     print ("!")
-    print ("!  RUNTIME PARAMETERS: None")
+    print ("! :Dependencies: physcon")
     print ("!")
-    print ("!  DEPENDENCIES: physcon")
+    print ("! :Generated:",datetime.datetime.now())
     print ("!")
-    print ("!  GENERATED:",datetime.datetime.now())
-    print ("!+")
     print ("!--------------------------------------------------------------------------")
-    print ("module kernel")
     print (" use physcon, only:pi")
     print (" implicit none")
     print (" character(len=%i), public :: kernelname = '%s'" %(len(name),name))
@@ -660,9 +680,9 @@ def printkernel_phantom(w,R,name):
     print (" real, parameter, public  :: cnormk = %s" %fmt(c3D))
     print (" real, parameter, public  :: wab0 = %s, gradh0 = -3.*wab0" %fmt(w0))
     print (" real, parameter, public  :: dphidh0 = %s" %fmtp(dpotdh0))
-    print (" real, parameter, public  :: cnormk_drag = %s " %fmt(c3Ddrag))
+    print (" real, parameter, public  :: cnormk_drag = %s" %fmt(c3Ddrag))
     var, relvar, reldev = getvar(w,R)
-    print (" real, parameter, public  :: hfact_default = %.1f " %(1.2/reldev[2]))
+    print (" real, parameter, public  :: hfact_default = %.1f" %(1.2/reldev[2]))
     #print " real, parameter, public  :: hfact_default = %s " %fmt(reldev[2])
     print (" real, parameter, public  :: av_factor = %s" %fmt(avratio))
     print ("\ncontains\n")
@@ -774,7 +794,7 @@ def printkernel_phantom(w,R,name):
     print ("pure subroutine kernel_softening(q2,q,potensoft,fsoft)")
     print (" real, intent(in)  :: q2,q")
     print (" real, intent(out) :: potensoft,fsoft")
-    print_decl(pot)
+    print_decl(fsoft)
     if isinstance(dw, Piecewise):
        for i, (de, c) in enumerate(dw.args):
            (pote, potc) = pot.args[i]
@@ -793,6 +813,30 @@ def printkernel_phantom(w,R,name):
        print ("    potensoft = %s" %fmtp(pot))
        print ("    fsoft     = %s" %fmtp(fsoft))
     print ("\nend subroutine kernel_softening\n")
+
+    print ("!------------------------------------------")
+    print ("! gradient acceleration kernel needed for")
+    print ("! use in Forward symplectic integrator")
+    print ("!------------------------------------------")
+    print ("pure subroutine kernel_grad_soft(q2,q,gsoft)")
+    print (" real, intent(in)  :: q2,q")
+    print (" real, intent(out) :: gsoft")
+    print_decl(gsoft)
+    if isinstance(dw, Piecewise):
+       for i, (de, c) in enumerate(dw.args):
+           (ge, gc) = gsoft.args[i]
+           if i == 0:
+              print (" if (%s) then" %fmt(c))
+           elif i == len(dw.args)-1 and c == True:
+              print (" else")
+           else:
+              print (" elseif (%s) then" %fmt(c))
+           print_defs(4,fmtp(ge))
+           print ("    gsoft = %s" %fmtp(ge))
+       print (" endif")
+    else:
+       print ("    gsoft = %s" %fmtp(gsoft))
+    print ("\nend subroutine kernel_grad_soft\n")
     print ("!------------------------------------------")
     print ("! double-humped version of the kernel for")
     print ("! use in drag force calculations")
@@ -956,8 +1000,8 @@ def f6(R):
 
 # define which kernel to use
 #f, name = sinq(R,3)
-#f, name = m5(R)
-f, name = w6(R)
+f, name = m4(R)
+#f, name = w6(R)
 
 #print_avdiss(f,R)
 #printvariances(f,R)

src/main/checksetup.f90

@@ -429,6 +429,10 @@ subroutine check_setup(nerror,nwarn,restart)
 !--check centre of mass
 !
  call get_centreofmass(xcom,vcom,npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass)
+!
+!--check Forward symplectic integration method imcompatiblity
+!
+ call check_vdep_extf (nwarn,iexternalforce)
 
  if (.not.h2chemistry .and. maxvxyzu >= 4 .and. icooling == 3 .and. iexternalforce/=iext_corotate .and. nptmass==0) then
     if (dot_product(xcom,xcom) >  1.e-2) then
@@ -522,11 +526,15 @@ subroutine check_setup_ptmass(nerror,nwarn,hmin)
  use part, only:nptmass,xyzmh_ptmass,ihacc,ihsoft,gr,iTeff,sinks_have_luminosity,&
                 ilum,iJ2,ispinx,ispinz,iReff
  use ptmass_radiation, only:isink_radiation
+ use ptmass, only:use_fourthorder
  integer, intent(inout) :: nerror,nwarn
  real,    intent(in)    :: hmin
  integer :: i,j,n
  real :: dx(3)
  real :: r,hsink,hsoft,J2
+ logical :: isoblate
+
+ isoblate = .false.
 
  if (gr .and. nptmass > 0) then
     print*,' ERROR: nptmass = ',nptmass, ' should be = 0 for GR'
@@ -615,6 +623,7 @@ subroutine check_setup_ptmass(nerror,nwarn,hmin)
     ! in order to specify the rotation direction
     !
     if (J2 > 0.) then
+       isoblate = .true.
        if (dot_product(xyzmh_ptmass(ispinx:ispinz,i),xyzmh_ptmass(ispinx:ispinz,i)) < tiny(0.)) then
           nerror = nerror + 1
           print*,'ERROR! non-zero J2 requires non-zero spin on sink particle ',i
@@ -625,6 +634,13 @@ subroutine check_setup_ptmass(nerror,nwarn,hmin)
        endif
     endif
  enddo
+
+ if (isoblate .and. use_fourthorder) then
+    nwarn = nwarn + 1
+    print*, 'WARNING: Substepping integration switched back to leapfrog due to oblateness'
+    use_fourthorder = .false.
+ endif
+
  !
  !  check that radiation properties are sensible
  !
@@ -999,4 +1015,17 @@ subroutine check_setup_radiation(npart,nerror,nwarn,radprop,rad)
 
 end subroutine check_setup_radiation
 
+subroutine check_vdep_extf(nwarn,iexternalforce)
+ use externalforces, only: is_velocity_dependent
+ use ptmass, only : use_fourthorder
+ integer, intent(inout) :: nwarn
+ integer, intent(in)    :: iexternalforce
+ if (is_velocity_dependent(iexternalforce) .and. use_fourthorder) then
+    print "(/,a,/)","Warning: velocity dependant external forces are not compatible with FSI switch back to Leapfrog..."
+    nwarn = nwarn + 1
+    use_fourthorder = .false.
+ endif
+
+end subroutine check_vdep_extf
+
 end module checksetup

src/main/config.F90

@@ -121,7 +121,6 @@ module dim
 #else
  logical, parameter :: do_radiation = .false.
 #endif
-
  ! rhosum
  integer, parameter :: maxrhosum = 39 + &
                                    maxdustlarge - 1 + &

src/main/cons2prim.f90

@@ -123,20 +123,20 @@ subroutine cons2primall(npart,xyzh,metrics,pxyzu,vxyzu,dens,eos_vars)
  use part,            only:isdead_or_accreted,massoftype,igas,rhoh,igasP,ics,ien_type,&
                            itemp,igamma
  use io,              only:fatal
- use eos,             only:ieos,gamma,done_init_eos,init_eos,get_spsound
+ use eos,             only:ieos,done_init_eos,init_eos,get_spsound
  integer, intent(in)    :: npart
  real,    intent(in)    :: pxyzu(:,:),xyzh(:,:),metrics(:,:,:,:)
  real,    intent(inout) :: vxyzu(:,:),dens(:)
  real,    intent(out)   :: eos_vars(:,:)
  integer :: i, ierr
- real    :: p_guess,rhoi,pondens,spsound,tempi,gammai
+ real    :: p_guess,rhoi,tempi,gammai
 
  if (.not.done_init_eos) call init_eos(ieos,ierr)
 
 !$omp parallel do default (none) &
 !$omp shared(xyzh,metrics,vxyzu,dens,pxyzu,npart,massoftype) &
-!$omp shared(ieos,gamma,eos_vars,ien_type) &
-!$omp private(i,ierr,spsound,pondens,p_guess,rhoi,tempi,gammai)
+!$omp shared(ieos,eos_vars,ien_type) &
+!$omp private(i,ierr,p_guess,rhoi,tempi,gammai)
  do i=1,npart
     if (.not.isdead_or_accreted(xyzh(4,i))) then
        ! get pressure, temperature and gamma from previous step as the initial guess