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 resrit.f90
28!>
29!> \brief This subroutine perform the solving of the transport equation
30!> of the turbulent heat fluxes.
31!-------------------------------------------------------------------------------
32
33!-------------------------------------------------------------------------------
34! Arguments
35!______________________________________________________________________________.
36!  mode           name          role
37!______________________________________________________________________________!
38!> \param[in]     nscal         total number of scalars
39!> \param[in]     iscal         number of the scalar used
40!> \param[in]     xcpp          \f$ C_p \f$
41!> \param[in,out] xut, xuta     calculated variables at cell centers
42!>                               (at current and previous time steps)
43!> \param[in]     dt            time step (per cell)
44!> \param[in]     gradv         mean velocity gradient
45!> \param[in]     gradt         mean scalar gradient
46!> \param[in]     grad_al       alpha scalar gradient
47!______________________________________________________________________________!
48
49subroutine resrit &
50 ( nscal  ,                                                       &
51   iscal  , xcpp   , xut    , xuta   ,                            &
52   dt     ,                                                       &
53   gradv  , gradt  , grad_al)
54
55!===============================================================================
56! Module files
57!===============================================================================
58
59use, intrinsic :: iso_c_binding
60
61use paramx
62use numvar
63use entsor
64use optcal
65use cstnum
66use cstphy
67use parall
68use period
69use field
70use mesh
71use cs_f_interfaces
72use cs_c_bindings
73
74!===============================================================================
75
76implicit none
77
78! Arguments
79
80integer          nscal , iscal
81
82double precision dt(ncelet)
83double precision xcpp(ncelet), xut(3,ncelet), xuta(3,ncelet)
84double precision gradv(3,3,ncelet)
85double precision gradt(3,ncelet)
86double precision grad_al(3,ncelet)
87
88! Local variables
89
90integer          iel
91integer          ii, ivar
92integer          iflmas, iflmab
93integer          iwarnp
94integer          imvisp
95integer          iescap
96integer          st_prv_id
97integer          ivisep, ifcvsl
98integer          isou, jsou
99integer          itt
100integer          icvflb
101integer          f_id
102integer          keyvar, iut
103integer          init
104integer          kturt, turb_flux_model
105
106integer          ivoid(1)
107
108double precision trrij
109double precision thets , thetv , thetp1
110double precision xttke , prdtl
111double precision grav(3)
112double precision xrij(3,3),phiith(3), phiitw(3)
113double precision xnal(3), xnoral
114double precision alpha, xttdrbt, xttdrbw
115double precision pk, gk, xxc1, xxc2, xxc3, imp_term
116double precision rctse, visls_0
117
118double precision rvoid(1)
119
120character(len=80) :: fname, name
121
122double precision, allocatable, dimension(:,:) :: viscce
123double precision, allocatable, dimension(:) :: viscb
124double precision, allocatable, dimension(:) :: viscf
125double precision, allocatable, dimension(:,:) :: weighf
126double precision, allocatable, dimension(:) :: weighb
127double precision, allocatable, dimension(:,:) :: smbrut
128double precision, allocatable, dimension(:,:,:) :: fimp
129double precision, allocatable, dimension(:) :: w1
130
131double precision, dimension(:,:), pointer :: coefav, cofafv, visten
132double precision, dimension(:,:,:), pointer :: coefbv, cofbfv
133double precision, dimension(:), pointer :: imasfl, bmasfl
134double precision, dimension(:), pointer :: crom, cpro_beta
135double precision, dimension(:), pointer :: cvar_al
136double precision, dimension(:), pointer :: cvar_ep
137double precision, dimension(:,:), pointer :: cvar_rij
138double precision, dimension(:), pointer :: cvar_tt, cvara_tt
139double precision, dimension(:), pointer :: viscl, visct, viscls, c_st_prv
140
141type(var_cal_opt) :: vcopt, vcopt_ut
142type(var_cal_opt), target :: vcopt_loc
143type(var_cal_opt), pointer :: p_k_value
144type(c_ptr) :: c_k_value
145
146!===============================================================================
147
148!===============================================================================
149! 1. Initialization
150!===============================================================================
151
152! Allocate work arrays
153allocate(w1(ncelet))
154allocate(viscce(6,ncelet))
155allocate(smbrut(3,ncelet))
156allocate(fimp(3,3,ncelet))
157allocate(viscf(nfac), viscb(nfabor))
158allocate(weighf(2,nfac))
159allocate(weighb(nfabor))
160
161if ((itytur.eq.2).or.(itytur.eq.5).or.(iturb.eq.60)) then
162  write(nfecra,*)'Utiliser un modele Rij avec ces modeles de thermiques'!FIXME
163  call csexit(1)
164endif
165
166call field_get_val_s(icrom, crom)
167call field_get_val_s(iviscl, viscl)
168call field_get_val_s(ivisct, visct)
169
170call field_get_val_s(ivarfl(iep), cvar_ep)
171call field_get_val_v(ivarfl(irij), cvar_rij)
172
173call field_get_key_int(ivarfl(iu), kimasf, iflmas)
174call field_get_key_int(ivarfl(iu), kbmasf, iflmab)
175call field_get_val_s(iflmas, imasfl)
176call field_get_val_s(iflmab, bmasfl)
177
178if (ibeta.ge.0) then
179  call field_get_val_s(ibeta, cpro_beta)
180endif
181
182call field_get_val_v(ivsten, visten)
183
184ivar = isca(iscal)
185
186call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt)
187
188call field_get_key_id("variable_id", keyvar)
189! Get name of the reference scalar
190call field_get_name(ivarfl(ivar), fname)
191! Index of the corresponding turbulent flux
192call field_get_id(trim(fname)//'_turbulent_flux', f_id)
193! Set pointer values of turbulent fluxes
194call field_get_key_int(f_id, keyvar, iut)
195
196call field_get_key_struct_var_cal_opt(ivarfl(iut), vcopt_ut)
197
198if (vcopt%iwarni.ge.1) then
199  call field_get_name(ivarfl(ivar), name)
200  write(nfecra,1000) trim(name)//'_turbulent_flux'
201endif
202
203! Get the turbulent flux model
204call field_get_key_id('turbulent_flux_model', kturt)
205call field_get_key_int(ivarfl(isca(iscal)), kturt, turb_flux_model)
206
207! S pour Source, V pour Variable
208thets  = thetst
209thetv  = vcopt%thetav
210
211call field_get_key_int(ivarfl(ivar), kstprv, st_prv_id)
212if (st_prv_id.ge.0) then
213  call field_get_val_s(st_prv_id, c_st_prv)
214else
215  c_st_prv=> null()
216endif
217
218call field_get_key_int (ivarfl(ivar), kivisl, ifcvsl)
219if (ifcvsl .ge. 0) then
220  call field_get_val_s(ifcvsl, viscls)
221else
222  call field_get_key_double(ivarfl(ivar), kvisl0, visls_0)
223endif
224
225do iel = 1, ncelet
226  do isou = 1, 3
227    smbrut(isou,iel) = 0.d0
228    do jsou = 1, 3
229      fimp(isou,jsou,iel) = 0.d0
230    enddo
231  enddo
232enddo
233
234! Find the corresponding variance of the scalar iscal
235itt = -1
236if ((abs(gx)+abs(gy)+abs(gz)).gt.0) then
237  grav(1) = gx
238  grav(2) = gy
239  grav(3) = gz
240  do ii = 1, nscal
241    if (iscavr(ii).eq.iscal) itt = ii
242  enddo
243endif
244
245!===============================================================================
246! 2. Mass source terms FIXME
247!===============================================================================
248
249if (st_prv_id.ge.0) then
250  do iel = 1, ncel
251    do isou = 1,3
252      smbrut(isou,iel) = fimp(isou,isou,iel)*xuta(isou,iel)
253      fimp(isou,isou,iel) = - thetv*fimp(isou,isou,iel)
254    enddo
255  enddo
256! If we do not extrapolate the source terms
257else
258  do iel = 1, ncel
259    do isou = 1, 3
260      ! User source term
261      smbrut(isou,iel) = smbrut(isou,iel) + fimp(isou,isou,iel)*xuta(isou,iel)
262      ! Diagonal
263      fimp(isou,isou,iel) = max(-fimp(isou,isou,iel),zero)
264    enddo
265  enddo
266endif
267
268!===============================================================================
269! 3. Unsteady term
270!===============================================================================
271
272do iel = 1, ncel
273  do isou = 1, 3
274    fimp(isou,isou,iel) = fimp(isou,isou,iel)                                  &
275                        + vcopt%istat*(crom(iel)/dt(iel))*cell_f_vol(iel)
276  enddo
277enddo
278
279!===============================================================================
280! 4. Right Hand Side of the thermal fluxes:
281!     rho*(Pit + Git + Phi*_it - eps_it)
282!===============================================================================
283
284if (itt.gt.0) then
285  call field_get_val_s(ivarfl(isca(itt)), cvar_tt)
286  call field_get_val_prev_s(ivarfl(isca(itt)), cvara_tt)
287endif
288if (turb_flux_model.eq.31) then
289  ! Name of the scalar
290  call field_get_name(ivarfl(isca(iscal)), fname)
291
292  ! Index of the corresponding turbulent flux
293  call field_get_id(trim(fname)//'_alpha', f_id)
294
295  call field_get_val_s(f_id, cvar_al)
296endif
297
298do iel = 1, ncel
299
300  xrij(1,1) = cvar_rij(1,iel)
301  xrij(2,2) = cvar_rij(2,iel)
302  xrij(3,3) = cvar_rij(3,iel)
303  xrij(1,2) = cvar_rij(4,iel)
304  xrij(2,3) = cvar_rij(5,iel)
305  xrij(1,3) = cvar_rij(6,iel)
306
307  xrij(2,1) = xrij(1,2)
308  xrij(3,1) = xrij(1,3)
309  xrij(3,2) = xrij(2,3)
310
311  if (ifcvsl.ge.0) then
312    prdtl = viscl(iel)*xcpp(iel)/viscls(iel)
313  else
314    prdtl = viscl(iel)*xcpp(iel)/visls_0
315  endif
316
317  trrij  = 0.5d0*(xrij(1,1) + xrij(2,2) + xrij(3,3))
318  ! --- Compute Durbin time scheme
319  xttke  = trrij/cvar_ep(iel)
320
321  if (turb_flux_model.eq.31) then
322    alpha = cvar_al(iel)
323    !FIXME Warning / rhebdfm**0.5 compared to F Dehoux
324    xttdrbt = xttke * sqrt( (1.d0-alpha)*prdtl/ rhebdfm + alpha )
325    xttdrbw = xttdrbt * sqrt(rhebdfm/prdtl) ! And so multiplied by (R/Prandt)^0.5
326
327    ! Compute the unite normal vector
328    xnoral = &
329      ( grad_al(1, iel)*grad_al(1, iel) &
330      + grad_al(2, iel)*grad_al(2, iel) &
331      + grad_al(3, iel)*grad_al(3, iel))
332    xnoral = sqrt(xnoral)
333
334    if (xnoral.le.epzero/cell_f_vol(iel)**(1.d0/3.d0)) then
335      xnal(1) = 0.d0
336      xnal(2) = 0.d0
337      xnal(3) = 0.d0
338    else
339      xnal(1) = grad_al(1, iel)/xnoral
340      xnal(2) = grad_al(2, iel)/xnoral
341      xnal(3) = grad_al(3, iel)/xnoral
342    endif
343
344    ! Production and buoyancy for TKE
345    pk = -( xrij(1,1)*gradv(1,1,iel) +&
346            xrij(1,2)*gradv(1,2,iel) +&
347            xrij(1,3)*gradv(1,3,iel) +&
348            xrij(2,1)*gradv(2,1,iel) +&
349            xrij(2,2)*gradv(2,2,iel) +&
350            xrij(2,3)*gradv(2,3,iel) +&
351            xrij(3,1)*gradv(3,1,iel) +&
352            xrij(3,2)*gradv(3,2,iel) +&
353            xrij(3,3)*gradv(3,3,iel) )
354
355    if (ibeta.ge.0) then
356      gk = cpro_beta(iel)*(xuta(1,iel)*gx + xuta(2,iel)*gy + xuta(3,iel)*gz)!FIXME make buoyant term coherent elsewhere
357    else
358      gk = 0.d0
359    endif
360    xxc1 = 1.d0+2.d0*(1.d0 - cvar_al(iel))*(pk+gk)/cvar_ep(iel)
361    xxc2 = 0.5d0*(1.d0+1.d0/prdtl)*(1.d0-0.3d0*(1.d0 - cvar_al(iel)) &
362      *(pk+gk)/cvar_ep(iel))
363    xxc3 = xxc2
364
365  else
366    xttdrbt = xttke
367    xttdrbw = xttke
368    alpha = 1.d0
369    xxc1 = 0.d0
370    xxc2 = 0.d0
371    xxc3 = 0.d0
372    xnal(1) = 0.d0
373    xnal(2) = 0.d0
374    xnal(3) = 0.d0
375  endif
376
377  do isou = 1, 3
378    phiith(isou) = -c1trit / xttdrbt * xuta(isou,iel)               &
379                 + c2trit*(xuta(1,iel)*gradv(1,isou,iel)            &
380                          +xuta(2,iel)*gradv(2,isou,iel)            &
381                          +xuta(3,iel)*gradv(3,isou,iel))           &
382                 + c4trit*(-xrij(isou,1)*gradt(1,iel)               &
383                           -xrij(isou,2)*gradt(2,iel)               &
384                           -xrij(isou,3)*gradt(3,iel))
385    if (itt.gt.0.and.ibeta.ge.0) then
386      phiith(isou) = phiith(isou)                                              &
387             + c3trit*(cpro_beta(iel)*grav(isou)*cvar_tt(iel))
388    endif
389
390    phiitw(isou) = -1.d0/xttdrbw * xxc1 *             & !FIXME full implicit
391                 ( xuta(1,iel)*xnal(1)*xnal(isou)     &
392                 + xuta(2,iel)*xnal(2)*xnal(isou)     &
393                 + xuta(3,iel)*xnal(3)*xnal(isou))
394
395    ! Pressure/thermal fluctuation correlation term
396    !----------------------------------------------
397    smbrut(isou,iel) = smbrut(isou,iel) +                                      &
398                cell_f_vol(iel)*crom(iel)*(  alpha *         phiith(isou)          &
399                                      + (1.d0 - alpha) * phiitw(isou))
400
401    imp_term = max(cell_f_vol(iel)*crom(iel)*(                                     &
402              alpha        * (c1trit/xttdrbt-c2trit*gradv(isou,isou,iel))      &
403            + (1.d0-alpha) * (xxc1*xnal(isou)*xnal(isou) / xttdrbw) &
404            ), 0.d0)
405
406    fimp(isou,isou,iel) = fimp(isou,isou,iel) + imp_term
407
408    ! Production terms
409    !-----------------
410    smbrut(isou,iel) = smbrut(isou,iel)                              &
411                     + cell_f_vol(iel)*crom(iel)                         &
412                       ! Production term due to the mean velcoity
413                       *( -gradv(1,isou,iel)*xuta(1,iel)             &
414                          -gradv(2,isou,iel)*xuta(2,iel)             &
415                          -gradv(3,isou,iel)*xuta(3,iel)             &
416                       ! Production term due to the mean temperature
417                         -xrij(1,isou)*gradt(1,iel)                  &
418                         -xrij(2,isou)*gradt(2,iel)                  &
419                         -xrij(3,isou)*gradt(3,iel)                  &
420                        )
421
422    ! Production term due to the gravity
423    if (itt.gt.0.and.ibeta.ge.0) then
424      smbrut(isou,iel) = smbrut(isou,iel)                            &
425                       + cell_f_vol(iel)*crom(iel)*(            &
426               -grav(isou)*cpro_beta(iel)*cvara_tt(iel))
427    endif
428
429    ! Dissipation (Wall term only because "h" term is zero
430    smbrut(isou,iel) = smbrut(isou,iel) -                                      &
431                cell_f_vol(iel)*crom(iel)*(1.d0 - alpha) / xttdrbw *               &
432                       ( xxc2 * xuta(isou, iel)                                &
433                       + xxc3*( xuta(1,iel)*xnal(1)*xnal(isou)                 &
434                              + xuta(2,iel)*xnal(2)*xnal(isou)                 &
435                              + xuta(3,iel)*xnal(3)*xnal(isou)))
436
437    ! TODO we can implicite more terms
438    imp_term = max(cell_f_vol(iel)*crom(iel)*(1.d0 - alpha) / xttdrbw *            &
439                  ( xxc2 + xxc3 * xnal(isou)*xnal(isou)), 0.d0)
440    fimp(isou,isou,iel) = fimp(isou,isou,iel) + imp_term
441
442  enddo
443enddo
444
445!===============================================================================
446! 5. Tensor diffusion
447!===============================================================================
448! Symmetric tensor diffusivity (GGDH)
449if (iand(vcopt_ut%idften, ANISOTROPIC_RIGHT_DIFFUSION).ne.0) then
450
451  call field_get_key_double(ivarfl(isca(iscal)), kvisl0, visls_0)
452
453  do iel = 1, ncel
454
455    if (ifcvsl.ge.0) then
456      prdtl = viscl(iel)*xcpp(iel)/viscls(iel)
457    else
458      prdtl = viscl(iel)*xcpp(iel)/visls_0
459    endif
460
461    do isou = 1, 6
462      if (isou.le.3) then
463        viscce(isou,iel) = 0.5d0*(viscl(iel)*(1.d0+1.d0/prdtl))    &
464                         + ctheta(iscal)*visten(isou,iel)/csrij
465      else
466        viscce(isou,iel) = ctheta(iscal)*visten(isou,iel)/csrij
467      endif
468    enddo
469  enddo
470
471  iwarnp = vcopt%iwarni
472
473  call vitens &
474  ( viscce , iwarnp ,             &
475    weighf , weighb ,             &
476    viscf  , viscb  )
477
478! Scalar diffusivity
479else
480
481  do iel = 1, ncel
482    rctse = ctheta(iscal) * visct(iel) / cmu
483    w1(iel) = viscl(iel) + vcopt%idifft*rctse
484  enddo
485
486  imvisp = vcopt%imvisf
487
488  call viscfa                    &
489  ( imvisp ,                     &
490   w1     ,                      &
491   viscf  , viscb  )
492
493endif
494
495!===============================================================================
496! 6. Vectorial solving of the turbulent thermal fluxes
497!===============================================================================
498
499if (st_prv_id.ge.0) then
500  thetp1 = 1.d0 + thets
501  do iel = 1, ncel
502    do isou = 1, 3
503      smbrut(isou,iel) = smbrut(isou,iel) + thetp1*c_st_prv(iel) !FIXME
504    enddo
505  enddo
506endif
507
508! Name of the scalar ivar
509call field_get_name(ivarfl(ivar), fname)
510
511! Index of the corresponding turbulent flux
512call field_get_id(trim(fname)//'_turbulent_flux', f_id)
513
514call field_get_coefa_v(f_id,coefav)
515call field_get_coefb_v(f_id,coefbv)
516call field_get_coefaf_v(f_id,cofafv)
517call field_get_coefbf_v(f_id,cofbfv)
518
519iescap = 0
520
521! We do not take into account transpose of grad
522ivisep = 0
523
524! all boundary convective flux with upwind
525icvflb = 0
526init   = 1
527
528vcopt_loc = vcopt
529
530vcopt_loc%istat  = -1
531vcopt_loc%idifft = -1
532vcopt_loc%iwgrec = 0
533vcopt_loc%thetav = thetv
534vcopt_loc%blend_st = 0 ! Warning, may be overwritten if a field
535
536p_k_value => vcopt_loc
537c_k_value = equation_param_from_vcopt(c_loc(p_k_value))
538
539call cs_equation_iterative_solve_vector                     &
540 ( idtvar , init   ,                                        &
541   f_id   , c_null_char ,                                   &
542   ivisep , iescap , c_k_value       ,                      &
543   xuta   , xuta   ,                                        &
544   coefav , coefbv , cofafv , cofbfv ,                      &
545   imasfl , bmasfl , viscf  ,                               &
546   viscb  , viscf  , viscb  , rvoid  ,                      &
547   rvoid  , viscce , weighf , weighb ,                      &
548   icvflb , ivoid  ,                                        &
549   fimp   , smbrut , xut    , rvoid )
550
551!===============================================================================
552! 7. Writings
553!===============================================================================
554
555! Free memory
556deallocate(viscce)
557deallocate(viscf, viscb)
558deallocate(smbrut)
559deallocate(fimp)
560deallocate(weighf, weighb)
561deallocate(w1)
562!--------
563! Formats
564!--------
565
566 1000 format(/,'           Solving variable ',A23           ,/)
567
568!----
569! End
570!----
571
572return
573end subroutine
574