1!------------------------------------------------------------------------------- 2 3! This file is part of Code_Saturne, a general-purpose CFD tool. 4! 5! Copyright (C) 1998-2021 EDF S.A. 6! 7! This program is free software; you can redistribute it and/or modify it under 8! the terms of the GNU General Public License as published by the Free Software 9! Foundation; either version 2 of the License, or (at your option) any later 10! version. 11! 12! This program is distributed in the hope that it will be useful, but WITHOUT 13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15! details. 16! 17! You should have received a copy of the GNU General Public License along with 18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 19! Street, Fifth Floor, Boston, MA 02110-1301, USA. 20 21!------------------------------------------------------------------------------- 22 23subroutine cfmsfp & 24!================ 25 26 ( nvar , nscal , iterns , ncepdp , ncesmp , & 27 icepdc , icetsm , itypsm , & 28 dt , vela , & 29 ckupdc , smacel , & 30 flumas , flumab ) 31 32!=============================================================================== 33! FONCTION : 34! ---------- 35 36! "MASS FLUX" AT THE FACES CALCULATION FOR THE CFL RESTRICTION CALCULATION 37! AND THE SOLVING OF THE PRESSURE 38 39!------------------------------------------------------------------------------- 40!ARGU ARGUMENTS 41!__________________.____._____.________________________________________________. 42! name !type!mode ! role ! 43!__________________!____!_____!________________________________________________! 44! nvar ! i ! <-- ! total number of variables ! 45! nscal ! i ! <-- ! total number of scalars ! 46! iterns ! i ! <-- ! Navier-Stokes iteration number ! 47! ncepdp ! i ! <-- ! number of cells with head loss ! 48! ncesmp ! i ! <-- ! number of cells with mass source term ! 49! icepdc(ncelet) ! te ! <-- ! numero des ncepdp cellules avec pdc ! 50! icetsm(ncesmp) ! te ! <-- ! numero des cellules a source de masse ! 51! itypsm ! te ! <-- ! type de source de masse pour les ! 52! dt(ncelet) ! ra ! <-- ! time step (per cell) ! 53! vela ! ra ! <-- ! variable value at time step beginning ! 54! ckupdc ! tr ! <-- ! work array for the head loss ! 55! (ncepdp,6) ! ! ! ! 56! smacel ! tr ! <-- ! variable value associated to the mass source ! 57! (ncesmp,*) ! ! ! term (for ivar=ipr, smacel is the mass flux ! 58! ! ! ! \f$ \Gamma^n \f$) ! 59! flumas(nfac) ! tr ! --> ! flux de masse aux faces internes ! 60! flumab(nfabor) ! tr ! --> ! flux de masse aux faces de bord ! 61!__________________!____!_____!________________________________________________! 62 63! Type: i (integer), r (real), s (string), a (array), l (logical), 64! and composite types (ex: ra real array) 65! mode: <-- input, --> output, <-> modifies data, --- work array 66!=============================================================================== 67 68!=============================================================================== 69! Module files 70!=============================================================================== 71 72use cfpoin 73use paramx 74use numvar 75use entsor 76use optcal 77use cstphy 78use cstnum 79use parall 80use period 81use ppppar 82use ppthch 83use ppincl 84use mesh 85use field 86use cs_f_interfaces 87use cs_c_bindings 88 89!=============================================================================== 90 91implicit none 92 93! Arguments 94 95integer nvar , nscal, iterns 96integer ncepdp , ncesmp 97 98integer icepdc(ncepdp) 99integer icetsm(ncesmp), itypsm(ncesmp,nvar) 100 101double precision dt(ncelet) 102double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar) 103double precision flumas(nfac), flumab(nfabor) 104double precision vela (3 ,ncelet) 105 106! Local variables 107 108integer ifac , iel, ischcp, idftnp, ircflp 109integer init , inc , iccocg, isstpp 110integer imrgrp, nswrgp, imligp, iwarnp, iconvp, idiffp, imvisp 111integer icvflb, f_id0 112integer isou , jsou 113integer iflmb0, itypfl 114integer itsqdm, imasac 115 116double precision epsrgp, climgp, thetap, blencp, relaxp 117double precision rom 118 119double precision, allocatable, dimension(:) :: w1 120double precision, dimension(:), pointer :: crom, brom 121double precision, dimension(:), pointer :: viscl, visct 122double precision, dimension(:,:), allocatable :: tsexp, gavinj 123double precision, allocatable, dimension(:,:) :: vel0 124double precision, dimension(:,:,:), allocatable :: tsimp 125double precision, allocatable, dimension(:,:,:), target :: viscf 126double precision, allocatable, dimension(:,:) :: viscce 127double precision, allocatable, dimension(:) :: viscb 128double precision, allocatable, dimension(:) :: secvif, secvib 129double precision, allocatable, dimension(:,:,:) :: coefbv 130 131double precision, dimension(:,:), pointer :: coefau, cofafu 132double precision, dimension(:,:,:), pointer :: coefbu, cofbfu 133 134type(var_cal_opt) :: vcopt_u, vcopt_p 135type(var_cal_opt), target :: vcopt_loc 136type(var_cal_opt), pointer :: p_k_value 137type(c_ptr) :: c_k_value 138 139double precision rvoid(1) 140 141!=============================================================================== 142 143!=============================================================================== 144! 1. INITIALISATION 145!=============================================================================== 146 147call field_get_coefa_v(ivarfl(iu), coefau) 148call field_get_coefb_v(ivarfl(iu), coefbu) 149call field_get_coefaf_v(ivarfl(iu), cofafu) 150call field_get_coefbf_v(ivarfl(iu), cofbfu) 151 152! Allocate work arrays 153allocate(w1(ncelet)) 154allocate(tsexp(3,ncelet)) 155allocate(gavinj(3,ncelet)) 156allocate(tsimp(3,3,ncelet)) 157allocate(coefbv(3,3,nfabor)) 158allocate(vel0(3,ncelet)) 159 160call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt_u) 161call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt_p) 162 163if (iand(vcopt_u%idften, ISOTROPIC_DIFFUSION).ne.0) then 164 allocate(viscf(1, 1, nfac), viscb(nfabor)) 165else if (iand(vcopt_u%idften, ANISOTROPIC_LEFT_DIFFUSION).ne.0) then 166 allocate(viscf(3, 3, nfac), viscb(nfabor)) 167 allocate(viscce(6,ncelet)) 168endif 169 170if (ivisse.eq.1) then 171 allocate(secvif(nfac),secvib(nfabor)) 172endif 173 174f_id0 = -1 175 176! Density 177 178call field_get_val_s(icrom, crom) 179call field_get_val_s(ibrom, brom) 180 181!=============================================================================== 182! 2. MASS FLUX AT THE FACES 183!=============================================================================== 184 185! 2.1 SOURCE TERMS OF THE MOMENTUM EQUATIONS 186! ========================================== 187 188! Some first tests (double expansion waves in a shock tube) 189! has shown that taking into account all the 190! momentum equation terms in the mass equation seems to be a 191! bad idea (in particular the convective term, but the diffusive 192! term, the transposed gradient, the mass and user source terms 193! were all null in the considered tests). 194! However, it may be due to a bug at that early stage of implementation 195! of the algorithm (but we didn't find it). 196! We thus recommand not to take into account the momentum source terms, 197! except the gravity term (because it is in balance with the pressure 198! gradient and because its effect is visible at equilibrium). 199! However, we keep here the implementation of the preliminary tests 200! (1.1.0.h version) with an overall test so that the correction is not 201! active (thus, there is no user question and there is always the 202! possibility to perform other tests in the future). 203! Note that, with these terms, the thoeretical analysis is harder 204! (Without these terms we are in the configuration Euler + gravity) 205 206! --- Initialization 207do iel = 1, ncel 208 do isou = 1, 3 209 tsexp(isou,iel) = 0.d0 210 do jsou = 1, 3 211 tsimp(isou,jsou,iel) = 0.d0 212 enddo 213 enddo 214enddo 215 216! Test on momentum source terms 217itsqdm = 0 218if (itsqdm.ne.0) then 219 220 if (ivisse.eq.1) then 221 222 call visecv & 223 ( secvif , secvib ) 224 225 endif 226 227 228 ! --- User source term 229 call ustsnv & 230 ( nvar , nscal , ncepdp , ncesmp , & 231 iu , & 232 icepdc , icetsm , itypsm , & 233 dt , & 234 ckupdc , smacel , tsexp , tsimp ) 235 236 ! C version 237 call user_source_terms(ivarfl(iu), tsexp, tsimp) 238 239 ! Convective of the momentum equation 240 ! in upwind and without reconstruction 241 iconvp = vcopt_u%iconv 242 243 init = 1 244 inc = 1 245 iccocg = 1 246 iflmb0 = 1 247 imrgrp = vcopt_u%imrgra 248 nswrgp = vcopt_u%nswrgr 249 imligp = vcopt_u%imligr 250 iwarnp = vcopt_u%iwarni 251 epsrgp = vcopt_u%epsrgr 252 climgp = vcopt_u%climgr 253 relaxp = vcopt_u%relaxv 254 thetap = vcopt_u%thetav 255 256 itypfl = 1 257 258 ! Mass flux calculation 259 call inimav & 260 !========== 261( ivarfl(iu) , itypfl , & 262 iflmb0 , init , inc , imrgrp , nswrgp , imligp , & 263 iwarnp , & 264 epsrgp , climgp , & 265 crom, brom, & 266 vela, & 267 coefau , coefbu , & 268 flumas , flumab ) 269 270 ! ---> Face diffusivity for the velocity 271 idiffp = vcopt_u%idiff 272 imvisp = vcopt_u%imvisf 273 if (idiffp.ge. 1) then 274 275 call field_get_val_s(iviscl, viscl) 276 call field_get_val_s(ivisct, visct) 277 278 if (itytur.eq.3) then 279 do iel = 1, ncel 280 w1(iel) = viscl(iel) 281 enddo 282 else 283 do iel = 1, ncel 284 w1(iel) = viscl(iel) + vcopt_u%idifft*visct(iel) 285 enddo 286 endif 287 288 ! Scalar diffusivity (Default) 289 if (iand(vcopt_u%idften, ISOTROPIC_DIFFUSION).ne.0) then 290 291 call viscfa & 292 !========== 293 ( imvisf , & 294 w1 , & 295 viscf , viscb ) 296 297 ! Tensorial diffusion of the velocity (in case of tensorial porosity) 298 else if (iand(vcopt_u%idften, ANISOTROPIC_LEFT_DIFFUSION).ne.0) then 299 300 do iel = 1, ncel 301 do isou = 1, 3 302 viscce(isou, iel) = w1(iel) 303 enddo 304 do isou = 4, 6 305 viscce(isou, iel) = 0.d0 306 enddo 307 enddo 308 309 call vistnv (imvisf, viscce, viscf, viscb) 310 311 endif 312 313 ! --- If no diffusion, viscosity is set to 0. 314 else 315 316 do ifac = 1, nfac 317 viscf(1,1,ifac) = 0.d0 318 enddo 319 do ifac = 1, nfabor 320 viscb(ifac) = 0.d0 321 enddo 322 323 endif 324 325 idftnp = vcopt_u%idften 326 327 ! no recontruction 328 ircflp = 0; 329 330 ! upwind 331 ischcp = 0; 332 blencp = 0; 333 isstpp = 0; 334 335 inc = 1; 336 337 icvflb = 0; 338 339 ! The added convective scalar mass flux is: 340 ! (thetap*Y_\face-imasac*Y_\celli)*mf. 341 ! When building the implicit part of the rhs, one 342 ! has to impose 1 on mass accumulation. 343 imasac = 1 344 345 vcopt_loc = vcopt_u 346 347 vcopt_loc%istat = -1 348 vcopt_loc%idifft = -1 349 vcopt_loc%iswdyn = -1 350 vcopt_loc%nswrsm = -1 351 vcopt_loc%iwgrec = 0 352 vcopt_loc%blend_st = 0 ! Warning, may be overwritten if a field 353 vcopt_loc%epsilo = -1 354 vcopt_loc%epsrsm = -1 355 356 p_k_value => vcopt_loc 357 c_k_value = equation_param_from_vcopt(c_loc(p_k_value)) 358 359 call cs_balance_vector & 360 (idtvar, ivarfl(iu), imasac, inc, ivisse, & 361 c_k_value, vela, vela, coefau, coefbu, cofafu, cofbfu, & 362 flumas, flumab, viscf, viscb, secvif, secvib, & 363 rvoid, rvoid, rvoid, & 364 icvflb, icvfli, tsexp) 365 366endif 367 368! End of the test on momentum source terms 369 370! Mass source term 371if (ncesmp.gt.0) then 372 373 ! The momentum balance is used in its conservative form here 374 ! so the mass source term is only composed of gamma*uinj 375 ! => array of previous velocity has to be set to zero 376 377 do iel = 1, ncel 378 do isou = 1,3 379 vel0(isou,iel) = 0.d0 380 enddo 381 enddo 382 383 call catsmv(ncesmp, iterns, icetsm, itypsm(:,iu), & 384 cell_f_vol, vel0, smacel(:,iu), smacel(:,ipr), & 385 tsexp, tsimp, gavinj) 386 387 do iel = 1, ncel 388 do isou = 1, 3 389 tsexp(isou,iel) = tsexp(isou,iel) + gavinj(isou,iel) 390 enddo 391 enddo 392 393endif 394 395! --- Volumic forces term (gravity) 396do iel = 1, ncel 397 rom = crom(iel) 398 tsexp(1,iel) = gx + tsexp(1,iel)/rom 399 tsexp(2,iel) = gy + tsexp(2,iel)/rom 400 tsexp(3,iel) = gz + tsexp(3,iel)/rom 401enddo 402 403! --- Calculation of the convective "velocities at the cell centers 404! (Calculation of u^n+dt*f^n) 405 406do iel = 1, ncel 407 do isou = 1, 3 408 tsexp(isou,iel) = dt(iel)*tsexp(isou,iel) 409 enddo 410enddo 411 412! Computation of the flux 413 414! volumic flux part based on dt*f^n 415 416! Initialization of the mass flux 417init = 1 418! As mentioned above, for homogeneous Neumann 419inc = 0 420iccocg = 1 421iflmb0 = 1 422! Reconstruction is useless here 423nswrgp = 0 ! FIXME 424imrgrp = vcopt_p%imrgra 425imligp = vcopt_p%imligr 426iwarnp = vcopt_p%iwarni 427epsrgp = vcopt_p%epsrgr 428climgp = vcopt_p%climgr 429 430! Velocity flux (crom, brom not used) 431itypfl = 0 432 433! No contribution of f to the boundary mass flux 434do ifac = 1, nfabor 435 do isou = 1, 3 436 do jsou = 1, 3 437 coefbv(isou,jsou,ifac) = 0.d0 438 enddo 439 enddo 440enddo 441 442call inimav & 443!========== 444( f_id0 , itypfl , & 445 iflmb0 , init , inc , imrgrp , nswrgp , imligp , & 446 iwarnp , & 447 epsrgp , climgp , & 448 crom, brom, & 449 tsexp, & 450 coefau , coefbv , & 451 flumas , flumab ) 452 453! volumic flux part based on velocity u^n 454init = 0 455! take into account Dirichlet velocity boundary conditions 456inc = 1 457 458call inimav & 459!========== 460( ivarfl(iu) , itypfl , & 461 iflmb0 , init , inc , imrgrp , nswrgp , imligp , & 462 iwarnp , & 463 epsrgp , climgp , & 464 crom, brom, & 465 vela, & 466 coefau , coefbu , & 467 flumas , flumab ) 468 469! Free memory 470deallocate(w1) 471deallocate(tsexp) 472deallocate(gavinj) 473deallocate(tsimp) 474deallocate(viscf, viscb) 475if (allocated(secvif)) deallocate(secvif, secvib) 476if (allocated(viscce)) deallocate(viscce) 477deallocate(coefbv) 478deallocate(vel0) 479 480!-------- 481! FORMATS 482!-------- 483 484 485!---- 486! FIN 487!---- 488 489return 490 491end subroutine 492