1!------------------------------------------------------------------------------- 2 3! This file is part of Code_Saturne, a general-purpose CFD tool. 4! 5! Copyright (C) 1998-2021 EDF S.A. 6! 7! This program is free software; you can redistribute it and/or modify it under 8! the terms of the GNU General Public License as published by the Free Software 9! Foundation; either version 2 of the License, or (at your option) any later 10! version. 11! 12! This program is distributed in the hope that it will be useful, but WITHOUT 13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15! details. 16! 17! You should have received a copy of the GNU General Public License along with 18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 19! Street, Fifth Floor, Boston, MA 02110-1301, USA. 20 21!------------------------------------------------------------------------------- 22 23!> \file calhyd.f90 24!> \brief Poisson equation resolution for hydrostatic pressure: 25!> \f$ \divs ( \grad P ) = \divs ( f ) \f$ 26!> 27!------------------------------------------------------------------------------ 28 29!------------------------------------------------------------------------------ 30! Arguments 31!------------------------------------------------------------------------------ 32! mode name role 33!------------------------------------------------------------------------------ 34!> \param[out] indhyd indicateur de mise a jour de phydr 35!> \param[in] iterns Navier-Stokes iteration number 36!> \param[in] fext external force generating hydrostatic pressure 37!> \param[in] dfext external force increment 38!> generating hydrostatic pressure 39!> \param[out] phydr hydrostatic pressure increment 40!> \param[in] flumas work array 41!> \param[in] flumab work array 42!> \param[in,out] viscf work array 43!> \param[in,out] viscb work array 44!> \param[in,out] dam work array 45!> \param[in,out] xam work array 46!> \param[in,out] dpvar work array 47!> \param[in,out] rhs work array 48!______________________________________________________________________________ 49 50subroutine calhyd & 51 ( indhyd , iterns , & 52 fext , dfext , & 53 phydr , flumas , flumab , & 54 viscf , viscb , & 55 dam , xam , & 56 dpvar , rhs ) & 57 bind(C, name='cs_hydrostatic_pressure_compute') 58 59!=============================================================================== 60! Module files 61!=============================================================================== 62 63use, intrinsic :: iso_c_binding 64use atincl, only: iatmst 65use paramx 66use numvar 67use entsor 68use cstnum 69use cstphy 70use optcal 71use parall 72use period 73use mesh 74use field_operator 75use cs_c_bindings 76 77!=============================================================================== 78 79implicit none 80 81! Arguments 82 83integer(c_int) :: indhyd 84integer(c_int), value :: iterns 85 86double precision fext(3,ncelet) 87double precision dfext(3,ncelet) 88double precision phydr(ncelet) 89double precision flumas(nfac), flumab(nfabor) 90double precision viscf(nfac), viscb(nfabor) 91double precision dam(ncelet), xam(nfac) 92double precision dpvar(ncelet) 93double precision rhs(ncelet) 94 95! Local variables 96 97character(len=80) :: chaine 98integer lchain 99integer f_id, iccocg, inc , init , isym 100integer iel , ical, ifac 101integer nswmpr 102integer isweep, niterf 103integer iphydp 104integer imrgrp, nswrgp, imligp, iwarnp, iwgrp 105integer idiffp, iconvp, ndircp, imvisp 106integer ibsize, iesize 107integer imucpp, f_id0 108 109double precision ressol, residu, rnorm , rnrmf , rnrmdf 110double precision epsrgp, climgp, extrap, epsilp 111double precision precre, precab, thetap 112double precision qimp, hint 113 114double precision rvoid(1) 115 116double precision, allocatable, dimension(:) :: rovsdt, div_fext, viscce 117double precision, allocatable, dimension(:,:) :: next_fext 118double precision, dimension(:), pointer :: coefap, coefbp, cofafp, cofbfp 119 120type(solving_info) sinfo 121type(var_cal_opt) :: vcopt_pr 122 123!=============================================================================== 124 125!=============================================================================== 126! 1. Initialisations 127!=============================================================================== 128 129! Allocate temporary arrays 130allocate(rovsdt(ncelet), div_fext(ncelet), viscce(ncelet)) 131allocate(next_fext(3, ncelet)) 132 133! Field id 134call field_get_id_try("hydrostatic_pressure", f_id) 135 136! Get variables calculation options 137call field_get_key_struct_var_cal_opt(f_id, vcopt_pr) 138call field_get_key_struct_solving_info(f_id, sinfo) 139 140if (iterns.le.1) then 141 sinfo%nbivar = 0 142endif 143 144call field_get_coefa_s (f_id, coefap) 145call field_get_coefb_s (f_id, coefbp) 146call field_get_coefaf_s(f_id, cofafp) 147call field_get_coefbf_s(f_id, cofbfp) 148 149chaine = 'hydrostatic_p' 150lchain = 16 151 152! Symmetric 153isym = 1 154 155! Matrix block size 156ibsize = 1 157iesize = 1 158 159! TEST DE VARIATION DE LA PRESSION HYDROSTATIQUE EN SORTIE 160 161! on regarde si les terme source ont varie 162! on ne passe dans calhyd que si on a des faces de sortie std 163! la precision pour les tests est a peu pres arbitraire. 164precre = sqrt(epzero) 165precab = 1.d2*epzero 166 167ical = 0 168do iel = 1, ncel 169 rnrmf = fext(1,iel)**2+fext(2,iel)**2+fext(3,iel)**2 170 rnrmdf = dfext(1,iel)**2+dfext(2,iel)**2+dfext(3,iel)**2 171 if ((rnrmdf.ge.precre*rnrmf).and.(rnrmdf.ge.precab)) then 172 ical = 1 173 endif 174enddo 175if (irangp.ge.0) then 176 call parcpt (ical) 177endif 178 179if (iatmst.eq.0) then 180 if (ical.eq.0) then 181 indhyd = 0 182 return 183 endif 184endif 185 186if (mod(ntcabs,ntlist).eq.0.or.vcopt_pr%iwarni.ge.1) write(nfecra,1000) 187 188 1000 format( & 189' Hydrostatic pressure computation: ',/, & 190' updating the Dirichlets at the end (CALHYD)',/) 191 192indhyd = 1 193 194f_id0 = -1 195 196do iel = 1, ncel 197 next_fext(1 ,iel) = fext(1 ,iel) * cell_is_active(iel) + dfext(1 ,iel) 198 next_fext(2 ,iel) = fext(2 ,iel) * cell_is_active(iel) + dfext(2 ,iel) 199 next_fext(3 ,iel) = fext(3 ,iel) * cell_is_active(iel) + dfext(3 ,iel) 200enddo 201 202! Parallel or periodicity synchronization 203call synvin(next_fext) 204 205!=============================================================================== 206! 2. Prepare matrix and boundary conditions 207!=============================================================================== 208 209! Boundary conditions 210 211! On resout avec des CL de flux nul partout 212ndircp = 0 213 214do ifac = 1, nfabor 215 ! LOCAL Neumann Boundary Conditions on the hydrostatic pressure 216 !-------------------------------------------------------------- 217 218 hint = 1.d0/distb(ifac) 219 qimp = 0.d0 220 221 call set_neumann_scalar & 222 !================== 223 ( coefap(ifac), cofafp(ifac), & 224 coefbp(ifac), cofbfp(ifac), & 225 qimp , hint ) 226 227enddo 228 229! Unsteady term 230do iel = 1, ncel 231 rovsdt(iel) = 0.d0 232enddo 233 234! Face diffusivity 235do iel = 1, ncel 236 viscce(iel) = 1.d0 237enddo 238 239imvisp = vcopt_pr%imvisf 240 241call viscfa & 242 ( imvisp , & 243 viscce , & 244 viscf , viscb ) 245 246iconvp = vcopt_pr%iconv 247idiffp = vcopt_pr%idiff 248thetap = vcopt_pr%thetav 249imucpp = 0 250 251call matrix & 252!========== 253 ( iconvp , idiffp , ndircp , isym , & 254 thetap , imucpp , & 255 coefbp , cofbfp , rovsdt , & 256 flumas , flumab , viscf , viscb , & 257 rvoid , dam , xam ) 258 259!=============================================================================== 260! 3. Compute right hand side 261!=============================================================================== 262 263init = 1 264inc = 0 265iccocg = 1 266nswrgp = vcopt_pr%nswrgr 267imligp = vcopt_pr%imligr 268iwarnp = vcopt_pr%iwarni 269epsrgp = vcopt_pr%epsrgr 270climgp = vcopt_pr%climgr 271 272! Compute div(f_ext^n+1) 273call projts & 274!========== 275 ( init , nswrgp , & 276 next_fext, & 277 cofbfp , & 278 flumas, flumab , & 279 viscf , viscb , & 280 viscce , viscce , viscce ) 281 282init = 1 283call divmas(init,flumas,flumab, div_fext) 284rnorm = sqrt(cs_gdot(ncel,div_fext,div_fext)) 285 286!=============================================================================== 287! 4. BOUCLES SUR LES NON ORTHOGONALITES (RESOLUTION) 288!=============================================================================== 289 290! --- Nombre de sweeps 291nswmpr = vcopt_pr%nswrsm 292 293! --- Mise a zero des variables 294! dpvar sera l'increment d'increment a chaque sweep 295! div_fext sera la divergence du flux de masse predit 296do iel = 1, ncel 297 dpvar(iel) = 0.d0 298enddo 299 300! Reconstruction loop (beginning) 301!-------------------------------- 302 303do isweep = 1, nswmpr 304 305 ! --- Update the right hand side and update the residual 306 ! rhs^{k+1} = - div(fext^n+1) - D(1, p_h^{k+1}) 307 !------------------------------------------------------------- 308 309 iccocg = 1 310 init = 1 ! re-init rhs to 0 if init = 1 311 inc = 1 312 imrgrp = vcopt_pr%imrgra 313 nswrgp = vcopt_pr%nswrgr 314 imligp = vcopt_pr%imligr 315 iwgrp = vcopt_pr%iwgrec 316 iwarnp = vcopt_pr%iwarni 317 epsrgp = vcopt_pr%epsrgr 318 climgp = vcopt_pr%climgr 319 extrap = 0.d0 320 iphydp = 1 321 322 call itrgrp & 323 !========== 324 ( f_id0 , init , inc , imrgrp , iccocg , nswrgp , imligp , iphydp , & 325 iwgrp , iwarnp , & 326 epsrgp , climgp , extrap , & 327 next_fext, & 328 phydr , & 329 coefap , coefbp , & 330 cofafp , cofbfp , & 331 viscf , viscb , & 332 viscce , & 333 rhs ) 334 335 336 do iel = 1, ncel 337 rhs(iel) = - div_fext(iel) - rhs(iel) 338 enddo 339 340 ! --- Convergence test 341 residu = sqrt(cs_gdot(ncel,rhs,rhs)) 342 if (vcopt_pr%iwarni.ge.2) then 343 write(nfecra,1400) chaine(1:16), isweep,residu, rnorm 344 endif 345 346 if (isweep.eq.1) then 347 sinfo%rnsmbr = residu 348 endif 349 350 if (residu.le.vcopt_pr%epsrsm*rnorm) then 351 ! Si convergence, sortie 352 353 goto 101 354 355 endif 356 357 ! Solving on the increment 358 !------------------------- 359 do iel = 1, ncel 360 dpvar(iel) = 0.d0 361 enddo 362 363 iwarnp = vcopt_pr%iwarni 364 epsilp = vcopt_pr%epsilo 365 366 ! Solver residual 367 ressol = residu 368 369 call sles_solve_native(f_id, '', & 370 isym, ibsize, iesize, dam, xam, & 371 epsilp, rnorm, niterf, ressol, rhs, dpvar) 372 373 ! Writing 374 sinfo%nbivar = sinfo%nbivar + niterf 375 376 ! Update the increment of pressure 377 !--------------------------------- 378 379 do iel = 1, ncel 380 phydr(iel) = phydr(iel) + dpvar(iel) 381 enddo 382 383enddo 384 385! --- Reconstruction loop (end) 386 387if (vcopt_pr%iwarni.ge.2) then 388 write(nfecra,1600) chaine(1:16), nswmpr 389endif 390 391 101 continue 392 393 ! For logging 394if (abs(rnorm).gt.0.d0) then 395 sinfo%resvar = residu/rnorm 396else 397 sinfo%resvar = 0.d0 398endif 399 400! Save convergence info 401call field_set_key_struct_solving_info(f_id, sinfo) 402 403! Free memory 404deallocate(rovsdt, div_fext, viscce) 405deallocate(next_fext) 406 407!=============================================================================== 408! 5. Free solver setup 409!=============================================================================== 410 411call sles_free_native(f_id, '') 412 413!-------- 414! Formats 415!-------- 416 417 1400 format(1X,A16,' : sweep = ',I5,' RHS residual = ',E14.6, ' NORM =',E14.6) 418 1600 format( & 419'@ ',/,& 420'@ @@ WARNING: ', A16,' HYDROSTATIC PRESSURE STEP ',/,& 421'@ ======== ',/,& 422'@ Maximum number of iterations ',I10 ,' reached ',/,& 423'@ ' ) 424 425!---- 426! End 427!---- 428 429return 430 431end subroutine 432