← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/base/navstv.f90

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2011-11-24 00:00:08 UTC
  • mfrom: (6.1.9 sid)
  • Revision ID: package-import@ubuntu.com-20111124000008-2vo99e38267942q5
Tags: 2.1.0-3
Install a missing file

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/base/navstv.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 navstv &
 
24
!================
 
25
 
 
26
 ( nvar   , nscal  , iterns , icvrge ,                            &
 
27
   isostd ,                                                       &
 
28
   dt     , tpucou , rtp    , rtpa   , propce , propfa , propfb , &
 
29
   tslagr , coefa  , coefb  , frcxt  ,                            &
 
30
   trava  , ximpa  , uvwk   )
 
31
 
 
32
!===============================================================================
 
33
! FONCTION :
 
34
! ----------
 
35
 
 
36
! Solving of NS equations for incompressible or slightly compressible flows for
 
37
! one time step. Both convection-diffusion and continuity steps are performed.
 
38
! The velocity components are solved together in once.
 
39
 
 
40
!-------------------------------------------------------------------------------
 
41
!ARGU                             ARGUMENTS
 
42
!__________________.____._____.________________________________________________.
 
43
! name             !type!mode ! role                                           !
 
44
!__________________!____!_____!________________________________________________!
 
45
! nvar             ! i  ! <-- ! total number of variables                      !
 
46
! nscal            ! i  ! <-- ! total number of scalars                        !
 
47
! iterns           ! e  ! <-- ! numero d'iteration sur navsto                  !
 
48
! icvrge           ! e  ! <-- ! indicateur de convergence du pnt fix           !
 
49
! isostd           ! te ! <-- ! indicateur de sortie standard                  !
 
50
!    (nfabor+1)    !    !     !  +numero de la face de reference               !
 
51
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
 
52
! tpucou(ncelet,3) ! ra ! <-- ! velocity-pressure coupling (interleaved)       !
 
53
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
 
54
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
 
55
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
 
56
! propfa(nfac, *)  ! ra ! <-- ! physical properties at interior face centers   !
 
57
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers   !
 
58
! coefa, coefb     ! ra ! <-- ! boundary conditions                            !
 
59
!  (nfabor, *)     !    !     !                                                !
 
60
! frcxt(ncelet,3)  ! tr ! <-- ! force exterieure generant la pression          !
 
61
!                  !    !     !  hydrostatique                                 !
 
62
! tslagr           ! tr ! <-- ! terme de couplage retour du                    !
 
63
!(ncelet,*)        !    !     !     lagrangien                                 !
 
64
! trava            ! tr ! <-- ! tableau de travail pour couplage               !
 
65
!(3,ncelet)        !    !     ! vitesse pression par point fixe                !
 
66
! ximpa            ! tr ! <-- ! tableau de travail pour couplage               !
 
67
!(3,3,ncelet)      !    !     ! vitesse pression par point fixe                !
 
68
! uvwk             ! tr ! <-- ! tableau de travail pour couplage u/p           !
 
69
!(3,ncelet)        !    !     ! sert a stocker la vitesse de                   !
 
70
!                  !    !     ! l'iteration precedente                         !
 
71
!__________________!____!_____!________________________________________________!
 
72
 
 
73
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
 
74
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
 
75
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
 
76
!            --- tableau de travail
 
77
 
 
78
!===============================================================================
 
79
 
 
80
!===============================================================================
 
81
! Module files
 
82
!===============================================================================
 
83
 
 
84
use paramx
 
85
use dimens, only: ndimfb
 
86
use numvar
 
87
use entsor
 
88
use cstphy
 
89
use cstnum
 
90
use optcal
 
91
use pointe
 
92
use albase
 
93
use parall
 
94
use period
 
95
use ppppar
 
96
use ppthch
 
97
use ppincl
 
98
use cplsat
 
99
use mesh
 
100
 
 
101
!===============================================================================
 
102
 
 
103
implicit none
 
104
 
 
105
! Arguments
 
106
 
 
107
integer          nvar   , nscal  , iterns , icvrge
 
108
 
 
109
integer          isostd(nfabor+1)
 
110
 
 
111
double precision dt(ncelet), tpucou(ncelet,3), rtp(ncelet,*), rtpa(ncelet,*)
 
112
double precision propce(ncelet,*)
 
113
double precision propfa(nfac,*), propfb(ndimfb,*)
 
114
double precision tslagr(ncelet,*)
 
115
double precision coefa(ndimfb,*), coefb(ndimfb,*)
 
116
double precision frcxt(ncelet,3)
 
117
double precision trava(ndim,ncelet),ximpa(ndim,ndim,ncelet)
 
118
double precision uvwk(ndim,ncelet)
 
119
 
 
120
! Local variables
 
121
 
 
122
integer          iccocg, inc, iel, iel1, iel2, ifac, imax
 
123
integer          ii    , inod
 
124
integer          isou, ivar, iitsm
 
125
integer          iclipr, iclipf
 
126
integer          icliup, iclivp, icliwp, init
 
127
integer          icluma, iclvma, iclwma
 
128
integer          iflmas, iflmab, ipcrom, ipbrom
 
129
integer          iflms1, iflmb1, iflmb0
 
130
integer          nswrgp, imligp, iwarnp
 
131
integer          nbrval, iappel, iescop, idtsca
 
132
integer          ndircp, icpt  , iecrw
 
133
integer          numcpl
 
134
double precision rnorm , rnorma, rnormi, vitnor
 
135
double precision dtsrom, unsrom, surf  , rhom, rovolsdt
 
136
double precision epsrgp, climgp, extrap, xyzmax(3)
 
