Merge pull request #536 from Yrisch/subsystems_SDAR

Add subsystem integration using regularization method like TTL or logH (Slow down will come in an other PR)

build/Makefile

@@ -535,6 +535,7 @@ SOURCES= physcon.f90 ${CONFIG} ${SRCKERNEL} io.F90 units.f90 \
          utils_deriv.f90 utils_implicit.f90 radiation_implicit.f90 ${SRCTURB} \
          ${SRCINJECT_DEPS} wind_equations.f90 wind.F90 ${SRCINJECT} \
          ${SRCKROME} memory.f90 ${SRCREADWRITE_DUMPS}  \
+         utils_subgroup.f90 utils_kepler.f90 subgroup.f90\
          quitdump.f90 ptmass.F90 \
          readwrite_infile.F90 dens.F90 force.F90 deriv.F90 energies.F90 sort_particles.f90 \
          utils_shuffleparticles.F90 evwrite.f90 substepping.F90 step_leapfrog.F90 writeheader.F90 ${SRCAN} step_supertimestep.F90 \

src/main/checksetup.f90

@@ -433,6 +433,10 @@ subroutine check_setup(nerror,nwarn,restart)
 !--check Forward symplectic integration method imcompatiblity
 !
  call check_vdep_extf (nwarn,iexternalforce)
+!
+!--check Regularization imcompatibility
+!
+ call check_regnbody (nerror)
 
  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
@@ -1028,4 +1032,14 @@ subroutine check_vdep_extf(nwarn,iexternalforce)
 
 end subroutine check_vdep_extf
 
+subroutine check_regnbody (nerror)
+ use ptmass, only:use_regnbody,use_fourthorder
+ integer, intent(inout) :: nerror
+ if (use_regnbody .and. .not.(use_fourthorder)) then
+    print "(/,a,/)","Error: TTL integration and regularization tools are not available without FSI. Turn off TTL..."
+    nerror = nerror + 1
+ endif
+end subroutine check_regnbody
+
+
 end module checksetup

src/main/energies.F90

@@ -71,7 +71,7 @@ subroutine compute_energies(t)
                           isdead_or_accreted,epot_sinksink,imacc,ispinx,ispiny,&
                           ispinz,mhd,gravity,poten,dustfrac,eos_vars,itemp,igasP,ics,&
                           nden_nimhd,eta_nimhd,iion,ndustsmall,graindens,grainsize,&
-                          iamdust,ndusttypes,rad,iradxi
+                          iamdust,ndusttypes,rad,iradxi,gtgrad,group_info,n_group
  use part,           only:pxyzu,fxyzu,fext
  use gravwaveutils,  only:calculate_strain,calc_gravitwaves
  use centreofmass,   only:get_centreofmass_accel
@@ -81,7 +81,8 @@ subroutine compute_energies(t)
  use externalforces, only:externalforce,externalforce_vdependent,was_accreted,accradius1
  use options,        only:iexternalforce,calc_erot,alpha,ieos,use_dustfrac
  use mpiutils,       only:reduceall_mpi
- use ptmass,         only:get_accel_sink_gas
+ use ptmass,         only:get_accel_sink_gas,use_regnbody
+ use subgroup,     only:get_pot_subsys
  use viscosity,      only:irealvisc,shearfunc
  use nicil,          only:nicil_update_nimhd,nicil_get_halldrift,nicil_get_ambidrift, &
                      use_ohm,use_hall,use_ambi,n_data_out,n_warn,eta_constant
@@ -601,7 +602,14 @@ subroutine compute_energies(t)
  emag = reduceall_mpi('+',emag)
  epot = reduceall_mpi('+',epot)
  erad = reduceall_mpi('+',erad)
- if (nptmass > 1) epot = epot + epot_sinksink
+ if (nptmass > 1) then
+    if (use_regnbody) then
+       call get_pot_subsys(n_group,group_info,xyzmh_ptmass,fxyz_ptmass,gtgrad,epot_sinksink)
+    endif
+    epot = epot + epot_sinksink
+ endif
+
+
 
  etot = ekin + etherm + emag + epot + erad
 

src/main/initial.F90

@@ -130,7 +130,8 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
                             nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,&
                             epot_sinksink,get_ntypes,isdead_or_accreted,dustfrac,ddustevol,&
                             nden_nimhd,dustevol,rhoh,gradh, &
