Skip to content

Commit 405ea0a

Browse files
hyeoksu-leehyeoksu-lee
committed
Fix a bug on using pi_inf in subgrid bubble model and enable subgrid bubble models for viscous flow
1 parent dae6563 commit 405ea0a

File tree

14 files changed

+218
-42
lines changed

14 files changed

+218
-42
lines changed

‎examples/3D_turb_mixing/case.py‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -22,9 +22,9 @@
2222
Lz = 59.0
2323

2424
# Number of grid cells
25-
Nx = 255
26-
Ny = 255
27-
Nz = 255
25+
Nx = 191
26+
Ny = 191
27+
Nz = 191
2828

2929
# Grid spacing
3030
dx = Lx/float(Nx)
@@ -75,7 +75,6 @@
7575
'num_fluids' : 1,
7676
'adv_alphan' : 'T',
7777
'time_stepper' : 3,
78-
'weno_vars' : 2,
7978
'weno_order' : 5,
8079
'weno_eps' : 1.E-16,
8180
'weno_Re_flux' : 'T',

‎src/common/include/inline_conversions.fpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -32,11 +32,11 @@
3232

3333
if (mpp_lim .and. (num_fluids > 1)) then
3434
c = (1d0/gamma + 1d0)* &
35-
(pres + pi_inf)/rho
35+
(pres + pi_inf/(gamma+1d0))/rho
3636
else
3737
c = &
3838
(1d0/gamma + 1d0)* &
39-
(pres + pi_inf)/ &
39+
(pres + pi_inf/(gamma+1d0))/ &
4040
(rho*(1d0 - adv(num_fluids)))
4141
end if
4242
else
@@ -48,6 +48,7 @@
4848
else
4949
c = sqrt(c)
5050
end if
51+
5152
end subroutine s_compute_speed_of_sound
5253
#:enddef
5354

‎src/common/m_variables_conversion.fpp‎

Lines changed: 64 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -233,19 +233,33 @@ contains
233233
real(kind(0d0)), intent(OUT), target :: gamma
234234
real(kind(0d0)), intent(OUT), target :: pi_inf
235235

236+
real(kind(0d0)), dimension(num_fluids) :: alpha_rho_K, alpha_K
236237
real(kind(0d0)), optional, dimension(2), intent(OUT) :: Re_K
237238

238239
real(kind(0d0)), optional, intent(OUT) :: G_K
239240
real(kind(0d0)), optional, dimension(num_fluids), intent(IN) :: G
240241

241-
integer :: i
242+
integer :: i, q
242243

243244
! Constraining the partial densities and the volume fractions within
244245
! their physical bounds to make sure that any mixture variables that
245246
! are derived from them result within the limits that are set by the
246247
! fluids physical parameters that make up the mixture
247-
! alpha_rho_K(1) = qK_vf(i)%sf(i,j,k)
248-
! alpha_K(1) = qK_vf(E_idx+i)%sf(i,j,k)
248+
do i = 1, num_fluids
249+
alpha_rho_K(i) = q_vf(i)%sf(j, k, l)
250+
alpha_K(i) = q_vf(advxb + i - 1)%sf(j, k, l)
251+
end do
252+
253+
if (mpp_lim) then
254+
255+
do i = 1, num_fluids
256+
alpha_rho_K(i) = max(0d0, alpha_rho_K(i))
257+
alpha_K(i) = min(max(0d0, alpha_K(i)), 1d0)
258+
end do
259+
260+
alpha_K = alpha_K/max(sum(alpha_K), 1d-16)
261+
262+
end if
249263

250264
! Performing the transfer of the density, the specific heat ratio
251265
! function as well as the liquid stiffness function, respectively
@@ -284,6 +298,26 @@ contains
284298
end if
285299
end if
286300

301+
#ifdef MFC_SIMULATION
302+
! Computing the shear and bulk Reynolds numbers from species analogs
303+
if (any(Re_size > 0)) then
304+
if (num_fluids == 1) then ! need to consider case with num_fluids >= 2
305+
do i = 1, 2
306+
307+
Re_K(i) = dflt_real; if (Re_size(i) > 0) Re_K(i) = 0d0
308+
309+
do q = 1, Re_size(i)
310+
Re_K(i) = (1-alpha_K(Re_idx(i, q)))/fluid_pp(Re_idx(i, q))%Re(i) &
311+
+ Re_K(i)
312+
end do
313+
314+
Re_K(i) = 1d0/max(Re_K(i), sgm_eps)
315+
316+
end do
317+
end if
318+
end if
319+
#endif
320+
287321
! Post process requires rho_sf/gamma_sf/pi_inf_sf to also be updated
288322
#ifdef MFC_POST_PROCESS
289323
rho_sf (j, k, l) = rho
@@ -467,13 +501,16 @@ contains
467501

