Skip to content

Commit e2c7f7c

Browse files
author
Diego Vaca
committed
p_inf from subroutine to function
1 parent 654bc3c commit e2c7f7c

File tree

1 file changed

+30
-35
lines changed

1 file changed

+30
-35
lines changed

‎src/simulation/m_bubbles_EL.fpp‎

Lines changed: 30 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -518,10 +518,10 @@ contains
518518
integer, intent(in) :: stage
519519

520520
real(wp) :: vaporflux, pliqint
521-
real(wp) :: preterm1, term2, paux, pint, Romega, term1_fac, Rb
521+
real(wp) :: preterm1, term2, pint, Romega, term1_fac, Rb
522522
real(wp) :: conc_v, R_m, gamma_m, fpb, fmass_n, fmass_v
523523
real(wp) :: fR, fV, fbeta_c, fbeta_t, fR0, fpbdt
524-
real(wp) :: pinf, aux1, aux2, velint, cson, rhol
524+
real(wp) :: pinf, velint, cson, rhol
525525
real(wp) :: gamma, pi_inf, qv
526526
real(wp), dimension(contxe) :: myalpha_rho, myalpha
527527
real(wp), dimension(2) :: Re
@@ -538,19 +538,18 @@ contains
538538
! Calculate velocity potentials (valid for one bubble per cell)
539539
!$acc parallel loop gang vector default(present) private(k, cell)
540540
do k = 1, nBubs
541-
call s_get_pinf(k, q_prim_vf, 2, paux, cell, preterm1, term2, Romega)
541+
pinf = f_pinf(k, q_prim_vf, 2, cell, preterm1, term2, Romega)
542542
fR0 = bub_R0(k)
543543
fR = intfc_rad(k, 2)
544544
fV = intfc_vel(k, 2)
545545
fpb = gas_p(k, 2)
546546
pint = f_cpbw_KM(fR0, fR, fV, fpb)
547547
pint = pint + 0.5_wp*fV**2._wp
548548
if (lag_params%cluster_type == 2) then
549-
bub_dphidt(k) = (paux - pint) + term2
549+
bub_dphidt(k) = (pinf - pint) + term2
550550
! Accounting for the potential induced by the bubble averaged over the control volume
551551
! Note that this is based on the incompressible flow assumption near the bubble.
552-
Rb = intfc_rad(k, 2)
553-
term1_fac = 3._wp/2._wp*(Rb*(Romega**2._wp - Rb**2._wp))/(Romega**3._wp - Rb**3._wp)
552+
term1_fac = 3._wp/2._wp*(fR*(Romega**2._wp - fR**2._wp))/(Romega**3._wp - fR**3._wp)
554553
bub_dphidt(k) = bub_dphidt(k)/(1._wp - term1_fac)
555554
end if
556555
end do
@@ -587,7 +586,7 @@ contains
587586
pliqint = f_cpbw_KM(fR0, fR, fV, fpb)
588587

589588
! Obtaining driving pressure
590-
call s_get_pinf(k, q_prim_vf, 1, pinf, cell, aux1, aux2)
589+
pinf = f_pinf(k, q_prim_vf, 1, cell, preterm1, term2, Romega)
591590

592591
! Obtain liquid density and computing speed of sound from pinf
593592
!$acc loop seq
@@ -804,32 +803,28 @@ contains
804803
!! @param bub_id Particle identifier
805804
!! @param q_prim_vf Primitive variables
806805
!! @param ptype 1: p at infinity, 2: averaged P at the bubble location
807-
!! @param f_pinfl Driving pressure
808806
!! @param cell Bubble cell
809807
!! @param Romega Control volume radius
810-
subroutine s_get_pinf(bub_id, q_prim_vf, ptype, f_pinfl, cell, preterm1, term2, Romega)
811-
#ifdef _CRAYFTN
812-
!DIR$ INLINEALWAYS s_get_pinf
813-
#else
808+
function f_pinf(bub_id, q_prim_vf, ptype, cell, preterm1, term2, Romega)
814809
!$acc routine seq
815-
#endif
816810
integer, intent(in) :: bub_id, ptype
817811
type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
818-
real(wp), intent(out) :: f_pinfl
819812
integer, dimension(3), intent(out) :: cell
820-
real(wp), intent(out), optional :: preterm1, term2, Romega
813+
real(wp), intent(out) :: preterm1, term2, Romega
814+
815+
real(wp) :: f_pinf
821816

