← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/base/itrmas.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/itrmas.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 itrmas &
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
 
   flumas , flumab ,                                              &
47
 
   dpdx   , dpdy   , dpdz   , dpdxa  , dpdya  , dpdza  ,          &
48
 
   rdevel , rtuser , ra     )
49
 
 
50
 
!===============================================================================
51
 
! FONCTION :
52
 
! ----------
53
 
 
54
 
! INCREMENTATION DU FLUX DE MASSE A PARTIR DE GRAD(P)
55
 
! P = PRESSION, INCREMENT DE PRESSION, DOUBLE INCREMENT DE PRESSION
56
 
! grad(P) = GRADIENT FACETTE POUR L'INSTANT
57
 
 
58
 
!  .          .               --->     -->
59
 
!  m   =  m     -  Visc   grad(P) . n
60
 
!   ij     ij          ij       ij   ij
61
 
 
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
 
! flumas(nfac)     ! tr ! <-- ! flux de masse aux faces internes               !
116
 
! flumab(nfabor    ! tr ! <-- ! flux de masse aux faces de bord                !
117
 
! fextx,y,z        ! tr ! <-- ! force exterieure generant la pression          !
118
 
!   (ncelet)       !    !     !  hydrostatique                                 !
119
 
! dpd.(ncelet      ! tr ! --- ! tableau de travail pour le grad de p           !
120
 
! dpdxa,y,z        ! tr ! --- ! tableau de travail pour le grad de p           !
121
 
!    (ncelet       !    !     !  avec decentrement amont                       !
122
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary development      !
123
 
! rtuser(nrtuse)   ! ra ! <-> ! user-reserved real work array                  !
124
 
!__________________!____!_____!________________________________________________!
125
 
 
126
 
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
127
 
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
128
 
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
129
 
!            --- tableau de travail
130
 
 
131
 
!-------------------------------------------------------------------------------
132
 
!===============================================================================
133
 
 
134
 
implicit none
135
 
 
136
 
!===============================================================================
137
 
! Common blocks
138
 
!===============================================================================
139
 
 
140
 
include "paramx.h"
141
 
include "pointe.h"
142
 
include "numvar.h"
143
 
include "period.h"
144
 
include "parall.h"
145
 
 
146
 
!===============================================================================
147
 
 
148
 
! Arguments
149
 
 
150
 
integer          idbia0 , idbra0
151
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
152
 
integer          nfml   , nprfml
153
 
integer          nnod   , lndfac , lndfbr , ncelbr
154
 
integer          nvar   , nscal  , nphas
155
 
integer          nideve , nrdeve , nituse , nrtuse
156
 
integer          init   , inc    , imrgra , iccocg
157
 
integer          nswrgp , imligp
158
 
integer          iwarnp , iphydp , nfecra
159
 
double precision epsrgp , climgp , extrap
160
 
 
161
 
integer          ifacel(2,nfac) , ifabor(nfabor)
162
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
163
 
integer          iprfml(nfml,nprfml)
164
 
integer          ipnfac(nfac+1), nodfac(lndfac)
165
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
166
 
integer          idevel(nideve), ituser(nituse)
167
 
integer          ia(*)
168
 
 
169
 
double precision xyzcen(ndim,ncelet)
170
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
171
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
172
 
double precision xyznod(ndim,nnod), volume(ncelet)
173
 
double precision pvar(ncelet), coefap(nfabor), coefbp(nfabor)
174
 
double precision viscf(nfac), viscb(nfabor)
175
 
double precision viselx(ncelet), visely(ncelet), viselz(ncelet)
176
 
double precision flumas(nfac), flumab(nfabor)
177
 
double precision fextx(ncelet),fexty(ncelet),fextz(ncelet)
178
 
double precision dpdx (ncelet),dpdy (ncelet),dpdz (ncelet)
179
 
double precision dpdxa(ncelet),dpdya(ncelet),dpdza(ncelet)
180
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
181
 
 
182
 
! Local variables
183
 
 
184
 
integer          idebia, idebra
185
 
integer          ifac, ii, jj, iij, iii
186
 
integer          idimte, itenso
187
 
double precision pfac,pip
188
 
double precision dpxf  , dpyf  , dpzf  ,surfan, dist
189
 
double precision dijpfx, dijpfy, dijpfz
190
 
double precision diipbx, diipby, diipbz
191
 
double precision dijx  , dijy  , dijz
192
 
 
193
 
!===============================================================================
194
 
 
195
 
 
196
 
!     C'est penible de passer viscf et viscel qui portent une info
197
 
!                similaire
198
 
 
199
 
 
200
 
 
201
 
!===============================================================================
202
 
! 1. INITIALISATION
203
 
!===============================================================================
204
 
 
205
 
idebia = idbia0
206
 
idebra = idbra0
207
 
 
208
 
if( init.ge.1 ) then
209
 
  do ifac = 1, nfac
210
 
    flumas(ifac) = 0.d0
211
 
  enddo
212
 
  do ifac = 1, nfabor
213
 
    flumab(ifac) = 0.d0
214
 
  enddo
215
 
elseif(init.ne.0) then
216
 
  write(nfecra,1000) init
217
 
  call csexit (1)
218
 
endif
219
 
 
220
 
! ---> TRAITEMENT DU PARALLELISME
221
 
 
222
 
if(irangp.ge.0) call parcom (pvar)
223
 
                !==========
224
 
 
225
 
! ---> TRAITEMENT DE LA PERIODICITE
226
 
 
227
 
if(iperio.eq.1) then
228
 
  idimte = 0
229
 
  itenso = 0
230
 
  call percom                                                     &
231
 
  !==========
232
 
  ( idimte , itenso ,                                             &
233
 
    pvar   , pvar   , pvar  ,                                     &
234
 
    pvar   , pvar   , pvar  ,                                     &
235
 
    pvar   , pvar   , pvar  )
236
 
endif
237
 
 
238
 
!===============================================================================
239
 
! 2.  INCREMENT DU FLUX DE MASSE SS TECHNIQUE DE RECONSTRUCTION
240
 
!===============================================================================
241
 
 
242
 
if( nswrgp.le.1 ) then
243
 
 
244
 
!     FLUX DE MASSE SUR LES FACETTES FLUIDES
245
 
 
246
 
  do ifac = 1, nfac
247
 
 
248
 
    ii = ifacel(1,ifac)
249
 
    jj = ifacel(2,ifac)
250
 
 
251
 
    flumas(ifac) = flumas(ifac) + viscf(ifac)*(pvar(ii) -pvar(jj))
252
 
 
253
 
  enddo
254
 
 
255
 
!     FLUX DE MASSE SUR LES FACETTES DE BORD
256
 
 
257
 
  do ifac = 1, nfabor
258
 
 
259
 
    ii = ifabor(ifac)
260
 
    pfac = inc*coefap(ifac) +coefbp(ifac)*pvar(ii)
261
 
 
262
 
    flumab(ifac) = flumab(ifac) +viscb(ifac)*( pvar(ii) -pfac )
263
 
 
264
 
  enddo
265
 
 
266
 
endif
267
 
 
268
 
 
269
 
!===============================================================================
270
 
! 3.  INCREMENTATION DU FLUX DE MASSE AVEC TECHNIQUE DE
271
 
!         RECONSTRUCTION SI LE MAILLAGE EST NON ORTHOGONAL
272
 
!===============================================================================
273
 
 
274
 
if( nswrgp.gt.1 ) then
275
 
 
276
 
!     CALCUL DU GRADIENT
277
 
 
278
 
  call grdcel                                                     &
279
 
  !==========
280
 
 ( idebia , idebra ,                                              &
281
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
282
 
   nnod   , lndfac , lndfbr , ncelbr , nphas  ,                   &
283
 
   nideve , nrdeve , nituse , nrtuse ,                            &
284
 
   ipr(1) , imrgra , inc    , iccocg , nswrgp , imligp , iphydp , &
285
 
   iwarnp , nfecra , epsrgp , climgp , extrap ,                   &
286
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
287
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
288
 
   idevel , ituser , ia     ,                                     &
289
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
290
 
   fextx  , fexty  , fextz  ,                                     &
291
 
   pvar   , coefap , coefbp ,                                     &
292
 
   dpdx   , dpdy   , dpdz   ,                                     &
293
 
!        ------   ------   ------
294
 
   dpdxa  , dpdya  , dpdza  ,                                     &
295
 
   rdevel , rtuser , ra     )
296
 
 
297
 
! ---> TRAITEMENT DU PARALLELISME
298
 
 
299
 
if(irangp.ge.0) then
300
 
  call parcom (viselx)
301
 
  !==========
302
 
  call parcom (visely)
303
 
  !==========
304
 
  call parcom (viselz)
305
 
  !==========
306
 
endif
307
 
 
308
 
! ---> TRAITEMENT DE LA PERIODICITE
309
 
 
310
 
if(iperio.eq.1) then
311
 
  idimte = 1
312
 
  itenso = 0
313
 
  call percom                                                     &
314
 
  !==========
315
 
  ( idimte , itenso ,                                             &
316
 
    viselx , viselx , viselx ,                                    &
317
 
    visely , visely , visely ,                                    &
318
 
    viselz , viselz , viselz )
319
 
endif
320
 
 
321
 
!     FLUX DE MASSE SUR LES FACETTES FLUIDES
322
 
 
323
 
  do ifac = 1, nfac
324
 
 
325
 
    ii = ifacel(1,ifac)
326
 
    jj = ifacel(2,ifac)
327
 
 
328
 
    dpxf = 0.5d0*(viselx(ii)*dpdx(ii) + viselx(jj)*dpdx(jj))
329
 
    dpyf = 0.5d0*(visely(ii)*dpdy(ii) + visely(jj)*dpdy(jj))
330
 
    dpzf = 0.5d0*(viselz(ii)*dpdz(ii) + viselz(jj)*dpdz(jj))
331
 
 
332
 
    iij = idijpf-1+3*(ifac-1)
333
 
    dijpfx = ra(iij+1)
334
 
    dijpfy = ra(iij+2)
335
 
    dijpfz = ra(iij+3)
336
 
 
337
 
!---> DIJ = IJ - (IJ.N) N
338
 
    dijx = (xyzcen(1,jj)-xyzcen(1,ii))-dijpfx
339
 
    dijy = (xyzcen(2,jj)-xyzcen(2,ii))-dijpfy
340
 
    dijz = (xyzcen(3,jj)-xyzcen(3,ii))-dijpfz
341
 
 
342
 
    dist   = ra(idist-1+ifac)
343
 
    surfan = ra(isrfan-1+ifac)
344
 
 
345
 
    flumas(ifac) = flumas(ifac)                                   &
346
 
     + viscf(ifac)*( pvar(ii) -pvar(jj) )                         &
347
 
     + ( dpxf * dijx                                              &
348
 
     +   dpyf * dijy                                              &
349
 
     +   dpzf * dijz )*surfan/dist
350
 
 
351
 
  enddo
352
 
 
353
 
!     FLUX DE MASSE SUR LES FACETTES DE BORD
354
 
 
355
 
  do ifac = 1, nfabor
356
 
 
357
 
    ii = ifabor(ifac)
358
 
    iii = idiipb-1+3*(ifac-1)
359
 
    diipbx = ra(iii+1)
360
 
    diipby = ra(iii+2)
361
 
    diipbz = ra(iii+3)
362
 
 
363
 
    pip = pvar(ii)                                                &
364
 
        + dpdx(ii)*diipbx                                         &
365
 
        + dpdy(ii)*diipby + dpdz(ii)*diipbz
366
 
    pfac = inc*coefap(ifac) +coefbp(ifac)*pip
367
 
 
368
 
    flumab(ifac) = flumab(ifac) +viscb(ifac)*( pip -pfac )
369
 
 
370
 
  enddo
371
 
 
372
 
endif
373
 
 
374
 
!--------
375
 
! FORMATS
376
 
!--------
377
 
 
378
 
#if defined(_CS_LANG_FR)
379
 
 
380
 
 1000 format('ITRMAS APPELE AVEC INIT = ',I10)
381
 
 
382
 
#else
383
 
 
384
 
 1000 format('ITRMAS CALLED WITH INIT = ',I10)
385
 
 
386
 
#endif
387
 
 
388
 
!----
389
 
! FIN
390
 
!----
391
 
 
392
 
return
393
 
 
394
 
end subroutine