← Back to branch summary

~ubuntu-branches/ubuntu/precise/code-saturne/precise

« back to all changes in this revision

Viewing changes to src/turb/reseps.f90

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2011-11-01 17:43:32 UTC
  • mto: (6.1.7 sid)
  • mto: This revision was merged to the branch mainline in revision 11.
  • Revision ID: package-import@ubuntu.com-20111101174332-tl4vk45no0x3emc3
Tags: upstream-2.1.0
Import upstream version 2.1.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/turb/reseps.f90
 
1
!-------------------------------------------------------------------------------
 
2
 
 
3
! This file is part of Code_Saturne, a general-purpose CFD tool.
 
4
!
 
5
! Copyright (C) 1998-2011 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
subroutine reseps &
 
24
!================
 
25
 
 
26
 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
 
27
   ivar   , isou   , ipp    ,                                     &
 
28
   icepdc , icetsm , itpsmp ,                                     &
 
29
   dt     , rtp    , rtpa   , propce , propfa , propfb ,          &
 
30
   coefa  , coefb  , grdvit , produc , gradro ,                   &
 
31
   ckupdc , smcelp , gamma  ,                                     &
 
32
   viscf  , viscb  ,                                              &
 
33
   tslagr ,                                                       &
 
34
   smbr   , rovsdt )
 
35
 
 
36
!===============================================================================
 
37
! FONCTION :
 
38
! ----------
 
39
 
 
40
! RESOLUTION DES EQUATIONS CONVECTION DIFFUSION TERME SOURCE
 
41
!   POUR EPSILON DANS LE CAS DU MODELE Rij-epsilon
 
42
 
 
43
!   On a ici ISOU   = 7 (7 ieme variable du Rij-epsilon)
 
44
 
 
45
!-------------------------------------------------------------------------------
 
46
!ARGU                             ARGUMENTS
 
47
!__________________.____._____.________________________________________________.
 
48
! name             !type!mode ! role                                           !
 
49
!__________________!____!_____!________________________________________________!
 
50
! nvar             ! i  ! <-- ! total number of variables                      !
 
51
! nscal            ! i  ! <-- ! total number of scalars                        !
 
52
! ncepdp           ! i  ! <-- ! number of cells with head loss                 !
 
53
! ncesmp           ! i  ! <-- ! number of cells with mass source term          !
 
54
! ivar             ! i  ! <-- ! variable number                                !
 
55
! isou             ! e  ! <-- ! numero de passage                              !
 
56
! ipp              ! e  ! <-- ! numero de variable pour sorties post           !
 
