Skip to content

Commit cc3572b

Browse files
sbryngelsonsbryngelson
authored
Merge pull request #86 from anshgupta1234/conversions_update
2 parents 44c89b2 + 3ef9f69 commit cc3572b

File tree

12 files changed

+280
-254
lines changed

12 files changed

+280
-254
lines changed

‎examples/2D_shockbubble/case.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -63,7 +63,7 @@
6363
'bc_x%end' : -6,
6464
'bc_y%beg' : -6,
6565
'bc_y%end' : -6,
66-
# =================================================== =======================
66+
# ==========================================================================
6767

6868
# Formatted Database Files Structure Parameters ============================
6969
'format' : 1,

‎src/common/inline_conversions.fpp‎

Lines changed: 50 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -1,56 +1,53 @@
1-
!> This procedure calculates the speed of sound
2-
!! @param pres Pressure
3-
!! @param rho Cell averaged density
4-
!! @param pi_inf Cell averaged liquid stiffness
5-
!! @param gamma Cell averaged specific heat ratio
6-
!! @param H Cell averaged enthalpy
7-
!! @param adv Advection Variables
8-
!! @param vel_sum Sum of all velocities
9-
!! @param q_prim_vf Primitive vars in 1 direction
10-
!! @param flg Helps determine which conditionals to be called.
11-
! flg >= 2: Check all conditionals
12-
! flg = 1: Check for alt_soundspeed, otherwise run the 3rd conditional block
13-
! flg = 0: Check for alt_soundspeed, otherwise use enthalpy
14-
!! @param c Speed of sound
15-
16-
#:def compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, q_prim_vf, j, k, l, flg, c)
17-
18-
if (alt_soundspeed .and. ${flg}$ >= 0) then
19-
blkmod1 = ((gammas(1) + 1d0)*${pres}$ + &
20-
pi_infs(1))/gammas(1)
21-
blkmod2 = ((gammas(2) + 1d0)*${pres}$ + &
22-
pi_infs(2))/gammas(2)
23-
${c}$ = (1d0/(${rho}$*(${adv}$(1)/blkmod1 + ${adv}$(2)/blkmod2)))
24-
elseif (model_eqns == 3 .and. ${flg}$ >= 2) then
25-
${c}$ = 0d0
1+
#:def s_compute_speed_of_sound()
2+
subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c)
3+
4+
real(kind(0d0)), intent(IN) :: pres
5+
real(kind(0d0)), intent(IN) :: rho, gamma, pi_inf
6+
real(kind(0d0)), intent(IN) :: H
7+
real(kind(0d0)), dimension(num_fluids), intent(IN) :: adv
8+
real(kind(0d0)), intent(IN) :: vel_sum
9+
real(kind(0d0)), intent(OUT) :: c
10+
11+
real(kind(0d0)) :: blkmod1, blkmod2
12+
13+
integer :: q
14+
15+
if (alt_soundspeed) then
16+
blkmod1 = ((gammas(1) + 1d0)*pres + &
17+
pi_infs(1))/gammas(1)
18+
blkmod2 = ((gammas(2) + 1d0)*pres + &
19+
pi_infs(2))/gammas(2)
20+
c = (1d0/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2)))
21+
elseif (model_eqns == 3) then
22+
c = 0d0
2623
!$acc loop seq
27-
do q = 1, num_fluids
28-
${c}$ = ${c}$ + ${q_prim_vf}$(${j}$, ${k}$, ${l}$, q + advxb - 1)*(1d0/gammas(q) + 1d0)* &
29-
(${q_prim_vf}$(${j}$, ${k}$, ${l}$, E_idx) + pi_infs(q)/(gammas(q) + 1d0))
30-
end do
31-
${c}$ = ${c}$/${rho}$
32-
33-
elseif (((model_eqns == 4) .or. (model_eqns == 2 .and. bubbles) .or. ${flg}$ == 1) &
34-
.and. ${flg}$ >= 1) then
35-
! Sound speed for bubble mmixture to order O(\alpha)
36-
37-
if (mpp_lim .and. (num_fluids > 1)) then
38-
${c}$ = (1d0/${gamma}$ + 1d0)* &
39-
(${pres}$ + ${pi_inf}$)/${rho}$
24+
do q = 1, num_fluids
25+
c = c + adv(q)*(1d0/gammas(q) + 1d0)* &
26+
(pres + pi_infs(q)/(gammas(q) + 1d0))
27+
end do
28+
c = c/rho
29+
30+
elseif (((model_eqns == 4) .or. (model_eqns == 2 .and. bubbles))) then
31+
! Sound speed for bubble mmixture to order O(\alpha)
32+
33+
if (mpp_lim .and. (num_fluids > 1)) then
34+
c = (1d0/gamma + 1d0)* &
35+
(pres + pi_inf)/rho
36+
else
37+
c = &
38+
(1d0/gamma + 1d0)* &
39+
(pres + pi_inf)/ &
40+
(rho*(1d0 - adv(num_fluids)))
41+
end if
42+
else
43+
c = ((H - 5d-1*vel_sum)/gamma)
44+
end if
45+
46+
if (mixture_err .and. c < 0d0) then
47+
c = 100.d0*sgm_eps
4048
else
41-
${c}$ = &
42-
(1d0/${gamma}$ + 1d0)* &
43-
(${pres}$ + ${pi_inf}$)/ &
44-
(${rho}$*(1d0 - ${adv}$(num_fluids)))
49+
c = sqrt(c)
4550
end if
46-
else
47-
${c}$ = ((${H}$ - 5d-1*${vel_sum}$)/${gamma}$)
48-
end if
49-
50-
if (mixture_err .and. ${c}$ < 0d0) then
51-
${c}$ = 100.d0*sgm_eps
52-
else
53-
${c}$ = sqrt(${c}$)
54-
end if
55-
56-
#:enddef compute_speed_of_sound
51+
end subroutine s_compute_speed_of_sound
52+
#:enddef
53+

