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!=============================================================================== 24! Function : 25! -------- 26 27!> \file condli.f90 28!> 29!> \brief Translation of the boundary conditions given by 30!> cs_user_boundary_conditions 31!> in a form that fits to the solver. 32!> 33!> The values at a boundary face \f$ \fib \f$ stored in the face center 34!> \f$ \centf \f$ of the variable \f$ P \f$ and its diffusive flux \f$ Q \f$ 35!> are written as: 36!> \f[ 37!> P_\centf = A_P^g + B_P^g P_\centi 38!> \f] 39!> and 40!> \f[ 41!> Q_\centf = -\left(A_P^f + B_P^f P_\centi\right) 42!> \f] 43!> where \f$ P_\centi \f$ is the value of the variable \f$ P \f$ at the 44!> neighboring cell. 45!> 46!> Warning: 47!> - if we consider an increment of a variable, the boundary conditions 48!> read: 49!> \f[ 50!> \delta P_\centf = B_P^g \delta P_\centi 51!> \f] 52!> and 53!> \f[ 54!> \delta Q_\centf = -B_P^f \delta P_\centi 55!> \f] 56!> 57!> - for a vector field such as the veclocity \f$ \vect{u} \f$ the boundary 58!> conditions may read: 59!> \f[ 60!> \vect{u}_\centf = \vect{A}_u^g + \tens{B}_u^g \vect{u}_\centi 61!> \f] 62!> and 63!> \f[ 64!> \vect{Q}_\centf = -\left(\vect{A}_u^f + \tens{B}_u^f \vect{u}_\centi\right) 65!> \f] 66!> where \f$ \tens{B}_u^g \f$ and \f$ \tens{B}_u^f \f$ are 3x3 tensor matrix 67!> which coupled veclocity components next to a boundary. 68!> 69!> Please refer to the 70!> <a href="../../theory.pdf#boundary"><b>boundary conditions</b></a> section 71!> of the theory guide for more informations, as well as the 72!> <a href="../../theory.pdf#condli"><b>condli</b></a> section. 73!------------------------------------------------------------------------------- 74 75!------------------------------------------------------------------------------- 76! Arguments 77!______________________________________________________________________________. 78! mode name role ! 79!______________________________________________________________________________! 80!> \param[in] nvar total number of variables 81!> \param[in] nscal total number of scalars 82!> \param[in] iterns iteration number on Navier-Stokes equations 83!> \param[in] isvhb indicator to save exchange coeffient 84!> at the walls 85!> \param[in] itrale ALE iteration number 86!> \param[in] itrale ALE iteration number 87!> \param[in] italim for ALE 88!> \param[in] itrfin for ALE 89!> \param[in] ineefl for ALE 90!> \param[in] itrfup for ALE 91!> \param[in] cofale work array for FSI 92!> \param[in] xprale work array for FSI 93!> \param[in,out] icodcl face boundary condition code: 94!> - 1 Dirichlet 95!> - 2 Radiative outlet 96!> - 3 Neumann 97!> - 4 sliding and 98!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 99!> - 5 smooth wall and 100!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 101!> - 6 rough wall and 102!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 103!> - 9 free inlet/outlet 104!> (input mass flux blocked to 0) 105!> - 10 Boundary value related to the next cell 106!> value by an affine function 107!> - 11 Generalized Dirichlet for vectors 108!> - 12 Dirichlet boundary value related to the 109!> next cell value by an affine function for 110!> the advection operator and Neumann for the 111!> diffusion operator 112!> - 13 Dirichlet for the advection operator and 113!> Neumann for the diffusion operator 114!> - 14 Generalized symmetry for vectors (used for 115!> Marangoni effects modeling) 116!> - 15 Neumann for the advection operator and 117!> homogeneous Neumann for the diffusion 118!> operator (walls with hydro. pressure for 119!> the compressible module) 120!> \param[in,out] isostd indicator for standard outlet 121!> and reference face index 122!> \param[in] dt time step (per cell) 123!> \param[in,out] rcodcl boundary condition values: 124!> - rcodcl(1) value of the dirichlet 125!> - rcodcl(2) value of the exterior exchange 126!> coefficient (infinite if no exchange) 127!> - rcodcl(3) value flux density 128!> (negative if gain) in w/m2 or roughness 129!> in m if icodcl=6 130!> -# for the velocity \f$ (\mu+\mu_T) 131!> \gradv \vect{u} \cdot \vect{n} \f$ 132!> -# for the pressure \f$ \Delta t 133!> \grad P \cdot \vect{n} \f$ 134!> -# for a scalar \f$ cp \left( K + 135!> \dfrac{K_T}{\sigma_T} \right) 136!> \grad T \cdot \vect{n} \f$ 137!> \param[out] visvdr viscosite dynamique ds les cellules 138!> de bord apres amortisst de v driest 139!> \param[out] hbord coefficients d'echange aux bords 140!> \param[out] theipb boundary temperature in \f$ \centip \f$ 141!> (more exaclty the energetic variable) 142!_______________________________________________________________________________ 143 144subroutine condli & 145 ( nvar , nscal , iterns , & 146 isvhb , & 147 itrale , italim , itrfin , ineefl , itrfup , & 148 cofale , xprale , & 149 icodcl , isostd , & 150 dt , rcodcl , & 151 visvdr , hbord , theipb ) 152 153!=============================================================================== 154! Module files 155!=============================================================================== 156 157use paramx 158use numvar 159use optcal 160use alaste 161use alstru 162use atincl, only: iautom, iprofm 163use coincl, only: fment, ientfu, ientgb, ientgf, ientox, tkent, qimp 164use ppcpfu, only: inmoxy 165use cstphy 166use cstnum 167use pointe 168use entsor 169use albase 170use parall 171use ppppar 172use ppthch 173use ppincl 174use cpincl 175use radiat 176use cplsat 177use mesh 178use field 179use field_operator 180use radiat 181use turbomachinery 182use darcy_module 183use cs_c_bindings 184 185!=============================================================================== 186 187implicit none 188 189! Arguments 190 191integer nvar , nscal , iterns, isvhb 192integer itrale , italim , itrfin , ineefl , itrfup, nbt2h 193 194double precision, pointer, dimension(:) :: xprale 195double precision, pointer, dimension(:,:) :: cofale 196integer, pointer, dimension(:,:) :: icodcl 197integer, dimension(nfabor+1) :: isostd 198double precision, pointer, dimension(:) :: dt 199double precision, pointer, dimension(:,:,:) :: rcodcl 200double precision, pointer, dimension(:) :: visvdr, hbord, theipb 201 202! Local variables 203 204integer ifac , iel , ivar 205integer isou , jsou , ii 206integer ihcp , iscal 207integer inc , iprev , iccocg 208integer isoent, isorti, ncpt, isocpt(2) 209integer iclsym, ipatur, ipatrg, isvhbl 210integer iscacp, ifcvsl 211integer itplus, itstar 212integer f_id, iut, ivt, iwt, ialt, iflmab 213integer kbfid, b_f_id 214integer keyvar 215integer dimrij, f_dim 216integer kturt, turb_flux_model, turb_flux_model_type 217 218double precision sigma , cpp , rkl , visls_0 219double precision hint , hext , pimp , dimp, cfl 220double precision pinf , ratio 221double precision hintt(6) 222double precision flumbf, visclc, visctc, distbf 223double precision srfbnf, normal(3) 224double precision rinfiv(3), pimpv(3), qimpv(3), hextv(3), cflv(3) 225double precision b_pvari(3) 226double precision visci(3,3), fikis, viscis, distfi 227double precision temp, exchange_coef 228double precision turb_schmidt 229double precision, allocatable, dimension(:) :: pimpts, hextts, qimpts, cflts 230double precision sigmae 231 232character(len=80) :: fname 233 234double precision, dimension(:,:), pointer :: disale 235double precision, dimension(:,:), pointer :: xyzno0 236double precision, allocatable, dimension(:,:) :: velipb 237double precision, pointer, dimension(:,:) :: rijipb 238double precision, allocatable, dimension(:,:) :: grad 239double precision, allocatable, dimension(:,:,:) :: gradv 240double precision, allocatable, dimension(:,:,:) :: gradts 241double precision, pointer, dimension(:,:) :: dttens, visten 242double precision, dimension(:), pointer :: tplusp, tstarp, yplbr 243double precision, pointer, dimension(:,:) :: forbr 244double precision, dimension(:,:), pointer :: coefaut, cofafut, cofarut 245double precision, dimension(:,:,:), pointer :: coefbut, cofbfut, cofbrut 246double precision, dimension(:), pointer :: bmasfl 247double precision, dimension(:), pointer :: bfconv, bhconv 248double precision, dimension(:,:), pointer :: coefau, cofafu, cfaale, claale 249double precision, dimension(:,:,:), pointer :: coefbu, cofbfu, cfbale, clbale 250double precision, dimension(:), pointer :: coefap, coefbp, cofafp, cofbfp 251double precision, dimension(:,:), pointer :: coefav, cofafv 252double precision, dimension(:,:,:), pointer :: coefbv, cofbfv 253double precision, dimension(:,:), pointer :: coefats, cofafts, cofadts 254double precision, dimension(:,:,:), pointer :: cofbdts, coefbts, cofbfts 255double precision, dimension(:,:), pointer :: vel, vela 256double precision, dimension(:), pointer :: crom 257 258double precision, dimension(:), pointer :: viscl, visct, viscls 259double precision, dimension(:), pointer :: cpro_cp, cpro_cv, cvar_s, cvara_s 260double precision, dimension(:,:), pointer :: cvar_v, cvara_v 261double precision, dimension(:), pointer :: bvar_s, btemp_s 262double precision, dimension(:,:), pointer :: bvar_v 263double precision, dimension(:), pointer :: cpro_visma_s 264double precision, dimension(:,:), pointer :: cvar_ts, cvara_ts, cpro_visma_v 265 266integer, allocatable, dimension(:) :: lbt2h 267double precision, allocatable, dimension(:) :: vbt2h 268 269! darcy arrays 270double precision, dimension(:), pointer :: permeability 271double precision, dimension(:,:), pointer :: tensor_permeability 272 273type(var_cal_opt) :: vcopt 274 275!=============================================================================== 276! Interfaces 277!=============================================================================== 278 279!=============================================================================== 280! Interfaces 281!=============================================================================== 282 283interface 284 285 subroutine clptur(nscal, isvhb, icodcl, rcodcl, velipb, rijipb, & 286 visvdr, hbord, theipb) 287 288 implicit none 289 integer :: nscal, isvhb 290 integer, pointer, dimension(:,:) :: icodcl 291 double precision, pointer, dimension(:,:,:) :: rcodcl 292 double precision, dimension(:,:) :: velipb 293 double precision, pointer, dimension(:,:) :: rijipb 294 double precision, pointer, dimension(:) :: visvdr, hbord, theipb 295 296 end subroutine clptur 297 298 subroutine clptrg(nscal, isvhb, icodcl, rcodcl, velipb, rijipb, & 299 visvdr, hbord, theipb) 300 301 implicit none 302 integer :: nscal, isvhb 303 integer, pointer, dimension(:,:) :: icodcl 304 double precision, pointer, dimension(:,:,:) :: rcodcl 305 double precision, dimension(:,:) :: velipb 306 double precision, pointer, dimension(:,:) :: rijipb 307 double precision, pointer, dimension(:) :: visvdr, hbord, theipb 308 309 end subroutine clptrg 310 311 subroutine clsyvt(nscal, icodcl, rcodcl, velipb, rijipb) 312 313 implicit none 314 integer :: nscal 315 integer, pointer, dimension(:,:) :: icodcl 316 double precision, pointer, dimension(:,:,:) :: rcodcl 317 double precision, dimension(:,:) :: velipb 318 double precision, pointer, dimension(:,:) :: rijipb 319 320 end subroutine clsyvt 321 322 subroutine cs_boundary_conditions_complete(nvar, itypfb, icodcl, rcodcl) & 323 bind(C, name='cs_boundary_conditions_complete') 324 325 use, intrinsic :: iso_c_binding 326 implicit none 327 integer(c_int), value :: nvar 328 integer(c_int), dimension(*), intent(inout) :: itypfb, icodcl 329 real(kind=c_double), dimension(*), intent(inout) :: rcodcl 330 331 end subroutine cs_boundary_conditions_complete 332 333 subroutine cs_syr_coupling_recv_boundary(nvar, bc_type, icodcl, rcodcl) & 334 bind(C, name = 'cs_syr_coupling_recv_boundary') 335 336 use, intrinsic :: iso_c_binding 337 implicit none 338 integer(kind=c_int), value :: nvar 339 integer(kind=c_int), dimension(*), intent(inout) :: bc_type 340 integer(kind=c_int), dimension(*), intent(inout) :: icodcl 341 real(kind=c_double), dimension(*), intent(inout) :: rcodcl 342 343 end subroutine cs_syr_coupling_recv_boundary 344 345 subroutine cs_syr_coupling_exchange_volume() & 346 bind(C, name = 'cs_syr_coupling_exchange_volume') 347 348 use, intrinsic :: iso_c_binding 349 implicit none 350 351 end subroutine cs_syr_coupling_exchange_volume 352 353 subroutine strpre(itrale, italim, ineefl, impale, xprale, cofale) 354 355 implicit none 356 integer :: itrale, italim, ineefl 357 integer, dimension(:) :: impale 358 double precision, pointer, dimension(:) :: xprale 359 double precision, pointer, dimension(:,:) :: cofale 360 361 end subroutine strpre 362 363 end interface 364 365!=============================================================================== 366! 0. User calls 367!=============================================================================== 368 369rijipb => null() 370 371call precli(nvar, icodcl, rcodcl) 372 373! - Interface Code_Saturne 374! ====================== 375 376! N.B. Zones de face de bord : on utilise provisoirement les zones des 377! physiques particulieres, meme sans physique particuliere 378! -> sera modifie lors de la restructuration des zones de bord 379 380call uiclim & 381 ( nozppm, & 382 iqimp, icalke, ientat, ientcp, inmoxy, ientox, & 383 ientfu, ientgb, ientgf, iprofm, iautom, & 384 itypfb, izfppp, icodcl, & 385 qimp, qimpat, qimpcp, dh, xintur, & 386 timpat, timpcp, tkent , fment, distch, nvar, rcodcl) 387 388if (ippmod(iphpar).eq.0.or.ippmod(igmix).ge.0.or.ippmod(icompf).ge.0) then 389 390 ! ON NE FAIT PAS DE LA PHYSIQUE PARTICULIERE 391 392 call stdtcl & 393 ( nbzppm , nozppm , & 394 iqimp , icalke , qimp , dh , xintur, & 395 itypfb , izfppp , & 396 rcodcl ) 397 398endif 399 400call cs_boundary_conditions_complete(nvar, itypfb, icodcl, rcodcl) 401 402! User-defined functions 403! ========================== 404 405call cs_f_user_boundary_conditions & 406 (nvar, nscal, icodcl, itrifb, itypfb, izfppp, dt, rcodcl ) 407 408call user_boundary_conditions(nvar, itypfb, icodcl, rcodcl) 409 410! Check consistency with GUI definitions 411 412call uiclve(nozppm) 413 414! BC'based coupling with other code_saturne instances. 415 416if (nbrcpl.gt.0) then 417 call cscfbr(nscal, icodcl, itypfb, dt, rcodcl) 418endif 419 420! Synthetic Eddy Method for L.E.S. 421 422call synthe(ttcabs, dt, rcodcl) 423 424! ALE method (mesh velocity BC and vertices displacement) 425 426if (iale.ge.1) then 427 428 call field_get_val_v(fdiale, disale) 429 call field_get_val_v_by_name("vtx_coord0", xyzno0) 430 431 do ii = 1, nnod 432 impale(ii) = 0 433 enddo 434 435 ! - Interface Code_Saturne 436 ! ====================== 437 438 call uialcl(ibfixe, igliss, ivimpo, ifresf, & 439 ialtyb, impale, disale, & 440 iuma, ivma, iwma, & 441 rcodcl) 442 443 ! TODO in the future version: remove dt, xyzno0, and disale 444 ! because they are avaliable as fields. 445 446 call usalcl(itrale, nvar, nscal, & 447 icodcl, itypfb, ialtyb, impale, & 448 dt, rcodcl, xyzno0, disale) 449 450 ! Au cas ou l'utilisateur aurait touche disale sans mettre impale=1, on 451 ! remet le deplacement initial 452 do ii = 1, nnod 453 if (impale(ii).eq.0) then 454 disale(1,ii) = xyznod(1,ii)-xyzno0(1,ii) 455 disale(2,ii) = xyznod(2,ii)-xyzno0(2,ii) 456 disale(3,ii) = xyznod(3,ii)-xyzno0(3,ii) 457 endif 458 enddo 459 460 ! En cas de couplage de structures, on calcule un deplacement predit 461 if (nbstru.gt.0.or.nbaste.gt.0) then 462 call strpre(itrale, italim, ineefl, impale, xprale, cofale) 463 endif 464 465endif 466 467! UNE FOIS CERTAINS CODES DE CONDITIONS LIMITES INITIALISES PAR 468! L'UTILISATEUR, ON PEUT COMPLETER CES CODES PAR LES COUPLAGES 469! AUX BORDS (TYPE SYRTHES), SAUF SI ON DOIT Y REPASSER ENSUITE 470! POUR CENTRALISER CE QUI EST RELATIF AU COUPLAGE AVEC SYRTHES 471! ON POSITIONNE ICI L'APPEL AU COUPLAGE VOLUMIQUE SYRTHES 472! UTILE POUR BENIFICER DE LA DERNIERE VITESSE CALCULEE SI ON 473! BOUCLE SUR U/P. 474! LE COUPLAGE VOLUMIQUE DOIT ETRE APPELE AVANT LE SURFACIQUE 475! POUR RESPECTER LE SCHEMA DE COMMUNICATION 476 477if (itrfin.eq.1 .and. itrfup.eq.1) then 478 479 call cs_syr_coupling_exchange_volume 480 481 call cs_syr_coupling_recv_boundary(nvar, itypfb, icodcl, rcodcl) 482 483 if (nfpt1t.gt.0) then 484 call cou1di(nfabor, iscalt, icodcl, rcodcl) 485 endif 486 487 ! coupling 1D thermal model with condensation modelling 488 ! to take into account the solid temperature evolution over time 489 if (nftcdt.gt.0) then 490 call cs_tagmri(nfabor, iscalt, icodcl, rcodcl) 491 endif 492 493endif 494 495! Radiative transfer: add contribution to energy BCs. 496if (iirayo.gt.0 .and. itrfin.eq.1 .and. itrfup.eq.1) then 497 call cs_rad_transfer_bcs(nvar, itypfb, icodcl, dt, rcodcl) 498endif 499 500! For internal coupling, set itypfb to wall function by default 501! if not set by the user 502call cs_internal_coupling_bcs(itypfb) 503 504! Convert temperature to enthalpy for Dirichlet conditions 505 506if (itherm.eq.2) then 507 508 nbt2h = 0 509 allocate(lbt2h(nfabor)) 510 allocate(vbt2h(nfabor)) 511 512 ivar = isca(iscalt) 513 514 ! Filter Dirichlet/imposed value faces 515 516 do ii = 1, nfabor 517 if (icodcl(ii,ivar).lt.0) then 518 nbt2h = nbt2h + 1 519 lbt2h(nbt2h) = ii - 1 ! 0-based numbering 520 icodcl(ii,ivar) = -icodcl(ii,ivar) 521 vbt2h(ii) = rcodcl(ii,ivar,1) 522 else 523 vbt2h(ii) = 0.d0 524 endif 525 enddo 526 527 call cs_ht_convert_t_to_h_faces_l(nbt2h, lbt2h, vbt2h, rcodcl(:,ivar,1)) 528 529endif 530 531!=============================================================================== 532! 1. initializations 533!=============================================================================== 534 535if (ippmod(idarcy).eq.1) then 536 if (darcy_anisotropic_permeability.eq.0) then 537 call field_get_val_s_by_name('permeability', permeability) 538 else 539 call field_get_id('permeability', f_id) 540 call field_get_val_v(f_id, tensor_permeability) 541 endif 542endif 543 544call field_get_key_id("variable_id", keyvar) 545 546! allocate temporary arrays 547allocate(velipb(nfabor,3)) 548if (irij.ge.1) then 549 call field_get_dim(ivarfl(irij), dimrij) ! dimension of Rij 550 allocate(pimpts(dimrij)) 551 allocate(hextts(dimrij)) 552 allocate(qimpts(dimrij)) 553 allocate(cflts(dimrij)) 554 do isou = 1 , dimrij 555 pimpts(isou) = 0 556 hextts(isou) = 0 557 qimpts(isou) = 0 558 cflts(isou) = 0 559 enddo 560endif 561! coefa and coefb are required to compute the cell gradients for the wall 562! turbulent boundary conditions. 563! so, their initial values are kept (note that at the first time step, they are 564! initialized to zero flux in inivar.f90) 565 566! velipb stores the velocity in i' of boundary cells 567 568! initialize variables to avoid compiler warnings 569 570rinfiv(1) = rinfin 571rinfiv(2) = rinfin 572rinfiv(3) = rinfin 573 574! pointers to y+, t+ and t* if saved 575 576yplbr => null() 577tplusp => null() 578tstarp => null() 579 580! initialization of the array storing yplus 581! which is computed in clptur.f90 and/or clptrg.f90 582call field_get_id_try('yplus', f_id) 583if (f_id.ge.0) then 584 call field_get_val_s(f_id, yplbr) 585 do ifac = 1, nfabor 586 yplbr(ifac) = 0.d0 587 enddo 588endif 589 590call field_get_id_try('tplus', itplus) 591if (itplus.ge.0) then 592 call field_get_val_s (itplus, tplusp) 593 do ifac = 1, nfabor 594 tplusp(ifac) = 0.d0 595 enddo 596endif 597 598call field_get_id_try('tstar', itstar) 599if (itstar.ge.0) then 600 call field_get_val_s (itstar, tstarp) 601 do ifac = 1, nfabor 602 tstarp(ifac) = 0.d0 603 enddo 604endif 605 606! map field arrays 607call field_get_val_v(ivarfl(iu), vel) 608call field_get_val_prev_v(ivarfl(iu), vela) 609 610! pointers to the mass fluxes 611call field_get_key_int(ivarfl(iu), kbmasf, iflmab) 612call field_get_val_s(iflmab, bmasfl) 613 614! pointers to specific fields 615 616if (iirayo.ge.1) then 617 call field_get_val_s_by_name("rad_convective_flux", bfconv) 618 call field_get_val_s_by_name("rad_exchange_coefficient", bhconv) 619endif 620 621if (idtten.ge.0) call field_get_val_v(idtten, dttens) 622 623if (iforbr.ge.0 .and. iterns.eq.1) call field_get_val_v(iforbr, forbr) 624 625! pointers to velocity bc coefficients 626call field_get_coefa_v(ivarfl(iu), coefau) 627call field_get_coefb_v(ivarfl(iu), coefbu) 628call field_get_coefaf_v(ivarfl(iu), cofafu) 629call field_get_coefbf_v(ivarfl(iu), cofbfu) 630 631! pointers to boundary variable values 632 633call field_get_key_id("boundary_value_id", kbfid) 634 635!=============================================================================== 636! 2. treatment of types of bcs given by itypfb 637!=============================================================================== 638 639if ( ippmod(iphpar).ge.1.and.ippmod(igmix).eq.-1 & 640 .and.ippmod(ieljou).eq.-1.and.ippmod(ielarc).eq.-1 & 641 .or.ippmod(icompf).ge.0.and.ippmod(igmix).ge.0) then 642 call pptycl & 643 ( nvar , .false., & 644 icodcl , itypfb , izfppp , & 645 dt , & 646 rcodcl ) 647endif 648 649if (iale.ge.1) then 650 call altycl & 651 ( itypfb , ialtyb , icodcl , impale , .false. , & 652 dt , & 653 rcodcl ) 654endif 655 656if (iturbo.ne.0) then 657 call mmtycl(itypfb, rcodcl) 658endif 659 660call typecl & 661 ( nvar , nscal , .false. , & 662 itypfb , itrifb , icodcl , isostd , & 663 rcodcl ) 664 665!=============================================================================== 666! 3. check the consistency of the bcs 667!=============================================================================== 668 669call vericl(nvar, nscal, itypfb, icodcl) 670 671!=============================================================================== 672! 4. variables 673!=============================================================================== 674 675! --- physical quantities 676call field_get_val_s(iviscl, viscl) 677call field_get_val_s(ivisct, visct) 678 679!=============================================================================== 680! 5. compute the temperature or the enthalpy in i' for boundary cells 681! (thanks to the formula: fi + grad(fi).ii') 682 683! for the coupling with syrthes 684! theipb is used by cs_syr_coupling_send_boundary after condli 685! for the coupling with the 1d wall thermal module 686! theipb is used by cou1do after condli 687! for the radiation module 688! theipb is used to compute the required flux in raypar 689 690! ceci pourrait en pratique etre hors de la boucle. 691 692!=============================================================================== 693 694! allocate a temporary array for the gradient reconstruction 695allocate(grad(3,ncelet)) 696 697! pour le couplage syrthes ou module thermique 1d 698! ----------------------------------------------- 699! ici, on fait une boucle "inutile" (on ne fait quelque chose 700! que pour icpsyr = 1). c'est pour preparer le traitement 701! eventuel de plusieurs temperatures (ie plusieurs couplages 702! syrthes a la fois ; noter cependant que meme dans ce cas, 703! une seule temperature sera recue de chaque couplage. en polyph, 704! il faudrait ensuite reconstruire les enthalpies ... 705! plus tard si necessaire). 706! ici, il ne peut y avoir qu'un seul scalaire avec icpsyr = 1 et 707! ce uniquement s'il y a effectivement couplage avec syrthes 708! (sinon, on s'est arrete dans verini) 709! dans le cas du couplage avec le module 1d, on utilise iscalt. 710 711! pour le rayonnement 712! ------------------- 713! on calcule la valeur en i' s'il y a une variable 714! thermique 715 716 717! on recherche l'unique scalaire qui convient 718! (ce peut etre t, h, ou e (en compressible)) 719 720! compute the boundary value of required scalars 721 722! Check for boundary values 723 724do ii = 1, nscal 725 726 ivar = isca(ii) 727 f_id = ivarfl(ivar) 728 729 call field_get_key_int(f_id, kbfid, b_f_id) 730 call field_get_dim(f_id, f_dim) 731 732 ! if thermal variable has no boundary but temperature does, use it 733 if (b_f_id .lt. 0 .and. ii.eq.iscalt .and. itherm.eq.2) then 734 b_f_id = itempb 735 endif 736 737 if (b_f_id .ge. 0) then 738 if (f_dim.eq.1) then 739 call field_get_val_s(b_f_id, bvar_s) 740 else 741 call field_get_val_v(b_f_id, bvar_v) 742 endif 743 else if (ii.eq.iscalt) then 744 bvar_s => null() ! no boundary field, but may need theipb 745 else 746 cycle ! nothing to do for this scalar 747 endif 748 749 call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 750 751 if (f_dim.eq.1) then 752 if (itbrrb.eq.1 .and. vcopt%ircflu.eq.1) then 753 754 call field_get_val_s(ivarfl(ivar), cvar_s) 755 756 inc = 1 757 iprev = 1 758 iccocg = 1 759 760 call field_gradient_scalar(ivarfl(ivar), iprev, 0, inc, & 761 iccocg, & 762 grad) 763 764 call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 765 766 if (b_f_id .ge. 0) then 767 do ifac = 1 , nfabor 768 iel = ifabor(ifac) 769 bvar_s(ifac) = cvar_s(iel) & 770 + vcopt%ircflu & 771 * ( & 772 + grad(1,iel)*diipb(1,ifac) & 773 + grad(2,iel)*diipb(2,ifac) & 774 + grad(3,iel)*diipb(3,ifac) & 775 ) 776 enddo 777 else 778 do ifac = 1 , nfabor 779 iel = ifabor(ifac) 780 theipb(ifac) = cvar_s(iel) & 781 + vcopt%ircflu & 782 * ( & 783 + grad(1,iel)*diipb(1,ifac) & 784 + grad(2,iel)*diipb(2,ifac) & 785 + grad(3,iel)*diipb(3,ifac) & 786 ) 787 enddo 788 endif 789 790 else 791 792 call field_get_val_prev_s(ivarfl(ivar), cvara_s) 793 794 if (b_f_id .ge. 0) then 795 do ifac = 1 , nfabor 796 iel = ifabor(ifac) 797 bvar_s(ifac) = cvara_s(iel) 798 enddo 799 else 800 do ifac = 1 , nfabor 801 iel = ifabor(ifac) 802 theipb(ifac) = cvara_s(iel) 803 enddo 804 endif 805 806 endif 807 808 ! Copy bvar_s to theipb if both theipb and bvar_s present 809 810 if (b_f_id .ge. 0 .and. ii.eq.iscalt) then 811 do ifac = 1 , nfabor 812 theipb(ifac) = bvar_s(ifac) 813 enddo 814 endif 815 816 elseif (b_f_id.ge.0) then 817 if (itbrrb.eq.1 .and. vcopt%ircflu.eq.1) then 818 call field_get_val_v(ivarfl(ivar), cvar_v) 819 820 allocate(gradv(3,3,ncelet)) 821 822 inc = 1 823 iprev = 1 824 call field_gradient_vector(ivarfl(ivar), iprev, 0, inc, gradv) 825 826 do ifac = 1 , nfabor 827 iel = ifabor(ifac) 828 do isou = 1, 3 829 bvar_v(isou,ifac) = cvar_v(isou,iel) & 830 + gradv(1,isou,iel)*diipb(1,ifac) & 831 + gradv(2,isou,iel)*diipb(2,ifac) & 832 + gradv(3,isou,iel)*diipb(3,ifac) 833 enddo 834 enddo 835 836 deallocate(gradv) 837 838 else 839 call field_get_val_prev_v(ivarfl(ivar), cvara_v) 840 841 do ifac = 1 , nfabor 842 iel = ifabor(ifac) 843 do isou = 1, 3 844 bvar_v(isou,ifac) = cvara_v(isou,iel) 845 enddo 846 enddo 847 endif 848 endif 849 850enddo !nscal 851 852!=============================================================================== 853! 6. compute the velocity and Reynolds stesses tensor in i' for boundary cells 854! (thanks to the formula: fi + grad(fi).ii') if there are symmetry or 855! wall with wall functions boundary conditions 856!=============================================================================== 857 858! ---> indicator for symmetries or wall with wall functions 859 860iclsym = 0 861ipatur = 0 862ipatrg = 0 863do ifac = 1, nfabor 864 if (icodcl(ifac,iu).eq.4) then 865 iclsym = 1 866 elseif (icodcl(ifac,iu).eq.5) then 867 ipatur = 1 868 elseif (icodcl(ifac,iu).eq.6) then 869 ipatrg = 1 870 endif 871 if (iclsym.ne.0.and.ipatur.ne.0.and.ipatrg.ne.0) goto 100 872enddo 873 874100 continue 875 876if (irangp.ge.0) then 877 call parcmx(iclsym) 878 call parcmx(ipatur) 879 call parcmx(ipatrg) 880endif 881 882! ---> compute the velocity in i' for boundary cells 883 884if (iclsym.ne.0.or.ipatur.ne.0.or.ipatrg.ne.0.or.iforbr.ge.0) then 885 886 if (ntcabs.gt.1) then 887 888 ! allocate a temporary array 889 allocate(gradv(3,3,ncelet)) 890 891 inc = 1 892 iprev = 1 893 894 call field_gradient_vector(ivarfl(iu), iprev, 0, inc, gradv) 895 896 call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt) 897 898 do isou = 1, 3 899 do ifac = 1, nfabor 900 iel = ifabor(ifac) 901 velipb(ifac,isou) = vela(isou,iel) & 902 + vcopt%ircflu & 903 * ( & 904 + gradv(1,isou,iel)*diipb(1,ifac) & 905 + gradv(2,isou,iel)*diipb(2,ifac) & 906 + gradv(3,isou,iel)*diipb(3,ifac) & 907 ) 908 enddo 909 enddo 910 911 deallocate(gradv) 912 913 ! nb: at the first time step, coefa and coefb are unknown, so the walue 914 ! in i is stored instead of the value in i' 915 else 916 917 do isou = 1, 3 918 919 do ifac = 1, nfabor 920 iel = ifabor(ifac) 921 velipb(ifac,isou) = vela(isou,iel) 922 enddo 923 924 enddo 925 926 endif 927 928endif 929 930! ---> compute rij in i' for boundary cells 931 932if ((iclsym.ne.0.or.ipatur.ne.0.or.ipatrg.ne.0).and.itytur.eq.3) then 933 934 ! allocate a work array to store rij values at boundary faces 935 allocate(rijipb(nfabor,6)) 936 937 if (ntcabs.gt.1.and.irijrb.eq.1) then 938 939 call field_get_val_v(ivarfl(irij), cvar_ts) 940 941 inc = 1 942 iprev = 1 943 944 call field_get_key_struct_var_cal_opt(ivarfl(irij), vcopt) 945 946 ! allocate a temporary array 947 allocate(gradts(6,3,ncelet)) 948 949 call field_gradient_tensor(ivarfl(irij), iprev, 0, inc, gradts) 950 951 do ifac = 1 , nfabor 952 iel = ifabor(ifac) 953 do isou = 1, 6 954 rijipb(ifac,isou) = cvar_ts(isou,iel) & 955 + vcopt%ircflu & 956 * ( gradts(isou,1,iel)*diipb(1,ifac) & 957 + gradts(isou,2,iel)*diipb(2,ifac) & 958 + gradts(isou,3,iel)*diipb(3,ifac)) 959 enddo 960 enddo 961 962 ! nb: at the first time step, coefa and coefb are unknown, so the walue 963 ! in i is stored instead of the value in i' 964 else 965 966 call field_get_val_prev_v(ivarfl(irij), cvara_ts) 967 968 do ifac = 1 , nfabor 969 iel = ifabor(ifac) 970 do isou = 1, 6 971 rijipb(ifac,isou) = cvara_ts(isou, iel) 972 enddo 973 enddo 974 endif 975 976endif 977 978! free memory 979deallocate(grad) 980 981!=============================================================================== 982! 7. turbulence at walls: 983! (u,v,w,k,epsilon,rij,temperature) 984!=============================================================================== 985! --- on a besoin de velipb et de rijipb (et theipb pour le rayonnement) 986 987! on initialise visvdr a -999.d0. 988! dans clptur, on amortit la viscosite turbulente sur les cellules 989! de paroi si on a active van driest. la valeur finale est aussi 990! stockee dans visvdr. 991! plus loin, dans distyp, la viscosite sur les cellules 992! de paroi sera amortie une seconde fois. on se sert alors de 993! visvdr pour lui redonner une valeur correcte. 994if(itytur.eq.4.and.idries.eq.1) then 995 do iel=1,ncel 996 visvdr(iel) = -999.d0 997 enddo 998endif 999 1000if (ipatur.ne.0) then 1001 1002 ! smooth wall laws 1003 call clptur & 1004 ( nscal , isvhb , icodcl , & 1005 rcodcl , & 1006 velipb , rijipb , visvdr , & 1007 hbord , theipb ) 1008 1009endif 1010 1011if (ipatrg.ne.0) then 1012 1013 ! rough wall laws 1014 call clptrg & 1015 ( nscal , isvhb , icodcl , & 1016 rcodcl , & 1017 velipb , rijipb , visvdr , & 1018 hbord , theipb ) 1019 1020endif 1021 1022!=============================================================================== 1023! 8. symmetry for vectors and tensors 1024! (u,v,w,rij) 1025!=============================================================================== 1026! on a besoin de velipb et de rijipb 1027 1028do ifac = 1, nfabor 1029 isympa(ifac) = 1 1030enddo 1031 1032if (iclsym.ne.0) then 1033 1034 call clsyvt & 1035 ( nscal , icodcl , & 1036 rcodcl , & 1037 velipb , rijipb ) 1038 1039endif 1040 1041!=============================================================================== 1042! 9. velocity: outlet, Dirichlet and Neumann and convective outlet 1043!=============================================================================== 1044 1045! ---> outlet: in case of incoming mass flux, the mass flux is set to zero. 1046 1047isoent = 0 1048isorti = 0 1049do ifac = 1, nfabor 1050 1051 if (icodcl(ifac,iu).eq.9) then 1052 1053 flumbf = bmasfl(ifac) 1054 1055 ! --- physical properties 1056 iel = ifabor(ifac) 1057 visclc = viscl(iel) 1058 visctc = visct(iel) 1059 1060 ! --- geometrical quantities 1061 distbf = distb(ifac) 1062 1063 if (itytur.eq.3) then 1064 hint = visclc /distbf 1065 else 1066 hint = (visclc + visctc)/distbf 1067 endif 1068 1069 isorti = isorti + 1 1070 1071 if (flumbf.lt.-epzero) then 1072 1073 ! Dirichlet boundary condition 1074 !----------------------------- 1075 1076 ! coupled solving of the velocity components 1077 1078 pimpv(1) = 0.d0 1079 pimpv(2) = 0.d0 1080 pimpv(3) = 0.d0 1081 1082 call set_dirichlet_vector & 1083 ( coefau(:,ifac) , cofafu(:,ifac) , & 1084 coefbu(:,:,ifac), cofbfu(:,:,ifac), & 1085 pimpv , hint , rinfiv ) 1086 1087 1088 isoent = isoent + 1 1089 1090 else 1091 1092 ! Neumann boundary conditions 1093 !---------------------------- 1094 1095 dimp = 0.d0 1096 1097 ! coupled solving of the velocity components 1098 1099 qimpv(1) = 0.d0 1100 qimpv(2) = 0.d0 1101 qimpv(3) = 0.d0 1102 1103 call set_neumann_vector & 1104 ( coefau(:,ifac) , cofafu(:,ifac) , & 1105 coefbu(:,:,ifac), cofbfu(:,:,ifac), & 1106 qimpv , hint ) 1107 1108 endif 1109 1110 endif 1111 1112enddo 1113 1114call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt) 1115 1116if (mod(ntcabs,ntlist).eq.0 .or. vcopt%iwarni.ge. 0) then 1117 isocpt(1) = isoent 1118 isocpt(2) = isorti 1119 if (irangp.ge.0) then 1120 ncpt = 2 1121 call parism(ncpt, isocpt) 1122 endif 1123 if (isocpt(2).gt.0 .and. (vcopt%iwarni.ge.2.or.isocpt(1).gt.0)) then 1124 write(nfecra,3010) isocpt(1), isocpt(2) 1125 endif 1126endif 1127 1128! ---> dirichlet and neumann 1129 1130do ifac = 1, nfabor 1131 1132 iel = ifabor(ifac) 1133 1134 ! --- physical propreties 1135 visclc = viscl(iel) 1136 visctc = visct(iel) 1137 1138 ! --- geometrical quantities 1139 distbf = distb(ifac) 1140 1141 if (itytur.eq.3) then 1142 hint = visclc /distbf 1143 else 1144 hint = ( visclc+visctc )/distbf 1145 endif 1146 1147 ! dirichlet boundary conditions 1148 !------------------------------ 1149 1150 if (icodcl(ifac,iu).eq.1) then 1151 1152 1153 pimpv(1) = rcodcl(ifac,iu,1) 1154 pimpv(2) = rcodcl(ifac,iv,1) 1155 pimpv(3) = rcodcl(ifac,iw,1) 1156 hextv(1) = rcodcl(ifac,iu,2) 1157 hextv(2) = rcodcl(ifac,iv,2) 1158 hextv(3) = rcodcl(ifac,iw,2) 1159 1160 call set_dirichlet_vector & 1161 ( coefau(:,ifac) , cofafu(:,ifac) , & 1162 coefbu(:,:,ifac), cofbfu(:,:,ifac), & 1163 pimpv , hint , hextv ) 1164 1165 1166 ! neumann boundary conditions 1167 !---------------------------- 1168 1169 elseif (icodcl(ifac,iu).eq.3) then 1170 1171 ! coupled solving of the velocity components 1172 1173 qimpv(1) = rcodcl(ifac,iu,3) 1174 qimpv(2) = rcodcl(ifac,iv,3) 1175 qimpv(3) = rcodcl(ifac,iw,3) 1176 1177 call set_neumann_vector & 1178 ( coefau(:,ifac) , cofafu(:,ifac) , & 1179 coefbu(:,:,ifac), cofbfu(:,:,ifac), & 1180 qimpv , hint ) 1181 1182 ! convective boundary conditions 1183 !------------------------------- 1184 1185 elseif (icodcl(ifac,iu).eq.2 .and. iterns.le.1) then 1186 1187 ! coupled solving of the velocity components 1188 1189 pimpv(1) = rcodcl(ifac,iu,1) 1190 cflv(1) = rcodcl(ifac,iu,2) 1191 pimpv(2) = rcodcl(ifac,iv,1) 1192 cflv(2) = rcodcl(ifac,iv,2) 1193 pimpv(3) = rcodcl(ifac,iw,1) 1194 cflv(3) = rcodcl(ifac,iw,2) 1195 1196 call set_convective_outlet_vector & 1197 ( coefau(:,ifac) , cofafu(:,ifac) , & 1198 coefbu(:,:,ifac), cofbfu(:,:,ifac), & 1199 pimpv , cflv , hint ) 1200 1201 ! imposed value for the convection operator, imposed flux for diffusion 1202 !---------------------------------------------------------------------- 1203 1204 elseif (icodcl(ifac,iu).eq.13) then 1205 1206 pimpv(1) = rcodcl(ifac,iu,1) 1207 pimpv(2) = rcodcl(ifac,iv,1) 1208 pimpv(3) = rcodcl(ifac,iw,1) 1209 1210 qimpv(1) = rcodcl(ifac,iu,3) 1211 qimpv(2) = rcodcl(ifac,iv,3) 1212 qimpv(3) = rcodcl(ifac,iw,3) 1213 1214 call set_dirichlet_conv_neumann_diff_vector & 1215 ( coefau(:,ifac) , cofafu(:,ifac) , & 1216 coefbu(:,:,ifac), cofbfu(:,:,ifac), & 1217 pimpv , qimpv ) 1218 1219 ! convective boundary for Marangoni effects (generalized symmetry condition) 1220 !--------------------------------------------------------------------------- 1221 1222 elseif (icodcl(ifac,iu).eq.14) then 1223 1224 pimpv(1) = rcodcl(ifac,iu,1) 1225 pimpv(2) = rcodcl(ifac,iv,1) 1226 pimpv(3) = rcodcl(ifac,iw,1) 1227 1228 qimpv(1) = rcodcl(ifac,iu,3) 1229 qimpv(2) = rcodcl(ifac,iv,3) 1230 qimpv(3) = rcodcl(ifac,iw,3) 1231 1232 normal(1) = surfbo(1,ifac)/surfbn(ifac) 1233 normal(2) = surfbo(2,ifac)/surfbn(ifac) 1234 normal(3) = surfbo(3,ifac)/surfbn(ifac) 1235 1236 1237 ! coupled solving of the velocity components 1238 1239 call set_generalized_sym_vector & 1240 ( coefau(:,ifac) , cofafu(:,ifac) , & 1241 coefbu(:,:,ifac), cofbfu(:,:,ifac), & 1242 pimpv , qimpv , hint , normal ) 1243 1244 ! Neumann on the normal component, Dirichlet on tangential components 1245 !-------------------------------------------------------------------- 1246 1247 elseif (icodcl(ifac,iu).eq.11) then 1248 1249 ! Dirichlet to impose on the tangential components 1250 pimpv(1) = rcodcl(ifac,iu,1) 1251 pimpv(2) = rcodcl(ifac,iv,1) 1252 pimpv(3) = rcodcl(ifac,iw,1) 1253 1254 ! Flux to impose on the normal component 1255 qimpv(1) = rcodcl(ifac,iu,3) 1256 qimpv(2) = rcodcl(ifac,iv,3) 1257 qimpv(3) = rcodcl(ifac,iw,3) 1258 1259 normal(1) = surfbo(1,ifac)/surfbn(ifac) 1260 normal(2) = surfbo(2,ifac)/surfbn(ifac) 1261 normal(3) = surfbo(3,ifac)/surfbn(ifac) 1262 1263 1264 ! coupled solving of the velocity components 1265 1266 call set_generalized_dirichlet_vector & 1267 ( coefau(:,ifac) , cofafu(:,ifac) , & 1268 coefbu(:,:,ifac), cofbfu(:,:,ifac), & 1269 pimpv , qimpv , hint , normal ) 1270 1271 1272 endif 1273 1274enddo 1275 1276!=============================================================================== 1277! 10. pressure: Dirichlet and Neumann and convective outlet 1278!=============================================================================== 1279 1280call field_get_coefa_s(ivarfl(ipr), coefap) 1281call field_get_coefb_s(ivarfl(ipr), coefbp) 1282call field_get_coefaf_s(ivarfl(ipr), cofafp) 1283call field_get_coefbf_s(ivarfl(ipr), cofbfp) 1284 1285call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt) 1286 1287if (ivofmt.gt.0) call field_get_val_s(icrom, crom) ! FIXME consistency with correction step 1288 1289do ifac = 1, nfabor 1290 1291 iel = ifabor(ifac) 1292 1293 ! --- geometrical quantities 1294 distbf = distb(ifac) 1295 1296 ! if a flux dt.grad p (w/m2) is set in cs_user_boundary 1297 if (iand(vcopt%idften, ISOTROPIC_DIFFUSION).ne.0) then 1298 hint = dt(iel)/distbf 1299 if (ivofmt.gt.0) hint = hint/crom(iel) 1300 if (ippmod(idarcy).eq.1) hint = permeability(iel)/distbf 1301 else if (iand(vcopt%idften, ORTHOTROPIC_DIFFUSION).ne.0) then 1302 hint = ( dttens(1, iel)*surfbo(1,ifac)**2 & 1303 + dttens(2, iel)*surfbo(2,ifac)**2 & 1304 + dttens(3, iel)*surfbo(3,ifac)**2 & 1305 ) / (surfbn(ifac)**2 * distbf) 1306 if (ivofmt.gt.0) hint = hint/crom(iel) 1307 ! symmetric tensor diffusivity 1308 else if (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 1309 if (ippmod(idarcy).eq.-1) then 1310 visci(1,1) = dttens(1,iel) 1311 visci(2,2) = dttens(2,iel) 1312 visci(3,3) = dttens(3,iel) 1313 visci(1,2) = dttens(4,iel) 1314 visci(2,1) = dttens(4,iel) 1315 visci(2,3) = dttens(5,iel) 1316 visci(3,2) = dttens(5,iel) 1317 visci(1,3) = dttens(6,iel) 1318 visci(3,1) = dttens(6,iel) 1319 else 1320 visci(1,1) = tensor_permeability(1,iel) 1321 visci(2,2) = tensor_permeability(2,iel) 1322 visci(3,3) = tensor_permeability(3,iel) 1323 visci(1,2) = tensor_permeability(4,iel) 1324 visci(2,1) = tensor_permeability(4,iel) 1325 visci(2,3) = tensor_permeability(5,iel) 1326 visci(3,2) = tensor_permeability(5,iel) 1327 visci(1,3) = tensor_permeability(6,iel) 1328 visci(3,1) = tensor_permeability(6,iel) 1329 endif 1330 1331 ! ||ki.s||^2 1332 viscis = ( visci(1,1)*surfbo(1,ifac) & 1333 + visci(1,2)*surfbo(2,ifac) & 1334 + visci(1,3)*surfbo(3,ifac))**2 & 1335 + ( visci(2,1)*surfbo(1,ifac) & 1336 + visci(2,2)*surfbo(2,ifac) & 1337 + visci(2,3)*surfbo(3,ifac))**2 & 1338 + ( visci(3,1)*surfbo(1,ifac) & 1339 + visci(3,2)*surfbo(2,ifac) & 1340 + visci(3,3)*surfbo(3,ifac))**2 1341 1342 ! if.ki.s 1343 fikis = ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,1) & 1344 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,1) & 1345 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,1) & 1346 )*surfbo(1,ifac) & 1347 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,2) & 1348 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,2) & 1349 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,2) & 1350 )*surfbo(2,ifac) & 1351 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,3) & 1352 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,3) & 1353 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,3) & 1354 )*surfbo(3,ifac) 1355 1356 distfi = distb(ifac) 1357 1358 ! take i" so that i"f= eps*||fi||*ki.n when j" is in cell rji 1359 ! nb: eps =1.d-1 must be consistent with vitens.f90 1360 fikis = max(fikis, 1.d-1*sqrt(viscis)*distfi) 1361 1362 hint = viscis/surfbn(ifac)/fikis 1363 if (ivofmt.gt.0) hint = hint/crom(iel) 1364 1365 endif 1366 1367 ! on doit remodifier la valeur du dirichlet de pression de maniere 1368 ! a retrouver p*. car dans typecl.f90 on a travaille avec la pression 1369 ! totale fournie par l'utilisateur : ptotale= p*+ rho.g.r 1370 ! en compressible, on laisse rcodcl tel quel 1371 1372 ! dirichlet boundary condition 1373 !----------------------------- 1374 1375 if (icodcl(ifac,ipr).eq.1) then 1376 1377 hext = rcodcl(ifac,ipr,2) 1378 pimp = rcodcl(ifac,ipr,1) 1379 1380 call set_dirichlet_scalar & 1381 ( coefap(ifac), cofafp(ifac), & 1382 coefbp(ifac), cofbfp(ifac), & 1383 pimp , hint , hext ) 1384 1385 endif 1386 1387 ! neumann boundary conditions 1388 !---------------------------- 1389 1390 if (icodcl(ifac,ipr).eq.3) then 1391 1392 dimp = rcodcl(ifac,ipr,3) 1393 1394 call set_neumann_scalar & 1395 ( coefap(ifac), cofafp(ifac), & 1396 coefbp(ifac), cofbfp(ifac), & 1397 dimp , hint ) 1398 1399 ! convective boundary conditions 1400 !------------------------------- 1401 1402 elseif (icodcl(ifac,ipr).eq.2) then 1403 1404 pimp = rcodcl(ifac,ipr,1) 1405 cfl = rcodcl(ifac,ipr,2) 1406 1407 call set_convective_outlet_scalar & 1408 ( coefap(ifac), cofafp(ifac), & 1409 coefbp(ifac), cofbfp(ifac), & 1410 pimp , cfl , hint ) 1411 1412 ! Boundary value proportional to boundary cell value 1413 !--------------------------------------------------- 1414 1415 elseif (icodcl(ifac,ipr).eq.10) then 1416 1417 pinf = rcodcl(ifac,ipr,1) 1418 ratio = rcodcl(ifac,ipr,2) 1419 1420 call set_affine_function_scalar & 1421 ( coefap(ifac), cofafp(ifac), & 1422 coefbp(ifac), cofbfp(ifac), & 1423 pinf , ratio , hint ) 1424 1425 ! Imposed value for the convection operator is proportional to boundary 1426 ! cell value, imposed flux for diffusion 1427 !---------------------------------------------------------------------- 1428 1429 elseif (icodcl(ifac,ipr).eq.12) then 1430 1431 pinf = rcodcl(ifac,ipr,1) 1432 ratio = rcodcl(ifac,ipr,2) 1433 dimp = rcodcl(ifac,ipr,3) 1434 1435 call set_affine_function_conv_neumann_diff_scalar & 1436 ( coefap(ifac), cofafp(ifac), & 1437 coefbp(ifac), cofbfp(ifac), & 1438 pinf , ratio , dimp ) 1439 1440 ! Imposed value for the convection operator, imposed flux for diffusion 1441 !---------------------------------------------------------------------- 1442 1443 elseif (icodcl(ifac,ipr).eq.13) then 1444 1445 pimp = rcodcl(ifac,ipr,1) 1446 dimp = rcodcl(ifac,ipr,3) 1447 1448 call set_dirichlet_conv_neumann_diff_scalar & 1449 ( coefap(ifac), cofafp(ifac), & 1450 coefbp(ifac), cofbfp(ifac), & 1451 pimp , dimp ) 1452 1453 ! Neumann for the convection operator, zero flux for diffusion 1454 !---------------------------------------------------------------------- 1455 1456 elseif (icodcl(ifac,ipr).eq.15) then 1457 1458 dimp = rcodcl(ifac,ipr,3) 1459 1460 call set_neumann_conv_h_neumann_diff_scalar & 1461 ( coefap(ifac), cofafp(ifac), & 1462 coefbp(ifac), cofbfp(ifac), & 1463 dimp , hint ) 1464 1465 endif 1466 1467enddo 1468 1469!=============================================================================== 1470! 11. void fraction (VOF): dirichlet and neumann and convective outlet 1471!=============================================================================== 1472 1473if (ivofmt.gt.0) then 1474 1475 call field_get_coefa_s(ivarfl(ivolf2), coefap) 1476 call field_get_coefb_s(ivarfl(ivolf2), coefbp) 1477 call field_get_coefaf_s(ivarfl(ivolf2), cofafp) 1478 call field_get_coefbf_s(ivarfl(ivolf2), cofbfp) 1479 1480 do ifac = 1, nfabor 1481 1482 ! hint is unused since there is no diffusion for the void fraction 1483 hint = 1.d0 1484 1485 ! dirichlet boundary condition 1486 !----------------------------- 1487 1488 if (icodcl(ifac,ivolf2).eq.1) then 1489 1490 pimp = rcodcl(ifac,ivolf2,1) 1491 hext = rcodcl(ifac,ivolf2,2) 1492 1493 call set_dirichlet_scalar & 1494 ( coefap(ifac), cofafp(ifac), & 1495 coefbp(ifac), cofbfp(ifac), & 1496 pimp , hint , hext ) 1497 1498 endif 1499 1500 ! neumann boundary conditions 1501 !---------------------------- 1502 1503 if (icodcl(ifac,ivolf2).eq.3) then 1504 1505 dimp = rcodcl(ifac,ivolf2,3) 1506 1507 call set_neumann_scalar & 1508 ( coefap(ifac), cofafp(ifac), & 1509 coefbp(ifac), cofbfp(ifac), & 1510 dimp , hint ) 1511 1512 ! convective boundary conditions 1513 !------------------------------- 1514 1515 elseif (icodcl(ifac,ivolf2).eq.2) then 1516 1517 pimp = rcodcl(ifac,ivolf2,1) 1518 cfl = rcodcl(ifac,ivolf2,2) 1519 1520 call set_convective_outlet_scalar & 1521 ( coefap(ifac), cofafp(ifac), & 1522 coefbp(ifac), cofbfp(ifac), & 1523 pimp , cfl , hint ) 1524 1525 endif 1526 1527 enddo 1528 1529endif 1530 1531!=============================================================================== 1532! 12. turbulent quantities: dirichlet and neumann and convective outlet 1533!=============================================================================== 1534 1535! ---> k-epsilon and k-omega 1536 1537if (itytur.eq.2.or.iturb.eq.60) then 1538 1539 do ii = 1, 2 1540 1541 ! pour le k-omega, on met les valeurs sigma_k2 et sigma_w2 car ce terme 1542 ! ne concerne en pratique que les entrees (pas de pb en paroi ou en flux 1543 ! nul) 1544 if (ii.eq.1.and.itytur.eq.2) then 1545 ivar = ik 1546 call field_get_key_double(ivarfl(ik), ksigmas, sigma) 1547 elseif (ii.eq.1.and.iturb.eq.60) then 1548 ivar = ik 1549 sigma = ckwsk2 !fixme it is not consistent with the model 1550 elseif (itytur.eq.2) then 1551 ivar = iep 1552 call field_get_key_double(ivarfl(iep), ksigmas, sigma) 1553 else 1554 ivar = iomg 1555 sigma = ckwsw2 !fixme it is not consistent with the model 1556 endif 1557 1558 call field_get_coefa_s(ivarfl(ivar), coefap) 1559 call field_get_coefb_s(ivarfl(ivar), coefbp) 1560 call field_get_coefaf_s(ivarfl(ivar), cofafp) 1561 call field_get_coefbf_s(ivarfl(ivar), cofbfp) 1562 1563 do ifac = 1, nfabor 1564 1565 iel = ifabor(ifac) 1566 1567 ! --- physical propreties 1568 visclc = viscl(iel) 1569 visctc = visct(iel) 1570 1571 ! --- geometrical quantities 1572 distbf = distb(ifac) 1573 1574 hint = (visclc+visctc/sigma)/distbf 1575 1576 ! dirichlet boundary condition 1577 !----------------------------- 1578 1579 if (icodcl(ifac,ivar).eq.1) then 1580 1581 pimp = rcodcl(ifac,ivar,1) 1582 hext = rcodcl(ifac,ivar,2) 1583 1584 call set_dirichlet_scalar & 1585 ( coefap(ifac), cofafp(ifac), & 1586 coefbp(ifac), cofbfp(ifac), & 1587 pimp , hint , hext ) 1588 1589 endif 1590 1591 ! neumann boundary conditions 1592 !---------------------------- 1593 1594 if (icodcl(ifac,ivar).eq.3) then 1595 1596 dimp = rcodcl(ifac,ivar,3) 1597 1598 call set_neumann_scalar & 1599 ( coefap(ifac), cofafp(ifac), & 1600 coefbp(ifac), cofbfp(ifac), & 1601 dimp , hint ) 1602 1603 ! convective boundary conditions 1604 !------------------------------- 1605 1606 elseif (icodcl(ifac,ivar).eq.2) then 1607 1608 pimp = rcodcl(ifac,ivar,1) 1609 cfl = rcodcl(ifac,ivar,2) 1610 1611 call set_convective_outlet_scalar & 1612 ( coefap(ifac), cofafp(ifac), & 1613 coefbp(ifac), cofbfp(ifac), & 1614 pimp , cfl , hint ) 1615 1616 ! imposed value for the convection operator, imposed flux for diffusion 1617 !---------------------------------------------------------------------- 1618 1619 elseif (icodcl(ifac,ivar).eq.13) then 1620 1621 pimp = rcodcl(ifac,ivar,1) 1622 dimp = rcodcl(ifac,ivar,3) 1623 1624 call set_dirichlet_conv_neumann_diff_scalar & 1625 ( coefap(ifac), cofafp(ifac), & 1626 coefbp(ifac), cofbfp(ifac), & 1627 pimp , dimp ) 1628 1629 1630 endif 1631 1632 enddo 1633 1634 enddo 1635 1636! ---> rij-epsilon 1637elseif (itytur.eq.3) then 1638 1639 ! --> rij 1640 ivar = irij 1641 call field_get_coefa_v(ivarfl(irij), coefats) 1642 call field_get_coefb_v(ivarfl(irij), coefbts) 1643 call field_get_coefaf_v(ivarfl(irij), cofafts) 1644 call field_get_coefbf_v(ivarfl(irij), cofbfts) 1645 call field_get_coefad_v(ivarfl(irij), cofadts) 1646 call field_get_coefbd_v(ivarfl(irij), cofbdts) 1647 1648 call field_get_key_struct_var_cal_opt(ivarfl(irij), vcopt) 1649 1650 if (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 1651 call field_get_val_v(ivsten, visten) 1652 endif 1653 1654 do ifac = 1, nfabor 1655 1656 iel = ifabor(ifac) 1657 1658 ! --- physical propreties 1659 visclc = viscl(iel) 1660 1661 ! --- geometrical quantities 1662 distbf = distb(ifac) 1663 1664 ! symmetric tensor diffusivity (daly harlow - ggdh) TODO 1665 if (iand(vcopt%idften, ANISOTROPIC_RIGHT_DIFFUSION).ne.0) then 1666 1667 visci(1,1) = visclc + visten(1,iel) 1668 visci(2,2) = visclc + visten(2,iel) 1669 visci(3,3) = visclc + visten(3,iel) 1670 visci(1,2) = visten(4,iel) 1671 visci(2,1) = visten(4,iel) 1672 visci(2,3) = visten(5,iel) 1673 visci(3,2) = visten(5,iel) 1674 visci(1,3) = visten(6,iel) 1675 visci(3,1) = visten(6,iel) 1676 1677 ! ||ki.s||^2 1678 viscis = ( visci(1,1)*surfbo(1,ifac) & 1679 + visci(1,2)*surfbo(2,ifac) & 1680 + visci(1,3)*surfbo(3,ifac))**2 & 1681 + ( visci(2,1)*surfbo(1,ifac) & 1682 + visci(2,2)*surfbo(2,ifac) & 1683 + visci(2,3)*surfbo(3,ifac))**2 & 1684 + ( visci(3,1)*surfbo(1,ifac) & 1685 + visci(3,2)*surfbo(2,ifac) & 1686 + visci(3,3)*surfbo(3,ifac))**2 1687 1688 ! if.ki.s 1689 fikis = ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,1) & 1690 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,1) & 1691 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,1) & 1692 )*surfbo(1,ifac) & 1693 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,2) & 1694 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,2) & 1695 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,2) & 1696 )*surfbo(2,ifac) & 1697 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,3) & 1698 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,3) & 1699 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,3) & 1700 )*surfbo(3,ifac) 1701 1702 distfi = distb(ifac) 1703 1704 ! take i" so that i"f= eps*||fi||*ki.n when j" is in cell rji 1705 ! nb: eps =1.d-1 must be consistent with vitens.f90 1706 fikis = max(fikis, 1.d-1*sqrt(viscis)*distfi) 1707 1708 hint = viscis/surfbn(ifac)/fikis 1709 1710 ! scalar diffusivity 1711 else 1712 visctc = visct(iel) 1713 hint = (visclc+visctc*csrij/cmu)/distbf 1714 endif 1715 1716 do isou = 1, dimrij 1717 ivar = irij + isou - 1 1718 1719 ! allocate(pimpts(dimrij)) 1720 ! allocate(hextts(dimrij)) 1721 ! allocate(qimpts(dimrij)) 1722 ! allocate(cflts(dimrij)) 1723 ! Dirichlet Boundary Condition 1724 !----------------------------- 1725 1726 if (icodcl(ifac,ivar).eq.1) then 1727 pimpts(isou) = rcodcl(ifac,ivar,1) 1728 hextts(isou) = rcodcl(ifac,ivar,2) 1729 1730 call set_dirichlet_tensor & 1731 ( coefats(:, ifac), cofafts(:,ifac), & 1732 coefbts(:,:,ifac), cofbfts(:,:,ifac), & 1733 pimpts , hint , hextts ) 1734 1735 ! Boundary conditions for the momentum equation 1736 cofadts(isou, ifac) = coefats(isou, ifac) 1737 cofbdts(isou, isou, ifac) = coefbts(isou, isou, ifac) 1738 endif 1739 1740 ! Neumann Boundary Conditions 1741 !---------------------------- 1742 1743 if (icodcl(ifac,ivar).eq.3) then 1744 qimpts(isou) = rcodcl(ifac,ivar,3) 1745 1746 call set_neumann_tensor & 1747 ( coefats(:,ifac), cofafts(:,ifac), & 1748 coefbts(:,:,ifac), cofbfts(:,:,ifac), & 1749 qimpts , hint ) 1750 1751 ! Boundary conditions for the momentum equation 1752 cofadts(isou, ifac) = coefats(isou, ifac) 1753 cofbdts(isou, isou, ifac) = coefbts(isou ,isou ,ifac) 1754 1755 ! Convective Boundary Conditions 1756 !------------------------------- 1757 1758 elseif (icodcl(ifac,ivar).eq.2) then 1759 pimpts(isou) = rcodcl(ifac,ivar,1) 1760 cflts(isou) = rcodcl(ifac,ivar,2) 1761 1762 call set_convective_outlet_tensor & 1763 ( coefats(:, ifac), cofafts(:, ifac), & 1764 coefbts(:,:, ifac), cofbfts(:,:, ifac), & 1765 pimpts , cflts , hint ) 1766 ! Boundary conditions for the momentum equation 1767 cofadts(isou, ifac) = coefats(isou, ifac) 1768 cofbdts(isou, isou, ifac) = coefbts(isou, isou, ifac) 1769 1770 ! Imposed value for the convection operator, imposed flux for diffusion 1771 !---------------------------------------------------------------------- 1772 1773 elseif (icodcl(ifac,ivar).eq.13) then 1774 1775 pimpts(isou) = rcodcl(ifac,ivar,1) 1776 qimpts(isou) = rcodcl(ifac,ivar,3) 1777 1778 call set_dirichlet_conv_neumann_diff_tensor & 1779 ( coefats(:, ifac), cofafts(:, ifac), & 1780 coefbts(:,:, ifac), cofbfts(:,:, ifac), & 1781 pimpts , qimpts ) 1782 1783 ! Boundary conditions for the momentum equation 1784 cofadts(isou, ifac) = coefats(isou, ifac) 1785 cofbdts(isou, isou, ifac) = coefbts(isou, isou, ifac) 1786 1787 endif 1788 1789 enddo 1790 enddo 1791 1792 ! --> epsilon 1793 1794 ivar = iep 1795 1796 call field_get_coefa_s(ivarfl(ivar), coefap) 1797 call field_get_coefb_s(ivarfl(ivar), coefbp) 1798 call field_get_coefaf_s(ivarfl(ivar), cofafp) 1799 call field_get_coefbf_s(ivarfl(ivar), cofbfp) 1800 call field_get_key_double(ivarfl(iep), ksigmas, sigmae) 1801 1802 call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 1803 1804 if (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 1805 call field_get_val_v(ivsten, visten) 1806 endif 1807 1808 do ifac = 1, nfabor 1809 1810 iel = ifabor(ifac) 1811 1812 ! --- Physical Propreties 1813 visclc = viscl(iel) 1814 visctc = visct(iel) 1815 1816 ! --- Geometrical quantities 1817 distbf = distb(ifac) 1818 1819 ! Symmetric tensor diffusivity (Daly Harlow -- GGDH) 1820 if (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 1821 1822 visci(1,1) = visclc + visten(1,iel)/sigmae 1823 visci(2,2) = visclc + visten(2,iel)/sigmae 1824 visci(3,3) = visclc + visten(3,iel)/sigmae 1825 visci(1,2) = visten(4,iel)/sigmae 1826 visci(2,1) = visten(4,iel)/sigmae 1827 visci(2,3) = visten(5,iel)/sigmae 1828 visci(3,2) = visten(5,iel)/sigmae 1829 visci(1,3) = visten(6,iel)/sigmae 1830 visci(3,1) = visten(6,iel)/sigmae 1831 1832 ! ||Ki.S||^2 1833 viscis = ( visci(1,1)*surfbo(1,ifac) & 1834 + visci(1,2)*surfbo(2,ifac) & 1835 + visci(1,3)*surfbo(3,ifac))**2 & 1836 + ( visci(2,1)*surfbo(1,ifac) & 1837 + visci(2,2)*surfbo(2,ifac) & 1838 + visci(2,3)*surfbo(3,ifac))**2 & 1839 + ( visci(3,1)*surfbo(1,ifac) & 1840 + visci(3,2)*surfbo(2,ifac) & 1841 + visci(3,3)*surfbo(3,ifac))**2 1842 1843 ! IF.Ki.S 1844 fikis = ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,1) & 1845 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,1) & 1846 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,1) & 1847 )*surfbo(1,ifac) & 1848 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,2) & 1849 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,2) & 1850 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,2) & 1851 )*surfbo(2,ifac) & 1852 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,3) & 1853 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,3) & 1854 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,3) & 1855 )*surfbo(3,ifac) 1856 1857 distfi = distb(ifac) 1858 1859 ! Take I" so that I"F= eps*||FI||*Ki.n when J" is in cell rji 1860 ! NB: eps =1.d-1 must be consistent with vitens.f90 1861 fikis = max(fikis, 1.d-1*sqrt(viscis)*distfi) 1862 1863 hint = viscis/surfbn(ifac)/fikis 1864 1865 ! Scalar diffusivity 1866 else 1867 hint = (visclc+visctc/sigmae)/distbf 1868 endif 1869 1870 ! Dirichlet Boundary Condition 1871 !----------------------------- 1872 1873 if (icodcl(ifac,ivar).eq.1) then 1874 1875 pimp = rcodcl(ifac,ivar,1) 1876 hext = rcodcl(ifac,ivar,2) 1877 1878 call set_dirichlet_scalar & 1879 ( coefap(ifac), cofafp(ifac), & 1880 coefbp(ifac), cofbfp(ifac), & 1881 pimp , hint , hext ) 1882 1883 endif 1884 1885 ! Neumann Boundary Conditions 1886 !---------------------------- 1887 1888 if (icodcl(ifac,ivar).eq.3) then 1889 1890 dimp = rcodcl(ifac,ivar,3) 1891 1892 call set_neumann_scalar & 1893 ( coefap(ifac), cofafp(ifac), & 1894 coefbp(ifac), cofbfp(ifac), & 1895 dimp , hint ) 1896 1897 ! Convective Boundary Conditions 1898 !------------------------------- 1899 1900 elseif (icodcl(ifac,ivar).eq.2) then 1901 1902 pimp = rcodcl(ifac,ivar,1) 1903 cfl = rcodcl(ifac,ivar,2) 1904 1905 call set_convective_outlet_scalar & 1906 ( coefap(ifac), cofafp(ifac), & 1907 coefbp(ifac), cofbfp(ifac), & 1908 pimp , cfl , hint ) 1909 1910 1911 ! Imposed value for the convection operator, imposed flux for diffusion 1912 !---------------------------------------------------------------------- 1913 1914 elseif (icodcl(ifac,ivar).eq.13) then 1915 1916 pimp = rcodcl(ifac,ivar,1) 1917 dimp = rcodcl(ifac,ivar,3) 1918 1919 call set_dirichlet_conv_neumann_diff_scalar & 1920 ( coefap(ifac), cofafp(ifac), & 1921 coefbp(ifac), cofbfp(ifac), & 1922 pimp , dimp ) 1923 1924 endif 1925 1926 enddo 1927 1928 ! --> alpha for the EBRSM 1929 1930 if (iturb.eq.32) then 1931 ivar = ial 1932 1933 call field_get_coefa_s(ivarfl(ivar), coefap) 1934 call field_get_coefb_s(ivarfl(ivar), coefbp) 1935 call field_get_coefaf_s(ivarfl(ivar), cofafp) 1936 call field_get_coefbf_s(ivarfl(ivar), cofbfp) 1937 1938 do ifac = 1, nfabor 1939 1940 iel = ifabor(ifac) 1941 1942 distbf = distb(ifac) 1943 1944 hint = 1.d0/distbf 1945 1946 ! Dirichlet Boundary Condition 1947 !----------------------------- 1948 1949 if (icodcl(ifac,ivar).eq.1) then 1950 1951 pimp = rcodcl(ifac,ivar,1) 1952 hext = rcodcl(ifac,ivar,2) 1953 1954 call set_dirichlet_scalar & 1955 ( coefap(ifac), cofafp(ifac), & 1956 coefbp(ifac), cofbfp(ifac), & 1957 pimp , hint , hext ) 1958 1959 endif 1960 1961 ! Neumann Boundary Conditions 1962 !---------------------------- 1963 1964 if (icodcl(ifac,ivar).eq.3) then 1965 1966 dimp = rcodcl(ifac,ivar,3) 1967 1968 call set_neumann_scalar & 1969 ( coefap(ifac), cofafp(ifac), & 1970 coefbp(ifac), cofbfp(ifac), & 1971 dimp , hint ) 1972 1973 ! Convective Boundary Conditions 1974 !------------------------------- 1975 1976 elseif (icodcl(ifac,ivar).eq.2) then 1977 1978 pimp = rcodcl(ifac,ivar,1) 1979 cfl = rcodcl(ifac,ivar,2) 1980 1981 call set_convective_outlet_scalar & 1982 ( coefap(ifac), cofafp(ifac), & 1983 coefbp(ifac), cofbfp(ifac), & 1984 pimp , cfl , hint ) 1985 1986 ! Imposed value for the convection operator, imposed flux for diffusion 1987 !---------------------------------------------------------------------- 1988 1989 elseif (icodcl(ifac,ivar).eq.13) then 1990 1991 pimp = rcodcl(ifac,ivar,1) 1992 dimp = rcodcl(ifac,ivar,3) 1993 1994 call set_dirichlet_conv_neumann_diff_scalar & 1995 ( coefap(ifac), cofafp(ifac), & 1996 coefbp(ifac), cofbfp(ifac), & 1997 pimp , dimp ) 1998 1999 2000 endif 2001 2002 enddo 2003 endif 2004 2005! ---> v2f type models (phi_bar and Bl-v2/k) 2006 2007elseif (itytur.eq.5) then 2008 2009 ! --> k, epsilon and phi 2010 do ii = 1, 3 2011 2012 if (ii.eq.1) then 2013 ivar = ik 2014 call field_get_key_double(ivarfl(ik), ksigmas, sigma) 2015 elseif (ii.eq.2) then 2016 ivar = iep 2017 call field_get_key_double(ivarfl(iep), ksigmas, sigma) 2018 else 2019 ivar = iphi 2020 call field_get_key_double(ivarfl(iphi), ksigmas, sigma) 2021 endif 2022 2023 call field_get_coefa_s(ivarfl(ivar), coefap) 2024 call field_get_coefb_s(ivarfl(ivar), coefbp) 2025 call field_get_coefaf_s(ivarfl(ivar), cofafp) 2026 call field_get_coefbf_s(ivarfl(ivar), cofbfp) 2027 2028 do ifac = 1, nfabor 2029 2030 iel = ifabor(ifac) 2031 2032 ! --- Physical Properties 2033 visclc = viscl(iel) 2034 visctc = visct(iel) 2035 2036 ! --- Geometrical quantities 2037 distbf = distb(ifac) 2038 2039 hint = (visclc+visctc/sigma)/distbf 2040 2041 ! Dirichlet Boundary Condition 2042 !----------------------------- 2043 2044 if (icodcl(ifac,ivar).eq.1) then 2045 2046 pimp = rcodcl(ifac,ivar,1) 2047 hext = rcodcl(ifac,ivar,2) 2048 2049 call set_dirichlet_scalar & 2050 ( coefap(ifac), cofafp(ifac), & 2051 coefbp(ifac), cofbfp(ifac), & 2052 pimp , hint , hext ) 2053 2054 endif 2055 2056 ! Neumann Boundary Conditions 2057 !---------------------------- 2058 2059 if (icodcl(ifac,ivar).eq.3) then 2060 2061 dimp = rcodcl(ifac,ivar,3) 2062 2063 call set_neumann_scalar & 2064 ( coefap(ifac), cofafp(ifac), & 2065 coefbp(ifac), cofbfp(ifac), & 2066 dimp , hint ) 2067 2068 ! Convective Boundary Conditions 2069 !------------------------------- 2070 2071 elseif (icodcl(ifac,ivar).eq.2) then 2072 2073 pimp = rcodcl(ifac,ivar,1) 2074 cfl = rcodcl(ifac,ivar,2) 2075 2076 call set_convective_outlet_scalar & 2077 ( coefap(ifac), cofafp(ifac), & 2078 coefbp(ifac), cofbfp(ifac), & 2079 pimp , cfl , hint ) 2080 2081 ! Imposed value for the convection operator, imposed flux for diffusion 2082 !---------------------------------------------------------------------- 2083 2084 elseif (icodcl(ifac,ivar).eq.13) then 2085 2086 pimp = rcodcl(ifac,ivar,1) 2087 dimp = rcodcl(ifac,ivar,3) 2088 2089 call set_dirichlet_conv_neumann_diff_scalar & 2090 ( coefap(ifac), cofafp(ifac), & 2091 coefbp(ifac), cofbfp(ifac), & 2092 pimp , dimp ) 2093 2094 2095 endif 2096 2097 enddo 2098 2099 enddo 2100 2101 if (iturb.eq.50) then 2102 2103 ! --> FB 2104 2105 ivar = ifb 2106 2107 call field_get_coefa_s(ivarfl(ivar), coefap) 2108 call field_get_coefb_s(ivarfl(ivar), coefbp) 2109 call field_get_coefaf_s(ivarfl(ivar), cofafp) 2110 call field_get_coefbf_s(ivarfl(ivar), cofbfp) 2111 2112 do ifac = 1, nfabor 2113 2114 ! --- Physical Properties 2115 visclc = 1.d0 2116 2117 ! --- Geometrical quantities 2118 distbf = distb(ifac) 2119 2120 hint = visclc/distbf 2121 2122 ! Dirichlet Boundary Condition 2123 !----------------------------- 2124 2125 if (icodcl(ifac,ivar).eq.1) then 2126 2127 pimp = rcodcl(ifac,ivar,1) 2128 hext = rcodcl(ifac,ivar,2) 2129 2130 call set_dirichlet_scalar & 2131 ( coefap(ifac), cofafp(ifac), & 2132 coefbp(ifac), cofbfp(ifac), & 2133 pimp , hint , hext ) 2134 2135 endif 2136 2137 ! Neumann Boundary Conditions 2138 !---------------------------- 2139 2140 if (icodcl(ifac,ivar).eq.3) then 2141 2142 dimp = rcodcl(ifac,ivar,3) 2143 2144 call set_neumann_scalar & 2145 ( coefap(ifac), cofafp(ifac), & 2146 coefbp(ifac), cofbfp(ifac), & 2147 dimp , hint ) 2148 2149 ! Convective Boundary Conditions 2150 !------------------------------- 2151 2152 elseif (icodcl(ifac,ivar).eq.2) then 2153 2154 pimp = rcodcl(ifac,ivar,1) 2155 cfl = rcodcl(ifac,ivar,2) 2156 2157 call set_convective_outlet_scalar & 2158 ( coefap(ifac), cofafp(ifac), & 2159 coefbp(ifac), cofbfp(ifac), & 2160 pimp , cfl , hint ) 2161 2162 ! Imposed value for the convection operator, imposed flux for diffusion 2163 !---------------------------------------------------------------------- 2164 2165 elseif (icodcl(ifac,ivar).eq.13) then 2166 2167 pimp = rcodcl(ifac,ivar,1) 2168 dimp = rcodcl(ifac,ivar,3) 2169 2170 call set_dirichlet_conv_neumann_diff_scalar & 2171 ( coefap(ifac), cofafp(ifac), & 2172 coefbp(ifac), cofbfp(ifac), & 2173 pimp , dimp ) 2174 2175 2176 endif 2177 2178 enddo 2179 2180 elseif (iturb.eq.51) then 2181 2182 ! --> alpha 2183 2184 ivar = ial 2185 2186 call field_get_coefa_s(ivarfl(ivar), coefap) 2187 call field_get_coefb_s(ivarfl(ivar), coefbp) 2188 call field_get_coefaf_s(ivarfl(ivar), cofafp) 2189 call field_get_coefbf_s(ivarfl(ivar), cofbfp) 2190 2191 do ifac = 1, nfabor 2192 2193 ! --- Physical Propreties 2194 visclc = 1.d0 2195 2196 ! --- Geometrical quantities 2197 distbf = distb(ifac) 2198 2199 hint = visclc/distbf 2200 2201 ! Dirichlet Boundary Condition 2202 !----------------------------- 2203 2204 if (icodcl(ifac,ivar).eq.1) then 2205 2206 pimp = rcodcl(ifac,ivar,1) 2207 hext = rcodcl(ifac,ivar,2) 2208 2209 call set_dirichlet_scalar & 2210 ( coefap(ifac), cofafp(ifac), & 2211 coefbp(ifac), cofbfp(ifac), & 2212 pimp , hint , hext ) 2213 2214 endif 2215 2216 ! Neumann Boundary Conditions 2217 !---------------------------- 2218 2219 if (icodcl(ifac,ivar).eq.3) then 2220 2221 dimp = rcodcl(ifac,ivar,3) 2222 2223 call set_neumann_scalar & 2224 ( coefap(ifac), cofafp(ifac), & 2225 coefbp(ifac), cofbfp(ifac), & 2226 dimp , hint ) 2227 2228 ! Convective Boundary Conditions 2229 !------------------------------- 2230 2231 elseif (icodcl(ifac,ivar).eq.2) then 2232 2233 pimp = rcodcl(ifac,ivar,1) 2234 cfl = rcodcl(ifac,ivar,2) 2235 2236 call set_convective_outlet_scalar & 2237 ( coefap(ifac), cofafp(ifac), & 2238 coefbp(ifac), cofbfp(ifac), & 2239 pimp , cfl , hint ) 2240 2241 ! Imposed value for the convection operator, imposed flux for diffusion 2242 !---------------------------------------------------------------------- 2243 2244 elseif (icodcl(ifac,ivar).eq.13) then 2245 2246 pimp = rcodcl(ifac,ivar,1) 2247 dimp = rcodcl(ifac,ivar,3) 2248 2249 call set_dirichlet_conv_neumann_diff_scalar & 2250 ( coefap(ifac), cofafp(ifac), & 2251 coefbp(ifac), cofbfp(ifac), & 2252 pimp , dimp ) 2253 2254 2255 endif 2256 2257 enddo 2258 2259 endif 2260 2261! ---> Spalart Allmaras 2262 2263elseif (iturb.eq.70) then 2264 2265 ivar = inusa 2266 2267 call field_get_coefa_s(ivarfl(ivar), coefap) 2268 call field_get_coefb_s(ivarfl(ivar), coefbp) 2269 call field_get_coefaf_s(ivarfl(ivar), cofafp) 2270 call field_get_coefbf_s(ivarfl(ivar), cofbfp) 2271 2272 do ifac = 1, nfabor 2273 2274 iel = ifabor(ifac) 2275 2276 ! --- Physical Propreties 2277 visclc = viscl(iel) 2278 2279 ! --- Geometrical quantities 2280 distbf = distb(ifac) 2281 2282 hint = visclc/distbf 2283 2284 ! Dirichlet Boundary Condition 2285 !----------------------------- 2286 2287 if (icodcl(ifac,ivar).eq.1) then 2288 2289 pimp = rcodcl(ifac,ivar,1) 2290 hext = rcodcl(ifac,ivar,2) 2291 2292 call set_dirichlet_scalar & 2293 ( coefap(ifac), cofafp(ifac), & 2294 coefbp(ifac), cofbfp(ifac), & 2295 pimp , hint , hext ) 2296 2297 endif 2298 2299 ! Neumann Boundary Conditions 2300 !---------------------------- 2301 2302 if (icodcl(ifac,ivar).eq.3) then 2303 2304 dimp = rcodcl(ifac,ivar,3) 2305 2306 call set_neumann_scalar & 2307 ( coefap(ifac), cofafp(ifac), & 2308 coefbp(ifac), cofbfp(ifac), & 2309 dimp , hint ) 2310 2311 ! Convective Boundary Conditions 2312 !------------------------------- 2313 2314 elseif (icodcl(ifac,ivar).eq.2) then 2315 2316 pimp = rcodcl(ifac,ivar,1) 2317 cfl = rcodcl(ifac,ivar,2) 2318 2319 call set_convective_outlet_scalar & 2320 ( coefap(ifac), cofafp(ifac), & 2321 coefbp(ifac), cofbfp(ifac), & 2322 pimp , cfl , hint ) 2323 2324 ! Imposed value for the convection operator, imposed flux for diffusion 2325 !---------------------------------------------------------------------- 2326 2327 elseif (icodcl(ifac,ivar).eq.13) then 2328 2329 pimp = rcodcl(ifac,ivar,1) 2330 dimp = rcodcl(ifac,ivar,3) 2331 2332 call set_dirichlet_conv_neumann_diff_scalar & 2333 ( coefap(ifac), cofafp(ifac), & 2334 coefbp(ifac), cofbfp(ifac), & 2335 pimp , dimp ) 2336 2337 2338 endif 2339 2340 enddo 2341 2342endif 2343 2344!=============================================================================== 2345! 13. Other scalars (except variances): 2346! Dirichlet and Neumann and convective outlet 2347!=============================================================================== 2348 2349if (nscal.ge.1) then 2350 2351 if(icp.ge.0) then 2352 call field_get_val_s(icp, cpro_cp) 2353 endif 2354 2355 if (ippmod(icompf).ge.0.and.icv.ge.0) then 2356 call field_get_val_s(icv, cpro_cv) 2357 endif 2358 2359 call field_get_key_id('turbulent_flux_model', kturt) 2360 2361 do ii = 1, nscal 2362 2363 ivar = isca(ii) 2364 2365 isvhbl = 0 2366 if (ii.eq.isvhb) then 2367 isvhbl = isvhb 2368 endif 2369 2370 call field_get_key_int (ivarfl(ivar), kivisl, ifcvsl) 2371 if (ifcvsl .ge. 0) then 2372 call field_get_val_s(ifcvsl, viscls) 2373 endif 2374 2375 call field_get_key_int(ivarfl(ivar), kscacp, iscacp) 2376 2377 2378 ! Get the turbulent flux model for the scalar 2379 call field_get_key_int(ivarfl(isca(ii)), kturt, turb_flux_model) 2380 turb_flux_model_type = turb_flux_model / 10 2381 2382 ! --- Indicateur de prise en compte de Cp ou non 2383 ! (selon si le scalaire (scalaire associe pour une fluctuation) 2384 ! doit etre ou non traite comme une temperature) 2385 ! Si le scalaire est une variance et que le 2386 ! scalaire associe n'est pas resolu, on suppose alors qu'il 2387 ! doit etre traite comme un scalaire passif (defaut IHCP = 0) 2388 ihcp = 0 2389 2390 iscal = ii 2391 if (iscavr(ii).gt.0) then 2392 iscal = iscavr(ii) 2393 endif 2394 2395 ! Reference diffusivity 2396 call field_get_key_double(ivarfl(isca(iscal)), kvisl0, visls_0) 2397 2398 if (iscacp.eq.1) then 2399 if(icp.ge.0) then 2400 ihcp = 2 2401 else 2402 ihcp = 1 2403 endif 2404 endif 2405 2406 call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 2407 2408 if (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0.or.turb_flux_model_type.eq.3) then 2409 if (iturb.ne.32.or.turb_flux_model_type.eq.3) then 2410 call field_get_val_v(ivsten, visten) 2411 else ! EBRSM and (GGDH or AFM) 2412 call field_get_val_v(ivstes, visten) 2413 endif 2414 endif 2415 2416 call field_get_key_double(ivarfl(isca(ii)), ksigmas, turb_schmidt) 2417 2418 ! Get boundary value (for post-processing) 2419 call field_get_key_int(ivarfl(ivar), kbfid, b_f_id) 2420 call field_get_dim(ivarfl(isca(iscal)), f_dim) 2421 2422 ! Scalar transported quantity 2423 if (f_dim.eq.1) then 2424 call field_get_coefa_s(ivarfl(ivar), coefap) 2425 call field_get_coefb_s(ivarfl(ivar), coefbp) 2426 call field_get_coefaf_s(ivarfl(ivar), cofafp) 2427 call field_get_coefbf_s(ivarfl(ivar), cofbfp) 2428 2429 ! if thermal variable has no boundary but temperature does, use it 2430 if (b_f_id .lt. 0 .and. ii.eq.iscalt .and. itherm.eq.2) then 2431 b_f_id = itempb 2432 endif 2433 2434 if (b_f_id .ge. 0) then 2435 call field_get_val_s(b_f_id, bvar_s) 2436 else 2437 bvar_s => null() 2438 endif 2439 2440 do ifac = 1, nfabor 2441 2442 iel = ifabor(ifac) 2443 2444 ! --- Physical Properties 2445 visctc = visct(iel) 2446 2447 ! --- Geometrical quantities 2448 distbf = distb(ifac) 2449 2450 ! --- Prise en compte de Cp ou CV 2451 ! (dans le Cas compressible ihcp=0) 2452 2453 cpp = 1.d0 2454 if (ihcp.eq.0) then 2455 cpp = 1.d0 2456 elseif (ihcp.eq.2) then 2457 cpp = cpro_cp(iel) 2458 elseif (ihcp.eq.1) then 2459 cpp = cp0 2460 endif 2461 2462 ! --- Viscosite variable ou non 2463 if (ifcvsl.lt.0) then 2464 rkl = visls_0 2465 else 2466 rkl = viscls(iel) 2467 endif 2468 2469 ! Scalar diffusivity 2470 if (iand(vcopt%idften, ISOTROPIC_DIFFUSION).ne.0) then 2471 if (ippmod(idarcy).eq.-1) then !FIXME 2472 hint = (rkl+vcopt%idifft*cpp*visctc/turb_schmidt)/distbf 2473 else ! idarcy = 1 2474 hint = rkl/distbf 2475 endif 2476 2477 ! Symmetric tensor diffusivity 2478 elseif (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 2479 2480 temp = vcopt%idifft*cpp*ctheta(ii)/csrij 2481 visci(1,1) = rkl + temp*visten(1,iel) 2482 visci(2,2) = rkl + temp*visten(2,iel) 2483 visci(3,3) = rkl + temp*visten(3,iel) 2484 visci(1,2) = temp*visten(4,iel) 2485 visci(2,1) = temp*visten(4,iel) 2486 visci(2,3) = temp*visten(5,iel) 2487 visci(3,2) = temp*visten(5,iel) 2488 visci(1,3) = temp*visten(6,iel) 2489 visci(3,1) = temp*visten(6,iel) 2490 2491 ! ||Ki.S||^2 2492 viscis = ( visci(1,1)*surfbo(1,ifac) & 2493 + visci(1,2)*surfbo(2,ifac) & 2494 + visci(1,3)*surfbo(3,ifac))**2 & 2495 + ( visci(2,1)*surfbo(1,ifac) & 2496 + visci(2,2)*surfbo(2,ifac) & 2497 + visci(2,3)*surfbo(3,ifac))**2 & 2498 + ( visci(3,1)*surfbo(1,ifac) & 2499 + visci(3,2)*surfbo(2,ifac) & 2500 + visci(3,3)*surfbo(3,ifac))**2 2501 2502 ! IF.Ki.S 2503 fikis = ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,1) & 2504 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,1) & 2505 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,1) & 2506 )*surfbo(1,ifac) & 2507 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,2) & 2508 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,2) & 2509 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,2) & 2510 )*surfbo(2,ifac) & 2511 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,3) & 2512 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,3) & 2513 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,3) & 2514 )*surfbo(3,ifac) 2515 2516 distfi = distb(ifac) 2517 2518 ! Take I" so that I"F= eps*||FI||*Ki.n when J" is in cell rji 2519 ! NB: eps =1.d-1 must be consistent with vitens.f90 2520 fikis = max(fikis, 1.d-1*sqrt(viscis)*distfi) 2521 2522 hint = viscis/surfbn(ifac)/fikis 2523 2524 endif 2525 2526 ! Dirichlet Boundary Condition 2527 !----------------------------- 2528 2529 if (icodcl(ifac,ivar).eq.1) then 2530 2531 pimp = rcodcl(ifac,ivar,1) 2532 hext = rcodcl(ifac,ivar,2) 2533 2534 call set_dirichlet_scalar & 2535 ( coefap(ifac), cofafp(ifac), & 2536 coefbp(ifac), cofbfp(ifac), & 2537 pimp , hint , hext ) 2538 2539 ! Store boundary value 2540 if (b_f_id .ge. 0) then 2541 bvar_s(ifac) = coefap(ifac) + coefbp(ifac) * bvar_s(ifac) 2542 endif 2543 endif 2544 2545 ! Neumann Boundary Conditions 2546 !---------------------------- 2547 2548 if (icodcl(ifac,ivar).eq.3) then 2549 2550 dimp = rcodcl(ifac,ivar,3) 2551 2552 call set_neumann_scalar & 2553 ( coefap(ifac), cofafp(ifac), & 2554 coefbp(ifac), cofbfp(ifac), & 2555 dimp , hint ) 2556 2557 ! Store boundary value 2558 if (b_f_id .ge. 0) then 2559 bvar_s(ifac) = coefap(ifac) + coefbp(ifac) * bvar_s(ifac) 2560 endif 2561 2562 ! Convective Boundary Conditions 2563 !------------------------------- 2564 2565 elseif (icodcl(ifac,ivar).eq.2 .and. iterns.le.1) then 2566 2567 pimp = rcodcl(ifac,ivar,1) 2568 cfl = rcodcl(ifac,ivar,2) 2569 2570 call set_convective_outlet_scalar & 2571 ( coefap(ifac), cofafp(ifac), & 2572 coefbp(ifac), cofbfp(ifac), & 2573 pimp , cfl , hint ) 2574 2575 ! Store boundary value 2576 if (b_f_id .ge. 0) then 2577 bvar_s(ifac) = coefap(ifac) + coefbp(ifac) * bvar_s(ifac) 2578 endif 2579 2580 ! Set total flux as a Robin condition 2581 !------------------------------------ 2582 2583 elseif (icodcl(ifac,ivar).eq.12) then 2584 2585 hext = rcodcl(ifac,ivar,2) 2586 dimp = rcodcl(ifac,ivar,3) 2587 2588 call set_total_flux & 2589 ( coefap(ifac), cofafp(ifac), & 2590 coefbp(ifac), cofbfp(ifac), & 2591 hext , dimp ) 2592 2593 ! Store boundary value 2594 if (b_f_id .ge. 0) then 2595 bvar_s(ifac) = coefap(ifac) + coefbp(ifac) * bvar_s(ifac) 2596 endif 2597 2598 ! Imposed value for the convection operator, imposed flux for diffusion 2599 !---------------------------------------------------------------------- 2600 2601 elseif (icodcl(ifac,ivar).eq.13) then 2602 2603 pimp = rcodcl(ifac,ivar,1) 2604 dimp = rcodcl(ifac,ivar,3) 2605 2606 call set_dirichlet_conv_neumann_diff_scalar & 2607 ( coefap(ifac), cofafp(ifac), & 2608 coefbp(ifac), cofbfp(ifac), & 2609 pimp , dimp ) 2610 2611 ! Store boundary value 2612 if (b_f_id .ge. 0) then 2613 bvar_s(ifac) = coefap(ifac) + coefbp(ifac) * bvar_s(ifac) 2614 endif 2615 2616 endif 2617 2618 ! Storage of the thermal exchange coefficient 2619 ! (conversion in case of energy or enthalpy) 2620 ! the exchange coefficient is in W/(m2 K) 2621 ! Useful for thermal coupling or radiative transfer 2622 if (icodcl(ifac,ivar).eq.1.or.icodcl(ifac,ivar).eq.3) then 2623 if (iirayo.ge.1 .and. ii.eq.iscalt.or.isvhbl.gt.0) then 2624 2625 ! Enthalpy 2626 if (itherm.eq.2) then 2627 ! If Cp is variable 2628 if (icp.ge.0) then 2629 exchange_coef = hint*cpro_cp(iel) 2630 else 2631 exchange_coef = hint*cp0 2632 endif 2633 2634 ! Total energy (compressible module) 2635 elseif (itherm.eq.3) then 2636 ! If Cv is variable 2637 if (icv.ge.0) then 2638 exchange_coef = hint*cpro_cv(iel) 2639 else 2640 exchange_coef = hint*cv0 2641 endif 2642 2643 ! Temperature 2644 elseif (iscacp.eq.1) then 2645 exchange_coef = hint 2646 endif 2647 endif 2648 2649 ! ---> Thermal coupling, store hint = lambda/d 2650 if (isvhbl.gt.0) hbord(ifac) = exchange_coef 2651 2652 ! ---> Radiative transfer 2653 if (iirayo.ge.1 .and. ii.eq.iscalt) then 2654 bhconv(ifac) = exchange_coef 2655 2656 ! The outgoing flux is stored (Q = h(Ti'-Tp): negative if 2657 ! gain for the fluid) in W/m2 2658 bfconv(ifac) = cofafp(ifac) + cofbfp(ifac)*theipb(ifac) 2659 endif 2660 2661 endif 2662 2663 ! Thermal heat flux boundary conditions 2664 if (turb_flux_model_type.eq.3) then 2665 2666 ! Name of the scalar ivar !TODO move outside of the loop 2667 call field_get_name(ivarfl(ivar), fname) 2668 2669 ! Index of the corresponding turbulent flux 2670 call field_get_id(trim(fname)//'_turbulent_flux', f_id) 2671 2672 call field_get_coefa_v(f_id,coefaut) 2673 call field_get_coefb_v(f_id,coefbut) 2674 call field_get_coefaf_v(f_id,cofafut) 2675 call field_get_coefbf_v(f_id,cofbfut) 2676 call field_get_coefad_v(f_id,cofarut) 2677 call field_get_coefbd_v(f_id,cofbrut) 2678 2679 ! --- Physical Properties 2680 visclc = viscl(iel) 2681 2682 ! --- Geometrical quantities 2683 distbf = distb(ifac) 2684 2685 if (ifcvsl.lt.0) then 2686 rkl = visls_0/cpp 2687 else 2688 rkl = viscls(iel)/cpp 2689 endif 2690 hintt(1) = 0.5d0*(visclc+rkl)/distbf & 2691 + visten(1,iel)*ctheta(iscal)/distbf/csrij !FIXME ctheta (iscal) 2692 hintt(2) = 0.5d0*(visclc+rkl)/distbf & 2693 + visten(2,iel)*ctheta(iscal)/distbf/csrij 2694 hintt(3) = 0.5d0*(visclc+rkl)/distbf & 2695 + visten(3,iel)*ctheta(iscal)/distbf/csrij 2696 hintt(4) = visten(4,iel)*ctheta(iscal)/distbf/csrij 2697 hintt(5) = visten(5,iel)*ctheta(iscal)/distbf/csrij 2698 hintt(6) = visten(6,iel)*ctheta(iscal)/distbf/csrij 2699 2700 ! Set pointer values of turbulent fluxes in icodcl 2701 call field_get_key_int(f_id, keyvar, iut) 2702 ivt = iut + 1 2703 iwt = iut + 2 2704 2705 ! Dirichlet Boundary Condition 2706 !----------------------------- 2707 2708 if (icodcl(ifac,iut).eq.1) then 2709 2710 pimpv(1) = rcodcl(ifac,iut,1) 2711 pimpv(2) = rcodcl(ifac,ivt,1) 2712 pimpv(3) = rcodcl(ifac,iwt,1) 2713 hextv(1) = rcodcl(ifac,iut,2) 2714 hextv(2) = rcodcl(ifac,ivt,2) 2715 hextv(3) = rcodcl(ifac,iwt,2) 2716 2717 call set_dirichlet_vector_aniso & 2718 ( coefaut(:,ifac) , cofafut(:,ifac) , & 2719 coefbut(:,:,ifac), cofbfut(:,:,ifac), & 2720 pimpv , hintt , hextv ) 2721 2722 ! Boundary conditions for thermal transport equation 2723 do isou = 1, 3 2724 cofarut(isou,ifac) = coefaut(isou,ifac) 2725 do jsou =1, 3 2726 cofbrut(isou,jsou,ifac) = coefbut(isou,jsou,ifac) 2727 enddo 2728 enddo 2729 2730 ! Neumann Boundary Conditions 2731 !---------------------------- 2732 2733 elseif (icodcl(ifac,iut).eq.3) then 2734 2735 qimpv(1) = rcodcl(ifac,iut,3) 2736 qimpv(2) = rcodcl(ifac,ivt,3) 2737 qimpv(3) = rcodcl(ifac,iwt,3) 2738 2739 call set_neumann_vector_aniso & 2740 ( coefaut(:,ifac) , cofafut(:,ifac) , & 2741 coefbut(:,:,ifac), cofbfut(:,:,ifac), & 2742 qimpv , hintt ) 2743 2744 ! Boundary conditions for thermal transport equation 2745 do isou = 1, 3 2746 cofarut(isou,ifac) = coefaut(isou,ifac) 2747 do jsou =1, 3 2748 cofbrut(isou,jsou,ifac) = coefbut(isou,jsou,ifac) 2749 enddo 2750 enddo 2751 2752 ! Convective Boundary Conditions 2753 !------------------------------- 2754 2755 elseif (icodcl(ifac,iut).eq.2) then 2756 2757 pimpv(1) = rcodcl(ifac,iut,1) 2758 cflv(1) = rcodcl(ifac,iut,2) 2759 pimpv(2) = rcodcl(ifac,ivt,1) 2760 cflv(2) = rcodcl(ifac,ivt,2) 2761 pimpv(3) = rcodcl(ifac,iwt,1) 2762 cflv(3) = rcodcl(ifac,iwt,2) 2763 2764 call set_convective_outlet_vector_aniso & 2765 ( coefaut(:,ifac) , cofafut(:,ifac) , & 2766 coefbut(:,:,ifac), cofbfut(:,:,ifac), & 2767 pimpv , cflv , hintt ) 2768 2769 ! Boundary conditions for thermal transport equation 2770 do isou = 1, 3 2771 cofarut(isou,ifac) = coefaut(isou,ifac) 2772 do jsou =1, 3 2773 cofbrut(isou,jsou,ifac) = coefbut(isou,jsou,ifac) 2774 enddo 2775 enddo 2776 2777 endif 2778 2779 endif 2780 2781 enddo 2782 2783 ! Vector transported quantity (dimension may be greater than 3) 2784 else 2785 2786 if (b_f_id .ge. 0) call field_get_val_v(b_f_id, bvar_v) 2787 2788 call field_get_coefa_v(ivarfl(ivar), coefav) 2789 call field_get_coefb_v(ivarfl(ivar), coefbv) 2790 call field_get_coefaf_v(ivarfl(ivar), cofafv) 2791 call field_get_coefbf_v(ivarfl(ivar), cofbfv) 2792 2793 do ifac = 1, nfabor 2794 2795 iel = ifabor(ifac) 2796 2797 ! --- Physical Properties 2798 visctc = visct(iel) 2799 2800 ! --- Geometrical quantities 2801 distbf = distb(ifac) 2802 2803 ! --- Prise en compte de Cp ou CV 2804 ! (dans le Cas compressible ihcp=0) 2805 2806 cpp = 1.d0 2807 if (ihcp.eq.0) then 2808 cpp = 1.d0 2809 elseif (ihcp.eq.2) then 2810 cpp = cpro_cp(iel) 2811 elseif (ihcp.eq.1) then 2812 cpp = cp0 2813 endif 2814 2815 ! --- Viscosite variable ou non 2816 if (ifcvsl.lt.0) then 2817 rkl = visls_0 2818 else 2819 rkl = viscls(iel) 2820 endif 2821 2822 ! Scalar diffusivity 2823 if (iand(vcopt%idften, ISOTROPIC_DIFFUSION).ne.0) then 2824 if (ippmod(idarcy).eq.-1) then !FIXME 2825 hint = (rkl+vcopt%idifft*cpp*visctc/turb_schmidt)/distbf 2826 else ! idarcy = 1 2827 hint = rkl/distbf 2828 endif 2829 2830 hintt(1) = hint 2831 hintt(2) = hint 2832 hintt(3) = hint 2833 hintt(4) = 0.d0 2834 hintt(5) = 0.d0 2835 hintt(6) = 0.d0 2836 2837 ! Symmetric tensor diffusivity 2838 elseif (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 2839 2840 temp = vcopt%idifft*cpp*ctheta(ii)/csrij 2841 hintt(1) = (rkl + temp*visten(1,iel))/distbf 2842 hintt(2) = (rkl + temp*visten(2,iel))/distbf 2843 hintt(3) = (rkl + temp*visten(3,iel))/distbf 2844 hintt(4) = temp*visten(4,iel) /distbf 2845 hintt(5) = temp*visten(5,iel) /distbf 2846 hintt(6) = temp*visten(6,iel) /distbf 2847 2848 endif 2849 2850 ! Dirichlet Boundary Condition 2851 !----------------------------- 2852 2853 if (icodcl(ifac,ivar).eq.1) then 2854 2855 pimpv(1) = rcodcl(ifac,ivar ,1) 2856 pimpv(2) = rcodcl(ifac,ivar+1,1) 2857 pimpv(3) = rcodcl(ifac,ivar+2,1) 2858 hextv(1) = rcodcl(ifac,ivar ,2) 2859 hextv(2) = rcodcl(ifac,ivar+1,2) 2860 hextv(3) = rcodcl(ifac,ivar+2,2) 2861 2862 call set_dirichlet_vector_aniso & 2863 ( coefav(:,ifac), cofafv(:,ifac), & 2864 coefbv(:,:,ifac), cofbfv(:,:,ifac), & 2865 pimpv , hintt , hextv) 2866 2867 ! Store boundary value 2868 if (b_f_id .ge. 0) then 2869 ! B_ij. Pj(I) 2870 do isou = 1, 3 2871 b_pvari(isou) = 0.d0 2872 do jsou = 1, 3 2873 b_pvari(isou) = b_pvari(isou) & 2874 + coefbv(jsou, isou, ifac) * bvar_v(jsou,ifac) 2875 enddo 2876 enddo 2877 do isou = 1, 3 2878 bvar_v(isou, ifac) = coefav(isou, ifac) + b_pvari(isou) 2879 enddo 2880 endif 2881 endif 2882 2883 ! Neumann Boundary Conditions 2884 !---------------------------- 2885 2886 if (icodcl(ifac,ivar).eq.3) then 2887 2888 qimpv(1) = rcodcl(ifac,ivar ,3) 2889 qimpv(2) = rcodcl(ifac,ivar+1,3) 2890 qimpv(3) = rcodcl(ifac,ivar+2,3) 2891 2892 call set_neumann_vector_aniso & 2893 ( coefav(:,ifac), cofafv(:,ifac), & 2894 coefbv(:,:,ifac), cofbfv(:,:,ifac), & 2895 qimpv , hintt ) 2896 2897 ! Store boundary value 2898 if (b_f_id .ge. 0) then 2899 ! B_ij. Pj(I) 2900 do isou = 1, 3 2901 b_pvari(isou) = 0.d0 2902 do jsou = 1, 3 2903 b_pvari(isou) = b_pvari(isou) & 2904 + coefbv(jsou, isou, ifac) * bvar_v(jsou,ifac) 2905 enddo 2906 enddo 2907 do isou = 1, 3 2908 bvar_v(isou, ifac) = coefav(isou, ifac) + b_pvari(isou) 2909 enddo 2910 endif 2911 2912 ! Convective Boundary Conditions 2913 !------------------------------- 2914 2915 elseif (icodcl(ifac,ivar).eq.2) then 2916 2917 pimpv(1) = rcodcl(ifac,ivar ,1) 2918 cflv(1) = rcodcl(ifac,ivar ,2) 2919 pimpv(2) = rcodcl(ifac,ivar+1,1) 2920 cflv(2) = rcodcl(ifac,ivar+1,2) 2921 pimpv(3) = rcodcl(ifac,ivar+2,1) 2922 cflv(3) = rcodcl(ifac,ivar+2,2) 2923 2924 call set_convective_outlet_vector_aniso & 2925 ( coefav(:,ifac), cofafv(:,ifac), & 2926 coefbv(:,:,ifac), cofbfv(:,:,ifac), & 2927 pimpv , cflv , hintt) 2928 2929 ! Store boundary value 2930 if (b_f_id .ge. 0) then 2931 ! B_ij. Pj(I) 2932 do isou = 1, 3 2933 b_pvari(isou) = 0.d0 2934 do jsou = 1, 3 2935 b_pvari(isou) = b_pvari(isou) & 2936 + coefbv(jsou, isou, ifac) * bvar_v(jsou,ifac) 2937 enddo 2938 enddo 2939 do isou = 1, 3 2940 bvar_v(isou, ifac) = coefav(isou, ifac) + b_pvari(isou) 2941 enddo 2942 endif 2943 2944 ! Imposed value for the convection operator, imposed flux for diffusion 2945 !---------------------------------------------------------------------- 2946 2947 elseif (icodcl(ifac,ivar).eq.13) then 2948 2949 pimpv = rcodcl(ifac,ivar ,1) 2950 qimpv = rcodcl(ifac,ivar ,3) 2951 pimpv = rcodcl(ifac,ivar+1,1) 2952 qimpv = rcodcl(ifac,ivar+1,3) 2953 pimpv = rcodcl(ifac,ivar+2,1) 2954 qimpv = rcodcl(ifac,ivar+2,3) 2955 2956 call set_dirichlet_conv_neumann_diff_vector & 2957 ( coefav(:,ifac) , cofafv(:,ifac) , & 2958 coefbv(:,:,ifac), cofbfv(:,:,ifac), & 2959 pimpv , qimpv ) 2960 2961 ! Store boundary value 2962 if (b_f_id .ge. 0) then 2963 ! B_ij. Pj(I) 2964 do isou = 1, 3 2965 b_pvari(isou) = 0.d0 2966 do jsou = 1, 3 2967 b_pvari(isou) = b_pvari(isou) & 2968 + coefbv(jsou, isou, ifac) * bvar_v(jsou,ifac) 2969 enddo 2970 enddo 2971 do isou = 1, 3 2972 bvar_v(isou, ifac) = coefav(isou, ifac) + b_pvari(isou) 2973 enddo 2974 endif 2975 2976 ! convective boundary for Marangoni effects 2977 !(generalized symmetry condition) 2978 !------------------------------------------ 2979 2980 elseif (icodcl(ifac,ivar).eq.14) then 2981 2982 pimpv(1) = rcodcl(ifac,ivar ,1) 2983 pimpv(2) = rcodcl(ifac,ivar+1,1) 2984 pimpv(3) = rcodcl(ifac,ivar+2,1) 2985 2986 qimpv(1) = rcodcl(ifac,ivar ,3) 2987 qimpv(2) = rcodcl(ifac,ivar+1,3) 2988 qimpv(3) = rcodcl(ifac,ivar+2,3) 2989 2990 normal(1) = surfbo(1,ifac)/surfbn(ifac) 2991 normal(2) = surfbo(2,ifac)/surfbn(ifac) 2992 normal(3) = surfbo(3,ifac)/surfbn(ifac) 2993 2994 ! coupled solving of the velocity components 2995 2996 call set_generalized_sym_vector_aniso & 2997 ( coefav(:,ifac) , cofafv(:,ifac) , & 2998 coefbv(:,:,ifac), cofbfv(:,:,ifac), & 2999 pimpv , qimpv , hintt, normal ) 3000 3001 ! Store boundary value 3002 if (b_f_id .ge. 0) then 3003 ! B_ij. Pj(I) 3004 do isou = 1, 3 3005 b_pvari(isou) = 0.d0 3006 do jsou = 1, 3 3007 b_pvari(isou) = b_pvari(isou) & 3008 + coefbv(jsou, isou, ifac) * bvar_v(jsou,ifac) 3009 enddo 3010 enddo 3011 do isou = 1, 3 3012 bvar_v(isou, ifac) = coefav(isou, ifac) + b_pvari(isou) 3013 enddo 3014 endif 3015 3016 ! Neumann on the normal component, Dirichlet on tangential components 3017 !-------------------------------------------------------------------- 3018 3019 elseif (icodcl(ifac,ivar).eq.11) then 3020 3021 ! Dirichlet to impose on the tangential components 3022 pimpv(1) = rcodcl(ifac,ivar ,1) 3023 pimpv(2) = rcodcl(ifac,ivar+1,1) 3024 pimpv(3) = rcodcl(ifac,ivar+2,1) 3025 3026 ! Flux to impose on the normal component 3027 qimpv(1) = rcodcl(ifac,ivar ,3) 3028 qimpv(2) = rcodcl(ifac,ivar+1,3) 3029 qimpv(3) = rcodcl(ifac,ivar+2,3) 3030 3031 normal(1) = surfbo(1,ifac)/surfbn(ifac) 3032 normal(2) = surfbo(2,ifac)/surfbn(ifac) 3033 normal(3) = surfbo(3,ifac)/surfbn(ifac) 3034 3035 ! coupled solving of the velocity components 3036 3037 call set_generalized_dirichlet_vector_aniso & 3038 ( coefav(:,ifac) , cofafv(:,ifac) , & 3039 coefbv(:,:,ifac), cofbfv(:,:,ifac), & 3040 pimpv , qimpv , hintt, normal ) 3041 3042 ! Store boundary value 3043 if (b_f_id .ge. 0) then 3044 ! B_ij. Pj(I) 3045 do isou = 1, 3 3046 b_pvari(isou) = 0.d0 3047 do jsou = 1, 3 3048 b_pvari(isou) = b_pvari(isou) & 3049 + coefbv(jsou, isou, ifac) * bvar_v(jsou,ifac) 3050 enddo 3051 enddo 3052 do isou = 1, 3 3053 bvar_v(isou, ifac) = coefav(isou, ifac) + b_pvari(isou) 3054 enddo 3055 endif 3056 3057 endif 3058 3059 enddo ! End of loop on faces 3060 3061 endif ! End of vector transported quantities 3062 3063 ! EB-GGDH/AFM/DFM alpha boundary conditions 3064 if (turb_flux_model.eq.11 .or. turb_flux_model.eq.21 .or. turb_flux_model.eq.31) then 3065 3066 ! Name of the scalar ivar 3067 call field_get_name(ivarfl(ivar), fname) 3068 3069 ! Index of the corresponding turbulent flux 3070 call field_get_id(trim(fname)//'_alpha', f_id) 3071 3072 call field_get_key_int(f_id, keyvar, ialt) 3073 3074 call field_get_coefa_s (f_id, coefap) 3075 call field_get_coefb_s (f_id, coefbp) 3076 call field_get_coefaf_s(f_id, cofafp) 3077 call field_get_coefbf_s(f_id, cofbfp) 3078 3079 do ifac = 1, nfabor 3080 3081 iel = ifabor(ifac) 3082 3083 distbf = distb(ifac) 3084 3085 hint = 1.d0/distbf 3086 3087 ! Dirichlet Boundary Condition 3088 !----------------------------- 3089 3090 if (icodcl(ifac,ialt).eq.1) then 3091 3092 pimp = rcodcl(ifac,ialt,1) 3093 hext = rcodcl(ifac,ialt,2) 3094 3095 call set_dirichlet_scalar & 3096 ( coefap(ifac), cofafp(ifac), & 3097 coefbp(ifac), cofbfp(ifac), & 3098 pimp , hint , hext ) 3099 3100 endif 3101 3102 ! Neumann Boundary Conditions 3103 !---------------------------- 3104 3105 if (icodcl(ifac,ialt).eq.3) then 3106 3107 dimp = rcodcl(ifac,ialt,3) 3108 3109 call set_neumann_scalar & 3110 ( coefap(ifac), cofafp(ifac), & 3111 coefbp(ifac), cofbfp(ifac), & 3112 dimp , hint ) 3113 3114 ! Radiative Boundary Conditions 3115 !------------------------------- 3116 3117 elseif (icodcl(ifac,ialt).eq.2) then 3118 3119 pimp = rcodcl(ifac,ialt,1) 3120 cfl = rcodcl(ifac,ialt,2) 3121 3122 call set_convective_outlet_scalar & 3123 ( coefap(ifac), cofafp(ifac), & 3124 coefbp(ifac), cofbfp(ifac), & 3125 pimp , cfl , hint ) 3126 3127 ! Imposed value for the convection operator, 3128 ! imposed flux for diffusion 3129 !------------------------------------------- 3130 3131 elseif (icodcl(ifac,ialt).eq.13) then 3132 3133 pimp = rcodcl(ifac,ialt,1) 3134 dimp = rcodcl(ifac,ialt,3) 3135 3136 call set_dirichlet_conv_neumann_diff_scalar & 3137 ( coefap(ifac), cofafp(ifac), & 3138 coefbp(ifac), cofbfp(ifac), & 3139 pimp , dimp ) 3140 3141 3142 endif 3143 enddo! End of loop on face 3144 endif 3145 3146 enddo! End of loop on scalars 3147 3148endif 3149 3150!=============================================================================== 3151! 14. Mesh velocity (ALE module): Dirichlet and Neumann and convective outlet 3152!=============================================================================== 3153 3154if (iale.eq.1) then 3155 3156 call field_get_coefa_v(ivarfl(iuma), claale) 3157 call field_get_coefb_v(ivarfl(iuma), clbale) 3158 call field_get_coefaf_v(ivarfl(iuma), cfaale) 3159 call field_get_coefbf_v(ivarfl(iuma), cfbale) 3160 3161 call field_get_key_struct_var_cal_opt(ivarfl(iuma), vcopt) 3162 3163 if (iand(vcopt%idften, ISOTROPIC_DIFFUSION).ne.0) then 3164 call field_get_val_s(ivisma, cpro_visma_s) 3165 elseif (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 3166 call field_get_val_v(ivisma, cpro_visma_v) 3167 endif 3168 3169 do ifac = 1, nfabor 3170 3171 iel = ifabor(ifac) 3172 distbf = distb(ifac) 3173 3174 if (iand(vcopt%idften, ISOTROPIC_DIFFUSION).ne.0) then 3175 3176 hintt(1) = cpro_visma_s(iel)/distbf 3177 hintt(2) = cpro_visma_s(iel)/distbf 3178 hintt(3) = cpro_visma_s(iel)/distbf 3179 hintt(4) = 0.d0 3180 hintt(5) = 0.d0 3181 hintt(6) = 0.d0 3182 3183 elseif (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 3184 3185 hintt(1) = cpro_visma_v(1,iel)/distbf 3186 hintt(2) = cpro_visma_v(2,iel)/distbf 3187 hintt(3) = cpro_visma_v(3,iel)/distbf 3188 hintt(4) = cpro_visma_v(4,iel)/distbf 3189 hintt(5) = cpro_visma_v(5,iel)/distbf 3190 hintt(6) = cpro_visma_v(6,iel)/distbf 3191 3192 endif 3193 3194 ! Dirichlet Boundary Conditions 3195 !------------------------------ 3196 3197 if (icodcl(ifac,iuma).eq.1) then 3198 3199 pimpv(1) = rcodcl(ifac,iuma,1) 3200 pimpv(2) = rcodcl(ifac,ivma,1) 3201 pimpv(3) = rcodcl(ifac,iwma,1) 3202 hextv(1) = rcodcl(ifac,iuma,2) 3203 hextv(2) = rcodcl(ifac,ivma,2) 3204 hextv(3) = rcodcl(ifac,iwma,2) 3205 3206 call set_dirichlet_vector_aniso & 3207 ( claale(:,ifac) , cfaale(:,ifac) , & 3208 clbale(:,:,ifac), cfbale(:,:,ifac), & 3209 pimpv , hintt , hextv ) 3210 3211 ! Neumann Boundary Conditions 3212 !---------------------------- 3213 3214 elseif (icodcl(ifac,iuma).eq.3) then 3215 3216 ! Coupled solving of the velocity components 3217 3218 qimpv(1) = rcodcl(ifac,iuma,3) 3219 qimpv(2) = rcodcl(ifac,ivma,3) 3220 qimpv(3) = rcodcl(ifac,iwma,3) 3221 3222 call set_neumann_vector_aniso & 3223 ( claale(:,ifac) , cfaale(:,ifac) , & 3224 clbale(:,:,ifac), cfbale(:,:,ifac), & 3225 qimpv , hintt ) 3226 3227 3228 ! Convective Boundary Conditions 3229 !------------------------------- 3230 3231 elseif (icodcl(ifac,iuma).eq.2) then 3232 3233 ! Coupled solving of the velocity components 3234 3235 pimpv(1) = rcodcl(ifac,iuma,1) 3236 cflv(1) = rcodcl(ifac,iuma,2) 3237 pimpv(2) = rcodcl(ifac,ivma,1) 3238 cflv(2) = rcodcl(ifac,ivma,2) 3239 pimpv(3) = rcodcl(ifac,iwma,1) 3240 cflv(3) = rcodcl(ifac,iwma,2) 3241 3242 call set_convective_outlet_vector_aniso & 3243 ( claale(:,ifac) , cfaale(:,ifac) , & 3244 clbale(:,:,ifac), cfbale(:,:,ifac), & 3245 pimpv , cflv , hintt ) 3246 3247 endif 3248 3249 enddo 3250 3251endif 3252 3253!=============================================================================== 3254! 15. Compute stresses at boundary (step 1 over 5) 3255!=============================================================================== 3256 3257if (iforbr.ge.0 .and. iterns.eq.1) then 3258 3259 ! Coupled solving of the velocity components 3260 do ifac = 1, nfabor 3261 srfbnf = surfbn(ifac) 3262 3263 ! The implicit term is added after having updated the velocity 3264 forbr(1,ifac) = ( cofafu(1,ifac) & 3265 + cofbfu(1,1,ifac) * velipb(ifac,1) & 3266 + cofbfu(1,2,ifac) * velipb(ifac,2) & 3267 + cofbfu(1,3,ifac) * velipb(ifac,3) )*srfbnf 3268 forbr(2,ifac) = ( cofafu(2,ifac) & 3269 + cofbfu(2,1,ifac) * velipb(ifac,1) & 3270 + cofbfu(2,2,ifac) * velipb(ifac,2) & 3271 + cofbfu(2,3,ifac) * velipb(ifac,3) )*srfbnf 3272 forbr(3,ifac) = ( cofafu(3,ifac) & 3273 + cofbfu(3,1,ifac) * velipb(ifac,1) & 3274 + cofbfu(3,2,ifac) * velipb(ifac,2) & 3275 + cofbfu(3,3,ifac) * velipb(ifac,3) )*srfbnf 3276 enddo 3277endif 3278 3279! Free memory 3280deallocate(velipb) 3281if (associated(rijipb)) deallocate(rijipb) 3282 3283!=============================================================================== 3284! 16. Update of boundary temperature when saved and not a variable. 3285!=============================================================================== 3286 3287if (itherm.eq.2) then 3288 3289 if (itempb.ge.0) then 3290 3291 call field_get_val_s(itempb, btemp_s) 3292 3293 ! If we also have a boundary value field for the thermal 3294 ! scalar, copy its values first. 3295 3296 ! If we do not have a boundary value field for the thermal scalar, 3297 ! then boundary values for the thermal scalar were directly 3298 ! saved to the boundary temperature field, so no copy is needed. 3299 3300 f_id = ivarfl(isca(iscalt)) 3301 call field_get_key_int(f_id, kbfid, b_f_id) 3302 3303 if (b_f_id .ge. 0) then 3304 call field_get_val_s(b_f_id, bvar_s) 3305 else 3306 allocate(bvar_s(nfabor)) 3307 do ii = 1, nfabor 3308 bvar_s(ii) = btemp_s(ii) 3309 enddo 3310 endif 3311 3312 call cs_ht_convert_h_to_t_faces(bvar_s, btemp_s) 3313 3314 if (b_f_id .lt. 0) then 3315 deallocate(bvar_s) 3316 endif 3317 3318 ! In case of assigned temperature values, overwrite computed 3319 ! wall temperature with prescribed one to avoid issues due to 3320 ! enthalpy -> temperature conversion precision 3321 ! (T -> H -> T at the boundary does not preserve T) 3322 do ii = 1, nbt2h 3323 ifac = lbt2h(ii) + 1 3324 btemp_s(ifac) = vbt2h(ifac) 3325 enddo 3326 3327 endif 3328 3329 deallocate(lbt2h) 3330 deallocate(vbt2h) 3331 3332endif 3333 3334!=============================================================================== 3335! 17. Formats 3336!=============================================================================== 3337 3338 3010 format( & 3339 'Incoming flow detained for ', I10 , & 3340 ' outlet faces on ',I10) 3341 3342!---- 3343! End 3344!---- 3345 3346return 3347end subroutine condli 3348 3349!=============================================================================== 3350! Local functions 3351!=============================================================================== 3352 3353!------------------------------------------------------------------------------- 3354! Arguments 3355!______________________________________________________________________________. 3356! mode name role ! 3357!______________________________________________________________________________! 3358!> \param[out] coefa explicit BC coefficient for gradients 3359!> \param[out] cofaf explicit BC coefficient for diffusive flux 3360!> \param[out] coefb implicit BC coefficient for gradients 3361!> \param[out] cofbf implicit BC coefficient for diffusive flux 3362!> \param[in] pimp Dirichlet value to impose 3363!> \param[in] hint Internal exchange coefficient 3364!> \param[in] hext External exchange coefficient (10^30 by default) 3365!_______________________________________________________________________________ 3366 3367subroutine set_dirichlet_scalar & 3368 ( coefa , cofaf, coefb , cofbf, pimp , hint, hext) 3369 3370!=============================================================================== 3371! Module files 3372!=============================================================================== 3373 3374use cstnum, only: rinfin 3375 3376!=============================================================================== 3377 3378implicit none 3379 3380! Arguments 3381 3382double precision coefa, cofaf, coefb, cofbf, pimp, hint, hext 3383 3384! Local variables 3385 3386double precision heq 3387 3388!=============================================================================== 3389 3390if (abs(hext).gt.rinfin*0.5d0) then 3391 3392 ! Gradient BCs 3393 coefa = pimp 3394 coefb = 0.d0 3395 3396 ! Flux BCs 3397 cofaf = -hint*pimp 3398 cofbf = hint 3399 3400else 3401 3402 ! Gradient BCs 3403 coefa = hext*pimp/(hint + hext) 3404 coefb = hint /(hint + hext) 3405 3406 ! Flux BCs 3407 heq = hint*hext/(hint + hext) 3408 cofaf = -heq*pimp 3409 cofbf = heq 3410 3411endif 3412 3413return 3414end subroutine set_dirichlet_scalar 3415 3416!=============================================================================== 3417 3418!------------------------------------------------------------------------------- 3419! Arguments 3420!______________________________________________________________________________. 3421! mode name role ! 3422!______________________________________________________________________________! 3423!> \param[out] coefa explicit BC coefficient for gradients 3424!> \param[out] cofaf explicit BC coefficient for diffusive flux 3425!> \param[out] coefb implicit BC coefficient for gradients 3426!> \param[out] cofbf implicit BC coefficient for diffusive flux 3427!> \param[in] pimpv Dirichlet value to impose 3428!> \param[in] hint Internal exchange coefficient 3429!> \param[in] hextv External exchange coefficient (10^30 by default) 3430!_______________________________________________________________________________ 3431 3432subroutine set_dirichlet_vector & 3433 ( coefa , cofaf, coefb , cofbf, pimpv , hint , hextv) 3434 3435!=============================================================================== 3436! Module files 3437!=============================================================================== 3438 3439use cstnum, only: rinfin 3440 3441!=============================================================================== 3442 3443implicit none 3444 3445! Arguments 3446 3447double precision coefa(3), cofaf(3) 3448double precision coefb(3,3), cofbf(3,3) 3449double precision pimpv(3) 3450double precision hint 3451double precision hextv(3) 3452 3453! Local variables 3454 3455integer isou , jsou 3456double precision heq 3457 3458!=============================================================================== 3459 3460do isou = 1, 3 3461 3462 if (abs(hextv(isou)).gt.rinfin*0.5d0) then 3463 ! Gradient BCs 3464 coefa(isou) = pimpv(isou) 3465 do jsou = 1, 3 3466 coefb(isou,jsou) = 0.d0 3467 enddo 3468 3469 ! Flux BCs 3470 cofaf(isou) = -hint*pimpv(isou) 3471 do jsou = 1, 3 3472 if (jsou.eq.isou) then 3473 cofbf(isou,jsou) = hint 3474 else 3475 cofbf(isou,jsou) = 0.d0 3476 endif 3477 enddo 3478 3479 else 3480 3481 heq = hint*hextv(isou)/(hint + hextv(isou)) 3482 3483 ! Gradient BCs 3484 coefa(isou) = hextv(isou)*pimpv(isou)/(hint + hextv(isou)) 3485 do jsou = 1, 3 3486 if (jsou.eq.isou) then 3487 coefb(isou,jsou) = hint/(hint + hextv(isou)) 3488 else 3489 coefb(isou,jsou) = 0.d0 3490 endif 3491 enddo 3492 3493 ! Flux BCs 3494 cofaf(isou) = -heq*pimpv(isou) 3495 do jsou = 1, 3 3496 if (jsou.eq.isou) then 3497 cofbf(isou,jsou) = heq 3498 else 3499 cofbf(isou,jsou) = 0.d0 3500 endif 3501 enddo 3502 3503 endif 3504 3505enddo 3506 3507return 3508end subroutine set_dirichlet_vector 3509 3510!=============================================================================== 3511 3512!------------------------------------------------------------------------------- 3513! Arguments 3514!______________________________________________________________________________. 3515! mode name role ! 3516!______________________________________________________________________________! 3517!> \param[out] coefa explicit BC coefficient for gradients 3518!> \param[out] cofaf explicit BC coefficient for diffusive flux 3519!> \param[out] coefb implicit BC coefficient for gradients 3520!> \param[out] cofbf implicit BC coefficient for diffusive flux 3521!> \param[in] pimpts Dirichlet value to impose 3522!> \param[in] hint Internal exchange coefficient 3523!> \param[in] hextts External exchange coefficient (10^30 by default) 3524!_______________________________________________________________________________ 3525 3526subroutine set_dirichlet_tensor & 3527 ( coefa , cofaf, coefb , cofbf, pimpts , hint , hextts) 3528 3529!=============================================================================== 3530! Module files 3531!=============================================================================== 3532 3533use cstnum, only: rinfin 3534 3535!=============================================================================== 3536 3537implicit none 3538 3539! Arguments 3540 3541double precision coefa(6), cofaf(6) 3542double precision coefb(6,6), cofbf(6,6) 3543double precision pimpts(6) 3544double precision hint 3545double precision hextts(6) 3546 3547! Local variables 3548 3549integer isou , jsou 3550double precision heq 3551 3552!=============================================================================== 3553 3554do isou = 1, 6 3555 3556 if (abs(hextts(isou)).gt.rinfin*0.5d0) then 3557 ! Gradient BCs 3558 coefa(isou) = pimpts(isou) 3559 do jsou = 1, 6 3560 coefb(isou,jsou) = 0.d0 3561 enddo 3562 3563 ! Flux BCs 3564 cofaf(isou) = -hint*pimpts(isou) 3565 do jsou = 1, 6 3566 if (jsou.eq.isou) then 3567 cofbf(isou,jsou) = hint 3568 else 3569 cofbf(isou,jsou) = 0.d0 3570 endif 3571 enddo 3572 3573 else 3574 3575 heq = hint*hextts(isou)/(hint + hextts(isou)) 3576 3577 ! Gradient BCs 3578 coefa(isou) = hextts(isou)*pimpts(isou)/(hint + hextts(isou)) 3579 do jsou = 1, 6 3580 if (jsou.eq.isou) then 3581 coefb(isou,jsou) = hint/(hint + hextts(isou)) 3582 else 3583 coefb(isou,jsou) = 0.d0 3584 endif 3585 enddo 3586 3587 ! Flux BCs 3588 cofaf(isou) = -heq*pimpts(isou) 3589 do jsou = 1, 6 3590 if (jsou.eq.isou) then 3591 cofbf(isou,jsou) = heq 3592 else 3593 cofbf(isou,jsou) = 0.d0 3594 endif 3595 enddo 3596 3597 endif 3598 3599enddo 3600 3601return 3602end subroutine set_dirichlet_tensor 3603 3604!=============================================================================== 3605 3606!------------------------------------------------------------------------------- 3607! Arguments 3608!______________________________________________________________________________. 3609! mode name role ! 3610!______________________________________________________________________________! 3611!> \param[out] coefa explicit BC coefficient for gradients 3612!> \param[out] cofaf explicit BC coefficient for diffusive flux 3613!> \param[out] coefb implicit BC coefficient for gradients 3614!> \param[out] cofbf implicit BC coefficient for diffusive flux 3615!> \param[in] pimpv Dirichlet value to impose 3616!> \param[in] hint Internal exchange coefficient 3617!> \param[in] hextv External exchange coefficient (10^30 by default) 3618!_______________________________________________________________________________ 3619 3620subroutine set_dirichlet_vector_aniso & 3621 ( coefa , cofaf, coefb , cofbf, pimpv , hint , hextv) 3622 3623!=============================================================================== 3624! Module files 3625!=============================================================================== 3626 3627use cstnum, only: rinfin 3628 3629!=============================================================================== 3630 3631implicit none 3632 3633! Arguments 3634 3635double precision coefa(3), cofaf(3) 3636double precision coefb(3,3), cofbf(3,3) 3637double precision pimpv(3) 3638double precision hint(6) 3639double precision hextv(3) 3640 3641! Local variables 3642 3643integer isou , jsou 3644 3645!=============================================================================== 3646 3647do isou = 1, 3 3648 3649 if (abs(hextv(isou)).gt.rinfin*0.5d0) then 3650 ! Gradient BCs 3651 coefa(isou) = pimpv(isou) 3652 do jsou = 1, 3 3653 coefb(isou,jsou) = 0.d0 3654 enddo 3655 3656 else 3657 3658 call csexit(1) 3659 3660 endif 3661 3662enddo 3663 3664! Flux BCs 3665cofaf(1) = -(hint(1)*pimpv(1) + hint(4)*pimpv(2) + hint(6)*pimpv(3)) 3666cofaf(2) = -(hint(4)*pimpv(1) + hint(2)*pimpv(2) + hint(5)*pimpv(3)) 3667cofaf(3) = -(hint(6)*pimpv(1) + hint(5)*pimpv(2) + hint(3)*pimpv(3)) 3668cofbf(1,1) = hint(1) 3669cofbf(2,2) = hint(2) 3670cofbf(3,3) = hint(3) 3671cofbf(1,2) = hint(4) 3672cofbf(2,1) = hint(4) 3673cofbf(2,3) = hint(5) 3674cofbf(3,2) = hint(5) 3675cofbf(1,3) = hint(6) 3676cofbf(3,1) = hint(6) 3677 3678return 3679end subroutine set_dirichlet_vector_aniso 3680 3681!=============================================================================== 3682 3683!------------------------------------------------------------------------------- 3684! Arguments 3685!______________________________________________________________________________. 3686! mode name role ! 3687!______________________________________________________________________________! 3688!> \param[out] coefa explicit BC coefficient for gradients 3689!> \param[out] cofaf explicit BC coefficient for diffusive flux 3690!> \param[out] coefb implicit BC coefficient for gradients 3691!> \param[out] cofbf implicit BC coefficient for diffusive flux 3692!> \param[in] dimp Flux value to impose 3693!> \param[in] hint Internal exchange coefficient 3694!_______________________________________________________________________________ 3695 3696subroutine set_neumann_scalar & 3697 ( coefa , cofaf, coefb , cofbf, dimp , hint) 3698 3699!=============================================================================== 3700! Module files 3701!=============================================================================== 3702 3703!=============================================================================== 3704 3705implicit none 3706 3707! Arguments 3708 3709double precision coefa, cofaf, coefb, cofbf, dimp, hint 3710 3711! Local variables 3712 3713!=============================================================================== 3714 3715! Gradient BCs 3716coefa = -dimp/max(hint, 1.d-300) 3717coefb = 1.d0 3718 3719! Flux BCs 3720cofaf = dimp 3721cofbf = 0.d0 3722 3723return 3724end subroutine set_neumann_scalar 3725 3726!=============================================================================== 3727 3728!------------------------------------------------------------------------------- 3729! Arguments 3730!______________________________________________________________________________. 3731! mode name role ! 3732!______________________________________________________________________________! 3733!> \param[out] coefa explicit BC coefficient for gradients 3734!> \param[out] cofaf explicit BC coefficient for diffusive flux 3735!> \param[out] coefb implicit BC coefficient for gradients 3736!> \param[out] cofbf implicit BC coefficient for diffusive flux 3737!> \param[in] qimpv Flux value to impose 3738!> \param[in] hint Internal exchange coefficient 3739!_______________________________________________________________________________ 3740 3741subroutine set_neumann_vector & 3742 ( coefa , cofaf, coefb , cofbf, qimpv , hint) 3743 3744!=============================================================================== 3745! Module files 3746!=============================================================================== 3747 3748!=============================================================================== 3749 3750implicit none 3751 3752! Arguments 3753 3754double precision coefa(3), cofaf(3) 3755double precision coefb(3,3), cofbf(3,3) 3756double precision qimpv(3) 3757double precision hint 3758 3759! Local variables 3760 3761integer isou , jsou 3762 3763!=============================================================================== 3764 3765do isou = 1, 3 3766 3767 ! Gradient BCs 3768 coefa(isou) = -qimpv(isou)/max(hint, 1.d-300) 3769 do jsou = 1, 3 3770 if (jsou.eq.isou) then 3771 coefb(isou,jsou) = 1.d0 3772 else 3773 coefb(isou,jsou) = 0.d0 3774 endif 3775 enddo 3776 3777 ! Flux BCs 3778 cofaf(isou) = qimpv(isou) 3779 do jsou = 1, 3 3780 cofbf(isou,jsou) = 0.d0 3781 enddo 3782 3783enddo 3784 3785return 3786end subroutine set_neumann_vector 3787 3788!=============================================================================== 3789 3790!------------------------------------------------------------------------------- 3791! Arguments 3792!______________________________________________________________________________. 3793! mode name role ! 3794!______________________________________________________________________________! 3795!> \param[out] coefa explicit BC coefficient for gradients 3796!> \param[out] cofaf explicit BC coefficient for diffusive flux 3797!> \param[out] coefb implicit BC coefficient for gradients 3798!> \param[out] cofbf implicit BC coefficient for diffusive flux 3799!> \param[in] qimpts Flux value to impose 3800!> \param[in] hint Internal exchange coefficient 3801!_______________________________________________________________________________ 3802 3803subroutine set_neumann_tensor & 3804 ( coefa , cofaf, coefb , cofbf, qimpts , hint) 3805 3806!=============================================================================== 3807! Module files 3808!=============================================================================== 3809 3810!=============================================================================== 3811 3812implicit none 3813 3814! Arguments 3815 3816double precision coefa(6), cofaf(6) 3817double precision coefb(6,6), cofbf(6,6) 3818double precision qimpts(6) 3819double precision hint 3820 3821! Local variables 3822 3823integer isou , jsou 3824 3825!=============================================================================== 3826 3827do isou = 1, 6 3828 3829 ! Gradient BCs 3830 coefa(isou) = -qimpts(isou)/max(hint, 1.d-300) 3831 do jsou = 1, 6 3832 if (jsou.eq.isou) then 3833 coefb(isou,jsou) = 1.d0 3834 else 3835 coefb(isou,jsou) = 0.d0 3836 endif 3837 enddo 3838 3839 ! Flux BCs 3840 cofaf(isou) = qimpts(isou) 3841 do jsou = 1, 6 3842 cofbf(isou,jsou) = 0.d0 3843 enddo 3844 3845enddo 3846 3847return 3848end subroutine set_neumann_tensor 3849 3850!=============================================================================== 3851 3852!------------------------------------------------------------------------------- 3853! Arguments 3854!______________________________________________________________________________. 3855! mode name role ! 3856!______________________________________________________________________________! 3857!> \param[out] coefa explicit BC coefficient for gradients 3858!> \param[out] cofaf explicit BC coefficient for diffusive flux 3859!> \param[out] coefb implicit BC coefficient for gradients 3860!> \param[out] cofbf implicit BC coefficient for diffusive flux 3861!> \param[in] qimpv Flux value to impose 3862!> \param[in] hint Internal exchange coefficient 3863!_______________________________________________________________________________ 3864 3865subroutine set_neumann_vector_aniso & 3866 ( coefa , cofaf, coefb , cofbf, qimpv , hint) 3867 3868!=============================================================================== 3869! Module files 3870!=============================================================================== 3871 3872!=============================================================================== 3873 3874implicit none 3875 3876! Arguments 3877 3878double precision coefa(3), cofaf(3) 3879double precision coefb(3,3), cofbf(3,3) 3880double precision qimpv(3) 3881double precision hint(6) 3882 3883! Local variables 3884 3885integer isou , jsou 3886double precision invh(6), invdet, m(6) 3887 3888!=============================================================================== 3889m(1) = hint(2)*hint(3) - hint(5)*hint(5) 3890m(2) = hint(1)*hint(3) - hint(6)*hint(6) 3891m(3) = hint(1)*hint(2) - hint(4)*hint(4) 3892m(4) = hint(5)*hint(6) - hint(4)*hint(3) 3893m(5) = hint(4)*hint(6) - hint(1)*hint(5) 3894m(6) = hint(4)*hint(5) - hint(2)*hint(6) 3895 3896invdet = 1.d0/(hint(1)*m(1) + hint(4)*m(4) + hint(6)*m(6)) 3897 3898invh(1) = m(1) * invdet 3899invh(2) = m(2) * invdet 3900invh(3) = m(3) * invdet 3901invh(4) = m(4) * invdet 3902invh(5) = m(5) * invdet 3903invh(6) = m(6) * invdet 3904 3905! Gradient BCs 3906coefa(1) = -(invh(1)*qimpv(1) + invh(4)*qimpv(2) + invh(6)*qimpv(3)) 3907coefa(2) = -(invh(4)*qimpv(1) + invh(2)*qimpv(2) + invh(5)*qimpv(3)) 3908coefa(3) = -(invh(6)*qimpv(1) + invh(5)*qimpv(2) + invh(3)*qimpv(3)) 3909coefb(1,1) = 1.d0 3910coefb(2,2) = 1.d0 3911coefb(3,3) = 1.d0 3912coefb(1,2) = 0.d0 3913coefb(2,1) = 0.d0 3914coefb(2,3) = 0.d0 3915coefb(3,2) = 0.d0 3916coefb(1,3) = 0.d0 3917coefb(3,1) = 0.d0 3918 3919do isou = 1, 3 3920 3921 ! Flux BCs 3922 cofaf(isou) = qimpv(isou) 3923 do jsou = 1, 3 3924 cofbf(isou,jsou) = 0.d0 3925 enddo 3926 3927enddo 3928 3929return 3930end subroutine set_neumann_vector_aniso 3931 3932!=============================================================================== 3933 3934!------------------------------------------------------------------------------- 3935! Arguments 3936!______________________________________________________________________________. 3937! mode name role ! 3938!______________________________________________________________________________! 3939!> \param[out] coefa explicit BC coefficient for gradients 3940!> \param[out] cofaf explicit BC coefficient for diffusive flux 3941!> \param[out] coefb implicit BC coefficient for gradients 3942!> \param[out] cofbf implicit BC coefficient for diffusive flux 3943!> \param[in] pimpv Dirichlet value to impose on the normal 3944!> component 3945!> \param[in] qimpv Flux value to impose on the 3946!> tangential components 3947!> \param[in] hint Internal exchange coefficient 3948!> \param[in] normal normal 3949!_______________________________________________________________________________ 3950 3951subroutine set_generalized_sym_vector & 3952 ( coefa , cofaf, coefb , cofbf, pimpv, qimpv, hint, normal) 3953 3954!=============================================================================== 3955! Module files 3956!=============================================================================== 3957 3958!=============================================================================== 3959 3960implicit none 3961 3962! Arguments 3963 3964double precision coefa(3), cofaf(3) 3965double precision coefb(3,3), cofbf(3,3) 3966double precision hint 3967double precision normal(3) 3968double precision pimpv(3), qimpv(3) 3969 3970! Local variables 3971 3972integer isou , jsou 3973 3974!=============================================================================== 3975 3976do isou = 1, 3 3977 3978 ! Gradient BCs 3979 coefa(isou) = - qimpv(isou)/max(hint, 1.d-300) 3980 ! "[1 -n(x)n] Qimp / hint" is divided into two 3981 do jsou = 1, 3 3982 coefa(isou) = coefa(isou) & 3983 + normal(isou)*normal(jsou)*(pimpv(jsou)+qimpv(jsou)/max(hint, 1.d-300)) 3984 if (jsou.eq.isou) then 3985 coefb(isou,jsou) = 1.d0 - normal(isou)*normal(jsou) 3986 else 3987 coefb(isou,jsou) = - normal(isou)*normal(jsou) 3988 endif 3989 enddo 3990 3991 ! Flux BCs 3992 cofaf(isou) = qimpv(isou) 3993 ! "[1 -n(x)n] Qimp" is divided into two 3994 do jsou = 1, 3 3995 cofaf(isou) = cofaf(isou) & 3996 - normal(isou)*normal(jsou)*(hint*pimpv(jsou)+qimpv(jsou)) 3997 cofbf(isou,jsou) = hint*normal(isou)*normal(jsou) 3998 enddo 3999 4000enddo 4001 4002return 4003end subroutine set_generalized_sym_vector 4004 4005!=============================================================================== 4006 4007!------------------------------------------------------------------------------- 4008! Arguments 4009!______________________________________________________________________________. 4010! mode name role ! 4011!______________________________________________________________________________! 4012!> \param[out] coefa explicit BC coefficient for gradients 4013!> \param[out] cofaf explicit BC coefficient for diffusive flux 4014!> \param[out] coefb implicit BC coefficient for gradients 4015!> \param[out] cofbf implicit BC coefficient for diffusive flux 4016!> \param[in] pimpv Dirichlet value to impose on the normal 4017!> component 4018!> \param[in] qimpv Flux value to impose on the 4019!> tangential components 4020!> \param[in] hint Internal exchange coefficient 4021!> \param[in] normal normal 4022!_______________________________________________________________________________ 4023 4024subroutine set_generalized_sym_vector_aniso & 4025 ( coefa , cofaf, coefb , cofbf, pimpv, qimpv, hint, normal) 4026 4027!=============================================================================== 4028! Module files 4029!=============================================================================== 4030 4031!=============================================================================== 4032 4033implicit none 4034 4035! Arguments 4036 4037double precision coefa(3), cofaf(3) 4038double precision coefb(3,3), cofbf(3,3) 4039double precision hint(6) 4040double precision normal(3) 4041double precision pimpv(3), qimpv(3) 4042 4043! Local variables 4044 4045integer isou , jsou 4046double precision invh(6), invdet, m(6), qshint(3), hintpv(3), hintnm(3) 4047 4048!=============================================================================== 4049 4050m(1) = hint(2)*hint(3) - hint(5)*hint(5) 4051m(2) = hint(1)*hint(3) - hint(6)*hint(6) 4052m(3) = hint(1)*hint(2) - hint(4)*hint(4) 4053m(4) = hint(5)*hint(6) - hint(4)*hint(3) 4054m(5) = hint(4)*hint(6) - hint(1)*hint(5) 4055m(6) = hint(4)*hint(5) - hint(2)*hint(6) 4056 4057invdet = 1.d0/(hint(1)*m(1) + hint(4)*m(4) + hint(6)*m(6)) 4058 4059invh(1) = m(1) * invdet 4060invh(2) = m(2) * invdet 4061invh(3) = m(3) * invdet 4062invh(4) = m(4) * invdet 4063invh(5) = m(5) * invdet 4064invh(6) = m(6) * invdet 4065 4066qshint(1) = invh(1)*qimpv(1) + invh(4)*qimpv(2) + invh(6)*qimpv(3) 4067qshint(2) = invh(4)*qimpv(1) + invh(2)*qimpv(2) + invh(5)*qimpv(3) 4068qshint(3) = invh(6)*qimpv(1) + invh(5)*qimpv(2) + invh(3)*qimpv(3) 4069 4070hintpv(1) = hint(1)*pimpv(1) + hint(4)*pimpv(2) + hint(6)*pimpv(3) 4071hintpv(2) = hint(4)*pimpv(1) + hint(2)*pimpv(2) + hint(5)*pimpv(3) 4072hintpv(3) = hint(6)*pimpv(1) + hint(5)*pimpv(2) + hint(3)*pimpv(3) 4073 4074hintnm(1) = hint(1)*normal(1) + hint(4)*normal(2) + hint(6)*normal(3) 4075hintnm(2) = hint(4)*normal(1) + hint(2)*normal(2) + hint(5)*normal(3) 4076hintnm(3) = hint(6)*normal(1) + hint(5)*normal(2) + hint(3)*normal(3) 4077 4078do isou = 1, 3 4079 4080 ! Gradient BCs 4081 coefa(isou) = - qshint(isou) 4082 ! "[1 -n(x)n] Qimp / hint" is divided into two 4083 do jsou = 1, 3 4084 coefa(isou) = coefa(isou) & 4085 + normal(isou)*normal(jsou)*(pimpv(jsou)+qshint(jsou)) 4086 if (jsou.eq.isou) then 4087 coefb(isou,jsou) = 1.d0 - normal(isou)*normal(jsou) 4088 else 4089 coefb(isou,jsou) = - normal(isou)*normal(jsou) 4090 endif 4091 enddo 4092 4093 ! Flux BCs 4094 cofaf(isou) = qimpv(isou) 4095 ! "[1 -n(x)n] Qimp" is divided into two 4096 do jsou = 1, 3 4097 cofaf(isou) = cofaf(isou) & 4098 - normal(isou)*normal(jsou)*(hintpv(jsou)+qimpv(jsou)) 4099 cofbf(isou,jsou) = hintnm(isou)*normal(jsou) 4100 enddo 4101 4102enddo 4103 4104return 4105end subroutine set_generalized_sym_vector_aniso 4106 4107!=============================================================================== 4108 4109!------------------------------------------------------------------------------- 4110! Arguments 4111!______________________________________________________________________________. 4112! mode name role ! 4113!______________________________________________________________________________! 4114!> \param[out] coefa explicit BC coefficient for gradients 4115!> \param[out] cofaf explicit BC coefficient for diffusive flux 4116!> \param[out] coefb implicit BC coefficient for gradients 4117!> \param[out] cofbf implicit BC coefficient for diffusive flux 4118!> \param[in] pimpv Dirichlet value to impose on the tangential 4119!> components 4120!> \param[in] qimpv Flux value to impose on the 4121!> normal component 4122!> \param[in] hint Internal exchange coefficient 4123!> \param[in] normal normal 4124!_______________________________________________________________________________ 4125 4126subroutine set_generalized_dirichlet_vector & 4127 ( coefa , cofaf, coefb , cofbf, pimpv, qimpv, hint, normal) 4128 4129!=============================================================================== 4130! Module files 4131!=============================================================================== 4132 4133!=============================================================================== 4134 4135implicit none 4136 4137! Arguments 4138 4139double precision coefa(3), cofaf(3) 4140double precision coefb(3,3), cofbf(3,3) 4141double precision hint 4142double precision normal(3) 4143double precision pimpv(3), qimpv(3) 4144 4145! Local variables 4146 4147integer isou , jsou 4148 4149!=============================================================================== 4150 4151do isou = 1, 3 4152 4153 ! Gradient BCs 4154 ! "[1 -n(x)n] Pimp" is divided into two 4155 coefa(isou) = pimpv(isou) 4156 do jsou = 1, 3 4157 coefa(isou) = coefa(isou) & 4158 - normal(isou)*normal(jsou)*(pimpv(jsou)+qimpv(jsou)/max(hint, 1.d-300)) 4159 coefb(isou,jsou) = normal(isou)*normal(jsou) 4160 enddo 4161 4162 ! Flux BCs 4163 ! "[1 -n(x)n] Pimp" is divided into two 4164 cofaf(isou) = -hint*pimpv(isou) 4165 do jsou = 1, 3 4166 cofaf(isou) = cofaf(isou) & 4167 + normal(isou)*normal(jsou)*(qimpv(jsou)+pimpv(jsou)*hint) 4168 if (jsou.eq.isou) then 4169 cofbf(isou,jsou) = hint*(1.d0-normal(isou)*normal(jsou)) 4170 else 4171 cofbf(isou,jsou) = -hint*normal(isou)*normal(jsou) 4172 endif 4173 enddo 4174 4175enddo 4176 4177return 4178end subroutine set_generalized_dirichlet_vector 4179 4180!=============================================================================== 4181 4182!------------------------------------------------------------------------------- 4183! Arguments 4184!______________________________________________________________________________. 4185! mode name role ! 4186!______________________________________________________________________________! 4187!> \param[out] coefa explicit BC coefficient for gradients 4188!> \param[out] cofaf explicit BC coefficient for diffusive flux 4189!> \param[out] coefb implicit BC coefficient for gradients 4190!> \param[out] cofbf implicit BC coefficient for diffusive flux 4191!> \param[in] pimpv Dirichlet value to impose on the tangential 4192!> components 4193!> \param[in] qimpv Flux value to impose on the 4194!> normal component 4195!> \param[in] hint Internal exchange coefficient 4196!> \param[in] normal normal 4197!_______________________________________________________________________________ 4198 4199subroutine set_generalized_dirichlet_vector_aniso & 4200 ( coefa , cofaf, coefb , cofbf, pimpv, qimpv, hint, normal) 4201 4202!=============================================================================== 4203! Module files 4204!=============================================================================== 4205 4206!=============================================================================== 4207 4208implicit none 4209 4210! Arguments 4211 4212double precision coefa(3), cofaf(3) 4213double precision coefb(3,3), cofbf(3,3) 4214double precision hint(6) 4215double precision normal(3) 4216double precision pimpv(3), qimpv(3) 4217 4218! Local variables 4219 4220integer isou , jsou 4221double precision invh(6), invdet, m(6), qshint(3), hintpv(3), hintnm(3) 4222 4223!=============================================================================== 4224 4225m(1) = hint(2)*hint(3) - hint(5)*hint(5) 4226m(2) = hint(1)*hint(3) - hint(6)*hint(6) 4227m(3) = hint(1)*hint(2) - hint(4)*hint(4) 4228m(4) = hint(5)*hint(6) - hint(4)*hint(3) 4229m(5) = hint(4)*hint(6) - hint(1)*hint(5) 4230m(6) = hint(4)*hint(5) - hint(2)*hint(6) 4231 4232invdet = 1.d0/(hint(1)*m(1) + hint(4)*m(4) + hint(6)*m(6)) 4233 4234invh(1) = m(1) * invdet 4235invh(2) = m(2) * invdet 4236invh(3) = m(3) * invdet 4237invh(4) = m(4) * invdet 4238invh(5) = m(5) * invdet 4239invh(6) = m(6) * invdet 4240 4241qshint(1) = invh(1)*qimpv(1) + invh(4)*qimpv(2) + invh(6)*qimpv(3) 4242qshint(2) = invh(4)*qimpv(1) + invh(2)*qimpv(2) + invh(5)*qimpv(3) 4243qshint(3) = invh(6)*qimpv(1) + invh(5)*qimpv(2) + invh(3)*qimpv(3) 4244 4245hintpv(1) = hint(1)*pimpv(1) + hint(4)*pimpv(2) + hint(6)*pimpv(3) 4246hintpv(2) = hint(4)*pimpv(1) + hint(2)*pimpv(2) + hint(5)*pimpv(3) 4247hintpv(3) = hint(6)*pimpv(1) + hint(5)*pimpv(2) + hint(3)*pimpv(3) 4248 4249hintnm(1) = hint(1)*normal(1) + hint(4)*normal(2) + hint(6)*normal(3) 4250hintnm(2) = hint(4)*normal(1) + hint(2)*normal(2) + hint(5)*normal(3) 4251hintnm(3) = hint(6)*normal(1) + hint(5)*normal(2) + hint(3)*normal(3) 4252 4253do isou = 1, 3 4254 4255 ! Gradient BCs 4256 ! "[1 -n(x)n] Pimp" is divided into two 4257 coefa(isou) = pimpv(isou) 4258 do jsou = 1, 3 4259 coefa(isou) = coefa(isou) & 4260 - normal(isou)*normal(jsou)*(pimpv(jsou)+qshint(jsou)) 4261 coefb(isou,jsou) = normal(isou)*normal(jsou) 4262 enddo 4263 4264 ! Flux BCs 4265 ! "[1 -n(x)n] Pimp" is divided into two 4266 cofaf(isou) = -hintpv(isou) 4267 do jsou = 1, 3 4268 cofaf(isou) = cofaf(isou) & 4269 + normal(isou)*normal(jsou)*(qimpv(jsou)+hintpv(jsou)) 4270 if (jsou.eq.isou) then 4271 cofbf(isou,jsou) = hint(isou)-hintnm(isou)*normal(jsou) 4272 else 4273 cofbf(isou,jsou) = -hintnm(isou)*normal(jsou) 4274 endif 4275 enddo 4276 4277enddo 4278 4279return 4280end subroutine set_generalized_dirichlet_vector_aniso 4281 4282!=============================================================================== 4283 4284!------------------------------------------------------------------------------- 4285! Arguments 4286!______________________________________________________________________________. 4287! mode name role ! 4288!______________________________________________________________________________! 4289!> \param[out] coefa explicit BC coefficient for gradients 4290!> \param[out] cofaf explicit BC coefficient for diffusive flux 4291!> \param[out] coefb implicit BC coefficient for gradients 4292!> \param[out] cofbf implicit BC coefficient for diffusive flux 4293!> \param[in] pimp Flux value to impose 4294!> \param[in] cfl Local Courant number used to convect 4295!> \param[in] hint Internal exchange coefficient 4296!_______________________________________________________________________________ 4297 4298subroutine set_convective_outlet_scalar & 4299 ( coefa , cofaf, coefb , cofbf, pimp , cfl , hint) 4300 4301!=============================================================================== 4302! Module files 4303!=============================================================================== 4304 4305!=============================================================================== 4306 4307implicit none 4308 4309! Arguments 4310 4311double precision coefa, cofaf, coefb, cofbf, pimp, cfl, hint 4312 4313! Local variables 4314 4315!=============================================================================== 4316 4317! Gradient BCs 4318coefb = cfl/(1.d0+cfl) 4319coefa = (1.d0-coefb)*pimp 4320 4321! Flux BCs 4322cofaf = -hint*coefa 4323cofbf = hint*(1.d0 - coefb) 4324 4325return 4326end subroutine 4327 4328!=============================================================================== 4329 4330!------------------------------------------------------------------------------- 4331! Arguments 4332!______________________________________________________________________________. 4333! mode name role ! 4334!______________________________________________________________________________! 4335!> \param[out] coefa explicit BC coefficient for gradients 4336!> \param[out] cofaf explicit BC coefficient for diffusive flux 4337!> \param[out] coefb implicit BC coefficient for gradients 4338!> \param[out] cofbf implicit BC coefficient for diffusive flux 4339!> \param[in] pimpv Dirichlet value to impose 4340!> \param[in] cflv Local Courant number used to convect 4341!> \param[in] hint Internal exchange coefficient 4342!_______________________________________________________________________________ 4343 4344subroutine set_convective_outlet_vector & 4345 ( coefa , cofaf, coefb , cofbf, pimpv , cflv , hint) 4346 4347!=============================================================================== 4348! Module files 4349!=============================================================================== 4350 4351!=============================================================================== 4352 4353implicit none 4354 4355! Arguments 4356 4357double precision coefa(3), cofaf(3) 4358double precision coefb(3,3), cofbf(3,3) 4359double precision pimpv(3), cflv(3) 4360double precision hint 4361 4362! Local variables 4363 4364integer isou , jsou 4365 4366!=============================================================================== 4367 4368do isou = 1, 3 4369 4370 ! Gradient BCs 4371 do jsou = 1, 3 4372 if (jsou.eq.isou) then 4373 coefb(isou,jsou) = cflv(isou)/(1.d0+cflv(isou)) 4374 else 4375 coefb(isou,jsou) = 0.d0 4376 endif 4377 enddo 4378 coefa(isou) = (1.d0-coefb(isou,isou))*pimpv(isou) 4379 4380 ! Flux BCs 4381 cofaf(isou) = -hint*coefa(isou) 4382 do jsou = 1, 3 4383 if (jsou.eq.isou) then 4384 cofbf(isou,jsou) = hint*(1.d0 - coefb(isou,jsou)) 4385 else 4386 cofbf(isou,jsou) = 0.d0 4387 endif 4388 enddo 4389 4390enddo 4391 4392return 4393end subroutine 4394 4395!=============================================================================== 4396 4397!------------------------------------------------------------------------------- 4398! Arguments 4399!______________________________________________________________________________. 4400! mode name role ! 4401!______________________________________________________________________________! 4402!> \param[out] coefa explicit BC coefficient for gradients 4403!> \param[out] cofaf explicit BC coefficient for diffusive flux 4404!> \param[out] coefb implicit BC coefficient for gradients 4405!> \param[out] cofbf implicit BC coefficient for diffusive flux 4406!> \param[in] pimpts Dirichlet value to impose 4407!> \param[in] cflts Local Courant number used to convect 4408!> \param[in] hint Internal exchange coefficient 4409!_______________________________________________________________________________ 4410 4411subroutine set_convective_outlet_tensor & 4412 ( coefa , cofaf, coefb , cofbf, pimpts , cflts , hint) 4413 4414!=============================================================================== 4415! Module files 4416!=============================================================================== 4417 4418!=============================================================================== 4419 4420implicit none 4421 4422! Arguments 4423 4424double precision coefa(6), cofaf(6) 4425double precision coefb(6,6), cofbf(6,6) 4426double precision pimpts(6), cflts(6) 4427double precision hint 4428 4429! Local variables 4430 4431integer isou , jsou 4432 4433!=============================================================================== 4434 4435do isou = 1, 6 4436 4437 ! Gradient BCs 4438 do jsou = 1, 6 4439 if (jsou.eq.isou) then 4440 coefb(isou,jsou) = cflts(isou)/(1.d0+cflts(isou)) 4441 else 4442 coefb(isou,jsou) = 0.d0 4443 endif 4444 enddo 4445 coefa(isou) = (1.d0-coefb(isou,isou))*pimpts(isou) 4446 4447 ! Flux BCs 4448 cofaf(isou) = -hint*coefa(isou) 4449 do jsou = 1, 6 4450 if (jsou.eq.isou) then 4451 cofbf(isou,jsou) = hint*(1.d0 - coefb(isou,jsou)) 4452 else 4453 cofbf(isou,jsou) = 0.d0 4454 endif 4455 enddo 4456 4457enddo 4458 4459return 4460end subroutine 4461 4462 4463!=============================================================================== 4464 4465!------------------------------------------------------------------------------- 4466! Arguments 4467!______________________________________________________________________________. 4468! mode name role ! 4469!______________________________________________________________________________! 4470!> \param[out] coefa explicit BC coefficient for gradients 4471!> \param[out] cofaf explicit BC coefficient for diffusive flux 4472!> \param[out] coefb implicit BC coefficient for gradients 4473!> \param[out] cofbf implicit BC coefficient for diffusive flux 4474!> \param[in] pimpv Dirichlet value to impose 4475!> \param[in] cflv Local Courant number used to convect 4476!> \param[in] hint Internal exchange coefficient 4477!_______________________________________________________________________________ 4478 4479subroutine set_convective_outlet_vector_aniso & 4480 ( coefa , cofaf, coefb , cofbf, pimpv , cflv , hint ) 4481 4482!=============================================================================== 4483! Module files 4484!=============================================================================== 4485 4486!=============================================================================== 4487 4488implicit none 4489 4490! Arguments 4491 4492double precision coefa(3), cofaf(3) 4493double precision coefb(3,3), cofbf(3,3) 4494double precision pimpv(3), cflv(3) 4495double precision hint(6) 4496 4497! Local variables 4498 4499integer isou , jsou 4500 4501!=============================================================================== 4502 4503do isou = 1, 3 4504 4505 ! Gradient BCs 4506 do jsou = 1, 3 4507 if (jsou.eq.isou) then 4508 coefb(isou,jsou) = cflv(isou)/(1.d0+cflv(isou)) 4509 else 4510 coefb(isou,jsou) = 0.d0 4511 endif 4512 enddo 4513 coefa(isou) = (1.d0-coefb(isou,isou))*pimpv(isou) 4514 4515enddo 4516 4517! Flux BCs 4518cofaf(1) = -(hint(1)*coefa(1) + hint(4)*coefa(2) + hint(6)*coefa(3)) 4519cofaf(2) = -(hint(4)*coefa(1) + hint(2)*coefa(2) + hint(5)*coefa(3)) 4520cofaf(3) = -(hint(6)*coefa(1) + hint(5)*coefa(2) + hint(3)*coefa(3)) 4521cofbf(1,1) = hint(1)*(1.d0 - coefb(1,1)) 4522cofbf(2,2) = hint(2)*(1.d0 - coefb(2,2)) 4523cofbf(3,3) = hint(3)*(1.d0 - coefb(3,3)) 4524cofbf(1,2) = hint(4)*(1.d0 - coefb(1,1)) 4525cofbf(2,1) = hint(4)*(1.d0 - coefb(1,1)) 4526cofbf(2,3) = hint(5)*(1.d0 - coefb(2,2)) 4527cofbf(3,2) = hint(5)*(1.d0 - coefb(2,2)) 4528cofbf(1,3) = hint(6)*(1.d0 - coefb(3,3)) 4529cofbf(3,1) = hint(6)*(1.d0 - coefb(3,3)) 4530 4531return 4532end subroutine 4533 4534 4535!=============================================================================== 4536 4537!------------------------------------------------------------------------------- 4538! Arguments 4539!______________________________________________________________________________. 4540! mode name role ! 4541!______________________________________________________________________________! 4542!> \param[out] coefa explicit BC coefficient for gradients 4543!> \param[out] cofaf explicit BC coefficient for diffusive flux 4544!> \param[out] coefb implicit BC coefficient for gradients 4545!> \param[out] cofbf implicit BC coefficient for diffusive flux 4546!> \param[in] pinf affine part 4547!> \param[in] ratio linear part 4548!> \param[in] hint internal exchange coefficient 4549!_______________________________________________________________________________ 4550 4551subroutine set_affine_function_scalar & 4552 ( coefa , cofaf, coefb, cofbf, pinf , ratio, hint ) 4553 4554!=============================================================================== 4555! Module files 4556!=============================================================================== 4557 4558!=============================================================================== 4559 4560implicit none 4561 4562! Arguments 4563 4564double precision coefa, cofaf, coefb, cofbf, pinf, ratio, hint 4565 4566!=============================================================================== 4567 4568! Gradient BCs 4569coefb = ratio 4570coefa = pinf 4571 4572! Flux BCs 4573cofaf = -hint*coefa 4574cofbf = hint*(1.d0 - coefb) 4575 4576return 4577end subroutine 4578 4579 4580!=============================================================================== 4581 4582!------------------------------------------------------------------------------- 4583! Arguments 4584!______________________________________________________________________________. 4585! mode name role ! 4586!______________________________________________________________________________! 4587!> \param[out] coefa explicit BC coefficient for gradients 4588!> \param[out] cofaf explicit BC coefficient for diffusive flux 4589!> \param[out] coefb implicit BC coefficient for gradients 4590!> \param[out] cofbf implicit BC coefficient for diffusive flux 4591!> \param[in] dimp Flux value to impose 4592!> \param[in] hint Internal exchange coefficient 4593!_______________________________________________________________________________ 4594 4595subroutine set_neumann_conv_h_neumann_diff_scalar & 4596 (coefa , cofaf, coefb, cofbf, dimp, hint) 4597 4598!=============================================================================== 4599! Module files 4600!=============================================================================== 4601 4602!=============================================================================== 4603 4604implicit none 4605 4606! Arguments 4607 4608double precision coefa, cofaf, coefb, cofbf, dimp, hint 4609 4610!=============================================================================== 4611 4612! Gradient BCs 4613call set_neumann_scalar(coefa, cofaf, coefb, cofbf, dimp, hint) 4614 4615! Flux BCs 4616cofaf = 0.d0 4617cofbf = 0.d0 4618 4619return 4620end subroutine 4621 4622 4623!=============================================================================== 4624 4625!------------------------------------------------------------------------------- 4626! Arguments 4627!______________________________________________________________________________. 4628! mode name role ! 4629!______________________________________________________________________________! 4630!> \param[out] coefa explicit BC coefficient for gradients 4631!> \param[out] cofaf explicit BC coefficient for diffusive flux 4632!> \param[out] coefb implicit BC coefficient for gradients 4633!> \param[out] cofbf implicit BC coefficient for diffusive flux 4634!> \param[in] pinf affine part 4635!> \param[in] ratio linear part 4636!> \param[in] dimp Flux value to impose 4637!_______________________________________________________________________________ 4638 4639subroutine set_affine_function_conv_neumann_diff_scalar & 4640 (coefa, cofaf, coefb, cofbf, pinf, ratio, dimp) 4641 4642!=============================================================================== 4643! Module files 4644!=============================================================================== 4645 4646!=============================================================================== 4647 4648implicit none 4649 4650! Arguments 4651 4652double precision coefa, cofaf, coefb, cofbf, pinf, ratio, dimp 4653 4654!=============================================================================== 4655 4656! Gradient BCs 4657coefb = ratio 4658coefa = pinf 4659 4660! Flux BCs 4661cofaf = dimp 4662cofbf = 0.d0 4663 4664return 4665end subroutine 4666 4667 4668!=============================================================================== 4669 4670!------------------------------------------------------------------------------- 4671! Arguments 4672!______________________________________________________________________________. 4673! mode name role ! 4674!______________________________________________________________________________! 4675!> \param[out] coefa explicit BC coefficient for gradients 4676!> \param[out] cofaf explicit BC coefficient for diffusive flux 4677!> \param[out] coefb implicit BC coefficient for gradients 4678!> \param[out] cofbf implicit BC coefficient for diffusive flux 4679!> \param[in] hext convective flux to be imposed 4680!> \param[in] dimp Flux value to impose 4681!_______________________________________________________________________________ 4682 4683subroutine set_total_flux & 4684 ( coefa, cofaf, coefb, cofbf, hext, dimp ) 4685 4686!=============================================================================== 4687! Module files 4688!=============================================================================== 4689 4690!=============================================================================== 4691 4692implicit none 4693 4694! Arguments 4695 4696double precision coefa, cofaf, coefb, cofbf, dimp, hext 4697 4698! Gradients BCs 4699coefa = 0.d0 4700coefb = 1.d0 4701 4702! Flux BCs 4703cofaf = dimp 4704cofbf = hext 4705 4706return 4707end subroutine 4708 4709 4710!=============================================================================== 4711 4712!------------------------------------------------------------------------------- 4713! Arguments 4714!______________________________________________________________________________. 4715! mode name role ! 4716!______________________________________________________________________________! 4717!> \param[out] coefa explicit BC coefficient for gradients 4718!> \param[out] cofaf explicit BC coefficient for diffusive flux 4719!> \param[out] coefb implicit BC coefficient for gradients 4720!> \param[out] cofbf implicit BC coefficient for diffusive flux 4721!> \param[in] pimp Dirichlet value to impose 4722!> \param[in] dimp Flux value to impose 4723!_______________________________________________________________________________ 4724 4725subroutine set_dirichlet_conv_neumann_diff_scalar & 4726 ( coefa, cofaf, coefb, cofbf, pimp, dimp ) 4727 4728!=============================================================================== 4729! Module files 4730!=============================================================================== 4731 4732!=============================================================================== 4733 4734implicit none 4735 4736! Arguments 4737 4738double precision coefa, cofaf, coefb, cofbf, pimp, dimp 4739 4740! Gradients BCs 4741coefa = pimp 4742coefb = 0.d0 4743 4744! Flux BCs 4745cofaf = dimp 4746cofbf = 0.d0 4747 4748return 4749end subroutine 4750 4751!=============================================================================== 4752 4753!------------------------------------------------------------------------------- 4754! Arguments 4755!______________________________________________________________________________. 4756! mode name role ! 4757!______________________________________________________________________________! 4758!> \param[out] coefa explicit BC coefficient for gradients 4759!> \param[out] cofaf explicit BC coefficient for diffusive flux 4760!> \param[out] coefb implicit BC coefficient for gradients 4761!> \param[out] cofbf implicit BC coefficient for diffusive flux 4762!> \param[in] pimpv Dirichlet value to impose 4763!> \param[in] qimpv Flux value to impose 4764!_______________________________________________________________________________ 4765 4766subroutine set_dirichlet_conv_neumann_diff_vector & 4767 ( coefa, cofaf, coefb, cofbf, pimpv, qimpv ) 4768 4769!=============================================================================== 4770! Module files 4771!=============================================================================== 4772 4773!=============================================================================== 4774 4775implicit none 4776 4777! Arguments 4778 4779double precision coefa(3), cofaf(3) 4780double precision coefb(3,3), cofbf(3,3) 4781double precision pimpv(3), qimpv(3) 4782 4783! Local variables 4784 4785integer isou , jsou 4786 4787!=============================================================================== 4788 4789do isou = 1, 3 4790 4791 ! Gradient BCs 4792 coefa(isou) = pimpv(isou) 4793 do jsou = 1, 3 4794 coefb(isou,jsou) = 0.d0 4795 enddo 4796 4797 ! Flux BCs 4798 cofaf(isou) = qimpv(isou) 4799 do jsou = 1, 3 4800 cofbf(isou,jsou) = 0.d0 4801 enddo 4802 4803enddo 4804 4805return 4806end subroutine 4807 4808!------------------------------------------------------------------------------- 4809! Arguments 4810!______________________________________________________________________________. 4811! mode name role ! 4812!______________________________________________________________________________! 4813!> \param[out] coefa explicit BC coefficient for gradients 4814!> \param[out] cofaf explicit BC coefficient for diffusive flux 4815!> \param[out] coefb implicit BC coefficient for gradients 4816!> \param[out] cofbf implicit BC coefficient for diffusive flux 4817!> \param[in] pimpts Dirichlet value to impose 4818!> \param[in] qimpts Flux value to impose 4819!_______________________________________________________________________________ 4820 4821subroutine set_dirichlet_conv_neumann_diff_tensor & 4822 ( coefa, cofaf, coefb, cofbf, pimpts, qimpts ) 4823 4824!=============================================================================== 4825! Module files 4826!=============================================================================== 4827 4828!=============================================================================== 4829 4830implicit none 4831 4832! Arguments 4833 4834double precision coefa(6), cofaf(6) 4835double precision coefb(6,6), cofbf(6,6) 4836double precision pimpts(6), qimpts(6) 4837 4838! Local variables 4839 4840integer isou , jsou 4841 4842!=============================================================================== 4843 4844do isou = 1, 6 4845 4846 ! BS test sur hextv ? if (abs(hextv(isou)).gt.rinfin*0.5d0) then 4847 4848 ! Gradient BCs 4849 coefa(isou) = pimpts(isou) 4850 do jsou = 1, 6 4851 coefb(isou,jsou) = 0.d0 4852 enddo 4853 4854 ! Flux BCs 4855 cofaf(isou) = qimpts(isou) 4856 do jsou = 1, 6 4857 cofbf(isou,jsou) = 0.d0 4858 enddo 4859 4860enddo 4861 4862return 4863end subroutine 4864 4865!=============================================================================== 4866! Function : 4867! -------- 4868 4869!> \brief Initialisation of boundary condition arrays. 4870!> 4871!------------------------------------------------------------------------------- 4872 4873!------------------------------------------------------------------------------- 4874! Arguments 4875!______________________________________________________________________________. 4876! mode name role ! 4877!______________________________________________________________________________! 4878!> \param[in] nvar total number of variables 4879!> \param[in] nscal total number of scalars 4880!> \param[in] itrale ALE iteration number 4881!> \param[in,out] icodcl face boundary condition code (see \ref condli) 4882!> \param[in,out] isostd indicator for standard outlet 4883!> and reference face index 4884!> \param[in] dt time step (per cell) 4885!> \param[in,out] rcodcl boundary condition values: 4886!> - rcodcl(1) value of the dirichlet 4887!> - rcodcl(2) value of the exterior exchange 4888!> coefficient (infinite if no exchange) 4889!> - rcodcl(3) value flux density 4890!> (negative if gain) in w/m2 or roughness 4891!> in m if icodcl=6 4892!> -# for the velocity \f$ (\mu+\mu_T) 4893!> \gradv \vect{u} \cdot \vect{n} \f$ 4894!> -# for the pressure \f$ \Delta t 4895!> \grad P \cdot \vect{n} \f$ 4896!> -# for a scalar \f$ cp \left( K + 4897!> \dfrac{K_T}{\sigma_T} \right) 4898!> \grad T \cdot \vect{n} \f$ 4899!_______________________________________________________________________________ 4900 4901subroutine condli_ini & 4902 ( nvar , nscal , & 4903 itrale , & 4904 icodcl , isostd , & 4905 dt , rcodcl ) 4906 4907!=============================================================================== 4908! Module files 4909!=============================================================================== 4910 4911use paramx 4912use numvar 4913use optcal 4914use alaste 4915use alstru 4916use atincl, only: iautom, iprofm 4917use coincl, only: fment, ientfu, ientgb, ientgf, ientox, tkent, qimp 4918use ppcpfu, only: inmoxy 4919use cstphy 4920use cstnum 4921use pointe 4922use entsor 4923use albase 4924use parall 4925use ppppar 4926use ppthch 4927use ppincl 4928use cpincl 4929use radiat 4930use cplsat 4931use mesh 4932use field 4933use field_operator 4934use radiat 4935use turbomachinery 4936use darcy_module 4937use cs_c_bindings 4938 4939!=============================================================================== 4940 4941implicit none 4942 4943! Arguments 4944 4945integer nvar , nscal 4946integer itrale 4947 4948integer icodcl(nfabor,nvar) 4949integer isostd(nfabor+1) 4950 4951double precision, pointer, dimension(:) :: dt 4952double precision rcodcl(nfabor,nvar,3) 4953 4954! Local variables 4955 4956integer iterns, ii 4957 4958double precision, dimension(:,:), pointer :: disale 4959double precision, dimension(:,:), pointer :: xyzno0 4960 4961!=============================================================================== 4962! 0. User calls 4963!=============================================================================== 4964 4965call precli(nvar, icodcl, rcodcl) 4966 4967! - Interface Code_Saturne 4968! ====================== 4969 4970! N.B. Zones de face de bord : on utilise provisoirement les zones des 4971! physiques particulieres, meme sans physique particuliere 4972! -> sera modifie lors de la restructuration des zones de bord 4973 4974call uiclim & 4975 ( nozppm, & 4976 iqimp, icalke, ientat, ientcp, inmoxy, ientox, & 4977 ientfu, ientgb, ientgf, iprofm, iautom, & 4978 itypfb, izfppp, icodcl, & 4979 qimp, qimpat, qimpcp, dh, xintur, & 4980 timpat, timpcp, tkent , fment, distch, nvar, rcodcl) 4981 4982! - Sous-programme utilisateur 4983! ========================== 4984 4985call cs_f_user_boundary_conditions & 4986 ( nvar , nscal , & 4987 icodcl , itrifb , itypfb , izfppp , & 4988 dt , & 4989 rcodcl ) 4990 4991call user_boundary_conditions(nvar, itypfb, icodcl, rcodcl) 4992 4993! -- Methode ALE (CL de vitesse de maillage et deplacement aux noeuds) 4994 4995if (iale.ge.1) then 4996 4997 call field_get_val_v(fdiale, disale) 4998 call field_get_val_v_by_name("vtx_coord0", xyzno0) 4999 5000 do ii = 1, nnod 5001 impale(ii) = 0 5002 enddo 5003 5004 ! - Interface Code_Saturne 5005 ! ====================== 5006 5007 call uialcl & 5008 ( ibfixe, igliss, ivimpo, ifresf, & 5009 ialtyb, & 5010 impale, & 5011 disale, & 5012 iuma, ivma, iwma, & 5013 rcodcl) 5014 5015 call usalcl & 5016 ( itrale , & 5017 nvar , nscal , & 5018 icodcl , itypfb , ialtyb , & 5019 impale , & 5020 dt , & 5021 rcodcl , xyzno0 , disale ) 5022 5023 ! Au cas ou l'utilisateur aurait touche disale sans mettre impale=1, on 5024 ! remet le deplacement initial 5025 do ii = 1, nnod 5026 if (impale(ii).eq.0) then 5027 disale(1,ii) = xyznod(1,ii)-xyzno0(1,ii) 5028 disale(2,ii) = xyznod(2,ii)-xyzno0(2,ii) 5029 disale(3,ii) = xyznod(3,ii)-xyzno0(3,ii) 5030 endif 5031 enddo 5032 5033endif 5034 5035! For internal coupling, set itypfb to wall function by default 5036! if not set by the user 5037call cs_internal_coupling_bcs(itypfb) 5038 5039!=============================================================================== 5040! 2. treatment of types of bcs given by itypfb 5041!=============================================================================== 5042 5043if (iale.ge.1) then 5044 call altycl & 5045 ( itypfb , ialtyb , icodcl , impale , .true. , & 5046 dt , & 5047 rcodcl ) 5048endif 5049 5050if (iturbo.ne.0) then 5051 call mmtycl(itypfb, rcodcl) 5052endif 5053 5054iterns = 1 5055 5056! Locate internal BC-based coupling 5057if (nbrcpl.gt.0) then 5058 call cscloc 5059 call cscfbr_init(icodcl, itypfb) 5060endif 5061 5062if ( ippmod(iphpar).ge.1.and.ippmod(igmix).eq.-1 & 5063 .and.ippmod(ieljou).eq.-1.and.ippmod(ielarc).eq.-1 & 5064 .or.ippmod(icompf).ge.0.and.ippmod(igmix).ge.0) then 5065 call pptycl & 5066 ( nvar , .true., & 5067 icodcl , itypfb , izfppp , & 5068 dt , & 5069 rcodcl ) 5070endif 5071 5072call typecl & 5073 ( nvar , nscal , .true. , & 5074 itypfb , itrifb , icodcl , isostd , & 5075 rcodcl ) 5076 5077!=============================================================================== 5078! 3. check the consistency of the bcs 5079!=============================================================================== 5080 5081! When called before time loop, some values are not yet available. 5082if (ntcabs .eq. ntpabs) return 5083 5084call vericl(nvar, nscal, itypfb, icodcl) 5085 5086!---- 5087! End 5088!---- 5089 5090return 5091end subroutine condli_ini 5092