diff --git a/src/main/part.F90 b/src/main/part.F90 index 652935668..8fd586af6 100644 --- a/src/main/part.F90 +++ b/src/main/part.F90 @@ -206,7 +206,7 @@ module part integer, parameter :: iJ2 = 19 ! 2nd gravity moment due to oblateness real, allocatable :: xyzmh_ptmass(:,:) real, allocatable :: vxyz_ptmass(:,:) - real, allocatable :: fxyz_ptmass(:,:),fxyz_ptmass_sinksink(:,:) + real, allocatable :: fxyz_ptmass(:,:),fxyz_ptmass_sinksink(:,:),fsink_old(:,:) real, allocatable :: dsdt_ptmass(:,:),dsdt_ptmass_sinksink(:,:) integer :: nptmass = 0 ! zero by default real :: epot_sinksink @@ -431,6 +431,7 @@ subroutine allocate_part 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) + call allocate_array('fsink_old', fsink_old, 4, maxptmass) call allocate_array('dsdt_ptmass', dsdt_ptmass, 3, maxptmass) call allocate_array('dsdt_ptmass_sinksink', dsdt_ptmass_sinksink, 3, maxptmass) call allocate_array('poten', poten, maxgrav) @@ -517,6 +518,7 @@ subroutine deallocate_part if (allocated(vxyz_ptmass)) deallocate(vxyz_ptmass) if (allocated(fxyz_ptmass)) deallocate(fxyz_ptmass) if (allocated(fxyz_ptmass_sinksink)) deallocate(fxyz_ptmass_sinksink) + if (allocated(fsink_old)) deallocate(fsink_old) if (allocated(dsdt_ptmass)) deallocate(dsdt_ptmass) if (allocated(dsdt_ptmass_sinksink)) deallocate(dsdt_ptmass_sinksink) if (allocated(poten)) deallocate(poten) diff --git a/src/main/ptmass.F90 b/src/main/ptmass.F90 index 18534adc5..eab603fea 100644 --- a/src/main/ptmass.F90 +++ b/src/main/ptmass.F90 @@ -123,7 +123,8 @@ module ptmass !+ !---------------------------------------------------------------- subroutine get_accel_sink_gas(nptmass,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi,phi, & - pmassi,fxyz_ptmass,dsdt_ptmass,fonrmax,dtphi2,extrapfac) + pmassi,fxyz_ptmass,dsdt_ptmass,fonrmax, & + dtphi2,extrapfac,fsink_old) #ifdef FINVSQRT use fastmath, only:finvsqrt #endif @@ -136,6 +137,7 @@ subroutine get_accel_sink_gas(nptmass,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi,phi, real, intent(in) :: xyzmh_ptmass(nsinkproperties,nptmass) real, optional, intent(in) :: pmassi,extrapfac real, optional, intent(inout) :: fxyz_ptmass(4,nptmass),dsdt_ptmass(3,nptmass) + real, optional, intent(in) :: fsink_old(4,nptmass) real, optional, intent(out) :: fonrmax,dtphi2 real :: ftmpxi,ftmpyi,ftmpzi real :: dx,dy,dz,rr2,ddr,dr3,f1,f2,pmassj,J2,shat(3),Rsink @@ -169,9 +171,9 @@ subroutine get_accel_sink_gas(nptmass,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi,phi, do j=1,nptmass if (extrap)then - dx = xi - (xyzmh_ptmass(1,j) + extrapfac*fxyz_ptmass(1,j)) - dy = yi - (xyzmh_ptmass(2,j) + extrapfac*fxyz_ptmass(2,j)) - dz = zi - (xyzmh_ptmass(3,j) + extrapfac*fxyz_ptmass(3,j)) + dx = xi - (xyzmh_ptmass(1,j) + extrapfac*fsink_old(1,j)) + dy = yi - (xyzmh_ptmass(2,j) + extrapfac*fsink_old(2,j)) + dz = zi - (xyzmh_ptmass(3,j) + extrapfac*fsink_old(3,j)) else dx = xi - xyzmh_ptmass(1,j) dy = yi - xyzmh_ptmass(2,j) @@ -280,7 +282,7 @@ 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) + iexternalforce,ti,merge_ij,merge_n,dsdt_ptmass,extrapfac,fsink_old) #ifdef FINVSQRT use fastmath, only:finvsqrt #endif @@ -297,6 +299,7 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin integer, intent(out) :: merge_ij(:),merge_n real, intent(out) :: dsdt_ptmass(3,nptmass) real, optional, intent(in) :: extrapfac + real, optional, intent(in) :: fsink_old(4,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 @@ -315,7 +318,7 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin merge_ij = 0 if (nptmass <= 1) return ! check if it is a force computed using Omelyan extrapolation method for FSI - if (present(extrapfac)) then + if (present(extrapfac) .and. present(fsink_old)) then extrap = .true. else extrap = .false. @@ -338,7 +341,7 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin !$omp parallel do default(none) & !$omp shared(nptmass,xyzmh_ptmass,fxyz_ptmass,merge_ij,r_merge2,dsdt_ptmass) & !$omp shared(iexternalforce,ti,h_soft_sinksink,potensoft0,hsoft1,hsoft21) & - !$omp shared(extrapfac,extrap) & + !$omp shared(extrapfac,extrap,fsink_old) & !$omp private(i,xi,yi,zi,pmassi,pmassj) & !$omp private(dx,dy,dz,rr2,rr2j,ddr,dr3,f1,f2) & !$omp private(fxi,fyi,fzi,phii,dsx,dsy,dsz) & @@ -349,9 +352,9 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin !$omp reduction(+:phitot,merge_n) do i=1,nptmass if (extrap)then - xi = xyzmh_ptmass(1,i) + extrapfac*fxyz_ptmass(1,i) - yi = xyzmh_ptmass(2,i) + extrapfac*fxyz_ptmass(2,i) - zi = xyzmh_ptmass(3,i) + extrapfac*fxyz_ptmass(3,i) + 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) else xi = xyzmh_ptmass(1,i) yi = xyzmh_ptmass(2,i) @@ -372,9 +375,9 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin do j=1,nptmass if (i==j) cycle if (extrap)then - dx = xi - (xyzmh_ptmass(1,j) + extrapfac*fxyz_ptmass(1,j)) - dy = yi - (xyzmh_ptmass(2,j) + extrapfac*fxyz_ptmass(2,j)) - dz = zi - (xyzmh_ptmass(3,j) + extrapfac*fxyz_ptmass(3,j)) + dx = xi - (xyzmh_ptmass(1,j) + extrapfac*fsink_old(1,j)) + dy = yi - (xyzmh_ptmass(2,j) + extrapfac*fsink_old(2,j)) + dz = zi - (xyzmh_ptmass(3,j) + extrapfac*fsink_old(3,j)) else dx = xi - xyzmh_ptmass(1,j) dy = yi - xyzmh_ptmass(2,j) @@ -513,7 +516,7 @@ end subroutine get_accel_sink_sink !+ !---------------------------------------------------------------- subroutine get_gradf_sink_gas(nptmass,dt,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi, & - pmassi,fxyz_ptmass) + pmassi,fxyz_ptmass,fsink_old) use kernel, only:kernel_softening,kernel_grad_soft,radkern integer, intent(in) :: nptmass real, intent(in) :: xi,yi,zi,hi,dt @@ -521,6 +524,7 @@ subroutine get_gradf_sink_gas(nptmass,dt,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi, & real, intent(in) :: xyzmh_ptmass(nsinkproperties,nptmass) real, intent(in) :: pmassi real, intent(inout) :: fxyz_ptmass(4,nptmass) + real, intent(in) :: fsink_old(4,nptmass) real :: gtmpxi,gtmpyi,gtmpzi real :: dx,dy,dz,rr2,ddr,dr3,g11,g12,g21,g22,pmassj real :: dfx,dfy,dfz,drdotdf @@ -535,9 +539,9 @@ subroutine get_gradf_sink_gas(nptmass,dt,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi, & dx = xi - xyzmh_ptmass(1,j) dy = yi - xyzmh_ptmass(2,j) dz = zi - xyzmh_ptmass(3,j) - dfx = fxi - fxyz_ptmass(1,j) - dfy = fyi - fxyz_ptmass(2,j) - dfz = fzi - fxyz_ptmass(3,j) + dfx = fxi - fsink_old(1,j) + dfy = fyi - fsink_old(2,j) + dfz = fzi - fsink_old(3,j) pmassj = xyzmh_ptmass(4,j) hsoft = xyzmh_ptmass(ihsoft,j) if (hsoft > 0.0) hsoft = max(hsoft,hi) @@ -623,12 +627,13 @@ end subroutine get_gradf_sink_gas ! get gradient correction of the force for FSI integrator (sink-gas) !+ !---------------------------------------------------------------- -subroutine get_gradf_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,dt) +subroutine get_gradf_sink_sink(nptmass,dt,xyzmh_ptmass,fxyz_ptmass,fsink_old) use kernel, only:kernel_softening,kernel_grad_soft,radkern - integer, intent(in) :: nptmass - real, intent(in) :: xyzmh_ptmass(nsinkproperties,nptmass) + integer, intent(in) :: nptmass + real, intent(in) :: xyzmh_ptmass(nsinkproperties,nptmass) real, intent(inout) :: fxyz_ptmass(4,nptmass) - real, intent(in) :: dt + real, intent(in) :: fsink_old(4,nptmass) + real, intent(in) :: dt real :: xi,yi,zi,pmassi,pmassj,fxi,fyi,fzi,gxi,gyi,gzi real :: ddr,dx,dy,dz,dfx,dfy,dfz,drdotdf,rr2,dr3,g1,g2 real :: hsoft1,hsoft21,q2i,qi,psoft,fsoft,gsoft @@ -647,7 +652,7 @@ subroutine get_gradf_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,dt) !--compute N^2 gradf on point mass particles due to each other ! !$omp parallel do default(none) & - !$omp shared(nptmass,xyzmh_ptmass,fxyz_ptmass) & + !$omp shared(nptmass,xyzmh_ptmass,fxyz_ptmass,fsink_old) & !$omp shared(h_soft_sinksink,hsoft21,dt) & !$omp private(i,xi,yi,zi,pmassi,pmassj) & !$omp private(dx,dy,dz,dfx,dfy,dfz,drdotdf,rr2,ddr,dr3,g1,g2) & @@ -659,9 +664,9 @@ subroutine get_gradf_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,dt) zi = xyzmh_ptmass(3,i) pmassi = xyzmh_ptmass(4,i) if (pmassi < 0.) cycle - fxi = fxyz_ptmass(1,i) - fyi = fxyz_ptmass(2,i) - fzi = fxyz_ptmass(3,i) + fxi = fsink_old(1,i) + fyi = fsink_old(2,i) + fzi = fsink_old(3,i) gxi = 0. gyi = 0. gzi = 0. @@ -670,14 +675,14 @@ subroutine get_gradf_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,dt) dx = xi - xyzmh_ptmass(1,j) dy = yi - xyzmh_ptmass(2,j) dz = zi - xyzmh_ptmass(3,j) - dfx = fxi - fxyz_ptmass(1,j) - dfy = fyi - fxyz_ptmass(2,j) - dfz = fzi - fxyz_ptmass(3,j) + dfx = fxi - fsink_old(1,j) + dfy = fyi - fsink_old(2,j) + dfz = fzi - fsink_old(3,j) pmassj = xyzmh_ptmass(4,j) if (pmassj < 0.) cycle rr2 = dx*dx + dy*dy + dz*dz + epsilon(rr2) - drdotdf = dx*dfx + dy*dfy + dz*dfz +epsilon(drdotdf) + drdotdf = dx*dfx + dy*dfy + dz*dfz ddr = 1./sqrt(rr2) gpref = pmassj*((dt**2)/24.) diff --git a/src/main/step_extern.F90 b/src/main/step_extern.F90 index bb54636c0..69fb57ccb 100644 --- a/src/main/step_extern.F90 +++ b/src/main/step_extern.F90 @@ -426,7 +426,8 @@ end subroutine step_extern_sph ! and external forces except ptmass. (4th order scheme) !+ !---------------------------------------------------------------- -subroutine step_extern_FSI(dtextforce,dtsph,time,npart,nptmass,xyzh,vxyzu,fext,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass) +subroutine step_extern_FSI(dtextforce,dtsph,time,npart,nptmass,xyzh,vxyzu,fext, & + xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fsink_old,dsdt_ptmass) use part, only: isdead_or_accreted,igas,massoftype use io, only:iverbose,id,master,iprint,warning,fatal use io_summary, only:summary_variable,iosumextr,iosumextt @@ -434,7 +435,7 @@ subroutine step_extern_FSI(dtextforce,dtsph,time,npart,nptmass,xyzh,vxyzu,fext,x integer, intent(in) :: npart,nptmass real, intent(inout) :: dtextforce real, intent(inout) :: xyzh(:,:),vxyzu(:,:),fext(:,:) - real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:),fxyz_ptmass(4,nptmass),dsdt_ptmass(3,nptmass) + real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:),fxyz_ptmass(4,nptmass),dsdt_ptmass(3,nptmass),fsink_old(4,nptmass) real :: dt,t_end_step,dtextforce_min real :: pmassi,timei logical :: done,last_step @@ -468,9 +469,10 @@ subroutine step_extern_FSI(dtextforce,dtsph,time,npart,nptmass,xyzh,vxyzu,fext,x call drift_4th(ck(1),dt,npart,nptmass,xyzh,xyzmh_ptmass,vxyzu,vxyz_ptmass,dsdt_ptmass) call get_force_4th(nptmass,npart,nsubsteps,pmassi,timei,dtextforce,& xyzh,fext,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass) ! Direct calculation of the force and force gradient - !call get_gradf_extrap_4th(nptmass,npart,nsubsteps,pmassi,timei,dtextforce,& - ! xyzh,fext,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt) ! extrapolation method Omelyan - call get_gradf_4th(nptmass,npart,pmassi,dt,xyzh,fext,xyzmh_ptmass,fxyz_ptmass) + fsink_old=fxyz_ptmass + call get_gradf_extrap_4th(nptmass,npart,nsubsteps,pmassi,timei,dtextforce,& + xyzh,fext,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,fsink_old) ! extrapolation method Omelyan + !call get_gradf_4th(nptmass,npart,pmassi,dt,xyzh,fext,xyzmh_ptmass,fxyz_ptmass,fsink_old) call kick_4th (dk(2),dt,npart,nptmass,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,fext,fxyz_ptmass,dsdt_ptmass) call drift_4th(ck(2),dt,npart,nptmass,xyzh,xyzmh_ptmass,vxyzu,vxyz_ptmass,dsdt_ptmass) call get_force_4th(nptmass,npart,nsubsteps,pmassi,timei,dtextforce,& @@ -518,18 +520,21 @@ subroutine drift_4th(ck,dt,npart,nptmass,xyzh,xyzmh_ptmass,vxyzu,vxyz_ptmass,dsd integer, intent(in) :: npart,nptmass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:),dsdt_ptmass(:,:) + real :: ckdt integer :: i + ckdt = ck*dt + ! Drift gas particles !$omp parallel do default(none) & - !$omp shared(npart,xyzh,vxyzu,dt,ck) & + !$omp shared(npart,xyzh,vxyzu,ckdt) & !$omp private(i) do i=1,npart if (.not.isdead_or_accreted(xyzh(4,i))) then - xyzh(1,i) = xyzh(1,i) + ck*dt*vxyzu(1,i) - xyzh(2,i) = xyzh(2,i) + ck*dt*vxyzu(2,i) - xyzh(3,i) = xyzh(3,i) + ck*dt*vxyzu(3,i) + xyzh(1,i) = xyzh(1,i) + ckdt*vxyzu(1,i) + xyzh(2,i) = xyzh(2,i) + ckdt*vxyzu(2,i) + xyzh(3,i) = xyzh(3,i) + ckdt*vxyzu(3,i) endif enddo !$omp end parallel do @@ -538,16 +543,16 @@ subroutine drift_4th(ck,dt,npart,nptmass,xyzh,xyzmh_ptmass,vxyzu,vxyz_ptmass,dsd if(nptmass>0) then if(id==master) then !$omp parallel do default(none) & - !$omp shared(nptmass,xyzmh_ptmass,vxyz_ptmass,dsdt_ptmass,ck,dt) & + !$omp shared(nptmass,xyzmh_ptmass,vxyz_ptmass,dsdt_ptmass,ckdt) & !$omp private(i) do i=1,nptmass if (xyzmh_ptmass(4,i) > 0.) then - xyzmh_ptmass(1,i) = xyzmh_ptmass(1,i) + ck*dt*vxyz_ptmass(1,i) - xyzmh_ptmass(2,i) = xyzmh_ptmass(2,i) + ck*dt*vxyz_ptmass(2,i) - xyzmh_ptmass(3,i) = xyzmh_ptmass(3,i) + ck*dt*vxyz_ptmass(3,i) - xyzmh_ptmass(ispinx,i) = xyzmh_ptmass(ispinx,i) + ck*dt*dsdt_ptmass(1,i) - xyzmh_ptmass(ispiny,i) = xyzmh_ptmass(ispiny,i) + ck*dt*dsdt_ptmass(2,i) - xyzmh_ptmass(ispinz,i) = xyzmh_ptmass(ispinz,i) + ck*dt*dsdt_ptmass(3,i) + 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) + xyzmh_ptmass(3,i) = xyzmh_ptmass(3,i) + ckdt*vxyz_ptmass(3,i) + xyzmh_ptmass(ispinx,i) = xyzmh_ptmass(ispinx,i) + ckdt*dsdt_ptmass(1,i) + xyzmh_ptmass(ispiny,i) = xyzmh_ptmass(ispiny,i) + ckdt*dsdt_ptmass(2,i) + xyzmh_ptmass(ispinz,i) = xyzmh_ptmass(ispinz,i) + ckdt*dsdt_ptmass(3,i) endif enddo !$omp end parallel do @@ -573,17 +578,20 @@ subroutine kick_4th(dk,dt,npart,nptmass,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,fext real, intent(inout) :: vxyzu(:,:),fext(:,:) real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:),fxyz_ptmass(:,:),dsdt_ptmass(:,:) integer :: i + real :: dkdt + + dkdt = dk*dt ! Kick gas particles !$omp parallel do default(none) & - !$omp shared(npart,fext,xyzh,vxyzu,dt,dk) & + !$omp shared(npart,fext,xyzh,vxyzu,dkdt) & !$omp private(i) do i=1,npart if (.not.isdead_or_accreted(xyzh(4,i))) then - vxyzu(1,i) = vxyzu(1,i) + dk*dt*fext(1,i) - vxyzu(2,i) = vxyzu(2,i) + dk*dt*fext(2,i) - vxyzu(3,i) = vxyzu(3,i) + dk*dt*fext(3,i) + vxyzu(1,i) = vxyzu(1,i) + dkdt*fext(1,i) + vxyzu(2,i) = vxyzu(2,i) + dkdt*fext(2,i) + vxyzu(3,i) = vxyzu(3,i) + dkdt*fext(3,i) endif enddo !$omp end parallel do @@ -592,16 +600,16 @@ subroutine kick_4th(dk,dt,npart,nptmass,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,fext if (nptmass>0) then if(id==master) then !$omp parallel do default(none) & - !$omp shared(nptmass,xyzmh_ptmass,fxyz_ptmass,vxyz_ptmass,dsdt_ptmass,dk,dt) & + !$omp shared(nptmass,xyzmh_ptmass,fxyz_ptmass,vxyz_ptmass,dsdt_ptmass,dkdt) & !$omp private(i) do i=1,nptmass if (xyzmh_ptmass(4,i) > 0.) then - vxyz_ptmass(1,i) = vxyz_ptmass(1,i) + dk*dt*fxyz_ptmass(1,i) - vxyz_ptmass(2,i) = vxyz_ptmass(2,i) + dk*dt*fxyz_ptmass(2,i) - vxyz_ptmass(3,i) = vxyz_ptmass(3,i) + dk*dt*fxyz_ptmass(3,i) - xyzmh_ptmass(ispinx,i) = xyzmh_ptmass(ispinx,i) + dk*dt*dsdt_ptmass(1,i) - xyzmh_ptmass(ispiny,i) = xyzmh_ptmass(ispiny,i) + dk*dt*dsdt_ptmass(2,i) - xyzmh_ptmass(ispinz,i) = xyzmh_ptmass(ispinz,i) + dk*dt*dsdt_ptmass(3,i) + vxyz_ptmass(1,i) = vxyz_ptmass(1,i) + dkdt*fxyz_ptmass(1,i) + vxyz_ptmass(2,i) = vxyz_ptmass(2,i) + dkdt*fxyz_ptmass(2,i) + vxyz_ptmass(3,i) = vxyz_ptmass(3,i) + dkdt*fxyz_ptmass(3,i) + xyzmh_ptmass(ispinx,i) = xyzmh_ptmass(ispinx,i) + dkdt*dsdt_ptmass(1,i) + xyzmh_ptmass(ispiny,i) = xyzmh_ptmass(ispiny,i) + dkdt*dsdt_ptmass(2,i) + xyzmh_ptmass(ispinz,i) = xyzmh_ptmass(ispinz,i) + dkdt*dsdt_ptmass(3,i) endif enddo !$omp end parallel do @@ -636,7 +644,7 @@ subroutine get_force_4th(nptmass,npart,nsubsteps,pmassi,timei,dtextforce,xyzh,fe integer :: i real :: dtf,dtextforcenew,dtsinkgas,dtphi2,fonrmax real :: fextx,fexty,fextz - real :: fonrmaxi,phii,dtphi2i,extrapfac + real :: fonrmaxi,phii,dtphi2i dtextforcenew = bignumber dtsinkgas = bignumber @@ -713,7 +721,7 @@ end subroutine get_force_4th !---------------------------------------------------------------- -subroutine get_gradf_4th(nptmass,npart,pmassi,dt,xyzh,fext,xyzmh_ptmass,fxyz_ptmass) +subroutine get_gradf_4th(nptmass,npart,pmassi,dt,xyzh,fext,xyzmh_ptmass,fxyz_ptmass,fsink_old) use dim, only:maxptmass use ptmass, only:get_gradf_sink_gas,get_gradf_sink_sink use mpiutils, only:reduce_in_place_mpi @@ -721,6 +729,7 @@ subroutine get_gradf_4th(nptmass,npart,pmassi,dt,xyzh,fext,xyzmh_ptmass,fxyz_ptm integer, intent(in) :: nptmass,npart real, intent(inout) :: xyzh(:,:),fext(:,:) real, intent(inout) :: xyzmh_ptmass(:,:),fxyz_ptmass(4,nptmass) + real, intent(in) :: fsink_old(4,nptmass) real, intent(inout) :: dt real, intent(in) :: pmassi real :: fextx,fexty,fextz @@ -729,23 +738,23 @@ subroutine get_gradf_4th(nptmass,npart,pmassi,dt,xyzh,fext,xyzmh_ptmass,fxyz_ptm if (nptmass>0) then if(id==master) then - call get_gradf_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,dt) + call get_gradf_sink_sink(nptmass,dt,xyzmh_ptmass,fxyz_ptmass,fsink_old) else fxyz_ptmass(:,:) = 0. endif endif !$omp parallel default(none) & - !$omp shared(npart,nptmass,xyzh,xyzmh_ptmass,fext,dt,pmassi) & + !$omp shared(npart,nptmass,xyzh,xyzmh_ptmass,fext,dt,pmassi,fsink_old) & !$omp private(fextx,fexty,fextz) & !$omp reduction(+:fxyz_ptmass) !$omp do do i=1,npart - fextx = 0. - fexty = 0. - fextz = 0. + fextx = fext(1,i) + fexty = fext(2,i) + fextz = fext(3,i) call get_gradf_sink_gas(nptmass,dt,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),& - xyzmh_ptmass,fextx,fexty,fextz,pmassi,fxyz_ptmass) + xyzmh_ptmass,fextx,fexty,fextz,pmassi,fxyz_ptmass,fsink_old) fext(1,i) = fext(1,i)+ fextx fext(2,i) = fext(2,i)+ fexty fext(3,i) = fext(3,i)+ fextz @@ -768,17 +777,18 @@ end subroutine get_gradf_4th subroutine get_gradf_extrap_4th(nptmass,npart,nsubsteps,pmassi,timei,dtextforce,xyzh, & - fext,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt) + fext,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,fsink_old) use options, only:iexternalforce use dim, only:maxptmass use part, only:epot_sinksink - use io, only:iverbose,master,id,iprint,warning,fatal + use io, only:master,id use ptmass, only:get_accel_sink_gas,get_accel_sink_sink,merge_sinks use timestep, only:bignumber use mpiutils, only:reduce_in_place_mpi integer, intent(in) :: nptmass,npart,nsubsteps real, intent(inout) :: xyzh(:,:),fext(:,:) real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:),fxyz_ptmass(4,nptmass),dsdt_ptmass(3,nptmass) + real, intent(in) :: fsink_old(4,nptmass) real, intent(inout) :: dtextforce real, intent(in) :: timei,pmassi,dt integer :: merge_ij(nptmass) @@ -788,14 +798,18 @@ subroutine get_gradf_extrap_4th(nptmass,npart,nsubsteps,pmassi,timei,dtextforce, real :: fextx,fexty,fextz,xi,yi,zi real :: fonrmaxi,phii,dtphi2i,extrapfac + extrapfac = (1/24.)*dt**2 + if (nptmass>0) then if (id==master) then call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,& - dtf,iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass,extrapfac) + dtf,iexternalforce,timei,merge_ij,merge_n, & + dsdt_ptmass,extrapfac,fsink_old) if (merge_n > 0) then call merge_sinks(timei,nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,merge_ij) call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,& - dtf,iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass,extrapfac) + dtf,iexternalforce,timei,merge_ij,merge_n, & + dsdt_ptmass,extrapfac,fsink_old) endif else fxyz_ptmass(:,:) = 0. @@ -805,7 +819,7 @@ subroutine get_gradf_extrap_4th(nptmass,npart,nsubsteps,pmassi,timei,dtextforce, !$omp parallel default(none) & - !$omp shared(npart,nptmass,xyzh,xyzmh_ptmass,fext,extrapfac) & + !$omp shared(npart,nptmass,xyzh,xyzmh_ptmass,fext,extrapfac,fsink_old) & !$omp private(fextx,fexty,fextz,xi,yi,zi) & !$omp private(fonrmaxi,dtphi2i,phii,pmassi,dtf) & !$omp reduction(+:fxyz_ptmass,dsdt_ptmass) @@ -815,10 +829,11 @@ subroutine get_gradf_extrap_4th(nptmass,npart,nsubsteps,pmassi,timei,dtextforce, fexty = 0. fextz = 0. xi = xyzh(1,i) + extrapfac*fext(1,i) - xi = xyzh(1,i) + extrapfac*fext(1,i) - xi = xyzh(1,i) + extrapfac*fext(1,i) + yi = xyzh(2,i) + extrapfac*fext(2,i) + zi = xyzh(3,i) + extrapfac*fext(3,i) call get_accel_sink_gas(nptmass,xi,yi,zi,xyzh(4,i),xyzmh_ptmass,& - fextx,fexty,fextz,phii,pmassi,fxyz_ptmass,dsdt_ptmass,fonrmaxi,dtphi2i,extrapfac) + fextx,fexty,fextz,phii,pmassi,fxyz_ptmass, & + dsdt_ptmass,fonrmaxi,dtphi2i,extrapfac,fsink_old) fext(1,i) = fextx fext(2,i) = fexty fext(3,i) = fextz diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 6edbd45a0..f643d4c5e 100644 --- a/src/main/step_leapfrog.F90 +++ b/src/main/step_leapfrog.F90 @@ -104,7 +104,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) use deriv, only:derivs use timestep, only:dterr,bignumber,tolv use mpiutils, only:reduceall_mpi - use part, only:nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,ibin_wake + use part, only:nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,fsink_old,ibin_wake use io_summary, only:summary_printout,summary_variable,iosumtvi,iowake, & iosumflrp,iosumflrps,iosumflrc use boundary_dyn, only:dynamic_bdy,update_xyzminmax @@ -250,7 +250,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) if (nptmass > 0 .or. iexternalforce > 0 .or. h2chemistry .or. cooling_in_step .or. idamp > 0) then if (use_fourthorder) then call step_extern_FSI(dtextforce,dtsph,t,npart,nptmass,xyzh,vxyzu,fext, & - xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass) + xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fsink_old,dsdt_ptmass) else call step_extern_lf(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,fext,fxyzu,t, & nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,nbinmax,ibin_wake)