← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/alge/gradrv.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:
src/alge/gradrv.f90
 
1
!-------------------------------------------------------------------------------
 
2
 
 
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.
 
20
 
 
21
!-------------------------------------------------------------------------------
 
22
 
 
23
subroutine gradrv &
 
24
!================
 
25
 ( ncelet , ncel   , nfac   , nfabor ,                            &
 
26
   imrgra , inc    , nswrgp ,                                     &
 
27
   iwarnp , nfecra , epsrgp , extrap ,                            &
 
28
   ifacel , ifabor , ivar   ,                                     &
 
29
   volume , surfac , surfbo , pond   , xyzcen , cdgfac , cdgfbo , &
 
30
   dijpf  , diipb  , dofij  ,                                     &
 
31
   coefav , coefbv , pvar   ,                                     &
 
32
   cocg   ,                                                       &
 
33
   gradv  )
 
34
 
 
35
!===============================================================================
 
36
! FONCTION :
 
37
! ----------
 
38
 
 
39
! COMPUTE THE GRADIENT OF A VECTOR WITH AN ITERATIVE TECHNIC IN ORDER TO HANDLE
 
40
! WITH NONE ORTHOGANALITIES (NSWRGP>1). WE DO NOT TAKE INTO ACCOUNT ANY VOLUMIC
 
41
! FORCE HERE.
 
42
 
 
43
! COMPUTE COCG AT THE FISRT CALL AND IF NEEDED.
 
44
 
 
45
!-------------------------------------------------------------------------------
 
46
! Arguments
 
47
!__________________.____._____.________________________________________________.
 
48
! name             !type!mode ! role                                           !
 
49
!__________________!____!_____!________________________________________________!
 
50
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
 
51
! ncel             ! i  ! <-- ! number of cells                                !
 
52
! nfac             ! i  ! <-- ! number of interior faces                       !
 
53
! nfabor           ! i  ! <-- ! number of boundary faces                       !
 
54
! ivar             ! e  ! <-- ! numero de la variable                          !
 
55
!                  !    !     !   destine a etre utilise pour la               !
 
56
!                  !    !     !   periodicite uniquement (pering)              !
 
57
!                  !    !     !   on pourra donner ivar=0 si la                !
 
58
!                  !    !     !   variable n'est ni une composante de          !
 
59
!                  !    !     !   la vitesse, ni une composante du             !
 
60
!                  !    !     !   tenseur des contraintes rij                  !
 
61
! imrgra           ! e  ! <-- ! methode de reconstruction du gradient          !
 
62
!                  !    !     !  0 reconstruction 97                           !
 
63
!                  !    !     !  1 moindres carres                             !
 
64
!                  !    !     !  2 moindres carres support etendu              !
 
65
!                  !    !     !    complet                                     !
 
66
!                  !    !     !  3 moindres carres avec selection du           !
 
67
!                  !    !     !    support etendu                              !
 
68
! inc              ! e  ! <-- ! indicateur = 0 resol sur increment             !
 
69
!                  !    !     !              1 sinon                           !
 
70
! nswrgp           ! e  ! <-- ! nombre de sweep pour reconstruction            !
 
71
!                  !    !     !             des gradients                      !
 
72
! imligp           ! e  ! <-- ! methode de limitation du gradient              !
 
73
!                  !    !     !  < 0 pas de limitation                         !
 
74
!                  !    !     !  = 0 a partir des gradients voisins            !
 
75
!                  !    !     !  = 1 a partir du gradient moyen                !
 
76
! iwarnp           ! i  ! <-- ! verbosity                                      !
 
77
! nfecra           ! e  ! <-- ! unite du fichier sortie std                    !
 
78
! epsrgp           ! r  ! <-- ! precision relative pour la                     !
 
79
!                  !    !     !  reconstruction des gradients 97               !
 
80
! climgp           ! r  ! <-- ! coef gradient*distance/ecart                   !
 
81
! extrap           ! r  ! <-- ! coef extrap gradient                           !
 
82
! volume(ncelet)   ! tr ! <-- ! volume des elements                            !
 
83
! surfac(3,nfac)   ! tr ! <-- ! surf vectorielle des surfaces interne          !
 
84
! surfbo           ! tr ! <-- ! surf vectorielle des surfaces de bord          !
 