822817
real(wp), dimension(3) :: scoord, psi
823818
real(wp) :: dc, vol, aux
824819
real(wp) :: volgas, term1, Rbeq, denom
825820
real(wp) :: charvol, charpres, charvol2, charpres2
826821
integer, dimension(3) :: cellaux
827822
integer :: i, j, k
828-
integer :: smearGrid, smearGrid_zDir
823+
integer :: smearGrid, smearGridz
829824
logical :: celloutside
830825

831826
scoord = mtn_s(bub_id, 1:3, 2)
832-
f_pinfl = 0._wp
827+
f_pinf = 0._wp
833828

834829
!< Find current bubble cell
835830
cell(:) = int(scoord(:))
@@ -892,19 +887,19 @@ contains
892887

893888
!< Perform bilinear interpolation
894889
if (p == 0) then !2D
895-
f_pinfl = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))
896-
f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))
897-
f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)
898-
f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)
890+
f_pinf = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))
891+
f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))
892+
f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)
893+
f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)
899894
else !3D
900-
f_pinfl = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))*(1._wp - psi(3))
901-
f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))*(1._wp - psi(3))
902-
f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)*(1._wp - psi(3))
903-
f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)*(1._wp - psi(3))
904-
f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 1)*(1._wp - psi(1))*(1._wp - psi(2))*psi(3)
905-
f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3) + 1)*psi(1)*(1._wp - psi(2))*psi(3)
906-
f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3) + 1)*psi(1)*psi(2)*psi(3)
907-
f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3) + 1)*(1._wp - psi(1))*psi(2)*psi(3)
895+
f_pinf = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))*(1._wp - psi(3))
896+
f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))*(1._wp - psi(3))
897+
f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)*(1._wp - psi(3))
898+
f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)*(1._wp - psi(3))
899+
f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 1)*(1._wp - psi(1))*(1._wp - psi(2))*psi(3)
900+
f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3) + 1)*psi(1)*(1._wp - psi(2))*psi(3)
901+
f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3) + 1)*psi(1)*psi(2)*psi(3)
902+
f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3) + 1)*(1._wp - psi(1))*psi(2)*psi(3)
908903
end if
909904

910905
!R_Omega
@@ -914,8 +909,8 @@ contains
914909
! Bubble dynamic closure from Maeda and Colonius (2018)
915910
! Include the cell that contains the bubble (mapCells+1+mapCells)
916911
smearGrid = 2*mapCells + 1
917-
smearGrid_zDir = smearGrid
918-
if (p == 0) smearGrid_zDir = 1
912+
smearGridz = smearGrid
913+
if (p == 0) smearGridz = 1
919914

920915
charvol = 0._wp
921916
charpres = 0._wp
@@ -928,7 +923,7 @@ contains
928923
!$acc loop seq
929924
do j = 1, smearGrid
930925
!$acc loop seq
931-
do k = 1, smearGrid_zDir
926+
do k = 1, smearGridz
932927
cellaux(1) = cell(1) + i - (mapCells + 1)
933928
cellaux(2) = cell(2) + j - (mapCells + 1)
934929
cellaux(3) = cell(3) + k - (mapCells + 1)
@@ -984,7 +979,7 @@ contains
984979
end do
985980
end do
986981

987-
f_pinfl = charpres2/charvol2
982+
f_pinf = charpres2/charvol2
988983
vol = charvol
989984
dc = (3._wp*abs(vol)/(4._wp*pi))**(1._wp/3._wp)
990985

@@ -1009,12 +1004,12 @@ contains
10091004

10101005
! Getting p_inf
10111006
if (ptype == 1) then
1012-
f_pinfl = f_pinfl + preterm1*term1 + term2
1007+
f_pinf = f_pinf + preterm1*term1 + term2
10131008
end if
10141009

10151010
end if
10161011

1017-
end subroutine s_get_pinf
1012+
end function f_pinf
10181013

10191014
!> This subroutine updates the Lagrange variables using the tvd RK time steppers.
10201015
!! The time derivative of the bubble variables must be stored at every stage to avoid precision errors.

0 commit comments

Comments
 (0)