Skip to content

Commit d5f6ccd

Browse files
author
Anand Radhakrishnan
committed
Major Refactoring Done, All 1D and 2D tests pass on GPUs
1 parent cb2d236 commit d5f6ccd

File tree

12 files changed

+375
-1193
lines changed

12 files changed

+375
-1193
lines changed

‎src/simulation/m_bubbles.f90‎

Lines changed: 190 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,11 @@ module m_bubbles
2222
real(kind(0.d0)) :: chi_vw !< Bubble wall properties (Ando 2010)
2323
real(kind(0.d0)) :: k_mw !< Bubble wall properties (Ando 2010)
2424
real(kind(0.d0)) :: rho_mw !< Bubble wall properties (Ando 2010)
25-
!$acc declare create(chi_vw, k_mw, rho_mw)
25+
!$acc declare create(chi_vw, k_mw, rho_mw)
26+
27+
integer, allocatable, dimension(:) :: rs, vs, ms, ps
28+
!$acc declare create(rs, vs, ms, ps)
29+
2630

2731
contains
2832

@@ -37,20 +41,51 @@ module m_bubbles
3741
!! @param bub_v_src Bubble velocity equation source
3842
!! @param bub_p_src Bubble pressure equation source
3943
!! @param bub_m_src Bubble mass equation source
40-
subroutine s_compute_bubble_source(idir, q_prim_vf, q_cons_vf, mydivu, &
41-
bub_adv_src, bub_r_src, bub_v_src, bub_p_src, bub_m_src)
44+
45+
subroutine s_initialize_bubbles_module()
46+
47+
integer :: i, j, k, l, q
48+
49+
allocate (rs(1:nb))
50+
allocate (vs(1:nb))
51+
if (.not. polytropic) then
52+
allocate (ps(1:nb))
53+
allocate (ms(1:nb))
54+
end if
55+
56+
do l = 1, nb
57+
rs(l) = bub_idx%rs(l)
58+
vs(l) = bub_idx%vs(l)
59+
if (.not. polytropic) then
60+
ps(l) = bub_idx%ps(l)
61+
ms(l) = bub_idx%ms(l)
62+
end if
63+
end do
64+
65+
!$acc update device(rs, vs)
66+
if (.not. polytropic) then
67+
!$acc update device(ps, ms)
68+
end if
69+
70+
end subroutine
71+
72+
73+
subroutine s_compute_bubble_source(bub_adv_src, bub_r_src, bub_v_src, bub_p_src, bub_m_src, divu, nbub, &
74+
q_cons_vf, q_prim_vf, t_step, id, rhs_vf)
4275

4376
type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf, q_cons_vf
44-
type(scalar_field), intent(IN) :: mydivu
45-
integer, intent(IN) :: idir
77+
type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf
78+
type(scalar_field), intent(IN) :: divu
79+
real(kind(0d0)), dimension(0:m, 0:n, 0:p), intent(INOUT) :: nbub
80+
integer, intent(IN) :: t_step, id
4681

4782
real(kind(0d0)), dimension(0:m, 0:n, 0:p), intent(INOUT) :: bub_adv_src
48-
real(kind(0d0)), dimension(0:m, 0:n, 0:p, 1:nb), intent(INOUT) :: bub_r_src, &
83+
real(kind(0d0)), dimension(0:m, 0:n, 0:p, 1:nb ), intent(INOUT) :: bub_r_src, &
4984
bub_v_src, &
5085
bub_p_src, &
5186
bub_m_src
5287

53-
real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: nbub !< Bubble number density
88+
!< Bubble number density
5489