‎src/common/m_derived_types.mod‎

-4.56 KB
Binary file not shown.

‎src/common/m_variables_conversion.fpp‎

Lines changed: 34 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,7 @@ module m_variables_conversion
7777

7878
!! In simulation, gammas and pi_infs is already declared in m_global_variables
7979
#ifndef MFC_SIMULATION
80-
real(kind(0d0)), allocatable, dimension(:) :: gammas, pi_infs
80+
real(kind(0d0)), allocatable, public, dimension(:) :: gammas, pi_infs
8181
!$acc declare create(gammas, pi_infs)
8282
#endif
8383

@@ -103,19 +103,26 @@ contains
103103
!> This procedure conditionally calculates the appropriate pressure
104104
!! @param energy Energy
105105
!! @param alf Void Fraction
106+
!! @param stress Shear Stress
107+
!! @param mom Momentum
106108
!! @param dyn_p Dynamic Pressure
107109
!! @param pi_inf Liquid Stiffness
108110
!! @param gamma Specific Heat Ratio
109111
!! @param pres Pressure to calculate
110-
subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, pres)
112+
subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, pres, stress, mom, G)
111113
!$acc routine seq
112114

113-
real(kind(0d0)), intent(IN) :: energy, alf
115+
real(kind(0d0)), intent(IN) :: energy, alf
116+
real(kind(0d0)), intent(IN), optional :: stress, mom, G
114117

115118
real(kind(0d0)), intent(IN) :: dyn_p
116119
real(kind(0d0)), intent(OUT) :: pres
117120

118-
real(kind(0d0)), intent(IN) :: pi_inf, gamma
121+
real(kind(0d0)), intent(IN) :: pi_inf, gamma, rho
122+
123+
real(kind(0d0)) :: E_e
124+
125+
integer :: s !< Generic loop iterator
119126

120127
! Depending on model_eqns and bubbles, the appropriate procedure
121128
! for computing pressure is targeted by the procedure pointer
@@ -131,6 +138,28 @@ contains
131138
)**(1/gamma + 1) - pi_inf
132139
end if
133140

141+
if (hypoelasticity .and. present(G)) then
142+
! calculate elastic contribution to Energy
143+
E_e = 0d0
144+
do s = stress_idx%beg, stress_idx%end
145+
if (G > 0) then
146+
E_e = E_e + ((stress/rho)**2d0)/(4d0*G)
147+
! Additional terms in 2D and 3D
148+
if ((s == stress_idx%beg + 1) .or. &
149+
(s == stress_idx%beg + 3) .or. &
150+
(s == stress_idx%beg + 4)) then
151+
E_e = E_e + ((stress/rho)**2d0)/(4d0*G)
152+
end if
153+
end if
154+
end do
155+
156+
pres = ( &
157+
energy - &
158+
0.5d0*(mom**2.d0)/rho - &
159+
pi_inf - E_e &
160+
)/gamma
161+
end if
162+
134163
end subroutine s_compute_pressure
135164

