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

« back to all changes in this revision

Viewing changes to users/elec/uselcl.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:
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 uselcl &
32
26
!================
33
27
 
34
 
 ( idbia0 , idbra0 ,                                              &
35
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
36
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
37
 
   nvar   , nscal  , nphas  ,                                     &
38
 
   nideve , nrdeve , nituse , nrtuse ,                            &
39
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml , maxelt , lstelt , &
40
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
 
28
 ( nvar   , nscal  ,                                              &
41
29
   icodcl , itrifb , itypfb , izfppp ,                            &
42
 
   idevel , ituser , ia     ,                                     &
43
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
44
30
   dt     , rtp    , rtpa   , propce , propfa , propfb ,          &
45
 
   coefa  , coefb  , rcodcl ,                                     &
46
 
   w1     , w2     , w3     , w4     , w5     , w6     , coefu  , &
47
 
   rdevel , rtuser , ra     )
 
31
   coefa  , coefb  , rcodcl )
48
32
 
49
33
!===============================================================================
50
34
! Purpose
93
77
!__________________.____._____.________________________________________________.
94
78
! name             !type!mode ! role                                           !
95
79
!__________________!____!_____!________________________________________________!
96
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
97
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
98
 
! ndim             ! i  ! <-- ! spatial dimension                              !
99
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
100
 
! ncel             ! i  ! <-- ! number of cells                                !
101
 
! nfac             ! i  ! <-- ! number of interior faces                       !
102
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
103
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
104
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
105
 
! nnod             ! i  ! <-- ! number of vertices                             !
106
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
107
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
108
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
109
80
! nvar             ! i  ! <-- ! total number of variables                      !
110
81
! nscal            ! i  ! <-- ! total number of scalars                        !
111
 
! nphas            ! i  ! <-- ! number of phases                               !
112
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
113
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
114
 
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
115
 
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
116
 
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
117
 
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
118
 
! iprfml           ! ia ! <-- ! property numbers per family                    !
119
 
!  (nfml, nprfml)  !    !     !                                                !
120
 
! maxelt           !  e ! <-- ! max number of cells and faces (int/boundary)   !
121
 
! lstelt(maxelt)   ! ia ! --- ! work array                                     !
122
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
123
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
124
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
125
 
