← Back to branch summary

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

« back to all changes in this revision

Viewing changes to users/rayt/usray3.f90

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2011-11-01 17:43:32 UTC
  • mto: (6.1.7 sid)
  • mto: This revision was merged to the branch mainline in revision 11.
  • Revision ID: package-import@ubuntu.com-20111101174332-tl4vk45no0x3emc3
Tags: upstream-2.1.0
Import upstream version 2.1.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/rayt/usray3.f90
2
2
 
3
3
!VERS
4
4
 
5
 
 
6
 
!     This file is part of the Code_Saturne Kernel, element of the
7
 
!     Code_Saturne CFD tool.
8
 
 
9
 
!     Copyright (C) 1998-2009 EDF S.A., France
10
 
 
11
 
!     contact: saturne-support@edf.fr
12
 
 
13
 
!     The Code_Saturne Kernel is free software; you can redistribute it
14
 
!     and/or modify it under the terms of the GNU General Public License
15
 
!     as published by the Free Software Foundation; either version 2 of
16
 
!     the License, or (at your option) any later version.
17
 
 
18
 
!     The Code_Saturne Kernel is distributed in the hope that it will be
19
 
!     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
20
 
!     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
21
 
!     GNU General Public License for more details.
22
 
 
23
 
!     You should have received a copy of the GNU General Public License
24
 
!     along with the Code_Saturne Kernel; if not, write to the
25
 
!     Free Software Foundation, Inc.,
26
 
!     51 Franklin St, Fifth Floor,
27
 
!     Boston, MA  02110-1301  USA
 
5
! This file is part of Code_Saturne, a general-purpose CFD tool.
 
6
!
 
7
! Copyright (C) 1998-2011 EDF S.A.
 
8
!
 
9
! This program is free software; you can redistribute it and/or modify it under
 
10
! the terms of the GNU General Public License as published by the Free Software
 
11
! Foundation; either version 2 of the License, or (at your option) any later
 
12
! version.
 
13
!
 
14
! This program is distributed in the hope that it will be useful, but WITHOUT
 
15
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 
16
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 
17
! details.
 
18
!
 
19
! You should have received a copy of the GNU General Public License along with
 
20
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
 
21
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
28
22
 
29
23
!-------------------------------------------------------------------------------
30
24
 
31
25
subroutine usray3 &
32
26
!================
33
27
 
34
 
 ( idbia0 , idbra0 ,                                              &
35
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
36
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
37
 
   nvar   , nscal  , iphas  , iappel ,                            &
38
 
   nideve , nrdeve , nituse , nrtuse ,                            &
39
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml , itypfb ,          &
40
 
   ipnfac , nodfac , ipnfbr , nodfbr , izfrdp ,                   &
41
 
   idevel , ituser , ia     ,                                     &
42
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
 
28
 ( nvar   , nscal  , iappel ,                                     &
 
29
   itypfb ,                                                       &
 
30
   izfrdp ,                                                       &
43
31
   dt     , rtp    , rtpa   , propce , propfa , propfb ,          &
44
 
   ck     , w1     , w2     , w3     , w4     , w5     ,  w6    , &
45
 
   rdevel , rtuser ,                                              &
46
 
   ra     )
 
32
   ck     )
47
33
 
48
34
!===============================================================================
49
 
! FONCTION :
50
 
! ----------
51
 
 
52
 
!   SOUS-PROGRAMME DU MODULE DE RAYONNEMENT :
53
 
!   -----------------------------------------
54
 
 
55
 
 
56
 
!     Coefficient d'absorption
57
 
!    -------------------------
58
 
 
59
 
!      Il est indispensable de renseigner la valeur du coefficient
60
 
!        d'absorption du fluide CK.
61
 
 
62
 
!      Pour un milieu transparent, le coefficient doit etre
63
 
!        fixe a 0.D0.
64
 
 
65
 
!  DANS LE CAS DU MODELE P-1 ON VERIFIE QUE LA LONGUEUR OPTIQUE
66
 
