← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/alge/inimav.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/alge/inimav.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 inimav &
 
24
!================
 
25
 
 
26
 ( nvar   , nscal  ,                                              &
 
27
   ivar   ,                                                       &
 
28
   iflmb0 , init   , inc    , imrgra , iccocg , nswrgu , imligu , &
 
29
   iwarnu , nfecra ,                                              &
 
30
   epsrgu , climgu , extrau ,                                     &
 
31
   rom    , romb   ,                                              &
 
32
   vel    ,                                                       &
 
33
   coefav , coefbv ,                                              &
 
34
   flumas , flumab )
 
35
 
 
36
!===============================================================================
 
37
! FONCTION :
 
38
! ----------
 
39
!                                                              -->
 
40
! INCREMENTATION DU FLUX DE MASSE A PARTIR DU CHAMP VECTORIEL ROM.U
 
41
!  .     .          -->   -->
 
42
!  m   = m   +(rom* U ) . n
 
43
!   ij    ij         ij    ij
 
44
 
 
45
 
 
46
! Pour la reconstruction, grad(rho u) est calcule avec des
 
47
! conditions aux limites approchees :
 
48
!    COEFA(rho u) = ROMB * COEFA(u)
 
49
!    COEFB(rho u) = COEFB (u)
 
50
 
 
51
! et pour le flux de masse au bord on ecrit
 