! nodfac(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
126
82
! icodcl           ! ia ! --> ! boundary condition code                        !
127
83
!  (nfabor, nvar)  !    !     ! = 1  -> Dirichlet                              !
128
84
!                  !    !     ! = 2  -> flux density                           !
132
88
!                  !    !     ! = 9  -> free inlet/outlet (velocity)           !
133
89
!                  !    !     !         inflowing possibly blocked             !
134
90
! itrifb(nfabor    ! ia ! <-- ! indirection for boundary faces ordering)       !
135
 
!  (nfabor, nphas) !    !     !                                                !
136
91
! itypfb           ! ia ! --> ! boundary face types                            !
137
 
!  (nfabor, nphas) !    !     !                                                !
138
 
! idevel(nideve)   ! ia ! <-- ! integer work array for temporary developpement !
139
 
! ituser(nituse    ! ia ! <-- ! user-reserved integer work array               !
140
 
! ia(*)            ! ia ! --- ! main integer work array                        !
141
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
142
 
!  (ndim, ncelet)  !    !     !                                                !
143
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
144
 
!  (ndim, nfac)    !    !     !                                                !
145
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
146
 
!  (ndim, nfavor)  !    !     !                                                !
147
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
148
 
!  (ndim, nfac)    !    !     !                                                !
149
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
150
 
!  (ndim, nfabor)  !    !     !                                                !
151
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
152
 
!  (ndim, nnod)    !    !     !                                                !
153
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
154
92
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
155
93
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
156
94
!  (ncelet, *)     !    !     !  (at current and preceding time steps)         !
169
107
!                  !    !     ! for velocities           ( vistl+visct)*gradu  !
170
108
!                  !    !     ! for pressure                         dt*gradp  !
171
109
!                  !    !     ! for scalars    cp*(viscls+visct/sigmas)*gradt  !
172
 
! w1,2,3,4,5,6     ! ra ! --- ! work arrays                                    !
173
 
!  (ncelet)        !    !     !  (computation of pressure gradient)            !
174
 
! coefu            ! ra ! --- ! work array                                     !
175
 
!  (nfabor, 3)     !    !     !  (computation of pressure gradient)            !
176
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary developpement    !
177
 
! rtuser(nituse    ! ra ! <-- ! user-reserved real work array                  !
178
 
! ra(*)            ! ra ! --- ! main real work array                           !
179
110
!__________________!____!_____!________________________________________________!
180
111
 
181
112
!     Type: i (integer), r (real), s (string), a (array), l (logical),
183
114
!     mode: <-- input, --> output, <-> modifies data, --- work array
184
115
!===============================================================================
185
116
 
 
117
!===============================================================================
 
118
! Module files
 
119
!===============================================================================
 
120
 
 
121
use paramx
 
122
use numvar
 
123
use optcal
 
124
use cstphy
 
125
use cstnum
 
126
use entsor
 
127
use ppppar
 
128
use ppthch
 
129
use ppincl
 
130
use elincl
 
131
use mesh
 
132
 
 
133
!===============================================================================
 
134
 
186
135
implicit none
187
136
 
188
 
!===============================================================================
189
 
!    Common blocks
190
 
!===============================================================================
191
 
 
192
 
include "paramx.h"
193
 
include "numvar.h"
194
 
include "optcal.h"
195
 
include "cstphy.h"
196
 
include "cstnum.h"
197
 
include "entsor.h"
198
 
include "ppppar.h"
199
 
include "ppthch.h"
200
 
include "ppincl.h"
201
 
include "elincl.h"
202
 
 
203
 
!===============================================================================
204
 
 
205
137
! Arguments
206
138
 
207
 
integer          idbia0 , idbra0
208
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
209
 
integer          nfml   , nprfml
210
 
integer          nnod   , lndfac , lndfbr , ncelbr
211
 
integer          nvar   , nscal  , nphas
212
 
integer          nideve , nrdeve , nituse , nrtuse
 
139
integer          nvar   , nscal
213
140
 
214
 
integer          ifacel(2,nfac) , ifabor(nfabor)
215
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
216
 
integer          iprfml(nfml,nprfml)
217
 
integer          maxelt, lstelt(maxelt)
218
 
integer          ipnfac(nfac+1), nodfac(lndfac)
219
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
220
141
integer          icodcl(nfabor,nvar)
221
 
integer          itrifb(nfabor,nphas), itypfb(nfabor,nphas)
 
142
integer          itrifb(nfabor), itypfb(nfabor)
222
143
integer          izfppp(nfabor)
223
 
integer          idevel(nideve), ituser(nituse), ia(*)
224
144
 
225
 
double precision xyzcen(ndim,ncelet)
226
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
227
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
228
 
double precision xyznod(ndim,nnod), volume(ncelet)
229
145
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
230
146
double precision propce(ncelet,*)
231
147
double precision propfa(nfac,*), propfb(nfabor,*)
232
148
double precision coefa(nfabor,*), coefb(nfabor,*)
233
149
double precision rcodcl(nfabor,nvar,3)
234
 
double precision w1(ncelet),w2(ncelet),w3(ncelet)
235
 
double precision w4(ncelet),w5(ncelet),w6(ncelet)
236
 
double precision coefu(nfabor,ndim)
237
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
238
150
 
239
151
! Local variables
240
152
 
241
 
integer          idebia, idebra
242
 
integer          ifac, ii, iphas , iel
 
153
integer          ifac, ii, iel
243
154
integer          idim
244
155
integer          izone,iesp
245
156
integer          ilelt, nlelt
249
160
double precision xkent, xeent
250
161
double precision z1   , z2
251
162
 
 
163
integer, allocatable, dimension(:) :: lstelt
 
164
 
252
165
!===============================================================================
253
166
 
254
167
! TEST_TO_REMOVE_FOR_USE_OF_SUBROUTINE_START
285
198
! 1.  Initialization
286
199
!===============================================================================
287
200
 
288
 
idebia = idbia0
289
 
idebra = idbra0
 
201
! Allocate a temporary array for boundary faces selection
 
202
allocate(lstelt(nfabor))
 
203
 
290
204
 
291
205
d2s3 = 2.d0/3.d0
292
206
 
298
212
!         Set the boundary condition for each face
299
213
!===============================================================================
300
214
 
301
 
iphas = 1
302
 
 
303
215
! --- For boundary faces of color 1 assign an inlet for all phases
304
216
!     ============================================================
305
217
!        and assign a cathode for "electric" variables.
312
224
 
313
225
  ifac = lstelt(ilelt)
314
226
 
315
 
  itypfb(ifac,iphas) = ientre
 
227
  itypfb(ifac) = ientre
316
228
 
317
229
!      - Zone Number (from 1 to n)
318
230
  izone = 1
320
232
!      - Zone localization for a given face
321
233
  izfppp(ifac) = izone
322
234
 
323
 
  rcodcl(ifac,iu(iphas),1) = 0.d0
324
 
  rcodcl(ifac,iv(iphas),1) = 0.d0
325
 
  rcodcl(ifac,iw(iphas),1) = 0.d0
 
235
  rcodcl(ifac,iu,1) = 0.d0
 
236
  rcodcl(ifac,iv,1) = 0.d0
 
237
  rcodcl(ifac,iw,1) = 0.d0
326
238
 
327
239
!         Turbulence
328
240
 
329
241
!     (ITYTUR est un indicateur qui vaut ITURB/10)
330
 
  if (itytur(iphas).eq.2 .or. itytur(iphas).eq.3                  &
331
 
       .or. iturb(iphas).eq.50 .or. iturb(iphas).eq.60) then
 
242
  if (itytur.eq.2 .or. itytur.eq.3                  &
 
243
       .or. iturb.eq.50 .or. iturb.eq.60            &
 
244
       .or. iturb.eq.70) then
332
245
 
333
 
    uref2 = rcodcl(ifac,iu(iphas),1)**2                           &
334
 
           +rcodcl(ifac,iv(iphas),1)**2                           &
335
 
           +rcodcl(ifac,iw(iphas),1)**2
 
246
    uref2 = rcodcl(ifac,iu,1)**2                           &
 
247
           +rcodcl(ifac,iv,1)**2                           &
 
248
           +rcodcl(ifac,iw,1)**2
336
249
    uref2 = max(uref2,1.d-12)
337
250
 
338
251
!   Turbulence example computed using equations valid for a pipe.
358
271
!     standard laws for a circular pipe
359
272
!     (their initialization is not needed here but is good practice).
360
273
 
361
 
    rhomoy = propfb(ifac,ipprob(irom(iphas)))
 
274
    rhomoy = propfb(ifac,ipprob(irom))
362
275
    ustar2 = 0.d0
363
276
    xkent  = epzero
364
277
    xeent  = epzero
365
278
 
366
279
    call keendb                                                   &
367
280
    !==========
368
 
     ( uref2, dhy, rhomoy, viscl0(iphas), cmu, xkappa,            &
 
281
     ( uref2, dhy, rhomoy, viscl0, cmu, xkappa,            &
369
282
       ustar2, xkent, xeent )
370
283
 
371
 
    if (itytur(iphas).eq.2) then
372
 
 
373
 
      rcodcl(ifac,ik(iphas),1)  = xkent
374
 
      rcodcl(ifac,iep(iphas),1) = xeent
375
 
 
376
 
    elseif(itytur(iphas).eq.3) then
377
 
 
378
 
      rcodcl(ifac,ir11(iphas),1) = d2s3*xkent
379
 
      rcodcl(ifac,ir22(iphas),1) = d2s3*xkent
380
 
      rcodcl(ifac,ir33(iphas),1) = d2s3*xkent
381
 
      rcodcl(ifac,ir12(iphas),1) = 0.d0
382
 
      rcodcl(ifac,ir13(iphas),1) = 0.d0
383
 
      rcodcl(ifac,ir23(iphas),1) = 0.d0
384
 
      rcodcl(ifac,iep(iphas),1)  = xeent
385
 
 
386
 
    elseif (iturb(iphas).eq.50) then
387
 
 
388
 
      rcodcl(ifac,ik(iphas),1)   = xkent
389
 
      rcodcl(ifac,iep(iphas),1)  = xeent
390
 
      rcodcl(ifac,iphi(iphas),1) = d2s3
391
 
      rcodcl(ifac,ifb(iphas),1)  = 0.d0
392
 
 
393
 
    elseif (iturb(iphas).eq.60) then
394
 
 
395
 
      rcodcl(ifac,ik(iphas),1)   = xkent
396
 
      rcodcl(ifac,iomg(iphas),1) = xeent/cmu/xkent
 
284
    if (itytur.eq.2) then
 
285
 
 
286
      rcodcl(ifac,ik,1)  = xkent
 
287
      rcodcl(ifac,iep,1) = xeent
 
288
 
 
289
    elseif(itytur.eq.3) then
 
290
 
 
291
      rcodcl(ifac,ir11,1) = d2s3*xkent
 
292
      rcodcl(ifac,ir22,1) = d2s3*xkent
 
293
      rcodcl(ifac,ir33,1) = d2s3*xkent
 
294
      rcodcl(ifac,ir12,1) = 0.d0
 
295
      rcodcl(ifac,ir13,1) = 0.d0
 
296
      rcodcl(ifac,ir23,1) = 0.d0
 
297
      rcodcl(ifac,iep,1)  = xeent
 
298
 
 
299
    elseif (iturb.eq.50) then
 
300
 
 
301
      rcodcl(ifac,ik,1)   = xkent
 
302
      rcodcl(ifac,iep,1)  = xeent
 
303
      rcodcl(ifac,iphi,1) = d2s3
 
304
      rcodcl(ifac,ifb,1)  = 0.d0
 
305
 
 
306
    elseif (iturb.eq.60) then
 
307
 
 
308
      rcodcl(ifac,ik,1)   = xkent
 
309
      rcodcl(ifac,iomg,1) = xeent/cmu/xkent
 
310
 
 
311
    elseif (iturb.eq.70) then
 
312
 
 
313
      rcodcl(ifac,inusa,1) = cmu*xkent**2/xeent
397
314
 
398
315
    endif
399
316
 
472
389
 
473
390
  ifac = lstelt(ilelt)
474
391
!
475
 
  itypfb(ifac,iphas)   = isolib
 
392
  itypfb(ifac)   = isolib
476
393
 
477
394
!      - Zone Number (from 1 to n)
478
395
  izone = 2
540
457
 
541
458
  ifac = lstelt(ilelt)
542
459
!
543
 
  itypfb(ifac,iphas)   = isolib
 
460
  itypfb(ifac)   = isolib
544
461
 
545
462
!      - Zone number (from 1 to n)
546
463
  izone = 3
611
528
 
612
529
  ifac = lstelt(ilelt)
613
530
!
614
 
  itypfb(ifac,iphas)   = iparoi
 
531
  itypfb(ifac)   = iparoi
615
532
 
616
533
!      - Zone number (from 1 to n)
617
534
  izone = 4
686
603
 
687
604
  ifac = lstelt(ilelt)
688
605
 
689
 
  itypfb(ifac,iphas)   = iparoi
 
606
  itypfb(ifac)   = iparoi
690
607
 
691
608
!      - Zone number (from 1 to n)
692
609
  izone = 5
763
680
 
764
681
!          SYMETRIES
765
682
 
766
 
  itypfb(ifac,iphas)   = isymet
 
683
  itypfb(ifac)   = isymet
767
684
 
768
685
!      - Zone number (from 1 to n)
769
686
  izone = 6
793
710
! FIN
794
711
!----
795
712
 
 
713
! Deallocate the temporary array
 
714
deallocate(lstelt)
 
715
 
796
716
return
797
717
end subroutine