1!-------------------------------------------------------------------------------
2
3!VERS
4
5! This file is part of Code_Saturne, a general-purpose CFD tool.
6!
7! Copyright (C) 1998-2021 EDF S.A.
8!
9! This program is free software; you can redistribute it and/or modify it under
10! the terms of the GNU General Public License as published by the Free Software
11! Foundation; either version 2 of the License, or (at your option) any later
12! version.
13!
14! This program is distributed in the hope that it will be useful, but WITHOUT
15! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
17! details.
18!
19! You should have received a copy of the GNU General Public License along with
20! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
21! Street, Fifth Floor, Boston, MA 02110-1301, USA.
22
23!-------------------------------------------------------------------------------
24
25!===============================================================================
26! Function:
27! ---------
28
29! Example of cs_f_user_boundary_conditions subroutine.f90 for inlet
30! with automatic inlet profile.
31
32! This example assumes the mesh is orthogonal at the inlet.
33!
34!-------------------------------------------------------------------------------
35
36!-------------------------------------------------------------------------------
37! Arguments
38!______________________________________________________________________________.
39!  mode           name          role                                           !
40!______________________________________________________________________________!
41!> \param[in]     nvar          total number of variables
42!> \param[in]     nscal         total number of scalars
43!> \param[out]    icodcl        boundary condition code:
44!>                               - 1 Dirichlet
45!>                               - 2 Radiative outlet
46!>                               - 3 Neumann
47!>                               - 4 sliding and
48!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
49!>                               - 5 smooth wall and
50!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
51!>                               - 6 rough wall and
52!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
53!>                               - 9 free inlet/outlet
54!>                                 (input mass flux blocked to 0)
55!>                               - 13 Dirichlet for the advection operator and
56!>                                    Neumann for the diffusion operator
57!> \param[in]     itrifb        indirection for boundary faces ordering
58!> \param[in,out] itypfb        boundary face types
59!> \param[in,out] izfppp        boundary face zone number
60!> \param[in]     dt            time step (per cell)
61!> \param[in,out] rcodcl        boundary condition values:
62!>                               - rcodcl(1) value of the dirichlet
63!>                               - rcodcl(2) value of the exterior exchange
64!>                                 coefficient (infinite if no exchange)
65!>                               - rcodcl(3) value flux density
66!>                                 (negative if gain) in w/m2
67!>                                 -# for the velocity \f$ (\mu+\mu_T)
68!>                                    \gradt \, \vect{u} \cdot \vect{n}  \f$
69!>                                 -# for the pressure \f$ \Delta t
70!>                                    \grad P \cdot \vect{n}  \f$
71!>                                 -# for a scalar \f$ cp \left( K +
72!>                                     \dfrac{K_T}{\sigma_T} \right)
73!>                                     \grad T \cdot \vect{n} \f$
74!_______________________________________________________________________________
75
76subroutine cs_f_user_boundary_conditions &
77 ( nvar   , nscal  ,                                              &
78   icodcl , itrifb , itypfb , izfppp ,                            &
79   dt     ,                                                       &
80   rcodcl )
81
82!===============================================================================
83
84!===============================================================================
85! Module files
86!===============================================================================
87
88use paramx
89use numvar
90use optcal
91use cstphy
92use cstnum
93use entsor
94use parall
95use period
96use ppppar
97use ppthch
98use coincl
99use cpincl
100use ppincl
101use ppcpfu
102use atincl
103use atsoil
104use ctincl
105use cs_fuel_incl
106use mesh
107use field
108use turbomachinery
109use cs_c_bindings
110
111!===============================================================================
112
113implicit none
114
115! Arguments
116
117integer          nvar   , nscal
118
119integer          icodcl(nfabor,nvar)
120integer          itrifb(nfabor), itypfb(nfabor)
121integer          izfppp(nfabor)
122
123double precision dt(ncelet)
124double precision rcodcl(nfabor,nvar,3)
125
126! Local variables
127
128!< [loc_var_dec]
129integer          ifac, iel, ii, ilelt, nlelt
130
131double precision xustar2, xdh, d2s3, rhomoy
132double precision acc(2), fmprsc, fmul, uref2, vnrm
133
134integer, allocatable, dimension(:) :: lstelt, mrkcel
135
136double precision, dimension(:), pointer :: bfpro_rom
137double precision, dimension(:), pointer :: cvar_k, cvar_ep, cvar_phi
138double precision, dimension(:), pointer :: cvar_omg
139double precision, dimension(:), pointer :: cvar_al, cvar_fb
140double precision, dimension(:), pointer :: cvar_nusa
141double precision, dimension(:), pointer :: cvar_scal
142double precision, dimension(:,:), pointer :: cvar_vel
143double precision, dimension(:,:), pointer :: cvar_rij
144!< [loc_var_dec]
145
146!===============================================================================
147! Initialization
148!===============================================================================
149
150!< [init]
151allocate(lstelt(nfabor))  ! temporary array for boundary faces selection
152
153! Map field arrays
154call field_get_val_v(ivarfl(iu), cvar_vel)
155call field_get_val_s(ibrom, bfpro_rom)
156
157if (itytur.eq.2) then
158
159  call field_get_val_s(ivarfl(ik), cvar_k)
160  call field_get_val_s(ivarfl(iep), cvar_ep)
161
162elseif (itytur.eq.3) then
163  call field_get_val_v(ivarfl(irij), cvar_rij)
164  call field_get_val_s(ivarfl(iep), cvar_ep)
165
166  if (iturb.eq.32) then
167    call field_get_val_s(ivarfl(ial), cvar_al)
168  endif
169
170elseif (itytur.eq.5) then
171
172  call field_get_val_s(ivarfl(ik), cvar_k)
173  call field_get_val_s(ivarfl(iep), cvar_ep)
174  call field_get_val_s(ivarfl(iphi), cvar_phi)
175
176  if (iturb.eq.50) then
177    call field_get_val_s(ivarfl(ifb), cvar_fb)
178  elseif (iturb.eq.51) then
179    call field_get_val_s(ivarfl(ial), cvar_al)
180  endif
181
182elseif (iturb.eq.60) then
183
184  call field_get_val_s(ivarfl(ik), cvar_k)
185  call field_get_val_s(ivarfl(iomg), cvar_omg)
186
187elseif (iturb.eq.70) then
188
189  call field_get_val_s(ivarfl(inusa), cvar_nusa)
190
191endif
192
193d2s3 = 2.d0/3.d0
194!< [init]
195
196!===============================================================================
197! Assign a pseudo-periodic channel type inlet to a set of boundary faces.
198
199! For each subset:
200! - use selection criteria to filter boundary faces of a given subset
201! - loop on faces from a subset
202!   - set the boundary condition for each face
203
204! A feedback loop is used so as to progressively reach a state similar
205! to that of a periodic channel at the inlet.
206!===============================================================================
207
208!< [example_1]
209call getfbr('INLET', nlelt, lstelt)
210!==========
211
212fmprsc = 1.d0 ! mean prescribed velocity
213
214if (ntcabs.eq.1) then
215
216  ! For the Rij-EBRSM model (and possibly V2f), we need a non-flat profile,
217  ! so as to ensure turbulent production, and avoid laminarization;
218  ! here, we simply divide the initial velocity by 10 for inlet
219  ! faces adjacent to the wall.
220
221  ! The loop below assumes wall conditions have been defined first
222  ! (in the GUI, or in this file, before the current test).
223
224  if (iturb.eq.32 .or. itytur.eq.5) then
225
226    allocate(mrkcel(ncelet))
227    do iel = 1, ncelet
228      mrkcel(iel) = 0
229    enddo
230
231    do ifac = 1, nfabor
232      if (itypfb(ifac) .eq. iparoi) then
233        iel = ifabor(ifac)
234        mrkcel(iel) = 1
235      endif
236    enddo
237
238  endif
239
240  do ilelt = 1, nlelt
241
242    ifac = lstelt(ilelt)
243    iel = ifabor(ifac)
244
245    itypfb(ifac) = ientre
246
247    rcodcl(ifac,iu,1) = - fmprsc * surfbo(1,ifac) / surfbn(ifac)
248    rcodcl(ifac,iv,1) = - fmprsc * surfbo(2,ifac) / surfbn(ifac)
249    rcodcl(ifac,iw,1) = - fmprsc * surfbo(3,ifac) / surfbn(ifac)
250
251    if (iturb.eq.32 .or. itytur.eq.5) then
252      if (mrkcel(iel) .eq. 1) then
253        rcodcl(ifac,iu,1) = fmprsc/10.d0
254      endif
255    endif
256
257    uref2 = rcodcl(ifac,iu,1)**2  &
258          + rcodcl(ifac,iv,1)**2  &
259          + rcodcl(ifac,iw,1)**2
260    uref2 = max(uref2,1.d-12)
261
262    !   Turbulence example computed using equations valid for a pipe.
263
264    !   We will be careful to specify a hydraulic diameter adapted
265    !     to the current inlet.
266
267    !   We will also be careful if necessary to use a more precise
268    !     formula for the dynamic viscosity use in the calculation of
269    !     the Reynolds number (especially if it is variable, it may be
270    !     useful to take the law from 'usphyv'. Here, we use by default
271    !     the 'viscl0" value.
272    !   Regarding the density, we have access to its value at boundary
273    !     faces (romb) so this value is the one used here (specifically,
274    !     it is consistent with the processing in 'usphyv', in case of
275    !     variable density)
276
277    !     Hydraulic diameter
278    xdh     = 1.d0
279
280    ! Calculation of turbulent inlet conditions using
281    !   the turbulence intensity and standard laws for a circular pipe
282    !   (their initialization is not needed here but is good practice)
283
284    rhomoy  = bfpro_rom(ifac)
285
286    call turbulence_bc_inlet_hyd_diam(ifac, uref2, xdh, rhomoy, viscl0,  &
287                                      rcodcl)
288
289    ! Handle scalars
290    if (nscal.gt.0) then
291      do ii = 1, nscal
292        rcodcl(ifac,isca(ii),1) = 1.d0
293      enddo
294    endif
295
296  enddo
297
298  if (iturb.eq.32 .or. itytur.eq.5) then
299    deallocate(mrkcel)
300  endif
301
302else
303
304! Subsequent time steps
305!----------------------
306
307  acc(1) = 0.d0
308  acc(2) = 0.d0
309
310  ! Estimate multiplier
311
312  do ilelt = 1, nlelt
313
314    ifac = lstelt(ilelt)
315    iel = ifabor(ifac)
316
317    vnrm = sqrt(cvar_vel(1,iel)**2 + cvar_vel(2,iel)**2 + cvar_vel(3,iel)**2)
318    acc(1) = acc(1) + vnrm*surfbn(ifac)
319    acc(2) = acc(2) + surfbn(ifac)
320
321  enddo
322
323  if (irangp.ge.0) then
324    call parrsm(2, acc)
325  endif
326
327  if (acc(1).le.epzero) then
328     fmul = 0.d0 ! zero velocity in bulk domain
329  else
330     fmul = fmprsc/(acc(1)/acc(2)) ! 1 / estimate flow multiplier
331  endif
332
333  ! Apply BC
334
335  do ilelt = 1, nlelt
336
337    ifac = lstelt(ilelt)
338    iel = ifabor(ifac)
339
340    itypfb(ifac) = ientre
341
342    vnrm = sqrt(cvar_vel(1,iel)**2 + cvar_vel(2,iel)**2 + cvar_vel(3,iel)**2)
343
344    rcodcl(ifac,iu,1) = - fmul * vnrm * surfbo(1,ifac) / surfbn(ifac)
345    rcodcl(ifac,iv,1) = - fmul * vnrm * surfbo(2,ifac) / surfbn(ifac)
346    rcodcl(ifac,iw,1) = - fmul * vnrm * surfbo(3,ifac) / surfbn(ifac)
347
348    if (itytur.eq.2) then
349
350      rcodcl(ifac,ik,1)  = cvar_k(iel)
351      rcodcl(ifac,iep,1) = cvar_ep(iel)
352
353    elseif (itytur.eq.3) then
354
355        rcodcl(ifac,ir11,1) = cvar_rij(1,iel)
356        rcodcl(ifac,ir22,1) = cvar_rij(2,iel)
357        rcodcl(ifac,ir33,1) = cvar_rij(3,iel)
358        rcodcl(ifac,ir12,1) = cvar_rij(4,iel)
359        rcodcl(ifac,ir13,1) = cvar_rij(6,iel)
360        rcodcl(ifac,ir23,1) = cvar_rij(5,iel)
361
362      rcodcl(ifac,iep,1)  = cvar_ep(iel)
363
364      if (iturb.eq.32) then
365        rcodcl(ifac,ial,1)  = cvar_al(iel)
366      endif
367
368    elseif (itytur.eq.5) then
369
370      rcodcl(ifac,ik,1)  = cvar_k(iel)
371      rcodcl(ifac,iep,1) = cvar_ep(iel)
372      rcodcl(ifac,iphi,1) = cvar_phi(iel)
373
374      if (iturb.eq.50) then
375        rcodcl(ifac,ifb,1)  = cvar_fb(iel)
376      elseif (iturb.eq.51) then
377        rcodcl(ifac,ial,1)  = cvar_al(iel)
378      endif
379
380    elseif (iturb.eq.60) then
381
382      rcodcl(ifac,ik,1)  = cvar_k(iel)
383      rcodcl(ifac,iomg,1) = cvar_omg(iel)
384
385    elseif (iturb.eq.70) then
386
387      rcodcl(ifac,inusa,1) = cvar_nusa(iel)
388
389    endif
390
391    ! Handle scalars (a correction similar to that of velocity is suggested
392    !                 rather than the simpler code below)
393    if (nscal.gt.0) then
394      do ii = 1, nscal
395        call field_get_val_s(ivarfl(isca(ii)), cvar_scal)
396        rcodcl(ifac,isca(ii),1) = cvar_scal(iel)
397      enddo
398    endif
399
400  enddo
401
402endif
403!< [example_1]
404
405!--------
406! Formats
407!--------
408
409!----
410! End
411!----
412
413deallocate(lstelt)  ! temporary array for boundary faces selection
414
415return
416end subroutine cs_f_user_boundary_conditions
417