137
double precision thetap, xdu, xdv, xdw
 
138
double precision xxp0 , xyp0 , xzp0
 
139
double precision rhofac, dtfac, ddepx , ddepy, ddepz
 
140
double precision xnrdis
 
141
double precision vitbox, vitboy, vitboz
 
142
 
 
143
double precision rvoid(1)
 
144
 
 
145
double precision, allocatable, dimension(:), target :: viscf, viscb
 
146
double precision, allocatable, dimension(:), target :: wvisfi, wvisbi
 
147
double precision, allocatable, dimension(:) :: drtp, smbr, rovsdt
 
148
double precision, allocatable, dimension(:,:) :: trav
 
149
double precision, allocatable, dimension(:) :: w1, w2, w3
 
150
double precision, allocatable, dimension(:) :: w4, w5, w6
 
151
double precision, allocatable, dimension(:) :: w7, w8, w9
 
152
double precision, allocatable, dimension(:) :: w10
 
153
double precision, allocatable, dimension(:,:) :: dfrcxt
 
154
double precision, allocatable, dimension(:,:) :: frchy, dfrchy
 
155
double precision, allocatable, dimension(:) :: esflum, esflub
 
156
 
 
157
double precision, pointer, dimension(:) :: viscfi => null(), viscbi => null()
 
158
 
 
159
double precision, dimension(:,:), allocatable :: gradp
 
160
double precision, dimension(:,:), allocatable :: vel
 
161
double precision, dimension(:,:), allocatable :: vela
 
162
double precision, dimension(:,:), allocatable :: tpucov
 
163
 
 
164
!===============================================================================
 
165
 
 
166
!===============================================================================
 
167
! 1.  INITIALISATION
 
168
!===============================================================================
 
169
 
 
170
! Allocate temporary arrays for the velocity-pressure resolution
 
171
allocate(viscf(nfac), viscb(nfabor))
 
172
allocate(drtp(ncelet), smbr(ncelet), rovsdt(ncelet))
 
173
allocate(trav(3,ncelet))
 
174
allocate(vela(3,ncelet))
 
175
allocate(vel(3,ncelet))
 
176
if (ipucou.eq.1 .or. ncpdct.gt.0) then
 
177
  allocate(tpucov(3,ncelet))
 
178
endif
 
179
 
 
180
! Allocate other arrays, depending on user options
 
181
!if (iphydr.eq.1) allocate(dfrcxt(ncelet,3))
 
182
allocate(dfrcxt(ncelet,3))
 
183
if (icalhy.eq.1) allocate(frchy(ncelet,ndim), dfrchy(ncelet,ndim))
 
184
if (iescal(iestot).gt.0) allocate(esflum(nfac), esflub(nfabor))
 
185
if (itytur.eq.3.and.irijnu.eq.1) then
 
186
  allocate(wvisfi(nfac), wvisbi(nfabor))
 
187
  viscfi => wvisfi(1:nfac)
 
188
  viscbi => wvisbi(1:nfabor)
 
189
else
 
190
  viscfi => viscf(1:nfac)
 
191
  viscbi => viscb(1:nfabor)
 
192
endif
 
193
 
 
194
! Allocate work arrays
 
195
allocate(w1(ncelet), w2(ncelet), w3(ncelet))
 
196
allocate(w4(ncelet), w5(ncelet), w6(ncelet))
 
197
allocate(w7(ncelet), w8(ncelet), w9(ncelet))
 
198
if (irnpnw.eq.1) allocate(w10(ncelet))
 
199
 
 
200
! Interleaved value of vel and vela and tpucou
 
201
do iel = 1, ncelet
 
202
  vel (1,iel) = rtp (iel,iu)
 
203
  vel (2,iel) = rtp (iel,iv)
 
204
  vel (3,iel) = rtp (iel,iw)
 
205
  vela(1,iel) = rtpa(iel,iu)
 
206
  vela(2,iel) = rtpa(iel,iv)
 
207
  vela(3,iel) = rtpa(iel,iw)
 
208
enddo
 
209
 
 
210
if(iwarni(iu).ge.1) then
 
211
  write(nfecra,1000)
 
212
endif
 
213
 
 
214
! Initialize variables to avoid compiler warnings
 
215
 
 
216
ivar = 0
 
217
iflmas = 0
 
218
ipcrom = 0
 
219
imax = 0
 
220
 
 
221
! Memoire
 
222
 
 
223
if(nterup.gt.1) then
 
224
 
 
225
  do iel = 1,ncelet
 
226
    do isou = 1, 3
 
227
    !     La boucle sur NCELET est une securite au cas
 
228
    !       ou on utiliserait UVWK par erreur a ITERNS = 1
 
229
      uvwk(isou,iel) = vel(isou,iel)
 
230
    enddo
 
231
  enddo
 
232
 
 
233
  ! Calcul de la norme L2 de la vitesse
 
234
  if(iterns.eq.1) then
 
235
    xnrmu0 = 0.d0
 
236
    do iel = 1, ncel
 
237
      xnrmu0 = xnrmu0 +(vela(1,iel)**2        &
 
238
                      + vela(2,iel)**2        &
 
239
                      + vela(3,iel)**2)       &
 
240
                      * volume(iel)
 
241
    enddo
 
242
    if(irangp.ge.0) then
 
243
      call parsom (xnrmu0)
 
244
      !==========
 
245
    endif
 
246
    ! En cas de couplage entre deux instances de Code_Saturne, on calcule
 
247
    ! la norme totale de la vitesse
 
248
    ! Necessaire pour que l'une des instances ne stoppe pas plus tot que les autres
 
