Skip to content

Commit cf81cd8

Browse files
sbryngelsonsbryngelson
authored
Merge pull request #62 from sbryngelson/master
2 parents 45be9fa + 75804a0 commit cf81cd8

File tree

2 files changed

+86
-206
lines changed

2 files changed

+86
-206
lines changed

‎src/simulation/inline_riemann.fpp‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
#:def arithmetic_avg()
2+
rho_avg = 5d-1*(rho_L + rho_R)
3+
vel_avg_rms = 0d0
4+
!$acc loop seq
5+
do i = 1, num_dims
6+
vel_avg_rms = vel_avg_rms + (5d-1*(vel_L(i) + vel_R(i)))**2d0
7+
end do
8+
9+
H_avg = 5d-1*(H_L + H_R)
10+
gamma_avg = 5d-1*(gamma_L + gamma_R)
11+
12+
13+
14+
#:enddef arithmetic_avg
15+
16+
17+
#:def roe_avg()
18+
rho_avg = sqrt(rho_L*rho_R)
19+
vel_avg_rms = 0d0
20+
!$acc loop seq
21+
do i = 1, num_dims
22+
vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(i) + sqrt(rho_R)*vel_R(i))**2d0/ &
23+
(sqrt(rho_L) + sqrt(rho_R))**2d0
24+
end do
25+
26+
H_avg = (sqrt(rho_L)*H_L + sqrt(rho_R)*H_R)/ &
27+
(sqrt(rho_L) + sqrt(rho_R))
28+
29+
gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/ &
30+
(sqrt(rho_L) + sqrt(rho_R))
31+
32+
rho_avg = sqrt(rho_L*rho_R)
33+
vel_avg_rms = (sqrt(rho_L)*vel_L(1) + sqrt(rho_R)*vel_R(1))**2d0/ &
34+
(sqrt(rho_L) + sqrt(rho_R))**2d0
35+
36+
#:enddef roe_avg
37+
38+
#:def compute_average_state()
39+
40+
if (avg_state == 1) then
41+
@:roe_avg()
42+
end if
43+
44+
if (avg_state == 2) then
45+
@:arithmetic_avg()
46+
end if
47+
48+
#:enddef compute_average_state

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 38 additions & 206 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,8 @@
1717
!! 1) Harten-Lax-van Leer (HLL)
1818
!! 2) Harten-Lax-van Leer-Contact (HLLC)
1919
!! 3) Exact
20+
#:include 'inline_riemann.fpp'
21+
2022
module m_riemann_solvers
2123

2224
! Dependencies =============================================================
@@ -34,10 +36,10 @@ module m_riemann_solvers
3436
implicit none
3537

3638
private; public :: s_initialize_riemann_solvers_module, &
37-
s_riemann_solver, &
38-
s_hll_riemann_solver, &
39-
s_hllc_riemann_solver, &
40-
s_finalize_riemann_solvers_module
39+
s_riemann_solver, &
40+
s_hll_riemann_solver, &
41+
s_hllc_riemann_solver, &
42+
s_finalize_riemann_solvers_module
4143

4244
abstract interface ! =======================================================
4345

@@ -108,20 +110,6 @@ module m_riemann_solvers
108110

109111
end subroutine s_abstract_riemann_solver
110112

111-
!> The abstract interface to the subroutines that are used to calculate
112-
!! the Roe and arithmetic average states. For more information refer to:
113-
!! 1) s_compute_roe_average_state
114-
!! 2) s_compute_arithmetic_average_state
115-
!! @param i First coordinate location index
116-
!! @param j Second coordinate location index
117-
!! @param k Third coordinate location index
118-
subroutine s_compute_abstract_average_state(qL_prim_rs_vf, qR_prim_rs_vf, i, j, k)
119-
import :: scalar_field, int_bounds_info, sys_size
120-
integer, intent(IN) :: i, j, k
121-
type(scalar_field), dimension(sys_size), intent(IN) :: qL_prim_rs_vf, qR_prim_rs_vf
122-
123-
end subroutine s_compute_abstract_average_state
124-
125113
!> The abstract interface to the subroutines that are utilized to compute
126114
!! the wave speeds of the Riemann problem either directly or by the means
127115
!! of pressure-velocity estimates. For more information please refer to:
@@ -328,11 +316,6 @@ module m_riemann_solvers
328316
!! Pointer to the procedure that is utilized to calculate either the HLL,
329317
!! HLLC or exact intercell fluxes, based on the choice of Riemann solver
330318

