Skip to content

Commit ae827a5

Browse files
Hyeoksu LeeHyeoksu Lee
Hyeoksu Lee
authored and
Hyeoksu Lee
committed
update mixture rules
1 parent ff8ff17 commit ae827a5

File tree

3 files changed

+67
-159
lines changed

3 files changed

+67
-159
lines changed

‎src/common/m_variables_conversion.fpp‎

Lines changed: 58 additions & 118 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,7 @@ module m_variables_conversion
3939
s_convert_primitive_to_conservative_variables, &
4040
s_convert_primitive_to_flux_variables, &
4141
s_compute_pressure, &
42+
s_compute_species_fraction, &
4243
#ifndef MFC_PRE_PROCESS
4344
s_compute_speed_of_sound, &
4445
s_compute_fast_magnetosonic_speed, &
@@ -250,70 +251,31 @@ contains
250251
real(wp), intent(out), target :: qv
251252

252253
real(wp), optional, dimension(2), intent(out) :: Re_K
253-
!! Partial densities and volume fractions
254254
real(wp), optional, intent(out) :: G_K
255255
real(wp), optional, dimension(num_fluids), intent(in) :: G
256-
257256
real(wp), dimension(num_fluids) :: alpha_rho_K, alpha_K !<
258257

259258
integer :: i, j !< Generic loop iterator
260259

261260
! Computing the density, the specific heat ratio function and the
262261
! liquid stiffness function, respectively
263-
264-
if (igr) then
265-
if (num_fluids == 1) then
266-
alpha_rho_K(1) = q_vf(contxb)%sf(k, l, r)
267-
alpha_K(1) = 1._wp
268-
else
269-
do i = 1, num_fluids - 1
270-
alpha_rho_K(i) = q_vf(i)%sf(k, l, r)
271-
alpha_K(i) = q_vf(advxb + i - 1)%sf(k, l, r)
272-
end do
273-
274-
alpha_rho_K(num_fluids) = q_vf(num_fluids)%sf(k, l, r)
275-
alpha_K(num_fluids) = 1._wp - sum(alpha_K(1:num_fluids - 1))
276-
end if
277-
else
278-
do i = 1, num_fluids
279-
alpha_rho_K(i) = q_vf(i)%sf(k, l, r)
280-
alpha_K(i) = q_vf(advxb + i - 1)%sf(k, l, r)
281-
end do
282-
end if
283-
284-
if (mpp_lim) then
285-
do i = 1, num_fluids
286-
alpha_rho_K(i) = max(0._wp, alpha_rho_K(i))
287-
alpha_K(i) = min(max(0._wp, alpha_K(i)), 1._wp)
288-
end do
289-
290-
alpha_K = alpha_K/max(sum(alpha_K), 1.e-16_wp)
291-
end if
262+
call s_compute_species_fraction(q_vf, k, l, r, alpha_rho_K, alpha_K)
292263

293264
! Calculating the density, the specific heat ratio function, the
294265
! liquid stiffness function, and the energy reference function,
295266
! respectively, from the species analogs
296267
rho = 0._wp; gamma = 0._wp; pi_inf = 0._wp; qv = 0._wp
297-
298-
if (bubbles_euler) then
299-
rho = alpha_rho_K(1)
300-
gamma = gammas(1)
301-
pi_inf = pi_infs(1)
302-
qv = qvs(1)
303-
else
304-
do i = 1, num_fluids
305-
rho = rho + alpha_rho_K(i)
306-
gamma = gamma + alpha_K(i)*gammas(i)
307-
pi_inf = pi_inf + alpha_K(i)*pi_infs(i)
308-
qv = qv + alpha_rho_K(i)*qvs(i)
309-
end do
310-
end if
268+
do i = 1, num_fluids
269+
rho = rho + alpha_rho_K(i)
270+
gamma = gamma + alpha_K(i)*gammas(i)
271+
pi_inf = pi_inf + alpha_K(i)*pi_infs(i)
272+
qv = qv + alpha_rho_K(i)*qvs(i)
273+
end do
311274

312275
#ifdef MFC_SIMULATION
313276
! Computing the shear and bulk Reynolds numbers from species analogs
314277
if (viscous) then
315278
do i = 1, 2
316-
317279
Re_K(i) = dflt_real; if (Re_size(i) > 0) Re_K(i) = 0._wp
318280