249
    ! (il faudrait quand meme verifier les options numeriques, ...)
 
250
    do numcpl = 1, nbrcpl
 
251
      call tbrcpl ( numcpl, 1, 1, xnrmu0, xnrdis )
 
252
      !==========
 
253
      xnrmu0 = xnrmu0 + xnrdis
 
254
    enddo
 
255
    xnrmu0 = sqrt(xnrmu0)
 
256
  endif
 
257
 
 
258
  ! On assure la periodicite ou le parallelisme de UVWK et la pression
 
259
  ! (cette derniere vaut la pression a l'iteration precedente)
 
260
  if(iterns.gt.1) then
 
261
    if (irangp.ge.0.or.iperio.eq.1) then
 
262
      call synvin(uvwk(1,1))
 
263
      !==========
 
264
      call synsca(rtpa(1,ipr))
 
265
      !==========
 
266
    endif
 
267
  endif
 
268
 
 
269
endif
 
270
 
 
271
 
 
272
!===============================================================================
 
273
! 2.  ETAPE DE PREDICTION DES VITESSES
 
274
!===============================================================================
 
275
 
 
276
iappel = 1
 
277
iflmas = ipprof(ifluma(iu))
 
278
iflmab = ipprob(ifluma(iu))
 
279
 
 
280
call predvv &
 
281
!==========
 
282
( iappel ,                                                       &
 
283
  nvar   , nscal  , iterns ,                                     &
 
284
  ncepdc , ncetsm ,                                              &
 
285
  icepdc , icetsm , itypsm ,                                     &
 
286
  dt     , rtp    , rtpa   , vel    , vela   ,                   &
 
287
  propce , propfa , propfb ,                                     &
 
288
  propfa(1,iflmas), propfb(1,iflmab),                            &
 
289
  tslagr , coefa  , coefb  , coefau , coefbu , cofafu , cofbfu , &
 
290
  ckupdc , smacel , frcxt ,                                      &
 
291
  trava  , ximpa  , uvwk   , dfrcxt , tpucov ,  trav  ,          &
 
292
  viscf  , viscb  , viscfi , viscbi ,                            &
 
293
  w1     , w7     , w8     , w9     , w10    )
 
294
 
 
295
! --- Sortie si pas de pression continuite
 
296
!       on met a jour les flux de masse, et on sort
 
297
 
 
298
if( iprco.le.0 ) then
 
299
 
 
300
  icliup = iclrtp(iu,icoef)
 
301
  iclivp = iclrtp(iv,icoef)
 
302
  icliwp = iclrtp(iw,icoef)
 
303
 
 
304
  iflmas = ipprof(ifluma(iu))
 
305
  iflmab = ipprob(ifluma(iu))
 
306
  ipcrom = ipproc(irom)
 
307
  ipbrom = ipprob(irom)
 
308
 
 
309
  init   = 1
 
310
  inc    = 1
 
311
  iccocg = 1
 
312
  iflmb0 = 1
 
313
  if (iale.eq.1) iflmb0 = 0
 
314
  nswrgp = nswrgr(iu)
 
315
  imligp = imligr(iu)
 
316
  iwarnp = iwarni(iu)
 
317
  epsrgp = epsrgr(iu)
 
318
  climgp = climgr(iu)
 
319
  extrap = extrag(iu)
 
320
 
 
321
  call inimav                                                     &
 
322
  !==========
 
323
 ( nvar   , nscal  ,                                              &
 
324
   iu     ,                                                       &
 
325
   iflmb0 , init   , inc    , imrgra ,iccocg , nswrgp , imligp ,  &
 
326
   iwarnp , nfecra ,                                              &
 
327
   epsrgp , climgp , extrap ,                                     &
 
328
   propce(1,ipcrom), propfb(1,ipbrom),                            &
 
329
   vel    ,                                                       &
 
330
   coefau , coefbu ,                                              &
 
331
   propfa(1,iflmas), propfb(1,iflmab)  )
 
332
 
 
333
 
 
334
  ! Ajout de la vitesse du solide dans le flux convectif,
 
335
  ! si le maillage est mobile (solide rigide)
 
336
  ! En turbomachine, on conna�t exactement la vitesse de maillage � ajouter
 
337
  if (imobil.eq.1) then
 
338
 
 
339
    iflmas = ipprof(ifluma(iu))
 
340
    iflmab = ipprob(ifluma(iu))
 
341
    ipcrom = ipproc(irom)
 
342
    ipbrom = ipprob(irom)
 
343
 
 
344
    do ifac = 1, nfac
 
345
      iel1 = ifacel(1,ifac)
 
346
      iel2 = ifacel(2,ifac)
 
347
      dtfac  = 0.5d0*(dt(iel1) + dt(iel2))
 
348
      rhofac = 0.5d0*(propce(iel1,ipcrom) + propce(iel2,ipcrom))
 
349
      vitbox = omegay*cdgfac(3,ifac) - omegaz*cdgfac(2,ifac)
 
350
      vitboy = omegaz*cdgfac(1,ifac) - omegax*cdgfac(3,ifac)
 
351
      vitboz = omegax*cdgfac(2,ifac) - omegay*cdgfac(1,ifac)
 
352
      propfa(ifac,iflmas) = propfa(ifac,iflmas) - rhofac*(        &
 
353
           vitbox*surfac(1,ifac) + vitboy*surfac(2,ifac) + vitboz*surfac(3,ifac) )
 
354
    enddo
 
355
    do ifac = 1, nfabor
 
356
      iel = ifabor(ifac)
 
357
      dtfac  = dt(iel)
 
358
      rhofac = propfb(ifac,ipbrom)
 
359
      vitbox = omegay*cdgfbo(3,ifac) - omegaz*cdgfbo(2,ifac)
 
360
      vitboy = omegaz*cdgfbo(1,ifac) - omegax*cdgfbo(3,ifac)
 
361
      vitboz = omegax*cdgfbo(2,ifac) - omegay*cdgfbo(1,ifac)
 
362
      propfb(ifac,iflmab) = propfb(ifac,iflmab) - rhofac*(        &
 
363
           vitbox*surfbo(1,ifac) + vitboy*surfbo(2,ifac) + vitboz*surfbo(3,ifac) )
 
364
    enddo
 
365
 
 
366
  endif
 
367
 
 
368
  return
 
369
 
 
370
endif
 
371
 
 
372
!===============================================================================
 
373
! 3.  ETAPE DE PRESSION/CONTINUITE ( VITESSE/PRESSION )
 
374
!===============================================================================
 
375
 
 
376
if(iwarni(iu).ge.1) then
 
377
  write(nfecra,1200)
 
378
endif
 
379
 
 
380
! --- Pas de temps scalaire ou pas
 
381
idtsca = 0
 
382
if ((ipucou.eq.1).or.(ncpdct.gt.0)) idtsca = 1
 
383
 
 
384
call resopv &
 
385
!==========
 
386
 ( nvar   , nscal  ,                                              &
 
387
   ncepdc , ncetsm ,                                              &
 
388
   icepdc , icetsm , itypsm ,                                     &
 
389
   isostd , idtsca ,                                              &
 
390
   dt     , rtp    , rtpa   , vel    , vela   ,                   &
 
391
   propce , propfa , propfb ,                                     &
 
392
   coefa  , coefb  , coefau , coefbu ,                            &
 
393
   ckupdc , smacel ,                                              &
 
394
   frcxt  , dfrcxt , tpucov , trav   ,                            &
 
395
   viscf  , viscb  , viscfi , viscbi ,                            &
 
396
   drtp   , smbr   , rovsdt , tslagr ,                            &
 
397
   frchy  , dfrchy , trava  )
 
398
 
 
399
 
 
400
!===============================================================================
 
401
! 4.  REACTUALISATION DU CHAMP DE VITESSE
 
402
!===============================================================================
 
403
 
 
404
iclipr = iclrtp(ipr,icoef)
 
405
iclipf = iclrtp(ipr,icoeff)
 
406
icliup = iclrtp(iu ,icoef)
 
407
iclivp = iclrtp(iv ,icoef)
 
408
icliwp = iclrtp(iw ,icoef)
 
409
 
 
410
iflmas = ipprof(ifluma(iu))
 
411
iflmab = ipprob(ifluma(iu))
 
412
ipcrom = ipproc(irom  )
 
413
ipbrom = ipprob(irom  )
 
414
 
 
415
 
 
416
!       IREVMC = 0 : Only the standard method is available for the coupled
 
417
!                    version of navstv.
 
418
 
 
419
!     The predicted velocity is corrected by the cell gradient of the
 
420
!     pressure increment.
 
421
 
 
422
!     GRADIENT DE L'INCREMENT TOTAL DE PRESSION
 
423
 
 
424
if (idtvar.lt.0) then
 
425
  do iel = 1, ncel
 
426
    drtp(iel) = (rtp(iel,ipr) -rtpa(iel,ipr)) / relaxv(ipr)
 
427
  enddo
 
428
else
 
429
  do iel = 1, ncel
 
430
    drtp(iel) = rtp(iel,ipr) -rtpa(iel,ipr)
 
431
  enddo
 
432
endif
 
433
 
 
434
! --->    TRAITEMENT DU PARALLELISME ET DE LA PERIODICITE
 
435
 
 
436
if (irangp.ge.0.or.iperio.eq.1) then
 
437
  call synsca(drtp)
 
438
  !==========
 
439
endif
 
440
 
 
441
 
 
442
iccocg = 1
 
443
inc = 0
 
444
if (iphydr.eq.1) inc = 1
 
445
nswrgp = nswrgr(ipr)
 
446
imligp = imligr(ipr)
 
447
iwarnp = iwarni(ipr)
 
448
epsrgp = epsrgr(ipr)
 
449
climgp = climgr(ipr)
 
450
extrap = extrag(ipr)
 
451
 
 
452
!Allocation
 
453
allocate(gradp (ncelet,3))
 
454
 
 
455
call grdpot &
 
456
!==========
 
457
( ipr , imrgra , inc    , iccocg , nswrgp , imligp , iphydr ,    &
 
458
 iwarnp , nfecra , epsrgp , climgp , extrap ,                   &
 
459
 rvoid  ,                                                       &
 
460
 dfrcxt(1,1),dfrcxt(1,2),dfrcxt(1,3),                           &
 
461
 drtp   , coefa(1,iclipf) , coefb(1,iclipr)  ,                  &
 
462
 gradp  )
 
463
 
 
464
!     REACTUALISATION DU CHAMP DE VITESSES
 
465
 
 
466
thetap = thetav(ipr)
 
467
do iel = 1, ncelet
 
468
  do isou = 1, 3
 
469
    trav(isou,iel) = gradp(iel,isou)
 
470
  enddo
 
471
enddo
 
472
 
 
473
!Free memory
 
474
deallocate(gradp)
 
475
 
 
476
 
 
477
!     REACTUALISATION DU CHAMP DE VITESSES
 
478
 
 
479
thetap = thetav(ipr)
 
480
if (iphydr.eq.0) then
 
481
  if (idtsca.eq.0) then
 
482
    do iel = 1, ncel
 
483
      dtsrom = -thetap*dt(iel)/propce(iel,ipcrom)
 
484
      do isou = 1, 3
 
485
        vel(isou,iel) = vel(isou,iel)+dtsrom*trav(isou,iel)
 
486
      enddo
 
487
    enddo
 
488
  else
 
489
    do iel = 1, ncel
 
490
      unsrom = -thetap/propce(iel,ipcrom)
 
491
      ! tpucov is an interleaved array
 
492
      do isou = 1, 3
 
493
        vel(isou,iel) = vel(isou,iel) + unsrom*tpucov(isou,iel)*trav(isou,iel)
 
494
      enddo
 
495
    enddo
 
496
  endif
 
497
else
 
498
  if (idtsca.eq.0) then
 
499
    do iel = 1, ncel
 
500
      dtsrom = thetap*dt(iel)/propce(iel,ipcrom)
 
501
      do isou = 1, 3
 
502
        vel(isou,iel) = vel(isou,iel)                           &
 
503
           +dtsrom*(dfrcxt(iel,isou)-trav(isou,iel) )
 
504
      enddo
 
505
    enddo
 
506
  else
 
507
    do iel = 1, ncel
 
508
      unsrom = thetap/propce(iel,ipcrom)
 
509
      ! tpucov is an interleaved array
 
510
      do isou = 1, 3
 
511
        vel(isou,iel) = vel(isou,iel)                         &
 
512
             +unsrom*tpucov(isou,iel)                         &
 
513
             *(dfrcxt(iel,isou)-trav(isou,iel))
 
514
      enddo
 
515
    enddo
 
516
  endif
 
517
!     mise a jour des forces exterieures pour le calcul des gradients
 
518
  do iel=1,ncel
 
519
    frcxt(iel,1) = frcxt(iel,1) + dfrcxt(iel,1)
 
520
    frcxt(iel,2) = frcxt(iel,2) + dfrcxt(iel,2)
 
521
    frcxt(iel,3) = frcxt(iel,3) + dfrcxt(iel,3)
 
522
  enddo
 
523
  if (irangp.ge.0.or.iperio.eq.1) then
 
524
    call synvec(frcxt(1,1), frcxt(1,2), frcxt(1,3))
 
525
    !==========
 
526
  endif
 
527
  !     mise a jour des Dirichlets de pression en sortie dans COEFA
 
528
  iclipr = iclrtp(ipr,icoef)
 
529
  iclipf = iclrtp(ipr,icoeff)
 
530
  do ifac = 1,nfabor
 
531
    if (isostd(ifac).eq.1)                              &
 
532
         coefa(ifac,iclipr) = coefa(ifac,iclipr)        &
 
533
         + coefa(ifac,iclipf)
 
534
  enddo
 
535
endif
 
536
 
 
537
 
 
538
!===============================================================================
 
539
! 5.  CALCUL D'UN ESTIMATEUR D'ERREUR DE L'ETAPE DE CORRECTION ET TOTAL
 
540
!===============================================================================
 
541
 
 
542
if(iescal(iescor).gt.0.or.iescal(iestot).gt.0) then
 
543
 
 
544
  ! ---> REPERAGE DES VARIABLES
 
545
 
 
546
  icliup = iclrtp(iu,icoef)
 
547
  iclivp = iclrtp(iv,icoef)
 
548
  icliwp = iclrtp(iw,icoef)
 
549
 
 
550
  ipcrom = ipproc(irom)
 
551
  ipbrom = ipprob(irom)
 
552
 
 
553
 
 
554
  ! ---> ECHANGE DES VITESSES ET PRESSION EN PERIODICITE ET PARALLELISME
 
555
 
 
556
  !    Pour les estimateurs IESCOR et IESTOT, la vitesse doit etre echangee.
 
557
 
 
558
  !    Pour l'estimateur IESTOT, la pression doit etre echangee aussi.
 
559
 
 
560
  !    Cela ne remplace pas l'echange du debut de pas de temps
 
561
  !     a cause de usproj qui vient plus tard et des calculs suite)
 
