Skip to content

Commit

Permalink
QBMM + Non-polytropic Fix (#391)
Browse files Browse the repository at this point in the history 
Co-authored-by: Anand <anand@lawn-128-61-11-198.lawn.gatech.edu>
  • Loading branch information
@anandrdbz
anandrdbz and Anand authored Apr 9, 2024
1 parent e9549c1 commit e3c4a06
Show file tree
Hide file tree
Showing 4 changed files with 190 additions and 165 deletions.
121 changes: 73 additions & 48 deletions src/simulation/m_qbmm.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -432,13 +432,13 @@ contains
integer :: idir
integer :: i, j, k, l, q

real(kind(0d0)) :: nb_q, nb_dot, R, R2, nR, nR2, nR_dot, nR2_dot, var
real(kind(0d0)) :: nb_q, nb_dot, R, R2, nR, nR2, nR_dot, nR2_dot, var, AX

if (idir == 1) then

!Non-polytropic qbmm needs to account for change in bubble radius due to a change in nb
if (.not. polytropic) then
!$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var)
!$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var, AX)
do i = 1, nb
do q = 1, nnode
do l = 0, p
Expand All @@ -451,29 +451,35 @@ contains
R = q_prim_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
R2 = q_prim_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)

nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j - 1, k, l) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l)
nR_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j - 1, k, l) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
nR2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j - 1, k, l) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)

rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dx(j)*R*nb_q**2)* &
(nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i))

if (R2 - R**2d0 > 0d0) then
var = R2 - R**2d0
else
var = verysmall
end if

if (q <= 2) then
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dx(j)*R*nb_q**2*dsqrt(var))* &
AX = R - dsqrt(var)
else
AX = R + dsqrt(var)
end if

nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j - 1, k, l) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l)
nR_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j - 1, k, l) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
nR2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j - 1, k, l) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)

rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dx(j)*AX*nb_q**2)* &
(nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i))

if (q <= 2) then
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dx(j)*AX*nb_q**2*dsqrt(var)*2d0)* &
(nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i))
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dx(j)*R*nb_q**2*dsqrt(var))* &
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dx(j)*AX*nb_q**2*dsqrt(var)*2d0)* &
(-2d0*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i))

else
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dx(j)*R*nb_q**2*dsqrt(var))* &
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dx(j)*AX*nb_q**2*dsqrt(var)*2d0)* &
(nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i))
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dx(j)*R*nb_q**2*dsqrt(var))* &
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dx(j)*AX*nb_q**2*dsqrt(var)*2d0)* &
(-2d0*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i))
end if

Expand Down Expand Up @@ -516,7 +522,7 @@ contains

!Non-polytropic qbmm needs to account for change in bubble radius due to a change in nb
if (.not. polytropic) then
!$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var)
!$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var, AX)
do i = 1, nb
do q = 1, nnode
do l = 0, p
Expand All @@ -529,29 +535,35 @@ contains
R = q_prim_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
R2 = q_prim_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)

nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k - 1, l) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l)
nR_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k - 1, l) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
nR2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k - 1, l) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)

rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dy(k)*R*nb_q**2)* &
(nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i))

if (R2 - R**2d0 > 0d0) then
var = R2 - R**2d0
else
var = verysmall
end if

if (q <= 2) then
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dy(k)*R*nb_q**2*dsqrt(var))* &
AX = R - dsqrt(var)
else
AX = R + dsqrt(var)
end if

nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k - 1, l) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l)
nR_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k - 1, l) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
nR2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k - 1, l) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)

rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dy(k)*AX*nb_q**2)* &
(nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i))

if (q <= 2) then
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dy(k)*AX*nb_q**2*dsqrt(var)*2d0)* &
(nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i))
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dy(k)*R*nb_q**2*dsqrt(var))* &
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dy(k)*AX*nb_q**2*dsqrt(var)*2d0)* &
(-2d0*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i))

else
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dy(k)*R*nb_q**2*dsqrt(var))* &
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dy(k)*AX*nb_q**2*dsqrt(var)*2d0)* &
(nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i))
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dy(k)*R*nb_q**2*dsqrt(var))* &
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dy(k)*AX*nb_q**2*dsqrt(var)*2d0)* &
(-2d0*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i))
end if