85
!   (3,nfabor)     !    !     !                                                !
 
86
! pond(nfac)       ! tr ! <-- ! ponderation geometrique (entre 0 et 1          !
 
87
! xyzcen           ! tr ! <-- ! point associes aux volumes de control          !
 
88
! (3,ncelet)       !    !     !                                                !
 
89
! cdgfac           ! tr ! <-- ! centre de gravite des faces internes           !
 
90
! (3,nfac)         !    !     !                                                !
 
91
! cdgfbo           ! tr ! <-- ! centre de gravite des faces de bord            !
 
92
! (3,nfabor)       !    !     !                                                !
 
93
! dijpf(3,nfac)    ! tr ! <-- ! vect i'j', i' (resp. j') projection            !
 
94
!                  !    !     !  du centre i (resp. j) sur la normale          !
 
95
!                  !    !     !  a la face interne                             !
 
96
! diipb            ! tr ! <-- ! vect ii', ii projection du centre i            !
 
97
!   (3,nfabor)     !    !     !  sur la normale a la face de bord              !
 
98
! cocg(ncelet,3    ! tr ! <-- ! couplage des composantes du gradient           !
 
99
!     ,3)          !    !     ! lors de la reconstruction                      !
 
100
! dofij(3,nfac)    ! tr ! --> ! vecteur of pour les faces internes             !
 
101
!                  !    !     ! o : intersection de ij et la face              !
 
102
!                  !    !     ! f : centre de la face                          !
 
103
! pond(nfac)       ! tr ! <-- ! ponderation geometrique (entre 0 et 1)         !
 
104
! pvar(3,ncelet)    ! tr ! <-- ! variable (vitesse)                             !
 
105
! coefav,coefbv    ! tr ! <-- ! tableaux des cond lim pour pvar                 !
 
106
!   (3,nfabor)     !    !     !  sur la normale a la face de bord              !
 
107
! gradv            ! tr ! --> ! gradient de vitesse                            !
 
108
!   (3,3,ncelet)   !    !     !                                                !
 
109
! gradva           ! tr ! --- ! tableau de travail pour le grad de vitesse     !
 
110
!   (3,3,ncelet)   !    !     !                                                !
 
111
!__________________!____!_____!________________________________________________!
 
112
 
 
113
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
 
114
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
 
115
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
 
116
!            --- tableau de travail
 
117
!===============================================================================
 
118
 
 
119
!===============================================================================
 
120
! Module files
 
121
!===============================================================================
 
122
 
 
123
use paramx
 
124
use cstphy
 
125
use cstnum
 
126
use albase
 
127
use cplsat
 
128
use parall
 
129
use period
 
130
 
 
131
!===============================================================================
 
132
 
 
133
implicit none
 
134
 
 
135
! Arguments
 
136
 
 
137
integer          ncelet , ncel   , nfac   , nfabor
 
138
integer          ivar   , imrgra , inc    , nswrgp
 
139
integer          imligp , iwarnp , nfecra
 
140
double precision epsrgp , climgp , extrap
 
141
 
 
142
integer          ifacel(2,nfac), ifabor(nfabor)
 
143
double precision volume(ncelet), surfac(3,nfac)
 
144
double precision surfbo(3,nfabor)
 
145
double precision pond(nfac)
 
146
double precision xyzcen(3,ncelet)
 
147
double precision cdgfac(3,nfac),cdgfbo(3,nfabor)
 
148
double precision dijpf(3,nfac), diipb(3,nfabor)
 
149
double precision cocg(3,3,ncelet)
 
150
double precision dofij(3,nfac)
 
151
 
 
152
double precision pvar(3,ncelet), coefav(3,nfabor), coefbv(3,3,nfabor)
 
153
double precision gradv (3,3,ncelet)
 
154
 
 
155
! Local variables
 
156
 
 
157
integer          imlini
 
158
 
 
159
double precision climin
 
160
 
 
161
integer          lbloc
 
162
parameter       (lbloc = 1024)
 
163
integer          iel, ifac, ii, jj, kk, ll, mm
 
164
integer          isou, jsou
 
165
integer          isqrt, isweep, nswmax
 
166
integer          ibloc,nbloc,irel
 
