← Back to branch summary

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

« back to all changes in this revision

Viewing changes to users/base/uskpdc.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/base/uskpdc.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 uskpdc &
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
 
   ncepdp , iphas  , iappel ,                                     &
40
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml , maxelt , lstelt , &
41
 
   ipnfac , nodfac , ipnfbr , nodfbr , icepdc ,                   &
42
 
   idevel , ituser , ia     ,                                     &
43
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
 
28
 ( nvar   , nscal  ,                                              &
 
29
   ncepdp , iappel ,                                              &
 
30
   icepdc , izcpdc ,                                              &
44
31
   dt     , rtpa   , rtp    , propce , propfa , propfb ,          &
45
 
   coefa  , coefb  , ckupdc ,                                     &
46
 
   rdevel , rtuser , ra     )
 
32
   coefa  , coefb  , ckupdc )
47
33
 
48
34
!===============================================================================
49
35
! FONCTION :
50
36
! ----------
51
37
 
52
38
!                    PERTES DE CHARGE (PDC)
53
 
!                       POUR LA PHASE IPHAS
54
39
 
55
40
! IAPPEL = 1 :
56
41
!             CALCUL DU NOMBRE DE CELLULES OU L'ON IMPOSE UNE PDC
103
88
!__________________.____._____.________________________________________________.
104
89
! name             !type!mode ! role                                           !
105
90
!__________________!____!_____!________________________________________________!
106
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
107
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
108
 
! ndim             ! i  ! <-- ! spatial dimension                              !
109
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
110
 
! ncel             ! i  ! <-- ! number of cells                                !
111
 
! nfac             ! i  ! <-- ! number of interior faces                       !
112
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
113
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
114
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
115
 
! nnod             ! i  ! <-- ! number of vertices                             !
116
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
117
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
118
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
119
91
! nvar             ! i  ! <-- ! total number of variables                      !
120
92
! nscal            ! i  ! <-- ! total number of scalars                        !
121
 
! nphas            ! i  ! <-- ! number of phases                               !
122
93
! ncepdp           ! i  ! <-- ! number of cells with head loss                 !
123
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
124
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
125
94
! iappel           ! e  ! <-- ! indique les donnes a renvoyer                  !
126
 
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
127
 
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
128
 
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
129
 
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
130
 
! iprfml           ! ia ! <-- ! property numbers per family                    !
131
 
!  (nfml, nprfml)  !    !     !                                                !
132
 
! maxelt           ! i  ! <-- ! max number of cells and faces (int/boundary)   !
133
 
! lstelt(maxelt)   ! ia ! --- ! work array                                     !
134
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
135
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
136
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
137
 
