← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/rayt/raypun.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/rayt/raypun.f90
1
1
!-------------------------------------------------------------------------------
2
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
 
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.
25
20
 
26
21
!-------------------------------------------------------------------------------
27
22
 
28
23
subroutine raypun &
29
24
!================
30
25
 
31
 
 ( idbia0 , idbra0 ,                                              &
32
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
33
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
34
 
   nvar   , nscal  , nphas  , iphas  ,                            &
35
 
   nideve , nrdeve , nituse , nrtuse ,                            &
36
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml , itypfb ,          &
37
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
38
 
   idevel , ituser , ia     ,                                     &
39
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
 
26
 ( nvar   , nscal  ,                                              &
 
27
   itypfb ,                                                       &
40
28
   dt     , rtp    , rtpa   , propce , propfa , propfb ,          &
41
29
   coefa  , coefb  ,                                              &
42
30
   cofrua , cofrub ,                                              &
43
31
   flurds , flurdb ,                                              &
44
 
   dtr    , viscf  , viscb  ,                                     &
45
 
   dam    , xam    ,                                              &
46
 
   drtp   , smbrs  , rovsdt ,                                     &
 
32
   viscf  , viscb  ,                                              &
 
33
   smbrs  , rovsdt ,                                              &
47
34
   theta4 , thetaa , sa     ,                                     &
48
35
   qx     , qy     , qz     ,                                     &
49
36
   qincid , eps    , tparoi ,                                     &
50
 
   w1     , w2     , w3     , w4     , w5     ,                   &
51
 
   w6     , w7     , w8     , w9     , ckmel  ,                   &
52
 
   rdevel , rtuser , ra     )
 
37
   ckmel  )
53
38
 
54
39
!===============================================================================
55
40
! FONCTION :
66
51
!__________________.____._____.________________________________________________.
67
52
! name             !type!mode ! role                                           !
68
53
!__________________!____!_____!________________________________________________!
69
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
70
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
71
 
! ndim             ! i  ! <-- ! spatial dimension                              !
72
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
73
 
! ncel             ! i  ! <-- ! number of cells                                !
74
 
! nfac             ! i  ! <-- ! number of interior faces                       !
75
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
76
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
77
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
78
 
! nnod             ! i  ! <-- ! number of vertices                             !
79
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
80
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
81
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
82
54
! nvar             ! i  ! <-- ! total number of variables                      !
83
55
! nscal            ! i  ! <-- ! total number of scalars                        !
84
 
! nphas            ! i  ! <-- ! number of phases                               !
85
 
! iphas            ! i  ! --> ! phase number                                   !
86
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
87
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
88
 
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
89
 
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
90
 
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
91
 
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
92
 
! iprfml           ! ia ! <-- ! property numbers per family                    !
93
 
!  (nfml, nprfml)  !    !     !                                                !
94
56
! itypfb           ! ia ! <-- ! boundary face types                            !
95
 
!  (nfabor, nphas) !    !     !                                                !
96
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
97
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
98
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
99
 