167
double precision pfac,pip,deltpx,deltpy,deltpz
 
168
double precision rnorx,rnory,rnorz,rnorm,residu
 
169
double precision pfaci
 
170
double precision dof,fmoyx,fmoyy,fmoyz
 
171
double precision aa(3,3,lbloc),unsdet,pfsx,pfsy,pfsz
 
172
double precision a11,a12,a13,a21,a22,a23,a31,a32,a33
 
173
double precision cocg11,cocg12,cocg13,cocg21,cocg22,cocg23
 
174
double precision cocg31,cocg32,cocg33
 
175
double precision unsvol, vecfac
 
176
 
 
177
double precision, dimension(:,:,:), allocatable :: gradva
 
178
 
 
179
integer          ipassv
 
180
data             ipassv /0/
 
181
save             ipassv
 
182
 
 
183
!===============================================================================
 
184
 
 
185
!===============================================================================
 
186
! 0. PREPARATION POUR PERIODICITE DE ROTATION
 
187
!===============================================================================
 
188
 
 
189
! Par defaut, on traitera le gradient comme un tenseur ...
 
190
!   (i.e. on suppose que c'est le gradient d'une grandeurs vectorielle)
 
191
 
 
192
if(iperio.eq.1) then
 
193
  call synvin(pvar)
 
194
  !==========
 
195
endif
 
196
 
 
197
 
 
198
!===============================================================================
 
199
! 1. INITIALISATION
 
200
!===============================================================================
 
201
 
 
202
isqrt = 1
 
203
nswmax = nswrgp
 
204
isweep = 1
 
205
 
 
206
! Allocate temporary array
 
207
allocate(gradva(3,3,ncelet))
 
208
 
 
209
!===============================================================================
 
210
! 2. CALCUL SANS RECONSTRUCTION
 
211
!===============================================================================
 
212
 
 
213
!     SI INITIALISATION PAR MOINDRES CARRES (IMRGRA=4), B EST DEJA REMPLI
 
214
!     SINON (IMRGRA=0) ON CALCULE UN GRADIENT SANS RECONSTRUCTION
 
215
 
 
216
if (imrgra.eq.0) then
 
217
 
 
218
 
 
219
  do iel = 1, ncelet
 
220
  ! BOUCLE SUR LES COMPOSANTES X, Y, Z
 
221
    do jsou = 1, 3
 
222
    ! BOUCLE SUR LES COMPOSANTES U, V, W
 
223
      do isou = 1 ,3
 
224
        gradva(isou,jsou,iel) = 0.d0
 
225
      enddo
 
226
    enddo
 
227
  enddo
 
228
 
 
229
  !     ASSEMBLAGE A PARTIR DES FACETTES FLUIDES
 
230
 
 
231
  do ifac = 1,nfac
 
232
  ! BOUCLE SUR LES COMPOSANTES U, V, W
 
233
    ii = ifacel(1,ifac)
 
234
    jj = ifacel(2,ifac)
 
235
    do isou = 1 ,3
 
236
      pfac   = pond(ifac)*pvar(isou,ii) +(1.d0-pond(ifac))*pvar(isou,jj)
 
237
      do jsou = 1, 3
 
238
        gradva(isou,jsou,ii) = gradva(isou,jsou,ii) + pfac*surfac(jsou,ifac)
 
239
        gradva(isou,jsou,jj) = gradva(isou,jsou,jj) - pfac*surfac(jsou,ifac)
 
240
      enddo
 
241
    enddo
 
242
  enddo
 
243
 
 
244
  !     ASSEMBLAGE A PARTIR DES FACETTES DE BORD
 
245
 
 
246
  do ifac = 1,nfabor
 
247
  ! BOUCLE SUR LES COMPOSANTES U, V, W
 
248
    ii = ifabor(ifac)
 
249
    do isou = 1 ,3
 
250
      pfac = inc*coefav(isou,ifac) + coefbv(isou,1,ifac)*pvar(1,ii)    &
 
251
                                   + coefbv(isou,2,ifac)*pvar(2,ii)    &
 
252
                                   + coefbv(isou,3,ifac)*pvar(3,ii)
 
253
      do jsou = 1, 3
 
254
        gradva(isou,jsou,ii) = gradva(isou,jsou,ii) + pfac*surfbo(jsou,ifac)
 