319281
do j = 1, Re_size(i)
@@ -322,7 +284,6 @@ contains
322284
end do
323285

324286
Re_K(i) = 1._wp/max(Re_K(i), sgm_eps)
325-
326287
end do
327288
end if
328289
#endif
@@ -353,52 +314,25 @@ contains
353314
& parallelism='[seq]', cray_inline=True)
354315

355316
real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K
356-
357317
real(wp), dimension(num_fluids), intent(inout) :: alpha_rho_K, alpha_K !<
358318
real(wp), dimension(2), intent(out) :: Re_K
359-
!! Partial densities and volume fractions
360-
361319
real(wp), optional, intent(out) :: G_K
362320
real(wp), optional, dimension(num_fluids), intent(in) :: G
363321

364322
integer :: i, j !< Generic loop iterators
365-
real(wp) :: alpha_K_sum
366323

367324
#ifdef MFC_SIMULATION
368325
! Constraining the partial densities and the volume fractions within
369326
! their physical bounds to make sure that any mixture variables that
370327
! are derived from them result within the limits that are set by the
371328
! fluids physical parameters that make up the mixture
372-
rho_K = 0._wp
373-
gamma_K = 0._wp
374-
pi_inf_K = 0._wp
375-
qv_K = 0._wp
376-
377-
alpha_K_sum = 0._wp
378-
379-
if (mpp_lim) then
380-
do i = 1, num_fluids
381-
alpha_rho_K(i) = max(0._wp, alpha_rho_K(i))
382-
alpha_K(i) = min(max(0._wp, alpha_K(i)), 1._wp)
383-
alpha_K_sum = alpha_K_sum + alpha_K(i)
384-
end do
385-
386-
alpha_K = alpha_K/max(alpha_K_sum, sgm_eps)
387-
end if
388-
389-
if (bubbles_euler) then
390-
rho_K = alpha_rho_K(1)
391-
gamma_K = gammas(1)
392-
pi_inf_K = pi_infs(1)
393-
qv_K = qvs(1)
394-
else
395-
do i = 1, num_fluids
396-
rho_K = rho_K + alpha_rho_K(i)
397-
gamma_K = gamma_K + alpha_K(i)*gammas(i)
398-
pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i)
399-
qv_K = qv_K + alpha_rho_K(i)*qvs(i)
400-
end do
401-
end if
329+
rho_K = 0._wp; gamma_K = 0._wp; pi_inf_K = 0._wp; qv_K = 0._wp
330+
do i = 1, num_fluids
331+
rho_K = rho_K + alpha_rho_K(i)
332+
gamma_K = gamma_K + alpha_K(i)*gammas(i)
333+
pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i)
334+
qv_K = qv_K + alpha_rho_K(i)*qvs(i)
335+
end do
402336

403337
if (present(G_K)) then
404338
G_K = 0._wp
@@ -437,7 +371,6 @@ contains
437371

438372
$:GPU_ENTER_DATA(copyin='[is1b,is1e,is2b,is2e,is3b,is3e]')
439373

440-
#ifdef MFC_SIMULATION
441374
@:ALLOCATE(gammas (1:num_fluids))
442375
@:ALLOCATE(gs_min (1:num_fluids))
443376
@:ALLOCATE(pi_infs(1:num_fluids))
@@ -446,16 +379,6 @@ contains
446379
@:ALLOCATE(qvs (1:num_fluids))
447380
@:ALLOCATE(qvps (1:num_fluids))
448381
@:ALLOCATE(Gs_vc (1:num_fluids))
449-
#else
450-
@:ALLOCATE(gammas (1:num_fluids))
451-
@:ALLOCATE(gs_min (1:num_fluids))
452-
@:ALLOCATE(pi_infs(1:num_fluids))
453-
@:ALLOCATE(ps_inf(1:num_fluids))
454-
@:ALLOCATE(cvs (1:num_fluids))
455-
@:ALLOCATE(qvs (1:num_fluids))
456-
@:ALLOCATE(qvps (1:num_fluids))
457-
@:ALLOCATE(Gs_vc (1:num_fluids))
458-
#endif
459382

460383
do i = 1, num_fluids
461384
gammas(i) = fluid_pp(i)%gamma
@@ -484,12 +407,7 @@ contains
484407
#endif
485408