52
!  FLUMAB = [ROMB*COEFA(u) + ROMB*COEFB(u)*Ui
 
53
!                + COEFB(u)*II'.grad(rho u) ].Sij
 
54
! ce qui utilise de petites approximations sur les
 
55
! non-orthogonalites (cf. notice)
 
56
 
 
57
!-------------------------------------------------------------------------------
 
58
! Arguments
 
59
!__________________.____._____.________________________________________________.
 
60
! name             !type!mode ! role                                           !
 
61
!__________________!____!_____!________________________________________________!
 
62
! nvar             ! i  ! <-- ! total number of variables                      !
 
63
! nscal            ! i  ! <-- ! total number of scalars                        !
 
64
! ivar             ! e  ! <-- ! variable                                       !
 
65
! iflmb0           ! e  ! <-- ! =1 : flux de masse annule sym-paroi            !
 
66
! init             ! e  ! <-- ! > 0 : initialisation du flux de masse          !
 
67
! inc              ! e  ! <-- ! indicateur = 0 resol sur increment             !
 
68
!                  !    !     !              1 sinon                           !
 
69
! imrgra           ! e  ! <-- ! indicateur = 0 gradrc 97                       !
 
70
!                  ! e  ! <-- !            = 1 gradmc 99                       !
 
71
! nswrgu           ! e  ! <-- ! nombre de sweep pour reconstruction            !
 
72
!                  !    !     !             des gradients                      !
 
73
! imligu           ! e  ! <-- ! methode de limitation du gradient              !
 
74
!                  !    !     !  < 0 pas de limitation                         !
 
75
!                  !    !     !  = 0 a partir des gradients voisins            !
 
76
!                  !    !     !  = 1 a partir du gradient moyen                !
 
77
! iwarnu           ! e  ! <-- ! niveau d'impression                            !
 
78
! nfecra           ! e  ! <-- ! unite du fichier sortie std                    !
 
79
! epsrgu           ! r  ! <-- ! precision relative pour la                     !
 
80
!                  !    !     !  reconstruction des gradients 97               !
 
81
! climgu           ! r  ! <-- ! coef gradient*distance/ecart                   !
 
82
! extrau           ! r  ! <-- ! coef extrap gradient                           !
 
83
! isympa           ! te ! <-- ! zero pour annuler le flux de masse             !
 
84
! (nfabor     )    !    !     !(symetries et parois avec cl couplees)          !
 
85
!                  !    !     ! un sinon                                       !
 
86
! ia(*)            ! ia ! --- ! main integer work array                        !
 
87
! rom(ncelet       ! tr ! <-- ! masse volumique aux cellules                   !
 
88
! romb(nfabor)     ! tr ! <-- ! masse volumique aux bords                      !
 
89
! vel(3,ncelet)    ! tr ! <-- ! vitesse                                        !
 
90
! coefav, b        ! tr ! <-- ! tableaux des cond lim pour ux, uy, uz          !
 
91
!   (3,nfabor)     !    !     !  sur la normale a la face de bord              !
 
92
! coefbv, q        ! tr ! <-- ! tableaux des cond lim pour ux, uy, uz          !
 
93
!   (3,3,nfabor)   !    !     !  sur la normale a la face de bord              !
 
94
! flumas(nfac)     ! tr ! <-- ! flux de masse aux faces internes               !
 
95
! flumab(nfabor    ! tr ! <-- ! flux de masse aux faces de bord                !
 
96
!__________________!____!_____!________________________________________________!
 
97
 
 
98
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
 
99
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
 
100
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
 
101
!            --- tableau de travail
 
102
!===============================================================================
 
103
 
 
104
!===============================================================================
 
105
! Module files
 
106
!===============================================================================
 
107
 
 
108
use paramx
 
109
use dimens, only: ndimfb
 
110
use pointe
 
111
use parall
 
112
use period
 
113
use mesh
 
114
 
 
115
!===============================================================================
 
116
 
 
117
implicit none
 
118
 
 
119
! Arguments
 
120
 
 
121
integer          nvar   , nscal
 
122
integer          ivar
 
123
integer          iflmb0 , init   , inc    , imrgra , iccocg
 
124
integer          nswrgu , imligu
 
125
integer          iwarnu , nfecra
 
126
double precision epsrgu , climgu , extrau
 
127
 
 
128
 
 
129
double precision rom(ncelet), romb(nfabor)
 
130
double precision vel(3,ncelet)
 
131
double precision coefav(3,ndimfb)
 
132
double precision coefbv(3,3,nfabor)
 
133
double precision flumas(nfac), flumab(nfabor)
 
134
 
 
135
! Local variables
 
136
 
 
137
integer          ifac, ii, jj, iel
 
138
integer          iappel, isou, jsou
 
139
double precision pfac, pip
 
140
double precision dofx,dofy,dofz,pnd
 
141
double precision diipbx, diipby, diipbz
 
142
logical          ilved
 
143
 
 
144
double precision, dimension(:,:), allocatable :: qdm, coefaq
 
145
double precision, dimension(:,:,:), allocatable :: grdqdm
 
146
 
 
147
allocate(qdm(3,ncelet))
 
148
allocate(coefaq(3,ndimfb))
 
149
 
 
150
!===============================================================================
 
151
 
 
152
!===============================================================================
 
153
! 1.  INITIALISATION
 
154
!===============================================================================
 
155
 
 
156
! ---> CALCUL DE LA QTE DE MOUVEMENT
 
157
 
 
158
 
 
159
if( init.eq.1 ) then
 
160
  do ifac = 1, nfac
 
161
    flumas(ifac) = 0.d0
 
162
  enddo
 
163
  do ifac = 1, nfabor
 
164
    flumab(ifac) = 0.d0
 
165
  enddo
 
166
 
 
167
elseif(init.ne.0) then
 
168
  write(nfecra,1000) init
 
169
  call csexit (1)
 
170
endif
 
171
 
 
172
do iel = 1, ncel
 
173
  do isou = 1, 3
 
174
    qdm(isou,iel) = rom(iel)*vel(isou,iel)
 
175
  enddo
 
176
enddo
 
177
 
 
178
! ---> TRAITEMENT DU PARALLELISME ET DE LA PERIODICITE
 
179
 
 
180
if (irangp.ge.0.or.iperio.eq.1) then
 
181
  call synvin(qdm)
 
182
  !==========
 
183
endif
 
184
 
 
185
do ifac =1, nfabor
 
186
  do isou = 1, 3
 
187
    coefaq(isou,ifac) = romb(ifac)*coefav(isou,ifac)
 
188
  enddo
 
189
enddo
 
190
!===============================================================================
 
191
! 2.  CALCUL DU FLUX DE MASSE SANS TECHNIQUE DE RECONSTRUCTION
 
192
!===============================================================================
 
193
 
 
194
if( nswrgu.le.1 ) then
 
195
 
 
196
!     FLUX DE MASSE SUR LES FACETTES FLUIDES
 
197
 
 
198
  do ifac = 1, nfac
 
199
    ii = ifacel(1,ifac)
 
200
    jj = ifacel(2,ifac)
 
201
    pnd = pond(ifac)
 
202
    ! Components U, V, W
 
203
    do isou = 1, 3
 
204
      flumas(ifac) = flumas(ifac) +                                        &
 
205
         (pnd*qdm(isou,ii)+(1.d0-pnd)*qdm(isou,jj)) *surfac(isou,ifac)
 
206
    enddo
 
207
  enddo
 
208
 
 
209
 
 
210
!     FLUX DE MASSE SUR LES FACETTES DE BORD
 
211
 
 
212
  do ifac = 1, nfabor
 
213
    ii = ifabor(ifac)
 
214
    ! Components U, V, W
 
215
    do isou = 1, 3
 
216
      pfac = inc*coefaq(isou,ifac)
 
217
 
 
218
      ! coefbv is a matrix
 
219
      do jsou = 1, 3
 
220
        pfac = pfac + romb(ifac)*coefbv(isou,jsou,ifac)*vel(jsou,ii)
 
221
      enddo
 
222
 
 
223
      flumab(ifac) = flumab(ifac) + pfac*surfbo(isou,ifac)
 
224
 
 
225
    enddo
 
226
  enddo
 
227
endif
 
228
 
 
229
 
 
230
!===============================================================================
 
231
! 4.  CALCUL DU FLUX DE MASSE AVEC TECHNIQUE DE RECONSTRUCTION
 
232
!        SI LE MAILLAGE EST NON ORTHOGONAL
 
233
!===============================================================================
 
234
 
 
235
if( nswrgu.gt.1 ) then
 
236
 
 
237
  allocate(grdqdm(3,3,ncelet))
 
238
 
 
239
 
 
240
!     CALCUL DU GRADIENT SUIVANT de QDM
 
241
!     =================================
 
242
  ! gradient vectoriel la periodicite est deja traitee
 
243
  ilved = .true.
 
244
 
 
245
  call grdvec &
 
246
  !==========
 
247
( ivar   , imrgra , inc    , iccocg , nswrgu , imligu ,          &
 
248
  iwarnu , nfecra , epsrgu , climgu , extrau ,                   &
 
249
  ilved  ,                                                       &
 
250
  qdm    , coefaq , coefbv ,                                     &
 
251
  grdqdm )
 
252
 
 
253
 
 
254
! ---> FLUX DE MASSE SUR LES FACETTES FLUIDES
 
255
 
 
256
  do ifac = 1, nfac
 
257
 
 
258
    ii = ifacel(1,ifac)
 
259
    jj = ifacel(2,ifac)
 
260
 
 
261
    pnd = pond(ifac)
 
262
 
 
263
    dofx = dofij(1,ifac)
 
264
    dofy = dofij(2,ifac)
 
265
    dofz = dofij(3,ifac)
 
266
 
 
267
! Termes suivant U, V, W
 
268
    do isou = 1, 3
 
269
      flumas(ifac) = flumas(ifac)                                   &
 
270
! Terme non reconstruit
 
271
         +( pnd*qdm(isou,ii) +(1.d0-pnd)*qdm(isou,jj)             &
 
272
!  --->
 
273
!  --->     ->    -->      ->
 
274
! (Grad(rho U ) . OFij ) . Sij
 
275
           +0.5d0*( grdqdm(isou,1,ii) +grdqdm(isou,1,jj) )*dofx   &
 
276
           +0.5d0*( grdqdm(isou,2,ii) +grdqdm(isou,2,jj) )*dofy   &
 
277
           +0.5d0*( grdqdm(isou,3,ii) +grdqdm(isou,3,jj) )*dofz   &
 
278
           )*surfac(isou,ifac)
 
279
    enddo
 
280
 
 
281
  enddo
 
282
 
 
283
! ---> FLUX DE MASSE SUR LES FACETTES DE BORD
 
284
  do ifac = 1, nfabor
 
285
 
 
286
    ii = ifabor(ifac)
 
287
    diipbx = diipb(1,ifac)
 
288
    diipby = diipb(2,ifac)
 
289
    diipbz = diipb(3,ifac)
 
290
 
 
291
! SUIVANT U, V, W
 
292
    do isou = 1, 3
 
293
 
 
294
      pfac = inc*coefaq(isou,ifac)
 
295
 
 
296
      ! coefu is a matrix
 
297
      do jsou = 1, 3
 
298
 
 
299
        pip =  romb(ifac)*vel(jsou,ii)                &
 
300
              +grdqdm(jsou,1,ii)*diipbx              &
 
301
              +grdqdm(jsou,2,ii)*diipby              &
 
302
              +grdqdm(jsou,3,ii)*diipbz
 
303
 
 
304
        pfac = pfac +coefbv(isou,jsou,ifac)*pip
 
305
      enddo
 
306
 
 
307
     flumab(ifac) = flumab(ifac) +pfac*surfbo(isou,ifac)
 
308
    enddo
 
309
 
 
310
  enddo
 
311
 
 
312
! DESALOCATION
 
313
deallocate(grdqdm)
 
314
deallocate(qdm, coefaq)
 
315
 
 
316
endif
 
317
 
 
318
!===============================================================================
 
319
! 6.  POUR S'ASSURER DE LA NULLITE DU FLUX DE MASSE AUX LIMITES
 
320
!       SYMETRIES PAROIS COUPLEES
 
321
!===============================================================================
 
322
 
 
323
if(iflmb0.eq.1) then
 
324
! FORCAGE DE FLUMAB a 0 pour la vitesse'
 
325
  do ifac = 1, nfabor
 
326
    if(isympa(ifac).eq.0) then
 
327
      flumab(ifac) = 0.d0
 
328
    endif
 
329
  enddo
 
330
endif
 
331
 
 
332
!--------
 
333
! FORMATS
 
334
!--------
 
335
 
 
336
#if defined(_CS_LANG_FR)
 
337
 
 
338
 1000 format('INIMAV APPELE AVEC INIT =',I10)
 
339
 
 
340
#else
 
341
 
 
342
 1000 format('INIMAV CALLED WITH INIT =',I10)
 
343
 
 
344
#endif
 
345
 
 
346
!----
 
347
! FIN
 
348
!----
 
349
 
 
350
return
 
351
 
 
352
end subroutine