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
23!> \file post_util.f90
24
25!===============================================================================
26! Function:
27! ---------
28
29!> \brief Compute thermal flux at boundary.
30
31!> If working with enthalpy, compute an enthalpy flux.
32
33!-------------------------------------------------------------------------------
34! Arguments
35!______________________________________________________________________________.
36!  mode           name          role                                           !
37!______________________________________________________________________________!
38!> \param[in]     nfbrps        number of boundary faces to postprocess
39!> \param[in]     lstfbr        list of boundary faces to postprocess
40!> \param[out]    bflux         boundary heat flux at selected faces
41!_______________________________________________________________________________
42
43subroutine post_boundary_thermal_flux &
44 ( nfbrps , lstfbr ,                                              &
45   bflux )
46
47!===============================================================================
48
49!===============================================================================
50! Module files
51!===============================================================================
52
53use paramx
54use pointe
55use entsor
56use cstnum
57use cstphy
58use optcal
59use numvar
60use parall
61use period
62use mesh
63use field
64use field_operator
65use cs_c_bindings
66
67!===============================================================================
68
69implicit none
70
71! Arguments
72
73integer, intent(in)                                :: nfbrps
74integer, dimension(nfbrps), intent(in)             :: lstfbr
75double precision, dimension(nfbrps), intent(out)   :: bflux
76
77! Local variables
78
79integer ::         f_id
80integer ::         iloc, ivar
81
82character(len=80) :: f_name
83integer(c_int), dimension(:), allocatable :: c_faces
84
85!===============================================================================
86! Interfaces
87!===============================================================================
88
89interface
90
91  subroutine cs_post_boundary_flux(scalar_name, n_loc_b_faces, b_face_ids,   &
92                                   b_face_flux)                              &
93    bind(C, name='cs_post_boundary_flux')
94    use, intrinsic :: iso_c_binding
95    implicit none
96    character(kind=c_char, len=1), dimension(*), intent(in) :: scalar_name
97    integer(c_int), value :: n_loc_b_faces
98    integer(c_int), dimension(*) :: b_face_ids
99    real(kind=c_double), dimension(*) :: b_face_flux
100  end subroutine cs_post_boundary_flux
101
102end interface
103
104!===============================================================================
105
106! Initialize variables to avoid compiler warnings
107
108if (iscalt.gt.0) then
109
110  ivar = isca(iscalt)
111  f_id = ivarfl(ivar)
112
113  call field_get_name (f_id, f_name)
114
115  allocate(c_faces(nfbrps))
116
117  do iloc = 1, nfbrps
118    c_faces(iloc) = lstfbr(iloc) - 1
119  enddo
120
121  call cs_post_boundary_flux(trim(f_name)//c_null_char, nfbrps, c_faces, bflux)
122
123  deallocate(c_faces)
124
125else ! if thermal variable is not available
126
127  do iloc = 1, nfbrps
128    bflux(iloc) = 0.d0
129  enddo
130
131endif
132
133!----
134! End
135!----
136
137return
138end subroutine post_boundary_thermal_flux
139
140!===============================================================================
141! Function:
142! ---------
143
144!> \brief Compute Nusselt number near boundary.
145
146!-------------------------------------------------------------------------------
147! Arguments
148!______________________________________________________________________________.
149!  mode           name          role                                           !
150!______________________________________________________________________________!
151!> \param[in]     nfbrps        number of boundary faces to postprocess
152!> \param[in]     lstfbr        list of boundary faces to postprocess
153!> \param[out]    bnussl        Nusselt near boundary
154!_______________________________________________________________________________
155
156subroutine post_boundary_nusselt &
157 ( nfbrps , lstfbr ,                                              &
158   bnussl )
159
160!===============================================================================
161
162!===============================================================================
163! Module files
164!===============================================================================
165
166use paramx
167use pointe
168use entsor
169use cstnum
170use cstphy
171use optcal
172use numvar
173use parall
174use period
175use mesh
176use field
177use field_operator
178use cs_c_bindings
179
180!===============================================================================
181
182implicit none
183
184! Arguments
185
186integer, intent(in)                                :: nfbrps
187integer, dimension(nfbrps), intent(in)             :: lstfbr
188double precision, dimension(nfbrps), intent(out)   :: bnussl
189
190! Local variables
191
192integer ::         inc, iccocg
193integer ::         iel, ifac, iloc, ivar
194integer ::         ifcvsl, itplus, itstar
195
196double precision :: xvsl  , srfbn , heq
197double precision :: diipbx, diipby, diipbz
198double precision :: numer, denom, visls_0
199
200double precision, allocatable, dimension(:) :: theipb
201double precision, allocatable, dimension(:,:) :: grad
202double precision, dimension(:), pointer :: cviscl
203double precision, dimension(:), pointer :: coefap, coefbp, cofafp, cofbfp
204double precision, dimension(:), pointer :: tplusp, tstarp
205double precision, dimension(:), pointer :: tscalp, hextp, hintp, dist_theipb
206
207type(var_cal_opt) :: vcopt
208
209logical(c_bool), dimension(:), pointer ::  cpl_faces
210
211!===============================================================================
212
213! pointers to T+ and T* if saved
214
215call field_get_id_try('tplus', itplus)
216call field_get_id_try('tstar', itstar)
217
218if (itstar.ge.0 .and. itplus.ge.0) then
219
220  ivar = isca(iscalt)
221
222  call field_get_val_prev_s(ivarfl(ivar), tscalp)
223
224  call field_get_val_s(itplus, tplusp)
225  call field_get_val_s(itstar, tstarp)
226
227  ! Boundary condition pointers for diffusion
228
229  call field_get_coefaf_s(ivarfl(ivar), cofafp)
230  call field_get_coefbf_s(ivarfl(ivar), cofbfp)
231
232  ! Boundary condition pointers for diffusion with coupling
233
234  call field_get_hext(ivarfl(ivar), hextp)
235  call field_get_hint(ivarfl(ivar), hintp)
236
237  ! Compute variable values at boundary faces
238
239  call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt)
240
241  allocate(theipb(nfabor))
242  theipb = 0.d0
243
244  do iloc = 1, nfbrps
245    ifac = lstfbr(iloc)
246    iel = ifabor(ifac)
247    theipb(ifac) = tscalp(iel)
248  enddo
249
250  ! Reconstructed fluxes
251  if (vcopt%ircflu .gt. 0 .and. itbrrb.eq.1) then
252
253    ! Compute gradient of temperature / enthalpy
254
255    allocate(grad(3,ncelet))
256
257    inc = 1
258    iccocg = 1
259
260    call field_gradient_scalar(ivarfl(ivar), 0, 0, inc, iccocg, grad)
261
262    ! Compute reconstructed temperature
263
264    do iloc = 1, nfbrps
265      ifac = lstfbr(iloc)
266      iel = ifabor(ifac)
267
268      diipbx = diipb(1,ifac)
269      diipby = diipb(2,ifac)
270      diipbz = diipb(3,ifac)
271
272      theipb(ifac) = theipb(ifac)                                               &
273                   + diipbx*grad(1,iel) + diipby*grad(2,iel) + diipbz*grad(3,iel)
274    enddo
275
276    deallocate(grad)
277  endif
278
279  if (vcopt%icoupl.gt.0) then
280    call field_get_coupled_faces(ivarfl(ivar), cpl_faces)
281    allocate(dist_theipb(nfabor))
282    call cs_ic_field_dist_data_by_face_id(ivarfl(ivar), 1, theipb, dist_theipb)
283  endif
284
285  ! Physical property pointers
286
287  call field_get_key_int (ivarfl(ivar), kivisl, ifcvsl)
288  if (ifcvsl .ge. 0) then
289    call field_get_val_s(ifcvsl, cviscl)
290    visls_0 = -1
291  else
292    call field_get_key_double(ivarfl(ivar), kvisl0, visls_0)
293  endif
294
295  ! Boundary condition pointers for gradients and advection
296
297  call field_get_coefa_s(ivarfl(ivar), coefap)
298  call field_get_coefb_s(ivarfl(ivar), coefbp)
299
300  ! Compute using reconstructed temperature value in boundary cells
301
302  do iloc = 1, nfbrps
303
304    ifac = lstfbr(iloc)
305    iel = ifabor(ifac)
306
307    if (ifcvsl.ge.0) then
308      xvsl = cviscl(iel)
309    else
310      xvsl = visls_0
311    endif
312
313    srfbn = surfbn(ifac)
314
315    numer = (cofafp(ifac) + cofbfp(ifac)*theipb(ifac)) * distb(ifac)
316    ! here numer = 0 if current face is coupled
317
318    if (vcopt%icoupl.gt.0.and.ntcabs.gt.ntpabs) then
319      ! FIXME exchange coefs not computed at start of calculation
320      if (cpl_faces(ifac)) then
321        heq = hextp(ifac) * hintp(ifac) / ((hextp(ifac) + hintp(ifac))*srfbn)
322        numer = heq*(theipb(ifac)-dist_theipb(ifac)) * distb(ifac)
323      endif
324    endif
325
326    denom = xvsl * tplusp(ifac)*tstarp(ifac)
327
328    if (abs(denom).gt.1e-30) then
329      bnussl(iloc) = numer / denom
330    else
331      bnussl(iloc) = 0.d0
332    endif
333
334  enddo
335
336  if (vcopt%icoupl.gt.0) then
337    deallocate(dist_theipb)
338  endif
339
340  deallocate(theipb)
341
342else ! default if not computable
343
344  do iloc = 1, nfbrps
345    bnussl(iloc) = -1.d0
346  enddo
347
348endif
349
350!--------
351! Formats
352!--------
353
354!----
355! End
356!----
357
358return
359end subroutine post_boundary_nusselt
360
361!===============================================================================
362! Function:
363! ---------
364
365!> \brief Compute stress at boundary.
366
367!-------------------------------------------------------------------------------
368! Arguments
369!______________________________________________________________________________.
370!  mode           name          role                                           !
371!______________________________________________________________________________!
372!> \param[in]     nfbrps        number of boundary faces to postprocess
373!> \param[in]     lstfbr        list of boundary faces to postprocess
374!> \param[out]    stress        stress at selected faces
375!_______________________________________________________________________________
376
377subroutine post_stress &
378 ( nfbrps , lstfbr ,                                                    &
379   stress )
380
381!===============================================================================
382
383!===============================================================================
384! Module files
385!===============================================================================
386
387use paramx
388use pointe
389use entsor
390use cstnum
391use optcal
392use numvar
393use parall
394use period
395use mesh
396use field
397
398!===============================================================================
399
400implicit none
401
402! Arguments
403
404integer, intent(in)                                 :: nfbrps
405integer, dimension(nfbrps), intent(in)              :: lstfbr
406double precision, dimension(3, nfbrps), intent(out) :: stress
407
408! Local variables
409
410integer          :: ifac  , iloc
411double precision :: srfbn
412double precision, dimension(:,:), pointer :: forbr
413
414!===============================================================================
415
416call field_get_val_v(iforbr, forbr)
417
418do iloc = 1, nfbrps
419  ifac = lstfbr(iloc)
420  srfbn = surfbn(ifac)
421  stress(1,iloc) = forbr(1,ifac)/srfbn
422  stress(2,iloc) = forbr(2,ifac)/srfbn
423  stress(3,iloc) = forbr(3,ifac)/srfbn
424enddo
425
426!--------
427! Formats
428!--------
429
430!----
431! End
432!----
433
434return
435end subroutine post_stress
436
437!===============================================================================
438! Function:
439! ---------
440
441!> \brief Extract stress normal to the boundary.
442
443!-------------------------------------------------------------------------------
444! Arguments
445!______________________________________________________________________________.
446!  mode           name          role                                           !
447!______________________________________________________________________________!
448!> \param[in]     nfbrps        number of boundary faces to postprocess
449!> \param[in]     lstfbr        list of boundary faces to postprocess
450!> \param[out]    effnrm        stress normal to wall at selected faces
451!_______________________________________________________________________________
452
453subroutine post_stress_normal &
454 ( nfbrps , lstfbr ,                                              &
455   effnrm )
456
457!===============================================================================
458
459!===============================================================================
460! Module files
461!===============================================================================
462
463use paramx
464use pointe
465use entsor
466use cstnum
467use optcal
468use numvar
469use parall
470use period
471use mesh
472use field
473
474!===============================================================================
475
476implicit none
477
478! Arguments
479
480integer, intent(in)                                 :: nfbrps
481integer, dimension(nfbrps), intent(in)              :: lstfbr
482double precision, dimension(nfbrps), intent(out)    :: effnrm
483
484! Local variables
485
486integer                        :: ifac  , iloc
487double precision               :: srfbn
488double precision, dimension(3) :: srfnor
489double precision, dimension(:,:), pointer :: forbr
490
491!===============================================================================
492
493call field_get_val_v(iforbr, forbr)
494
495do iloc = 1, nfbrps
496  ifac = lstfbr(iloc)
497  srfbn = surfbn(ifac)
498  srfnor(1) = surfbo(1,ifac) / srfbn
499  srfnor(2) = surfbo(2,ifac) / srfbn
500  srfnor(3) = surfbo(3,ifac) / srfbn
501  effnrm(iloc) =  (  forbr(1,ifac)*srfnor(1)                                 &
502                   + forbr(2,ifac)*srfnor(2)                                 &
503                   + forbr(3,ifac)*srfnor(3)) / srfbn
504enddo
505
506!--------
507! Formats
508!--------
509
510!----
511! End
512!----
513
514return
515end subroutine post_stress_normal
516
517!===============================================================================
518! Function:
519! ---------
520
521!> \brief Compute tangential stress at boundary.
522
523!-------------------------------------------------------------------------------
524! Arguments
525!______________________________________________________________________________.
526!  mode           name          role                                           !
527!______________________________________________________________________________!
528!> \param[in]     nfbrps        number of boundary faces to postprocess
529!> \param[in]     lstfbr        list of boundary faces to postprocess
530!> \param[out]    stress        stress at selected faces
531!_______________________________________________________________________________
532
533subroutine post_stress_tangential &
534 ( nfbrps , lstfbr ,                                              &
535   stress )
536
537!===============================================================================
538
539!===============================================================================
540! Module files
541!===============================================================================
542
543use paramx
544use pointe
545use entsor
546use cstnum
547use optcal
548use numvar
549use parall
550use period
551use mesh
552use field
553
554!===============================================================================
555
556implicit none
557
558! Arguments
559
560integer, intent(in)                                 :: nfbrps
561integer, dimension(nfbrps), intent(in)              :: lstfbr
562double precision, dimension(3, nfbrps), intent(out) :: stress
563
564! Local variables
565
566integer                        :: ifac  , iloc
567double precision               :: srfbn, fornor
568double precision, dimension(3) :: srfnor
569double precision, dimension(:,:), pointer :: forbr
570
571!===============================================================================
572
573call field_get_val_v(iforbr, forbr)
574
575do iloc = 1, nfbrps
576  ifac = lstfbr(iloc)
577  srfbn = surfbn(ifac)
578  srfnor(1) = surfbo(1,ifac) / srfbn
579  srfnor(2) = surfbo(2,ifac) / srfbn
580  srfnor(3) = surfbo(3,ifac) / srfbn
581  fornor =    forbr(1,ifac)*srfnor(1)                                 &
582            + forbr(2,ifac)*srfnor(2)                                 &
583            + forbr(3,ifac)*srfnor(3)
584  stress(1,iloc) = (forbr(1,ifac) - fornor*srfnor(1)) / srfbn
585  stress(2,iloc) = (forbr(2,ifac) - fornor*srfnor(2)) / srfbn
586  stress(3,iloc) = (forbr(3,ifac) - fornor*srfnor(3)) / srfbn
587enddo
588
589!--------
590! Formats
591!--------
592
593!----
594! End
595!----
596
597return
598end subroutine post_stress_tangential
599