!    DU MILIEU EST AU MINIMUM DE L'ORDRE DE L'UNITE
67
 
 
68
 
!      ATTENTION :
69
 
!      =========
70
 
!        Pour les physiques particulieres (Combustion, charbon...)
71
 
 
72
 
!        il est INTERDIT de fournir le COEFFICIENT D'ABSORPTION ici.
73
 
!               ========
74
 
 
75
 
!        Voir le sous-programme PPCABS
76
 
 
 
35
! Purpose:
 
36
! --------
 
37
 
 
38
! Absorption coefficient for radiative module
 
39
! ----------------------
 
40
 
 
41
! It is necessary to define the value of the fluid's absorption
 
42
! coefficient Ck.
 
43
 
 
44
! For a transparent medium, the coefficient should be set to 0.d0.
 
45
 
 
46
! In the case of the P-1 model, we check that the optical length is at
 
47
! least of the order of 1.
 
48
 
 
49
! Caution:
 
50
! ========
 
51
!   For specific physics (Combustion, coal, ...),
 
52
 
 
53
!   it is Forbidden to define the absorption coefficient here.
 
54
!         =========
 
55
 
 
56
!   See subroutine ppcabs.
77
57
 
78
58
!-------------------------------------------------------------------------------
79
 
!ARGU                             ARGUMENTS
 
59
! Arguments
80
60
!__________________.____._____.________________________________________________.
81
61
! name             !type!mode ! role                                           !
82
62
!__________________!____!_____!________________________________________________!
83
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
84
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
85
 
! ndim             ! i  ! <-- ! spatial dimension                              !
86
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
87
 
! ncel             ! i  ! <-- ! number of cells                                !
88
 
! nfac             ! i  ! <-- ! number of interior faces                       !
89
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
90
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
91
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
92
 
! nnod             ! i  ! <-- ! number of vertices                             !
93
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
94
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
95
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
96
63
! nvar             ! i  ! <-- ! total number of variables                      !
97
64
! nscal            ! i  ! <-- ! total number of scalars                        !
98
 
! iphas            ! i  ! <-- ! phase number                                   !
99
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
100
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
101
 
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
102
 
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
103
 
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
104
 
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
105
 
! iprfml           ! ia ! <-- ! property numbers per family                    !
106
 
!  (nfml, nprfml)  !    !     !                                                !
107
65
! itypfb           ! ia ! <-- ! boundary face types                            !
108
 
!  (nfabor, nphas) !    !     !                                                !
109
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
110
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
111
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
112
 