255
      enddo
 
256
    enddo
 
257
  enddo
 
258
 
 
259
 
 
260
  ! GRAVEL = GRADIENT
 
261
 
 
262
  do iel = 1, ncel
 
263
    unsvol = 1.d0/volume(iel)
 
264
    do jsou = 1, 3
 
265
      do isou = 1, 3
 
266
        gradv(isou,jsou,iel) = gradva(isou,jsou,iel)*unsvol
 
267
      enddo
 
268
    enddo
 
269
  enddo
 
270
 
 
271
!     TRAITEMENT DU PARALLELISME ET DE LA PERIODICITE
 
272
 
 
273
  if(irangp.ge.0.or.iperio.eq.1) then
 
274
    call syntin (gradv)
 
275
    !==========
 
276
  endif
 
277
 
 
278
endif
 
279
 
 
280
if( nswrgp.le.1 ) then
 
281
  ! Free memory
 
282
  deallocate(gradva)
 
283
 
 
284
  return
 
285
endif
 
286
 
 
287
 
 
288
!     On incremente IPASS quand on calcule COCG pour la premiere fois
 
289
ipassv = ipassv + 1
 
290
 
 
291
!===============================================================================
 
292
! 3. RECONSTRUCTION DES GRADIENTS POUR LES MAILLAGES TORDUS
 
293
!===============================================================================
 
294
 
 
295
!     RESOLUTION SEMI-IMPLICITE SUR TOUT LE MAILLAGE
 
296
!     gradv = GRADIENT
 
297
 
 
298
 
 
299
if(ipassv.eq.1 .or. iale.eq.1 .or. imobil.eq.1) then
 
300
 
 
301
! ---> CALCUL DE COCG
 
302
 
 
303
  do iel =1,ncelet
 
304
    do jj = 1, 3
 
305
      do ii = 1, 3
 
306
        cocg(ii,jj,iel)= 0.d0
 
307
      enddo
 
308
    enddo
 
309
  enddo
 
310
  do iel=1,ncel
 
311
    cocg(1,1,iel) = volume(iel)
 
312
    cocg(2,2,iel) = volume(iel)
 
313
    cocg(3,3,iel) = volume(iel)
 
314
  enddo
 
315
 
 
316
! ---> AJOUT DES CONTRIBUTIONS DES FACES INTERNES
 
317
! ATTENTION, LA MATRICE 3X3 COCG(iel) EST ICI LA TRANSPOSEE
 
318
!            DE CELLE UTILISEE DANS LE GRADIENT SCALAIRE ET NE PREND PAS EN
 
319
!            COMPTE LES TERMES DE BORD
 
320
 
 
321
  do ifac=1,nfac
 
322
    ii = ifacel(1,ifac)
 
323
    jj = ifacel(2,ifac)
 
324
 
 
325
    do isou =1,3
 
326
      !---> DOF = OF
 
327
      dof = dofij(isou,ifac)
 
328
 
 
329
      pfaci  = -dof*0.5d0
 
330
      do jsou =1,3
 
331
        vecfac = pfaci*surfac(jsou,ifac)
 
332
        cocg(isou,jsou,ii) = cocg(isou,jsou,ii) + vecfac
 
333
        cocg(isou,jsou,jj) = cocg(isou,jsou,jj) - vecfac
 
334
      enddo
 
335
    enddo
 
336
 
 
337
  enddo
 
338
 
 
339
 
 
340
! prise en compte explicite des faces de bords: pas de recalcul de COCG
 
341
 
 
342
! ---> ON INVERSE POUR TOUTE LES CELLULES : LE COCG POUR TOUTES LES CELLULES
 
343
!      RESTE ENSUITE LE MEME TANT QUE LE MAILLAGE NE CHANGE PAS
 
344
  nbloc = ncel/lbloc
 
345
  if (nbloc.gt.0) then
 
346
    do ibloc = 1, nbloc
 
347
      do ii =1, lbloc
 
348
        iel = (ibloc-1)*lbloc+ii
 
349
 
 
350
        cocg11 = cocg(1,1,iel)
 
351
        cocg12 = cocg(1,2,iel)
 
352
        cocg13 = cocg(1,3,iel)
 
