Télécharger supgef.eso

Retour à la liste

Numérotation des lignes :

supgef
  1. C SUPGEF SOURCE PV 15/04/10 21:15:34 8474
  2. SUBROUTINE SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ,
  3. & NES,NP,NPG,IAXI,XCOOR,
  4. & LE,KDEB,KFIN,LRV,IDCEN,CMD,IKOMP,
  5. & TN,IPADT,UN,IPADU,NPTD,ANUK,
  6. & WT,WS,HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL,
  7. & DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
  8. C--------------------------------------------------------------------
  9. C Calcul des fonctions tests associées à la formulation Eléments
  10. C Finis Petrov-Galerkin afin de stabiliser les termes de convection.
  11. C On en profite pour calculer des temps caractéristiques.
  12. C--------------------------------------------------------------------
  13. C Cette subroutine étant écrite en fortran pur, certains arguments
  14. C servent uniquement au dimensionnement des tableaux.
  15. C Afin de doper les calculs, on decoupe la boucle sur les éléments
  16. C par paquets afin de bénéficier à fond de la vectorisation.
  17. C--------------------------------------------------------------------
  18. C
  19. C--------------------------
  20. C Paramètre Entrée/Sortie :
  21. C--------------------------
  22. C
  23. C /S FN : Fonction de forme associé à la transformation géométrique
  24. C /S GR : Gradient de FN dans l'élément de référence
  25. C /S PG : Poids d'intégration associé à chaque point de Gauss
  26. C /S XYZ : Coordonnées des sommets de l'élément
  27. C /S AJ : Jacobien
  28. C /S HR : Gradient des fonctions de forme dans l'élément courant
  29. C /S PGSQ : Produit du poids d'intégration par detJ
  30. C /S RPG : Rayon associé à chaque point de Gauss (cas axisymétrique)
  31. C
  32. C E/ NES : Dimension en espace du problème traité (2 en 2D et 3 en 3D)
  33. C E/ NP : Nombre de points support de ddl pour l'élément considéré
  34. C E/ NPG : Nombre de points de Gauss
  35. C E/ IAXI : Précise l'axe de symétrie dans le cas axisymétrique
  36. C (en fait l'axe de symétrie est toujours y -> IAXI=2)
  37. C E/ XCOOR : Coordonnées de l'ensemble des points (voir SMCOORD.INC)
  38. C Petrov-Galerkin choisi pour stabiliser la convection
  39. C Petrov-Galerkin choisi pour stabiliser la convection
  40. C /S HK : Gradient de FN dans l'élément courant
  41. C /S PGSK : Produit du poids d'intégration par detJ
  42. C /S RPGK : Rayon associé à chaque point de Gauss (cas axisymétrique)
  43. C /S AIRE : Volume de chaque élément
  44. C /S AJK : Jacobien
  45. C /S UIL : Vecteur vitesse aux points de Gauss
  46. C /S DUIL : Gradient du champ de vitesse
  47. C (la divergence est obtenue par sommation)
  48. C E/ KDEB : Indice du début de la fenetre
  49. C E/ KFIN : Indice de fin de la fenetre
  50. C E/ LRV : Dimension de la fenetre
  51. C E/ TN : Inconnue transportée au pas de temps précédent
  52. C E/ IPADT : Correspondance numérotation locale et globale pour les
  53. C points SOMMET : J=LECT(I) : le point numéro I est rangé
  54. C en Jème position pour le CHAMP TN
  55. C E/ UN : Champ de vitesse transportant aux sommets des éléments
  56. C E/ IPADU : Correspondance numérotation locale et globale pour les
  57. C points SOMMET : J=LECT(I) : le point numéro I est rangé
  58. C en Jème position pour le CHAMP UN
  59. C E/ NPTD : Nombre de points sous-tendant le spg associé à UN
  60. C E/ ANUK : Champ contenant le coefficient de diffusion (diffusivité
  61. C pour l'énergie, rapport viscosité sur densité pour qdm)
  62. C E/ IDCEN : Contient le schéma associé au terme convectif
  63. C (1=CENTREE 2=SUPGDC 3=SUPG 4=TVISQUEUX 5=CNG)
  64. C E/ DTM1 : Pas de temps imposé
  65. C E/S DT : Pas de temps
  66. C /S DTT1 : Temps caractéristique associé à la diffusion/convection
  67. C /S DTT2 : Temps caractéristique associé à la diffusion
  68. C /S DIAEL : Longueur carac. de l'élément le plus pénalisant pour dt
  69. C /S NUEL : Numéro de l'élément le plus pénalisant pour dt
  70. C
  71. C Variable locale
  72. C
  73. C DDUIL : gradient de chaque composante de la vitesse au pt Gauss
  74. C pour un élément.
  75. C
  76. C RTNL : valeur de TN au point de gauss
  77. C
  78. C--------------------------------------------------------------------
  79. IMPLICIT INTEGER(I-N)
  80. IMPLICIT REAL*8 (A-H,O-Z)
  81.  
  82. -INC PPARAM
  83. -INC CCOPTIO
  84. -INC CCREEL
  85. DIMENSION IPADT(*),IPADU(*),LE(NP,*)
  86.  
  87. DIMENSION XCOOR(*)
  88. DIMENSION UN(NPTD,IDIM),TN(*)
  89. DIMENSION ANUK(LRV)
  90.  
  91. DIMENSION FN(NP,NPG),GR(IDIM,NP,NPG),PG(NPG),XYZ(IDIM,NP)
  92. DIMENSION HR(IDIM,NP,NPG),PGSQ(NPG),RPG(NPG),AJ(IDIM,IDIM,NPG)
  93. DIMENSION WT(LRV,NP,NPG),WS(LRV,NP,NPG),HK(LRV,IDIM,NP,NPG)
  94. DIMENSION PGSK(LRV,NPG),RPGK(LRV,NPG),AIRE(LRV)
  95. DIMENSION AJK(LRV,IDIM,IDIM,NPG)
  96. C
  97. DIMENSION UIL(LRV,IDIM,NPG),DUIL(LRV,IDIM,NPG)
  98. DIMENSION UPIL(3),GRAD(3)
  99. DATA ACCT,ACC2/1.,1./ , ipas/0/
  100. C
  101. REAL*8 DDUIL(3,3),RTNL
  102. C
  103. C
  104. C- Pour les éléments de la fenetre KDEB:KFIN initialisations des
  105. C- données associées à l'élément fini utilisé : fonctions de forme,
  106. C- gradients des fonctions de forme, aire de l'élément, poids des
  107. C- points de Gauss et rayon des points de Gauss pour le cas axi.
  108. C
  109. C? write(6,*)'UN(NPTD,IDIM) ',nptd,idim
  110. C? write(6,1002)un
  111. XPETI=1.D-30
  112. CALL INITD(UPIL,3,0.D0)
  113.  
  114. DO 5020 K=KDEB,KFIN
  115. KP = K-KDEB+1
  116.  
  117. DO 20 I=1,NP
  118. J = LE(I,K)
  119. DO 10 N=1,IDIM
  120. XYZ(N,I) = XCOOR((J-1)*(IDIM+1)+N)
  121. 10 CONTINUE
  122. 20 CONTINUE
  123.  
  124. CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,
  125. & NES,IDIM,NP,NPG,IAXI,AIRE(KP),AJ,SGN)
  126.  
  127. DO 8934 N=1,IDIM
  128. DO 8934 M=1,IDIM
  129. DO 8934 L=1,NPG
  130. AJK(KP,M,N,L)=AJ(M,N,L)
  131. 8934 CONTINUE
  132.  
  133. CALL CALJTR(GR,XYZ,NES,IDIM,NP,NPG,AJ)
  134.  
  135. C
  136. DO 5030 L=1,NPG
  137.  
  138. DO 40 N=1,IDIM
  139. DO 30 I=1,NP
  140. HK(KP,N,I,L) = HR(N,I,L)
  141. 30 CONTINUE
  142. 40 CONTINUE
  143.  
  144. PGSK(KP,L) = PGSQ(L)
  145. RPGK(KP,L) = RPG(L)
  146. C
  147. C- Calcul en chaque élément, pour chaque point de Gauss
  148. C- UIL : Vitesse
  149. C
  150. DO 230 N=1,IDIM
  151. DUIL(KP,N,L) = XPETI
  152. UIL(KP,N,L) = XPETI
  153. DO 210 I=1,NP
  154. NF = IPADU(LE(I,K))
  155. UIL(KP,N,L) = UIL(KP,N,L) + UN(NF,N)*FN(I,L)
  156. DUIL(KP,N,L) = DUIL(KP,N,L) + UN(NF,N)*HK(KP,N,I,L)
  157. 210 CONTINUE
  158. 230 CONTINUE
  159. C
  160. C Dans le cas d'une formulation conservative du transport,
  161. C La vitesse de transport réelle est d(uT)/d(T). Cela revient
  162. C à rajouter le terme T (dU/dT) à la vitesse U.
  163. C Le terme dU/dT peut s'écrite pour chaque composante :
  164. C dx/dT dU/dx + dy/dT dU/dy
  165. C Lorsque dx/dT ou dy/dT est infini, cela signifie que les variables
  166. C sont indépendantes de T. Le terme correspondant doit alors être
  167. C éliminé du caclul.
  168. C
  169. IF (IDCEN .EQ. 2 .OR. IDCEN .EQ. 3) THEN
  170. IF (IKOMP .EQ. 1 .OR. IKOMP .EQ. 2) THEN
  171. C
  172. C Evaluation du gradient de TN.
  173. C
  174. DO 171 N=1,IDIM
  175. GRAD(N)=0.D0
  176. DO 172 I=1,NP
  177. NT = IPADT(LE(I,K))
  178. GRAD(N) = GRAD(N) + TN(NT) *HK(KP,N,I,L)
  179. 172 CONTINUE
  180. IF (ABS(GRAD(N)) .LT. 1.D-10) THEN
  181. C La variable x ou y est indépendante de T. On met le
  182. C gradient à l'infini pour le déconnecter des calculs
  183. C suivants.
  184. GRAD(N) = 1.D+100
  185. ENDIF
  186. 171 CONTINUE
  187. * WRITE(*,*) 'GRADT', GRAD(1),GRAD(2)
  188. C
  189. C On évalue le gradient de chaque composante de la vitesse.
  190. C
  191. DO 231 N=1,IDIM
  192. DO 231 M=1,IDIM
  193. DDUIL(N,M) = XPETI
  194. DO 211 I=1,NP
  195. NF = IPADU(LE(I,K))
  196. DDUIL(N,M) = DDUIL(N,M) + UN(NF,N)*HK(KP,M,I,L)
  197. 211 CONTINUE
  198. 231 CONTINUE
  199. C
  200. C Correction en axisymetrique
  201. C
  202. C? IF (IAXI .EQ. 1) THEN
  203. C? DO 212 I=1,NP
  204. C? NF = IPADU(LE(I,K))
  205. C? DDUIL(1,1)=DDUIL(1,1)+(UN(NF,1)*FN(I,L)/RPGK(KP,L))
  206. C? 212 CONTINUE
  207. C? ENDIF
  208. C
  209. C On évalue la variable scalaire TN au point de gauss L
  210. C
  211. RTNL=0.D0
  212. DO 173 I=1,NP
  213. NT = IPADT(LE(I,K))
  214. RTNL = RTNL + TN(NT) * FN(I,L)
  215. 173 CONTINUE
  216. * WRITE(*,*) 'RTNL', RTNL
  217. C
  218. C Modification de la vitesse U
  219. C
  220. IF (IKOMP .EQ. 1) THEN
  221. c decentrement de UN*TN en UN + TN*(dUN/dTN).
  222. DO 174 N=1,IDIM
  223. DO 174 M=1,IDIM
  224. UIL(KP,N,L) = UIL(KP,N,L)+RTNL*(DDUIL(N,M)/GRAD(M))
  225. 174 CONTINUE
  226. ELSE
  227. c décentrement de JN=UN*TN en dJN/dTN
  228. DO 175 N=1,IDIM
  229. DO 175 M=1,IDIM
  230. UIL(KP,N,L) = DDUIL(N,M)/GRAD(M)
  231. 175 CONTINUE
  232. ENDIF
  233. C
  234. ENDIF
  235. ENDIF
  236. C
  237. C- Calcul pour chaque point de Gauss de chaque élément de :
  238. C- /L UML : Norme de la vitesse aux points de Gauss
  239. C- BM : module d'un temps caractéristique associé à la convection
  240. C- XMB : caractéristique géométrique 1/2 He
  241. C
  242.  
  243. C
  244. BMI=0.D0
  245. UL=0.D0
  246. DO 310 N=1,IDIM
  247. UHAT=0.D0
  248. DO 311 M=1,IDIM
  249. UHAT=UHAT+AJ(M,N,L)*UIL(KP,M,L)
  250. 311 CONTINUE
  251.  
  252. BMI=BMI+UHAT*UHAT
  253. UL=UL+UIL(KP,N,L)*UIL(KP,N,L)
  254. 310 CONTINUE
  255. BM=SQRT(BMI) + XPETI
  256. UML=SQRT(UL)
  257. XMB=UML/BM
  258.  
  259.  
  260. C
  261. C
  262. C- Calcul en chaque élément, pour chaque point de Gauss de
  263. C- GRAD : vecteur unitaire porté par le gradient du champ scalaire
  264. C- UP : projection de la vitesse sur la direction donnée par GRAD
  265. C- UPIL : vecteur UP*GRAD aux points de Gauss
  266. C- pour l'option SUPGDC
  267. C
  268. IF (IDCEN.EQ.2) THEN
  269. DO 170 N=1,IDIM
  270. GRAD(N)=0.D0
  271. DO 170 I=1,NP
  272. NT = IPADT(LE(I,K))
  273. GRAD(N) = GRAD(N) + TN(NT) *HK(KP,N,I,L)
  274. 170 CONTINUE
  275.  
  276. AX=0.D0
  277. DO 2301 M=1,IDIM
  278. AX = AX + GRAD(M)*GRAD(M)
  279. 2301 CONTINUE
  280. AX = SQRT(AX) + XPETI
  281.  
  282. UPL=0.D0
  283. DO 2302 N=1,IDIM
  284. GRAD(N) = GRAD(N) / AX
  285. UPL = UPL + GRAD(N) * UIL(KP,N,L)
  286. 2302 CONTINUE
  287.  
  288. DO 2303 N=1,IDIM
  289. UPIL(N) = GRAD(N) * UPL
  290. 2303 CONTINUE
  291.  
  292. C
  293. BPI=0.D0
  294. DO 410 N=1,IDIM
  295. UPHAT=0.D0
  296. DO 411 M=1,IDIM
  297. UPHAT=UPHAT+AJ(M,N,L)*UPIL(M)
  298. 411 CONTINUE
  299.  
  300. BPI=BPI+UPHAT*UPHAT
  301. 410 CONTINUE
  302. BP=SQRT(BPI) + XPETI
  303. XPB=UPL/BP
  304.  
  305. ENDIF
  306. C
  307. C-----------------------------
  308. C- DECENTREMENT suivant IDCEN
  309. C-----------------------------
  310. C On calcule dans chaque cas TO1 et TO2 ainsi que le tenseur
  311. C associé à la viscosité numérique afin d'évaluer le pas de
  312. C temps de stabilité pour les schémas explicites.
  313. C
  314. C----------
  315. C CENTREE :
  316. C----------
  317. IF (IDCEN.EQ.1) THEN
  318. SI1 = 0.D0
  319. SI2 = 0.D0
  320. TO1 = 0.D0
  321. TO2 = 0.D0
  322. C---------
  323. C SUPGDC : Base théorique dans : A New FE formulation for computational
  324. C fluid dynamics, II Beyond SUPG, HUGHES et al., in Comp.Meth.Appl.Mech.
  325. C Eng., vol 54, pp 341-355 (1986).
  326. C---------
  327. ELSEIF (IDCEN.EQ.2) THEN
  328. SI1 = 1.D0
  329. SI2 = 1.D0
  330. C
  331. C- Approximation "doublement asymptotique" basée sur la vitesse moyenne
  332. C- HMK : Distance basé sur la vitesse moyenne
  333. C- ALFA : Peclet de maille basé sur la vitesse moyenne
  334. C
  335. ALFA = UML * XMB / (ANUK(KP)+XPETI) / 3.D0
  336. AKSI = MIN(1.D0,ALFA)
  337. CCT = AKSI / BM * ACCT * CMD
  338. C
  339. C- Approximation "doublement asymptotique" basée sur la projection de
  340. C- la vitesse sur le gradient du champ scalaire
  341. C- HMK : Distance basée sur la vitesse projetée
  342. C- ALFA : Peclet de maille basé sur la vitesse projetée
  343. C
  344. ALFA = UPL * XPB / (ANUK(KP)+XPETI) / 3.D0
  345. AKSI = MIN(1.D0,ALFA)
  346. CCP = AKSI / BP
  347. C
  348. CPT = CCP - CCT
  349. CC2 = MAX(0.D0,CPT) * ACC2 * CMD
  350. C
  351. TO1 = CCT
  352. TO2 = CC2
  353. C-------
  354. C SUPG :
  355. C-------
  356. ELSEIF (IDCEN.EQ.3) THEN
  357. SI1 = 1.D0
  358. SI2 = 1.D0
  359. C
  360. C- Approximation "doublement asymptotique" basée sur la vitesse moyenne
  361. C- HMK : Distance basé sur la vitesse moyenne
  362. C- ALFA : Peclet de maille basé sur la vitesse moyenne
  363. C
  364. ALFA = UML * XMB / (ANUK(KP)+XPETI) / 3.D0
  365. AKSI = MIN(1.D0,ALFA)
  366. CCT = AKSI / BM * ACCT * CMD
  367. C
  368. TO1 = CCT
  369. TO2 = 0.D0
  370. C-------------------
  371. C Tenseur Visqueux :
  372. C-------------------
  373. ELSEIF (IDCEN.EQ.4) THEN
  374. SI1 =-1.D0
  375. SI2 = 1.D0
  376. DT19 = DTM1 * 0.5D0
  377. TO1 = DT19
  378. TO2 = 0.D0
  379. C-----------------------------
  380. C Crank Nicholson généralisé :
  381. C-----------------------------
  382. ELSEIF (IDCEN.EQ.5) THEN
  383. SI1 =-1.D0
  384. SI2 = 1.D0
  385. DT19 = DTM1/6.D0
  386. TO1 = DT19
  387. TO2 = 0.D0
  388. ENDIF
  389. C
  390. C---------------------------
  391. C Pas de temps de stabilité
  392. C---------------------------
  393. C
  394. C La viscosité utilisée dans l'évaluation des pas de temps de stabilité
  395. C explicite ajoute à la viscosité physique la viscosité numérique :
  396. C DT1 : Correspond à une CFL de 1 et un Peclet de maille de 1
  397. C CFL=1 -> dx=Vdt |
  398. C |-> dt=2D/V2
  399. C Pe=1 -> dx=2D/V |
  400. C DT2 : Correspond au pas de temps de stabilité lié à la diffusion
  401. C FOU=1 -> dt=0.5 dx2/D
  402. C
  403. IF (IDCEN.NE.5) THEN
  404. DT0 = DT
  405. C DT1 = 2.D0 /
  406. C & ( UIL(KP,1,L)*UIL(KP,1,L)/(ANUK(KP)+XPETI)
  407. C & + UIL(KP,2,L)*UIL(KP,2,L)/(ANUK(KP)+XPETI) )
  408. C DT2 = 0.5D0 /
  409. C & ( ANUK(KP)*AL2(KP)
  410. C & + ANUK(KP)*AH2(KP) )
  411. * faut il provoquer une erreur?
  412. if (abs(uml).lt.xpetit) then
  413. dt1=xgrand
  414. else
  415. DT1 = 2.D0 * ANUK(KP) / (UML * UML)
  416. endif
  417. DT2 = DT1
  418. C
  419. IF (DT1.LT.DT) DT=DT1
  420. IF (DT2.LT.DT) DT=DT2
  421. IF (DT.NE.DT0) THEN
  422. DTT1 = DT1
  423. DTT2 = DT2
  424. DIAEL = XMB * 2.D0
  425. NUEL = K
  426. ENDIF
  427. ENDIF
  428. C
  429. C---------------------------------------------------
  430. C Fonction test pour la formulation Petrov-Galerkin
  431. C---------------------------------------------------
  432. C Ce qui est diffusion numérique en explicite se transforme en
  433. C ajoutant de la viscosité numérique (Tenseurs visqueux et CNG).
  434. C WS : Fonction test pour la partie explicite
  435. C
  436. IF(IDIM.EQ.2)THEN
  437. DO 2050 I=1,NP
  438. WT(KP,I,L) = FN(I,L) + SI1 *
  439. & (TO1*(UIL(KP,1,L)*HK(KP,1,I,L)+UIL(KP,2,L)*HK(KP,2,I,L))
  440. & +TO2*(UPIL(1)*HK(KP,1,I,L)+UPIL(2)*HK(KP,2,I,L)))
  441. WS(KP,I,L) = FN(I,L) + SI2 *
  442. & (TO1*(UIL(KP,1,L)*HK(KP,1,I,L)+UIL(KP,2,L)*HK(KP,2,I,L))
  443. & +TO2*(UPIL(1)*HK(KP,1,I,L)+UPIL(2)*HK(KP,2,I,L)))
  444. 2050 CONTINUE
  445. C
  446. ELSEIF(IDIM.EQ.3)THEN
  447. DO 2051 I=1,NP
  448. WT(KP,I,L) = FN(I,L) + SI1 *
  449. & (TO1*(UIL(KP,1,L)*HK(KP,1,I,L)+UIL(KP,2,L)*HK(KP,2,I,L)+
  450. & UIL(KP,3,L)*HK(KP,3,I,L) )
  451. & +TO2*(UPIL(1)*HK(KP,1,I,L)+UPIL(2)*HK(KP,2,I,L)+
  452. & UPIL(3)*HK(KP,3,I,L) ))
  453.  
  454. WS(KP,I,L) = FN(I,L) + SI2 *
  455. & (TO1*(UIL(KP,1,L)*HK(KP,1,I,L)+UIL(KP,2,L)*HK(KP,2,I,L)+
  456. & UIL(KP,3,L)*HK(KP,3,I,L))
  457. & + TO2*(UPIL(1)*HK(KP,1,I,L)+UPIL(2)*HK(KP,2,I,L)+
  458. & UPIL(3)*HK(KP,3,I,L) ))
  459. 2051 CONTINUE
  460. ENDIF
  461. C
  462. C- Si on est en conservatif, on rétablit les valeurs de la vitesse
  463. C de transport qui ont été modifiées car elles sont utilisées
  464. C dans d'autres subroutines ultérieures.
  465. C
  466. IF (IKOMP .EQ. 1) THEN
  467. DO 235 N=1,IDIM
  468. UIL(KP,N,L) = XPETI
  469. DO 215 I=1,NP
  470. NF = IPADU(LE(I,K))
  471. UIL(KP,N,L) = UIL(KP,N,L) + UN(NF,N)*FN(I,L)
  472. 215 CONTINUE
  473. 235 CONTINUE
  474. ENDIF
  475. C
  476. 5030 CONTINUE
  477.  
  478. 5020 CONTINUE
  479.  
  480. RETURN
  481. 1001 FORMAT(20(1X,I5))
  482. 1002 FORMAT(10(1X,1PE11.4))
  483. END
  484.  
  485.  
  486.  
  487.  
  488.  
  489.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales