← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/alge/matrxv.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/matrxv.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 matrxv &
 
24
!================
 
25
 
 
26
 ( ncelet , ncel   , nfac   , nfabor ,                            &
 
27
   iconvp , idiffp , ndircp , isym   , nfecra ,                   &
 
28
   thetap ,                                                       &
 
29
   ifacel , ifabor ,                                              &
 
30
   coefbu , cofbfu , fimp   , flumas , flumab , viscf  , viscb  , &
 
31
   da     , xa     )
 
32
 
 
33
!===============================================================================
 
34
! FONCTION :
 
35
! ----------
 
36
 
 
37
! CONSTRUCTION DE LA MATRICE DE CONVECTION UPWIND/DIFFUSION/TS DE LA VITESSE
 
38
 
 
39
!     IL EST INTERDIT DE MODIFIER FIMP   ICI
 
40
 
 
41
 
 
42
!-------------------------------------------------------------------------------
 
43
!ARGU                             ARGUMENTS
 
44
!__________________.____._____.________________________________________________.
 
45
! name             !type!mode ! role                                           !
 
46
!__________________!____!_____!________________________________________________!
 
47
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
 
48
! ncel             ! i  ! <-- ! number of cells                                !
 
49
! nfac             ! i  ! <-- ! number of interior faces                       !
 
50
! nfabor           ! i  ! <-- ! number of boundary faces                       !
 
51
! iconvp           ! e  ! <-- ! indicateur = 1 convection, 0 sinon             !
 
52
! idiffp           ! e  ! <-- ! indicateur = 1 diffusion , 0 sinon             !
 
53
! ndircp           ! e  ! <-- ! indicateur = 0 si decalage diagonale           !
 
54
! isym             ! e  ! <-- ! indicateur = 1 matrice symetrique              !
 
55
!                  !    !     !              2 matrice non symetrique          !
 
56
! thetap           ! r  ! <-- ! coefficient de ponderation pour le             !
 
