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