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