136165
!> This subroutine is designed for the gamma/pi_inf model
@@ -707,7 +736,7 @@ contains
707736
end do
708737
call s_compute_pressure(qK_cons_vf(E_idx)%sf(j, k, l), &
709738
qK_cons_vf(alf_idx)%sf(j, k, l), &
710-
dyn_pres_K, pi_inf_K, gamma_K, pres)
739+
dyn_pres_K, pi_inf_K, gamma_K, rho_K, pres)
711740

712741
qK_prim_vf(E_idx)%sf(j, k, l) = pres
713742

‎src/post_process/m_derived_variables.f90‎ renamed to ‎src/post_process/m_derived_variables.fpp‎

Lines changed: 71 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,9 @@
77
!! Currently, the available derived variables include the unadvected
88
!! volume fraction, specific heat ratio, liquid stiffness, speed of
99
!! sound, vorticity and the numerical Schlieren function.
10+
11+
#:include 'inline_conversions.fpp'
12+
1013
module m_derived_variables
1114

1215
! Dependencies =============================================================
@@ -15,6 +18,8 @@ module m_derived_variables
1518
use m_global_parameters !< Global parameters for the code
1619

1720
use m_mpi_proxy !< Message passing interface (MPI) module proxy
21+
22+
use m_variables_conversion
1823
! ==========================================================================
1924

2025
implicit none
@@ -28,7 +33,8 @@ module m_derived_variables
2833
s_derive_vorticity_component, &
2934
s_derive_qm, &
3035
s_derive_numerical_schlieren_function, &
31-
s_finalize_derived_variables_module
36+
s_finalize_derived_variables_module, &
37+
s_compute_speed_of_sound
3238