353
        cocg21 = cocg(2,1,iel)
 
354
        cocg22 = cocg(2,2,iel)
 
355
        cocg23 = cocg(2,3,iel)
 
356
        cocg31 = cocg(3,1,iel)
 
357
        cocg32 = cocg(3,2,iel)
 
358
        cocg33 = cocg(3,3,iel)
 
359
 
 
360
        a11=cocg22*cocg33-cocg32*cocg23
 
361
        a12=cocg32*cocg13-cocg12*cocg33
 
362
        a13=cocg12*cocg23-cocg22*cocg13
 
363
        a21=cocg31*cocg23-cocg21*cocg33
 
364
        a22=cocg11*cocg33-cocg31*cocg13
 
365
        a23=cocg21*cocg13-cocg11*cocg23
 
366
        a31=cocg21*cocg32-cocg31*cocg22
 
367
        a32=cocg31*cocg12-cocg11*cocg32
 
368
        a33=cocg11*cocg22-cocg21*cocg12
 
369
 
 
370
        unsdet = 1.d0/(cocg11*a11 +cocg21*a12+cocg31*a13)
 
371
 
 
372
        aa(1,1,ii) = a11*unsdet
 
373
        aa(1,2,ii) = a12*unsdet
 
374
        aa(1,3,ii) = a13*unsdet
 
375
        aa(2,1,ii) = a21*unsdet
 
376
        aa(2,2,ii) = a22*unsdet
 
377
        aa(2,3,ii) = a23*unsdet
 
378
        aa(3,1,ii) = a31*unsdet
 
379
        aa(3,2,ii) = a32*unsdet
 
380
        aa(3,3,ii) = a33*unsdet
 
381
 
 
382
      enddo
 
383
 
 
384
      do ii = 1, lbloc
 
385
        iel = (ibloc-1)*lbloc+ii
 
386
        cocg(1,1,iel) = aa(1,1,ii)
 
387
        cocg(1,2,iel) = aa(1,2,ii)
 
388
        cocg(1,3,iel) = aa(1,3,ii)
 
389
        cocg(2,1,iel) = aa(2,1,ii)
 
390
        cocg(2,2,iel) = aa(2,2,ii)
 
391
        cocg(2,3,iel) = aa(2,3,ii)
 
392
        cocg(3,1,iel) = aa(3,1,ii)
 
393
        cocg(3,2,iel) = aa(3,2,ii)
 
394
        cocg(3,3,iel) = aa(3,3,ii)
 
395
      enddo
 
396
 
 
397
    enddo
 
398
 
 
399
  endif
 
400
 
 
401
  irel = mod(ncel,lbloc)
 
402
  if (irel.gt.0) then
 
403
    ibloc = nbloc + 1
 
404
    do ii = 1, irel
 
405
      iel = (ibloc-1)*lbloc+ii
 
406
 
 
407
      cocg11 = cocg(1,1,iel)
 
408
      cocg12 = cocg(1,2,iel)
 
409
      cocg13 = cocg(1,3,iel)
 
410
      cocg21 = cocg(2,1,iel)
 
411
      cocg22 = cocg(2,2,iel)
 
412
      cocg23 = cocg(2,3,iel)
 
413
      cocg31 = cocg(3,1,iel)
 
414
      cocg32 = cocg(3,2,iel)
 
415
      cocg33 = cocg(3,3,iel)
 
416
 
 
417
      a11=cocg22*cocg33-cocg32*cocg23
 
418
      a12=cocg32*cocg13-cocg12*cocg33
 
419
      a13=cocg12*cocg23-cocg22*cocg13
 
420
      a21=cocg31*cocg23-cocg21*cocg33
 
421
      a22=cocg11*cocg33-cocg31*cocg13
 
422
      a23=cocg21*cocg13-cocg11*cocg23
 
423
      a31=cocg21*cocg32-cocg31*cocg22
 
424
      a32=cocg31*cocg12-cocg11*cocg32
 
425
      a33=cocg11*cocg22-cocg21*cocg12
 
426
 
 
427
      unsdet = 1.d0/(cocg11*a11 +cocg21*a12+cocg31*a13)
 
428
 
 
429
      aa(1,1,ii) = a11*unsdet
 