562
 
 
563
 
 
564
  ! --- Vitesse
 
565
 
 
566
  if (irangp.ge.0.or.iperio.eq.1) then
 
567
    call synvin(vel)
 
568
    !==========
 
569
  endif
 
570
 
 
571
  !  -- Pression
 
572
 
 
573
  if(iescal(iestot).gt.0) then
 
574
 
 
575
    if (irangp.ge.0.or.iperio.eq.1) then
 
576
      call synsca(rtp(1,ipr))
 
577
      !==========
 
578
    endif
 
579
 
 
580
  endif
 
581
 
 
582
 
 
583
  ! ---> CALCUL DU FLUX DE MASSE DEDUIT DE LA VITESSE REACTUALISEE
 
584
 
 
585
  init   = 1
 
586
  inc    = 1
 
587
  iccocg = 1
 
588
  iflmb0 = 1
 
589
  if (iale.eq.1) iflmb0 = 0
 
590
  nswrgp = nswrgr(iu)
 
591
  imligp = imligr(iu)
 
592
  iwarnp = iwarni(iu)
 
593
  epsrgp = epsrgr(iu)
 
594
  climgp = climgr(iu)
 
595
  extrap = extrag(iu)
 
596
 
 
597
  call inimav                                                     &
 
598
  !==========
 