468502
subroutine s_convert_species_to_mixture_variables_bubbles_acc(rho_K, &
469503
gamma_K, pi_inf_K, &
470-
alpha_K, alpha_rho_K, k, l, r)
504+
alpha_K, alpha_rho_K, Re_K, k, l, r)
471505
!$acc routine seq
472506

473507
real(kind(0d0)), intent(INOUT) :: rho_K, gamma_K, pi_inf_K
474508

475509
real(kind(0d0)), dimension(num_fluids), intent(IN) :: alpha_rho_K, alpha_K !<
476510
!! Partial densities and volume fractions
511+
512+
real(kind(0d0)), dimension(2), intent(OUT) :: Re_K
513+
477514
integer, intent(IN) :: k, l, r
478515
integer :: i, j !< Generic loop iterators
479516

@@ -499,6 +536,25 @@ contains
499536
gamma_K = gammas(1)
500537
pi_inf_K = pi_infs(1)
501538
end if
539+
540+
if (any(Re_size > 0)) then
541+
if (num_fluids == 1) then ! need to consider case with num_fluids >= 2
542+
543+
do i = 1, 2
544+
Re_K(i) = dflt_real
545+
546+
if (Re_size(i) > 0) Re_K(i) = 0d0
547+
548+
do j = 1, Re_size(i)
549+
Re_K(i) = (1d0-alpha_K(Re_idx(i, j)))/Res(i, j) &
550+
+ Re_K(i)
551+
end do
552+
553+
Re_K(i) = 1d0/max(Re_K(i), sgm_eps)
554+
555+
end do
556+
end if
557+
end if
502558
#endif
503559

504560
end subroutine s_convert_species_to_mixture_variables_bubbles_acc
@@ -706,7 +762,7 @@ contains
706762
alpha_rho_K, Re_K, j, k, l, G_K, Gs)
707763
else if (bubbles) then
708764
call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, pi_inf_K, &
709-
alpha_K, alpha_rho_K, j, k, l)
765+
alpha_K, alpha_rho_K, Re_K, j, k, l)
710766
else
711767
call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, &
712768
alpha_K, alpha_rho_K, Re_K, j, k, l)
@@ -1004,9 +1060,9 @@ contains
10041060
call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, &
10051061
alpha_K, alpha_rho_K, Re_K, &
10061062
j, k, l, G_K, Gs)
1007-
! else if (bubbles) then
1008-
! call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, &
1009-
! pi_inf_K, alpha_K, alpha_rho_K, j, k, l)
1063+
else if (bubbles) then
1064+
call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, &
1065+
pi_inf_K, alpha_K, alpha_rho_K, Re_K, j, k, l)
10101066
else
10111067
call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, &
10121068
alpha_K, alpha_rho_K, Re_K, j, k, l)

‎src/pre_process/m_assign_variables.f90‎

Lines changed: 14 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -266,7 +266,8 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l,
266266
end if
267267