57
!                  !    !     ! theta-schema (on ne l'utilise pour le          !
 
58
!                  !    !     ! moment que pour u,v,w et les scalaire          !
 
59
!                  !    !     ! - thetap = 0.5 correspond a un schema          !
 
60
!                  !    !     !   totalement centre en temps (mixage           !
 
61
!                  !    !     !   entre crank-nicolson et adams-               !
 
62
!                  !    !     !   bashforth)                                   !
 
63
! ifacel(2,nfac)   ! te ! <-- ! no des elts voisins d'une face intern          !
 
64
! ifabor(nfabor)   ! te ! <-- ! no de l'elt voisin d'une face de bord          !
 
65
! coefbu(3,3,nfabor! tr ! <-- ! tab b des cl pour la vitesse                   !
 
66
! cofbfu(3,3,nfabor! tr ! <-- ! tab b des cl de flux pour la vitesse           !
 
67
! flumas(nfac)     ! tr ! <-- ! flux de masse aux faces internes               !
 
68
! flumab(nfabor    ! tr ! <-- ! flux de masse aux faces de bord                !
 
69
! viscf(nfac)      ! tr ! <-- ! visc*surface/dist aux faces internes           !
 
70
! viscb(nfabor     ! tr ! <-- ! visc*surface/dist aux faces de bord            !
 
71
! da (3,3,ncelet)  ! tr ! --> ! partie diagonale de la matrice                 !
 
72
! xa (*,nfac)      ! tr ! --> ! extra interleaved diagonale de la matrice      !
 
73
!__________________!____!_____!________________________________________________!
 
74
 
 
75
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
 
76
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
 
77
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
 
78
!            --- tableau de travail
 
79
!===============================================================================
 
80
 
 
81
!===============================================================================
 
82
! Module files
 
83
!===============================================================================
 
84
 
 
85
use parall
 
86
 
 
87
!===============================================================================
 
88
 
 
89
implicit none
 
90
 
 
91
! Arguments
 
92
 
 
93
integer          ncelet , ncel   , nfac   , nfabor
 
94
integer          iconvp , idiffp , ndircp , isym
 
95
integer          nfecra
 
96
double precision thetap
 
97
 
 
98
integer          ifacel(2,nfac), ifabor(nfabor)
 
99
double precision coefbu(3,3,nfabor), fimp(3,3,ncelet), cofbfu(3,3,nfabor)
 
100
double precision flumas(nfac), flumab(nfabor)
 
101
double precision viscf(nfac), viscb(nfabor)
 
102
double precision da(3,3,ncelet),xa(isym,nfac)
 
103
 
 
104
! Local variables
 
105
 
 
106
integer          ifac,ii,jj,iel, isou, jsou
 
107
double precision flui,fluj,epsi
 
108
 
 
109
!===============================================================================
 
110
 
 
111
!===============================================================================
 
112
! 1. INITIALISATION
 
113
!===============================================================================
 
114
 
 
115
if(isym.ne.1.and.isym.ne.2) then
 
116
  write(nfecra,1000) isym
 
117
  call csexit (1)
 
118
endif
 
119
 
 
120
epsi = 1.d-7
 
121
 
 
122
do iel = 1, ncel
 
123
  do isou = 1, 3
 
124
    do jsou = 1, 3
 
125
      da(isou,jsou,iel) = fimp(isou,jsou,iel)
 
126
    enddo
 
127
  enddo
 
128
enddo
 
129
if(ncelet.gt.ncel) then
 
130
  do iel = ncel+1, ncelet
 
131
    do isou = 1, 3
 
132
      do jsou = 1, 3
 
133
        da(isou,jsou,iel) = 0.d0
 
134
      enddo
 
135
    enddo
 
136
  enddo
 
137
endif
 
138
 
 
139
if(isym.eq.2) then
 
140
  do ifac = 1, nfac
 
141
    xa(1,ifac) = 0.d0
 
142
    xa(2,ifac) = 0.d0
 
143
  enddo
 
144
else
 
145
  do ifac = 1, nfac
 
146
    xa(1,ifac) = 0.d0
 
147
  enddo
 
148
endif
 
149
 
 
150
!===============================================================================
 
151
! 2.    CALCUL DES TERMES EXTRADIAGONAUX
 
152
!===============================================================================
 
153
 
 
154
if(isym.eq.2) then
 
155
 
 
156
  do ifac = 1,nfac
 
157
    flui = 0.5d0*( flumas(ifac) -abs(flumas(ifac)) )
 
158
    fluj =-0.5d0*( flumas(ifac) +abs(flumas(ifac)) )
 
159
    xa(1,ifac) = thetap*(iconvp*flui -idiffp*viscf(ifac))
 
160
    xa(2,ifac) = thetap*(iconvp*fluj -idiffp*viscf(ifac))
 
161
  enddo
 
162
 
 
163
else
 
164
 
 
165
  do ifac = 1,nfac
 
166
    flui = 0.5d0*( flumas(ifac) -abs(flumas(ifac)) )
 
167
    xa(1,ifac) = thetap*(iconvp*flui -idiffp*viscf(ifac))
 
168
  enddo
 
169
 
 
170
endif
 
171
 
 
172
!===============================================================================
 
173
! 3.     CONTRIBUTION DES TERMES X-TRADIAGONAUX A LA DIAGONALE
 
174
!===============================================================================
 
175
 
 
176
if(isym.eq.2) then
 
177
 
 
178
  do ifac = 1,nfac
 
179
    ii = ifacel(1,ifac)
 
180
    jj = ifacel(2,ifac)
 
181
    do isou = 1, 3
 
182
      da(isou,isou,ii) = da(isou,isou,ii) -xa(2,ifac)
 
183
      da(isou,isou,jj) = da(isou,isou,jj) -xa(1,ifac)
 
184
    enddo
 
185
  enddo
 
186
 
 
187
 
 
188
else
 
189
 
 
190
  do ifac = 1,nfac
 
191
    ii = ifacel(1,ifac)
 
192
    jj = ifacel(2,ifac)
 
193
    do isou = 1, 3
 
194
      da(isou,isou,ii) = da(isou,isou,ii) -xa(1,ifac)
 
195
      da(isou,isou,jj) = da(isou,isou,jj) -xa(1,ifac)
 
196
    enddo
 
197
  enddo
 
198
 
 
199
endif
 
200
 
 
201
!===============================================================================
 
202
! 4.     CONTRIBUTION DES FACETTES DE BORDS A LA DIAGONALE
 
203
!===============================================================================
 
204
 
 
205
do ifac = 1,nfabor
 
206
  ii = ifabor(ifac)
 
207
  flui = 0.5d0*( flumab(ifac) -abs(flumab(ifac)) )
 
208
  fluj =-0.5d0*( flumab(ifac) +abs(flumab(ifac)) )
 
209
  do isou = 1, 3
 
210
    do jsou = 1, 3
 
211
      if(isou.eq.jsou) then
 
212
        da(isou,jsou,ii) = da(isou,jsou,ii) + thetap*(                    &
 
213
                       iconvp*(-fluj + flui*coefbu(isou,jsou,ifac) )      &
 
214
                      +idiffp*viscb(ifac)*(1.d0-cofbfu(isou,jsou,ifac))   &
 
215
                             )
 
216
      else
 
217
        da(isou,jsou,ii) = da(isou,jsou,ii) + thetap*(                    &
 
218
                       iconvp*( flui*coefbu(isou,jsou,ifac) )             &
 
219
                      +idiffp*viscb(ifac)*(-cofbfu(isou,jsou,ifac))       &
 
220
                             )
 
221
      endif
 
222
    enddo
 
223
  enddo
 
224
enddo
 
225
 
 
226
 
 
227
!===============================================================================
 
228
! 5.  NON PRESENCE DE PTS DIRICHLET --> LEGER RENFORCEMENT DE LA
 
229
!     DIAGONALE POUR DECALER LE SPECTRE DES VALEURS PROPRES
 
230
!===============================================================================
 
231
!     (si IDIRCL=0, on a force NDIRCP a valoir au moins 1 pour ne pas
 
232
!      decaler la diagonale)
 
233
 
 
234
if ( ndircp.le.0 ) then
 
235
  do iel=1,ncel
 
236
    do isou = 1, 3
 
237
      da(isou,isou,iel) = (1.d0+epsi)*da(isou,isou,iel)
 
238
    enddo
 
239
  enddo
 
240
endif
 
241
 
 
242
!--------
 
243
! FORMATS
 
244
!--------
 
245
 
 
246
#if defined(_CS_LANG_FR)
 
247
 
 
248
 1000 format(                                                           &
 
249
'@                                                            ',/,&
 
250
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
 
251
'@                                                            ',/,&
 
252
'@ @@ ATTENTION : ARRET DANS matrxv                           ',/,&
 
253
'@    =========                                               ',/,&
 
254
'@     APPEL DE matrxv              AVEC ISYM   = ',I10        ,/,&
 
255
'@                                                            ',/,&
 
256
'@  Le calcul ne peut pas etre execute.                       ',/,&
 
257
'@                                                            ',/,&
 
258
'@  Contacter l''assistance.                                  ',/,&
 
259
'@                                                            ',/,&
 
260
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
 
261
'@                                                            ',/)
 
262
 
 
263
#else
 
264
 
 
265
 1000 format(                                                           &
 
266
'@'                                                            ,/,&
 
267
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
 
268
'@'                                                            ,/,&
 
269
'@ @@ WARNING: ABORT IN matrxv'                                ,/,&
 
270
'@    ========'                                                ,/,&
 
271
'@     matrxv CALLED                WITH ISYM   = ',I10        ,/,&
 
272
'@'                                                            ,/,&
 
273
'@  The calculation will not be run.'                          ,/,&
 
274
'@'                                                            ,/,&
 
275
'@  Contact support.'                                          ,/,&
 
276
'@'                                                            ,/,&
 
277
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
 
278
'@'                                                            ,/)
 
279
 
 
280
#endif
 
281
 
 
282
!----
 
283
! FIN
 
284
!----
 
285
 
 
286
return
 
287
 
 
288
end subroutine