486409
if (bubbles_euler) then
487-
#ifdef MFC_SIMULATION
488410
@:ALLOCATE(bubrs_vc(1:nb))
489-
#else
490-
@:ALLOCATE(bubrs_vc(1:nb))
491-
#endif
492-
493411
do i = 1, nb
494412
bubrs_vc(i) = bub_idx%rs(i)
495413
end do
@@ -675,27 +593,7 @@ contains
675593
do j = ibounds(1)%beg, ibounds(1)%end
676594
dyn_pres_K = 0._wp
677595

678-
if (igr) then
679-
if (num_fluids == 1) then
680-
alpha_rho_K(1) = qK_cons_vf(contxb)%sf(j, k, l)
681-
alpha_K(1) = 1._wp
682-
else
683-
$:GPU_LOOP(parallelism='[seq]')
684-
do i = 1, num_fluids - 1
685-
alpha_rho_K(i) = qK_cons_vf(i)%sf(j, k, l)
686-
alpha_K(i) = qK_cons_vf(advxb + i - 1)%sf(j, k, l)
687-
end do
688-
689-
alpha_rho_K(num_fluids) = qK_cons_vf(num_fluids)%sf(j, k, l)
690-
alpha_K(num_fluids) = 1._wp - sum(alpha_K(1:num_fluids - 1))
691-
end if
692-
else
693-
$:GPU_LOOP(parallelism='[seq]')
694-
do i = 1, num_fluids
695-
alpha_rho_K(i) = qK_cons_vf(i)%sf(j, k, l)
696-
alpha_K(i) = qK_cons_vf(advxb + i - 1)%sf(j, k, l)
697-
end do
698-
end if
596+
call s_compute_species_fraction(qK_cons_vf, j, k, l, alpha_rho_K, alpha_K)
699597

700598
if (model_eqns /= 4) then
701599
#ifdef MFC_SIMULATION
@@ -1296,6 +1194,7 @@ contains
12961194
do i = advxb, advxe
12971195
alpha_K(i - E_idx) = qK_prim_vf(j, k, l, i)
12981196
end do
1197+
12991198
$:GPU_LOOP(parallelism='[seq]')
13001199
do i = 1, num_vels
13011200
vel_K(i) = qK_prim_vf(j, k, l, contxe + i)
@@ -1389,6 +1288,47 @@ contains
13891288
#endif
13901289
end subroutine s_convert_primitive_to_flux_variables
13911290

1291+
subroutine s_compute_species_fraction(q_vf, k, l, r, alpha_rho_K, alpha_K)
1292+
$:GPU_ROUTINE(function_name='s_compute_species_fraction', &
1293+
& parallelism='[seq]', cray_inline=True)
1294+
type(scalar_field), dimension(sys_size), intent(in) :: q_vf
1295+
integer, intent(in) :: k, l, r
1296+
real(wp), dimension(num_fluids), intent(out) :: alpha_rho_K, alpha_K
1297+
integer :: i
1298+
1299+
if (num_fluids == 1) then
1300+
alpha_rho_K(1) = q_vf(contxb)%sf(k, l, r)
1301+
if (igr .or. bubbles_euler) then
1302+
alpha_K(1) = 1._wp
1303+
else
1304+
alpha_K(1) = q_vf(advxb)%sf(k, l, r)
1305+
end if
1306+
else
1307+
if (igr) then
1308+
do i = 1, num_fluids - 1
1309+
alpha_rho_K(i) = q_vf(i)%sf(k, l, r)
1310+
alpha_K(i) = q_vf(advxb + i - 1)%sf(k, l, r)
1311+
end do
1312+
alpha_rho_K(num_fluids) = q_vf(num_fluids)%sf(k, l, r)
1313+
alpha_K(num_fluids) = 1._wp - sum(alpha_K(1:num_fluids - 1))
1314+
else
1315+
do i = 1, num_fluids
1316+
alpha_rho_K(i) = q_vf(i)%sf(k, l, r)
1317+
alpha_K(i) = q_vf(advxb + i - 1)%sf(k, l, r)
1318+
end do
1319+
end if
1320+
end if
1321+
1322+
if (mpp_lim) then
1323+
do i = 1, num_fluids
1324+
alpha_rho_K(i) = max(0._wp, alpha_rho_K(i))
1325+
alpha_K(i) = min(max(0._wp, alpha_K(i)), 1._wp)
1326+
end do
1327+
alpha_K = alpha_K/max(sum(alpha_K), 1.e-16_wp)
1328+
end if
1329+
1330+
end subroutine s_compute_species_fraction
1331+
13921332
impure subroutine s_finalize_variables_conversion_module()
13931333