331-
procedure(s_compute_abstract_average_state), &
332-
pointer :: s_compute_average_state => null() !<
333-
!! Pointer to the subroutine utilized to calculate either the Roe or the
334-
!! arithmetic average state variables, based on the chosen average state
335-
336319
procedure(s_compute_abstract_wave_speeds), &
337320
pointer :: s_compute_wave_speeds => null() !<
338321
!! Pointer to the subroutine that is utilized to compute the wave speeds of
@@ -600,40 +583,8 @@ contains
600583
end if
601584
end do
602585
end if
603-
604-
if (avg_state == 2) then
605-
rho_avg = 5d-1*(rho_L + rho_R)
606-
607-
!$acc loop seq
608-
do i = 1, num_dims
609-
vel_avg(i) = 5d-1*(vel_L(i) + vel_R(i))
610-
end do
611-
612-
H_avg = 5d-1*(H_L + H_R)
613-
614-
gamma_avg = 5d-1*(gamma_L + gamma_R)
615-
elseif (avg_state == 1) then
616-
rho_avg = sqrt(rho_L*rho_R)
617-
618-
!$acc loop seq
619-
do i = 1, num_dims
620-
vel_avg(i) = (sqrt(rho_L)*vel_L(i) + sqrt(rho_R)*vel_R(i))/ &
621-
(sqrt(rho_L) + sqrt(rho_R))
622-
end do
623-
624-
H_avg = (sqrt(rho_L)*H_L + sqrt(rho_R)*H_R)/ &
625-
(sqrt(rho_L) + sqrt(rho_R))
626-
627-
gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/ &
628-
(sqrt(rho_L) + sqrt(rho_R))
629-
end if
630-
631-
vel_avg_rms = 0d0
632-
633-
!$acc loop seq
634-
do i = 1, num_dims
635-
vel_avg_rms = vel_avg_rms + vel_avg(i)**2d0
636-
end do
586+
587+
@:compute_average_state()
637588

