Skip to content

Commit

Permalink
Lagrange bubbles fixing PR part 5
Browse files Browse the repository at this point in the history
  • Loading branch information
Diego Vaca committed Dec 10, 2024
1 parent 1fb0954 commit 733fb17
Show file tree
Hide file tree
Showing 11 changed files with 191 additions and 192 deletions.
25 changes: 14 additions & 11 deletions src/common/m_constants.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,25 +46,28 @@ module m_constants
! RKCK constants
integer, parameter :: num_ts_rkck = 6 !< Number of time-stages in the RKCK stepper
! RKCK 4th/5th time stepper coefficients based on Cash J. and Karp A. (1990)
real(kind(0d0)), parameter :: rkck_c1 = 0.0d0, rkck_c2 = 0.2d0, rkck_c3 = 0.3d0, rkck_c4 = 0.6d0, & ! c1 c2 c3 c4 c5 c6
real(kind(0d0)), parameter :: rkck_c1 = 0.0d0, rkck_c2 = 0.2d0, rkck_c3 = 0.3d0, rkck_c4 = 0.6d0, & ! c1 c2 c3 c4 c5 c6
rkck_c5 = 1.0d0, rkck_c6 = 0.875d0
real(kind(0.d0)), dimension(6), parameter :: rkck_coef1 = (/0.2d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0/) ! a21
real(kind(0.d0)), dimension(6), parameter :: rkck_coef2 = (/3.0d0/40.0d0, 9.0d0/40.0d0, 0.0d0, 0.0d0, & ! a31 a32
real(kind(0.d0)), dimension(6), parameter :: rkck_coef1 = (/0.2d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0/) ! a21
real(kind(0.d0)), dimension(6), parameter :: rkck_coef2 = (/3.0d0/40.0d0, 9.0d0/40.0d0, 0.0d0, 0.0d0, & ! a31 a32
0.0d0, 0.0d0/)
real(kind(0.d0)), dimension(6), parameter :: rkck_coef3 = (/0.3d0, -0.9d0, 1.2d0, 0.0d0, 0.0d0, 0.0d0/) ! a41 a42 a43
real(kind(0.d0)), dimension(6), parameter :: rkck_coef4 = (/-11.0d0/54.0d0, 2.5d0, -70.0d0/27.0d0, & ! a51 a52 a53 a54
real(kind(0.d0)), dimension(6), parameter :: rkck_coef3 = (/0.3d0, -0.9d0, 1.2d0, 0.0d0, 0.0d0, 0.0d0/) ! a41 a42 a43
real(kind(0.d0)), dimension(6), parameter :: rkck_coef4 = (/-11.0d0/54.0d0, 2.5d0, -70.0d0/27.0d0, & ! a51 a52 a53 a54
35.d0/27.d0, 0.0d0, 0.0d0/)
real(kind(0.d0)), dimension(6), parameter :: rkck_coef5 = (/1631.0d0/55296.0d0, 175.0d0/512.0d0, & ! a61 a62 a63 a64 a65
575.d0/13824.d0, 44275.d0/110592.d0, 253.d0/4096.d0, 0.0d0/)
real(kind(0.d0)), dimension(6), parameter :: rkck_coef6 = (/37.d0/378.d0, 0.0d0, 250.d0/621.d0, & ! b1 b2 b3 b4 b5 b6
real(kind(0.d0)), dimension(6), parameter :: rkck_coef5 = (/1631.0d0/55296.0d0, 175.0d0/512.0d0, & ! a61 a62 a63 a64 a65
575.d0/13824.d0, 44275.d0/110592.d0, &
253.d0/4096.d0, 0.0d0/)
real(kind(0.d0)), dimension(6), parameter :: rkck_coef6 = (/37.d0/378.d0, 0.0d0, 250.d0/621.d0, & ! b1 b2 b3 b4 b5 b6
125.0d0/594.0d0, 0.0d0, 512.0d0/1771.0d0/)
real(kind(0.d0)), dimension(6), parameter :: rkck_coefE = (/37.d0/378.d0 - 2825.0d0/27648.0d0, 0.0d0, & ! er1 er2 er3 er4 er5 er6 (4th/5th error)
250.d0/621.d0 - 18575.0d0/48384.0d0, 125.0d0/594.0d0 - 13525.0d0/55296.0d0, &
real(kind(0.d0)), dimension(6), parameter :: rkck_coefE = (/37.d0/378.d0 - 2825.0d0/27648.0d0, 0.0d0, & ! er1 er2 er3 er4 er5 er6 (4th/5th error)
250.d0/621.d0 - 18575.0d0/48384.0d0, &
125.0d0/594.0d0 - 13525.0d0/55296.0d0, &
-277.0d0/14336.0d0, 512.0d0/1771.0d0 - 0.25d0/)
! Adaptive rkck constants
real(kind(0d0)), parameter :: verysmall_dt = 1.d-14 !< Very small dt, stop stepper
real(kind(0d0)), parameter :: SAFETY = 0.9d0 !< Next dt will be maximum 0.9*dt if truncation error is above tolerance.
real(kind(0d0)), parameter :: RNDDEC = 1.0d05 !< Need to round the relative truncation error (5th decimal) to avoid the inclusion of random decimals when dividing by a very small number (rkck_tolerance)
real(kind(0d0)), parameter :: RNDDEC = 1.0d05 !< Need to round the relative truncation error (5th decimal) to avoid the inclusion of random decimals
! after dividing calculated tolerance by a very small number (rkck_tolerance)
real(kind(0d0)), parameter :: PSHRNK = -0.25d0 !< Factor to reduce dt when truncation error above tolerance
real(kind(0d0)), parameter :: SHRNKDT = 0.5d0 !< Factor to reduce dt due to negative bubble radius
real(kind(0d0)), parameter :: ERRCON = 1.89d-4 !< Limit to slightly increase dt when truncation error is between ERRCON and 1
Expand Down
78 changes: 37 additions & 41 deletions src/simulation/m_bubbles_EL.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@ contains
end subroutine s_initialize_bubbles_EL_module
!> The purpose of this procedure is to start lagrange bubble parameters applying nondimensionalization if needed
!> The purpose of this procedure is to start the lagrange bubble parameters applying nondimensionalization if needed
subroutine s_start_lagrange_inputs()
integer :: id_bubbles, id_host
Expand Down Expand Up @@ -197,7 +197,7 @@ contains
Web = 1d0/ss
Re_inv = mul0
! Need improvements to accept polytropic gas compression, isothermal and adiabatic thermal models, and
! Need improvement to accept polytropic gas compression, isothermal and adiabatic thermal models, and
! the Gilmore and RP bubble models.
polytropic = .false. ! Forcing no polytropic model
thermal = 3 ! Forcing constant transfer coefficient model based on Preston et al., 2007
Expand Down Expand Up @@ -350,7 +350,7 @@ contains
rhol, gamma, pi_inf, qv, Re)
dynP = 0d0
do i = 1, num_dims
dynP = dynP + 0.5*q_cons_vf(contxe + i)%sf(cell(1), cell(2), cell(3))**2/rhol
dynP = dynP + 0.5d0*q_cons_vf(contxe + i)%sf(cell(1), cell(2), cell(3))**2/rhol
end do
pliq = (q_cons_vf(E_idx)%sf(cell(1), cell(2), cell(3)) - dynP - pi_inf)/gamma
if (pliq < 0) print *, "Negative pressure", proc_rank, &
Expand All @@ -359,14 +359,14 @@ contains
! Initial particle pressure
gas_p(bub_id, 1) = pliq + 2d0*(1d0/Web)/bub_R0(bub_id)
if ((1d0/Web) /= 0d0) then
pcrit = pv - 4d0*(1d0/Web)/(3.d0*sqrt(3d0*gas_p(bub_id, 1)*bub_R0(bub_id)**3/(2d0*(1d0/Web))))
pcrit = pv - 4d0*(1d0/Web)/(3d0*sqrt(3d0*gas_p(bub_id, 1)*bub_R0(bub_id)**3d0/(2d0*(1d0/Web))))
pref = gas_p(bub_id, 1)
else
pcrit = 0d0
end if

! Initial particle mass
volparticle = 4d0/3d0*pi*bub_R0(bub_id)**3 ! volume
volparticle = 4d0/3d0*pi*bub_R0(bub_id)**3d0 ! volume
gas_mv(bub_id, 1) = pv*volparticle*(1d0/(R_v*Tw))*(massflag) ! vapermass
gas_mg(bub_id) = (gas_p(bub_id, 1) - pv*(massflag))*volparticle*(1d0/(R_n*Tw)) ! gasmass
if (gas_mg(bub_id) <= 0d0) then
Expand Down Expand Up @@ -473,9 +473,9 @@ contains
indomain = particle_in_domain(inputvals(1:3))
if (indomain .and. (id > 0)) then
bub_id = bub_id + 1
nBubs = bub_id ! local number of bubbles
lag_id(bub_id, 1) = id !global ID
lag_id(bub_id, 2) = bub_id !local ID
nBubs = bub_id ! local number of bubbles
lag_id(bub_id, 1) = id ! global ID
lag_id(bub_id, 2) = bub_id ! local ID
mtn_pos(bub_id, 1:3, 1) = inputvals(1:3)
mtn_posPrev(bub_id, 1:3, 1) = inputvals(4:6)
mtn_vel(bub_id, 1:3, 1) = inputvals(7:9)
Expand Down Expand Up @@ -540,14 +540,14 @@ contains
fV = intfc_vel(k, 2)
fpb = gas_p(k, 2)
pint = f_cpbw_KM(fR0, fR, fV, fpb)
pint = pint + 0.5d0*fV**2
pint = pint + 0.5d0*fV**2d0
if (lag_params%cluster_type == 2) then
bub_dphidt(k) = (paux - pint) + term2
! Accounting for the potential induced by the bubble averaged over the control volume
! Note that this is based on the incompressible flow assumption near the bubble.
Rb = intfc_rad(k, 2)
term1_fac = 3d0/2d0*(Rb*(Romega**2d0 - Rb**2d0))/(Romega**3d0 - Rb**3d0)
bub_dphidt(k) = bub_dphidt(k)/(1 - term1_fac)
bub_dphidt(k) = bub_dphidt(k)/(1d0 - term1_fac)
end if
end do
end if
Expand Down Expand Up @@ -594,10 +594,10 @@ contains
end do
call s_convert_species_to_mixture_variables_acc(rhol, gamma, pi_inf, qv, myalpha, &
myalpha_rho, Re, cell(1), cell(2), cell(3))
call s_compute_cson_from_pinf(k, q_prim_vf, pinf, cell, rhol, gamma, pi_inf, cson)
call s_compute_cson_from_pinf(q_prim_vf, pinf, cell, rhol, gamma, pi_inf, cson)

! Velocity correction due to massflux
velint = fV - gas_dmvdt(k, stage)/(4d0*pi*fR**2*rhol)
velint = fV - gas_dmvdt(k, stage)/(4d0*pi*fR**2d0*rhol)

! Interphase acceleration and update vars
intfc_dveldt(k, stage) = f_rddot_KM(fpbdt, pinf, pliqint, rhol, fR, velint, fR0, cson)
Expand Down Expand Up @@ -641,20 +641,18 @@ contains
!! @param gamma Liquid specific heat ratio
!! @param pi_inf Liquid stiffness
!! @param cson Calculated speed of sound
subroutine s_compute_cson_from_pinf(bub_id, q_prim_vf, pinf, cell, rhol, gamma, pi_inf, cson)
subroutine s_compute_cson_from_pinf(q_prim_vf, pinf, cell, rhol, gamma, pi_inf, cson)
#ifdef _CRAYFTN
!DIR$ INLINEALWAYS s_compute_cson_from_pinf
#else
!$acc routine seq
#endif
integer, intent(in) :: bub_id
type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
real(kind(0d0)), intent(in) :: pinf, rhol, gamma, pi_inf
integer, dimension(3), intent(in) :: cell
real(kind(0d0)), intent(out) :: cson

real(kind(0d0)) :: E, H
real(kind(0d0)), dimension(3) :: scoord
real(kind(0d0)), dimension(num_dims) :: vel
integer :: i

Expand Down Expand Up @@ -817,7 +815,7 @@ contains
real(kind(0d0)), intent(out), optional :: preterm1, term2, Romega

real(kind(0d0)), dimension(3) :: scoord, psi
real(kind(0d0)) :: dc, vol, aux, dist
real(kind(0d0)) :: dc, vol, aux
real(kind(0d0)) :: volgas, term1, Rbeq, denom
real(kind(0d0)) :: charvol, charpres, charvol2, charpres2
integer, dimension(3) :: cellaux
Expand Down Expand Up @@ -998,22 +996,22 @@ contains

end if

if (lag_params%pressure_corrector .and. present(preterm1)) then
if (lag_params%pressure_corrector) then

!Valid if only one bubble exists per cell
volgas = intfc_rad(bub_id, 2)**3
denom = intfc_rad(bub_id, 2)**2
term1 = bub_dphidt(bub_id)*intfc_rad(bub_id, 2)**2
term2 = intfc_vel(bub_id, 2)*intfc_rad(bub_id, 2)**2
volgas = intfc_rad(bub_id, 2)**3d0
denom = intfc_rad(bub_id, 2)**2d0
term1 = bub_dphidt(bub_id)*intfc_rad(bub_id, 2)**2d0
term2 = intfc_vel(bub_id, 2)*intfc_rad(bub_id, 2)**2d0

Rbeq = volgas**(1d0/3d0) !surrogate bubble radius
aux = dc**3 - Rbeq**3
aux = dc**3d0 - Rbeq**3d0
term2 = term2/denom
term2 = 3d0/2d0*term2**2*Rbeq**3*(1d0 - Rbeq/dc)/aux
preterm1 = 3d0/2d0*Rbeq*(dc**2 - Rbeq**2)/(aux*denom)
term2 = 3d0/2d0*term2**2d0*Rbeq**3d0*(1d0 - Rbeq/dc)/aux
preterm1 = 3d0/2d0*Rbeq*(dc**2d0 - Rbeq**2d0)/(aux*denom)

!Control volume radius
if (present(Romega)) Romega = dc
if (ptype == 2) Romega = dc

! Getting p_inf
if (ptype == 1) then
Expand Down Expand Up @@ -1167,7 +1165,7 @@ contains

lag_largestep = 0d0
remove_id = 0
!$acc parallel loop gang vector default(present) reduction(+: lag_largestep, remove_id) private(k) copyin(RKstep)
!$acc parallel loop gang vector default(present) reduction(+: lag_largestep) reduction(MAX: remove_id) private(k) copyin(RKstep)
do k = 1, nBubs

radiusOld = intfc_rad(k, 2)
Expand All @@ -1188,7 +1186,7 @@ contains
print *, 'Negative bubble radius encountered'
lag_largestep = lag_largestep + 1d0
if (dt < 2d0*verysmall_dt) then
remove_id = k
remove_id = max(remove_id, k)
end if
end if

Expand Down Expand Up @@ -1382,9 +1380,9 @@ contains
!! @param scoord Calculated particle coordinates
subroutine s_locate_cell(pos, cell, scoord)

real(kind(0d0)), dimension(3) :: pos
real(kind(0d0)), dimension(3), optional :: scoord
integer, dimension(3) :: cell
real(kind(0d0)), dimension(3), intent(in) :: pos
real(kind(0d0)), dimension(3), intent(out) :: scoord
integer, dimension(3), intent(inout) :: cell

integer :: i

Expand Down Expand Up @@ -1420,16 +1418,14 @@ contains
! In other words, the coordinate of the center of the cell is x_cc(cell).

!coordinates in computational space
if (present(scoord)) then
scoord(1) = cell(1) + (pos(1) - x_cb(cell(1) - 1))/dx(cell(1))
scoord(2) = cell(2) + (pos(2) - y_cb(cell(2) - 1))/dy(cell(2))
scoord(3) = 0d0
if (p > 0) scoord(3) = cell(3) + (pos(3) - z_cb(cell(3) - 1))/dz(cell(3))
cell(:) = int(scoord(:))
do i = 1, num_dims
if (scoord(i) < 0.0d0) cell(i) = cell(i) - 1
end do
end if
scoord(1) = cell(1) + (pos(1) - x_cb(cell(1) - 1))/dx(cell(1))
scoord(2) = cell(2) + (pos(2) - y_cb(cell(2) - 1))/dy(cell(2))
scoord(3) = 0d0
if (p > 0) scoord(3) = cell(3) + (pos(3) - z_cb(cell(3) - 1))/dz(cell(3))
cell(:) = int(scoord(:))
do i = 1, num_dims
if (scoord(i) < 0.0d0) cell(i) = cell(i) - 1
end do

end subroutine s_locate_cell

Expand Down Expand Up @@ -1847,7 +1843,7 @@ contains

integer :: k

!$acc parallel loop gang vector default(present) reduction(MAX: Rmax_glb) reduction(MAX: Rmin_glb) private(k)
!$acc parallel loop gang vector default(present) reduction(MAX: Rmax_glb) reduction(MIN: Rmin_glb) private(k)
do k = 1, nBubs
Rmax_glb = max(Rmax_glb, intfc_rad(k, 1)/bub_R0(k))
Rmin_glb = min(Rmin_glb, intfc_rad(k, 1)/bub_R0(k))
Expand Down
32 changes: 16 additions & 16 deletions src/simulation/m_bubbles_EL_kernels.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,12 +62,12 @@ contains
!$acc parallel loop gang vector default(present) private(l, s_coord, cell)
do l = 1, nBubs

volpart = 4.0d0/3.0d0*pi*lbk_rad(l, 2)**3
volpart = 4.0d0/3.0d0*pi*lbk_rad(l, 2)**3d0
s_coord(1:3) = lbk_s(l, 1:3, 2)
call s_get_cell(s_coord, cell)

strength_vol = volpart
strength_vel = 4.0d0*pi*lbk_rad(l, 2)**2*lbk_vel(l, 2)
strength_vel = 4.0d0*pi*lbk_rad(l, 2)**2d0*lbk_vel(l, 2)

if (num_dims == 2) then
Vol = dx(cell(1))*dy(cell(2))*lag_params%charwidth
Expand Down Expand Up @@ -128,15 +128,15 @@ contains
do l = 1, nBubs
nodecoord(1:3) = 0
center(1:3) = 0.0d0
volpart = 4.0d0/3.0d0*pi*lbk_rad(l, 2)**3
volpart = 4.0d0/3.0d0*pi*lbk_rad(l, 2)**3d0
s_coord(1:3) = lbk_s(l, 1:3, 2)
center(1:2) = lbk_pos(l, 1:2, 2)
if (p > 0) center(3) = lbk_pos(l, 3, 2)
call s_get_cell(s_coord, cell)
call s_compute_stddsv(cell, volpart, stddsv)

strength_vol = volpart
strength_vel = 4.0d0*pi*lbk_rad(l, 2)**2*lbk_vel(l, 2)
strength_vel = 4.0d0*pi*lbk_rad(l, 2)**2d0*lbk_vel(l, 2)

!$acc loop collapse(3) private(cellaux, nodecoord)
do i = 1, smearGrid
Expand Down Expand Up @@ -218,25 +218,25 @@ contains
real(kind(0.d0)), intent(out) :: func

real(kind(0.d0)) :: distance
real(kind(0.d0)) :: theta, dtheta, L2, dz, Lz2
real(kind(0.d0)) :: theta, dtheta, L2, dzp, Lz2
integer :: Nr, Nr_count

distance = sqrt((center(1) - nodecoord(1))**2 + (center(2) - nodecoord(2))**2 + (center(3) - nodecoord(3))**2)
distance = sqrt((center(1) - nodecoord(1))**2d0 + (center(2) - nodecoord(2))**2d0 + (center(3) - nodecoord(3))**2d0)

if (num_dims == 3) then
!< 3D gaussian function
func = exp(-0.5d0*(distance/stddsv)**2)/(DSQRT(2.0d0*pi)*stddsv)**3
func = exp(-0.5d0*(distance/stddsv)**2d0)/(DSQRT(2.0d0*pi)*stddsv)**3d0
else
if (cyl_coord) then
!< 2D cylindrical function:
! We smear particles in the azimuthal direction for given r
theta = 0d0
Nr = ceiling(2d0*PI*nodecoord(2)/(y_cb(cellaux(2)) - y_cb(cellaux(2) - 1)))
dtheta = 2d0*PI/Nr
Nr = ceiling(2d0*pi*nodecoord(2)/(y_cb(cellaux(2)) - y_cb(cellaux(2) - 1)))
dtheta = 2d0*pi/Nr
L2 = center(2)**2d0 + nodecoord(2)**2d0 - 2d0*center(2)*nodecoord(2)*cos(theta)
distance = DSQRT((center(1) - nodecoord(1))**2d0 + L2)
! Factor 2d0 is for symmetry (upper half of the 2D field (+r) is considered)
func = dtheta/2d0/PI*exp(-0.5d0*(distance/stddsv)**2)/(DSQRT(2.0d0*pi)*stddsv)**3
func = dtheta/2d0/pi*exp(-0.5d0*(distance/stddsv)**2d0)/(DSQRT(2.0d0*pi)*stddsv)**3d0
Nr_count = 0
do while (Nr_count < Nr - 1)
Nr_count = Nr_count + 1
Expand All @@ -246,7 +246,7 @@ contains
distance = DSQRT((center(1) - nodecoord(1))**2d0 + L2)
! nodecoord(2)*dtheta is the azimuthal width of the cell
func = func + &
dtheta/2d0/PI*exp(-0.5d0*(distance/stddsv)**2)/(DSQRT(2.0d0*pi)*stddsv)**(3*(strength_idx + 1))
dtheta/2d0/pi*exp(-0.5d0*(distance/stddsv)**2d0)/(DSQRT(2.0d0*pi)*stddsv)**(3d0*(strength_idx + 1))
end do
else

Expand All @@ -255,16 +255,16 @@ contains
theta = 0d0
Nr = ceiling(lag_params%charwidth/(y_cb(cellaux(2)) - y_cb(cellaux(2) - 1)))
Nr_count = 1 - mapCells
dz = y_cb(cellaux(2) + 1) - y_cb(cellaux(2))
Lz2 = (center(3) - (dz*(5d-1 + Nr_count) - lag_params%charwidth/2d0))**2d0
dzp = y_cb(cellaux(2) + 1) - y_cb(cellaux(2))
Lz2 = (center(3) - (dzp*(5d-1 + Nr_count) - lag_params%charwidth/2d0))**2d0
distance = DSQRT((center(1) - nodecoord(1))**2d0 + (center(2) - nodecoord(2))**2d0 + Lz2)
func = dz/lag_params%charwidth*exp(-0.5d0*(distance/stddsv)**2)/(DSQRT(2.0d0*pi)*stddsv)**3
func = dzp/lag_params%charwidth*exp(-0.5d0*(distance/stddsv)**2d0)/(DSQRT(2.0d0*pi)*stddsv)**3d0
do while (Nr_count < Nr - 1 + (mapCells - 1))
Nr_count = Nr_count + 1
Lz2 = (center(3) - (dz*(5d-1 + Nr_count) - lag_params%charwidth/2d0))**2d0
Lz2 = (center(3) - (dzp*(5d-1 + Nr_count) - lag_params%charwidth/2d0))**2d0
distance = DSQRT((center(1) - nodecoord(1))**2d0 + (center(2) - nodecoord(2))**2d0 + Lz2)
func = func + &
dz/lag_params%charwidth*exp(-0.5d0*(distance/stddsv)**2)/(DSQRT(2.0d0*pi)*stddsv)**(3*(strength_idx + 1))
dzp/lag_params%charwidth*exp(-0.5d0*(distance/stddsv)**2d0)/(DSQRT(2.0d0*pi)*stddsv)**(3d0*(strength_idx + 1))
end do
end if
end if
Expand Down
2 changes: 1 addition & 1 deletion src/simulation/m_global_parameters.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -699,7 +699,7 @@ contains
lag_params%write_bubbles = .false.
lag_params%write_bubbles_stats = .false.
lag_params%nBubs_glb = dflt_int
lag_params%epsilonb = dflt_real
lag_params%epsilonb = 1.0d0
lag_params%charwidth = dflt_real
lag_params%valmaxvoid = dflt_real
lag_params%c0 = dflt_real
Expand Down
Loading

0 comments on commit 733fb17

Please sign in to comment.