3339
real(kind(0d0)), allocatable, dimension(:, :, :) :: gm_rho_sf !<
3440
!! Gradient magnitude (gm) of the density for each cell of the computational
@@ -570,14 +576,14 @@ subroutine s_derive_vorticity_component(i, q_prim_vf, q_sf) ! ----------
570576
571577
end subroutine s_derive_vorticity_component ! --------------------------
572578
573-
!> This subroutine gets as inputs the primitive variables. From those
574-
!! inputs, it proceeds to calculate the value of the Q_M
575-
!! function, which are subsequently stored in the derived flow
576-
!! quantity storage variable, q_sf.
577-
!! @param q_prim_vf Primitive variables
578-
!! @param q_sf Q_M
579-
subroutine s_derive_qm(q_prim_vf,q_sf)
580-
type(scalar_field), &
579+
!> This subroutine gets as inputs the primitive variables. From those
580+
!! inputs, it proceeds to calculate the value of the Q_M
581+
!! function, which are subsequently stored in the derived flow
582+
!! quantity storage variable, q_sf.
583+
!! @param q_prim_vf Primitive variables
584+
!! @param q_sf Q_M
585+
subroutine s_derive_qm(q_prim_vf,q_sf)
586+
type(scalar_field), &
581587
dimension(sys_size), &
582588
intent(IN) :: q_prim_vf
583589
@@ -590,71 +596,73 @@ subroutine s_derive_qm(q_prim_vf,q_sf)
590596
real(kind(0d0)), &
591597
dimension(1:3, 1:3) :: q_jacobian_sf, S, S2, O, O2
592598
593-
real(kind(0d0)) :: trS, trS2, trO2, Q, IIS
599+
real(kind(0d0)) :: trS, trS2, trO2, Q, IIS
594600
integer :: j, k, l, r, jj, kk !< Generic loop iterators
595601
596602
do l = -offset_z%beg, p + offset_z%end
597603
do k = -offset_y%beg, n + offset_y%end
598604
do j = -offset_x%beg, m + offset_x%end
599605
600-
! Get velocity gradient tensor
606+
! Get velocity gradient tensor
601607
q_jacobian_sf(:, :) = 0d0
602-
608+
603609
do r = -fd_number, fd_number
604-
do jj = 1, 3
605-
! d()/dx
606-
q_jacobian_sf(jj, 1) = &
607-
q_jacobian_sf(jj, 1)+ &
608-
fd_coeff_x(r, j)* &
609-
q_prim_vf(mom_idx%beg+jj-1)%sf(r + j, k, l)
610-
! d()/dy
611-
q_jacobian_sf(jj, 2) = &
612-
q_jacobian_sf(jj, 2)+ &
613-
fd_coeff_y(r, k)* &
614-
q_prim_vf(mom_idx%beg+jj-1)%sf(j, r + k, l)
615-
! d()/dz
616-
q_jacobian_sf(jj, 3) = &
617-
q_jacobian_sf(jj, 3)+ &
618-
fd_coeff_z(r, l)* &
619-
q_prim_vf(mom_idx%beg+jj-1)%sf(j, k, r + l)
620-
end do
621-
end do
622-
623-
! Decompose J into asymmetric matrix, S, and a skew-symmetric matrix, O
624-
do jj = 1, 3
625-
do kk = 1, 3
626-
S(jj, kk) = 0.5D0* &
627-
(q_jacobian_sf(jj, kk) + q_jacobian_sf(kk, jj))
628-
O(jj, kk) = 0.5D0* &
629-
(q_jacobian_sf(jj, kk) - q_jacobian_sf(kk, jj))
630-
end do
631-
end do
632-
633-
! Compute S2 = S*S'
634-
do jj = 1, 3
635-
do kk = 1, 3
636-
O2(jj, kk) = O(jj,1)*O(kk,1)+ &
637-
O(jj,2)*O(kk,2)+ &
638-
O(jj,3)*O(kk,3)
639-
S2(jj, kk) = S(jj,1)*S(kk,1)+ &
640-
S(jj,2)*S(kk,2)+ &
641-
S(jj,3)*S(kk,3)
642-
end do
643-
end do
644-
645-
! Compute Q
646-
Q = 0.5*((O2(1,1)+O2(2,2)+O2(3,3))- &
647-
(S2(1,1)+S2(2,2)+S2(3,3)))
648-
trS = S(1,1)+S(2,2)+S(3,3)
649-
IIS = 0.5*((S(1,1)+S(2,2)+S(3,3))**2- &
650-
(S2(1,1)+S2(2,2)+S2(3,3)))
651-
q_sf(j, k, l) = Q+IIS
610+
do jj = 1, 3
611+
! d()/dx
612+
q_jacobian_sf(jj, 1) = &
613+
q_jacobian_sf(jj, 1)+ &
614+
fd_coeff_x(r, j)* &
615+
q_prim_vf(mom_idx%beg+jj-1)%sf(r + j, k, l)
616+
! d()/dy
617+
q_jacobian_sf(jj, 2) = &
618+
q_jacobian_sf(jj, 2)+ &
619+
fd_coeff_y(r, k)* &
620+
q_prim_vf(mom_idx%beg+jj-1)%sf(j, r + k, l)
621+
! d()/dz
622+
q_jacobian_sf(jj, 3) = &
623+
q_jacobian_sf(jj, 3)+ &
624+
fd_coeff_z(r, l)* &
625+
q_prim_vf(mom_idx%beg+jj-1)%sf(j, k, r + l)
626+
end do
627+
end do
628+
629+
! Decompose J into asymmetric matrix, S, and a skew-symmetric matrix, O
630+
do jj = 1, 3
631+
do kk = 1, 3
632+
S(jj, kk) = 0.5D0* &
633+
(q_jacobian_sf(jj, kk) + q_jacobian_sf(kk, jj))
634+
O(jj, kk) = 0.5D0* &
635+
(q_jacobian_sf(jj, kk) - q_jacobian_sf(kk, jj))
636+
end do
637+
end do
638+
639+
! Compute S2 = S*S'
640+
do jj = 1, 3
641+
do kk = 1, 3
642+
O2(jj, kk) = O(jj,1)*O(kk,1)+ &
643+
O(jj,2)*O(kk,2)+ &
644+
O(jj,3)*O(kk,3)
645+
S2(jj, kk) = S(jj,1)*S(kk,1)+ &
646+
S(jj,2)*S(kk,2)+ &
647+
S(jj,3)*S(kk,3)
648+
end do
649+
end do
650+
651+
! Compute Q
652+
Q = 0.5*((O2(1,1)+O2(2,2)+O2(3,3))- &
653+
(S2(1,1)+S2(2,2)+S2(3,3)))
654+
trS = S(1,1)+S(2,2)+S(3,3)
655+
IIS = 0.5*((S(1,1)+S(2,2)+S(3,3))**2- &
656+
(S2(1,1)+S2(2,2)+S2(3,3)))
657+
q_sf(j, k, l) = Q+IIS
652658

653659
end do
654660
end do
655-
end do
661+
end do
662+
663+
end subroutine s_derive_qm
656664

657-
end subroutine s_derive_qm
665+
@:s_compute_speed_of_sound()
658666

659667
!> This subroutine gets as inputs the conservative variables
660668
!! and density. From those inputs, it proceeds to calculate
@@ -833,4 +841,4 @@ subroutine s_finalize_derived_variables_module() ! -------------------
833841
834842
end subroutine s_finalize_derived_variables_module ! -----------------
835843
836-
end module m_derived_variables
844+
end module m_derived_variables

0 commit comments

Comments
 (0)