! nodfbr(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
138
95
! icepdc(ncepdp    ! te ! <-- ! numero des ncepdp cellules avec pdc            !
139
 
! idevel(nideve)   ! ia ! <-> ! integer work array for temporary development   !
140
 
! ituser(nituse)   ! ia ! <-> ! user-reserved integer work array               !
141
 
! ia(*)            ! ia ! --- ! main integer work array                        !
142
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
143
 
!  (ndim, ncelet)  !    !     !                                                !
144
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
145
 
!  (ndim, nfac)    !    !     !                                                !
146
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
147
 
!  (ndim, nfabor)  !    !     !                                                !
148
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
149
 
!  (ndim, nfac)    !    !     !                                                !
150
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
151
 
!  (ndim, nfabor)  !    !     !                                                !
152
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
153
 
!  (ndim, nnod)    !    !     !                                                !
154
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
 
96
! izcpdc(ncelet)   ! ia ! <-- ! cells zone for head loss definition            !
155
97
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
156
98
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
157
99
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
162
104
!  (nfabor, *)     !    !     !                                                !
163
105
! ckupdc           ! tr ! <-- ! tableau de travail pour pdc                    !
164
106
!  (ncepdp,6)      !    !     !                                                !
165
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary development      !
166
 
! rtuser(nrtuse)   ! ra ! <-> ! user-reserved real work array                  !
167
 
! ra(*)            ! ra ! --- ! main real work array                           !
168
107
!__________________!____!_____!________________________________________________!
169
108
 
170
109
!     Type: i (integer), r (real), s (string), a (array), l (logical),
172
111
!     mode: <-- input, --> output, <-> modifies data, --- work array
173
112
!===============================================================================
174
113
 
 
114
!===============================================================================
 
115
! Module files
 
116
!===============================================================================
 
117
 
 
118
use paramx
 
119
use numvar
 
120
use optcal
 
121
use cstnum
 
122
use parall
 
123
use period
 
124
use mesh
 
125
 
 
126
!===============================================================================
 
127
 
175
128
implicit none
176
129
 
177
 
!===============================================================================
178
 
! Common blocks
179
 
!===============================================================================
180
 
 
181
 
include "paramx.h"
182
 
include "pointe.h"
183
 
include "numvar.h"
184
 
include "optcal.h"
185
 
include "cstnum.h"
186
 
include "parall.h"
187
 
include "period.h"
188
 
 
189
 
!===============================================================================
190
 
 
191
130
! Arguments
192
131
 
193
 
integer          idbia0 , idbra0
194
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
195
 
integer          nfml   , nprfml
196
 
integer          nnod   , lndfac , lndfbr , ncelbr
197
 
integer          nvar   , nscal  , nphas
 
132
integer          nvar   , nscal
198
133
integer          ncepdp
199
 
integer          nideve , nrdeve , nituse , nrtuse
200
134
integer          iappel
201
135
 
202
 
integer          ifacel(2,nfac) , ifabor(nfabor)
203
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
204
 
integer          iprfml(nfml,nprfml)
205
 
integer          maxelt, lstelt(maxelt)
206
 
integer          ipnfac(nfac+1), nodfac(lndfac)
207
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
208
136
integer          icepdc(ncepdp)
209
 
integer          idevel(nideve), ituser(nituse), ia(*)
 
137
integer          izcpdc(ncel)
210
138
 
211
 
double precision xyzcen(ndim,ncelet)
212
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
213
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
214
 
double precision xyznod(ndim,nnod), volume(ncelet)
215
139
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
216
140
double precision propce(ncelet,*)
217
141
double precision propfa(nfac,*), propfb(nfabor,*)
218
142
double precision coefa(nfabor,*), coefb(nfabor,*)
219
143
double precision ckupdc(ncepdp,6)
220
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
221
144
 
222
145
! Local variables
223
146
 
224
 
integer          idebia, idebra
225
 
integer          iel, ielpdc, iphas, ikpdc
 
147
integer          iel, ielpdc, ikpdc
226
148
integer          ilelt, nlelt
 
149
integer          izone
227
150
integer          iutile
228
151
 
229
152
double precision alpha, cosalp, sinalp, vit, ck1, ck2
230
153
 
 
154
integer, allocatable, dimension(:) :: lstelt
 
155
 
231
156
!===============================================================================
232
157
 
233
158
! TEST_TO_REMOVE_FOR_USE_OF_SUBROUTINE_START
241
166
 
242
167
!===============================================================================
243
168
 
244
 
idebia = idbia0
245
 
idebra = idbra0
 
169
! Allocate a temporary array for cells selection
 
170
allocate(lstelt(ncel))
 
171
 
246
172
 
247
173
if(iappel.eq.1.or.iappel.eq.2) then
248
174
 
290
216
 
291
217
!       Ce test permet de desactiver l'exemple
292
218
  if(1.eq.0) then
293
 
    if(iphas.eq.1) then
294
 
      ielpdc = 0
295
 
 
296
 
      CALL GETCEL('X <= 6.0 and X >= 4.0 and Y >= 2.0 and'//      &
297
 
                  'Y <= 8.0',NLELT,LSTELT)
298
 
 
299
 
      do ilelt = 1, nlelt
300
 
        iel = lstelt(ilelt)
301
 
        ielpdc = ielpdc + 1
302
 
        if (iappel.eq.2) icepdc(ielpdc) = iel
303
 
      enddo
304
 
 
305
 
    else
306
 
      ielpdc = 0
307
 
    endif
 
219
 
 
220
    izone = 0
 
221
    ielpdc = 0
 
222
 
 
223
    call getcel('X <= 6.0 and X >= 4.0 and Y >= 2.0 and'//      &
 
224
         'Y <= 8.0',nlelt,lstelt)
 
225
 
 
226
    izone = izone + 1
 
227
 
 
228
    do ilelt = 1, nlelt
 
229
      iel = lstelt(ilelt)
 
230
      izcpdc(iel) = izone
 
231
      ielpdc = ielpdc + 1
 
232
      if (iappel.eq.2) icepdc(ielpdc) = iel
 
233
    enddo
 
234
 
308
235
  endif
309
236
 
310
237
 
364
291
    enddo
365
292
  enddo
366
293
 
367
 
  if(iphas.eq.1) then
368
 
 
369
 
! --- Tenseur diagonal
370
 
!   Exemple de pertes de charges dans la direction x
371
 
 
372
 
    iutile = 0
373
 
    if (iutile.eq.0) return
374
 
 
375
 
    do ielpdc = 1, ncepdp
376
 
      iel=icepdc(ielpdc)
377
 
      vit = sqrt(  rtpa(iel,iu(iphas))**2                         &
378
 
                 + rtpa(iel,iv(iphas))**2                         &
379
 
                 + rtpa(iel,iw(iphas))**2)
380
 
      ckupdc(ielpdc,1) = 10.d0*vit
381
 
      ckupdc(ielpdc,2) =  0.d0*vit
382
 
      ckupdc(ielpdc,3) =  0.d0*vit
383
 
    enddo
384
 
 
385
 
! --- Tenseur 3x3
386
 
!   Exemple de pertes de charges a ALPHA = 45 degres x,y
387
 
!      la direction x resiste par ck1 et y par ck2
388
 
!      ck2 nul represente des ailettes comme ceci :  ///////
389
 
!      dans le repere de calcul X Y
390
 
 
391
 
!                 Y|    /y
392
 
!                  |  /
393
 
!                  |/
394
 
!                  \--------------- X
395
 
!                   \ / ALPHA
396
 
!                    \
397
 
!                     \ x
398
 
 
399
 
    iutile = 0
400
 
    if (iutile.eq.0) return
401
 
 
402
 
    alpha  = pi/4.d0
403
 
    cosalp = cos(alpha)
404
 
    sinalp = sin(alpha)
405
 
    ck1 = 10.d0
406
 
    ck2 =  0.d0
407
 
 
408
 
    do ielpdc = 1, ncepdp
409
 
      iel=icepdc(ielpdc)
410
 
      vit = sqrt(  rtpa(iel,iu(iphas))**2                         &
411
 
                 + rtpa(iel,iv(iphas))**2                         &
412
 
                 + rtpa(iel,iw(iphas))**2)
413
 
      ckupdc(ielpdc,1) = (cosalp**2*ck1 + sinalp**2*ck2)*vit
414
 
      ckupdc(ielpdc,2) = (sinalp**2*ck1 + cosalp**2*ck2)*vit
415
 
      ckupdc(ielpdc,3) =  0.d0
416
 
      ckupdc(ielpdc,4) = cosalp*sinalp*(-ck1+ck2)*vit
417
 
      ckupdc(ielpdc,5) =  0.d0
418
 
      ckupdc(ielpdc,6) =  0.d0
419
 
    enddo
420
 
 
421
 
  endif
422
 
 
 
294
  ! --- Tenseur diagonal
 
295
  !   Exemple de pertes de charges dans la direction x
 
296
 
 
297
  iutile = 0
 
298
  if (iutile.eq.0) return
 
299
 
 
300
  do ielpdc = 1, ncepdp
 
301
    iel=icepdc(ielpdc)
 
302
    vit = sqrt(  rtpa(iel,iu)**2                         &
 
303
         + rtpa(iel,iv)**2                         &
 
304
         + rtpa(iel,iw)**2)
 
305
    ckupdc(ielpdc,1) = 10.d0*vit
 
306
    ckupdc(ielpdc,2) =  0.d0*vit
 
307
    ckupdc(ielpdc,3) =  0.d0*vit
 
308
  enddo
 
309
 
 
310
  ! --- Tenseur 3x3
 
311
  !   Exemple de pertes de charges a ALPHA = 45 degres x,y
 
312
  !      la direction x resiste par ck1 et y par ck2
 
313
  !      ck2 nul represente des ailettes comme ceci :  ///////
 
314
  !      dans le repere de calcul X Y
 
315
 
 
316
  !                 Y|    /y
 
317
  !                  |  /
 
318
  !                  |/
 
319
  !                  \--------------- X
 
320
  !                   \ / ALPHA
 
321
  !                    \
 
322
  !                     \ x
 
323
 
 
324
  iutile = 0
 
325
  if (iutile.eq.0) return
 
326
 
 
327
  alpha  = pi/4.d0
 
328
  cosalp = cos(alpha)
 
329
  sinalp = sin(alpha)
 
330
  ck1 = 10.d0
 
331
  ck2 =  0.d0
 
332
 
 
333
  do ielpdc = 1, ncepdp
 
334
    iel=icepdc(ielpdc)
 
335
    vit = sqrt(  rtpa(iel,iu)**2                         &
 
336
         + rtpa(iel,iv)**2                         &
 
337
         + rtpa(iel,iw)**2)
 
338
    ckupdc(ielpdc,1) = (cosalp**2*ck1 + sinalp**2*ck2)*vit
 
339
    ckupdc(ielpdc,2) = (sinalp**2*ck1 + cosalp**2*ck2)*vit
 
340
    ckupdc(ielpdc,3) =  0.d0
 
341
    ckupdc(ielpdc,4) = cosalp*sinalp*(-ck1+ck2)*vit
 
342
    ckupdc(ielpdc,5) =  0.d0
 
343
    ckupdc(ielpdc,6) =  0.d0
 
344
  enddo
423
345
 
424
346
!-------------------------------------------------------------------------------
425
347
 
426
348
endif
427
349
 
 
350
! Deallocate the temporary array
 
351
deallocate(lstelt)
 
352
 
428
353
return
429
 
 
430
354
end subroutine