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 predvv.f90 28!> 29!> \brief This subroutine performs the velocity prediction step of the Navier 30!> Stokes equations for incompressible or slightly compressible flows for 31!> the coupled velocity components solver. 32!> 33!> - at the first call, the predicted velocities are computed and also 34!> an estimator on the predicted velocity is computed. 35!> 36!> - at the second call, a global estimator on Navier Stokes is computed. 37!> This second call is done after the correction step 38!> (\ref cs_pressure_correction). 39!> 40!> Please refer to the 41!> <a href="../../theory.pdf#predvv"><b>predvv</b></b></a> section 42!> of the theory guide for more informations. 43!------------------------------------------------------------------------------- 44 45!------------------------------------------------------------------------------- 46! Arguments 47!______________________________________________________________________________. 48! mode name role ! 49!______________________________________________________________________________! 50!> \param[in] iappel call number (1 or 2) 51!> \param[in] nvar total number of variables 52!> \param[in] nscal total number of scalars 53!> \param[in] iterns index of the iteration on Navier-Stokes 54!> \param[in] ncepdp number of cells with head loss 55!> \param[in] ncesmp number of cells with mass source term 56!> \param[in] icepdc index of cells with head loss 57!> \param[in] icetsm index of cells with mass source term 58!> \param[in] itypsm type of mass source term for the variables 59!> \param[in] dt time step (per cell) 60!> \param[in] vel velocity 61!> \param[in] vela velocity at the previous time step 62!> \param[in] velk velocity at the previous sub iteration (or vela) 63!> \param[in,out] da_uu velocity matrix 64!> \param[in] tslagr coupling term for the Lagrangian module 65!> \param[in] coefav boundary condition array for the variable 66!> (explicit part) 67!> \param[in] coefbv boundary condition array for the variable 68!> (implicit part) 69!> \param[in] cofafv boundary condition array for the diffusion 70!> of the variable (explicit part) 71!> \param[in] cofbfv boundary condition array for the diffusion 72!> of the variable (implicit part) 73!> \param[in] ckupdc work array for the head loss 74!> \param[in] smacel variable value associated to the mass source 75!> term (for ivar=ipr, smacel is the mass flux 76!> \f$ \Gamma^n \f$) 77!> \param[in] frcxt external forces making hydrostatic pressure 78!> \param[in] trava working array for the velocity-pressure coupling 79!> \param[out] dfrcxt variation of the external forces 80! making the hydrostatic pressure 81!> \param[in] tpucou non scalar time step in case of 82!> velocity pressure coupling 83!> \param[in] trav right hand side for the normalizing 84!> the residual 85!> \param[in] viscf visc*surface/dist aux faces internes 86!> \param[in] viscb visc*surface/dist aux faces de bord 87!> \param[in] viscfi same as viscf for increments 88!> \param[in] viscbi same as viscb for increments 89!> \param[in] secvif secondary viscosity at interior faces 90!> \param[in] secvib secondary viscosity at boundary faces 91!> \param[in] w1 working array 92!_______________________________________________________________________________ 93 94subroutine predvv & 95 ( iappel , & 96 nvar , nscal , iterns , & 97 ncepdp , ncesmp , & 98 icepdc , icetsm , itypsm , & 99 dt , vel , vela , velk , da_uu , & 100 tslagr , coefav , coefbv , cofafv , cofbfv , & 101 ckupdc , smacel , frcxt , & 102 trava , dfrcxt , tpucou , trav , & 103 viscf , viscb , viscfi , viscbi , secvif , secvib , & 104 w1 ) 105 106!=============================================================================== 107 108!=============================================================================== 109! Module files 110!=============================================================================== 111 112use, intrinsic :: iso_c_binding 113 114use paramx 115use numvar 116use entsor 117use cstphy 118use cstnum 119use optcal 120use parall 121use period 122use lagran 123use ppppar 124use ppthch 125use ppincl 126use cplsat 127use mesh 128use rotation 129use turbomachinery 130use cs_f_interfaces 131use cs_c_bindings 132use cfpoin 133use field 134use field_operator 135use cavitation 136use vof 137use atincl, only: kopint, iatmst, ps 138 139!=============================================================================== 140 141implicit none 142 143! Arguments 144 145integer iappel 146integer nvar , nscal , iterns 147integer ncepdp , ncesmp 148 149integer icepdc(ncepdp) 150integer icetsm(ncesmp), itypsm(ncesmp,nvar) 151 152double precision dt(ncelet) 153double precision tslagr(ncelet,*) 154double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar) 155double precision frcxt(3,ncelet), dfrcxt(3,ncelet) 156double precision trava(ndim,ncelet) 157double precision tpucou(6, ncelet) 158double precision trav(3,ncelet) 159double precision viscf(*), viscb(nfabor) 160double precision viscfi(*), viscbi(nfabor) 161double precision secvif(nfac), secvib(nfabor) 162double precision w1(ncelet) 163double precision coefav(3 ,nfabor) 164double precision cofafv(3 ,nfabor) 165double precision coefbv(3,3,nfabor) 166double precision cofbfv(3,3,nfabor) 167 168double precision vel (3, ncelet) 169double precision velk (3, ncelet) 170double precision vela (3, ncelet) 171double precision da_uu (6, ncelet) 172 173! Local variables 174 175integer f_id , iel , ielpdc, ifac , isou , itypfl, n_fans 176integer iccocg, inc , iprev , init , ii , jj 177integer imrgrp, nswrgp, imligp, iwarnp 178integer idftnp 179integer nswrsp, imvisp 180integer iescap 181integer iflmb0 182integer idtva0, icvflb, f_oi_id 183integer jsou , ivisep, imasac 184integer f_dim , iflwgr 185integer iflmas, iflmab 186integer ivoid(1) 187 188double precision rnorm , vitnor 189double precision romvom, drom , rom 190double precision epsrgp, climgp 191double precision vit1 , vit2 , vit3, xkb, pip, pfac 192double precision cpdc11, cpdc22, cpdc33, cpdc12, cpdc13, cpdc23 193double precision d2s3 , thetap, thetp1, thets 194double precision diipbx, diipby, diipbz 195double precision dvol 196double precision tensor(6) 197double precision rscp, tref 198 199double precision rvoid(1) 200 201! Working arrays 202double precision, allocatable, dimension(:,:) :: eswork 203double precision, allocatable, dimension(:,:), target :: grad 204double precision, allocatable, dimension(:,:), target :: hl_exp 205double precision, dimension(:,:), allocatable :: smbr 206double precision, dimension(:,:,:), allocatable :: fimp 207double precision, dimension(:,:,:), allocatable :: fimpcp 208double precision, dimension(:,:), allocatable :: gavinj 209double precision, dimension(:,:), allocatable :: tsexp 210double precision, dimension(:,:,:), allocatable :: tsimp 211double precision, allocatable, dimension(:,:) :: viscce 212double precision, dimension(:,:), allocatable :: vect 213double precision, dimension(:), pointer :: crom, croma, cromaa, pcrom 214double precision, dimension(:), pointer :: brom_eos, crom_eos, brom, broma 215double precision, dimension(:), allocatable, target :: cproa_rho_tc 216double precision, dimension(:), pointer :: coefa_k, coefb_k 217double precision, dimension(:), pointer :: coefa_p, coefb_p 218double precision, dimension(:,:), allocatable :: rij 219double precision, dimension(:,:), pointer :: coefap 220double precision, dimension(:,:,:), pointer :: coefbp 221double precision, dimension(:,:), allocatable :: coefat 222double precision, dimension(:,:,:), allocatable :: coefbt 223double precision, dimension(:,:), allocatable :: tflmas, tflmab 224double precision, allocatable, dimension(:,:), target :: divt 225double precision, dimension(:,:), pointer :: forbr, c_st_vel 226double precision, dimension(:), pointer :: cvar_pr, cvara_k 227double precision, dimension(:,:), pointer :: cvara_rij 228double precision, dimension(:), pointer :: viscl, visct, c_estim 229double precision, dimension(:,:), pointer :: lapla, lagr_st_vel 230double precision, dimension(:,:), pointer :: cpro_gradp 231double precision, dimension(:,:), pointer :: cpro_divr 232double precision, dimension(:,:), pointer :: cpro_pred_vel 233double precision, dimension(:), pointer :: cpro_wgrec_s, wgrec_crom 234double precision, dimension(:,:), pointer :: cpro_wgrec_v 235double precision, dimension(:), pointer :: imasfl, bmasfl 236double precision, dimension(:), pointer :: imasfl_prev, bmasfl_prev 237double precision, dimension(:), pointer :: cpro_beta, cvar_t 238double precision, dimension(:), allocatable, target :: cpro_rho_tc 239double precision, dimension(:), pointer :: cpro_rho_mass 240 241type(var_cal_opt) :: vcopt_p, vcopt_u, vcopt 242type(var_cal_opt), target :: vcopt_loc 243type(var_cal_opt), pointer :: p_k_value 244type(c_ptr) :: c_k_value 245 246!=============================================================================== 247 248!=============================================================================== 249! 1. Initialization 250!=============================================================================== 251 252! Id of the mass flux 253call field_get_key_int(ivarfl(iu), kimasf, iflmas) 254call field_get_key_int(ivarfl(iu), kbmasf, iflmab) 255 256! Pointers to the mass fluxes 257call field_get_val_s(iflmas, imasfl) 258call field_get_val_s(iflmab, bmasfl) 259 260! Pointers to properties 261! Density at time n+1,iteration iterns+1 262call field_get_val_s(icrom, crom_eos) 263call field_get_val_s(ibrom, brom_eos) 264 265! Density at time (n) 266if (irovar.eq.1) then 267 call field_get_val_prev_s(icrom, croma) 268 call field_get_val_prev_s(ibrom, broma) 269else 270 croma => crom_eos 271 broma => brom_eos 272endif 273 274! Density at time (n-1) if needed 275if ((idilat.gt.1.or.ivofmt.gt.0.or.ippmod(icompf).eq.3).or.irovar.eq.1) then 276 call field_get_val_prev2_s(icrom, cromaa) 277endif 278 279! Density for the unsteady term (at time n) 280! Compressible algorithm (mass equation is already solved) 281! or Low Mach compressible algos with mass flux prediction 282if (( (ippmod(icompf).ge.0.and.ippmod(icompf).ne.3) & 283 .or.(idilat.gt.1.and.ipredfl.eq.1)) & 284 .and.irovar.eq.1) then 285 pcrom => croma 286 287! VOF algorithm and Low Mach compressible algos: density at time n-1 288else if ((idilat.gt.1.or.ivofmt.gt.0.or.ippmod(icompf).eq.3) & 289 .and.irovar.eq.1) then 290 if (itpcol.eq.0.and.iterns.eq.1) then 291 pcrom => cromaa 292 else 293 pcrom => croma 294 endif 295 296! Weakly variable density algo. (idilat <=1) or constant density 297else 298 pcrom => crom_eos 299endif 300 301! Allocate temporary arrays 302allocate(smbr(3,ncelet)) 303allocate(fimp(3,3,ncelet)) 304allocate(fimpcp(3,3,ncelet)) 305allocate(tsexp(3,ncelet)) 306allocate(tsimp(3,3,ncelet)) 307call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt_u) 308call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt_p) 309 310! Density for other terms such as buoyancy term 311! 2nd order in time 312if (vcopt_u%thetav .lt. 1.d0) then 313 ! map the density pointer: 314 ! 1/4(n-1) + 1/2(n) + 1/4(n+1) 315 ! here replaced by (n) 316 crom => croma 317 brom => broma 318! 1st order in time 319else 320 crom => crom_eos 321 brom => brom_eos 322endif 323 324! Interpolation of rho^n-1/2 (stored in pcrom) 325! Interpolation of the mass flux at (n+1/2) 326! NB: the mass flux (n+1) is overwritten because not used after. 327! The mass flux for (n->n+1) will be recomputed in cs_pressure_correction 328! FIXME irovar=1 and if dt varies, use theta(rho) = theta(u)*... 329if (vcopt_u%thetav .lt. 1.d0 .and. iappel.eq.1 .and. iterns.gt.1 & 330 .and. itpcol .eq. 0) then 331 allocate(cproa_rho_tc(ncelet)) 332 333 ! Pointer to the previous mass fluxes 334 call field_get_val_prev_s(iflmas, imasfl_prev) 335 call field_get_val_prev_s(iflmab, bmasfl_prev) 336 337 if (irovar.eq.1) then 338 ! remap the density pointer: n-1/2 339 do iel = 1, ncelet 340 cproa_rho_tc(iel) = vcopt_u%thetav * croma(iel) & 341 + (1.d0 - vcopt_u%thetav) * cromaa(iel) 342 enddo 343 pcrom => cproa_rho_tc 344 endif 345 346 ! Inner mass flux interpolation: n-1/2->n+1/2 347 do ifac = 1, nfac 348 imasfl(ifac) = vcopt_u%thetav * imasfl(ifac) & 349 + (1.d0 - vcopt_u%thetav) * imasfl_prev(ifac) 350 enddo 351 352 ! Boundary mass flux interpolation: n-1/2->n+1/2 353 do ifac = 1, nfabor 354 bmasfl(ifac) = vcopt_u%thetav * bmasfl(ifac) & 355 + (1.d0 - vcopt_u%thetav) * bmasfl_prev(ifac) 356 enddo 357 358endif 359 360idftnp = vcopt_u%idften 361imvisp = vcopt_u%imvisf 362 363if (iand(idftnp, ANISOTROPIC_LEFT_DIFFUSION).ne.0) allocate(viscce(6,ncelet)) 364 365! Allocate a temporary array for the prediction-stage error estimator 366if (iescal(iespre).gt.0) then 367 allocate(eswork(3,ncelet)) 368endif 369 370if (iappel.eq.2) then 371 if (iforbr.ge.0 .and. iterns.eq.1 .or. ivofmt.gt.0) then 372 call field_get_val_s(ivarfl(ipr), cvar_pr) 373 endif 374 if(iforbr.ge.0 .and. iterns.eq.1 & 375 .and. (itytur.eq.2 .or. itytur.eq.5 .or. iturb.eq.60) .and. igrhok.eq.1) then 376 call field_get_val_s(ivarfl(ik), cvara_k) 377 endif 378 if (itytur.eq.3.and.iterns.eq.1) then 379 call field_get_val_v(ivarfl(irij), cvara_rij) 380 endif 381else 382 if (iforbr.ge.0 .and. iterns.eq.1 .or. ivofmt.gt.0) then 383 call field_get_val_s(ivarfl(ipr), cvar_pr) 384 endif 385 if(iforbr.ge.0 .and. iterns.eq.1 & 386 .and. (itytur.eq.2 .or. itytur.eq.5 .or. iturb.eq.60) .and. igrhok.eq.1) then 387 call field_get_val_prev_s(ivarfl(ik), cvara_k) 388 endif 389 if (itytur.eq.3.and.iterns.eq.1) then 390 call field_get_val_v(ivarfl(irij), cvara_rij) 391 endif 392endif 393 394if (iforbr.ge.0 .and. iterns.eq.1) call field_get_val_v(iforbr, forbr) 395 396! Theta relatif aux termes sources explicites 397thets = thetsn 398if (isno2t.gt.0) then 399 call field_get_key_int(ivarfl(iu), kstprv, f_id) 400 call field_get_val_v(f_id, c_st_vel) 401else 402 c_st_vel => null() 403endif 404 405!------------------------------------------------------------------------------- 406! ---> Get user source terms 407 408do iel = 1, ncel 409 do isou = 1, 3 410 tsexp(isou,iel) = 0.d0 411 do jsou = 1, 3 412 tsimp(jsou,isou,iel) = 0.d0 413 enddo 414 enddo 415enddo 416 417! The computation of explicit and implicit source terms is performed 418! at the first iteration only. 419! If iphydr=1 or if we have buoyant scalars then we need to update source terms 420call uitsnv (vel, tsexp, tsimp) 421 422call ustsnv & 423 ( nvar , nscal , ncepdp , ncesmp , & 424 iu , & 425 icepdc , icetsm , itypsm , & 426 dt , & 427 ckupdc , smacel , tsexp , tsimp ) 428 429! C version 430call user_source_terms(ivarfl(iu), tsexp, tsimp) 431 432n_fans = cs_fan_n_fans() 433if (n_fans .gt. 0) then 434 if (ntcabs.eq.ntpabs+1) then 435 call debvtl(imasfl, bmasfl, crom, brom) 436 endif 437 call tsvvtl(tsexp) 438endif 439 440! Skip first time step after restart if previous values have not been read. 441if (vcopt_u%ibdtso.lt.0) then 442 vcopt_u%ibdtso = iabs(vcopt_u%ibdtso) 443 call field_set_key_struct_var_cal_opt(ivarfl(iu), vcopt_u) 444endif 445 446! Nudging towards optimal interpolation for velocity 447if (ippmod(iatmos).ge.0) then 448 call field_get_key_int(ivarfl(iu), kopint, f_oi_id) 449 if (f_oi_id.ge.0) then 450 call cs_at_data_assim_source_term(ivarfl(iu), tsexp, tsimp) 451 endif 452 if (iatmst.ge.1) then 453 call cs_at_source_term_for_inlet(tsexp) 454 endif 455endif 456 457! Coupling between two Code_Saturne 458if (nbrcpl.gt.0) then 459 !vectorial interleaved exchange 460 call csccel(iu, vela, coefav, coefbv, tsexp) 461endif 462 463if (vcopt_u%ibdtso.gt.1.and.ntcabs.gt.ntinit & 464 .and.(idtvar.eq.0.or.idtvar.eq.1)) then 465 ! TODO: remove test on ntcabs and implemente a "proper" condition for 466 ! initialization. 467 f_id = ivarfl(iu) 468 call cs_backward_differentiation_in_time(f_id, tsexp, tsimp) 469endif 470 471!=============================================================================== 472! 2. Potential forces (pressure gradient and gravity) 473!=============================================================================== 474 475!------------------------------------------------------------------------------- 476! ---> Pressure gradient 477 478call field_get_id_try("pressure_gradient", f_id) 479if (f_id.ge.0) then 480 call field_get_val_v(f_id, cpro_gradp) 481else 482 allocate(grad(3,ncelet)) 483 cpro_gradp => grad 484endif 485 486iccocg = 1 487inc = 1 488 489! Take the latest pressure field 490iprev = 0 491 492! Namely for the VOF algorithm: consistency of the gradient 493! with the diffusive flux scheme of the correction step 494if (vcopt_p%iwgrec.eq.1) then 495 496 ! retrieve density used in diffusive flux scheme (correction step) 497 if (irovar.eq.1.and.(idilat.gt.1.or.ivofmt.gt.0.or.ippmod(icompf).eq.3)) then 498 call field_get_id("density_mass", f_id) 499 call field_get_val_s(f_id, cpro_rho_mass) 500 501 ! Time interpolated density 502 if (vcopt_u%thetav.lt.1.d0.and.iterns.gt.1) then 503 allocate(cpro_rho_tc(ncelet)) 504 505 do iel = 1, ncelet 506 cpro_rho_tc(iel) = vcopt_u%thetav * cpro_rho_mass(iel) & 507 + (1.d0 - vcopt_u%thetav) * croma(iel) 508 enddo 509 510 wgrec_crom => cpro_rho_tc 511 else 512 wgrec_crom => cpro_rho_mass 513 endif 514 515 ! Weakly variable density algo. (idilat <=1) or constant density 516 else 517 wgrec_crom => crom_eos 518 endif 519 520 ! Id weighting field for gradient 521 call field_get_key_int(ivarfl(ipr), kwgrec, iflwgr) 522 call field_get_dim(iflwgr, f_dim) 523 if (f_dim.gt.1) then 524 call field_get_val_v(iflwgr, cpro_wgrec_v) 525 do iel = 1, ncel 526 ! FIXME should take head losses into account, 527 ! not compatible either with ipucou=1... 528 cpro_wgrec_v(1,iel) = dt(iel) / wgrec_crom(iel) 529 cpro_wgrec_v(2,iel) = dt(iel) / wgrec_crom(iel) 530 cpro_wgrec_v(3,iel) = dt(iel) / wgrec_crom(iel) 531 cpro_wgrec_v(4,iel) = 0.d0 532 cpro_wgrec_v(5,iel) = 0.d0 533 cpro_wgrec_v(6,iel) = 0.d0 534 enddo 535 call syntis(cpro_wgrec_v) 536 537 else 538 call field_get_val_s(iflwgr, cpro_wgrec_s) 539 do iel = 1, ncel 540 cpro_wgrec_s(iel) = dt(iel) / wgrec_crom(iel) 541 enddo 542 call synsca(cpro_wgrec_s) 543 endif 544 if (allocated(cpro_rho_tc)) deallocate(cpro_rho_tc) 545endif 546 547call grdpor(inc) 548 549! Pressure gradient 550call field_gradient_potential(ivarfl(ipr), iprev, 0, inc, & 551 iccocg, iphydr, & 552 frcxt, cpro_gradp) 553 554! Calcul des efforts aux parois (partie 2/5), si demande 555! La pression a la face est calculee comme dans gradrc/gradmc 556! et on la transforme en pression totale 557! On se limite a la premiere iteration (pour faire simple par 558! rapport a la partie issue de condli, hors boucle) 559if (iforbr.ge.0 .and. iterns.eq.1) then 560 call field_get_coefa_s (ivarfl(ipr), coefa_p) 561 call field_get_coefb_s (ivarfl(ipr), coefb_p) 562 do ifac = 1, nfabor 563 iel = ifabor(ifac) 564 diipbx = diipb(1,ifac) 565 diipby = diipb(2,ifac) 566 diipbz = diipb(3,ifac) 567 pip = cvar_pr(iel) & 568 + diipbx*cpro_gradp(1,iel) & 569 + diipby*cpro_gradp(2,iel) & 570 + diipbz*cpro_gradp(3,iel) 571 pfac = coefa_p(ifac) + coefb_p(ifac)*pip 572 pfac = pfac & 573 + ro0*(gx*(cdgfbo(1,ifac)-xyzp0(1)) & 574 + gy*(cdgfbo(2,ifac)-xyzp0(2)) & 575 + gz*(cdgfbo(3,ifac)-xyzp0(3)) ) & 576 - pred0 577 do isou = 1, 3 578 forbr(isou,ifac) = forbr(isou,ifac) + pfac*surfbo(isou,ifac) 579 enddo 580 enddo 581endif 582 583if (iappel.eq.1) then 584 585 ! Initialization 586 ! NB: at the second call, trav contains the temporal increment 587 do iel = 1, ncel 588 trav(1,iel) = 0.d0 589 trav(2,iel) = 0.d0 590 trav(3,iel) = 0.d0 591 enddo 592endif 593 594! FIXME : "rho g" will be second order only if extrapolated 595if (iphydr.eq.1) then 596 do iel = 1, ncel 597 trav(1,iel) = trav(1,iel)+(frcxt(1 ,iel) - cpro_gradp(1,iel)) * cell_f_vol(iel) 598 trav(2,iel) = trav(2,iel)+(frcxt(2 ,iel) - cpro_gradp(2,iel)) * cell_f_vol(iel) 599 trav(3,iel) = trav(3,iel)+(frcxt(3 ,iel) - cpro_gradp(3,iel)) * cell_f_vol(iel) 600 enddo 601else if (ippmod(icompf).ge.0) then 602 do iel = 1, ncel 603 rom = crom(iel) 604 trav(1,iel) = trav(1,iel)+(rom*gx - cpro_gradp(1,iel)) * cell_f_vol(iel) 605 trav(2,iel) = trav(2,iel)+(rom*gy - cpro_gradp(2,iel)) * cell_f_vol(iel) 606 trav(3,iel) = trav(3,iel)+(rom*gz - cpro_gradp(3,iel)) * cell_f_vol(iel) 607 enddo 608 ! Boussinesq approximation 609else if (idilat.eq.0) then 610 611 !FIXME make it dependant on the scalar and use is_buoyant field 612 call field_get_val_s(ibeta, cpro_beta) 613 call field_get_val_s(ivarfl(isca(iscalt)), cvar_t) 614 615 ! Delta rho = - rho_0 beta (T-T0) 616 tref = t0 617 ! for atmospheric flows, variable is potential temperature 618 if (ippmod(iatmos).gt.0) then 619 rscp = rair/cp0 620 tref = t0 * (ps / p0)**rscp 621 endif 622 623 do iel = 1, ncel 624 drom = - crom(iel) * cpro_beta(iel) * (cvar_t(iel) - tref) 625 trav(1,iel) = trav(1,iel) + (drom * gx - cpro_gradp(1,iel)) * cell_f_vol(iel) 626 trav(2,iel) = trav(2,iel) + (drom * gy - cpro_gradp(2,iel)) * cell_f_vol(iel) 627 trav(3,iel) = trav(3,iel) + (drom * gz - cpro_gradp(3,iel)) * cell_f_vol(iel) 628 enddo 629 630else 631 do iel = 1, ncel 632 if (ischtp.eq.2 .and. itpcol.eq.1) then 633 drom = (3.d0/2.d0*croma(iel) - 1.d0/2.d0*cromaa(iel) - ro0) 634 else 635 drom = (crom(iel)-ro0) 636 endif 637 trav(1,iel) = trav(1,iel)+(drom*gx - cpro_gradp(1,iel) ) * cell_f_vol(iel) 638 trav(2,iel) = trav(2,iel)+(drom*gy - cpro_gradp(2,iel) ) * cell_f_vol(iel) 639 trav(3,iel) = trav(3,iel)+(drom*gz - cpro_gradp(3,iel) ) * cell_f_vol(iel) 640 enddo 641endif 642 643! Free memory 644if (allocated(grad)) deallocate(grad) 645 646 647! Pour IAPPEL = 1 (ie appel standard sans les estimateurs) 648! TRAV rassemble les termes sources qui seront recalcules 649! a toutes les iterations sur navsto 650! Si on n'itere pas sur navsto et qu'on n'extrapole pas les 651! termes sources, TRAV contient tous les termes sources 652! jusqu'au basculement dans SMBR 653! A ce niveau, TRAV contient -grad P et rho g 654! P est suppose pris a n+1/2 655! rho est eventuellement interpole a n+1/2 656 657 658!------------------------------------------------------------------------------- 659! ---> Initialize trava array and source terms at the first call (iterns=1) 660 661! trava contains all source terms needed from the first sub iteration 662! (iterns=1) for the other iterations. 663! When there is only one iteration, we build source terms directly in trav 664! array. 665! Explicit source terms will be used at the next time step in case of 666! extrapolation (if there is only one or many iteration on navtsv) 667 668! At the first iteration on navstv 669if (iterns.eq.1) then 670 671 ! Si on extrapole les T.S. : -theta*valeur precedente 672 if (isno2t.gt.0) then 673 ! S'il n'y a qu'une iter : TRAV incremente 674 if (nterup.eq.1) then 675 do iel = 1, ncel 676 do ii = 1, ndim 677 trav (ii,iel) = trav (ii,iel) - thets*c_st_vel(ii,iel) 678 enddo 679 enddo 680 ! S'il y a plusieurs iter : TRAVA initialise 681 else 682 do iel = 1, ncel 683 do ii = 1, ndim 684 trava(ii,iel) = - thets*c_st_vel(ii,iel) 685 enddo 686 enddo 687 endif 688 ! Et on initialise le terme source pour le remplir ensuite 689 do iel = 1, ncel 690 do ii = 1, ndim 691 c_st_vel(ii,iel) = 0.d0 692 enddo 693 enddo 694 695 ! Si on n'extrapole pas les T.S. 696 else 697 ! S'il y a plusieurs iter : TRAVA initialise 698 ! sinon TRAVA n'existe pas 699 if(nterup.gt.1) then 700 do ii = 1, ndim 701 do iel = 1, ncel 702 trava(ii,iel) = 0.d0 703 enddo 704 enddo 705 endif 706 endif 707 708endif 709 710!------------------------------------------------------------------------------- 711! Initialization of the implicit terms 712 713if (iappel.eq.1) then 714 715 do iel = 1, ncel 716 do isou = 1, 3 717 fimp(isou,isou,iel) = vcopt_u%istat*pcrom(iel)/dt(iel)*cell_f_vol(iel) 718 do jsou = 1, 3 719 if(jsou.ne.isou) fimp(jsou,isou,iel) = 0.d0 720 enddo 721 enddo 722 enddo 723 724! Le remplissage de FIMP est toujours indispensable, 725! meme si on peut se contenter de n'importe quoi pour IAPPEL=2. 726else 727 do iel = 1, ncel 728 do isou = 1, 3 729 do jsou = 1, 3 730 fimp(jsou,isou,iel) = 0.d0 731 enddo 732 enddo 733 enddo 734endif 735 736!------------------------------------------------------------------------------- 737! ---> 2/3 RHO * GRADIENT DE K SI k-epsilon ou k-omega 738! NB : ON NE PREND PAS LE GRADIENT DE (RHO K), MAIS 739! CA COMPLIQUERAIT LA GESTION DES CL ... 740! On peut se demander si l'extrapolation en temps sert a 741! quelquechose 742 743! Ce terme explicite est calcule une seule fois, 744! a la premiere iter sur navsto : il est stocke dans un champ si on 745! doit l'extrapoler en temps ; il va dans TRAVA si on n'extrapole 746! pas mais qu'on itere sur navsto. Il va dans TRAV si on 747! n'extrapole pas et qu'on n'itere pas sur navsto. 748if( (itytur.eq.2 .or. itytur.eq.5 .or. iturb.eq.60) & 749 .and. igrhok.eq.1 .and. iterns.eq.1) then 750 751 ! Allocate a work array for the gradient calculation 752 allocate(grad(3,ncelet)) 753 754 iccocg = 1 755 iprev = 1 756 inc = 1 757 758 call field_gradient_scalar(ivarfl(ik), iprev, 0, inc, iccocg, grad) 759 760 d2s3 = 2.d0/3.d0 761 762 ! Si on extrapole les termes source en temps 763 if (isno2t.gt.0) then 764 ! Calcul de rho^n grad k^n si rho non extrapole 765 ! rho^n grad k^n si rho extrapole 766 767 do iel = 1, ncel 768 romvom = -croma(iel)*cell_f_vol(iel)*d2s3 769 do isou = 1, 3 770 c_st_vel(isou,iel) = c_st_vel(isou,iel)+grad(isou,iel)*romvom 771 enddo 772 enddo 773 ! Si on n'extrapole pas les termes sources en temps : TRAV ou TRAVA 774 else 775 if(nterup.eq.1) then 776 do iel = 1, ncel 777 romvom = -crom(iel)*cell_f_vol(iel)*d2s3 778 do isou = 1, 3 779 trav(isou,iel) = trav(isou,iel) + grad(isou,iel) * romvom 780 enddo 781 enddo 782 else 783 do iel = 1, ncel 784 romvom = -crom(iel)*cell_f_vol(iel)*d2s3 785 do isou = 1, 3 786 trava(isou,iel) = trava(isou,iel) + grad(isou,iel) * romvom 787 enddo 788 enddo 789 endif 790 endif 791 792 ! Calcul des efforts aux parois (partie 3/5), si demande 793 if (iforbr.ge.0) then 794 call field_get_coefa_s (ivarfl(ik), coefa_k) 795 call field_get_coefb_s (ivarfl(ik), coefb_k) 796 do ifac = 1, nfabor 797 iel = ifabor(ifac) 798 diipbx = diipb(1,ifac) 799 diipby = diipb(2,ifac) 800 diipbz = diipb(3,ifac) 801 xkb = cvara_k(iel) + diipbx*grad(1,iel) & 802 + diipby*grad(2,iel) + diipbz*grad(3,iel) 803 xkb = coefa_k(ifac)+coefb_k(ifac)*xkb 804 xkb = d2s3*crom(iel)*xkb 805 do isou = 1, 3 806 forbr(isou,ifac) = forbr(isou,ifac) + xkb*surfbo(isou,ifac) 807 enddo 808 enddo 809 endif 810 811 ! Free memory 812 deallocate(grad) 813 814endif 815 816 817!------------------------------------------------------------------------------- 818! ---> Transpose of velocity gradient in the diffusion term 819 820! These terms are taken into account in cs_balance_vector. 821! We only compute here the secondary viscosity. 822 823if (ivisse.eq.1) then 824 825 call visecv(secvif, secvib) 826 827endif 828 829!------------------------------------------------------------------------------- 830! ---> Head losses 831! (if iphydr=1 this term has already been taken into account) 832 833! ---> Explicit part 834if ((ncepdp.gt.0).and.(iphydr.ne.1)) then 835 836 ! Les termes diagonaux sont places dans TRAV ou TRAVA, 837 ! La prise en compte de velk a partir de la seconde iteration 838 ! est faite directement dans cs_equation_iterative_solve_vector. 839 if (iterns.eq.1) then 840 841 allocate(hl_exp(3, ncepdp)) 842 843 call tspdcv(ncepdp, icepdc, vela, ckupdc, hl_exp) 844 845 ! If we have inner iterations, we use trava, otherwise trav 846 if(nterup.gt.1) then 847 do ielpdc = 1, ncepdp 848 iel = icepdc(ielpdc) 849 trava(1,iel) = trava(1,iel) + hl_exp(1,ielpdc) 850 trava(2,iel) = trava(2,iel) + hl_exp(2,ielpdc) 851 trava(3,iel) = trava(3,iel) + hl_exp(3,ielpdc) 852 enddo 853 else 854 do ielpdc = 1, ncepdp 855 iel = icepdc(ielpdc) 856 trav(1,iel) = trav(1,iel) + hl_exp(1,ielpdc) 857 trav(2,iel) = trav(2,iel) + hl_exp(2,ielpdc) 858 trav(3,iel) = trav(3,iel) + hl_exp(3,ielpdc) 859 enddo 860 endif 861 862 deallocate(hl_exp) 863 endif 864 865endif 866 867! ---> Implicit part 868 869! At the second call, fimp is not needed anymore 870if (iappel.eq.1) then 871 if (ncepdp.gt.0) then 872 ! The theta-scheme for the head loss is the same as the other terms 873 thetap = vcopt_u%thetav 874 do ielpdc = 1, ncepdp 875 iel = icepdc(ielpdc) 876 romvom = crom(iel)*cell_f_vol(iel)*thetap 877 878 ! Diagonal part 879 do isou = 1, 3 880 fimp(isou,isou,iel) = fimp(isou,isou,iel) + romvom*ckupdc(isou,ielpdc) 881 enddo 882 ! Extra-diagonal part 883 cpdc12 = ckupdc(4,ielpdc) 884 cpdc23 = ckupdc(5,ielpdc) 885 cpdc13 = ckupdc(6,ielpdc) 886 887 fimp(1,2,iel) = fimp(1,2,iel) + romvom*cpdc12 888 fimp(2,1,iel) = fimp(2,1,iel) + romvom*cpdc12 889 fimp(1,3,iel) = fimp(1,3,iel) + romvom*cpdc13 890 fimp(3,1,iel) = fimp(3,1,iel) + romvom*cpdc13 891 fimp(2,3,iel) = fimp(2,3,iel) + romvom*cpdc23 892 fimp(3,2,iel) = fimp(3,2,iel) + romvom*cpdc23 893 enddo 894 endif 895endif 896 897 898!------------------------------------------------------------------------------- 899! ---> Coriolis force 900! (if iphydr=1 then this term is already taken into account) 901 902! ---> Explicit part 903 904if ((icorio.eq.1.or.iturbo.eq.1) .and. iphydr.ne.1) then 905 906 ! A la premiere iter sur navsto, on ajoute la partie issue des 907 ! termes explicites 908 if (iterns.eq.1) then 909 910 ! Si on n'itere pas sur navsto : TRAV 911 if (nterup.eq.1) then 912 913 ! Reference frame rotation 914 do iel = 1, ncel 915 romvom = -2.d0*crom(iel)*cell_f_vol(iel) 916 call add_coriolis_v(0, romvom, vela(:,iel), trav(:,iel)) 917 enddo 918 ! Turbomachinery frozen rotors rotation 919 if (iturbo.eq.1) then 920 do iel = 1, ncel 921 if (irotce(iel).gt.0) then 922 romvom = -crom(iel)*cell_f_vol(iel) 923 call add_coriolis_v(irotce(iel), romvom, vela(:,iel), trav(:,iel)) 924 endif 925 enddo 926 endif 927 928 ! Si on itere sur navsto : TRAVA 929 else 930 931 ! Reference frame rotation 932 do iel = 1, ncel 933 romvom = -2.d0*crom(iel)*cell_f_vol(iel) 934 call add_coriolis_v(0, romvom, vela(:,iel), trava(:,iel)) 935 enddo 936 ! Turbomachinery frozen rotors rotation 937 if (iturbo.eq.1) then 938 do iel = 1, ncel 939 if (irotce(iel).gt.0) then 940 romvom = -crom(iel)*cell_f_vol(iel) 941 call add_coriolis_v(irotce(iel), romvom, vela(:,iel), trava(:,iel)) 942 endif 943 enddo 944 endif 945 946 endif 947 endif 948endif 949 950! ---> Implicit part 951 952! At the second call, fimp is not needed anymore 953if (iappel.eq.1) then 954 if (icorio.eq.1 .or. iturbo.eq.1) then 955 ! The theta-scheme for the Coriolis term is the same as the other terms 956 thetap = vcopt_u%thetav 957 958 ! Reference frame rotation 959 do iel = 1, ncel 960 romvom = -2.d0*crom(iel)*cell_f_vol(iel)*thetap 961 call add_coriolis_t(0, romvom, fimp(:,:,iel)) 962 enddo 963 ! Turbomachinery frozen rotors rotation 964 if (iturbo.eq.1) then 965 do iel = 1, ncel 966 if (irotce(iel).gt.0) then 967 romvom = -crom(iel)*cell_f_vol(iel)*thetap 968 call add_coriolis_t(irotce(iel), romvom, fimp(:,:,iel)) 969 endif 970 enddo 971 endif 972 endif 973endif 974 975!------------------------------------------------------------------------------- 976! ---> - Divergence of tensor Rij 977! ---> - Non linear part of Rij for non-liear Eddy Viscosity Models 978 979if((itytur.eq.3.or.iturb.eq.23).and.iterns.eq.1) then 980 981 allocate(rij(6,ncelet)) 982 allocate(coefat(6,nfabor)) 983 allocate(coefbt(6,6,nfabor)) 984 985 ! Reynolds Stress Models 986 if (itytur.eq.3) then 987 988 do iel = 1, ncelet 989 rij(1,iel) = cvara_rij(1,iel) 990 rij(2,iel) = cvara_rij(2,iel) 991 rij(3,iel) = cvara_rij(3,iel) 992 rij(4,iel) = cvara_rij(4,iel) 993 rij(5,iel) = cvara_rij(5,iel) 994 rij(6,iel) = cvara_rij(6,iel) 995 enddo 996 ! --- Boundary conditions on the components of the tensor Rij 997 call field_get_coefad_v(ivarfl(irij),coefap) 998 coefat = coefap 999 1000 do ifac = 1, nfabor 1001 do ii = 1, 6 1002 do jj = 1, 6 1003 coefbt(jj,ii,ifac) = 0.d0 1004 enddo 1005 enddo 1006 enddo 1007 1008 call field_get_coefbd_v(ivarfl(irij),coefbp) 1009 coefbt = coefbp 1010 1011 ! Baglietto et al. quadratic k-epislon model 1012 else if(iturb.eq.23) then 1013 1014 ! --- Compute the non linear part of Rij 1015 call cnlevm (rij) 1016 1017 ! --- Boundary conditions : Homogeneous Neumann 1018 do ifac = 1, nfabor 1019 do ii = 1, 6 1020 coefat(ii,ifac) = 0.d0 1021 do jj = 1, 6 1022 coefbt(jj,ii,ifac) = 1.d0 1023 enddo 1024 enddo 1025 enddo 1026 1027 end if 1028 1029 ! Flux computation options 1030 f_id = -1 1031 init = 1 1032 inc = 1 1033 iflmb0 = 0 1034 if(itytur.eq.3) then 1035 call field_get_key_struct_var_cal_opt(ivarfl(irij), vcopt) 1036 else 1037 call field_get_key_struct_var_cal_opt(ivarfl(ik), vcopt) 1038 end if 1039 imrgrp = vcopt%imrgra 1040 nswrgp = vcopt%nswrgr 1041 imligp = vcopt%imligr 1042 iwarnp = vcopt%iwarni 1043 epsrgp = vcopt%epsrgr 1044 climgp = vcopt%climgr 1045 itypfl = 1 1046 1047 allocate(tflmas(3,nfac)) 1048 allocate(tflmab(3,nfabor)) 1049 1050 call divrij & 1051 ( f_id , itypfl , & 1052 iflmb0 , init , inc , imrgrp , nswrgp , imligp , & 1053 iwarnp , & 1054 epsrgp , climgp , & 1055 crom , brom , & 1056 rij , & 1057 coefat , coefbt , & 1058 tflmas , tflmab ) 1059 1060 deallocate(rij) 1061 deallocate(coefat, coefbt) 1062 1063 ! Calcul des efforts aux bords (partie 5/5), si necessaire 1064 1065 if (iforbr.ge.0) then 1066 do ifac = 1, nfabor 1067 do isou = 1, 3 1068 forbr(isou,ifac) = forbr(isou,ifac) + tflmab(isou,ifac) 1069 enddo 1070 enddo 1071 endif 1072 1073 1074 call field_get_id_try("reynolds_stress_divergence", f_id) 1075 if (f_id.ge.0) then 1076 call field_get_val_v(f_id, cpro_divr) 1077 else 1078 allocate(divt(3,ncelet)) 1079 cpro_divr => divt 1080 endif 1081 1082 init = 1 1083 call divmat(init,tflmas,tflmab, cpro_divr) 1084 1085 deallocate(tflmas, tflmab) 1086 1087 ! (if iphydr=1 then this term is already taken into account) 1088 if (iphydr.ne.1.or.igprij.ne.1) then 1089 1090 ! If extrapolation of source terms 1091 if (isno2t.gt.0) then 1092 do iel = 1, ncel 1093 do isou = 1, 3 1094 c_st_vel(isou,iel) = c_st_vel(isou,iel) - cpro_divr(isou,iel) 1095 enddo 1096 enddo 1097 1098 ! No extrapolation of source terms 1099 else 1100 1101 ! No inner iteration 1102 if (nterup.eq.1) then 1103 do iel = 1, ncel 1104 do isou = 1, 3 1105 trav(isou,iel) = trav(isou,iel) - cpro_divr(isou,iel) 1106 enddo 1107 enddo 1108 ! Inner iterations 1109 else 1110 do iel = 1, ncel 1111 do isou = 1, 3 1112 trava(isou,iel) = trava(isou,iel) - cpro_divr(isou,iel) 1113 enddo 1114 enddo 1115 endif 1116 endif 1117 endif 1118 1119endif 1120 1121!------------------------------------------------------------------------------- 1122! ---> Face diffusivity for the velocity 1123 1124if (vcopt_u%idiff.ge. 1) then 1125 1126 call field_get_val_s(iviscl, viscl) 1127 call field_get_val_s(ivisct, visct) 1128 1129 if (itytur.eq.3) then 1130 do iel = 1, ncel 1131 w1(iel) = viscl(iel) 1132 enddo 1133 else 1134 do iel = 1, ncel 1135 w1(iel) = viscl(iel) + vcopt_u%idifft*visct(iel) 1136 enddo 1137 endif 1138 1139 ! Scalar diffusivity (Default) 1140 if (iand(idftnp, ISOTROPIC_DIFFUSION).ne.0) then 1141 1142 call viscfa & 1143 ( imvisp , & 1144 w1 , & 1145 viscf , viscb ) 1146 1147 ! When using Rij-epsilon model with the option irijnu=1, the face 1148 ! viscosity for the Matrix (viscfi and viscbi) is increased 1149 if(itytur.eq.3.and.irijnu.eq.1) then 1150 1151 do iel = 1, ncel 1152 w1(iel) = viscl(iel) + vcopt_u%idifft*visct(iel) 1153 enddo 1154 1155 call viscfa & 1156 ( imvisp , & 1157 w1 , & 1158 viscfi , viscbi ) 1159 endif 1160 1161 ! Tensorial diffusion of the velocity (in case of tensorial porosity) 1162 else if (iand(idftnp, ANISOTROPIC_LEFT_DIFFUSION).ne.0) then 1163 1164 do iel = 1, ncel 1165 do isou = 1, 3 1166 viscce(isou, iel) = w1(iel) 1167 enddo 1168 do isou = 4, 6 1169 viscce(isou, iel) = 0.d0 1170 enddo 1171 enddo 1172 1173 call vistnv & 1174 ( imvisp , & 1175 viscce , & 1176 viscf , viscb ) 1177 1178 ! When using Rij-epsilon model with the option irijnu=1, the face 1179 ! viscosity for the Matrix (viscfi and viscbi) is increased 1180 if(itytur.eq.3.and.irijnu.eq.1) then 1181 1182 do iel = 1, ncel 1183 w1(iel) = viscl(iel) + vcopt_u%idifft*visct(iel) 1184 enddo 1185 1186 do iel = 1, ncel 1187 do isou = 1, 3 1188 viscce(isou, iel) = w1(iel) 1189 enddo 1190 do isou = 4, 6 1191 viscce(isou, iel) = 0.d0 1192 enddo 1193 enddo 1194 1195 call vistnv & 1196 ( imvisp , & 1197 viscce , & 1198 viscfi , viscbi ) 1199 1200 endif 1201 endif 1202 1203! --- If no diffusion, viscosity is set to 0. 1204else 1205 1206 do ifac = 1, nfac 1207 viscf(ifac) = 0.d0 1208 enddo 1209 do ifac = 1, nfabor 1210 viscb(ifac) = 0.d0 1211 enddo 1212 1213 if(itytur.eq.3.and.irijnu.eq.1) then 1214 do ifac = 1, nfac 1215 viscfi(ifac) = 0.d0 1216 enddo 1217 do ifac = 1, nfabor 1218 viscbi(ifac) = 0.d0 1219 enddo 1220 endif 1221 1222endif 1223 1224!------------------------------------------------------------------------------- 1225! ---> Take external forces partially equilibrated with the pressure gradient 1226! into account (only for the first call, the second one is dedicated 1227! to error estimators) 1228 1229if (iappel.eq.1.and.iphydr.eq.1) then 1230 1231 ! External forces at previous time step: 1232 ! frcxt was initialised to 0 1233 ! NB: frcxt was used in typecl, and will be updated 1234 ! at the end of navstv 1235 1236 ! External force variation between time step n and n+1 1237 ! (used in the correction step) 1238 !----------------------------- 1239 1240 ! Boussinesq approximation 1241 if (idilat.eq.0) then 1242 1243 !FIXME make it dependant on the scalar and use is_buoyant field 1244 call field_get_val_s(ibeta, cpro_beta) 1245 call field_get_val_s(ivarfl(isca(iscalt)), cvar_t) 1246 1247 ! Delta rho = - rho_0 beta (T-T0) 1248 tref = t0 1249 ! for atmospheric flows, variable is potential temperature 1250 if (ippmod(iatmos).gt.0) then 1251 rscp = rair/cp0 1252 tref = t0 * (ps / p0)**rscp 1253 endif 1254 1255 do iel = 1, ncel 1256 drom = - crom(iel) * cpro_beta(iel) * (cvar_t(iel) - tref) * cell_is_active(iel) 1257 dfrcxt(1, iel) = drom*gx - frcxt(1, iel) * cell_is_active(iel) 1258 dfrcxt(2, iel) = drom*gy - frcxt(2, iel) * cell_is_active(iel) 1259 dfrcxt(3, iel) = drom*gz - frcxt(3, iel) * cell_is_active(iel) 1260 enddo 1261 1262 else 1263 do iel = 1, ncel 1264 if (ischtp.eq.2 .and. itpcol.eq.1) then 1265 drom = (3.d0/2.d0*croma(iel) - 1.d0/2.d0*cromaa(iel) - ro0) & 1266 * cell_is_active(iel) 1267 else 1268 drom = (crom(iel)-ro0) * cell_is_active(iel) 1269 endif 1270 1271 dfrcxt(1, iel) = drom*gx - frcxt(1, iel) * cell_is_active(iel) 1272 dfrcxt(2, iel) = drom*gy - frcxt(2, iel) * cell_is_active(iel) 1273 dfrcxt(3, iel) = drom*gz - frcxt(3, iel) * cell_is_active(iel) 1274 enddo 1275 endif 1276 1277 ! Add head losses 1278 if (ncepdp.gt.0) then 1279 do ielpdc = 1, ncepdp 1280 iel=icepdc(ielpdc) 1281 vit1 = vela(1,iel) * cell_is_active(iel) 1282 vit2 = vela(2,iel) * cell_is_active(iel) 1283 vit3 = vela(3,iel) * cell_is_active(iel) 1284 cpdc11 = ckupdc(1,ielpdc) 1285 cpdc22 = ckupdc(2,ielpdc) 1286 cpdc33 = ckupdc(3,ielpdc) 1287 cpdc12 = ckupdc(4,ielpdc) 1288 cpdc23 = ckupdc(5,ielpdc) 1289 cpdc13 = ckupdc(6,ielpdc) 1290 dfrcxt(1 ,iel) = dfrcxt(1 ,iel) & 1291 - crom(iel)*(cpdc11*vit1+cpdc12*vit2+cpdc13*vit3) 1292 dfrcxt(2 ,iel) = dfrcxt(2 ,iel) & 1293 - crom(iel)*(cpdc12*vit1+cpdc22*vit2+cpdc23*vit3) 1294 dfrcxt(3 ,iel) = dfrcxt(3 ,iel) & 1295 - crom(iel)*(cpdc13*vit1+cpdc23*vit2+cpdc33*vit3) 1296 enddo 1297 endif 1298 1299 ! Add Coriolis force 1300 if (icorio.eq.1 .or. iturbo.eq.1) then 1301 1302 ! Reference frame rotation 1303 do iel = 1, ncel 1304 rom = -2.d0*crom(iel) * cell_is_active(iel) 1305 call add_coriolis_v(0, rom, vela(:,iel), dfrcxt(:,iel)) 1306 enddo 1307 ! Turbomachinery frozen rotors rotation 1308 if (iturbo.eq.1) then 1309 do iel = 1, ncel 1310 if (irotce(iel).gt.0) then 1311 rom = -crom(iel) * cell_is_active(iel) 1312 call add_coriolis_v(irotce(iel), rom, vela(:,iel), dfrcxt(:,iel)) 1313 endif 1314 enddo 1315 endif 1316 endif 1317 1318 ! Add -div( rho R) as external force 1319 if (itytur.eq.3.and.igprij.eq.1) then 1320 do iel = 1, ncel 1321 dvol = 0.d0 1322 ! If it is not a solid cell 1323 if (cell_is_active(iel).eq.1) dvol = 1.d0 / cell_f_vol(iel) 1324 do isou = 1, 3 1325 dfrcxt(isou, iel) = dfrcxt(isou, iel) - cpro_divr(isou, iel)*dvol 1326 enddo 1327 enddo 1328 endif 1329 1330 ! ---> Use user source terms 1331 1332 if (igpust.eq.1) then 1333 do iel = 1, ncel 1334 dvol = 0.d0 1335 ! If it is not a solid cell 1336 if (cell_is_active(iel).eq.1) dvol = 1.d0 / cell_f_vol(iel) 1337 1338 do isou = 1, 3 1339 dfrcxt(isou, iel) = dfrcxt(isou, iel) + tsexp(isou, iel)*dvol 1340 enddo 1341 enddo 1342 1343 endif 1344 1345 if (irangp.ge.0.or.iperio.eq.1) then 1346 call synvin(dfrcxt) 1347 endif 1348 1349endif 1350 1351!=============================================================================== 1352! 3. Solving of the 3x3xNcel coupled system 1353!=============================================================================== 1354 1355 1356! ---> AU PREMIER APPEL, 1357! MISE A ZERO DE L'ESTIMATEUR POUR LA VITESSE PREDITE 1358! S'IL DOIT ETRE CALCULE 1359 1360if (iappel.eq.1) then 1361 if (iestim(iespre).ge.0) then 1362 call field_get_val_s(iestim(iespre), c_estim) 1363 do iel = 1, ncel 1364 c_estim(iel) = 0.d0 1365 enddo 1366 endif 1367endif 1368 1369! ---> AU DEUXIEME APPEL, 1370! MISE A ZERO DE L'ESTIMATEUR TOTAL POUR NAVIER-STOKES 1371! (SI ON FAIT UN DEUXIEME APPEL, ALORS IL DOIT ETRE CALCULE) 1372 1373if (iappel.eq.2) then 1374 call field_get_val_s(iestim(iestot), c_estim) 1375 do iel = 1, ncel 1376 c_estim(iel) = 0.d0 1377 enddo 1378endif 1379 1380!------------------------------------------------------------------------------- 1381! ---> Use user source terms 1382 1383! ---> Explicit contribution due to implicit terms 1384 1385if (iterns.eq.1) then 1386 if (nterup.gt.1) then 1387 do iel = 1, ncel 1388 do isou = 1, 3 1389 do jsou = 1, 3 1390 trava(isou,iel) = trava(isou,iel) & 1391 + tsimp(jsou,isou,iel)*vela(jsou,iel) 1392 enddo 1393 enddo 1394 enddo 1395 else 1396 do iel = 1, ncel 1397 do isou = 1, 3 1398 do jsou = 1, 3 1399 trav(isou,iel) = trav(isou,iel) & 1400 + tsimp(jsou,isou,iel)*vela(jsou,iel) 1401 enddo 1402 enddo 1403 enddo 1404 endif 1405endif 1406 1407! Explicit user source terms are added 1408if ((iphydr.ne.1.or.igpust.ne.1)) then 1409 ! If source terms are time-extrapolated, they are stored in fields 1410 if (isno2t.gt.0) then 1411 if (iterns.eq.1) then 1412 do iel = 1, ncel 1413 do isou = 1, 3 1414 c_st_vel(isou,iel) = c_st_vel(isou,iel) + tsexp(isou,iel) 1415 enddo 1416 enddo 1417 endif 1418 1419 else 1420 ! Alwways in the current work array because this may be updated 1421 ! during inner iterations 1422 do iel = 1, ncel 1423 do isou = 1, 3 1424 trav(isou,iel) = trav(isou,iel) + tsexp(isou,iel) 1425 enddo 1426 enddo 1427 endif 1428endif 1429 1430! ---> Implicit terms 1431if (iappel.eq.1) then 1432 ! If source terms are time-extrapolated 1433 if (isno2t.gt.0) then 1434 thetap = vcopt_u%thetav 1435 do iel = 1, ncel 1436 do isou = 1, 3 1437 do jsou = 1, 3 1438 fimp(jsou,isou,iel) = fimp(jsou,isou,iel) & 1439 - tsimp(jsou,isou,iel)*thetap 1440 enddo 1441 enddo 1442 enddo 1443 else 1444 do iel = 1, ncel 1445 do isou = 1, 3 1446 do jsou = 1, 3 1447 fimp(jsou,isou,iel) = fimp(jsou,isou,iel) & 1448 + max(-tsimp(jsou,isou,iel),zero) 1449 enddo 1450 enddo 1451 enddo 1452 endif 1453endif 1454 1455!------------------------------------------------------------------------------- 1456! ---> Mass source terms 1457 1458if (ncesmp.gt.0) then 1459 1460! On calcule les termes Gamma (uinj - u) 1461! -Gamma u a la premiere iteration est mis dans 1462! TRAV ou TRAVA selon qu'on itere ou non sur navsto 1463! Gamma uinj a la premiere iteration est placee dans W1 1464! ROVSDT a chaque iteration recoit Gamma 1465 allocate(gavinj(3,ncelet)) 1466 if (nterup.eq.1) then 1467 call catsmv(ncesmp, iterns, icetsm, itypsm(1,iu), & 1468 cell_f_vol, vela, smacel(:,iu), smacel(:,ipr), & 1469 trav, fimp, gavinj) 1470 else 1471 call catsmv(ncesmp, iterns, icetsm, itypsm(1,iu), & 1472 cell_f_vol, vela, smacel(:,iu), smacel(:,ipr), & 1473 trava, fimp, gavinj) 1474 endif 1475 1476 ! At the first inner iteration, the explicit part "Gamma u^{in}" is added 1477 if (iterns.eq.1) then 1478 ! If source terms are extrapolated, stored in fields 1479 if(isno2t.gt.0) then 1480 do iel = 1, ncel 1481 do isou = 1, 3 1482 c_st_vel(isou,iel) = c_st_vel(isou,iel) + gavinj(isou,iel) 1483 enddo 1484 enddo 1485 1486 else 1487 ! If no inner iteration: in trav 1488 if (nterup.eq.1) then 1489 do iel = 1,ncel 1490 do isou = 1, 3 1491 trav(isou,iel) = trav(isou,iel) + gavinj(isou,iel) 1492 enddo 1493 enddo 1494 ! Otherwise, in trava 1495 else 1496 do iel = 1,ncel 1497 do isou = 1, 3 1498 trava(isou,iel) = trava(isou,iel) + gavinj(isou,iel) 1499 enddo 1500 enddo 1501 endif 1502 endif 1503 endif 1504 1505 deallocate(gavinj) 1506 1507endif 1508 1509! ---> Right Hand Side initialization 1510 1511! If source terms are extrapolated in time 1512if (isno2t.gt.0) then 1513 thetp1 = 1.d0 + thets 1514 ! If no inner iteration: trav 1515 if (nterup.eq.1) then 1516 do iel = 1, ncel 1517 do isou = 1, 3 1518 smbr(isou,iel) = trav(isou,iel) + thetp1*c_st_vel(isou,iel) 1519 enddo 1520 enddo 1521 1522 else 1523 do iel = 1, ncel 1524 do isou = 1, 3 1525 smbr(isou,iel) = trav(isou,iel) + trava(isou,iel) & 1526 + thetp1*c_st_vel(isou,iel) 1527 enddo 1528 enddo 1529 endif 1530 1531! No time extrapolation 1532else 1533 ! No inner iteration 1534 if (nterup.eq.1) then 1535 do iel = 1, ncel 1536 do isou = 1, 3 1537 smbr(isou,iel) = trav(isou,iel) 1538 enddo 1539 enddo 1540 ! Inner iterations 1541 else 1542 do iel = 1, ncel 1543 do isou = 1, 3 1544 smbr(isou,iel) = trav(isou,iel) + trava(isou,iel) 1545 enddo 1546 enddo 1547 endif 1548endif 1549 1550! ---> LAGRANGIEN : COUPLAGE RETOUR 1551 1552! L'ordre 2 sur les termes issus du lagrangien necessiterait de 1553! decomposer TSLAGR(IEL,ISOU) en partie implicite et 1554! explicite, comme c'est fait dans ustsnv. 1555! Pour le moment, on n'y touche pas. 1556 1557if (iilagr.eq.2 .and. ltsdyn.eq.1) then 1558 1559 call field_get_val_v_by_name('velocity_st_lagr', lagr_st_vel) 1560 1561 do iel = 1, ncel 1562 do isou = 1, 3 1563 smbr(isou,iel) = smbr(isou,iel) + lagr_st_vel(isou,iel) 1564 enddo 1565 enddo 1566 ! At the second call, fimp is unused 1567 if(iappel.eq.1) then 1568 do iel = 1, ncel 1569 do isou = 1, 3 1570 fimp(isou,isou,iel) = fimp(isou,isou,iel) + max(-tslagr(iel,itsli),zero) 1571 enddo 1572 enddo 1573 endif 1574 1575endif 1576 1577! ---> Electric Arc (Laplace Force) 1578! (No 2nd order in time yet) 1579if (ippmod(ielarc).ge.1) then 1580 call field_get_val_v_by_name('laplace_force', lapla) 1581 do iel = 1, ncel 1582 smbr(1,iel) = smbr(1,iel) + cell_f_vol(iel) * lapla(1,iel) 1583 smbr(2,iel) = smbr(2,iel) + cell_f_vol(iel) * lapla(2,iel) 1584 smbr(3,iel) = smbr(3,iel) + cell_f_vol(iel) * lapla(3,iel) 1585 enddo 1586endif 1587 1588! Solver parameters 1589 1590if (ippmod(icompf).ge.0) then 1591 ! impose boundary convective flux at some faces (face indicator icvfli) 1592 icvflb = 1 1593else 1594 ! all boundary convective flux with upwind 1595 icvflb = 0 1596endif 1597 1598if (staggered.eq.1) then 1599 do iel = 1, ncel 1600 smbr(1,iel) = 0 1601 smbr(2,iel) = 0 1602 smbr(3,iel) = 0 1603 enddo 1604endif 1605 1606if (iappel.eq.1) then 1607 1608 ! Store fimp as the velocity matrix is stored in it in codtiv call 1609 do iel = 1, ncel 1610 do isou = 1, 3 1611 do jsou = 1, 3 1612 fimpcp(jsou,isou,iel) = fimp(jsou,isou,iel) 1613 enddo 1614 enddo 1615 enddo 1616 1617 iescap = iescal(iespre) 1618 1619 vcopt_loc = vcopt_u 1620 1621 vcopt_loc%istat = -1 1622 vcopt_loc%idifft = -1 1623 vcopt_loc%iwgrec = 0 1624 vcopt_loc%blend_st = 0 ! Warning, may be overwritten if a field 1625 1626 p_k_value => vcopt_loc 1627 c_k_value = equation_param_from_vcopt(c_loc(p_k_value)) 1628 1629 ! Warning: in case of convergence estimators, eswork give the estimator 1630 ! of the predicted velocity 1631 call cs_equation_iterative_solve_vector & 1632 ( idtvar , iterns , & 1633 ivarfl(iu) , c_null_char , & 1634 ivisse , iescap , c_k_value , & 1635 vela , velk , & 1636 coefav , coefbv , cofafv , cofbfv , & 1637 imasfl , bmasfl , viscfi , & 1638 viscbi , viscf , viscb , secvif , & 1639 secvib , rvoid , rvoid , rvoid , & 1640 icvflb , icvfli , & 1641 fimp , smbr , vel , eswork ) 1642 1643 ! Store inverse of the velocity matrix for the correction step 1644 ! if needed (otherwise vitenp is used in cs_pressure_correction) 1645 if (rcfact.eq.1) then 1646 1647 do iel = 1, ncel 1648 1649 tensor(1) = fimp(1,1,iel)/crom(iel) 1650 tensor(2) = fimp(2,2,iel)/crom(iel) 1651 tensor(3) = fimp(3,3,iel)/crom(iel) 1652 tensor(4) = fimp(1,2,iel)/crom(iel) 1653 tensor(5) = fimp(2,3,iel)/crom(iel) 1654 tensor(6) = fimp(1,3,iel)/crom(iel) 1655 1656 call symmetric_matrix_inverse(tensor, da_uu(:, iel)) 1657 1658 do ii = 1, 6 1659 da_uu(ii,iel) = cell_f_vol(iel)*da_uu(ii,iel) 1660 enddo 1661 1662 enddo 1663 1664 call syntis(da_uu) 1665 1666 endif 1667 1668 1669 ! Velocity-pression coupling: compute the vector T, stored in tpucou, 1670 ! cs_equation_iterative_solve_vector is called, only one sweep is done, 1671 ! and tpucou is initialized by 0, so that the advection/diffusion added 1672 ! by cs_balance_vector is 0. 1673 ! nswrsp = -1 indicated that only one sweep is required and inc=0 1674 ! for boundary contitions on the weight matrix. 1675 if (ipucou.eq.1) then 1676 1677 ! Allocate temporary arrays for the velocity-pressure resolution 1678 allocate(vect(3,ncelet)) 1679 1680 nswrsp = -1 1681 do iel = 1, ncel 1682 do isou = 1, 3 1683 smbr(isou,iel) = cell_f_vol(iel) 1684 enddo 1685 enddo 1686 do iel = 1, ncelet 1687 do isou = 1, 3 1688 vect(isou,iel) = 0.d0 1689 enddo 1690 enddo 1691 iescap = 0 1692 1693 ! We do not take into account transpose of grad 1694 ivisep = 0 1695 1696 vcopt_loc%nswrsm = nswrsp 1697 1698 p_k_value => vcopt_loc 1699 c_k_value = equation_param_from_vcopt(c_loc(p_k_value)) 1700 1701 call cs_equation_iterative_solve_vector & 1702 ( idtvar , iterns , & 1703 ivarfl(iu) , c_null_char , & 1704 ivisep , iescap , c_k_value , & 1705 vect , vect , & 1706 coefav , coefbv , cofafv , cofbfv , & 1707 imasfl , bmasfl , viscfi , & 1708 viscbi , viscf , viscb , secvif , & 1709 secvib , rvoid , rvoid , rvoid , & 1710 icvflb , ivoid , & 1711 fimpcp , smbr , vect , rvoid ) 1712 1713 do iel = 1, ncelet 1714 rom = crom(iel) 1715 do isou = 1, 3 1716 tpucou(isou,iel) = rom*vect(isou,iel) 1717 enddo 1718 do isou = 4, 6 1719 tpucou(isou,iel) = 0.d0 1720 enddo 1721 enddo 1722 1723 ! Free memory 1724 deallocate(vect) 1725 1726 endif 1727 1728 ! ---> The estimator on the predicted velocity is summed up over the components 1729 if (iestim(iespre).ge.0) then 1730 call field_get_val_s(iestim(iespre), c_estim) 1731 do iel = 1, ncel 1732 do isou = 1, 3 1733 c_estim(iel) = c_estim(iel) + eswork(isou,iel) 1734 enddo 1735 enddo 1736 endif 1737 1738 1739! ---> End of the construction of the total estimator: 1740! RHS resiudal of (U^{n+1}, P^{n+1}) + rho*volume*(U^{n+1} - U^n)/dt 1741else if (iappel.eq.2) then 1742 1743 inc = 1 1744 ! Pas de relaxation en stationnaire 1745 idtva0 = 0 1746 imasac = 0 1747 1748 vcopt_loc = vcopt_u 1749 1750 vcopt_loc%istat = -1 1751 vcopt_loc%idifft = -1 1752 vcopt_loc%iswdyn = -1 1753 vcopt_loc%nswrsm = -1 1754 vcopt_loc%iwgrec = 0 1755 vcopt_loc%blend_st = 0 ! Warning, may be overwritten if a field 1756 vcopt_loc%epsilo = -1 1757 vcopt_loc%epsrsm = -1 1758 1759 p_k_value => vcopt_loc 1760 c_k_value = equation_param_from_vcopt(c_loc(p_k_value)) 1761 1762 call cs_balance_vector & 1763 ( idtva0 , ivarfl(iu) , imasac , inc , ivisse , & 1764 c_k_value , vel , vel , coefav , coefbv , cofafv , & 1765 cofbfv , imasfl , bmasfl , viscf , viscb , secvif , & 1766 secvib , rvoid , rvoid , rvoid , & 1767 icvflb , icvfli , smbr ) 1768 1769 call field_get_val_s(iestim(iestot), c_estim) 1770 do iel = 1, ncel 1771 do isou = 1, 3 1772 c_estim(iel) = c_estim(iel) + (smbr(isou,iel)/volume(iel))**2 1773 enddo 1774 enddo 1775endif 1776 1777! ---> Finilaze estimators + Printings 1778 1779call field_get_id_try("predicted_velocity", f_id) 1780if (f_id.ge.0) then 1781 call field_get_val_v(f_id, cpro_pred_vel) 1782 do iel = 1, ncel 1783 do isou = 1, 3 1784 cpro_pred_vel(isou, iel) = vel(isou, iel) 1785 enddo 1786 enddo 1787endif 1788 1789if (iappel.eq.1) then 1790 1791 ! ---> Estimator on the predicted velocity: 1792 ! square root (norm) or square root of the sum times the volume (L2 norm) 1793 if (iestim(iespre).ge.0) then 1794 call field_get_val_s(iestim(iespre), c_estim) 1795 if (iescal(iespre).eq.1) then 1796 do iel = 1, ncel 1797 c_estim(iel) = sqrt(c_estim(iel)) 1798 enddo 1799 else if (iescal(iespre).eq.2) then 1800 do iel = 1, ncel 1801 c_estim(iel) = sqrt(c_estim(iel)*volume(iel)) 1802 enddo 1803 endif 1804 endif 1805 1806 ! ---> Norm printings 1807 if (vcopt_u%iwarni.ge.2) then 1808 rnorm = -1.d0 1809 do iel = 1, ncel 1810 vitnor = sqrt(vel(1,iel)**2+vel(2,iel)**2+vel(3,iel)**2) 1811 rnorm = max(rnorm,vitnor) 1812 enddo 1813 1814 if (irangp.ge.0) call parmax (rnorm) 1815 1816 write(nfecra,1100) rnorm 1817 1818 do iel = 1, ncel 1819 vitnor = sqrt(vel(1,iel)**2+vel(2,iel)**2+vel(3,iel)**2) 1820 rnorm = min(rnorm,vitnor) 1821 enddo 1822 1823 if (irangp.ge.0) call parmin (rnorm) 1824 1825 write(nfecra,1200) rnorm 1826 1827 endif 1828 1829! ---> Estimator on the whole Navier-Stokes: 1830! square root (norm) or square root of the sum times the volume (L2 norm) 1831else if (iappel.eq.2) then 1832 1833 call field_get_val_s(iestim(iestot), c_estim) 1834 if (iescal(iestot).eq.1) then 1835 do iel = 1, ncel 1836 c_estim(iel) = sqrt(c_estim(iel)) 1837 enddo 1838 else if (iescal(iestot).eq.2) then 1839 do iel = 1, ncel 1840 c_estim(iel) = sqrt(c_estim(iel)*volume(iel)) 1841 enddo 1842 endif 1843 1844endif 1845 1846! Free memory 1847!------------ 1848deallocate(smbr) 1849deallocate(fimp) 1850deallocate(fimpcp) 1851deallocate(tsexp) 1852deallocate(tsimp) 1853if (allocated(viscce)) deallocate(viscce) 1854if (allocated(divt)) deallocate(divt) 1855if (allocated(cproa_rho_tc)) deallocate(cproa_rho_tc) 1856 1857!-------- 1858! Formats 1859!-------- 1860 1861 1100 format(/, & 1862 1X,'Maximum velocity after prediction ',E12.4) 1863 1864 1200 format(/, & 1865 1X,'Minimum velocity after prediction ',E12.4) 1866 1867!---- 1868! End 1869!---- 1870 1871return 1872 1873end subroutine 1874