diff --git a/src/main/initial.F90 b/src/main/initial.F90 index dbce5a3dc..b9544f08f 100644 --- a/src/main/initial.F90 +++ b/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 diff --git a/src/main/part.F90 b/src/main/part.F90 index b2ecb91dc..f2553f98e 100644 --- a/src/main/part.F90 +++ b/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) diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 323564b89..a6b7faf26 100644 --- a/src/main/step_leapfrog.F90 +++ b/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 diff --git a/src/main/substepping.F90 b/src/main/substepping.F90 index 1565570da..0fbed5f3b 100644 --- a/src/main/substepping.F90 +++ b/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