599
 ( nvar   , nscal  ,                                              &
 
600
   iu     ,                                                       &
 
601
   iflmb0 , init   , inc    , imrgra , iccocg , nswrgp , imligp , &
 
602
   iwarnp , nfecra ,                                              &
 
603
   epsrgp , climgp , extrap ,                                     &
 
604
   propce(1,ipcrom), propfb(1,ipbrom),                            &
 
605
   vel    ,                                                       &
 
606
   coefau , coefbu ,                                              &
 
607
   esflum , esflub )
 
608
 
 
609
 
 
610
  ! ---> CALCUL DE L'ESTIMATEUR CORRECTION : DIVERGENCE DE ROM * U (N + 1)
 
611
  !                                          - GAMMA
 
612
 
 
613
  if(iescal(iescor).gt.0) then
 
614
    init = 1
 
615
    call divmas(ncelet,ncel,nfac,nfabor,init,nfecra,            &
 
616
         ifacel,ifabor,esflum,esflub,w1)
 
617
 
 
618
    if (ncetsm.gt.0) then
 
619
      do iitsm = 1, ncetsm
 
620
        iel = icetsm(iitsm)
 
621
        w1(iel) = w1(iel)-volume(iel)*smacel(iitsm,ipr)
 
622
      enddo
 
623
    endif
 