430
      aa(1,2,ii) = a12*unsdet
 
431
      aa(1,3,ii) = a13*unsdet
 
432
      aa(2,1,ii) = a21*unsdet
 
433
      aa(2,2,ii) = a22*unsdet
 
434
      aa(2,3,ii) = a23*unsdet
 
435
      aa(3,1,ii) = a31*unsdet
 
436
      aa(3,2,ii) = a32*unsdet
 
437
      aa(3,3,ii) = a33*unsdet
 
438
 
 
439
    enddo
 
440
 
 
441
    do ii = 1, irel
 
442
      iel = (ibloc-1)*lbloc+ii
 
443
      cocg(1,1,iel) = aa(1,1,ii)
 
444
      cocg(1,2,iel) = aa(1,2,ii)
 
445
      cocg(1,3,iel) = aa(1,3,ii)
 
446
      cocg(2,1,iel) = aa(2,1,ii)
 
447
      cocg(2,2,iel) = aa(2,2,ii)
 
448
      cocg(2,3,iel) = aa(2,3,ii)
 
449
      cocg(3,1,iel) = aa(3,1,ii)
 
450
      cocg(3,2,iel) = aa(3,2,ii)
 
451
      cocg(3,3,iel) = aa(3,3,ii)
 
452
    enddo
 
453
 
 
454
  endif
 
455
 
 
456
endif
 
457
 
 
458
! ---> CALCUL DU RESIDU DE NORMALISATION
 
459
 
 
460
! TODO a opimized form prodcsc9
 
461
! Euclidian norm
 
462
call prodsc(9*ncelet,9*ncel,isqrt,gradva(1,1,1),gradva(1,1,1),       &
 
463
                              rnorm )
 
464
 
 
465
! TODO non-dimensional criterium
 
466
if (volmax.gt.1.d0) rnorm = rnorm / volmax
 
467
if( rnorm.le.epzero ) then
 
468
  ! Free memory
 
469
  deallocate(gradva)
 
470
  return
 
471
endif
 
472
 
 
473
!  LE VECTEUR OijFij EST CALCULE DANS CLDIJP
 
474
 
 
475
! ---> DEBUT DES ITERATIONS
 
476
 
 
477
 100  continue
 
478
 
 
479
isweep = isweep +1
 
480
 
 
481
 
 
482
!     CALCUL DU SECOND MEMBRE
 
483
 
 
484
 
 
485
do iel = 1, ncel
 
486
  do jsou =1, 3
 
487
    do isou = 1, 3
 
488
      gradva(isou,jsou,iel) = -gradv(isou,jsou,iel)*volume(iel)
 
489
    enddo
 
490
  enddo
 
491
enddo
 
492
 
 
493
 
 
494
!     ASSEMBLAGE A PARTIR DES FACETTES FLUIDES
 
495
 
 
496
do ifac = 1,nfac
 
497
  ii = ifacel(1,ifac)
 
498
  jj = ifacel(2,ifac)
 
499
  do isou = 1, 3
 
500
    vecfac = pond(ifac)*pvar(isou,ii) +(1.d0-pond(ifac))*pvar(isou,jj)         &
 
501
           +0.5d0*(  (gradv(isou,1,ii)+gradv(isou,1,jj))*dofij(1,ifac)   &
 
502
                    +(gradv(isou,2,ii)+gradv(isou,2,jj))*dofij(2,ifac)   &
 
503
                    +(gradv(isou,3,ii)+gradv(isou,3,jj))*dofij(3,ifac))
 
504
    do jsou = 1, 3
 
505
      gradva(isou,jsou,ii) = gradva(isou,jsou,ii)             &
 
506
                             + vecfac*surfac(jsou,ifac)
 
507
      gradva(isou,jsou,jj) = gradva(isou,jsou,jj)             &
 
508
                             - vecfac*surfac(jsou,ifac)
 
509
    enddo
 
510
  enddo
 
511
enddo
 
512
 
 
513
 
 
514
!     ASSEMBLAGE A PARTIR DES FACETTES DE BORD
 
515
 
 
516
do ifac = 1,nfabor
 
517
  ii = ifabor(ifac)
 
518
 
 
519
  do isou = 1,3
 
520
    pfac = inc*coefav(isou,ifac)
 