5590
real(kind(0d0)) :: tmp1, tmp2, tmp3, tmp4, &
5691
c_gas, c_liquid, &
@@ -63,97 +98,172 @@ subroutine s_compute_bubble_source(idir, q_prim_vf, q_cons_vf, mydivu, &
6398

6499
real(kind(0d0)), dimension(nb) :: Rtmp, Vtmp
65100
real(kind(0d0)) :: myR, myV, alf, myP, myRho, R2Vav
101+
real(kind(0d0)), dimension(num_fluids) :: myalpha, myalpha_rho
66102

67103
real(kind(0d0)), dimension(2) :: Re !< Reynolds number
68104

69-
integer :: j, k, l, q, s !< Loop variables
105+
integer :: i, j, k, l, q, ii !< Loop variables
70106
integer :: ndirs !< Number of coordinate directions
71107

72-
ndirs = 1; if (n > 0) ndirs = 2; if (p > 0) ndirs = 3
73108

74-
if (idir == ndirs) then
75-
bub_adv_src = 0.d0; bub_r_src = 0.d0; bub_v_src = 0.d0
76-
bub_p_src = 0.d0; bub_m_src = 0.d0
109+
!$acc parallel loop collapse(3) gang vector default(present) private(Rtmp, Vtmp)
110+
do l = 0, p
111+
do k = 0, n
112+
do j = 0, m
113+
bub_adv_src(j, k, l) = 0d0
77114

78-
! advection source
79-
do j = 0, m; do k = 0, n; do l = 0, p
80-
! = 3 \alpha \bar{R^2 V} / \bar{R^3} = 4 pi nbub \bar{R^2 V}
81-
do q = 1, nb
82-
Rtmp(q) = q_prim_vf(bub_idx%rs(q))%sf(j, k, l)
83-
Vtmp(q) = q_prim_vf(bub_idx%vs(q))%sf(j, k, l)
84-
end do
115+
!$acc loop seq
116+
do q = 1, nb
117+
bub_r_src(j, k, l, q) = 0d0
118+
bub_v_src(j, k, l, q) = 0d0
119+
bub_p_src(j, k, l, q) = 0d0
120+
bub_m_src(j, k, l, q) = 0d0
121+
end do
122+
end do
123+
end do
124+
end do
85125

86-
! Computes n_bub number bubble density
87-
call s_comp_n_from_prim(q_prim_vf(alf_idx)%sf(j, k, l), &
88-
Rtmp, nbub(j, k, l))
89126

90-
call s_quad((Rtmp**2.d0)*Vtmp, R2Vav)
91-
bub_adv_src(j, k, l) = 4.d0*pi*nbub(j, k, l)*R2Vav
92-
end do; end do; end do
93127

94-
! bubble radius and radial velocity source
95-
do q = 1, nb; do j = 0, m; do k = 0, n; do l = 0, p
96-
bub_r_src(j, k, l, q) = q_cons_vf(bub_idx%vs(q))%sf(j, k, l)
128+
!$acc parallel loop collapse(3) gang vector default(present) private(Rtmp, Vtmp)
129+
do l = 0, p
130+
do k = 0, n
131+
do j = 0, m
97132

98-
call s_convert_to_mixture_variables(q_cons_vf, myRho, n_tait, B_tait, Re, j, k, l)
133+
!$acc loop seq
134+
do q = 1, nb
135+
Rtmp(q) = q_prim_vf(rs(q))%sf(j, k, l)
136+
Vtmp(q) = q_prim_vf(vs(q))%sf(j, k, l)
137+
end do
99138

100-
n_tait = 1.d0/n_tait + 1.d0 !make this the usual little 'gamma'
139+
call s_comp_n_from_prim(q_prim_vf(alf_idx)%sf(j, k, l), &
140+
Rtmp, nbub(j, k, l))
101141

102-
myRho = q_prim_vf(1)%sf(j, k, l)
103-
myP = q_prim_vf(E_idx)%sf(j, k, l)
104-
alf = q_prim_vf(alf_idx)%sf(j, k, l)
105-
myR = q_prim_vf(bub_idx%rs(q))%sf(j, k, l)
106-
myV = q_prim_vf(bub_idx%vs(q))%sf(j, k, l)
142+
call s_quad((Rtmp**2.d0)*Vtmp, R2Vav)
107143

108-
if (.not. polytropic) then
109-
pb = q_prim_vf(bub_idx%ps(q))%sf(j, k, l)
110-
mv = q_prim_vf(bub_idx%ms(q))%sf(j, k, l)
111-
call s_bwproperty(pb, q)
112-
vflux = f_vflux(myR, myV, mv, q)
113-
pbdot = f_bpres_dot(vflux, myR, myV, pb, mv, q)
114-
115-
bub_p_src(j, k, l, q) = nbub(j, k, l)*pbdot
116-
bub_m_src(j, k, l, q) = nbub(j, k, l)*vflux*4.d0*pi*(myR**2.d0)
117-
else
118-
pb = 0d0; mv = 0d0; vflux = 0d0; pbdot = 0d0
119-
end if
144+
bub_adv_src(j, k, l) = 4.d0*pi*nbub(j, k, l)*R2Vav
120145

121-
if (bubble_model == 1) then
122-
! Gilmore bubbles
123-
Cpinf = myP - pref
124-
Cpbw = f_cpbw(R0(q), myR, myV, pb)
125-
myH = f_H(Cpbw, Cpinf, n_tait, B_tait)
126-
c_gas = f_cgas(Cpinf, n_tait, B_tait, myH)
127-
Cpinf_dot = f_cpinfdot(myRho, myP, alf, n_tait, B_tait, bub_adv_src(j, k, l), mydivu%sf(j, k, l))
128-
myHdot = f_Hdot(Cpbw, Cpinf, Cpinf_dot, n_tait, B_tait, myR, myV, R0(q), pbdot)
129-
rddot = f_rddot(Cpbw, myR, myV, myH, myHdot, c_gas, n_tait, B_tait)
130-
else if (bubble_model == 2) then
131-
! Keller-Miksis bubbles
132-
Cpinf = myP
133-
Cpbw = f_cpbw_KM(R0(q), myR, myV, pb)
134-
! c_gas = dsqrt( n_tait*(Cpbw+B_tait) / myRho)
135-
c_liquid = DSQRT(n_tait*(myP + B_tait)/(myRho*(1.d0 - alf)))
136-
rddot = f_rddot_KM(pbdot, Cpinf, Cpbw, myRho, myR, myV, R0(q), c_liquid)
137-
else if (bubble_model == 3) then
138-
! Rayleigh-Plesset bubbles
139-
Cpbw = f_cpbw_KM(R0(q), myR, myV, pb)
140-
rddot = f_rddot_RP(myP, myRho, myR, myV, R0(q), Cpbw)
141-
end if
146+
end do
147+
end do
148+
end do
142149

143-
bub_v_src(j, k, l, q) = nbub(j, k, l)*rddot
150+
!$acc parallel loop collapse(3) gang vector default(present) private(myalpha_rho, myalpha)
151+
do l = 0, p
152+
do k = 0, n
153+
do j = 0, m
154+
!$acc loop seq
155+
do q = 1, nb
144156

145-
if (alf < 1.d-11) then
146-
bub_adv_src(j, k, l) = 0d0
147-
bub_r_src(j, k, l, q) = 0d0
148-
bub_v_src(j, k, l, q) = 0d0
149-
if (.not. polytropic) then
150-
bub_p_src(j, k, l, q) = 0d0
151-
bub_m_src(j, k, l, q) = 0d0
152-
end if
153-
end if
157+
bub_r_src(j, k, l, q) = q_cons_vf(vs(q))%sf(j, k, l)
154158

155-
end do; end do; end do; end do
156-
end if
159+
!$acc loop seq
160+
do ii = 1, num_fluids
161+
myalpha_rho(ii) = q_cons_vf(ii)%sf(j, k, l)
162+
myalpha(ii) = q_cons_vf(advxb + ii - 1)%sf(j, k, l)
163+
end do
164+
165+
myRho = 0d0
166+
n_tait = 0d0
167+
B_tait = 0d0
168+
169+
if (mpp_lim .and. (num_fluids > 2)) then
170+
!$acc loop seq
171+
do ii = 1, num_fluids
172+
myRho = myRho + myalpha_rho(ii)
173+
n_tait = n_tait + myalpha(ii)*gammas(ii)
174+
B_tait = B_tait + myalpha(ii)*pi_infs(ii)
175+
end do
176+
else if (num_fluids > 2) then
177+
!$acc loop seq
178+
do ii = 1, num_fluids - 1
179+
myRho = myRho + myalpha_rho(ii)
180+
n_tait = n_tait + myalpha(ii)*gammas(ii)
181+
B_tait = B_tait + myalpha(ii)*pi_infs(ii)
182+
end do
183+
else
184+
myRho = myalpha_rho(1)
185+
n_tait = gammas(1)
186+
B_tait = pi_infs(1)
187+
end if
188+
189+
n_tait = 1.d0/n_tait + 1.d0 !make this the usual little 'gamma'
190+
191+
myRho = q_prim_vf(1)%sf(j, k, l)
192+
myP = q_prim_vf(E_idx)%sf(j, k, l)
193+
alf = q_prim_vf(alf_idx)%sf(j, k, l)
194+
myR = q_prim_vf(rs(q))%sf(j, k, l)
195+
myV = q_prim_vf(vs(q))%sf(j, k, l)
196+
197+
if (.not. polytropic) then
198+
pb = q_prim_vf(ps(q))%sf(j, k, l)
199+
mv = q_prim_vf(ms(q))%sf(j, k, l)
200+
call s_bwproperty(pb, q)
201+
vflux = f_vflux(myR, myV, mv, q)
202+
pbdot = f_bpres_dot(vflux, myR, myV, pb, mv, q)
203+
204+
bub_p_src(j, k, l, q) = nbub(j, k, l)*pbdot
205+
bub_m_src(j, k, l, q) = nbub(j, k, l)*vflux*4.d0*pi*(myR**2.d0)
206+
else
207+
pb = 0d0; mv = 0d0; vflux = 0d0; pbdot = 0d0
208+
end if
209+
210+
if (bubble_model == 1) then
211+
! Gilmore bubbles
212+
Cpinf = myP - pref
213+
Cpbw = f_cpbw(R0(q), myR, myV, pb)
214+
myH = f_H(Cpbw, Cpinf, n_tait, B_tait)
215+
c_gas = f_cgas(Cpinf, n_tait, B_tait, myH)
216+
Cpinf_dot = f_cpinfdot(myRho, myP, alf, n_tait, B_tait, bub_adv_src(j, k, l), divu%sf(j, k, l))
217+
myHdot = f_Hdot(Cpbw, Cpinf, Cpinf_dot, n_tait, B_tait, myR, myV, R0(q), pbdot)
218+
rddot = f_rddot(Cpbw, myR, myV, myH, myHdot, c_gas, n_tait, B_tait)
219+
else if (bubble_model == 2) then
220+
! Keller-Miksis bubbles
221+
Cpinf = myP
222+
Cpbw = f_cpbw_KM(R0(q), myR, myV, pb)
223+
! c_gas = dsqrt( n_tait*(Cpbw+B_tait) / myRho)
224+
c_liquid = DSQRT(n_tait*(myP + B_tait)/(myRho*(1.d0 - alf)))
225+
rddot = f_rddot_KM(pbdot, Cpinf, Cpbw, myRho, myR, myV, R0(q), c_liquid)
226+
else if (bubble_model == 3) then
227+
! Rayleigh-Plesset bubbles
228+
Cpbw = f_cpbw_KM(R0(q), myR, myV, pb)
229+
rddot = f_rddot_RP(myP, myRho, myR, myV, R0(q), Cpbw)
230+
end if
231+
232+
bub_v_src(j, k, l, q) = nbub(j, k, l)*rddot
233+
234+
if (alf < 1.d-11) then
235+
bub_adv_src(j, k, l) = 0d0
236+
bub_r_src(j, k, l, q) = 0d0
237+
bub_v_src(j, k, l, q) = 0d0
238+
if (.not. polytropic) then
239+
bub_p_src(j, k, l, q) = 0d0
240+
bub_m_src(j, k, l, q) = 0d0
241+
end if
242+
end if
243+
end do
244+
end do
245+
end do
246+
end do
247+
248+
!$acc parallel loop collapse(3) gang vector default(present)
249+
do l = 0, p
250+
do q = 0, n
251+
do i = 0, m
252+
rhs_vf(alf_idx)%sf(i, q, l) = rhs_vf(alf_idx)%sf(i, q, l) + bub_adv_src(i, q, l)
253+
if (num_fluids > 1) rhs_vf(advxb)%sf(i, q, l) = &
254+
rhs_vf(advxb)%sf(i, q, l) - bub_adv_src(i, q, l)
255+
!$acc loop seq
256+
do k = 1, nb
257+
rhs_vf(rs(k))%sf(i, q, l) = rhs_vf(rs(k))%sf(i, q, l) + bub_r_src(i, q, l, k)
258+
rhs_vf(vs(k))%sf(i, q, l) = rhs_vf(vs(k))%sf(i, q, l) + bub_v_src(i, q, l, k)
259+
if (polytropic .neqv. .true.) then
260+
rhs_vf(ps(k))%sf(i, q, l) = rhs_vf(ps(k))%sf(i, q, l) + bub_p_src(i, q, l, k)
261+
rhs_vf(ms(k))%sf(i, q, l) = rhs_vf(ms(k))%sf(i, q, l) + bub_m_src(i, q, l, k)
262+
end if
263+
end do
264+
end do
265+
end do
266+
end do
157267

158268
end subroutine s_compute_bubble_source
159269

‎src/simulation/m_cbc.f90‎

Lines changed: 1 addition & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -132,14 +132,12 @@ end subroutine s_compute_abstract_L
132132

133133
integer :: dj
134134

135-
integer :: momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe
136135
integer :: bcxb, bcxe, bcyb, bcye, bczb, bcze
137-
real(kind(0d0)), allocatable, dimension(:) :: gammas, pi_infs
138136

139137
!$acc declare create(q_prim_rsx_vf, q_prim_rsy_vf, q_prim_rsz_vf, F_rsx_vf, F_src_rsx_vf,flux_rsx_vf, flux_src_rsx_vf, &
140138
!$acc F_rsy_vf, F_src_rsy_vf,flux_rsy_vf, flux_src_rsy_vf, F_rsz_vf, F_src_rsz_vf,flux_rsz_vf, flux_src_rsz_vf,alpha_rho,vel,adv,mf,Re, &
141139
!$acc dalpha_rho_ds,dvel_ds,dadv_ds,lambda,L,ds,fd_coef_x,fd_coef_y,fd_coef_z, &
142-
!$acc pi_coef_x,pi_coef_y,pi_coef_z, momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe, gammas, pi_infs, bcxb, bcxe, bcyb, bcye, bczb, bcze, is1, is2, is3, dj)
140+
!$acc pi_coef_x,pi_coef_y,pi_coef_z, bcxb, bcxe, bcyb, bcye, bczb, bcze, is1, is2, is3, dj)
143141

