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!===============================================================================
25! Function:
26! --------
27!> \file visdyn.f90
28!> \brief Calculation of turbulent viscosity for
29!>        a dynamic Smagorinsky LES model
30!>
31!> \f[ smago = \dfrac{L_{ij}M_{ij}}{M_{ij}M_{ij}} \f]
32!>
33!> \f[ \mu_T = \rho smago L^2  \sqrt{2 S_{ij}S_{ij}} \f]
34!> \f[ S_{ij} = \dfrac{\der{u_i}{x_j} + \der{u_j}{x_i}}{2}\f]
35!>
36!> We have at edge faces types at previous time step
37!>   (except at first time step, when tables itypfb and itrifb
38!>   have not been filled).
39!>
40!> Please refer to the
41!> <a href="../../theory.pdf#dynsmago"><b>dynamic Smagorinsky model</b></a>
42!> section of the theory guide for more informations.
43!-------------------------------------------------------------------------------
44
45!-------------------------------------------------------------------------------
46! Arguments
47!______________________________________________________________________________.
48!  mode           name          role
49!______________________________________________________________________________!
50!> \param[in]     nvar          total number of variables
51!> \param[in]     nscal         total number of scalars
52!> \param[in]     ncepdp        number of cells with head loss
53!> \param[in]     ncesmp        number of cells with mass source term
54!> \param[in]     icepdc        number of ncepdp cells with losses
55!> \param[in]     icetsm        number of cells with mass source
56!> \param[in]     itypsm        type of mass source for the variable
57!>                               (cf. cs_user_mass_source_terms)
58!> \param[in]     dt            time step (per cell)
59!> \param[in]     ckupdc        work array for head losses
60!> \param[in]     smacel        value of variables associated to the
61!>                               mass source
62!>                               for ivar = ipr, smacel = mass flux
63!> \param[out]    gradv         the computed velocity gradients
64!______________________________________________________________________________!
65
66subroutine visdyn &
67 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
68   icepdc , icetsm , itypsm ,                                     &
69   dt     ,                                                       &
70   ckupdc , smacel, gradv )
71
72!===============================================================================
73! Module files
74!===============================================================================
75
76use paramx
77use numvar
78use cstnum
79use optcal
80use cstphy
81use ppincl
82use entsor
83use parall
84use period
85use mesh
86use field
87use field_operator
88use cs_c_bindings
89
90!===============================================================================
91
92implicit none
93
94! Arguments
95
96integer          nvar   , nscal
97integer          ncepdp , ncesmp
98
99integer          icepdc(ncepdp)
100integer          icetsm(ncesmp), itypsm(ncesmp,nvar)
101
102double precision dt(ncelet)
103double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar)
104double precision, intent(inout) :: gradv(3,3,ncelet)
105
106! Local variables
107
108integer          ii, jj, iel, inc, ivar
109integer          iprev
110integer          iclipc
111integer          iccocg, iclip
112integer          key_turb_diff, key_sgs_sca_coef
113integer          t_dif_id, sca_dync_id
114integer          nswrgp, imligp, iwarnp
115
116double precision coef, radeux, deux, delta, deltaf
117double precision s11, s22, s33, s11f, s22f, s33f
118double precision dudy, dudz, dvdx, dvdz, dwdx, dwdy
119double precision dudyf, dudzf, dvdxf, dvdzf, dwdxf, dwdyf
120double precision xfil, xa, xb, xfil2, xsmgmx, xsmgmn
121double precision xl11, xl22, xl33, xl12, xl13, xl23
122double precision xm11, xm22, xm33, xm12, xm13, xm23
123double precision smagma, smagmi, smagmy
124double precision scal1, scal2, scal3
125double precision epsrgp, climgp
126
127double precision, allocatable, dimension(:) :: w1, w2, w3
128double precision, allocatable, dimension(:) :: w4, w5, w6
129double precision, allocatable, dimension(:) :: w7, w8, w9
130double precision, allocatable, dimension(:) :: s_n, sf_n
131double precision, allocatable, dimension(:) :: w0, xrof, xro
132double precision, allocatable, dimension(:,:) :: grads, scami, scamif
133double precision, allocatable, dimension(:,:) :: xmij, w61, w62, gradsf
134double precision, dimension(:,:), pointer :: coefau
135double precision, dimension(:,:,:), pointer :: coefbu
136double precision, dimension(:), pointer :: crom, coefas, coefbs, cvar_sca
137double precision, dimension(:), pointer :: cpro_turb_diff
138double precision, dimension(:,:), pointer :: vel
139double precision, dimension(:), pointer :: visct, cpro_smago, cpro_sca_dync
140
141type(var_cal_opt) :: vcopt
142
143!===============================================================================
144
145!===============================================================================
146! 1.  Iniatilization
147!===============================================================================
148
149! Map field arrays
150call field_get_val_v(ivarfl(iu), vel)
151
152call field_get_coefa_v(ivarfl(iu), coefau)
153call field_get_coefb_v(ivarfl(iu), coefbu)
154
155call field_get_val_s(ivisct, visct)
156call field_get_val_s(icrom, crom)
157
158call field_get_val_s(ismago, cpro_smago)
159
160! For the calculation of the viscosity of the sub-mesh
161xfil   = xlesfl
162xfil2  = xlesfd
163xa     = ales
164xb     = bles
165deux   = 2.d0
166radeux = sqrt(deux)
167xsmgmx = smagmx
168xsmgmn = smagmn
169
170! Allocate some work arrays
171
172allocate(w0(ncelet), w1(ncelet))
173allocate(xmij(6,ncelet))
174allocate(xro(ncelet), xrof(ncelet))
175
176! Take into account variable density case: Favre filtering
177! Constant density case: Reynolds filtering
178if(irovar.eq.1) then
179  do iel = 1, ncel
180    xro(iel) = crom(iel)
181  enddo
182else
183  xro(:) = 1.d0
184endif
185
186! In case of constant density, xrof always 1.d0
187call les_filter(1, xro, xrof)
188!===============================================================================
189! 2.  Calculation of velocity gradient and of
190!       S11**2+S22**2+S33**2+2*(S12**2+S13**2+S23**2)
191!===============================================================================
192
193! Allocate temporary arrays for gradients calculation
194allocate(s_n(ncelet), sf_n(ncelet))
195allocate(w61(6,ncelet), w62(6,ncelet))
196
197inc = 1
198iprev = 0
199
200call field_gradient_vector(ivarfl(iu), iprev, 0, inc, gradv)
201
202do iel = 1, ncel
203
204  ! gradv(iel, xyz, uvw)
205  s11   = gradv(1, 1, iel)
206  s22   = gradv(2, 2, iel)
207  s33   = gradv(3, 3, iel)
208  dudy  = gradv(2, 1, iel)
209  dudz  = gradv(3, 1, iel)
210  dvdx  = gradv(1, 2, iel)
211  dvdz  = gradv(3, 2, iel)
212  dwdx  = gradv(1, 3, iel)
213  dwdy  = gradv(2, 3, iel)
214
215
216
217! In the case of constant density, s11+s22+s33 is zero
218  xmij(1,iel) = s11 - irovar*1.0d0/3.0d0*(s11+s22+s33)
219  xmij(2,iel) = s22 - irovar*1.0d0/3.0d0*(s11+s22+s33)
220  xmij(3,iel) = s33 - irovar*1.0d0/3.0d0*(s11+s22+s33)
221  xmij(4,iel) = 0.5d0*(dudy+dvdx)
222  xmij(5,iel) = 0.5d0*(dudz+dwdx)
223  xmij(6,iel) = 0.5d0*(dvdz+dwdy)
224
225  s_n(iel) = radeux*sqrt(                                      &
226             xmij(1,iel)**2 + xmij(2,iel)**2 + xmij(3,iel)**2  &
227           + 2.d0*(xmij(4,iel)**2 + xmij(5,iel)**2             &
228           + xmij(6,iel)**2) )
229  w62(:,iel) = xro(iel)*xmij(:,iel)
230enddo
231
232! w62 temperarily contains rho*S
233call les_filter(6, w62, w61)
234
235! w61 <rho*S>/<rho>, sf_n is ||<rho*S>/<rho>||
236do iel = 1, ncel
237  w61(:,iel) = w61(:,iel)/xrof(iel)
238  sf_n(iel) = radeux*sqrt(                                      &
239              w61(1,iel)**2 + w61(2,iel)**2 + w61(3,iel)**2     &
240            + 2.d0*(w61(4,iel)**2 + w61(5,iel)**2               &
241            + w61(6,iel)**2) )
242enddo
243
244!     Here XMIJ contains Sij
245!         S_n contains ||S||
246!            sqrt(2)*sqrt(S11^2+S22^2+S33^2+2(S12^2+S13^2+S23^2))
247!         Sf_n               contains ||SF||
248!            sqrt(2)*sqrt(S11F^2+S22F^2+S33F^2+2(S12F^2+S13F^2+S23F^2))
249
250!===============================================================================
251! 3.  Calculation of Mij
252!===============================================================================
253
254do iel = 1, ncel
255  w0(iel) = xfil *(xa*volume(iel))**xb
256enddo
257
258! Reuse xmij as temporary array
259
260do iel = 1, ncel
261  delta = w0(iel)
262  do ii = 1, 6
263    xmij(ii,iel) = -deux*xro(iel)*delta**2*s_n(iel)*xmij(ii,iel)
264  enddo
265enddo
266
267! w62 now contains <-2*rho*delta**2*||S||*S>
268call les_filter(6, xmij, w62)
269
270! Now compute final xmij value: M_ij = alpha_ij - beta_ij
271
272do iel = 1, ncel
273  delta = w0(iel)
274  deltaf = xfil2*delta
275  do ii = 1, 6
276    xmij(ii,iel) = -deux*xrof(iel)*deltaf**2*sf_n(iel)*w61(ii,iel) - w62(ii,iel)
277  enddo
278enddo
279
280deallocate(w61, w62)
281
282!===============================================================================
283! 4.  Calculation of the dynamic Smagorinsky constant
284!===============================================================================
285
286! Allocate work arrays
287allocate(w2(ncelet), w3(ncelet), w4(ncelet))
288allocate(w5(ncelet), w6(ncelet), w7(ncelet))
289allocate(w8(ncelet), w9(ncelet))
290
291! Filtering the velocity and its square
292
293! U**2
294do iel = 1,ncel
295  w0(iel) = xro(iel)*vel(1,iel)*vel(1,iel)
296enddo
297call les_filter(1, w0, w1)
298
299! V**2
300do iel = 1,ncel
301  w0(iel) = xro(iel)*vel(2,iel)*vel(2,iel)
302enddo
303call les_filter(1, w0, w2)
304
305! W**2
306do iel = 1,ncel
307  w0(iel) = xro(iel)*vel(3,iel)*vel(3,iel)
308enddo
309call les_filter(1, w0, w3)
310
311! UV
312do iel = 1,ncel
313  w0(iel) = xro(iel)*vel(1,iel)*vel(2,iel)
314enddo
315call les_filter(1, w0, w4)
316
317! UW
318do iel = 1,ncel
319  w0(iel) = xro(iel)*vel(1,iel)*vel(3,iel)
320enddo
321call les_filter(1, w0, w5)
322
323! VW
324do iel = 1,ncel
325  w0(iel) = xro(iel)*vel(2,iel)*vel(3,iel)
326enddo
327call les_filter(1, w0, w6)
328
329! U
330do iel = 1,ncel
331  w0(iel) = xro(iel)*vel(1,iel)
332enddo
333call les_filter(1, w0, w7)
334do iel = 1, ncel
335  w7(iel) = w7(iel)/xrof(iel)
336enddo
337
338! V
339do iel = 1,ncel
340  w0(iel) = xro(iel)*vel(2,iel)
341enddo
342call les_filter(1, w0, w8)
343do iel = 1, ncel
344  w8(iel) = w8(iel)/xrof(iel)
345enddo
346
347! W
348do iel = 1,ncel
349  w0(iel) = xro(iel)*vel(3,iel)
350enddo
351call les_filter(1, w0, w9)
352do iel = 1, ncel
353  w9(iel) = w9(iel)/xrof(iel)
354enddo
355
356do iel = 1, ncel
357
358  ! Calculation of Lij
359  xl11 = w1(iel) - xrof(iel) * w7(iel) * w7(iel)
360  xl22 = w2(iel) - xrof(iel) * w8(iel) * w8(iel)
361  xl33 = w3(iel) - xrof(iel) * w9(iel) * w9(iel)
362  xl12 = w4(iel) - xrof(iel) * w7(iel) * w8(iel)
363  xl13 = w5(iel) - xrof(iel) * w7(iel) * w9(iel)
364  xl23 = w6(iel) - xrof(iel) * w8(iel) * w9(iel)
365
366  xm11 = xmij(1,iel)
367  xm22 = xmij(2,iel)
368  xm33 = xmij(3,iel)
369  xm12 = xmij(4,iel)
370  xm13 = xmij(5,iel)
371  xm23 = xmij(6,iel)
372  ! Calculation of Mij :: Lij
373  w1(iel) = xm11 * xl11 + 2.d0* xm12 * xl12 + 2.d0* xm13 * xl13  &
374                        +       xm22 * xl22 + 2.d0* xm23 * xl23  &
375                        +                           xm33 * xl33
376  ! Calculation of Mij :: Mij
377  w2(iel) = xm11 * xm11 + 2.d0* xm12 * xm12 + 2.d0* xm13 * xm13  &
378                        +       xm22 * xm22 + 2.d0* xm23 * xm23  &
379                        +                           xm33 * xm33
380
381enddo
382
383deallocate(xmij)
384
385if (irangp.ge.0.or.iperio.eq.1) then
386  call synsca(w1)
387  call synsca(w2)
388endif
389
390! By default we make a local average of numerator and of
391! denominator, then only we make the quotient.
392! The user can do otherwise in ussmag.
393
394call les_filter(1, w1, w3)
395
396call les_filter(1, w2, w4)
397
398do iel = 1, ncel
399  if(abs(w4(iel)).le.epzero) then
400    cpro_smago(iel) = xsmgmx
401  else
402    cpro_smago(iel) = w3(iel)/w4(iel)
403  endif
404enddo
405
406call ussmag                                                       &
407 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
408   icepdc , icetsm , itypsm ,                                     &
409   dt     ,                                                       &
410   ckupdc , smacel ,                                              &
411   w1     , w2     )
412
413iclipc = 0
414do iel = 1, ncel
415  if(cpro_smago(iel).ge.xsmgmx) then
416    cpro_smago(iel) = xsmgmx
417    iclipc = iclipc + 1
418  elseif(cpro_smago(iel).le.xsmgmn) then
419    cpro_smago(iel) = xsmgmn
420    iclipc = iclipc + 1
421  endif
422enddo
423
424!===============================================================================
425! 5.  Calculation of (dynamic) viscosity
426!===============================================================================
427
428do iel = 1, ncel
429  coef = cpro_smago(iel)
430  delta  = xfil * (xa*volume(iel))**xb
431  visct(iel) = crom(iel) * coef * delta**2 * s_n(iel)
432enddo
433
434call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt)
435
436!     Some printings
437if (vcopt%iwarni.ge.1) then
438
439  smagma = -1.0d12
440  smagmi =  1.0d12
441  smagmy =  0.d0
442  do iel = 1, ncel
443    smagma = max(smagma,cpro_smago(iel))
444    smagmi = min(smagmi,cpro_smago(iel))
445    smagmy = smagmy + cpro_smago(iel)*volume(iel)
446  enddo
447  if (irangp.ge.0) then
448    call parmax(smagma)
449    call parmin(smagmi)
450    call parsom(smagmy)
451    call parcpt(iclipc)
452  endif
453  smagmy = smagmy / voltot
454  write(nfecra,1000) iclipc
455  write(nfecra,2001)
456  write(nfecra,2002) smagmy, smagmi, smagma
457  write(nfecra,2003)
458
459endif
460
461!===============================================================================
462! 6.  Clipping of the turbulent viscosity
463!===============================================================================
464
465! Pour la LES en modele dynamique on clippe la viscosite turbulente de maniere
466! a ce que mu_t soit positif,
467
468iclipc = 0
469do iel = 1, ncel
470  if (visct(iel).lt.0.d0) then
471    visct(iel) = 0.d0
472    iclipc = iclipc + 1
473  endif
474enddo
475if (vcopt%iwarni.ge.1) then
476  if (irangp.ge.0) then
477    call parcpt(iclipc)
478  endif
479  write(nfecra,2004) iclipc
480endif
481
482!===============================================================================
483! 7.  Scalar turbulent model
484!===============================================================================
485! In case of gaz combustion, the SGS scalar flux constant and the turbulent
486! diffusivity are only evaluated with the mixture fraction, then applied
487! automatically to the other scalar equations
488
489call field_get_key_id("turbulent_diffusivity_id", key_turb_diff)
490call field_get_key_id("sgs_scalar_flux_coef_id", key_sgs_sca_coef)
491
492do jj = 1, nscal
493  ivar = isca(jj)
494
495  ! For variance of a scalar, the turbulent diffsivity must not be computed
496  if (iscavr(jj).gt.0) cycle
497
498  ! For any scalar other than the mixture fraction in diffusion flames,
499  ! Dt is not computed either. TODO Soot may be an exception
500  if (ippmod(icod3p).ge.0) then
501    if (jj.ne.ifm)  cycle
502  endif
503
504  call field_get_key_int(ivarfl(ivar), key_turb_diff, t_dif_id)
505  call field_get_key_int(ivarfl(ivar), key_sgs_sca_coef, sca_dync_id)
506
507  if (t_dif_id.ge.0.and.sca_dync_id.ge.0) then
508    call field_get_val_s(t_dif_id, cpro_turb_diff)
509    call field_get_val_s(sca_dync_id, cpro_sca_dync)
510
511    !================================================================
512    ! 7.1.  Compute the Mi for scalar
513    !================================================================
514
515    allocate(grads(3, ncelet))
516    allocate(gradsf(3, ncelet))
517
518    inc = 1
519    iprev = 0
520    call field_get_coefa_s(ivarfl(ivar), coefas)
521    call field_get_coefb_s(ivarfl(ivar), coefbs)
522
523    call field_get_val_s(ivarfl(ivar),cvar_sca)
524    call field_gradient_scalar(ivarfl(ivar), iprev, imrgra, inc, 0, grads)
525
526    ! compute grad (<rho.Y>/<rho>)
527    allocate(scami(3,ncelet))
528    allocate(scamif(3,ncelet))
529
530    do iel = 1, ncel
531      w0(iel) = cvar_sca(iel)*xro(iel)
532    enddo
533    call les_filter(1, w0, w4)
534    do iel = 1, ncel
535      w4(iel) = w4(iel)/xrof(iel)
536    enddo
537
538    call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt)
539    inc = 1
540    nswrgp = vcopt%nswrgr
541    epsrgp = vcopt%epsrgr
542    imligp = vcopt%imligr
543    iwarnp = vcopt%iwarni
544    climgp = vcopt%climgr
545    iccocg = 1
546
547    call gradient_s                                                   &
548      (-1     , imrgra , inc    , iccocg , nswrgp , imligp ,          &
549      iwarnp  , epsrgp , climgp ,                                     &
550      w4      , coefas , coefbs ,                                     &
551      gradsf   )
552
553    do iel = 1, ncel
554      delta  = xfil * (xa*volume(iel))**xb
555      do ii = 1, 3
556        scami(ii,iel) = -xro(iel)*delta**2*s_n(iel)*grads(ii,iel)
557      enddo
558    enddo
559
560    call les_filter(3, scami, scamif)
561    do iel = 1, ncel
562      deltaf = xfil2*xfil * (xa*volume(iel))**xb
563      do ii = 1, 3
564        scami(ii,iel) = -deltaf**2*xrof(iel)*sf_n(iel)*gradsf(ii,iel)  &
565                      - scaMif(ii,iel)
566      enddo
567    enddo
568
569    !================================================================
570    ! 7.2.  Compute the Li for scalar
571    !================================================================
572    ! rho*U*Y
573    do iel = 1, ncel
574      w0(iel) = xro(iel)*vel(1,iel)*cvar_sca(iel)
575    enddo
576    call les_filter(1, w0, w1)
577
578    ! rho*V*Y
579    do iel = 1, ncel
580      w0(iel) = xro(iel)*vel(2,iel)*cvar_sca(iel)
581    enddo
582    call les_filter(1, w0, w2)
583
584    ! rho*W*Y
585    do iel = 1, ncel
586      w0(iel) = xro(iel)*vel(3,iel)*cvar_sca(iel)
587    enddo
588    call les_filter(1, w0, w3)
589
590    do iel = 1, ncel
591      scal1 = w1(iel) - xrof(iel)*w7(iel)*w4(iel)
592      scal2 = w2(iel) - xrof(iel)*w8(iel)*w4(iel)
593      scal3 = w3(iel) - xrof(iel)*w9(iel)*w4(iel)
594
595      w1(iel) = scal1*scami(1,iel) + scal2*scami(2,iel) + scal3*scami(3,iel)
596      w2(iel) = scami(1,iel)**2 + scami(2,iel)**2 + scami(3,iel)**2
597    enddo
598
599    if (irangp.ge.0.or.iperio.eq.1) then
600      call synsca(w1)
601      call synsca(w2)
602    endif
603
604    call les_filter(1, w1, w3)
605    call les_filter(1, w2, w4)
606
607    !================================================================
608    ! 7.3.  Compute the SGS flux coefficient and SGS diffusivity
609    !       Cs >= 0, Dt >=0
610    !================================================================
611
612    do iel = 1, ncel
613      if(abs(w4(iel)).le.epzero) then
614        cpro_sca_dync(iel) = 0.d0
615      else
616        cpro_sca_dync(iel) = max(w3(iel)/w4(iel), 0.d0)
617      endif
618
619      delta  = xfil * (xa*volume(iel))**xb
620      cpro_turb_diff(iel) = crom(iel) * cpro_sca_dync(iel) * delta**2 * s_n(iel)
621    enddo
622
623    deallocate(scami, scamif)
624    deallocate(grads, gradsf)
625  endif
626
627enddo
628
629
630! Free memory
631deallocate(s_n, sf_n)
632deallocate(w9, w8)
633deallocate(w7, w6, w5)
634deallocate(w4, w3, w2)
635deallocate(w1, w0)
636deallocate(xro, xrof)
637
638!----
639! Formats
640!----
641
642 1000 format(                                                           &
643' Nb of clipping of the Smagorinsky constant',              I10,/)
644 2001 format(                                                           &
645' --- Informations on the squared Smagorinsky constant'        ,/,&
646' --------------------------------'                            ,/,&
647' Mean value  Min value  Max value'                            ,/,&
648' --------------------------------'                              )
649 2002 format(                                                           &
650 e12.4    ,      e12.4,      e12.4                               )
651 2003 format(                                                           &
652' --------------------------------'                            ,/)
653 2004 format(                                                           &
654' Nb of clippings for the turbulent viscosity (mu_t>0):',   i10,/)
655
656!----
657! End
658!----
659
660return
661end subroutine
662