-                            Bevol,Bxyz,dustprop,filfac,ddustprop,ndustsmall,iboundary,eos_vars,dvdx
+                            Bevol,Bxyz,dustprop,filfac,ddustprop,ndustsmall,iboundary,eos_vars,dvdx, &
+                            n_group,n_ingroup,n_sing,nmatrix,group_info
  use part,             only:pxyzu,dens,metrics,rad,radprop,drad,ithick
  use densityforce,     only:densityiterate
  use linklist,         only:set_linklist
@@ -150,7 +151,8 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
  use nicil_sup,        only:use_consistent_gmw
  use ptmass,           only:init_ptmass,get_accel_sink_gas,get_accel_sink_sink, &
                             h_acc,r_crit,r_crit2,rho_crit,rho_crit_cgs,icreate_sinks, &
-                            r_merge_uncond,r_merge_cond,r_merge_uncond2,r_merge_cond2,r_merge2
+                            r_merge_uncond,r_merge_cond,r_merge_uncond2,r_merge_cond2,r_merge2, &
+                            use_regnbody
  use timestep,         only:time,dt,dtextforce,C_force,dtmax,dtmax_user,idtmax_n
  use timing,           only:get_timings
  use timestep_ind,     only:ibinnow,maxbins,init_ibin,istepfrac
@@ -210,6 +212,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
  use checkconserved,   only:get_conserv,etot_in,angtot_in,totmom_in,mdust_in
  use fileutils,        only:make_tags_unique
  use damping,          only:idamp
+ use subgroup,       only:group_identify
  character(len=*), intent(in)  :: infile
  character(len=*), intent(out) :: logfile,evfile,dumpfile
  logical,          intent(in), optional :: noread
@@ -495,10 +498,15 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
  endif
  if (nptmass > 0) then
     if (id==master) write(iprint,"(a,i12)") ' nptmass       = ',nptmass
-
     ! compute initial sink-sink forces and get timestep
-    call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,dtsinksink,&
+    if (use_regnbody) then
+       call group_identify(nptmass,n_group,n_ingroup,n_sing,xyzmh_ptmass,vxyz_ptmass,group_info,nmatrix)
+       call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,dtsinksink,&
+                             iexternalforce,time,merge_ij,merge_n,dsdt_ptmass,group_info=group_info)
+    else
+       call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,dtsinksink,&
                              iexternalforce,time,merge_ij,merge_n,dsdt_ptmass)
+    endif
     dtsinksink = C_force*dtsinksink
     if (id==master) write(iprint,*) 'dt(sink-sink) = ',dtsinksink
     dtextforce = min(dtextforce,dtsinksink)

src/main/part.F90

@@ -291,6 +291,20 @@ module part
 !
  real(kind=4), allocatable :: luminosity(:)
 !
+!-- Regularisation algorithm allocation
+!
+ integer,         allocatable :: group_info(:,:)
+ integer(kind=1), allocatable :: nmatrix(:,:)
+ integer, parameter :: igarg  = 1 ! idx of the particle member of a group
+ integer, parameter :: igcum  = 2 ! cumulative sum of the indices to find the starting and ending point of a group
+ integer, parameter :: igid   = 3 ! id of the group, correspond to the root of the group in the dfs/union find construction
+ ! needed for group identification and sorting
+ integer  :: n_group = 0
+ integer  :: n_ingroup = 0
+ integer  :: n_sing = 0
+ ! Gradient of the time transformation function
+ real, allocatable :: gtgrad(:,:)
+!
 !--derivatives (only needed if derivs is called)
 !
  real, allocatable         :: fxyzu(:,:)
@@ -480,6 +494,9 @@ subroutine allocate_part
     call allocate_array('abundance', abundance, nabundances, maxp_h2)
  endif
  call allocate_array('T_gas_cool', T_gas_cool, maxp_krome)
+ call allocate_array('group_info', group_info, 3, maxptmass)
+ call allocate_array("nmatrix", nmatrix, maxptmass, maxptmass)
+ call allocate_array("gtgrad", gtgrad, 3, maxptmass)
 
 
 end subroutine allocate_part
@@ -560,6 +577,9 @@ subroutine deallocate_part
  if (allocated(ibelong))      deallocate(ibelong)
  if (allocated(istsactive))   deallocate(istsactive)
  if (allocated(ibin_sts))     deallocate(ibin_sts)
