|
if (bubbles) then |
|
if (qbmm) then |
|
! advection source |
|
|
|
! bubble sources |
|
|
|
!$acc parallel loop collapse(3) gang vector default(present) |
|
do l = 0, p |
|
do q = 0, n |
|
do i = 0, m |
|
|
|
rhs_vf(alf_idx)%sf(i, q, l) = rhs_vf(alf_idx)%sf(i, q, l) + mom_sp(2)%sf(i, q, l) |
|
j = bubxb |
|
!$acc loop seq |
|
do k = 1, nb |
|
rhs_vf(j)%sf(i, q, l) = & |
|
rhs_vf(j)%sf(i, q, l) + mom_3d(0, 0, k)%sf(i, q, l) |
|
rhs_vf(j + 1)%sf(i, q, l) = & |
|
rhs_vf(j + 1)%sf(i, q, l) + mom_3d(1, 0, k)%sf(i, q, l) |
|
rhs_vf(j + 2)%sf(i, q, l) = & |
|
rhs_vf(j + 2)%sf(i, q, l) + mom_3d(0, 1, k)%sf(i, q, l) |
|
rhs_vf(j + 3)%sf(i, q, l) = & |
|
rhs_vf(j + 3)%sf(i, q, l) + mom_3d(2, 0, k)%sf(i, q, l) |
|
rhs_vf(j + 4)%sf(i, q, l) = & |
|
rhs_vf(j + 4)%sf(i, q, l) + mom_3d(1, 1, k)%sf(i, q, l) |
|
rhs_vf(j + 5)%sf(i, q, l) = & |
|
rhs_vf(j + 5)%sf(i, q, l) + mom_3d(0, 2, k)%sf(i, q, l) |
|
j = j + 6 |
|
end do |
|
end do |
|
end do |
|
end do |
|
else |
|
|
|
!$acc parallel loop collapse(3) gang vector default(present) |
|
do l = 0, p |
|
do k = 0, n |
|
do j = 0, m |
|
divu%sf(j, k, l) = 0d0 |
|
divu%sf(j, k, l) = & |
|
5d-1/dx(j)*(q_prim_qp%vf(contxe + id)%sf(j + 1, k, l) - & |
|
q_prim_qp%vf(contxe + id)%sf(j - 1, k, l)) |
|
|
|
end do |
|
end do |
|
end do |
|
|
|
!$acc parallel loop collapse(3) gang vector default(present) private(Rtmp, Vtmp) |
|
do l = 0, p |
|
do k = 0, n |
|
do j = 0, m |
|
bub_adv_src(j, k, l) = 0d0 |
|
|
|
!$acc loop seq |
|
do q = 1, nb |
|
bub_r_src(j, k, l, q) = 0d0 |
|
bub_v_src(j, k, l, q) = 0d0 |
|
bub_p_src(j, k, l, q) = 0d0 |
|
bub_m_src(j, k, l, q) = 0d0 |
|
end do |
|
end do |
|
end do |
|
end do |
|
|
|
ndirs = 1; if (n > 0) ndirs = 2; if (p > 0) ndirs = 3 |
|
if (id == ndirs) then |
|
|
|
!$acc parallel loop collapse(3) gang vector default(present) private(Rtmp, Vtmp) |
|
do l = 0, p |
|
do k = 0, n |
|
do j = 0, m |
|
|
|
!$acc loop seq |
|
do q = 1, nb |
|
Rtmp(q) = q_prim_qp%vf(rs(q))%sf(j, k, l) |
|
Vtmp(q) = q_prim_qp%vf(vs(q))%sf(j, k, l) |
|
end do |
|
|
|
call s_comp_n_from_prim(q_prim_qp%vf(alf_idx)%sf(j, k, l), & |
|
Rtmp, nbub(j, k, l)) |
|
|
|
call s_quad((Rtmp**2.d0)*Vtmp, R2Vav) |
|
|
|
bub_adv_src(j, k, l) = 4.d0*pi*nbub(j, k, l)*R2Vav |
|
|
|
end do |
|
end do |
|
end do |
|
|
|
!$acc parallel loop collapse(3) gang vector default(present) private(myalpha_rho, myalpha) |
|
do l = 0, p |
|
do k = 0, n |
|
do j = 0, m |
|
!$acc loop seq |
|
do q = 1, nb |
|
|
|
bub_r_src(j, k, l, q) = q_cons_qp%vf(vs(q))%sf(j, k, l) |
|
|
|
!$acc loop seq |
|
do ii = 1, num_fluids |
|
myalpha_rho(ii) = q_cons_qp%vf(ii)%sf(j, k, l) |
|
myalpha(ii) = q_cons_qp%vf(advxb + ii - 1)%sf(j, k, l) |
|
end do |
|
|
|
myRho = 0d0 |
|
n_tait = 0d0 |
|
B_tait = 0d0 |
|
|
|
if (mpp_lim .and. (num_fluids > 2)) then |
|
!$acc loop seq |
|
do ii = 1, num_fluids |
|
myRho = myRho + myalpha_rho(ii) |
|
n_tait = n_tait + myalpha(ii)*gammas(ii) |
|
B_tait = B_tait + myalpha(ii)*pi_infs(ii) |
|
end do |
|
else if (num_fluids > 2) then |
|
!$acc loop seq |
|
do ii = 1, num_fluids - 1 |
|
myRho = myRho + myalpha_rho(ii) |
|
n_tait = n_tait + myalpha(ii)*gammas(ii) |
|
B_tait = B_tait + myalpha(ii)*pi_infs(ii) |
|
end do |
|
else |
|
myRho = myalpha_rho(1) |
|
n_tait = gammas(1) |
|
B_tait = pi_infs(1) |
|
end if |
|
|
|
n_tait = 1.d0/n_tait + 1.d0 !make this the usual little 'gamma' |
|
|
|
myRho = q_prim_qp%vf(1)%sf(j, k, l) |
|
myP = q_prim_qp%vf(E_idx)%sf(j, k, l) |
|
alf = q_prim_qp%vf(alf_idx)%sf(j, k, l) |
|
myR = q_prim_qp%vf(rs(q))%sf(j, k, l) |
|
myV = q_prim_qp%vf(vs(q))%sf(j, k, l) |
|
|
|
if (.not. polytropic) then |
|
pb = q_prim_qp%vf(ps(q))%sf(j, k, l) |
|
mv = q_prim_qp%vf(ms(q))%sf(j, k, l) |
|
call s_bwproperty(pb, q) |
|
vflux = f_vflux(myR, myV, mv, q) |
|
pbdot = f_bpres_dot(vflux, myR, myV, pb, mv, q) |
|
|
|
bub_p_src(j, k, l, q) = nbub(j, k, l)*pbdot |
|
bub_m_src(j, k, l, q) = nbub(j, k, l)*vflux*4.d0*pi*(myR**2.d0) |
|
else |
|
pb = 0d0; mv = 0d0; vflux = 0d0; pbdot = 0d0 |
|
end if |
|
|
|
if (bubble_model == 1) then |
|
! Gilmore bubbles |
|
Cpinf = myP - pref |
|
Cpbw = f_cpbw(R0(q), myR, myV, pb) |
|
myH = f_H(Cpbw, Cpinf, n_tait, B_tait) |
|
c_gas = f_cgas(Cpinf, n_tait, B_tait, myH) |
|
Cpinf_dot = f_cpinfdot(myRho, myP, alf, n_tait, B_tait, bub_adv_src(j, k, l), divu%sf(j, k, l)) |
|
myHdot = f_Hdot(Cpbw, Cpinf, Cpinf_dot, n_tait, B_tait, myR, myV, R0(q), pbdot) |
|
rddot = f_rddot(Cpbw, myR, myV, myH, myHdot, c_gas, n_tait, B_tait) |
|
else if (bubble_model == 2) then |
|
! Keller-Miksis bubbles |
|
Cpinf = myP |
|
Cpbw = f_cpbw_KM(R0(q), myR, myV, pb) |
|
! c_gas = dsqrt( n_tait*(Cpbw+B_tait) / myRho) |
|
c_liquid = DSQRT(n_tait*(myP + B_tait)/(myRho*(1.d0 - alf))) |
|
rddot = f_rddot_KM(pbdot, Cpinf, Cpbw, myRho, myR, myV, R0(q), c_liquid) |
|
else if (bubble_model == 3) then |
|
! Rayleigh-Plesset bubbles |
|
Cpbw = f_cpbw_KM(R0(q), myR, myV, pb) |
|
rddot = f_rddot_RP(myP, myRho, myR, myV, R0(q), Cpbw) |
|
end if |
|
|
|
bub_v_src(j, k, l, q) = nbub(j, k, l)*rddot |
|
|
|
if (alf < 1.d-11) then |
|
bub_adv_src(j, k, l) = 0d0 |
|
bub_r_src(j, k, l, q) = 0d0 |
|
bub_v_src(j, k, l, q) = 0d0 |
|
if (.not. polytropic) then |
|
bub_p_src(j, k, l, q) = 0d0 |
|
bub_m_src(j, k, l, q) = 0d0 |
|
end if |
|
end if |
|
end do |
|
end do |
|
end do |
|
end do |
|
end if |
|
|
|
!$acc parallel loop collapse(3) gang vector default(present) |
|
do l = 0, p |
|
do q = 0, n |
|
do i = 0, m |
|
rhs_vf(alf_idx)%sf(i, q, l) = rhs_vf(alf_idx)%sf(i, q, l) + bub_adv_src(i, q, l) |
|
if (num_fluids > 1) rhs_vf(advxb)%sf(i, q, l) = & |
|
rhs_vf(advxb)%sf(i, q, l) - bub_adv_src(i, q, l) |
|
!$acc loop seq |
|
do k = 1, nb |
|
rhs_vf(rs(k))%sf(i, q, l) = rhs_vf(rs(k))%sf(i, q, l) + bub_r_src(i, q, l, k) |
|
rhs_vf(vs(k))%sf(i, q, l) = rhs_vf(vs(k))%sf(i, q, l) + bub_v_src(i, q, l, k) |
|
if (polytropic .neqv. .true.) then |
|
rhs_vf(ps(k))%sf(i, q, l) = rhs_vf(ps(k))%sf(i, q, l) + bub_p_src(i, q, l, k) |
|
rhs_vf(ms(k))%sf(i, q, l) = rhs_vf(ms(k))%sf(i, q, l) + bub_m_src(i, q, l, k) |
|
end if |
|
end do |
|
end do |
|
end do |
|
end do |
|
end if |
|
end if |
This code used to be in the
m_bubblesmodules, but now it's also ins_compute_rhs? Also,s_compute_rhsis 1863 lines long.MFC/src/simulation/m_rhs.f90
Lines 1113 to 1321 in ecbcb72