Skip to content

Commit

Permalink
add an additionnal array for FSI gradient force computation + few twe…
Browse files Browse the repository at this point in the history 
…aks on FSI routines
  • Loading branch information
@Yrisch
Yrisch committed Apr 10, 2024
1 parent 0c34f8c commit 2912d0a
Show file tree
Hide file tree
Showing 4 changed files with 97 additions and 75 deletions.
4 changes: 3 additions & 1 deletion src/main/part.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
63 changes: 34 additions & 29 deletions src/main/ptmass.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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.
Expand All @@ -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) &
Expand All @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -513,14 +516,15 @@ 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
real, intent(inout) :: 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
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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) &
Expand All @@ -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.
Expand All @@ -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.)
Expand Down
Loading

0 comments on commit 2912d0a

Please sign in to comment.