+ if (allocated(group_info))   deallocate(group_info)
+ if (allocated(nmatrix))      deallocate(nmatrix)
+ if (allocated(gtgrad))       deallocate(gtgrad)
 
 end subroutine deallocate_part
 

src/main/ptmass.F90

@@ -65,11 +65,12 @@ module ptmass
  real,    public :: f_acc  = 0.8
  real,    public :: h_soft_sinkgas  = 0.0
  real,    public :: h_soft_sinksink = 0.0
- real,    public :: r_merge_uncond  = 0.0    ! sinks will unconditionally merge if they touch
- real,    public :: r_merge_cond    = 0.0    ! sinks will merge if bound within this radius
- real,    public :: f_crit_override = 0.0    ! 1000.
+ real,    public :: r_merge_uncond  = 0.0     ! sinks will unconditionally merge if they touch
+ 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
  real, public, parameter :: dk2(3) = (/0.5,0.5,0.0/)
@@ -299,14 +300,15 @@ end subroutine get_accel_sink_gas
 !+
 !----------------------------------------------------------------
 subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksink,&
-            iexternalforce,ti,merge_ij,merge_n,dsdt_ptmass,extrapfac,fsink_old)
+            iexternalforce,ti,merge_ij,merge_n,dsdt_ptmass,extrapfac,fsink_old,group_info)
 #ifdef FINVSQRT
  use fastmath,       only:finvsqrt
 #endif
  use externalforces, only:externalforce
  use extern_geopot,  only:get_geopot_force
  use kernel,         only:kernel_softening,radkern
  use vectorutils,    only:unitvec
+ use part,           only:igarg,igid
  integer,           intent(in)  :: nptmass
  real,              intent(in)  :: xyzmh_ptmass(nsinkproperties,nptmass)
  real,              intent(out) :: fxyz_ptmass(4,nptmass)
@@ -317,15 +319,16 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin
  real,              intent(out) :: dsdt_ptmass(3,nptmass)
  real,    optional, intent(in)  :: extrapfac
  real,    optional, intent(in)  :: fsink_old(4,nptmass)
+ integer, optional, intent(in)  :: group_info(3,nptmass)
  real    :: xi,yi,zi,pmassi,pmassj,fxi,fyi,fzi,phii
  real    :: ddr,dx,dy,dz,rr2,rr2j,dr3,f1,f2
  real    :: hsoft1,hsoft21,q2i,qi,psoft,fsoft
  real    :: fextx,fexty,fextz,phiext !,hsofti
  real    :: fterm,pterm,potensoft0,dsx,dsy,dsz
  real    :: J2i,rsinki,shati(3)
  real    :: J2j,rsinkj,shatj(3)
- integer :: i,j
- logical :: extrap
+ integer :: k,l,i,j,gidi,gidj
+ logical :: extrap,subsys
 
  dtsinksink = huge(dtsinksink)
  fxyz_ptmass(:,:) = 0.
@@ -340,6 +343,12 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin
  else
     extrap = .false.
  endif
+
+ if (present(group_info)) then
+    subsys = .true.
+ else
+    subsys = .false.
+ endif
  !
  !--get self-contribution to the potential if sink-sink softening is used
  !
@@ -356,19 +365,26 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin
  !--compute N^2 forces on point mass particles due to each other
  !
  !$omp parallel do default(none) &
- !$omp shared(nptmass,xyzmh_ptmass,fxyz_ptmass,merge_ij,r_merge2,dsdt_ptmass) &
+ !$omp shared(nptmass,xyzmh_ptmass,fxyz_ptmass,merge_ij,r_merge2,dsdt_ptmass,group_info,subsys) &
  !$omp shared(iexternalforce,ti,h_soft_sinksink,potensoft0,hsoft1,hsoft21) &
  !$omp shared(extrapfac,extrap,fsink_old) &
- !$omp private(i,xi,yi,zi,pmassi,pmassj) &
+ !$omp private(i,j,xi,yi,zi,pmassi,pmassj) &
+ !$omp private(gidi,gidj) &
  !$omp private(dx,dy,dz,rr2,rr2j,ddr,dr3,f1,f2) &
  !$omp private(fxi,fyi,fzi,phii,dsx,dsy,dsz) &
  !$omp private(fextx,fexty,fextz,phiext) &
  !$omp private(q2i,qi,psoft,fsoft) &
  !$omp private(fterm,pterm,J2i,J2j,shati,shatj,rsinki,rsinkj) &
  !$omp reduction(min:dtsinksink) &
  !$omp reduction(+:phitot,merge_n)
