1!-------------------------------------------------------------------------------
2
3! This file is part of Code_Saturne, a general-purpose CFD tool.
4!
5! Copyright (C) 1998-2021 EDF S.A.
6!
7! This program is free software; you can redistribute it and/or modify it under
8! the terms of the GNU General Public License as published by the Free Software
9! Foundation; either version 2 of the License, or (at your option) any later
10! version.
11!
12! This program is distributed in the hope that it will be useful, but WITHOUT
13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
15! details.
16!
17! You should have received a copy of the GNU General Public License along with
18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
19! Street, Fifth Floor, Boston, MA 02110-1301, USA.
20
21!-------------------------------------------------------------------------------
22
23subroutine cfrusb &
24!================
25
26 ( ifac   ,                                                       &
27   bc_en  , bc_pr  , bc_vel    )
28
29!===============================================================================
30! FUNCTION :
31! ---------
32
33! Rusanov flux at the boundary for Euler + Energy
34
35! d rho   /dt + div rho u             = 0
36! d rho u /dt + div rho u u + grad  P = 0
37! d E     /dt + div rho u E + div u P = 0
38
39!-------------------------------------------------------------------------------
40! Arguments
41!__________________.____._____.________________________________________________.
42! name             !type!mode ! role                                           !
43!__________________!____!_____!________________________________________________!
44! ifac             ! i  ! <-- ! face number                                    !
45!__________________!____!_____!________________________________________________!
46
47!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
48!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
49!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
50!            --- tableau de travail
51!===============================================================================
52
53!===============================================================================
54! Module files
55!===============================================================================
56
57use paramx
58use pointe, only:rvoid1
59use numvar
60use optcal
61use cstphy
62use cstnum
63use parall
64use entsor
65use ppppar
66use ppthch
67use ppincl
68use cfpoin
69use mesh
70use field
71use cs_cf_bindings
72
73!===============================================================================
74
75implicit none
76
77! Arguments
78
79integer          ifac
80
81double precision bc_en(nfabor), bc_pr(nfabor), bc_vel(3,nfabor)
82
83! Local variables
84
85integer          iel, l_size
86integer          ien
87
88double precision inv_surfbn, rnx, rny, rnz
89double precision b_vel_n, c_vel_n, b_masfl, c_masfl, b_c, c_c, b_c2(1), c_c2(1)
90double precision rrus, r_b_masfl, r_b_vel_n
91double precision, dimension(:,:), pointer :: cofacv
92double precision, dimension(:), pointer :: coface
93double precision, dimension(:), pointer :: crom, brom
94double precision, dimension(:,:), pointer :: vel
95double precision, dimension(:), pointer :: cvar_pr, cvar_en, cpro_cp, cpro_cv
96double precision, dimension(:), pointer :: rvoid
97
98double precision c_cp(1), c_cv(1)
99
100!===============================================================================
101
102rvoid => null()
103
104!===============================================================================
105! 0. INITIALISATION
106!===============================================================================
107
108! Map field arrays
109call field_get_val_v(ivarfl(iu), vel)
110
111ien = isca(ienerg)
112
113call field_get_val_s(icrom, crom)
114call field_get_val_s(ibrom, brom)
115
116call field_get_val_s(ivarfl(ipr), cvar_pr)
117call field_get_val_s(ivarfl(ien), cvar_en)
118
119call field_get_coefac_v(ivarfl(iu), cofacv)
120call field_get_coefac_s(ivarfl(ien), coface)
121
122iel   = ifabor(ifac)
123
124! Initialize local specific heat values and handle uniform cases
125
126if (icp.ge.0) then
127  call field_get_val_s(icp, cpro_cp)
128  c_cp(1) = cpro_cp(iel)
129else
130  cpro_cp => rvoid1
131  c_cp(1) = 0.d0
132endif
133
134if (icv.ge.0) then
135  call field_get_val_s(icv, cpro_cv)
136  c_cv(1) = cpro_cv(iel)
137else
138  cpro_cv => rvoid1
139  c_cv(1) = 0.d0
140endif
141
142!===============================================================================
143! 1. COMPUTE VALUES NEEDED FOR RUSANOV SCHEME
144!===============================================================================
145
146inv_surfbn = 1. / surfbn(ifac)
147rnx = surfbo(1,ifac) * inv_surfbn
148rny = surfbo(2,ifac) * inv_surfbn
149rnz = surfbo(3,ifac) * inv_surfbn
150
151b_vel_n = bc_vel(1,ifac)*rnx + bc_vel(2,ifac)*rny + bc_vel(3,ifac)*rnz
152c_vel_n = vel(1,iel)*rnx + vel(2,iel)*rny + vel(3,iel)*rnz
153b_masfl  = brom(ifac)*b_vel_n
154c_masfl  = crom(iel)*c_vel_n
155
156l_size = 1
157
158call cs_cf_thermo_c_square(c_cp, c_cv, &
159                           bc_pr(ifac:ifac), brom(ifac:ifac), &
160                           rvoid, rvoid, rvoid, b_c2, l_size)
161call cs_cf_thermo_c_square(c_cp, c_cv, &
162                           cvar_pr(iel:iel), crom(iel:iel), &
163                           rvoid, rvoid, rvoid, c_c2, l_size)
164
165b_c    = sqrt(b_c2(1))
166c_c    = sqrt(c_c2(1))
167rrus  = max(abs(b_vel_n)+b_c, abs(c_vel_n)+c_c)
168
169! boundary mass flux computed with Rusanov scheme
170r_b_masfl  = 0.5d0*(b_masfl+c_masfl) - 0.5d0*rrus*(brom(ifac)-crom(iel))
171
172! Rusanov normal velocity (computed using boundary density)
173r_b_vel_n = r_b_masfl / brom(ifac)
174
175! Update velocity boundary condition
176bc_vel(1,ifac) = bc_vel(1,ifac) + (r_b_vel_n - b_vel_n)*rnx
177bc_vel(2,ifac) = bc_vel(2,ifac) + (r_b_vel_n - b_vel_n)*rny
178bc_vel(3,ifac) = bc_vel(3,ifac) + (r_b_vel_n - b_vel_n)*rnz
179
180!===============================================================================
181! 2. CONVECTIVE RUSANOV FLUX
182!===============================================================================
183
184! Tag the faces where a Rusanov flux is computed
185! The tag will be used in bilsc2 to retrieve the faces where a Rusanov flux
186! has to be imposed
187icvfli(ifac) = 1
188
189! Momentum flux (the centered pressure contribution is directly taken into account
190! in the pressure BC)
191cofacv(1,ifac) = suffbn(ifac)*                                                  &
192                 0.5d0*( b_masfl*bc_vel(1,ifac) + c_masfl*vel(1,iel)            &
193                        -rrus*(brom(ifac)*bc_vel(1,ifac)                        &
194                        -crom(iel)*vel(1,iel)) )
195
196cofacv(2,ifac) = suffbn(ifac)*                                                  &
197                 0.5d0*( b_masfl*bc_vel(2,ifac) + c_masfl*vel(2,iel)            &
198                        -rrus*( brom(ifac)*bc_vel(2,ifac)                       &
199                        -crom(iel)*vel(2,iel)) )
200
201cofacv(3,ifac) = suffbn(ifac)*                                                  &
202                 0.5d0*( b_masfl*bc_vel(3,ifac) + c_masfl*vel(3,iel)            &
203                        -rrus*(brom(ifac)*bc_vel(3,ifac)                        &
204                        -crom(iel)*vel(3,iel)) )
205
206! BC for the pressure gradient in the momentum balance
207bc_pr(ifac) = 0.5d0 * (bc_pr(ifac) + cvar_pr(iel))
208
209! Total energy flux
210coface(ifac) = suffbn(ifac)*                                                    &
211               0.5d0*( b_masfl*bc_en(ifac) + c_masfl*cvar_en(iel)               &
212                      +b_vel_n*bc_pr(ifac) + c_vel_n*cvar_pr(iel)               &
213                      -rrus*(brom(ifac)*bc_en(ifac)                             &
214                      -crom(iel)*cvar_en(iel)) )
215
216return
217
218end subroutine
219