624
 
 
625
    if(iescal(iescor).eq.2) then
 
626
      iescop = ipproc(iestim(iescor))
 
627
      do iel = 1, ncel
 
628
        propce(iel,iescop) =  abs(w1(iel))
 
629
      enddo
 
630
    elseif(iescal(iescor).eq.1) then
 
631
      iescop = ipproc(iestim(iescor))
 
632
      do iel = 1, ncel
 
633
        propce(iel,iescop) =  abs(w1(iel)) / volume(iel)
 
634
      enddo
 
635
    endif
 
636
  endif
 
637
 
 
638
 
 
639
  ! ---> CALCUL DE L'ESTIMATEUR TOTAL
 
640
 
 
641
  if(iescal(iestot).gt.0) then
 
642
 
 
643
    !   INITIALISATION DE TRAV AVEC LE TERME INSTATIONNAIRE
 
644
 
 
645
    do iel = 1, ncel
 
646
      rovolsdt = propce(iel,ipcrom)*volume(iel)/dt(iel)
 
647
      do isou = 1, 3
 
648
        trav(isou,iel) = rovolsdt *                               &
 
649
                 ( vela(isou,iel)- vel(isou,iel) )
 
650
      enddo
 
651
    enddo
 
652
 
 
653
    !   APPEL A PREDUV AVEC RTP ET RTP AU LIEU DE RTP ET RTPA
 
654
    !                  AVEC LE FLUX DE MASSE RECALCULE
 
655
    iappel = 2
 
656
    call predvv &
 
657
    !==========
 
658
 ( iappel ,                                                       &
 
659
   nvar   , nscal  , iterns ,                                     &
 
660
   ncepdc , ncetsm ,                                              &
 
661
   icepdc , icetsm , itypsm ,                                     &
 
662
   dt     , rtp    , rtp    , vel    , vel    ,                   &
 
663
   propce , propfa , propfb ,                                     &
 
664
   esflum , esflub ,                                              &
 
665
   tslagr , coefa  , coefb  , coefau , coefbu , cofafu , cofbfu , &
 
666
   ckupdc , smacel , frcxt  ,                                     &
 
667
   trava  , ximpa  , uvwk   , dfrcxt , tpucov , trav   ,          &
 
668
   viscf  , viscb  , viscfi , viscbi ,                            &
 
669
   w1     , w7     , w8     , w9     , w10    )
 
670
 
 
671
  endif
 
672
 
 
673
endif
 
674
 
 
675
!===============================================================================
 
676
! 6.  TRAITEMENT DU POINT FIXE SUR LE SYSTEME VITESSE/PRESSION
 
677
!===============================================================================
 
678
 
 
679
if(nterup.gt.1) then
 
680
! TEST DE CONVERGENCE DE L'ALGORITHME ITERATIF
 
681
! On initialise ICVRGE a 1 et on le met a 0 si une des phases n'est
 
682
! pas convergee
 
683
 
 
684
  icvrge = 1
 
685
 
 
686
  do iel = 1,ncel
 
687
    xdu = vel(1,iel) - uvwk(1,iel)
 
688
    xdv = vel(2,iel) - uvwk(2,iel)
 
689
    xdw = vel(3,iel) - uvwk(3,iel)
 
690
    xnrmu = xnrmu +(xdu**2 + xdv**2 + xdw**2)     &
 
691
                                * volume(iel)
 
692
  enddo
 
693
  ! --->    TRAITEMENT DU PARALLELISME
 
694
 
 
695
  if(irangp.ge.0) call parsom (xnrmu)
 
696
                              !==========
 
697
  ! -- >    TRAITEMENT DU COUPLAGE ENTRE DEUX INSTANCES DE CODE_SATURNE
 
698
  do numcpl = 1, nbrcpl
 
699
    call tbrcpl ( numcpl, 1, 1, xnrmu, xnrdis )
 
700
    !==========
 
701
    xnrmu = xnrmu + xnrdis
 
702
  enddo
 
703
  xnrmu = sqrt(xnrmu)
 
704
 
 
705
  ! Indicateur de convergence du point fixe
 
706
  if(xnrmu.ge.epsup*xnrmu0) icvrge = 0
 
707
 
 
708
endif
 
709
 
 
710
! ---> RECALAGE DE LA PRESSION SUR UNE PRESSION A MOYENNE NULLE
 
711
!  On recale si on n'a pas de Dirichlet. Or le nombre de Dirichlets
 
712
!  calcule dans typecl.F est NDIRCL si IDIRCL=1 et NDIRCL-1 si IDIRCL=0
 
713
!  (ISTAT vaut toujours 0 pour la pression)
 
714
 
 
715
if (idircl(ipr).eq.1) then
 