268268
! Computing Mixture Variables from Original Primitive Variables
269-
call s_convert_species_to_mixture_variables( &
269+
! call s_convert_species_to_mixture_variables( &
270+
call s_convert_to_mixture_variables( &
270271
q_prim_vf, j, k, l, &
271272
orig_rho, &
272273
orig_gamma, &
@@ -333,17 +334,22 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l,
333334
end if
334335
else
335336
q_prim_vf(bub_idx%rs(i))%sf(j, k, l) = muR
336-
q_prim_vf(bub_idx%vs(i))%sf(j, k, l) = muV
337+
if (.not. vel_profile) then
338+
q_prim_vf(bub_idx%vs(i))%sf(j, k, l) = muV
339+
else
340+
q_prim_vf(bub_idx%vs(i))%sf(j, k, l) = muV*tanh(y_cc(k))
341+
end if
337342
if (.not. polytropic) then
338343
q_prim_vf(bub_idx%ps(i))%sf(j, k, l) = patch_icpp(patch_id)%p0
339344
q_prim_vf(bub_idx%ms(i))%sf(j, k, l) = patch_icpp(patch_id)%m0
340345
end if
341346
end if
342347
end do
343348
end if
344-
349+
345350
! Density and the specific heat ratio and liquid stiffness functions
346-
call s_convert_species_to_mixture_variables( &
351+
! call s_convert_species_to_mixture_variables( &
352+
call s_convert_to_mixture_variables( &
347353
q_prim_vf, j, k, l, &
348354
patch_icpp(patch_id)%rho, &
349355
patch_icpp(patch_id)%gamma, &
@@ -412,7 +418,8 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l,
412418
end if
413419

414420
! Density and the specific heat ratio and liquid stiffness functions
415-
call s_convert_species_to_mixture_variables( &
421+
! call s_convert_species_to_mixture_variables( &
422+
call s_convert_to_mixture_variables( &
416423
q_prim_vf, j, k, l, &
417424
patch_icpp(smooth_patch_id)%rho, &
418425
patch_icpp(smooth_patch_id)%gamma, &
@@ -475,7 +482,8 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l,
475482
end if
476483

477484
! Density and the specific heat ratio and liquid stiffness functions
478-
call s_convert_species_to_mixture_variables(q_prim_vf, j, k, l, &
485+
! call s_convert_species_to_mixture_variables(q_prim_vf, j, k, l, &
486+
call s_convert_to_mixture_variables(q_prim_vf, j, k, l, &
479487
rho, gamma, pi_inf)
480488

481489
! Velocity

‎src/pre_process/m_initial_condition.fpp‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -377,10 +377,14 @@ contains
377377
! Set fluid flow properties
378378
gam = 1.+1./fluid_pp(1)%gamma
379379
pi_inf = fluid_pp(1)%pi_inf*(gam-1.)/gam
380-
rho1 = patch_icpp(1)%alpha_rho(1)/patch_icpp(1)%alpha(1)
380+
if (bubbles .and. num_fluids == 1) then
381+
rho1 = patch_icpp(1)%alpha_rho(1)/(1d0-patch_icpp(1)%alpha(1))
382+
else
383+
rho1 = patch_icpp(1)%alpha_rho(1)/patch_icpp(1)%alpha(1)
384+
end if
381385
c1 = sqrt((gam*(patch_icpp(1)%pres+pi_inf))/rho1)
382386
mach = 1./c1
383-
387+
384388
! Assign mean profiles
385389
do j=0,n
386390
u_mean(j)=tanh(y_cc(j))

‎src/pre_process/m_patches.f90‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -505,7 +505,7 @@ subroutine s_rectangle(patch_id, patch_id_fp, q_prim_vf) ! ---------------------
505505
real(kind(0d0)) :: pi_inf, gamma, lit_gamma !< Equation of state parameters
506506

507507
integer :: i, j !< generic loop iterators
508-
508+
509509
pi_inf = fluid_pp(1)%pi_inf
510510
gamma = fluid_pp(1)%gamma
511511
lit_gamma = (1d0 + gamma)/gamma

‎src/simulation/m_bubbles.fpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -195,8 +195,9 @@ contains
195195
n_tait = gammas(1)
196196
B_tait = pi_infs(1)
197197
end if
198-
198+
199199
n_tait = 1.d0/n_tait + 1.d0 !make this the usual little 'gamma'
200+
B_tait = B_tait*(n_tait-1)/n_tait ! make this the usual pi_inf
200201

201202
myRho = q_prim_vf(1)%sf(j, k, l)
202203
myP = q_prim_vf(E_idx)%sf(j, k, l)

‎src/simulation/m_cbc.fpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -750,7 +750,7 @@ contains
750750
end do
751751

752752
if (bubbles) then
753-
call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, adv, alpha_rho, 0, k, r)
753+
call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, adv, alpha_rho, Re_cbc, 0, k, r)
754754

755755
else
756756
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, adv, alpha_rho, Re_cbc, 0, k, r)
@@ -767,7 +767,7 @@ contains
767767

768768
! Compute mixture sound speed
769769
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_K_sum, c)
770-
770+
! write(44,*) cbc_dir,r,k,adv,c
771771
! ============================================================
772772

773773
! First-Order Spatial Derivatives of Primitive Variables =====

‎src/simulation/m_data_output.fpp‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -255,7 +255,7 @@ contains
255255
end do
256256

257257
if (bubbles) then
258-
call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, alpha, alpha_rho, j, k, l)
258+
call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, alpha, alpha_rho, Re, j, k, l)
259259
else
260260
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, alpha, alpha_rho, Re, j, k, l)
261261
end if
@@ -270,7 +270,7 @@ contains
270270
end do
271271

272272
pres = q_prim_vf(E_idx)%sf(j, k, l)
273-
273+
274274
E = gamma*pres + pi_inf + 5d-1*rho*vel_sum
275275

276276
H = (E + pres)/rho
@@ -302,7 +302,6 @@ contains
302302
end if
303303

304304
if (any(Re_size > 0)) then
305-
306305
if (grid_geometry == 3) then
307306
vcfl_sf(j, k, l) = maxval(dt/Re) &
308307
/min(dx(j), dy(k), fltr_dtheta)**2d0

‎src/simulation/m_rhs.fpp‎

Lines changed: 48 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -769,6 +769,15 @@ contains
769769
qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, &
770770
id)
771771

772+
if (bubbles) then
773+
iv%beg = bubxb; iv%end = bubxe
774+
call s_reconstruct_cell_boundary_values( &
775+
q_prim_qp%vf(iv%beg:iv%end), &
776+
qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, &
777+
qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, &
778+
id)
779+
end if
780+
772781
iv%beg = mom_idx%beg; iv%end = mom_idx%end
773782
if (weno_Re_flux) then
774783
call s_reconstruct_cell_boundary_values_visc_deriv( &
@@ -810,6 +819,28 @@ contains
810819
call nvtxStartRange("RHS-Riemann")
811820

812821
! Computing Riemann Solver Flux and Source Flux =================
822+
! do l = 0, p
823+
! do k = 0, n
824+
! do j = 0, m
825+
! do q = 1, sys_size
826+
! if(isnan(q_prim_qp%vf(q)%sf(j,k,l))) then
827+
! write(*,*) t_step,id,q,j,k,l,'q_prim_qp'
828+
! end if
829+
! end do
830+
! do q = momxb,momxe
831+
! if(isnan(dqR_prim_dx_n(id)%vf(q)%sf(j,k,l))) then
832+
! write(*,*) t_step,id,q,j,k,l,'dqR_prim_dx_n'
833+
! end if
834+
! if(isnan(dqR_prim_dy_n(id)%vf(q)%sf(j,k,l))) then
835+
! write(*,*) t_step,id,q,j,k,l,'dqR_prim_dy_n'
836+
! end if
837+
! if(isnan(dqR_prim_dz_n(id)%vf(q)%sf(j,k,l))) then
838+
! write(*,*) t_step,id,q,j,k,l,'dqR_prim_dz_n'
839+
! end if
840+
! end do
841+
! end do
842+
! end do
843+
! end do
813844

814845
call s_riemann_solver(qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, &
815846
dqR_prim_dx_n(id)%vf, &
@@ -827,7 +858,22 @@ contains
827858
flux_gsrc_n(id)%vf, &
828859
id, ix, iy, iz)
829860
call nvtxEndRange
830-
861+
! do l = 0, p
862+
! do k = 0, n
863+
! do j = 0, m
864+
! do q = 1, sys_size
865+
! if(isnan(flux_n(id)%vf(q)%sf(j,k,l))) then
866+
! write(*,*) t_step,id,q,j,k,l,'flux_n'
867+
! end if
868+
! end do
869+
! do q = momxb,advxe
870+
! if(isnan(flux_src_n(id)%vf(q)%sf(j,k,l))) then
871+
! write(*,*) t_step,id,q,j,k,l,'flux_src_n'
872+
! end if
873+
! end do
874+
! end do
875+
! end do
876+
! end do
831877
! ===============================================================
832878

833879

@@ -2586,4 +2632,4 @@ contains
25862632

25872633
end subroutine s_finalize_rhs_module ! ---------------------------------
25882634

2589-
end module m_rhs
2635+
end module m_rhs

0 commit comments

Comments
 (0)