638589
if (mixture_err) then
639590
if ((H_avg - 5d-1*vel_avg_rms) < 0d0) then
@@ -1210,35 +1161,8 @@ contains
12101161
12111162
H_L = (E_L + pres_L)/rho_L
12121163
H_R = (E_R + pres_R)/rho_R
1213-
if (avg_state == 2) then
1214-
1215-
rho_avg = 5d-1*(rho_L + rho_R)
1216-
vel_avg_rms = 0d0
1217-
!$acc loop seq
1218-
do i = 1, num_dims
1219-
vel_avg_rms = vel_avg_rms + (5d-1*(vel_L(i) + vel_R(i)))**2d0
1220-
end do
1221-
1222-
H_avg = 5d-1*(H_L + H_R)
1223-
1224-
gamma_avg = 5d-1*(gamma_L + gamma_R)
12251164
1226-
elseif (avg_state == 1) then
1227-
1228-
rho_avg = sqrt(rho_L*rho_R)
1229-
vel_avg_rms = 0d0
1230-
!$acc loop seq
1231-
do i = 1, num_dims
1232-
vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(i) + sqrt(rho_R)*vel_R(i))**2d0/ &
1233-
(sqrt(rho_L) + sqrt(rho_R))**2d0
1234-
end do
1235-
1236-
H_avg = (sqrt(rho_L)*H_L + sqrt(rho_R)*H_R)/ &
1237-
(sqrt(rho_L) + sqrt(rho_R))
1238-
1239-
gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/ &
1240-
(sqrt(rho_L) + sqrt(rho_R))
1241-
end if
1165+
@:compute_average_state()
12421166
12431167
if (mixture_err) then
12441168
if ((H_avg - 5d-1*vel_avg_rms) < 0d0) then
@@ -1545,39 +1469,8 @@ contains
15451469
15461470
H_L = (E_L + pres_L)/rho_L
15471471
H_R = (E_R + pres_R)/rho_R
1548-
if (avg_state == 2) then
1549-
1550-
rho_avg = 5d-1*(rho_L + rho_R)
1551-
!$acc loop seq
1552-
do i = 1, num_dims
1553-
vel_avg(i) = 5d-1*(vel_L(i) + vel_R(i))
1554-
end do
15551472
1556-
H_avg = 5d-1*(H_L + H_R)
1557-
1558-
gamma_avg = 5d-1*(gamma_L + gamma_R)
1559-
1560-
elseif (avg_state == 1) then
1561-
1562-
rho_avg = sqrt(rho_L*rho_R)
1563-
!$acc loop seq
1564-
do i = 1, num_dims
1565-
vel_avg(i) = (sqrt(rho_L)*vel_L(i) + sqrt(rho_R)*vel_R(i))/ &
1566-
(sqrt(rho_L) + sqrt(rho_R))
1567-
end do
1568-
1569-
H_avg = (sqrt(rho_L)*H_L + sqrt(rho_R)*H_R)/ &
1570-
(sqrt(rho_L) + sqrt(rho_R))
1571-
1572-
gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/ &
1573-
(sqrt(rho_L) + sqrt(rho_R))
1574-
end if
1575-
1576-
vel_avg_rms = 0d0
1577-
!$acc loop seq
1578-
do i = 1, num_dims
1579-
vel_avg_rms = vel_avg_rms + vel_avg(i)**2d0
1580-
end do
1473+
@:compute_average_state()
15811474
15821475
if (mixture_err) then
15831476
if ((H_avg - 5d-1*vel_avg_rms) < 0d0) then
@@ -1667,6 +1560,7 @@ contains
16671560
(pres_L - pres_R)/ &
16681561
(rho_avg*c_avg))
16691562
end if
1563+
16701564
! follows Einfeldt et al.
16711565
! s_M/P = min/max(0.,s_L/R)
16721566
s_M = min(0d0, s_L); s_P = max(0d0, s_R)
@@ -1692,44 +1586,32 @@ contains
16921586
16931587
! Momentum flux.
16941588
! f = \rho u u + p I, q = \rho u, q_star = \xi * \rho*(s_star, v, w)
1695-
if (bubbles .neqv. .true.) then
1589+
!$acc loop seq
1590+
do i = 1, num_dims
1591+
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) = &
1592+
xi_M*(rho_L*(vel_L(dir_idx(1))* &
1593+
vel_L(dir_idx(i)) + &
1594+
s_M*(xi_L*(dir_flg(dir_idx(i))*s_S + &
1595+
(1d0 - dir_flg(dir_idx(i)))* &
1596+
vel_L(dir_idx(i))) - vel_L(dir_idx(i)))) + &
1597+
dir_flg(dir_idx(i))*pres_L) &
1598+
+ xi_P*(rho_R*(vel_R(dir_idx(1))* &
1599+
vel_R(dir_idx(i)) + &
1600+
s_P*(xi_R*(dir_flg(dir_idx(i))*s_S + &
1601+
(1d0 - dir_flg(dir_idx(i)))* &
1602+
vel_R(dir_idx(i))) - vel_R(dir_idx(i)))) + &
1603+
dir_flg(dir_idx(i))*pres_R)
1604+
end do
1605+
1606+
if (bubbles) then
1607+
! Put p_tilde in
16961608
!$acc loop seq
16971609
do i = 1, num_dims
1698-
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) = &
1699-
xi_M*(rho_L*(vel_L(dir_idx(1))* &
1700-
vel_L(dir_idx(i)) + &
1701-
s_M*(xi_L*(dir_flg(dir_idx(i))*s_S + &
1702-
(1d0 - dir_flg(dir_idx(i)))* &
1703-
vel_L(dir_idx(i))) - vel_L(dir_idx(i)))) + &
1704-
dir_flg(dir_idx(i))*(pres_L)) &
1705-
+ xi_P*(rho_R*(vel_R(dir_idx(1))* &
1706-
vel_R(dir_idx(i)) + &
1707-
s_P*(xi_R*(dir_flg(dir_idx(i))*s_S + &
1708-
(1d0 - dir_flg(dir_idx(i)))* &
1709-
vel_R(dir_idx(i))) - vel_R(dir_idx(i)))) + &
1710-
dir_flg(dir_idx(i))*(pres_R))
1711-
! if (j==0) print*, 'flux_rs_vf', flux_rs_vf(cont_idx%end+dir_idx(i))%sf(j,k,l)
1610+
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) = &
1611+
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) + &
1612+
xi_M*(dir_flg(dir_idx(i))*(-1d0*ptilde_L)) &
1613+
+ xi_P*(dir_flg(dir_idx(i))*(-1d0*ptilde_R))
17121614
end do
1713-
else
1714-
! Include p_tilde
1715-
!$acc loop seq
1716-
do i = 1, num_dims
1717-
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) = &
1718-
xi_M*(rho_L*(vel_L(dir_idx(1))* &
1719-
vel_L(dir_idx(i)) + &
1720-
s_M*(xi_L*(dir_flg(dir_idx(i))*s_S + &
1721-
(1d0 - dir_flg(dir_idx(i)))* &
1722-
vel_L(dir_idx(i))) - vel_L(dir_idx(i)))) + &
1723-
dir_flg(dir_idx(i))*(pres_L - ptilde_L)) &
1724-
+ xi_P*(rho_R*(vel_R(dir_idx(1))* &
1725-
vel_R(dir_idx(i)) + &
1726-
s_P*(xi_R*(dir_flg(dir_idx(i))*s_S + &
1727-
(1d0 - dir_flg(dir_idx(i)))* &
1728-
vel_R(dir_idx(i))) - vel_R(dir_idx(i)))) + &
1729-
dir_flg(dir_idx(i))*(pres_R - ptilde_R))
1730-
! if (j==0) print*, 'flux_rs_vf', flux_rs_vf(cont_idx%end+dir_idx(i))%sf(j,k,l)
1731-
end do
1732-
17331615
end if
17341616
17351617
flux_rs${XYZ}$_vf(j, k, l, E_idx) = 0.d0
@@ -1906,9 +1788,6 @@ contains
19061788
end if
19071789
end do
19081790
1909-
!call s_comp_n_from_prim(qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + num_fluids), R0_L, nbub_L)
1910-
!call s_comp_n_from_prim(qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + num_fluids), R0_R, nbub_R)
1911-
19121791
nbub_L_denom = 0d0
19131792
nbub_R_denom = 0d0
19141793
@@ -2001,21 +1880,7 @@ contains
20011880
end do
20021881
20031882
elseif (avg_state == 1) then
2004-
2005-
rho_avg = sqrt(rho_L*rho_R)
2006-
2007-
vel_avg_rms = 0d0
2008-
!$acc loop seq
2009-
do i = 1, num_dims
2010-
vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(i) + sqrt(rho_R)*vel_R(i))**2d0/ &
2011-
(sqrt(rho_L) + sqrt(rho_R))**2d0
2012-
end do
2013-
2014-
H_avg = (sqrt(rho_L)*H_L + sqrt(rho_R)*H_R)/ &
2015-
(sqrt(rho_L) + sqrt(rho_R))
2016-
2017-
gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/ &
2018-
(sqrt(rho_L) + sqrt(rho_R))
1883+
call s_mpi_abort()
20191884
end if
20201885
20211886
if (mixture_err) then
@@ -2369,42 +2234,8 @@ contains
23692234
23702235
H_L = (E_L + pres_L)/rho_L
23712236
H_R = (E_R + pres_R)/rho_R
2372-
if (avg_state == 2) then
2373-
2374-
rho_avg = 5d-1*(rho_L + rho_R)
2375-
vel_avg_rms = (5d-1*(vel_L(1) + vel_R(1)))**2d0
2376-
if (num_dims >= 2) then
2377-
vel_avg_rms = vel_avg_rms + (5d-1*(vel_L(2) + vel_R(2)))**2d0
2378-
end if
2379-
if (num_dims == 3) then
2380-
vel_avg_rms = vel_avg_rms + (5d-1*(vel_L(3) + vel_R(3)))**2d0
2381-
end if
2382-
2383-
H_avg = 5d-1*(H_L + H_R)
2384-
2385-
gamma_avg = 5d-1*(gamma_L + gamma_R)
2386-
2387-
elseif (avg_state == 1) then
23882237
2389-
rho_avg = sqrt(rho_L*rho_R)
2390-
vel_avg_rms = (sqrt(rho_L)*vel_L(1) + sqrt(rho_R)*vel_R(1))**2d0/ &
2391-
(sqrt(rho_L) + sqrt(rho_R))**2d0
2392-
2393-
if (num_dims >= 2) then
2394-
vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(2) + sqrt(rho_R)*vel_R(2))**2d0/ &
2395-
(sqrt(rho_L) + sqrt(rho_R))**2d0
2396-
end if
2397-
if (num_dims == 3) then
2398-
vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(3) + sqrt(rho_R)*vel_R(3))**2d0/ &
2399-
(sqrt(rho_L) + sqrt(rho_R))**2d0
2400-
end if
2401-
2402-
H_avg = (sqrt(rho_L)*H_L + sqrt(rho_R)*H_R)/ &
2403-
(sqrt(rho_L) + sqrt(rho_R))
2404-
2405-
gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/ &
2406-
(sqrt(rho_L) + sqrt(rho_R))
2407-
end if
2238+
@:compute_average_state()
24082239
24092240
if (mixture_err) then
24102241
if ((H_avg - 5d-1*vel_avg_rms) < 0d0) then
@@ -2468,6 +2299,7 @@ contains
24682299
(s_R - vel_R(idx1))) &
24692300
/(rho_L*(s_L - vel_L(idx1)) - &
24702301
rho_R*(s_R - vel_R(idx1)))
2302+
24712303
elseif (wave_speeds == 2) then
24722304
pres_SL = 5d-1*(pres_L + pres_R + rho_avg*c_avg* &
24732305
(vel_L(idx1) - &
@@ -4602,7 +4434,7 @@ contains
46024434

46034435
! Disassociating the procedural pointers to the procedures that were
46044436
! utilized to compute the average state and estimate the wave speeds
4605-
s_compute_average_state => null(); s_compute_wave_speeds => null()
4437+
s_compute_wave_speeds => null()
46064438

46074439
! Disassociating procedural pointer to the subroutine which was
46084440
! utilized to calculate the viscous source flux

0 commit comments

Comments
 (0)