! nodfbr(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
113
 
! izfrdp(nfabor    ! te ! <-- ! numero de zone pour les faces de bord          !
114
 
! idevel(nideve)   ! ia ! <-> ! integer work array for temporary development   !
115
 
! ituser(nituse)   ! ia ! <-> ! user-reserved integer work array               !
116
 
! ia(*)            ! ia ! --- ! main integer work array                        !
117
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
118
 
!  (ndim, ncelet)  !    !     !                                                !
119
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
120
 
!  (ndim, nfac)    !    !     !                                                !
121
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
122
 
!  (ndim, nfabor)  !    !     !                                                !
123
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
124
 
!  (ndim, nfac)    !    !     !                                                !
125
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
126
 
!  (ndim, nfabor)  !    !     !                                                !
127
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
128
 
!  (ndim, nnod)    !    !     !                                                !
129
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
 
66
! izfrdp(nfabor    ! ia ! <-- ! zone number for boundary faces                 !
130
67
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
131
68
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
132
69
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
133
70
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
134
71
! propfa(nfac, *)  ! ra ! <-- ! physical properties at interior face centers   !
135
72
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers   !
136
 
! ck (ncelet)      ! tr ! --> ! coefficient d'absorption du milieu             !
137
 
!                  !    !     ! (nul si transparent)                           !
138
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary development      !
139
 
! w1...6(ncelet    ! tr ! --- ! tableau de travail                             !
140
 
! rtuser(nrtuse)   ! ra ! <-> ! user-reserved real work array                  !
141
 
! ra(*)            ! ra ! --- ! main real work array                           !
 
73
! ck(ncelet)       ! ra ! --> ! medium's absorption coefficient                !
 
74
!                  !    !     ! (zero if transparent)                          !
142
75
!__________________!____!_____!________________________________________________!
143
76
 
144
77
!     Type: i (integer), r (real), s (string), a (array), l (logical),
146
79
!     mode: <-- input, --> output, <-> modifies data, --- work array
147
80
!===============================================================================
148
81
 
 
82
!===============================================================================
 
83
! Module files
 
84
!===============================================================================
 
85
 
 
86
use paramx
 
87
use numvar
 
88
use entsor
 
89
use optcal
 
90
use cstphy
 
91
use cstnum
 
92
use parall
 
93
use period
 
94
use ppppar
 
95
use ppthch
 
96
use cpincl
 
97
use ppincl
 
98
use radiat
 
99
use ihmpre
 
100
use mesh
 
101
 
 
102
!===============================================================================
 
103
 
149
104
implicit none
150
105
 
151
 
!===============================================================================
152
 
! Common blocks
153
 
!===============================================================================
154
 
 
155
 
include "paramx.h"
156
 
include "numvar.h"
157
 
include "entsor.h"
158
 
include "optcal.h"
159
 
include "cstphy.h"
160
 
include "cstnum.h"
161
 
include "pointe.h"
162
 
include "parall.h"
163
 
include "period.h"
164
 
include "ppppar.h"
165
 
include "ppthch.h"
166
 
include "cpincl.h"
167
 
include "ppincl.h"
168
 
include "radiat.h"
169
 
include "ihmpre.h"
170
 
 
171
 
!===============================================================================
172
 
 
173
106
! Arguments
174
107
 
175
 
integer          idbia0 , idbra0
176
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
177
 
integer          nfml   , nprfml
178
 
integer          nnod   , lndfac , lndfbr , ncelbr
179
 
integer          nvar   , nscal  , iphas  , iappel
180
 
integer          nideve , nrdeve , nituse , nrtuse
181
 
 
182
 
integer          ifacel(2,nfac) , ifabor(nfabor)
183
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
184
 
integer          iprfml(nfml,nprfml) , itypfb(nfabor)
185
 
integer          ipnfac(nfac+1), nodfac(lndfac)
186
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr),izfrdp(nfabor)
187
 
integer          idevel(nideve), ituser(nituse), ia(*)
188
 
 
189
 
double precision xyzcen(ndim,ncelet)
190
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
191
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
192
 
double precision xyznod(ndim,nnod), volume(ncelet)
 
108
integer          nvar   , nscal  , iappel
 
109
 
 
110
integer          itypfb(nfabor)
 
111
integer          izfrdp(nfabor)
 
112
 
193
113
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
194
114
double precision propce(ncelet,*)
195
115
double precision propfa(nfac,*), propfb(nfabor,*)
196
116
 
197
 
double precision w1(ncelet), w2(ncelet), w3(ncelet)
198
 
double precision w4(ncelet), w5(ncelet), w6(ncelet)
199
117
 
200
118
double precision ck(ncelet)
201
119
 
202
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
203
 
 
204
 
 
205
120
! Local variables
206
121
 
207
 
integer          idebia , idebra , iel, ifac, iok
 
122
integer          iel, ifac, iok
208
123
double precision vv, sf, xlc, xkmin, pp
209
124
 
210
125
!===============================================================================
213
128
!===============================================================================
214
129
 
215
130
!===============================================================================
216
 
! 0.  CE TEST PERMET A L'UTILISATEUR D'ETRE CERTAIN QUE C'EST
217
 
!       SA VERSION DU SOUS PROGRAMME QUI EST UTILISEE
218
 
!       ET NON CELLE DE LA BIBLIOTHEQUE
 
131
! 0.  This test allows the user to be certain his version if the subroutine
 
132
!     is being used, and not that from the library.
219
133
!===============================================================================
220
134
 
221
135
if (iihmpr.eq.1) then
244
158
 
245
159
 
246
160
!===============================================================================
247
 
! 0 - GESTION MEMOIRE
 
161
! 0 - Memory management
248
162
!===============================================================================
249
163
 
250
 
idebia = idbia0
251
 
idebra = idbra0
252
164
 
253
 
! Indicateur d'arret (pour savoir si des faces ont ete oubliees)
 
165
! Stop flag (to determine if faces were forgotten)
254
166
iok = 0
255
167
 
256
168
!===============================================================================
257
 
!   COEFFICIENT D'ABSORPTION DU MILIEU (m-1)
258
 
 
259
 
!     DANS LE CAS DES PHYSIQUES PARTICULIERES (COMBUSTION GAZ/CHARBON/ELEC/FIOUL)
260
 
 
261
 
!         CK NE DOIT PAS ETRE FOURNI
262
 
       !==========
263
 
!         (il est determine automatiquement, eventuellement a partir
264
 
!          du fichier parametrique)
265
 
 
266
 
!     DANS LES AUTRES CAS
267
 
!         CK DOIT ETRE COMPLETE (IL EST NUL PAR DEFAUT)
268
 
 
269
 
 
 
169
! Absorption coefficient of the medium (m-1)
 
170
 
 
171
! In the case of specific physics (gas/coal/fuel combustion, elec)
 
172
 
 
173
! Ck must not be defined here
 
174
!============
 
175
! (it is determined automatically, possibly from the parametric file)
 
176
 
 
177
! In other cases, Ck must be defined (it is zero by default)
270
178
!===============================================================================
271
179
 
272
 
 
273
 
 
274
 
 
275
 
if( ippmod(iphpar).le.1 ) then
 
180
if (ippmod(iphpar).le.1) then
276
181
 
277
182
  do iel = 1, ncel
278
183
    ck(iel) = 0.d0
279
184
  enddo
280
185
 
281
 
!--> MODELE P1 : Controle standard des valeurs du coefficient
282
 
!                d'absorption. Ce coefficient doit assurer une
283
 
!                longueur optique au minimum de l'ordre de l'unite.
 
186
  !--> P1 model: standard control of absorption coefficient values.
 
187
  !              this coefficient must ensure an optical length
 
188
  !              at least of the order of 1.
284
189
 
285
190
  if (iirayo.eq.2) then
286
191
    sf = 0.d0
287
192
    vv = 0.d0
288
193
 
289
 
!         Calcul de la longueur caract�ristique du domaine de calcul
 
194
    ! Compute characteristic length of calculation domain
290
195
 
291
196
    do ifac = 1,nfabor
292
 
      sf = sf + sqrt(                                             &
293
 
                surfbo(1,ifac)**2 +                               &
294
 
                surfbo(2,ifac)**2 +                               &
295
 
                surfbo(3,ifac)**2 )
 
197
      sf = sf + sqrt(surfbo(1,ifac)**2 + surfbo(2,ifac)**2 + surfbo(3,ifac)**2)
296
198
    enddo
297
199
    if (irangp.ge.0) then
298
200
      call parsom(sf)
309
211
 
310
212
    xlc = 3.6d0 * vv / sf
311
213
 
312
 
!         Clipping pour la variable CK
 
214
    ! Clipping for variable CK
313
215
 
314
216
    xkmin = 1.d0 / xlc
315
217
 
321
223
      endif
322
224
    enddo
323
225
 
324
 
!     Arret en fin de pas de temps si epaisseur optique trop grande
325
 
!       (ISTPP1 = 1 permet d'arreter proprement a la fin du pas de temps
326
 
!        courant).
 
226
    ! Stop at the end of the time step if the optical thickness is too big
 
227
    ! (istpp1 = 1 allows stopping cleanly at the end of the current time step).
327
228
    pp = xnp1mx/100.0d0
328
229
    if (dble(iok).gt.pp*dble(ncel)) then
329
 
      write(nfecra,3000) xkmin, dble(iok)/dble(ncel)*100.d0,      &
330
 
                         xnp1mx
 
230
      write(nfecra,3000) xkmin, dble(iok)/dble(ncel)*100.d0, xnp1mx
331
231
      istpp1 = 1
332
 
!            CALL CSEXIT (1)
333
 
       !==========
 
232
      ! call csexit(1)
 
233
      !==========
334
234
    endif
335
235
  endif
336
236
 
337
237
endif
338
 
! -------
339
 
! FORMATS
340
 
! -------
341
 
 
342
 
 3000 format(                                                           &
 
238
 
 
239
! -------
 
240
! Formats
 
241
! -------
 
242
 
 
243
 3000 format(                                                     &
343
244
'@                                                            ',/,&
344
245
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
345
246
'@                                                            ',/,&
346
 
'@ @@ ATTENTION : RAYONNEMENT APPROXIMATION P-1 (USRAY3)      ',/,&
347
 
'@    =========                                               ',/,&
348
 
'@                                                            ',/,&
349
 
'@    LA LONGUEUR OPTIQUE DU MILIEU SEMI-TRANSPARENT          ',/,&
350
 
'@      DOIT AU MOINS ETRE DE L''ORDRE DE L''UNITE POUR ETRE  ',/,&
351
 
'@      DANS LE DOMAINE D''APPLICATION DE L''APPROXIMATION P-1',/,&
352
 
'@    CELA NE SEMBLE PAS ETRE LE CAS ICI.                     ',/,&
353
 
'@                                                            ',/,&
354
 
'@    LE COEFFICIENT D''ABSORPTION MINIMUM POUR ASSURER CETTE ',/,&
355
 
'@      LONGUEUR OPTIQUE EST XKMIN = ',E10.4                   ,/,&
356
 
'@    CETTE VALEUR N''EST PAS ATTEINTE POUR ', E10.4,'%       ',/,&
357
 
'@      DES CELLULES DU MAILLAGE.                             ',/,&
358
 
'@    LE POURCENTAGE DE CELLULES DU MAILLAGE POUR LESQUELLES  ',/,&
359
 
'@      ON ADMET QUE CETTE CONDITION SOIT VIOLEE EST IMPOSE   ',/,&
360
 
'@      PAR DEFAUT OU DANS USINI1 A XNP1MX = ', E10.4,'%      ',/,&
361
 
'@                                                            ',/,&
362
 
'@    Le calcul est interrompu.                               ',/,&
363
 
'@                                                            ',/,&
364
 
'@    Verifier les valeurs du coefficient d''absorption CK    ',/,&
365
 
'@      dans PPCABS, USRAY3 ou Fichier thermochimie.          ',/,&
 
247
'@ @@ WARNING:    P1 radiation approximation (usray3)         ',/,&
 
248
'@    ========                                                ',/,&
 
249
'@                                                            ',/,&
 
250
'@    The optical length of the semi-transparent medium       ',/,&
 
251
'@      must be at least of the order of one to be in the     ',/,&
 
252
'@      domain of validity of the P-1 approximation.          ',/,&
 
253
'@    This does not seem to be the case here.                 ',/,&
 
254
'@                                                            ',/,&
 
255
'@    The minimum absorption coefficient to ensure this       ',/,&
 
256
'@      optical length is XKmin = ', e10.4                     ,/,&
 
257
'@    This value is not reached for ', e10.4,'%               ',/,&
 
258
'@      of the meshe''s cells.                                ',/,&
 
259
'@    The percentage of mesh cells for which we allow this    ',/,&
 
260
'@      condition not to be rspected is set by default in     ',/,&
 
261
'@      usini1 to xnp1mx = ', e10.4,'%                        ',/,&
 
262
'@                                                            ',/,&
 
263
'@    The calculation is interrupted.                         ',/,&
 
264
'@                                                            ',/,&
 
265
'@    Check the values of the absorption coefficient Ck       ',/,&
 
266
'@      in ppcabs, usray3 or the thermochemistry data file.   ',/,&
366
267
'@                                                            ',/,&
367
268
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
368
269
'@                                                            ',/)
369
270
 
370
 
 
371
271
end subroutine