Skip to content

Commit 31420aa

Browse files
sbryngelsonsbryngelson
authored
Merge pull request #43 from wilfonba/refactoring
2 parents a39962c + e5bfdf5 commit 31420aa

16 files changed

+3219
-5836
lines changed

‎runTimeTests.sh‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
#!bash
2+
3+
for i in (0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
4+
./mfc.sh run ../MonopoleTests/1C0780C8Big/case.py -n 6 -g 6 -t pre_process simulation

‎src/simulation/m_bubbles.f90‎

Lines changed: 209 additions & 90 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, &
@@ -62,98 +97,182 @@ subroutine s_compute_bubble_source(idir, q_prim_vf, q_cons_vf, mydivu, &
6297
real(kind(0d0)) :: n_tait, B_tait
6398

6499
real(kind(0d0)), dimension(nb) :: Rtmp, Vtmp
65-
real(kind(0d0)) :: myR, myV, alf, myP, myRho, R2Vav
100+
real(kind(0d0)) :: myR, myV, alf, myP, myRho, R2Vav, R3
101+
real(kind(0d0)), dimension(num_fluids) :: myalpha, myalpha_rho
102+
real(kind(0d0)) :: start, finish
66103

67104
real(kind(0d0)), dimension(2) :: Re !< Reynolds number
68105

69-
integer :: j, k, l, q, s !< Loop variables
106+
integer :: i, j, k, l, q, ii !< Loop variables
70107
integer :: ndirs !< Number of coordinate directions
71108

72-
ndirs = 1; if (n > 0) ndirs = 2; if (p > 0) ndirs = 3
73-
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
77-
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)
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
114+
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
125+
126+
!$acc parallel loop collapse(3) gang vector default(present) private(Rtmp, Vtmp)
127+
do l = 0, p
128+
do k = 0, n
129+
do j = 0, m
130+
131+
!$acc loop seq
132+
do q = 1, nb
133+
Rtmp(q) = q_prim_vf(rs(q))%sf(j, k, l)
134+
Vtmp(q) = q_prim_vf(vs(q))%sf(j, k, l)
135+
end do
136+
137+
R3 = 0d0
138+
139+
!$acc loop seq
140+
do q = 1, nb
141+
R3 = R3 + weight(q)*Rtmp(q)**3.d0
142+
end do
143+
144+
nbub(j, k, l) = (3.d0/(4.d0*pi))*q_prim_vf(alf_idx)%sf(j, k, l)/R3
145+
146+
R2Vav = 0d0
147+
148+
!$acc loop seq
149+
do q = 1, nb
150+
R2Vav = R2Vav + weight(q)*Rtmp(q)**2.d0*Vtmp(q)
151+
end do
152+
153+
bub_adv_src(j, k, l) = 4.d0*pi*nbub(j, k, l)*R2Vav
154+
155+
end do
156+
end do
157+
end do
158+
159+
!$acc parallel loop collapse(3) gang vector default(present) private(myalpha_rho, myalpha)
160+
do l = 0, p
161+
do k = 0, n
162+
do j = 0, m
163+
!$acc loop seq
164+
do q = 1, nb
165+
166+
bub_r_src(j, k, l, q) = q_cons_vf(vs(q))%sf(j, k, l)
167+
168+
!$acc loop seq
169+
do ii = 1, num_fluids
170+
myalpha_rho(ii) = q_cons_vf(ii)%sf(j, k, l)
171+
myalpha(ii) = q_cons_vf(advxb + ii - 1)%sf(j, k, l)
84172
end do
85173

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

158277
end subroutine s_compute_bubble_source
159278

0 commit comments

Comments
 (0)