(sink gr) fext_ptmass is only defined once in part.f90

src/main/initial.F90

@@ -149,7 +149,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
  use boundary_dyn,     only:dynamic_bdy,init_dynamic_bdy
  use substepping,      only:combine_forces_gr
 #ifdef GR
- use part,             only:metricderivs,metricderivs_ptmass,metrics_ptmass,pxyzu_ptmass
+ use part,             only:metricderivs,metricderivs_ptmass,metrics_ptmass,pxyzu_ptmass,fext_ptmass
  use cons2prim,        only:prim2consall
  use eos,              only:ieos
  use extern_gr,        only:get_grforce_all,get_tmunu_all,get_tmunu_all_exact
@@ -240,9 +240,6 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
 #ifndef GR
  real            :: dtf,fextv(3)
 #endif
-#ifdef GR 
- real            :: fext_ptmass(4,nptmass)
-#endif 
  integer         :: itype,iposinit,ipostmp,ntypes,nderivinit
  logical         :: iexist,read_input_files
  character(len=len(dumpfile)) :: dumpfileold

src/main/part.F90

@@ -193,6 +193,7 @@ module part
  real, allocatable :: pxyzu_ptmass(:,:) !pxyz_ptmass(maxvxyzu,maxgr)
  real, allocatable :: metrics_ptmass(:,:,:,:) !metrics(0:3,0:3,2,maxgr)
  real, allocatable :: metricderivs_ptmass(:,:,:,:) !metricderivs(0:3,0:3,3,maxgr)
+ real, allocatable :: fext_ptmass(:,:)
 !
 !--sink particles
 !
@@ -484,6 +485,7 @@ subroutine allocate_part
  call allocate_array('metrics_ptmass', metrics_ptmass, 4, 4, 2, maxptmassgr)
  call allocate_array('metricderivs_ptmass', metricderivs_ptmass, 4, 4, 3, maxptmassgr)
  call allocate_array('xyzmh_ptmass', xyzmh_ptmass, nsinkproperties, maxptmass)
+ call allocate_array('fext_ptmass', fext_ptmass, 4, maxptmassgr)
  call allocate_array('vxyz_ptmass', vxyz_ptmass, 3, maxptmass)
  call allocate_array('fxyz_ptmass', fxyz_ptmass, 4, maxptmass)
  call allocate_array('fxyz_ptmass_sinksink', fxyz_ptmass_sinksink, 4, maxptmass)
@@ -581,6 +583,7 @@ subroutine deallocate_part
  if (allocated(metrics_ptmass))  deallocate(metrics_ptmass)
  if (allocated(metricderivs_ptmass))  deallocate(metricderivs_ptmass)
  if (allocated(xyzmh_ptmass)) deallocate(xyzmh_ptmass)
+ if (allocated(fext_ptmass))  deallocate(fext_ptmass)
  if (allocated(vxyz_ptmass))  deallocate(vxyz_ptmass)
  if (allocated(fxyz_ptmass))  deallocate(fxyz_ptmass)
  if (allocated(fxyz_ptmass_sinksink)) deallocate(fxyz_ptmass_sinksink)

src/main/step_leapfrog.F90

@@ -114,7 +114,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew)
  use timestep_ind,   only:get_dt,nbinmax,decrease_dtmax,dt_too_small
  use timestep_sts,   only:sts_get_dtau_next,use_sts,ibin_sts,sts_it_n
  use part,           only:ibin,ibin_old,twas,iactive,ibin_wake
- use part,           only:metricderivs,metricderivs_ptmass
+ use part,           only:metricderivs,metricderivs_ptmass,fext_ptmass
  use metric_tools,   only:imet_minkowski,imetric
  use cons2prim,      only:cons2primall,cons2primall_sink
  use extern_gr,      only:get_grforce_all
@@ -141,7 +141,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew)
  real               :: v2mean,hdti
  real               :: dtsinksink
  real               :: fonrmax,poti,dtphi2
- real               :: fext_gas(4,npart),fext_ptmass(4,nptmass)
+ real               :: fext_gas(4,npart)
  integer            :: merge_ij(nptmass)
  integer            :: merge_n
  real(kind=4)       :: t1,t2,tcpu1,tcpu2