521
 
 
522
    ! coefbv is a matrix
 
523
    do jsou = 1, 3
 
524
      pip =  pvar(jsou,ii)                                   &
 
525
         + gradv(jsou,1,ii)*diipb(1,ifac)                 &
 
526
         + gradv(jsou,2,ii)*diipb(2,ifac)                 &
 
527
         + gradv(jsou,3,ii)*diipb(3,ifac)
 
528
      pfac = pfac + coefbv(isou,jsou,ifac)*pip
 
529
    enddo
 
530
 
 
531
    do jsou = 1, 3
 
532
      gradva(isou,jsou,ii) = gradva(isou,jsou,ii) +pfac*surfbo(jsou,ifac)
 
533
    enddo
 
534
  enddo
 
535
enddo
 
536
 
 
537
 
 
538
!     INCREMENTATION DU GRADIENT
 
539
 
 
540
do iel =1,ncel
 
541
 
 
542
  do jsou = 1, 3
 
543
    do isou = 1, 3
 
544
 
 
545
      gradv(isou,jsou,iel) = gradv(isou,jsou,iel)             &
 
546
                        + gradva(isou,1,iel) *cocg(1,jsou,iel)  &
 
547
                        + gradva(isou,2,iel) *cocg(2,jsou,iel)  &
 
548
                        + gradva(isou,3,iel) *cocg(3,jsou,iel)
 
549
 
 
550
    enddo
 
551
  enddo
 
552
 
 
553
enddo
 
554
 
 
555
 
 
556
!     TRAITEMENT DU PARALLELISME
 
557
 
 
558
if(irangp.ge.0.or.iperio.eq.1) then
 
559
  call syntin (gradv)
 
560
endif
 
561
 
 
562
! ---> TEST DE CONVERGENCE
 
563
 
 
564
! Norme Euclidienne
 
565
call prodsc(9*ncelet,9*ncel,isqrt,gradva(1,1,1),gradva(1,1,1),       &
 
566
                              residu         )
 
567
 
 
568
! TODO a non-diemensional criterium
 
569
if (volmax.gt.1.d0) residu = residu / volmax
 
570
 
 
571
if( residu.le.epsrgp*rnorm) then
 
572
  if( iwarnp.ge.2 ) then
 
573
    write (nfecra,1000) isweep,residu/rnorm,rnorm,ivar
 
574
  endif
 
575
  goto 101
 
576
elseif( isweep.ge.nswmax ) then
 
577
  if( iwarnp.ge.0) then
 
578
     write (nfecra,1000)isweep,residu/rnorm,rnorm,ivar
 
579
     write (nfecra,1100)
 
580
  endif
 
581
  goto 101
 
582
else
 
583
  goto 100
 
584
endif
 
585
 
 
586
 101  continue
 
587
 
 
588
 
 
589
! Free memory
 
590
deallocate(gradva)
 
591
 
 
592
!--------
 
593
! FORMATS
 
594
!--------
 
595
#if defined(_CS_LANG_FR)
 
596
 
 
597
 1000 format(1X,'GRDVEC ISWEEP = ',I4,' RESIDU NORME: ',E11.4,          &
 
598
         ' NORME: ',E11.4,/,1X,'PARAMETRE IVAR = ',I4 )
 
599
 1100 format(                                                           &
 
600
'@                                                            ',/,&
 
601
'@ @@ ATTENTION :         NON CONVERGENCE DE GRDVEC           ',/,&
 
602
'@    =========                                               ',/,&
 
603
'@                                                            '  )
 
604
 
 
605
#else
 
606
 
 
607
 1000 format(1X,'GRDVEC ISWEEP = ',I4,' NORMED RESIDUAL: ',E11.4,       &
 
608
         ' NORM: ',E11.4,/,1X,'PARAMETER IVAR = ',I4 )
 
609
 1100 format(                                                           &
 
610
'@'                                                            ,/,&
 
611
'@ @@ WARNING:            NON CONVERGENCE OF GRDVEC'           ,/,&
 
612
'@    ========'                                                ,/,&
 
613
'@'                                                              )
 
614
 
 
615
#endif
 
616
 
 
617
!----
 
618
! FIN
 
619
!----
 
620
 
 
621
return
 
622
 
 
623
end subroutine