Expand All @@ -567,7 +579,7 @@ contains
if (.not. polytropic) then
if (grid_geometry == 3) then
!Non-polytropic qbmm needs to account for change in bubble radius due to a change in nb
!$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var)
!$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var, AX)
do i = 1, nb
do q = 1, nnode
do l = 0, p
Expand All @@ -580,28 +592,35 @@ contains
R = q_prim_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
R2 = q_prim_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)

nb_dot = q_prim_vf(contxe + idir)%sf(j, k, l)*(flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l))
nR_dot = q_prim_vf(contxe + idir)%sf(j, k, l)*(flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l))
nR2_dot = q_prim_vf(contxe + idir)%sf(j, k, l)*(flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l))

rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dz(l)*y_cc(k)*R*nb_q**2)* &
(nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i))
if (R2 - R**2d0 > 0d0) then
var = R2 - R**2d0
else
var = verysmall
end if

if (q <= 2) then
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dz(l)*y_cc(k)*R*nb_q**2*dsqrt(var))* &
AX = R - dsqrt(var)
else
AX = R + dsqrt(var)
end if

nb_dot = q_prim_vf(contxe + idir)%sf(j, k, l)*(flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l))
nR_dot = q_prim_vf(contxe + idir)%sf(j, k, l)*(flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l))
nR2_dot = q_prim_vf(contxe + idir)%sf(j, k, l)*(flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l))

rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dz(l)*y_cc(k)*AX*nb_q**2)* &
(nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i))

if (q <= 2) then
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dz(l)*y_cc(k)*AX*nb_q**2*dsqrt(var)*2d0)* &
(nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i))
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dz(l)*y_cc(k)*R*nb_q**2*dsqrt(var))* &
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dz(l)*y_cc(k)*AX*nb_q**2*dsqrt(var)*2d0)* &
(-2d0*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i))

else
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dz(l)*y_cc(k)*R*nb_q**2*dsqrt(var))* &
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dz(l)*y_cc(k)*AX*nb_q**2*dsqrt(var)*2d0)* &
(nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i))
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dz(l)*y_cc(k)*R*nb_q**2*dsqrt(var))* &
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dz(l)*y_cc(k)*AX*nb_q**2*dsqrt(var)*2d0)* &
(-2d0*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i))
end if
end do
Expand All @@ -611,7 +630,7 @@ contains
end do
else
!Non-polytropic qbmm needs to account for change in bubble radius due to a change in nb
!$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var)
!$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var, AX)
do i = 1, nb
do q = 1, nnode
do l = 0, p
Expand All @@ -624,29 +643,35 @@ contains
R = q_prim_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
R2 = q_prim_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)

nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l)
nR_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
nR2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)

rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dz(l)*R*nb_q**2)* &
(nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i))

if (R2 - R**2d0 > 0d0) then
var = R2 - R**2d0
else
var = verysmall
end if

if (q <= 2) then
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dz(l)*R*nb_q**2*dsqrt(var))* &
AX = R - dsqrt(var)
else
AX = R + dsqrt(var)
end if

nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l)
nR_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
nR2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)

rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dz(l)*AX*nb_q**2)* &
(nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i))

if (q <= 2) then
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dz(l)*AX*nb_q**2*dsqrt(var)*2d0)* &
(nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i))
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dz(l)*R*nb_q**2*dsqrt(var))* &
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0*gam/(dz(l)*AX*nb_q**2*dsqrt(var)*2d0)* &
(-2d0*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i))

else
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dz(l)*R*nb_q**2*dsqrt(var))* &
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dz(l)*AX*nb_q**2*dsqrt(var)*2d0)* &
(nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i))
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dz(l)*R*nb_q**2*dsqrt(var))* &
rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0*gam/(dz(l)*AX*nb_q**2*dsqrt(var)*2d0)* &
(-2d0*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i))
end if

Expand Down
Loading

0 comments on commit e3c4a06

Please sign in to comment.