← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/base/itrgrp.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/itrgrp.f90
1
 
!-------------------------------------------------------------------------------
2
 
 
3
 
!     This file is part of the Code_Saturne Kernel, element of the
4
 
!     Code_Saturne CFD tool.
5
 
 
6
 
!     Copyright (C) 1998-2009 EDF S.A., France
7
 
 
8
 
!     contact: saturne-support@edf.fr
9
 
 
10
 
!     The Code_Saturne Kernel is free software; you can redistribute it
11
 
!     and/or modify it under the terms of the GNU General Public License
12
 
!     as published by the Free Software Foundation; either version 2 of
13
 
!     the License, or (at your option) any later version.
14
 
 
15
 
!     The Code_Saturne Kernel is distributed in the hope that it will be
16
 
!     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
17
 
!     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18
 
!     GNU General Public License for more details.
19
 
 
20
 
!     You should have received a copy of the GNU General Public License
21
 
!     along with the Code_Saturne Kernel; if not, write to the
22
 
!     Free Software Foundation, Inc.,
23
 
!     51 Franklin St, Fifth Floor,
24
 
!     Boston, MA  02110-1301  USA
25
 
 
26
 
!-------------------------------------------------------------------------------
27
 
 
28
 
subroutine itrgrp &
29
 
!================
30
 
 
31
 
 ( idbia0 , idbra0 ,                                              &
32
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
33
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
34
 
   nvar   , nscal  , nphas  ,                                     &
35
 
   nideve , nrdeve , nituse , nrtuse ,                            &
36
 
   init   , inc    , imrgra , iccocg , nswrgp , imligp ,          &
37
 
   iphydp , iwarnp , nfecra ,                                     &
38
 
   epsrgp , climgp , extrap ,                                     &
39
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
40
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
41
 
   idevel , ituser , ia     ,                                     &
42
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
43
 
   fextx  , fexty  , fextz  ,                                     &
44
 
   pvar   , coefap , coefbp , viscf  , viscb  ,                   &
45
 
   viselx , visely , viselz ,                                     &
46
 
   diverg ,                                                       &
47
 
   dpdx   , dpdy   , dpdz   , dpdxa  , dpdya  , dpdza  ,          &
48
 
   rdevel , rtuser , ra     )
49
 
 
50
 
!===============================================================================
51
 
! FONCTION :
52
 
! ----------
53
 
 
54
 
! INCREMENTATION DE LA DIVERGENCE DU FLUX DE MASSE
55
 
! A PARTIR DE GRAD(P)
56
 
! grad(P) = GRADIENT FACETTE
57
 
 
58
 
!  .          .        --         --->     -->
59
 
!  m   =  m     - \    Visc   grad(P) . n
60
 
!   ij     ij     /__      ij       ij   ij
61
 
!                  j
62
 
 
63
 
!-------------------------------------------------------------------------------
64
 
! Arguments
65
 
!__________________.____._____.________________________________________________.
66
 
! name             !type!mode ! role                                           !
67
 
!__________________!____!_____!________________________________________________!
68
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
69
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
70
 
! ndim             ! i  ! <-- ! spatial dimension                              !
71
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
72
 
! ncel             ! i  ! <-- ! number of cells                                !
73
 
! nfac             ! i  ! <-- ! number of interior faces                       !
74
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
75
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
76
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
77
 
! nnod             ! i  ! <-- ! number of vertices                             !
78
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
79
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
80
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
81
 
! nvar             ! i  ! <-- ! total number of variables                      !
82
 
! nscal            ! i  ! <-- ! total number of scalars                        !
83
 
! nphas            ! i  ! <-- ! number of phases                               !
84
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
85
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
86
 
! init             ! e  ! <-- ! > 0 : initialisation du flux de masse          !
87
 
! inc              ! e  ! <-- ! indicateur = 0 resol sur increment             !
88
 
!                  !    !     !              1 sinon                           !
89
 
! imrgra           ! e  ! <-- ! indicateur = 0 gradrc 97                       !
90
 
!                  ! e  ! <-- !            = 1 gradmc 99                       !
91
 
! iccocg           ! e  ! <-- ! indicateur = 1 pour recalcul de cocg           !
92
 
!                  !    !     !              0 sinon                           !
93
 
! nswrgp           ! e  ! <-- ! nombre de sweep pour reconstruction            !
94
 
!                  !    !     !             des gradients                      !
95
 
! imligp           ! e  ! <-- ! methode de limitation du gradient              !
96
 
!                  !    !     !  < 0 pas de limitation                         !
97
 
