← Back to branch summary

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

« back to all changes in this revision

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