@@ -154,9 +154,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew)
 !
 ! set initial quantities
 !
-
  fext_gas = 0.
- fext_ptmass = 0.
  timei  = t
  hdtsph = 0.5*dtsph
  dterr  = bignumber

src/main/substepping.F90

@@ -1077,7 +1077,7 @@ subroutine predict_gr(xyzh,vxyzu,ntypes,pxyzu,fext,npart,nptmass,dt,timei,hdt, &
  use externalforces, only:externalforce,accrete_particles,update_externalforce
  use part,           only:maxphase,isdead_or_accreted,iamboundary,igas,iphase,iamtype,&
                              massoftype,rhoh,ien_type,eos_vars,igamma,itemp,igasP,&
-                             aprmassoftype,apr_level,epot_sinksink
+                             aprmassoftype,apr_level,epot_sinksink,fext_ptmass
  use io_summary,     only:summary_variable,iosumextr,iosumextt,summary_accrete
  use timestep,       only:bignumber,xtol,ptol
  use eos,            only:equationofstate,ieos
@@ -1110,17 +1110,16 @@ subroutine predict_gr(xyzh,vxyzu,ntypes,pxyzu,fext,npart,nptmass,dt,timei,hdt, &
  real :: bin_info(6,nptmass),dsdt_ptmass(3,nptmass)
  real :: dtphi2,dtsinksink,fonrmax
  integer :: merge_ij(nptmass),merge_n
- real    :: fext_gas(4,npart),fext_sinks(4,nptmass)
+ real    :: fext_gas(4,npart)
 
  pmassi  = massoftype(igas)
  itype   = igas
  fext_gas = 0.
- fext_sinks = 0.
 
  !----------------------------------------------
  ! calculate acceleration sink-sink
  !----------------------------------------------
- call get_accel_sink_sink(nptmass,xyzh_ptmass,fext_sinks,epot_sinksink,dtsinksink,&
+ call get_accel_sink_sink(nptmass,xyzh_ptmass,fext_ptmass,epot_sinksink,dtsinksink,&
                          iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass)
  !----------------------------------------------
  ! predictor during substeps for gas particles
@@ -1133,7 +1132,7 @@ subroutine predict_gr(xyzh,vxyzu,ntypes,pxyzu,fext,npart,nptmass,dt,timei,hdt, &
  !$omp shared(nptmass,xyzh_ptmass,vxyz_ptmass,fxyz_ptmass) &
  !$omp shared(maxphase,maxp,eos_vars) &
  !$omp shared(bin_info,dtphi2,poti,fonrmax) &
- !$omp shared(fext_gas,fext_sinks) &
+ !$omp shared(fext_gas,fext_ptmass) &
  !$omp shared(dt,hdt,xtol,ptol,aprmassoftype,apr_level) &
  !$omp shared(ieos,pxyzu,dens,metrics,metricderivs,ien_type) &
  !$omp shared(pxyzu_ptmass,metrics_ptmass,metricderivs_ptmass) &
@@ -1187,7 +1186,7 @@ subroutine predict_gr(xyzh,vxyzu,ntypes,pxyzu,fext,npart,nptmass,dt,timei,hdt, &
           pprev = pxyz
           ! calculate force between sink-gas particles
           call get_accel_sink_gas(nptmass,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),xyzh_ptmass, &
-                                  fext_gas(1,i),fext_gas(2,i),fext_gas(3,i),poti,pmassi,fext_sinks,&
+                                  fext_gas(1,i),fext_gas(2,i),fext_gas(3,i),poti,pmassi,fext_ptmass,&
                                   dsdt_ptmass,fonrmax,dtphi2,bin_info)
 
 
@@ -1257,7 +1256,7 @@ subroutine predict_gr(xyzh,vxyzu,ntypes,pxyzu,fext,npart,nptmass,dt,timei,hdt, &
  enddo predictor
 !$omp end parallel do
 
- call predict_gr_sink(xyzh_ptmass,vxyz_ptmass,ntypes,pxyzu_ptmass,fxyz_ptmass,fext_sinks,nptmass,&
+ call predict_gr_sink(xyzh_ptmass,vxyz_ptmass,ntypes,pxyzu_ptmass,fxyz_ptmass,fext_ptmass,nptmass,&
                       dt,timei,hdt,metrics_ptmass,metricderivs_ptmass,dtextforcenew,pitsmax,perrmax,&
                       xitsmax,xerrmax)
 
@@ -1438,7 +1437,8 @@ subroutine accrete_gr(xyzh,vxyzu,dens,fext,metrics,metricderivs,nlive,naccreted,
  use externalforces, only:externalforce,accrete_particles,update_externalforce
  use options,        only:iexternalforce
  use part,           only:maxphase,isdead_or_accreted,iamboundary,igas,iphase,iamtype,&
-                             massoftype,rhoh,igamma,itemp,igasP,aprmassoftype,apr_level
+                          massoftype,rhoh,igamma,itemp,igasP,aprmassoftype,apr_level,&
+                          fext_ptmass
  use io_summary,     only:summary_variable,iosumextr,iosumextt,summary_accrete
  use timestep,       only:bignumber,C_force
  use eos,            only:equationofstate,ieos
@@ -1470,17 +1470,16 @@ subroutine accrete_gr(xyzh,vxyzu,dens,fext,metrics,metricderivs,nlive,naccreted,
  real :: bin_info(6,nptmass),dsdt_ptmass(3,nptmass)
  real :: dtphi2,dtsinksink,fonrmax,poti
  integer :: merge_ij(nptmass),merge_n
- real    :: fext_gas(4,npart),fext_sinks(4,nptmass)
+ real    :: fext_gas(4,npart)
 
  pmassi     = massoftype(igas)
  itype      = igas
  fext_gas   = 0.
- fext_sinks = 0.
 
  !----------------------------------------------
  ! calculate acceleration sink-sink
  !----------------------------------------------
- call get_accel_sink_sink(nptmass,xyzmh_ptmass,fext_sinks,epot_sinksink,dtsinksink,&
+ call get_accel_sink_sink(nptmass,xyzmh_ptmass,fext_ptmass,epot_sinksink,dtsinksink,&
                          iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass)
 
  dtextforce_min = min(dtextforce_min,C_force*dtsinksink)
@@ -1492,7 +1491,7 @@ subroutine accrete_gr(xyzh,vxyzu,dens,fext,metrics,metricderivs,nlive,naccreted,
  !$omp shared(maxphase,maxp,aprmassoftype,apr_level) &
  !$omp shared(dtsinksink,epot_sinksink,merge_ij,merge_n,dsdt_ptmass) &
  !$omp shared(bin_info,dtphi2,poti,fonrmax) &
- !$omp shared(fext_gas,fext_sinks) &
+ !$omp shared(fext_gas,fext_ptmass) &
  !$omp private(i,accreted) &
  !$omp shared(ieos,dens,pxyzu,iexternalforce,C_force) &
  !$omp private(pri,pondensi,spsoundi,tempi,dtf) &
@@ -1518,7 +1517,7 @@ subroutine accrete_gr(xyzh,vxyzu,dens,fext,metrics,metricderivs,nlive,naccreted,
        pri = pondensi*dens(i)
 
        call get_accel_sink_gas(nptmass,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),xyzmh_ptmass, &
-                                  fext_gas(1,i),fext_gas(2,i),fext_gas(3,i),poti,pmassi,fext_sinks,&
+                                  fext_gas(1,i),fext_gas(2,i),fext_gas(3,i),poti,pmassi,fext_ptmass,&
                                   dsdt_ptmass,fonrmax,dtphi2,bin_info)
 
        call get_grforce(xyzh(:,i),metrics(:,:,:,i),metricderivs(:,:,:,i),vxyzu(1:3,i),dens(i),vxyzu(4,i),pri,fext(1:3,i),dtf)
@@ -1549,7 +1548,7 @@ subroutine accrete_gr(xyzh,vxyzu,dens,fext,metrics,metricderivs,nlive,naccreted,
  enddo accreteloop
  !$omp enddo
  !$omp end parallel
- call accrete_gr_sink(xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fext_sinks,&
+ call accrete_gr_sink(xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fext_ptmass,&
                       metrics_ptmass,metricderivs_ptmass,nlive_sinks,naccreted_sinks,pxyzu_ptmass,&
                       accretedmass,hdt,nptmass,dtextforce_min,timei,dtsinksink)
 end subroutine accrete_gr