716
  ndircp = ndircl(ipr)
 
717
else
 
718
  ndircp = ndircl(ipr)-1
 
719
endif
 
720
if(ndircp.le.0) then
 
721
  call prmoy0 &
 
722
  !==========
 
723
( ncelet , ncel   , nfac   , nfabor ,                         &
 
724
  volume , rtp(1,ipr) )
 
725
endif
 
726
 
 
727
! Calcul de la pression totale IPRTOT : (definie comme propriete )
 
728
! En compressible, la pression resolue est deja la pression totale
 
729
 
 
730
if (ippmod(icompf).lt.0) then
 
731
  xxp0   = xyzp0(1)
 
732
  xyp0   = xyzp0(2)
 
733
  xzp0   = xyzp0(3)
 
734
  do iel=1,ncel
 
735
    propce(iel,ipproc(iprtot))= rtp(iel,ipr)           &
 
736
         + ro0*( gx*(xyzcen(1,iel)-xxp0)               &
 
737
         + gy*(xyzcen(2,iel)-xyp0)                     &
 
738
         + gz*(xyzcen(3,iel)-xzp0) )                   &
 
739
         + p0 - pred0
 
740
  enddo
 
741
endif
 
742
 
 
743
!===============================================================================
 
744
! 7.  IMPRESSIONS
 
745
!===============================================================================
 
746
 
 
747
iflmas = ipprof(ifluma(iu))
 
748
iflmab = ipprob(ifluma(iu))
 
749
ipcrom = ipproc(irom)
 
750
ipbrom = ipprob(irom)
 
751
 
 
752
if (iwarni(iu).ge.1) then
 
753
 
 
754
  write(nfecra,2000)
 
755
 
 
756
  rnorm = -1.d0
 
757
  do iel = 1, ncel
 
758
    rnorm  = max(rnorm,abs(rtp(iel,ipr)))
 
759
  enddo
 
760
  if (irangp.ge.0) call parmax (rnorm)
 
761
                   !==========
 
762
  write(nfecra,2100)rnorm
 
763
 
 
764
  rnorm = -1.d0
 
765
  do iel = 1, ncel
 
766
    vitnor =                                                    &
 
767
     sqrt(vel(1,iel)**2+vel(2,iel)**2+vel(3,iel)**2)
 
768
    if(vitnor.ge.rnorm) then
 
769
      rnorm = vitnor
 
770
      imax  = iel
 
771
    endif
 
772
  enddo
 
773
 
 
774
  xyzmax(1) = xyzcen(1,imax)
 
775
  xyzmax(2) = xyzcen(2,imax)
 
776
  xyzmax(3) = xyzcen(3,imax)
 
777
 
 
778
  if (irangp.ge.0) then
 
779
    nbrval = 3
 
780
    call parmxl (nbrval, rnorm, xyzmax)
 
781
    !==========
 
782
  endif
 
783
 
 
784
  write(nfecra,2200) rnorm,xyzmax(1),xyzmax(2),xyzmax(3)
 
785
 
 
786
 
 
787
  ! Pour la periodicite et le parallelisme, rom est echange dans phyvar
 
788
 
 
789
 
 
790
  rnorma = -grand
 
791
  rnormi =  grand
 
792
  do ifac = 1, nfac
 
793
    iel1 = ifacel(1,ifac)
 
794
    iel2 = ifacel(2,ifac)
 
795
    surf = surfan(ifac)
 
796
    rhom = (propce(iel1,ipcrom)+propce(iel2,ipcrom))*0.5d0
 
797
    rnorm = propfa(ifac,iflmas)/(surf*rhom)
 
798
    rnorma = max(rnorma,rnorm)
 
799
    rnormi = min(rnormi,rnorm)
 
800
  enddo
 
801
  if (irangp.ge.0) then
 
802
    call parmax (rnorma)
 
803
    !==========
 
804
    call parmin (rnormi)
 
805
    !==========
 
806
  endif
 
807
  write(nfecra,2300)rnorma, rnormi
 
808
 
 
809
  rnorma = -grand
 
810
  rnormi =  grand
 
811
  do ifac = 1, nfabor
 
812
    rnorm = propfb(ifac,iflmab)/(surfbn(ifac)*propfb(ifac,ipbrom))
 
813
    rnorma = max(rnorma,rnorm)
 
814
    rnormi = min(rnormi,rnorm)
 
815
  enddo
 
816
  if (irangp.ge.0) then
 
817
    call parmax (rnorma)
 
818
    !==========
 
819
    call parmin (rnormi)
 
820
    !==========
 
821
  endif
 
822
  write(nfecra,2400)rnorma, rnormi
 
823
 
 
824
  rnorm = 0.d0
 
825
  do ifac = 1, nfabor
 
826
    rnorm = rnorm + propfb(ifac,iflmab)
 
827
  enddo
 
828
 
 
829
  if (irangp.ge.0) call parsom (rnorm)
 
830
                   !==========
 
831
 
 
832
  write(nfecra,2500)rnorm
 
833
 
 
834
  write(nfecra,2001)
 
835
 
 
836
  if(nterup.gt.1) then
 
837
    if(icvrge.eq.0) then
 
838
      write(nfecra,2600) iterns
 
839
      write(nfecra,2601) xnrmu, xnrmu0, epsup
 
840
      write(nfecra,2001)
 
841
      if(iterns.eq.nterup) then
 
842
        write(nfecra,2603)
 
843
        write(nfecra,2001)
 
844
      endif
 
845
    else
 
846
      write(nfecra,2602) iterns
 
847
      write(nfecra,2601) xnrmu, xnrmu0, epsup
 