- do i=1,nptmass
-    if (extrap) then
+ do k=1,nptmass
+    if (subsys) then
+       i = group_info(igarg,k)
+       gidi = group_info(igid,k)
+    else
+       i = k
+    endif
+    if (extrap)then
        xi     = xyzmh_ptmass(1,i) + extrapfac*fsink_old(1,i)
        yi     = xyzmh_ptmass(2,i) + extrapfac*fsink_old(2,i)
        zi     = xyzmh_ptmass(3,i) + extrapfac*fsink_old(3,i)
@@ -389,7 +405,14 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin
     dsx    = 0.
     dsy    = 0.
     dsz    = 0.
-    do j=1,nptmass
+    do l=1,nptmass
+       if (subsys) then
+          j = group_info(igarg,l)
+          gidj = group_info(igid,l)
+          if (gidi==gidj) cycle
+       else
+          j = l
+       endif
        if (i==j) cycle
        if (extrap) then
           dx     = xi - (xyzmh_ptmass(1,j) + extrapfac*fsink_old(1,j))
@@ -552,17 +575,35 @@ end subroutine ptmass_boundary_crossing
 !  (called from inside a parallel section)
 !+
 !----------------------------------------------------------------
-subroutine ptmass_drift(nptmass,ckdt,xyzmh_ptmass,vxyz_ptmass)
+subroutine ptmass_drift(nptmass,ckdt,xyzmh_ptmass,vxyz_ptmass,group_info,n_ingroup)
+ use part,only:igarg
  integer, intent(in)    :: nptmass
  real,    intent(in)    :: ckdt
  real,    intent(inout) :: xyzmh_ptmass(nsinkproperties,nptmass)
  real,    intent(inout) :: vxyz_ptmass(3,nptmass)
- integer :: i
+ integer, optional, intent(in)    :: n_ingroup
+ integer, optional, intent(in)    :: group_info(:,:)
+ integer :: i,k,istart_ptmass
+ logical :: woutsub
+
+ if (present(n_ingroup)) then
+    istart_ptmass = n_ingroup + 1
+    woutsub = .true.
+ else
+    istart_ptmass = 1
+    woutsub = .false.
+ endif
 
  !$omp parallel do schedule(static) default(none) &
  !$omp shared(nptmass,ckdt,xyzmh_ptmass,vxyz_ptmass) &
- !$omp private(i)
- do i=1,nptmass
+ !$omp shared(n_ingroup,group_info,woutsub,istart_ptmass) &
+ !$omp private(i,k)
+ do k=istart_ptmass,nptmass
+    if (woutsub) then
+       i = group_info(igarg,k)
+    else
+       i = k
+    endif
     if (xyzmh_ptmass(4,i) > 0.) then
        xyzmh_ptmass(1,i) = xyzmh_ptmass(1,i) + ckdt*vxyz_ptmass(1,i)
        xyzmh_ptmass(2,i) = xyzmh_ptmass(2,i) + ckdt*vxyz_ptmass(2,i)

src/main/step_leapfrog.F90

@@ -105,6 +105,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew)
  use mpiutils,       only:reduceall_mpi
  use part,           only:nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass, &
                           dsdt_ptmass,fsink_old,ibin_wake,dptmass
+ use part,           only:n_group,n_ingroup,n_sing,gtgrad,group_info,nmatrix
  use io_summary,     only:summary_printout,summary_variable,iosumtvi,iowake, &
                           iosumflrp,iosumflrps,iosumflrc
  use boundary_dyn,   only:dynamic_bdy,update_xyzminmax
@@ -247,9 +248,11 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew)
     endif
  else
     if (nptmass > 0 .or. iexternalforce > 0 .or. h2chemistry .or. cooling_in_step .or. idamp > 0) then
+
        call substep(npart,ntypes,nptmass,dtsph,dtextforce,t,xyzh,vxyzu,&
-                                 fext,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,&
-                                 dptmass,fsink_old,nbinmax,ibin_wake)
+                    fext,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,&
+                    dptmass,fsink_old,nbinmax,ibin_wake,gtgrad,group_info, &
+                    nmatrix,n_group,n_ingroup,n_sing)
     else
        call substep_sph(dtsph,npart,xyzh,vxyzu)
     endif