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 resrij2.f90
28!>
29!> \brief This subroutine prepares the solving of the coupled Reynolds stress
30!> components in \f$ R_{ij} - \varepsilon \f$ RANS (LRR) turbulence model.
31!>
32!> \remark
33!> - cvar_var(1,*) for \f$ R_{11} \f$
34!> - cvar_var(2,*) for \f$ R_{22} \f$
35!> - cvar_var(3,*) for \f$ R_{33} \f$
36!> - cvar_var(4,*) for \f$ R_{12} \f$
37!> - cvar_var(5,*) for \f$ R_{23} \f$
38!> - cvar_var(6,*) for \f$ R_{13} \f$
39!-------------------------------------------------------------------------------
40
41!-------------------------------------------------------------------------------
42! Arguments
43!______________________________________________________________________________.
44!  mode           name          role
45!______________________________________________________________________________!
46!> \param[in]     ivar          variable number
47!> \param[in]     gradv         work array for the velocity grad term
48!>                                 only for iturb=31
49!> \param[in]     produc        work array for production
50!> \param[in]     gradro        work array for grad rom
51!>                              (without rho volume) only for iturb=30
52!> \param[in]     viscf         visc*surface/dist at internal faces
53!> \param[in]     viscb         visc*surface/dist at edge faces
54!> \param[out]    viscce        Daly Harlow diffusion term
55!> \param[in]     smbr          working array
56!> \param[in]     rovsdt        working array
57!> \param[out]    weighf        working array
58!> \param[out]    weighb        working array
59!_______________________________________________________________________________
60
61subroutine resrij2 &
62 ( ivar   ,                                                       &
63   gradv  , produc , gradro ,                                     &
64   viscf  , viscb  , viscce ,                                     &
65   smbr   , rovsdt ,                                              &
66   weighf , weighb )
67
68!===============================================================================
69
70!===============================================================================
71! Module files
72!===============================================================================
73
74use, intrinsic :: iso_c_binding
75
76use paramx
77use numvar
78use entsor
79use optcal
80use cstphy
81use cstnum
82use parall
83use period
84use lagran
85use mesh
86use field
87use cs_f_interfaces
88use rotation
89use turbomachinery
90use cs_c_bindings
91
92!===============================================================================
93
94implicit none
95
96! Arguments
97
98integer          ivar
99
100double precision gradv(3, 3, ncelet)
101double precision produc(6, ncelet)
102double precision gradro(3, ncelet)
103double precision viscf(nfac), viscb(nfabor), viscce(6, ncelet)
104double precision smbr(6, ncelet)
105double precision rovsdt(6, 6, ncelet)
106double precision weighf(2, nfac), weighb(nfabor)
107
108! Local variables
109
110integer          iel
111integer          isou, jsou
112integer          ii    , jj    , kk
113integer          iflmas, iflmab
114integer          iwarnp
115integer          imvisp
116integer          st_prv_id
117integer          t2v(3,3)
118integer          iv2t(6), jv2t(6)
119integer          key_t_ext_id, f_id, rot_id
120integer          iroext
121
122double precision trprod, trrij
123double precision deltij(6), dij(3,3)
124double precision thets , thetv
125double precision matrot(3,3)
126double precision d1s2, d1s3, d2s3
127double precision ccorio
128double precision rctse
129double precision, dimension(3,3) :: cvara_r
130
131character(len=80) :: label
132double precision, allocatable, dimension(:,:), target :: buoyancy
133double precision, allocatable, dimension(:) :: w1
134double precision, allocatable, dimension(:,:) :: w2
135double precision, dimension(:), pointer :: imasfl, bmasfl
136double precision, dimension(:), pointer :: crom, cromo
137double precision, dimension(:,:), pointer :: visten
138double precision, dimension(:), pointer :: cvara_ep
139double precision, dimension(:,:), pointer :: cvar_var, cvara_var
140double precision, dimension(:), pointer :: viscl
141double precision, dimension(:,:), pointer :: c_st_prv
142double precision, dimension(:,:), pointer :: cpro_buoyancy
143
144integer iii, jjj
145
146double precision impl_drsm(6,6)
147double precision implmat2add(3,3)
148double precision impl_lin_cst, impl_id_cst
149double precision aiksjk, aikrjk, aii ,aklskl, aikakj
150double precision xaniso(3,3), xstrai(3,3), xrotac(3,3), xprod(3,3)
151double precision sym_strain(6)
152double precision matrn(6), oo_matrn(6)
153double precision eigen_vals(3)
154double precision ceps_impl
155double precision eigen_max
156double precision pij, phiij1, phiij2, epsij
157
158type(var_cal_opt) :: vcopt
159
160!===============================================================================
161
162!===============================================================================
163! Initialization
164!===============================================================================
165
166! Time extrapolation?
167call field_get_key_id("time_extrapolated", key_t_ext_id)
168
169call field_get_key_int(icrom, key_t_ext_id, iroext)
170
171! Allocate work arrays
172allocate(w1(ncelet))
173allocate(w2(6,ncelet))
174
175! Generating the tensor to vector (t2v) and vector to tensor (v2t) mask arrays
176
177! a) t2v
178t2v(1,1) = 1; t2v(1,2) = 4; t2v(1,3) = 6;
179t2v(2,1) = 4; t2v(2,2) = 2; t2v(2,3) = 5;
180t2v(3,1) = 6; t2v(3,2) = 5; t2v(3,3) = 3;
181
182! b) i index of v2t
183iv2t(1) = 1; iv2t(2) = 2; iv2t(3) = 3;
184iv2t(4) = 1; iv2t(5) = 2; iv2t(6) = 1;
185
186! c) j index of v2t
187jv2t(1) = 1; jv2t(2) = 2; jv2t(3) = 3;
188jv2t(4) = 2; jv2t(5) = 3; jv2t(6) = 3;
189
190! d) kronecker symbol
191dij(:,:) = 0.0d0;
192dij(1,1) = 1.0d0; dij(2,2) = 1.0d0; dij(3,3) = 1.0d0;
193
194call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt)
195
196if (vcopt%iwarni.ge.1) then
197  call field_get_label(ivarfl(ivar), label)
198  write(nfecra,1000) label
199endif
200
201call field_get_val_s(icrom, crom)
202call field_get_val_s(iviscl, viscl)
203
204call field_get_val_prev_s(ivarfl(iep), cvara_ep)
205
206call field_get_val_v(ivarfl(ivar), cvar_var)
207call field_get_val_prev_v(ivarfl(ivar), cvara_var)
208
209call field_get_key_int(ivarfl(iu), kimasf, iflmas)
210call field_get_key_int(ivarfl(iu), kbmasf, iflmab)
211call field_get_val_s(iflmas, imasfl)
212call field_get_val_s(iflmab, bmasfl)
213
214d1s2 = 1.d0/2.d0
215d1s3 = 1.d0/3.d0
216d2s3 = 2.d0/3.d0
217
218do isou = 1, 3
219  deltij(isou) = 1.0d0
220enddo
221do isou = 4, 6
222  deltij(isou) = 0.0d0
223enddo
224
225!     S as Source, V as Variable
226thets  = thetst
227thetv  = vcopt%thetav
228
229call field_get_key_int(ivarfl(ivar), kstprv, st_prv_id)
230if (st_prv_id .ge. 0) then
231  call field_get_val_v(st_prv_id, c_st_prv)
232else
233  c_st_prv=> null()
234endif
235
236if (st_prv_id.ge.0.and.iroext.gt.0) then
237  call field_get_val_prev_s(icrom, cromo)
238else
239  call field_get_val_s(icrom, cromo)
240endif
241
242! Coefficient of the "Coriolis-type" term
243if (icorio.eq.1) then
244  ! Relative velocity formulation
245  ccorio = 2.d0
246elseif (iturbo.eq.1) then
247  ! Mixed relative/absolute velocity formulation
248  ccorio = 1.d0
249else
250  ccorio = 0.d0
251endif
252
253!===============================================================================
254! Production, Pressure-Strain correlation, dissipation
255!===============================================================================
256
257do iel = 1, ncel
258
259  ! Initalize implicit matrices at 0
260  do isou = 1, 6
261    do jsou = 1, 6
262      impl_drsm(isou, jsou) = 0.0d0
263    end do
264  end do
265  do isou = 1, 3
266    do jsou = 1, 3
267      implmat2add(isou, jsou) = 0.0d0
268    end do
269  end do
270
271  impl_lin_cst = 0.0d0
272  impl_id_cst  = 0.0d0
273
274  ! Pij
275  xprod(1,1) = produc(1, iel)
276  xprod(1,2) = produc(4, iel)
277  xprod(1,3) = produc(6, iel)
278  xprod(2,2) = produc(2, iel)
279  xprod(2,3) = produc(5, iel)
280  xprod(3,3) = produc(3, iel)
281
282  xprod(2,1) = xprod(1,2)
283  xprod(3,1) = xprod(1,3)
284  xprod(3,2) = xprod(2,3)
285
286  trprod = d1s2 * (xprod(1,1) + xprod(2,2) + xprod(3,3) )
287  trrij  = d1s2 * (cvara_var(1,iel) + cvara_var(2,iel) + cvara_var(3,iel))
288
289  !-----> aII = aijaij
290  aii    = 0.d0
291  aklskl = 0.d0
292  aiksjk = 0.d0
293  aikrjk = 0.d0
294  aikakj = 0.d0
295
296  ! aij
297  xaniso(1,1) = cvara_var(1,iel)/trrij - d2s3
298  xaniso(2,2) = cvara_var(2,iel)/trrij - d2s3
299  xaniso(3,3) = cvara_var(3,iel)/trrij - d2s3
300  xaniso(1,2) = cvara_var(4,iel)/trrij
301  xaniso(1,3) = cvara_var(6,iel)/trrij
302  xaniso(2,3) = cvara_var(5,iel)/trrij
303  xaniso(2,1) = xaniso(1,2)
304  xaniso(3,1) = xaniso(1,3)
305  xaniso(3,2) = xaniso(2,3)
306
307  ! Sij
308  xstrai(1,1) = gradv(1, 1, iel)
309  xstrai(1,2) = d1s2*(gradv(2, 1, iel)+gradv(1, 2, iel))
310  xstrai(1,3) = d1s2*(gradv(3, 1, iel)+gradv(1, 3, iel))
311  xstrai(2,1) = xstrai(1,2)
312  xstrai(2,2) = gradv(2, 2, iel)
313  xstrai(2,3) = d1s2*(gradv(3, 2, iel)+gradv(2, 3, iel))
314  xstrai(3,1) = xstrai(1,3)
315  xstrai(3,2) = xstrai(2,3)
316  xstrai(3,3) = gradv(3, 3, iel)
317
318  ! omegaij
319  xrotac(1,1) = 0.d0
320  xrotac(1,2) = d1s2*(gradv(2, 1, iel)-gradv(1, 2, iel))
321  xrotac(1,3) = d1s2*(gradv(3, 1, iel)-gradv(1, 3, iel))
322  xrotac(2,1) = -xrotac(1,2)
323  xrotac(2,2) = 0.d0
324  xrotac(2,3) = d1s2*(gradv(3, 2, iel)-gradv(2, 3, iel))
325  xrotac(3,1) = -xrotac(1,3)
326  xrotac(3,2) = -xrotac(2,3)
327  xrotac(3,3) = 0.d0
328
329  do ii=1,3
330    do jj = 1,3
331      ! aii = aij.aij
332      aii    = aii+xaniso(ii,jj)*xaniso(ii,jj)
333      ! aklskl = aij.Sij
334      aklskl = aklskl + xaniso(ii,jj)*xstrai(ii,jj)
335    enddo
336  enddo
337
338  ! Computation of implicit components
339
340  sym_strain(1) = xstrai(1,1)
341  sym_strain(2) = xstrai(2,2)
342  sym_strain(3) = xstrai(3,3)
343  sym_strain(4) = xstrai(1,2)
344  sym_strain(5) = xstrai(2,3)
345  sym_strain(6) = xstrai(1,3)
346
347  do isou = 1, 6
348    matrn(isou) = cvara_var(isou,iel)/trrij
349    oo_matrn(isou) = 0.0d0
350  end do
351
352  ! Inversing the matrix
353  call symmetric_matrix_inverse(matrn, oo_matrn)
354  do isou = 1, 6
355    oo_matrn(isou) = oo_matrn(isou)/trrij
356  end do
357
358  ! Computing the maximal eigenvalue (in terms of norm!) of S
359  call calc_symtens_eigvals(sym_strain, eigen_vals)
360  eigen_max = maxval(abs(eigen_vals))
361
362  ! Constant for the dissipation
363  ceps_impl = d1s3 * cvara_ep(iel)
364
365  ! Identity constant
366  impl_id_cst = -d1s3*crij2*min(trprod,0.0d0)
367
368  ! Linear constant
369  impl_lin_cst = eigen_max *     ( &
370                 (1.d0-crij2)    )   ! Production + Phi2
371
372  do jsou = 1, 3
373    do isou = 1 ,3
374      iii = t2v(isou,jsou)
375      implmat2add(isou,jsou) = (1.d0-crij2)*xrotac(isou,jsou)    &
376                             + impl_lin_cst*deltij(iii)          &
377                             + impl_id_cst*d1s2*oo_matrn(iii)    &
378                             + ceps_impl*oo_matrn(iii)
379    end do
380  end do
381
382  impl_drsm(:,:) = 0.0d0
383  call reduce_symprod33_to_66(implmat2add, impl_drsm)
384
385  ! Rotating frame of reference => "absolute" vorticity
386  if (icorio.eq.1) then
387    do ii = 1, 3
388      do jj = 1, 3
389        xrotac(ii,jj) = xrotac(ii,jj) + matrot(ii,jj)
390      enddo
391    enddo
392  endif
393
394  do isou = 1, 6
395    iii = iv2t(isou)
396    jjj = jv2t(isou)
397    aiksjk = 0
398    aikrjk = 0
399    aikakj = 0
400    do kk = 1,3
401      ! aiksjk = aik.Sjk+ajk.Sik
402      aiksjk =   aiksjk + xaniso(iii,kk)*xstrai(jjj,kk)              &
403               + xaniso(jjj,kk)*xstrai(iii,kk)
404      ! aikrjk = aik.Omega_jk + ajk.omega_ik
405      aikrjk =   aikrjk + xaniso(iii,kk)*xrotac(jjj,kk)              &
406               + xaniso(jjj,kk)*xrotac(iii,kk)
407      ! aikakj = aik*akj
408      aikakj = aikakj + xaniso(iii,kk)*xaniso(kk,jjj)
409    enddo
410
411    ! Explicit terms
412    pij = (1.d0 - crij2)*xprod(iii,jjj)
413    phiij1 = -cvara_ep(iel)*crij1*xaniso(iii,jjj)
414    phiij2 = d2s3*crij2*trprod*deltij(isou)
415    epsij = -d2s3*cvara_ep(iel)*deltij(isou)
416
417    if (st_prv_id.ge.0) then
418      c_st_prv(isou,iel) = c_st_prv(isou,iel) &
419        + cromo(iel)*cell_f_vol(iel)*(pij+phiij1+phiij2+epsij)
420    else
421      smbr(isou,iel) = smbr(isou,iel) &
422        + cromo(iel)*cell_f_vol(iel)*(pij+phiij1+phiij2+epsij)
423      ! Implicit terms
424      rovsdt(isou,isou,iel) = rovsdt(isou,isou,iel) &
425        + cell_f_vol(iel)/trrij*crom(iel)*(crij1*cvara_ep(iel))
426
427      ! Careful ! Inversion of the order of the coefficients since
428      ! rovsdt matrix is then used by a c function for the linear solving
429      do jsou = 1, 6
430        rovsdt(jsou,isou,iel) = rovsdt(jsou,isou,iel) + cell_f_vol(iel) &
431                                *crom(iel) * impl_drsm(isou,jsou)
432      end do
433
434    endif
435
436  enddo
437
438enddo
439
440!===============================================================================
441! Coriolis terms in the Phi1 and production
442!===============================================================================
443
444if (icorio.eq.1 .or. iturbo.eq.1) then
445
446  do iel = 1, ncel
447
448    rot_id = icorio
449    if (iturbo.eq.1) rot_id = irotce(iel)
450
451    if (rot_id .lt. 1) cycle
452
453    call coriolis_t(rot_id, 1.d0, matrot)
454
455    cvara_r(1,1) = cvara_var(1,iel)
456    cvara_r(2,2) = cvara_var(2,iel)
457    cvara_r(3,3) = cvara_var(3,iel)
458    cvara_r(1,2) = cvara_var(4,iel)
459    cvara_r(2,3) = cvara_var(5,iel)
460    cvara_r(1,3) = cvara_var(6,iel)
461    cvara_r(2,1) = cvara_var(4,iel)
462    cvara_r(3,2) = cvara_var(5,iel)
463    cvara_r(3,1) = cvara_var(6,iel)
464
465    ! Compute Gij: (i,j) component of the Coriolis production
466    do isou = 1, 6
467      ii = iv2t(isou)
468      jj = jv2t(isou)
469
470      w2(isou,iel) = 0.d0
471      do kk = 1, 3
472        w2(isou,iel) = w2(isou,iel) - ccorio*(  matrot(ii,kk)*cvara_r(jj,kk) &
473                                    + matrot(jj,kk)*cvara_r(ii,kk) )
474      enddo
475    enddo
476  enddo
477
478  ! Coriolis contribution in the Phi1 term: (1-C2/2)Gij
479  if (icorio.eq.1) then
480    do iel = 1, ncel
481      do isou = 1, 6
482        w2(isou,iel) = crom(iel)*cell_f_vol(iel)*(1.d0 - 0.5d0*crij2)*w2(isou,iel)
483      enddo
484    enddo
485  endif
486
487  ! If source terms are extrapolated
488  if (st_prv_id.ge.0) then
489    do iel = 1, ncel
490      do isou = 1, 6
491        c_st_prv(isou,iel) = c_st_prv(isou,iel) + w2(isou,iel)
492      enddo
493    enddo
494  ! Otherwise, directly in smbr
495  else
496    do iel = 1, ncel
497      do isou = 1, 6
498        smbr(isou,iel) = smbr(isou,iel) + w2(isou,iel)
499      enddo
500    enddo
501  endif
502
503endif
504
505!===============================================================================
506! Wall echo terms
507!===============================================================================
508
509if (irijec.eq.1) then !todo
510
511  do iel = 1, ncel
512    do isou = 1, 6
513      w2(isou,iel) = 0.d0
514    enddo
515  enddo
516
517  call rijech2(produc, w2)
518
519  ! If we extrapolate the source terms: c_st_prv
520  if (st_prv_id.ge.0) then
521    do iel = 1, ncel
522      do isou = 1, 6
523        c_st_prv(isou,iel) = c_st_prv(isou,iel) + w2(isou,iel)
524      enddo
525    enddo
526  ! Otherwise smbr
527  else
528    do iel = 1, ncel
529      do isou = 1, 6
530        smbr(isou,iel) = smbr(isou,iel) + w2(isou,iel)
531      enddo
532    enddo
533  endif
534
535endif
536
537!===============================================================================
538! Buoyancy source term
539!===============================================================================
540
541if (igrari.eq.1) then
542
543  call field_get_id_try("rij_buoyancy", f_id)
544  if (f_id.ge.0) then
545    call field_get_val_v(f_id, cpro_buoyancy)
546  else
547    ! Allocate a work array
548    allocate(buoyancy(6,ncelet))
549    cpro_buoyancy => buoyancy
550  endif
551
552  call rijthe2(gradro, cpro_buoyancy)
553
554  ! If we extrapolate the source terms: previous ST
555  if (st_prv_id.ge.0) then
556    do iel = 1, ncel
557      do isou = 1, 6
558        c_st_prv(isou,iel) = c_st_prv(isou,iel) + cpro_buoyancy(isou,iel) * cell_f_vol(iel)
559      enddo
560    enddo
561  ! Otherwise smbr
562  else
563    do iel = 1, ncel
564      do isou = 1, 6
565        smbr(isou,iel) = smbr(isou,iel) + cpro_buoyancy(isou,iel) * cell_f_vol(iel)
566      enddo
567    enddo
568  endif
569
570  ! Free memory
571  if (allocated(buoyancy)) deallocate(buoyancy)
572
573endif
574
575!===============================================================================
576! Diffusion term (Daly Harlow: generalized gradient hypothesis method)
577!===============================================================================
578
579! Symmetric tensor diffusivity (GGDH)
580if (iand(vcopt%idften, ANISOTROPIC_RIGHT_DIFFUSION).ne.0) then
581
582  call field_get_val_v(ivsten, visten)
583
584  do iel = 1, ncel
585    viscce(1,iel) = visten(1,iel) + viscl(iel)
586    viscce(2,iel) = visten(2,iel) + viscl(iel)
587    viscce(3,iel) = visten(3,iel) + viscl(iel)
588    viscce(4,iel) = visten(4,iel)
589    viscce(5,iel) = visten(5,iel)
590    viscce(6,iel) = visten(6,iel)
591  enddo
592
593  iwarnp = vcopt%iwarni
594
595  call vitens(viscce, iwarnp, weighf, weighb, viscf, viscb)
596
597! Scalar diffusivity
598else
599
600  do iel = 1, ncel
601    trrij = 0.5d0 * (cvara_var(1,iel) + cvara_var(2,iel) + cvara_var(3,iel))
602    rctse = crom(iel) * csrij * trrij**2 / cvara_ep(iel)
603    w1(iel) = viscl(iel) + vcopt%idifft*rctse
604  enddo
605
606  imvisp = vcopt%imvisf
607
608  call viscfa(imvisp, w1, viscf, viscb)
609
610endif
611
612! Free memory
613
614deallocate(w1, w2)
615
616!--------
617! Formats
618!--------
619
620 1000 format(/,'           Solving variable ', a8          ,/)
621
622!----
623! End
624!----
625
626return
627
628end subroutine
629