144142
contains
145143

@@ -312,14 +310,6 @@ subroutine s_initialize_cbc_module() ! ---------------------------------
312310
! Allocating the cell-width distribution in the s-direction
313311
allocate (ds(0:buff_size))
314312

315-
allocate (gammas(1:num_fluids), pi_infs(1:num_fluids))
316-
317-
do i = 1, num_fluids
318-
gammas(i) = fluid_pp(i)%gamma
319-
pi_infs(i) = fluid_pp(i)%pi_inf
320-
end do
321-
322-
!$acc update device(gammas, pi_infs)
323313

324314
! Allocating/Computing CBC Coefficients in x-direction =============
325315
if (all((/bc_x%beg, bc_x%end/) <= -5)) then
@@ -439,17 +429,6 @@ subroutine s_initialize_cbc_module() ! ---------------------------------
439429
! Associating the procedural pointer to the appropriate subroutine
440430
! that will be utilized in the conversion to the mixture variables
441431

442-
momxb = mom_idx%beg
443-
momxe = mom_idx%end
444-
advxb = adv_idx%beg
445-
advxe = adv_idx%end
446-
contxb = cont_idx%beg
447-
contxe = cont_idx%end
448-
bubxb = bub_idx%beg
449-
bubxe = bub_idx%end
450-
451-
!$acc update device(momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe)
452-
453432
bcxb = bc_x%beg
454433
bcxe = bc_x%end
455434

0 commit comments

Comments
 (0)