848
      write(nfecra,2001)
 
849
    endif
 
850
  endif
 
851
 
 
852
endif
 
853
 
 
854
! Free memory
 
855
deallocate(viscf, viscb)
 
856
deallocate(drtp, smbr, rovsdt)
 
857
deallocate(trav)
 
858
if (allocated(dfrcxt)) deallocate(dfrcxt)
 
859
if (allocated(frchy))  deallocate(frchy, dfrchy)
 
860
if (allocated(esflum)) deallocate(esflum, esflub)
 
861
if (allocated(wvisfi)) deallocate(wvisfi, wvisbi)
 
862
deallocate(w1, w2, w3)
 
863
deallocate(w4, w5, w6)
 
864
deallocate(w7, w8, w9)
 
865
if (allocated(w10)) deallocate(w10)
 
866
 
 
867
! Interleaved values of vel and vela
 
868
 
 
869
do iel = 1, ncelet
 
870
  rtp (iel,iu) = vel (1,iel)
 
871
  rtp (iel,iv) = vel (2,iel)
 
872
  rtp (iel,iw) = vel (3,iel)
 
873
  rtpa(iel,iu) = vela(1,iel)
 
874
  rtpa(iel,iv) = vela(2,iel)
 
875
  rtpa(iel,iw) = vela(3,iel)
 
876
  if (ipucou.eq.1 .or. ncpdct.gt.0) then
 
877
    tpucou(iel,1) = tpucov(1,iel)
 
878
    tpucou(iel,2) = tpucov(2,iel)
 
879
    tpucou(iel,3) = tpucov(3,iel)
 
880
  endif
 
881
enddo
 
882
 
 
883
! Free memory
 
884
!--------------
 
885
deallocate(vel)
 
886
deallocate(vela)
 
887
if(ipucou.eq.1 .or. ncpdct.gt.0) deallocate(tpucov)
 
888
 
 
889
!--------
 
890
! FORMATS
 
891
!--------
 
892
#if defined(_CS_LANG_FR)
 
893
 
 
894
 1000 format(/,                                                   &
 
895
'   ** RESOLUTION POUR LA VITESSE                             ',/,&
 
896
'      --------------------------                             ',/)
 
897
 1200 format(/,                                                   &
 
898
'   ** RESOLUTION POUR LA PRESSION CONTINUITE                 ',/,&
 
899
'      --------------------------------------                 ',/)
 
900
 2000 format(/,' APRES PRESSION CONTINUITE',/,                    &
 
901
'-------------------------------------------------------------'  )
 
902
 2100 format(                                                           &
 
903
' Pression max.',E12.4   ,' (max. de la valeur absolue)       ',/)
 
904
 2200 format(                                                           &
 
905
' Vitesse  max.',E12.4   ,' en',3E11.3                         ,/)
 
906
 2300 format(                                                           &
 
907
' Vitesse  en face interne max.',E12.4   ,' ; min.',E12.4        )
 
908
 2400 format(                                                           &
 
909
' Vitesse  en face de bord max.',E12.4   ,' ; min.',E12.4        )
 
910
 2500 format(                                                           &
 
911
' Bilan de masse   au bord   ',E14.6                             )
 
912
 2600 format(                                                           &
 
913
' Informations Point fixe a l''iteration :',I10                ,/)
 
914
 2601 format('norme = ',E12.4,' norme 0 = ',E12.4,' toler  = ',E12.4 ,/)
 
915
 2602 format(                                                           &
 
916
' Convergence du point fixe a l''iteration ',I10               ,/)
 
917
 2603 format(                                                           &
 
918
' Non convergence du couplage vitesse pression par point fixe  ' )
 
919
 2001 format(                                                           &
 
920
'-------------------------------------------------------------',/)
 
921
 
 
922
#else
 
923
 
 
924
 1000 format(/,                                                   &
 
925
'   ** SOLVING VELOCITY'                                       ,/,&
 
926
'      ----------------'                                       ,/)
 
927
 1200 format(/,                                                   &
 
928
'   ** SOLVING CONTINUITY PRESSURE'                            ,/,&
 
929
'      ---------------------------'                            ,/)
 
930
 2000 format(/,' AFTER CONTINUITY PRESSURE',/,                    &
 
931
'-------------------------------------------------------------'  )
 
932
 2100 format(                                                           &
 
933
' Max. pressure',E12.4   ,' (max. absolute value)'             ,/)
 
934
 2200 format(                                                           &
 
935
' Max. velocity',E12.4   ,' en',3E11.3                         ,/)
 
936
 2300 format(                                                           &
 
937
' Max. velocity at interior face',E12.4   ,' ; min.',E12.4       )
 
938
 2400 format(                                                           &
 
939
' Max. velocity at boundary face',E12.4   ,' ; min.',E12.4       )
 
940
 2500 format(                                                           &
 
941
' Mass balance  at boundary  ',E14.6                             )
 
942
 2600 format(                                                           &
 
943
' Fixed point informations at iteration:',I10                  ,/)
 
944
 2601 format('norm = ',E12.4,' norm 0 = ',E12.4,' toler  = ',E12.4   ,/)
 
945
 2602 format(                                                           &
 
946
' Fixed point convergence at iteration ',I10                   ,/)
 
947
 2603 format(                                                           &
 
948
' Non convergence of fixed point for velocity pressure coupling' )
 
949
 2001 format(                                                           &
 
950
'-------------------------------------------------------------',/)
 
951
 
 
952
#endif
 
953
 
 
954
!----
 
955
! FIN
 
956
!----
 
957
 
 
958
return
 
959
 
 
960
end subroutine