13941334
! Deallocating the density, the specific heat ratio function and the

‎src/simulation/m_bubbles_EE.fpp‎

Lines changed: 6 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -246,30 +246,16 @@ contains
246246
n_tait = 0._wp
247247
B_tait = 0._wp
248248

249-
if (mpp_lim .and. (num_fluids > 2)) then
250-
$:GPU_LOOP(parallelism='[seq]')
251-
do ii = 1, num_fluids
252-
myRho = myRho + myalpha_rho(ii)
253-
n_tait = n_tait + myalpha(ii)*gammas(ii)
254-
B_tait = B_tait + myalpha(ii)*pi_infs(ii)
255-
end do
256-
else if (num_fluids > 2) then
257-
$:GPU_LOOP(parallelism='[seq]')
258-
do ii = 1, num_fluids - 1
259-
myRho = myRho + myalpha_rho(ii)
260-
n_tait = n_tait + myalpha(ii)*gammas(ii)
261-
B_tait = B_tait + myalpha(ii)*pi_infs(ii)
262-
end do
263-
else
264-
myRho = myalpha_rho(1)
265-
n_tait = gammas(1)
266-
B_tait = pi_infs(1)/pi_fac
267-
end if
249+
$:GPU_LOOP(parallelism='[seq]')
250+
do ii = 1, num_fluids
251+
myRho = myRho + myalpha_rho(ii)
252+
n_tait = n_tait + myalpha(ii)*gammas(ii)
253+
B_tait = B_tait + myalpha(ii)*pi_infs(ii)/pi_fac
254+
end do
268255

269256
n_tait = 1._wp/n_tait + 1._wp !make this the usual little 'gamma'
270257
B_tait = B_tait*(n_tait - 1)/n_tait ! make this the usual pi_inf
271258

272-
myRho = q_prim_vf(1)%sf(j, k, l)
273259
myP = q_prim_vf(E_idx)%sf(j, k, l)
274260
alf = q_prim_vf(alf_idx)%sf(j, k, l)
275261
myR = q_prim_vf(rs(q))%sf(j, k, l)

‎src/simulation/m_sim_helpers.fpp‎

Lines changed: 3 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -109,27 +109,7 @@ contains
109109

110110
integer :: i
111111

112-
if (igr) then
113-
if (num_fluids == 1) then
114-
alpha_rho(1) = q_prim_vf(contxb)%sf(j, k, l)
115-
alpha(1) = 1._wp
116-
else
117-
$:GPU_LOOP(parallelism='[seq]')
118-
do i = 1, num_fluids - 1
119-
alpha_rho(i) = q_prim_vf(i)%sf(j, k, l)
120-
alpha(i) = q_prim_vf(advxb + i - 1)%sf(j, k, l)
121-
end do
122-
123-
alpha_rho(num_fluids) = q_prim_vf(num_fluids)%sf(j, k, l)
124-
alpha(num_fluids) = 1._wp - sum(alpha(1:num_fluids - 1))
125-
end if
126-
else
127-
$:GPU_LOOP(parallelism='[seq]')
128-
do i = 1, num_fluids
129-
alpha_rho(i) = q_prim_vf(i)%sf(j, k, l)
130-
alpha(i) = q_prim_vf(advxb + i - 1)%sf(j, k, l)
131-
end do
132-
end if
112+
call s_compute_species_fraction(q_prim_vf, j, k, l, alpha_rho, alpha)
133113

134114
if (elasticity) then
135115
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, &
@@ -138,6 +118,8 @@ contains
138118
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re)
139119
end if
140120

121+
if (num_fluids == 1 .and. bubbles_euler) alpha = q_prim_vf(alf_idx)%sf(j, k, l)
122+
141123
if (igr) then
142124
$:GPU_LOOP(parallelism='[seq]')
143125
do i = 1, num_vels

0 commit comments

Comments
 (0)