! nodfbr(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
100
 
! idevel(nideve)   ! ia ! <-> ! integer work array for temporary development   !
101
 
! ituser(nituse)   ! ia ! <-> ! user-reserved integer work array               !
102
 
! ia(*)            ! ia ! --- ! main integer work array                        !
103
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
104
 
!  (ndim, ncelet)  !    !     !                                                !
105
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
106
 
!  (ndim, nfac)    !    !     !                                                !
107
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
108
 
!  (ndim, nfabor)  !    !     !                                                !
109
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
110
 
!  (ndim, nfac)    !    !     !                                                !
111
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
112
 
!  (ndim, nfabor)  !    !     !                                                !
113
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
114
 
!  (ndim, nnod)    !    !     !                                                !
115
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
116
57
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
117
58
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
118
59
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
125
66
!(nfabor)          !    !     !    faces de bord pour la luminance             !
126
67
! flurds,flurdb    ! tr ! --- ! pseudo flux de masse (faces internes           !
127
68
!(nfac)(nfabor)    !    !     !    et faces de bord )                          !
128
 
! dtr(ncelet)      ! tr ! --- ! dt*cdtvar                                      !
129
69
! viscf(nfac)      ! tr ! --- ! visc*surface/dist aux faces internes           !
130
70
! viscb(nfabor     ! tr ! --- ! visc*surface/dist aux faces de bord            !
131
 
! dam(ncelet       ! tr ! --- ! tableau de travail pour matrice                !
132
 
! xam(nfac,*)      ! tr ! --- ! tableau de travail pour matrice                !
133
 
! drtp(ncelet      ! tr ! --- ! tableau de travail pour increment              !
134
71
! smbrs(ncelet     ! tr ! --- ! tableau de travail pour sec mem                !
135
72
! rovsdt(ncelet    ! tr ! --- ! tableau de travail pour terme instat           !
136
73
! theta4(ncelet    ! tr ! --- ! pseudo temperature radiative                   !
141
78
! qincid(nfabor    ! tr ! --> ! densite de flux radiatif aux bords             !
142
79
! eps (nfabor)     ! tr ! <-- ! emissivite des facettes de bord                !
143
80
! tparoi(nfabor    ! tr ! <-- ! temperature de paroi en kelvin                 !
144
 
! w1...9(ncelet    ! tr ! --- ! tableau de travail                             !
145
81
! ckmel(ncelet)    ! tr ! <-- ! coeff d'absorption du melange                  !
146
82
!                  !    !     !   gaz-particules de charbon                    !
147
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary development      !
148
 
! rtuser(nrtuse)   ! ra ! <-> ! user-reserved real work array                  !
149
 
! ra(*)            ! ra ! --- ! main real work array                           !
150
83
!__________________!____!_____!________________________________________________!
151
84
 
152
85
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
153
86
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
154
87
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
155
88
!            --- tableau de travail
156
 
!-------------------------------------------------------------------------------
 
89
!===============================================================================
 
90
 
 
91
!===============================================================================
 
92
! Module files
 
93
!===============================================================================
 
94
 
 
95
use paramx
 
96
use numvar
 
97
use entsor
 
98
use optcal
 
99
use cstphy
 
100
use cstnum
 
101
use ppppar
 
102
use ppthch
 
103
use cpincl
 
104
use ppincl
 
105
use radiat
 
106
use parall
 
107
use period
 
108
use mesh
 
109
 
157
110
!===============================================================================
158
111
 
159
112
implicit none
160
113
 
161
 
!===============================================================================
162
 
! Common blocks
163
 
!===============================================================================
164
 
 
165
 
include "paramx.h"
166
 
include "numvar.h"
167
 
include "entsor.h"
168
 
include "optcal.h"
169
 
include "cstphy.h"
170
 
include "cstnum.h"
171
 
include "pointe.h"
172
 
include "ppppar.h"
173
 
include "ppthch.h"
174
 
include "cpincl.h"
175
 
include "ppincl.h"
176
 
include "radiat.h"
177
 
include "parall.h"
178
 
include "period.h"
179
 
 
180
 
 
181
 
!===============================================================================
182
 
 
183
114
! Arguments
184
115
 
185
 
integer          idbia0 , idbra0
186
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
187
 
integer          nfml   , nprfml
188
 
integer          nnod   , lndfac , lndfbr , ncelbr
189
 
integer          nvar   , nscal  , nphas  , iphas
190
 
integer          nideve , nrdeve , nituse , nrtuse
191
 
 
192
 
integer          ifacel(2,nfac) , ifabor(nfabor)
193
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
194
 
integer          iprfml(nfml,nprfml) , itypfb(nfabor,nphas)
195
 
integer          ipnfac(nfac+1), nodfac(lndfac)
196
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
197
 
integer          idevel(nideve), ituser(nituse)
198
 
integer          ia(*)
199
 
 
200
 
double precision xyzcen(ndim,ncelet)
201
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
202
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
203
 
double precision xyznod(ndim,nnod), volume(ncelet)
 
116
integer          nvar   , nscal
 
117
 
 
118
integer          itypfb(nfabor)
 
119
 
204
120
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
205
121
double precision propce(ncelet,*)
206
122
double precision propfa(nfac,*), propfb(nfabor,*)
209
125
double precision cofrua(nfabor), cofrub(nfabor)
210
126
double precision flurds(nfac), flurdb(nfabor)
211
127
 
212
 
double precision dtr(ncelet)
213
128
double precision viscf(nfac), viscb(nfabor)
214
 
double precision dam(ncelet), xam(nfac,2)
215
 
double precision drtp(ncelet), smbrs(ncelet)
 
129
double precision smbrs(ncelet)
216
130
double precision rovsdt(ncelet)
217
131
 
218
132
double precision theta4(ncelet), thetaa(ncelet)
220
134
double precision qx(ncelet), qy(ncelet), qz(ncelet)
221
135
double precision qincid(nfabor), tparoi(nfabor), eps(nfabor)
222
136
 
223
 
double precision w1(ncelet), w2(ncelet), w3(ncelet)
224
 
double precision w4(ncelet), w5(ncelet), w6(ncelet)
225
 
double precision w7(ncelet), w8(ncelet), w9(ncelet)
226
137
double precision ckmel(ncelet)
227
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
228
 
 
229
138
 
230
139
! Local variables
231
140
 
232
141
character*80     cnom
233
142
 
234
 
integer          idebia, idebra
235
143
integer          ifac  , iel
236
144
integer          iconv1, idiff1, ndirc1, ireso1
237
145
integer          nitmap, nswrsp, nswrgp, iwarnp
239
147
integer          ncymap, nitmgp
240
148
integer          inum
241
149
integer          idtva0, ivar0
242
 
integer          inc, iccocg, iphydp
243
 
integer          idimte , itenso
 
150
integer          inc, iccocg
244
151
double precision epsrgp, blencp, climgp, epsilp, extrap, epsrsp
245
152
double precision aa, aaa, aaaa, relaxp, thetap
246
153
 
 
154
double precision rvoid(1)
 
155
 
 
156
double precision, allocatable, dimension(:,:) :: grad
247
157
 
248
158
!===============================================================================
249
159
 
251
161
! 0. GESTION MEMOIRE
252
162
!===============================================================================
253
163
 
254
 
idebia = idbia0
255
 
idebra = idbra0
256
 
 
257
164
!===============================================================================
258
165
! 1. PARAMETRAGE DU SOLVEUR ET INITIALISATION
259
166
!===============================================================================
302
209
!--> Remise a zero des tableaux avant resolution
303
210
 
304
211
do iel = 1,ncel
305
 
  drtp(iel)   = zero
306
212
  theta4(iel) = zero
307
213
  thetaa(iel) = zero
308
214
enddo
325
231
 
326
232
call viscfa                                                       &
327
233
!==========
328
 
   ( idebia , idebra ,                                            &
329
 
     ndim   , ncelet , ncel   , nfac   , nfabor , nfml   ,        &
330
 
     nprfml , nnod   , lndfac , lndfbr , ncelbr ,                 &
331
 
     nideve , nrdeve , nituse , nrtuse , imvisf ,                 &
332
 
     ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                 &
333
 
     ipnfac , nodfac , ipnfbr , nodfbr ,                          &
334
 
     idevel , ituser , ia     ,                                   &
335
 
     xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod ,        &
336
 
     volume , ckmel  , viscf  , viscb  ,                          &
337
 
     rdevel , rtuser , ra     )
 
234
   ( imvisf ,                                                     &
 
235
     ckmel  , viscf  , viscb  )
338
236
 
339
237
!===============================================================================
340
238
! 3.  RESOLUTION
351
249
 
352
250
call codits                                                       &
353
251
!==========
354
 
 ( idebia , idebra ,                                              &
355
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml , nprfml ,   &
356
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
357
 
   nvar   , nscal  , nphas  ,                                     &
358
 
   nideve , nrdeve , nituse , nrtuse ,                            &
 
252
 ( nvar   , nscal  ,                                              &
359
253
   idtva0 , ivar0  , iconv1 , idiff1 , ireso1 , ndirc1 , nitmap , &
360
254
   imrgra , nswrsp , nswrgp , imligp , ircflp ,                   &
361
255
   ischcp , isstpp , iescap ,                                     &
362
256
   imgr1  , ncymap , nitmgp , inum   , iwarnp ,                   &
363
257
   blencp , epsilp , epsrsp , epsrgp , climgp , extrap ,          &
364
258
   relaxp , thetap ,                                              &
365
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
366
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
367
 
   idevel , ituser , ia     ,                                     &
368
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod ,          &
369
 
   volume ,                                                       &
370
259
   thetaa , thetaa , cofrua , cofrub , cofrua , cofrub ,          &
371
260
   flurds , flurdb ,                                              &
372
261
   viscf  , viscb  , viscf  , viscb  ,                            &
373
262
   rovsdt , smbrs  , theta4 ,                                     &
374
 
   dam    , xam    , drtp   ,                                     &
375
 
   w1     , w2     , w3     , w4     , w5     ,                   &
376
 
   w6     , w7     , w8     , w9     ,                            &
377
 
   rdevel , rtuser , ra     )
 
263
   rvoid  )
378
264
 
379
265
!===============================================================================
380
266
! 4. Vecteur densite de flux radiatif
381
267
!===============================================================================
382
268
 
 
269
! Allocate a temporary array for gradient computation
 
270
allocate(grad(ncelet,3))
 
271
 
383
272
!    En periodique et parallele, echange avant calcul du gradient
384
 
 
385
 
!    Parallele
386
 
if (irangp.ge.0) then
387
 
  call parcom (theta4)
388
 
  !==========
389
 
endif
390
 
 
391
 
!    Periodique
392
 
if (iperio.eq.1) then
393
 
  idimte = 0
394
 
  itenso = 0
395
 
  call percom                                                     &
396
 
  !==========
397
 
  ( idimte , itenso ,                                             &
398
 
    theta4 , theta4 , theta4 ,                                    &
399
 
    theta4 , theta4 , theta4 ,                                    &
400
 
    theta4 , theta4 , theta4)
 
273
if (irangp.ge.0.or.iperio.eq.1) then
 
274
  call synsca(theta4)
 
275
  !==========
401
276
endif
402
277
 
403
278
!     Calcul de la densite du flux radiatif QX, QY, QZ
411
286
extrap  = 0.d0
412
287
nswrgp  = 100
413
288
ivar0   = 0
414
 
iphydp  = 0
415
289
 
416
290
call grdcel                                                       &
417
291
!==========
418
 
   ( idebia , idebra ,                                            &
419
 
     ndim   , ncelet , ncel   , nfac   , nfabor , nfml,           &
420
 
     nprfml ,                                                     &
421
 
     nnod   , lndfac , lndfbr , ncelbr , nphas  ,                 &
422
 
     nideve , nrdeve , nituse , nrtuse ,                          &
423
 
     ivar0  , imrgra , inc    , iccocg , nswrgp , imligp,         &
424
 
     iphydp ,                                                     &
 
292
   ( ivar0  , imrgra , inc    , iccocg , nswrgp , imligp,         &
425
293
     iwarnp , nfecra , epsrgp , climgp , extrap ,                 &
426
 
     ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                 &
427
 
     ipnfac , nodfac , ipnfbr , nodfbr ,                          &
428
 
     idevel , ituser , ia     ,                                   &
429
 
     xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod ,        &
430
 
     volume ,                                                     &
431
 
     w7     , w7     , w7     ,                                   &
432
294
     theta4 , cofrua , cofrub ,                                   &
433
 
     w1     , w2     , w3     ,                                   &
434
 
     w4     , w5     , w6     ,                                   &
435
 
     rdevel , rtuser , ra     )
 
295
     grad   )
436
296
 
437
297
aa = - stephn * 4.d0 / 3.d0
438
298
 
439
299
do iel = 1,ncel
440
300
  aaa = aa * ckmel(iel)
441
 
  qx(iel) = w1(iel) * aaa
442
 
  qy(iel) = w2(iel) * aaa
443
 
  qz(iel) = w3(iel) * aaa
 
301
  qx(iel) = grad(iel,1) * aaa
 
302
  qy(iel) = grad(iel,2) * aaa
 
303
  qz(iel) = grad(iel,3) * aaa
444
304
enddo
445
305
 
 
306
! Free memory
 
307
deallocate(grad)
 
308
 
446
309
!===============================================================================
447
310
! 5. Terme Source Radiatif d'absorption et densite de flux incident
448
311
!===============================================================================
460
323
 
461
324
  iel = ifabor(ifac)
462
325
 
463
 
  if (itypfb(ifac,iphas).eq.iparoi .or.                           &
464
 
      itypfb(ifac,iphas).eq.iparug ) then
 
326
  if (itypfb(ifac).eq.iparoi .or.                           &
 
327
      itypfb(ifac).eq.iparug ) then
465
328
 
466
329
!--> Premiere version plus chere et legerement plus precise
467
330
 
468
331
    aaaa = tparoi(ifac)**4
469
332
 
470
 
    aaa  = 1.5d0 * ra(idistb-1+ifac) / ckmel(iel)                 &
 
333
    aaa  = 1.5d0 * distb(ifac) / ckmel(iel)                 &
471
334
           * ( 2.d0 /(2.d0-eps(ifac)) -1.d0 )
472
335
    aa   = ( aaa * aaaa + theta4(iel) ) / (1.d0 + aaa)
473
336
 
485
348
               + ( qx(iel) * surfbo(1,ifac) +                     &
486
349
                   qy(iel) * surfbo(2,ifac) +                     &
487
350
                   qz(iel) * surfbo(3,ifac) ) /                   &
488
 
                   (0.5d0 * ra(isrfbn-1+ifac) )
 
351
                   (0.5d0 * surfbn(ifac) )
489
352
  endif
490
353
 
491
354
enddo