57
! icepdc(ncelet    ! te ! <-- ! numero des ncepdp cellules avec pdc            !
 
58
! icetsm(ncesmp    ! te ! <-- ! numero des cellules a source de masse          !
 
59
! itpsmp           ! te ! <-- ! type de source de masse pour la                !
 
60
! (ncesmp)         !    !     !  variables (cf. ustsma)                        !
 
61
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
 
62
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
 
63
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
 
64
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
 
65
! propfa(nfac, *)  ! ra ! <-- ! physical properties at interior face centers   !
 
66
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers   !
 
67
! coefa, coefb     ! ra ! <-- ! boundary conditions                            !
 
68
!  (nfabor, *)     !    !     !                                                !
 
69
! grdvit           ! tr ! --- ! tableau de travail pour terme grad             !
 
70
!  (ncelet,3,3)    !    !     !    de vitesse     uniqt pour iturb=31          !
 
71
! produc           ! tr ! <-- ! tableau de travail pour production             !
 
72
!  (6,ncelet)      !    !     ! (sans rho volume) uniqt pour iturb=30          !
 
73
! gradro(ncelet,3) ! tr ! <-- ! tableau de travail pour grad rom               !
 
74
! ckupdc           ! tr ! <-- ! tableau de travail pour pdc                    !
 
75
!  (ncepdp,6)      !    !     !                                                !
 
76
! smcelp(ncesmp    ! tr ! <-- ! valeur de la variable associee a la            !
 
77
!                  !    !     !  source de masse                               !
 
78
! gamma(ncesmp)    ! tr ! <-- ! valeur du flux de masse                        !
 
79
! viscf(nfac)      ! tr ! --- ! visc*surface/dist aux faces internes           !
 
80
! viscb(nfabor     ! tr ! --- ! visc*surface/dist aux faces de bord            !
 
81
! tslagr           ! tr ! <-- ! terme de couplage retour du                    !
 
82
!  (ncelet,*)      !    !     !   lagrangien                                   !
 
83
! smbr(ncelet      ! tr ! --- ! tableau de travail pour sec mem                !
 
84
! rovsdt(ncelet    ! tr ! --- ! tableau de travail pour terme instat           !
 
85
!__________________!____!_____!________________________________________________!
 
86
 
 
87
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
 
88
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
 
89
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
 
90
!            --- tableau de travail
 
91
!===============================================================================
 
92
 
 
93
!===============================================================================
 
94
! Module files
 
95
!===============================================================================
 
96
 
 
97
use paramx
 
98
use dimens, only: ndimfb
 
99
use numvar
 
100
use entsor
 
101
use optcal
 
102
use cstnum
 
103
use cstphy
 
104
use parall
 
105
use period
 
106
use lagran
 
107
use mesh
 
108
 
 
109
!===============================================================================
 
110
 
 
111
implicit none
 
112
 
 
113
! Arguments
 
114
 
 
115
integer          nvar   , nscal
 
116
integer          ncepdp , ncesmp
 
117
integer          ivar   , isou   , ipp
 
118
 
 
119
integer          icepdc(ncepdp)
 
120
integer          icetsm(ncesmp), itpsmp(ncesmp)
 
121
 
 
122
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
 
123
double precision propce(ncelet,*)
 
124
double precision propfa(nfac,*), propfb(ndimfb,*)
 
125
double precision coefa(ndimfb,*), coefb(ndimfb,*)
 
126
double precision produc(6,ncelet), grdvit(ncelet,3,3)
 
127
double precision gradro(ncelet,3)
 
128
double precision ckupdc(ncepdp,6)
 
129
double precision smcelp(ncesmp), gamma(ncesmp)
 
130
double precision viscf(nfac), viscb(nfabor)
 
131
double precision tslagr(ncelet,*)
 
132
double precision smbr(ncelet), rovsdt(ncelet)
 
133
 
 
134
! Local variables
 
135
 
 
136
integer          init  , ifac  , iel   , inc   , iccocg
 
137
integer          ii    ,jj     , iiun
 
138
integer          iclvar, iclvaf
 
139
integer          ipcrom, ipcroo, ipcvis, ipcvst, iflmas, iflmab
 
140
integer          nswrgp, imligp, iwarnp
 
141
integer          iconvp, idiffp, ndircp, ireslp
 
142
integer          nitmap, nswrsp, ircflp, ischcp, isstpp, iescap
 
143
integer          imgrp , ncymxp, nitmfp
 
144
integer          iptsta
 
145
 
 
146
double precision blencp, epsilp, epsrgp, climgp, extrap, relaxp
 
147
double precision epsrsp
 
148
double precision trprod , trrij ,csteps, rctse
 
149
double precision grdpx,grdpy,grdpz,grdsn
 
150
double precision surfn2
 
151
double precision tseps , kseps , ceps2
 
152
double precision tuexpe, thets , thetv , thetap, thetp1
 
153
 
 
154
double precision rvoid(1)
 
155
 
 
156
double precision, allocatable, dimension(:,:) :: grad
 
157
double precision, allocatable, dimension(:) :: w1, w2, w3
 
158
double precision, allocatable, dimension(:) :: w4, w5, w6
 
159
double precision, allocatable, dimension(:) :: w7, w8, w9
 
160
 
 
161
!===============================================================================
 
162
 
 
163
!===============================================================================
 
164
! 1. INITIALISATION
 
165
!===============================================================================
 
166
 
 
167
! Allocate work arrays
 
168
allocate(w1(ncelet), w2(ncelet), w3(ncelet))
 
169
allocate(w4(ncelet), w5(ncelet), w6(ncelet))
 
170
allocate(w8(ncelet), w9(ncelet))
 
171
 
 
172
 
 
173
if(iwarni(ivar).ge.1) then
 
174
  write(nfecra,1000) nomvar(ipp)
 
175
endif
 
176
 
 
177
ipcrom = ipproc(irom  )
 
178
ipcvis = ipproc(iviscl)
 
179
ipcvst = ipproc(ivisct)
 
180
iflmas = ipprof(ifluma(iu))
 
181
iflmab = ipprob(ifluma(iu))
 
182
 
 
183
iclvar = iclrtp(ivar,icoef)
 
184
iclvaf = iclrtp(ivar,icoeff)
 
185
 
 
186
!     Constante Ce2, qui vaut CE2 pour ITURB=30 et CSSGE2 pour ITRUB=31
 
187
if (iturb.eq.30) then
 
188
  ceps2 = ce2
 
189
else
 
190
  ceps2 = cssge2
 
191
endif
 
192
 
 
193
!     S pour Source, V pour Variable
 
194
thets  = thetst
 
195
thetv  = thetav(ivar )
 
196
 
 
197
ipcroo = ipcrom
 
198
if(isto2t.gt.0.and.iroext.gt.0) then
 
199
  ipcroo = ipproc(iroma)
 
200
endif
 
201
if(isto2t.gt.0) then
 
202
  iptsta = ipproc(itstua)
 
203
else
 
204
  iptsta = 0
 
205
endif
 
206
 
 
207
do iel = 1, ncel
 
208
  smbr(iel) = 0.d0
 
209
enddo
 
210
do iel = 1, ncel
 
211
  rovsdt(iel) = 0.d0
 
212
enddo
 
213
 
 
214
!===============================================================================
 
215
! 2. TERMES SOURCES  UTILISATEURS
 
216
!===============================================================================
 
217
 
 
218
call ustsri                                                       &
 
219
!==========
 
220
 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
 
221
   ivar   ,                                                       &
 
222
   icepdc , icetsm , itpsmp ,                                     &
 
223
   dt     , rtpa   , propce , propfa , propfb ,                   &
 
224
   coefa  , coefb  , ckupdc , smcelp , gamma  , grdvit , produc , &
 
225
   smbr   , rovsdt )
 
226
 
 
227
!     Si on extrapole les T.S.
 
228
if(isto2t.gt.0) then
 
229
  do iel = 1, ncel
 
230
!       Sauvegarde pour echange
 
231
    tuexpe = propce(iel,iptsta+isou-1)
 
232
!       Pour la suite et le pas de temps suivant
 
233
    propce(iel,iptsta+isou-1) = smbr(iel)
 
234
!       Second membre du pas de temps precedent
 
235
!       On suppose -ROVSDT > 0 : on implicite
 
236
!          le terme source utilisateur (le reste)
 
237
    smbr(iel) = rovsdt(iel)*rtpa(iel,ivar) - thets*tuexpe
 
238
!       Diagonale
 
239
    rovsdt(iel) = - thetv*rovsdt(iel)
 
240
  enddo
 
241
else
 
242
  do iel = 1, ncel
 
243
    smbr(iel)   = rovsdt(iel)*rtpa(iel,ivar) + smbr(iel)
 
244
    rovsdt(iel) = max(-rovsdt(iel),zero)
 
245
  enddo
 
246
endif
 
247
 
 
248
!===============================================================================
 
249
! 3.  TERMES SOURCES LAGRANGIEN : COUPLAGE RETOUR
 
250
!===============================================================================
 
251
 
 
252
!     Ordre 2 non pris en compte
 
253
 if (iilagr.eq.2 .and. ltsdyn.eq.1) then
 
254
 
 
255
   do iel = 1,ncel
 
256
! Ts sur eps
 
257
     tseps = -0.5d0 * ( tslagr(iel,itsr11)                        &
 
258
                      + tslagr(iel,itsr22)                        &
 
259
                      + tslagr(iel,itsr33) )
 
260
! rapport k/eps
 
261
     kseps = 0.5d0 * ( rtpa(iel,ir11)                           &
 
262
                     + rtpa(iel,ir22)                           &
 
263
                     + rtpa(iel,ir33) )                         &
 
264
                     / rtpa(iel,iep)
 
265
 
 
266
     smbr(iel)   = smbr(iel) + ce4 *tseps *rtpa(iel,iep) /kseps
 
267
     rovsdt(iel) = rovsdt(iel) + max( (-ce4*tseps/kseps) , zero)
 
268
   enddo
 
269
 
 
270
 endif
 
271
 
 
272
!===============================================================================
 
273
! 4. TERME SOURCE DE MASSE
 
274
!===============================================================================
 
275
 
 
276
 
 
277
if (ncesmp.gt.0) then
 
278
 
 
279
!       Entier egal a 1 (pour navsto : nb de sur-iter)
 
280
  iiun = 1
 
281
 
 
282
!       On incremente SMBR par -Gamma RTPA et ROVSDT par Gamma (*theta)
 
283
  call catsma                                                     &
 
284
  !==========
 
285
 ( ncelet , ncel   , ncesmp , iiun   , isto2t , thetv ,    &
 
286
   icetsm , itpsmp ,                                              &
 
287
   volume , rtpa(1,ivar) , smcelp , gamma  ,                      &
 
288
   smbr   ,  rovsdt , w1 )
 
289
 
 
290
!       Si on extrapole les TS on met Gamma Pinj dans PROPCE
 
291
  if(isto2t.gt.0) then
 
292
    do iel = 1, ncel
 
293
      propce(iel,iptsta+isou-1) =                                 &
 
294
      propce(iel,iptsta+isou-1) + w1(iel)
 
295
    enddo
 
296
!       Sinon on le met directement dans SMBR
 
297
  else
 
298
    do iel = 1, ncel
 
299
      smbr(iel) = smbr(iel) + w1(iel)
 
300
    enddo
 
301
  endif
 
302
 
 
303
endif
 
304
 
 
305
!===============================================================================
 
306
! 5. TERME D'ACCUMULATION DE MASSE -(dRO/dt)*VOLUME
 
307
!    ET TERME INSTATIONNAIRE
 
308
!===============================================================================
 
309
 
 
310
! ---> Calcul de mij
 
311
 
 
312
init = 1
 
313
call divmas(ncelet,ncel,nfac,nfabor,init,nfecra,                  &
 
314
               ifacel,ifabor,propfa(1,iflmas),propfb(1,iflmab),w1)
 
315
 
 
316
! ---> Ajout au second membre
 
317
 
 
318
do iel = 1, ncel
 
319
  smbr(iel) = smbr(iel)                                           &
 
320
              + iconv(ivar)*w1(iel)*rtpa(iel,ivar)
 
321
enddo
 
322
 
 
323
! ---> Ajout dans la diagonale de la matrice
 
324
!     Extrapolation ou non, meme forme par coherence avec bilsc2
 
325
 
 
326
do iel = 1, ncel
 
327
  rovsdt(iel) = rovsdt(iel)                                       &
 
328
           + istat(ivar)*(propce(iel,ipcrom)/dt(iel))*volume(iel) &
 
329
           - iconv(ivar)*w1(iel)*thetv
 
330
enddo
 
331
 
 
332
 
 
333
!===============================================================================
 
334
! 6. PRODUCTION RHO * Ce1 * epsilon / k * P
 
335
!    DISSIPATION RHO*Ce2.epsilon/k*epsilon
 
336
!===============================================================================
 
337
 
 
338
 
 
339
! ---> Calcul de k pour la suite du sous-programme
 
340
!       on utilise un tableau de travail puisqu'il y en a...
 
341
do iel = 1, ncel
 
342
  w8(iel) = 0.5d0 * (rtpa(iel,ir11) + rtpa(iel,ir22) + rtpa(iel,ir33))
 
343
enddo
 
344
! ---> Calcul de la trace de la production, suivant qu'on est en
 
345
!     Rij standard ou en SSG (utilisation de PRODUC ou GRDVIT)
 
346
if (iturb.eq.30) then
 
347
  do iel = 1, ncel
 
348
    w9(iel) = 0.5d0*(produc(1,iel)+produc(2,iel)+produc(3,iel))
 
349
  enddo
 
350
else
 
351
  do iel = 1, ncel
 
352
    w9(iel) = -( rtpa(iel,ir11)*grdvit(iel,1,1) +               &
 
353
                 rtpa(iel,ir12)*grdvit(iel,1,2) +               &
 
354
                 rtpa(iel,ir13)*grdvit(iel,1,3) +               &
 
355
                 rtpa(iel,ir12)*grdvit(iel,2,1) +               &
 
356
                 rtpa(iel,ir22)*grdvit(iel,2,2) +               &
 
357
                 rtpa(iel,ir23)*grdvit(iel,2,3) +               &
 
358
                 rtpa(iel,ir13)*grdvit(iel,3,1) +               &
 
359
                 rtpa(iel,ir23)*grdvit(iel,3,2) +               &
 
360
                 rtpa(iel,ir33)*grdvit(iel,3,3) )
 
361
  enddo
 
362
endif
 
363
 
 
364
 
 
365
!     Terme explicite (Production)
 
366
 
 
367
do iel = 1, ncel
 
368
!     Demi-traces
 
369
  trprod = w9(iel)
 
370
  trrij  = w8(iel)
 
371
  w1(iel)   =             propce(iel,ipcroo)*volume(iel)*         &
 
372
       ce1*rtpa(iel,iep)/trrij*trprod
 
373
enddo
 
374
 
 
375
!     Si on extrapole les T.S : PROPCE
 
376
if(isto2t.gt.0) then
 
377
  do iel = 1, ncel
 
378
    propce(iel,iptsta+isou-1) =                                   &
 
379
    propce(iel,iptsta+isou-1) + w1(iel)
 
380
  enddo
 
381
!     Sinon : SMBR
 
382
else
 
383
  do iel = 1, ncel
 
384
    smbr(iel) = smbr(iel) + w1(iel)
 
385
  enddo
 
386
endif
 
387
 
 
388
!     Terme implicite (Dissipation)
 
389
do iel = 1, ncel
 
390
  trrij  = w8(iel)
 
391
  smbr(iel) = smbr(iel) - propce(iel,ipcrom)*volume(iel)*         &
 
392
                         ceps2*rtpa(iel,iep)**2/trrij
 
393
enddo
 
394
 
 
395
! ---> Matrice
 
396
 
 
397
if(isto2t.gt.0) then
 
398
  thetap = thetv
 
399
else
 
400
  thetap = 1.d0
 
401
endif
 
402
do iel = 1, ncel
 
403
  trrij  = w8(iel)
 
404
  rovsdt(iel) = rovsdt(iel) + ceps2*propce(iel,ipcrom)*volume(iel)&
 
405
                     *rtpa(iel,iep)/trrij*thetap
 
406
enddo
 
407
 
 
408
!===============================================================================
 
409
! 7. TERMES DE GRAVITE
 
410
!===============================================================================
 
411
 
 
412
if(igrari.eq.1) then
 
413
 
 
414
  ! Allocate a work array
 
415
  allocate(w7(ncelet))
 
416
 
 
417
  do iel = 1, ncel
 
418
    w7(iel) = 0.d0
 
419
  enddo
 
420
 
 
421
  call rijthe                                                     &
 
422
  !==========
 
423
 ( nvar   , nscal  ,                                              &
 
424
   ivar   , isou   , ipp    ,                                     &
 
425
   rtp    , rtpa   , propce , propfa , propfb ,                   &
 
426
   coefa  , coefb  , gradro , w7     )
 
427
 
 
428
  !     Si on extrapole les T.S. : PROPCE
 
429
  if(isto2t.gt.0) then
 
430
    do iel = 1, ncel
 
431
      propce(iel,iptsta+isou-1) = propce(iel,iptsta+isou-1) + w7(iel)
 
432
    enddo
 
433
  !     Sinon SMBR
 
434
  else
 
435
    do iel = 1, ncel
 
436
      smbr(iel) = smbr(iel) + w7(iel)
 
437
    enddo
 
438
  endif
 
439
 
 
440
  ! Free memory
 
441
  deallocate(w7)
 
442
 
 
443
endif
 
444
 
 
445
 
 
446
!===============================================================================
 
447
! 8.a TERMES DE DIFFUSION  A.grad(Eps) : PARTIE EXTRADIAGONALE EXPLICITE
 
448
!     RIJ STANDARD
 
449
!===============================================================================
 
450
 
 
451
if (iturb.eq.30) then
 
452
 
 
453
! ---> Calcul du grad(Eps)
 
454
 
 
455
  ! Allocate a temporary array for the gradient calculation
 
456
  allocate(grad(ncelet,3))
 
457
 
 
458
  iccocg = 1
 
459
  inc = 1
 
460
 
 
461
  nswrgp = nswrgr(ivar )
 
462
  imligp = imligr(ivar )
 
463
  iwarnp = iwarni(ivar )
 
464
  epsrgp = epsrgr(ivar )
 
465
  climgp = climgr(ivar )
 
466
  extrap = extrag(ivar )
 
467
 
 
468
  call grdcel                                                     &
 
469
  !==========
 
470
 ( ivar   , imrgra , inc    , iccocg , nswrgp , imligp ,          &
 
471
   iwarnp , nfecra , epsrgp , climgp , extrap ,                   &
 
472
   rtpa(1,ivar )   , coefa(1,iclvar) , coefb(1,iclvar) ,          &
 
473
   grad   )
 
474
 
 
475
! ---> Calcul des termes extradiagonaux de A.grad(Eps)
 
476
 
 
477
  do iel = 1 , ncel
 
478
    trrij  = w8(iel)
 
479
    csteps  = propce(iel,ipcroo) * crijep *trrij / rtpa(iel,iep)
 
480
    w4(iel) = csteps * (rtpa(iel,ir12)*grad(iel,2) + rtpa(iel,ir13)*grad(iel,3))
 
481
    w5(iel) = csteps * (rtpa(iel,ir12)*grad(iel,1) + rtpa(iel,ir23)*grad(iel,3))
 
482
    w6(iel) = csteps * (rtpa(iel,ir13)*grad(iel,1) + rtpa(iel,ir23)*grad(iel,2))
 
483
  enddo
 
484
 
 
485
! ---> Assemblage de { A.grad(Eps) } .S aux faces
 
486
 
 
487
  call vectds                                                     &
 
488
  !==========
 
489
( w4     , w5     , w6     ,                                      &
 
490
  viscf  , viscb  )
 
491
 
 
492
  init = 1
 
493
  call divmas(ncelet,ncel,nfac,nfabor,init,nfecra,                &
 
494
                                   ifacel,ifabor,viscf,viscb,w4)
 
495
 
 
496
!     Si on extrapole les termes sources
 
497
  if(isto2t.gt.0) then
 
498
    do iel = 1, ncel
 
499
      propce(iel,iptsta+isou-1) =                                 &
 
500
           propce(iel,iptsta+isou-1) + w4(iel)
 
501
    enddo
 
502
!     Sinon
 
503
  else
 
504
    do iel = 1, ncel
 
505
      smbr(iel) = smbr(iel) + w4(iel)
 
506
    enddo
 
507
  endif
 
508
 
 
509
 
 
510
!===============================================================================
 
511
! 8.b TERMES DE DIFFUSION  A.grad(Eps) : PARTIE DIAGONALE
 
512
!     RIJ STANDARD
 
513
!===============================================================================
 
514
!     Implicitation de (grad(eps).n)n en gradient facette
 
515
!     Si IDIFRE=1, terme correctif explicite
 
516
!        grad(eps)-(grad(eps).n)n calcule en gradient cellule
 
517
!     Les termes de bord sont uniquement pris en compte dans la partie
 
518
!        en (grad(eps).n)n
 
519
!     (W1,W2,W3) contient toujours le gradient de la variable traitee
 
520
 
 
521
!     La synchronisation des halos du gradient de epsilon a ete faite dans
 
522
!       grdcel. Pas utile de recommencer.
 
523
 
 
524
  if (idifre.eq.1) then
 
525
 
 
526
    do iel = 1, ncel
 
527
      trrij  = w8(iel)
 
528
      csteps = propce(iel,ipcroo) * crijep *trrij/rtpa(iel,iep)
 
529
      w4(iel)=csteps*rtpa(iel,ir11)
 
530
      w5(iel)=csteps*rtpa(iel,ir22)
 
531
      w6(iel)=csteps*rtpa(iel,ir33)
 
532
    enddo
 
533
 
 
534
! --->  TRAITEMENT DU PARALLELISME ET DE LA PERIODICITE
 
535
 
 
536
    if (irangp.ge.0.or.iperio.eq.1) then
 
537
      call syndia(w4, w5, w6)
 
538
      !==========
 
539
    endif
 
540
 
 
541
 
 
542
    do ifac = 1, nfac
 
543
 
 
544
      ii=ifacel(1,ifac)
 
545
      jj=ifacel(2,ifac)
 
546
 
 
547
      surfn2 = surfan(ifac)**2
 
548
 
 
549
      grdpx=0.5d0*(grad(ii,1)+grad(jj,1))
 
550
      grdpy=0.5d0*(grad(ii,2)+grad(jj,2))
 
551
      grdpz=0.5d0*(grad(ii,3)+grad(jj,3))
 
552
      grdsn=grdpx*surfac(1,ifac)+grdpy*surfac(2,ifac)             &
 
553
           +grdpz*surfac(3,ifac)
 
554
      grdpx=grdpx-grdsn*surfac(1,ifac)/surfn2
 
555
      grdpy=grdpy-grdsn*surfac(2,ifac)/surfn2
 
556
      grdpz=grdpz-grdsn*surfac(3,ifac)/surfn2
 
557
 
 
558
      viscf(ifac)= 0.5d0*(                                        &
 
559
           (w4(ii)+w4(jj))*grdpx*surfac(1,ifac)                   &
 
560
           +(w5(ii)+w5(jj))*grdpy*surfac(2,ifac)                  &
 
561
           +(w6(ii)+w6(jj))*grdpz*surfac(3,ifac))
 
562
 
 
563
    enddo
 
564
 
 
565
    ! Free memory
 
566
    deallocate(grad)
 
567
 
 
568
    do ifac = 1, nfabor
 
569
      viscb(ifac) = 0.d0
 
570
    enddo
 
571
 
 
572
    init = 1
 
573
    call divmas(ncelet,ncel,nfac,nfabor,init,nfecra,              &
 
574
       ifacel,ifabor,viscf,viscb,w1)
 
575
 
 
576
!     Si on extrapole les termes sources
 
577
    if(isto2t.gt.0) then
 
578
      do iel = 1, ncel
 
579
        propce(iel,iptsta+isou-1) =                               &
 
580
             propce(iel,iptsta+isou-1) + w1(iel)
 
581
      enddo
 
582
!     Sinon
 
583
    else
 
584
      do iel = 1, ncel
 
585
        smbr(iel) = smbr(iel) + w1(iel)
 
586
      enddo
 
587
    endif
 
588
 
 
589
  endif
 
590
 
 
591
! ---> Viscosite orthotrope pour partie implicite
 
592
 
 
593
  if( idiff(ivar).ge. 1 ) then
 
594
    do iel = 1, ncel
 
595
      trrij  = w8(iel)
 
596
      rctse = propce(iel,ipcrom) * crijep * trrij/rtpa(iel,iep)
 
597
      w1(iel) = propce(iel,ipcvis) + idifft(ivar)*rctse*rtpa(iel,ir11)
 
598
      w2(iel) = propce(iel,ipcvis) + idifft(ivar)*rctse*rtpa(iel,ir22)
 
599
      w3(iel) = propce(iel,ipcvis) + idifft(ivar)*rctse*rtpa(iel,ir33)
 
600
    enddo
 
601
 
 
602
    call visort                                                   &
 
603
    !==========
 
604
 ( imvisf ,                                                       &
 
605
   w1     , w2     , w3     ,                                     &
 
606
   viscf  , viscb  )
 
607
 
 
608
  else
 
609
 
 
610
    do ifac = 1, nfac
 
611
      viscf(ifac) = 0.d0
 
612
    enddo
 
613
    do ifac = 1, nfabor
 
614
      viscb(ifac) = 0.d0
 
615
    enddo
 
616
 
 
617
  endif
 
618
 
 
619
else
 
620
!===============================================================================
 
621
! 8.c TERMES DE DIFFUSION
 
622
!     RIJ SSG
 
623
!===============================================================================
 
624
! ---> Viscosite
 
625
 
 
626
  if( idiff(ivar).ge. 1 ) then
 
627
    do iel = 1, ncel
 
628
      w1(iel) = propce(iel,ipcvis)                                &
 
629
           + idifft(ivar)*propce(iel,ipcvst)/sigmae
 
630
    enddo
 
631
 
 
632
    call viscfa                                                   &
 
633
   !==========
 
634
 ( imvisf ,                                                       &
 
635
   w1     ,                                                       &
 
636
   viscf  , viscb  )
 
637
 
 
638
  else
 
639
 
 
640
    do ifac = 1, nfac
 
641
      viscf(ifac) = 0.d0
 
642
    enddo
 
643
    do ifac = 1, nfabor
 
644
      viscb(ifac) = 0.d0
 
645
    enddo
 
646
 
 
647
  endif
 
648
 
 
649
endif
 
650
 
 
651
 
 
652
!===============================================================================
 
653
! 9. RESOLUTION
 
654
!===============================================================================
 
655
 
 
656
if(isto2t.gt.0) then
 
657
  thetp1 = 1.d0 + thets
 
658
  do iel = 1, ncel
 
659
    smbr(iel) = smbr(iel) + thetp1*propce(iel,iptsta+isou-1)
 
660
  enddo
 
661
endif
 
662
 
 
663
iconvp = iconv (ivar)
 
664
idiffp = idiff (ivar)
 
665
ndircp = ndircl(ivar)
 
666
ireslp = iresol(ivar)
 
667
nitmap = nitmax(ivar)
 
668
nswrsp = nswrsm(ivar)
 
669
nswrgp = nswrgr(ivar)
 
670
imligp = imligr(ivar)
 
671
ircflp = ircflu(ivar)
 
672
ischcp = ischcv(ivar)
 
673
isstpp = isstpc(ivar)
 
674
iescap = 0
 
675
imgrp  = imgr  (ivar)
 
676
ncymxp = ncymax(ivar)
 
677
nitmfp = nitmgf(ivar)
 
678
!MO      IPP    =
 
679
iwarnp = iwarni(ivar)
 
680
blencp = blencv(ivar)
 
681
epsilp = epsilo(ivar)
 
682
epsrsp = epsrsm(ivar)
 
683
epsrgp = epsrgr(ivar)
 
684
climgp = climgr(ivar)
 
685
extrap = extrag(ivar)
 
686
relaxp = relaxv(ivar)
 
687
 
 
688
call codits                                                       &
 
689
!==========
 
690
 ( nvar   , nscal  ,                                              &
 
691
   idtvar , ivar   , iconvp , idiffp , ireslp , ndircp , nitmap , &
 
692
   imrgra , nswrsp , nswrgp , imligp , ircflp ,                   &
 
693
   ischcp , isstpp , iescap ,                                     &
 
694
   imgrp  , ncymxp , nitmfp , ipp    , iwarnp ,                   &
 
695
   blencp , epsilp , epsrsp , epsrgp , climgp , extrap ,          &
 
696
   relaxp , thetv  ,                                              &
 
697
   rtpa(1,ivar)    , rtpa(1,ivar)    ,                            &
 
698
                     coefa(1,iclvar) , coefb(1,iclvar) ,          &
 
699
                     coefa(1,iclvaf) , coefb(1,iclvaf) ,          &
 
700
                     propfa(1,iflmas), propfb(1,iflmab),          &
 
701
   viscf  , viscb  , viscf  , viscb  ,                            &
 
702
   rovsdt , smbr   , rtp(1,ivar)     ,                            &
 
703
   rvoid  )
 
704
 
 
705
!===============================================================================
 
706
! 10. IMPRESSIONS
 
707
!===============================================================================
 
708
 
 
709
! Free memory
 
710
deallocate(w1, w2, w3)
 
711
deallocate(w4, w5, w6)
 
712
deallocate(w8, w9)
 
713
 
 
714
!--------
 
715
! FORMATS
 
716
!--------
 
717
 
 
718
#if defined(_CS_LANG_FR)
 
719
 
 
720
 1000 format(/,'           RESOLUTION POUR LA VARIABLE ',A8,/)
 
721
 
 
722
#else
 
723
 
 
724
 1000 format(/,'           SOLVING VARIABLE ',A8           ,/)
 
725
 
 
726
#endif
 
727
 
 
728
!12345678 : MAX: 12345678901234 MIN: 12345678901234 NORM: 12345678901234
 
729
!----
 
730
! FIN
 
731
!----
 
732
 
 
733
return
 
734
 
 
735
end subroutine