!                  !    !     !  = 0 a partir des gradients voisins            !
98
 
!                  !    !     !  = 1 a partir du gradient moyen                !
99
 
! iwarnp           ! i  ! <-- ! verbosity                                      !
100
 
! iphydp           ! e  ! <-- ! indicateur de prise en compte de la            !
101
 
!                  !    !     ! pression hydrostatique                         !
102
 
! nfecra           ! e  ! <-- ! unite du fichier sortie std                    !
103
 
! epsrgp           ! r  ! <-- ! precision relative pour la                     !
104
 
!                  !    !     !  reconstruction des gradients 97               !
105
 
! climgp           ! r  ! <-- ! coef gradient*distance/ecart                   !
106
 
! extrap           ! r  ! <-- ! coef extrap gradient                           !
107
 
! pvar  (ncelet    ! tr ! <-- ! variable (pression)                            !
108
 
! coefap, b        ! tr ! <-- ! tableaux des cond lim pour pvar                !
109
 
!   (nfabor)       !    !     !  sur la normale a la face de bord              !
110
 
! viscf (nfac)     ! tr ! <-- ! "viscosite" face interne(dt*surf/dist          !
111
 
! viscb (nfabor    ! tr ! <-- ! "viscosite" face de bord(dt*surf/dist          !
112
 
! viselx(ncelet    ! tr ! <-- ! "viscosite" par cellule  dir x                 !
113
 
! visely(ncelet    ! tr ! <-- ! "viscosite" par cellule  dir y                 !
114
 
! viselz(ncelet    ! tr ! <-- ! "viscosite" par cellule  dir z                 !
115
 
! diverg(ncelet    ! tr ! <-- ! divergence du flux de masse                    !
116
 
! fextx,y,z        ! tr ! <-- ! force exterieure generant la pression          !
117
 
!   (ncelet)       !    !     !  hydrostatique                                 !
118
 
! dpd.(ncelet      ! tr ! --- ! tableau de travail pour le grad de p           !
119
 
! dpdxa,y,z        ! tr ! --- ! tableau de travail pour le grad de p           !
120
 
!    (ncelet       !    !     !  avec decentrement amont                       !
121
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary development      !
122
 
! rtuser(nrtuse)   ! ra ! <-> ! user-reserved real work array                  !
123
 
!__________________!____!_____!________________________________________________!
124
 
 
125
 
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
126
 
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
127
 
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
128
 
!            --- tableau de travail
129
 
 
130
 
!-------------------------------------------------------------------------------
131
 
!===============================================================================
132
 
 
133
 
implicit none
134
 
 
135
 
!===============================================================================
136
 
! Common blocks
137
 
!===============================================================================
138
 
 
139
 
include "paramx.h"
140
 
include "pointe.h"
141
 
include "numvar.h"
142
 
include "period.h"
143
 
include "parall.h"
144
 
 
145
 
!===============================================================================
146
 
 
147
 
! Arguments
148
 
 
149
 
integer          idbia0 , idbra0
150
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
151
 
integer          nfml   , nprfml
152
 
integer          nnod   , lndfac , lndfbr , ncelbr
153
 
integer          nvar   , nscal  , nphas
154
 
integer          nideve , nrdeve , nituse , nrtuse
155
 
integer          init   , inc    , imrgra , iccocg
156
 
integer          nswrgp , imligp
157
 
integer          iwarnp , iphydp , nfecra
158
 
double precision epsrgp , climgp , extrap
159
 
 
160
 
integer          ifacel(2,nfac) , ifabor(nfabor)
161
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
162
 
integer          iprfml(nfml,nprfml)
163
 
integer          ipnfac(nfac+1), nodfac(lndfac)
164
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
165
 
integer          idevel(nideve), ituser(nituse)
166
 
integer          ia(*)
167
 
 
168
 
double precision xyzcen(ndim,ncelet)
169
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
170
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
171
 
double precision xyznod(ndim,nnod), volume(ncelet)
172
 
double precision pvar(ncelet), coefap(nfabor), coefbp(nfabor)
173
 
double precision viscf(nfac), viscb(nfabor)
174
 
double precision viselx(ncelet), visely(ncelet), viselz(ncelet)
175
 
double precision diverg(ncelet)
176
 
double precision fextx(ncelet),fexty(ncelet),fextz(ncelet)
177
 
double precision dpdx (ncelet),dpdy (ncelet),dpdz (ncelet)
178
 
double precision dpdxa(ncelet),dpdya(ncelet),dpdza(ncelet)
179
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
180
 
 
181
 
! Local variables
182
 
 
183
 
integer          idebia, idebra
184
 
integer          ifac, ii, jj, iij, iii, ivar
185
 
integer          idimte, itenso
186
 
double precision pfac,pip
187
 
double precision dpxf  , dpyf  , dpzf  , flumas, flumab
188
 
double precision dijpfx, dijpfy, dijpfz
189
 
double precision diipbx, diipby, diipbz
190
 
double precision dist  , surfan
191
 
double precision dijx  , dijy  , dijz
192
 
 
193
 
!===============================================================================
194
 
 
195
 
!===============================================================================
196
 
! 1. INITIALISATION
197
 
!===============================================================================
198
 
 
199
 
idebia = idbia0
200
 
idebra = idbra0
201
 
 
202
 
if( init.ge.1 ) then
203
 
  do ii = 1, ncelet
204
 
    diverg(ii) = 0.d0
205
 
  enddo
206
 
elseif(init.eq.0.and.ncelet.gt.ncel) then
207
 
  do ii = ncel+1, ncelet
208
 
    diverg(ii) = 0.d0
209
 
  enddo
210
 
elseif(init.ne.0) then
211
 
  write(nfecra,1000) init
212
 
  call csexit (1)
213
 
endif
214
 
 
215
 
! ---> TRAITEMENT DU PARALLELISME
216
 
 
217
 
if(irangp.ge.0) call parcom (pvar)
218
 
                !==========
219
 
 
220
 
! ---> TRAITEMENT DE LA PERIODICITE
221
 
 
222
 
if(iperio.eq.1) then
223
 
  idimte = 0
224
 
  itenso = 0
225
 
  call percom                                                     &
226
 
  !==========
227
 
  ( idimte , itenso ,                                             &
228
 
    pvar   , pvar   , pvar  ,                                     &
229
 
    pvar   , pvar   , pvar  ,                                     &
230
 
    pvar   , pvar   , pvar  )
231
 
endif
232
 
 
233
 
!===============================================================================
234
 
! 2.  INCREMENT DU FLUX DE MASSE SS TECHNIQUE DE RECONSTRUCTION
235
 
!===============================================================================
236
 
 
237
 
if( nswrgp.le.1 ) then
238
 
 
239
 
!     FLUX DE MASSE SUR LES FACETTES FLUIDES
240
 
 
241
 
  do ifac = 1, nfac
242
 
 
243
 
    ii = ifacel(1,ifac)
244
 
    jj = ifacel(2,ifac)
245
 
 
246
 
    flumas = viscf(ifac)*(pvar(ii) -pvar(jj))
247
 
    diverg(ii) = diverg(ii) + flumas
248
 
    diverg(jj) = diverg(jj) - flumas
249
 
 
250
 
  enddo
251
 
 
252
 
!     FLUX DE MASSE SUR LES FACETTES DE BORD
253
 
 
254
 
  do ifac = 1, nfabor
255
 
 
256
 
    ii = ifabor(ifac)
257
 
    pfac = inc*coefap(ifac) +coefbp(ifac)*pvar(ii)
258
 
 
259
 
    flumab = viscb(ifac)*( pvar(ii) -pfac )
260
 
    diverg(ii) = diverg(ii) + flumab
261
 
 
262
 
  enddo
263
 
 
264
 
endif
265
 
 
266
 
 
267
 
!===============================================================================
268
 
! 3.  INCREMENTATION DU FLUX DE MASSE AVEC TECHNIQUE DE
269
 
!         RECONSTRUCTION SI LE MAILLAGE EST NON ORTHOGONAL
270
 
!===============================================================================
271
 
 
272
 
if( nswrgp.gt.1 ) then
273
 
 
274
 
!     CALCUL DU GRADIENT
275
 
 
276
 
!     IVAR ne sert a GRDCEL que si la variable est une composante de la vitesse
277
 
!     ou de Rij pour la periodicite. Ici la variable est soit la pression
278
 
!     soit phi, donc on peut mettre IVAR=0
279
 
  ivar = 0
280
 
 
281
 
  call grdcel                                                     &
282
 
  !==========
283
 
 ( idebia , idebra ,                                              &
284
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
285
 
   nnod   , lndfac , lndfbr , ncelbr , nphas  ,                   &
286
 
   nideve , nrdeve , nituse , nrtuse ,                            &
287
 
   ivar   , imrgra , inc    , iccocg , nswrgp , imligp , iphydp , &
288
 
 
289
 
   iwarnp , nfecra , epsrgp , climgp , extrap ,                   &
290
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
291
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
292
 
   idevel , ituser , ia     ,                                     &
293
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
294
 
   fextx  , fexty  , fextz  ,                                     &
295
 
   pvar   , coefap , coefbp ,                                     &
296
 
   dpdx   , dpdy   , dpdz   ,                                     &
297
 
!        ------   ------   ------
298
 
   dpdxa  , dpdya  , dpdza  ,                                     &
299
 
   rdevel , rtuser , ra     )
300
 
 
301
 
! ---> TRAITEMENT DU PARALLELISME
302
 
 
303
 
if(irangp.ge.0) then
304
 
  call parcom (viselx)
305
 
  !==========
306
 
  call parcom (visely)
307
 
  !==========
308
 
  call parcom (viselz)
309
 
  !==========
310
 
endif
311
 
 
312
 
! ---> TRAITEMENT DE LA PERIODICITE
313
 
 
314
 
if(iperio.eq.1) then
315
 
  idimte = 1
316
 
  itenso = 0
317
 
  call percom                                                     &
318
 
  !==========
319
 
  ( idimte , itenso ,                                             &
320
 
    viselx , viselx , viselx ,                                    &
321
 
    visely , visely , visely ,                                    &
322
 
    viselz , viselz , viselz )
323
 
endif
324
 
 
325
 
!     FLUX DE MASSE SUR LES FACETTES FLUIDES
326
 
 
327
 
  do ifac = 1, nfac
328
 
 
329
 
    ii = ifacel(1,ifac)
330
 
    jj = ifacel(2,ifac)
331
 
 
332
 
    iij = idijpf-1+3*(ifac-1)
333
 
 
334
 
    dijpfx = ra(iij+1)
335
 
    dijpfy = ra(iij+2)
336
 
    dijpfz = ra(iij+3)
337
 
    dist   = ra(idist-1+ifac)
338
 
    surfan = ra(isrfan-1+ifac)
339
 
 
340
 
!---> DIJ = IJ - (IJ.N) N
341
 
    dijx = (xyzcen(1,jj)-xyzcen(1,ii))-dijpfx
342
 
    dijy = (xyzcen(2,jj)-xyzcen(2,ii))-dijpfy
343
 
    dijz = (xyzcen(3,jj)-xyzcen(3,ii))-dijpfz
344
 
 
345
 
    dpxf = 0.5d0*(viselx(ii)*dpdx(ii) + viselx(jj)*dpdx(jj))
346
 
    dpyf = 0.5d0*(visely(ii)*dpdy(ii) + visely(jj)*dpdy(jj))
347
 
    dpzf = 0.5d0*(viselz(ii)*dpdz(ii) + viselz(jj)*dpdz(jj))
348
 
 
349
 
    flumas = viscf(ifac)*( pvar(ii) -pvar(jj) )                   &
350
 
     + ( dpxf * dijx                                              &
351
 
     +   dpyf * dijy                                              &
352
 
     +   dpzf * dijz )*surfan/dist
353
 
    diverg(ii) = diverg(ii) + flumas
354
 
    diverg(jj) = diverg(jj) - flumas
355
 
 
356
 
  enddo
357
 
 
358
 
!     FLUX DE MASSE SUR LES FACETTES DE BORD
359
 
 
360
 
  do ifac = 1, nfabor
361
 
 
362
 
    ii = ifabor(ifac)
363
 
 
364
 
    iii = idiipb-1+3*(ifac-1)
365
 
    diipbx = ra(iii+1)
366
 
    diipby = ra(iii+2)
367
 
    diipbz = ra(iii+3)
368
 
 
369
 
    pip = pvar(ii) +                                              &
370
 
      dpdx(ii)*diipbx+dpdy(ii)*diipby+dpdz(ii)*diipbz
371
 
    pfac = inc*coefap(ifac) +coefbp(ifac)*pip
372
 
 
373
 
    flumab = viscb(ifac)*( pip -pfac )
374
 
    diverg(ii) = diverg(ii) + flumab
375
 
 
376
 
  enddo
377
 
 
378
 
endif
379
 
 
380
 
!--------
381
 
! FORMATS
382
 
!--------
383
 
 
384
 
#if defined(_CS_LANG_FR)
385
 
 
386
 
 1000 format('ITRGRP APPELE AVEC INIT = ',I10)
387
 
 
388
 
#else
389
 
 
390
 
 1000 format('ITRGRP CALLED WITH INIT = ',I10)
391
 
 
392
 
#endif
393
 
 
394
 
!----
395
 
! FIN
396
 
!----
397
 
 
398
 
return
399
 
 
400
 
end subroutine