Télécharger gresp.procedur

Retour à la liste

Numérotation des lignes :

  1. * GRESP PROCEDUR GOUNAND 12/02/01 21:15:14 7266
  2. ************************************************************************
  3. * NOM : GRESP
  4. * DESCRIPTION : Résout un système par une méthode de projection
  5. * algébrique incrémentale
  6. *
  7. * LANGAGE : GIBIANE-CAST3M
  8. * AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  9. * mél : gounand@semt2.smts.cea.fr
  10. **********************************************************************
  11. * VERSION : v1, 22/11/2007, version initiale
  12. * HISTORIQUE : v1, 22/11/2007, création
  13. * HISTORIQUE : 21/12/2011 : correction d'un bug dans le
  14. * préconditionnement, nettoyage KOPS
  15. * HISTORIQUE :
  16. ************************************************************************
  17. *
  18. 'DEBPROC' GRESP ;
  19. 'ARGUMENT' mat*'MATRIK' ;
  20. 'ARGUMENT' ccli*'CHPOINT' ;
  21. 'ARGUMENT' smbr*'CHPOINT' ;
  22. 'ARGUMENT' rv*'TABLE' ;
  23. *
  24. * Nom des inconnues : vitesses et pressions
  25. * Dans rvg . 'METHINV', on stocke l'éventuelle table d'inversion
  26. * pour le laplacien de pression
  27. *
  28. debug = FAUX ;
  29. rvg = rv . 'GPROJ' ;
  30. rvm = rv . 'METHINV' ;
  31. 'SI' ('EXISTE' rvg 'dblproj') ;
  32. dblproj = rvg . 'dblproj' ;
  33. 'SINON' ;
  34. dblproj = VRAI ;
  35. 'FINSI' ;
  36. *
  37. 'SI' ('EXISTE' rvg 'NOPREC') ;
  38. noprec = rvg . 'NOPREC' ;
  39. 'SINON' ;
  40. noprec = faux ;
  41. 'FINSI' ;
  42. 'SI' debug ;
  43. 'MESSAGE' ('CHAINE' 'noprec=' noprec) ;
  44. 'FINSI' ;
  45. *
  46. 'SI' ('NON' noprec) ;
  47. 'SI' ('NON' ('EXISTE' rvg 'preco')) ;
  48. rvg . 'preco' = 'TABLE' ;
  49. 'FINSI' ;
  50. rvgp = rvg . 'preco' ;
  51. * Précision relative utilisée pour le préconditionnement du calcul de la
  52. * matrice de pression
  53. precrel = 1.D-12 ;
  54. * precrel = 1.D-2 ;
  55. 'FINSI' ;
  56. * Projection incrémentale donc pas besoin de XINIT
  57. 'SI' ('EXIS' rvm 'XINIT' ) ;
  58. 'OUBLIER' rvm 'XINIT' ;
  59. 'FINSI' ;
  60. *fd 'OUBLIER' rvm 'XINIT' ;
  61. * On a mis chaine sinon plantage qd la pression s'appelle 'PRES' ;
  62. ngvit = 'CHAINE' rvg . 'NOMVIT' ;
  63. ngpre = 'CHAINE' rvg . 'NOMPRES' ;
  64. 'SI' ('NON' ('EXISTE' rvg 'METHINV')) ;
  65. rvg . 'METHINV' = rvm ;
  66. 'FINSI' ;
  67. rvgm = rvg . 'METHINV' ;
  68. *
  69. dim = 'VALEUR' 'DIME' ;
  70. *fd
  71. * Si l'utilisateur a transmis les noms de composantes on ne fait rien.
  72. * Sinon (il utilise EXEC) on construit la liste nivit = 'MOTS' '1UN' '2UN'
  73. 'SI' ('EXIS' rvg 'COMPVIT') ;
  74. nivit = rvg . 'COMPVIT' ;
  75. nivit2 = nivit ;
  76. nipre = rvg . 'COMPPRES' ;
  77. nipre2 = nipre ;
  78. 'SINON' ;
  79. nivit = 'MOTS' ;
  80. 'REPETER' idim dim ;
  81. niviti = 'CHAINE' &idim ngvit ;
  82. nivit = 'ET' nivit ('MOTS' niviti) ;
  83. 'FIN' idim ;
  84. nipre = 'MOTS' ngpre ;
  85. nivit2 = 'EXTRAIRE' ('MOTS' 'UX' 'UY' 'UZ') ('LECT' 1 'PAS' 1 dim) ;
  86. nipre2 = 'MOTS' 'SCAL' ;
  87. 'FINSI' ;
  88. * Construit la liste des noms d'inconnues ni vitesse ni pression
  89. nitout = 'EXTR' mat 'COMP' ;
  90. niaut = 'MOTS' ;
  91. 'REPETER' itout ('DIME' nitout) ;
  92. ninc = 'EXTRAIRE' nitout &itout ;
  93. lexis = FAUX ;
  94. 'REPETER' ivit ('DIME' nivit) ;
  95. lexis = 'OU' lexis ('EGA' ninc ('EXTRAIRE' nivit &ivit)) ;
  96. 'FIN' ivit ;
  97. 'REPETER' ipre ('DIME' nipre) ;
  98. lexis = 'OU' lexis ('EGA' ninc ('EXTRAIRE' nipre &ipre)) ;
  99. 'FIN' ipre ;
  100. 'SI' ('NON' lexis) ;
  101. niaut = 'ET' niaut ('MOTS' ninc) ;
  102. 'FINSI' ;
  103. 'FIN' itout ;
  104. *
  105. * Initialisation de la solution du système total
  106. *
  107. sol = 0. ;
  108. *
  109. * On résout d'abord la partie ni vitesse ni pression
  110. *
  111. 'SI' ('>' ('DIME' niaut) 0) ;
  112. mtot = 'EXTRAIRE' mat niaut niaut ;
  113. ftot = 'EXCO' niaut smbr niaut 'NOID' ;
  114. cltot = 'EXCO' niaut ccli niaut 'NOID' ;
  115. rvm . 'MATASS' =mtot ;
  116. rvm . 'MAPREC' =mtot ;
  117. ct = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvm ;
  118. sol = '+' ct sol ;
  119. 'FINSI' ;
  120. ********************************************************************
  121. * *
  122. * On résout la partie vitesse-pression par projection incrémentale *
  123. * *
  124. ********************************************************************
  125. nivp = 'ET' nivit nipre ;
  126. * Réduction des arguments de la proc à la vitesse-pression
  127. *mat = 'KOPS' 'EXTRINCO' mat nivp nivp ;
  128. mat = 'EXTRAIRE' mat nivp nivp ;
  129. ccli = 'EXCO' nivp ccli nivp 'NOID' ;
  130. smbr = 'EXCO' nivp smbr nivp 'NOID' ;
  131. *
  132. un = rv . 'INCO' . ngvit ;
  133. pn = rv . 'INCO' . ngpre ;
  134. un2 = 'NOMC' nivit2 nivit un ;
  135. pn2 = 'NOMC' nipre2 nipre pn ;
  136. incovp = un2 '+' pn2 ;
  137. vpetit = -1.D200 ;
  138. * Comme on va faire une méthode incrémentale, il faut imposer
  139. * les conditions aux limites de Dirichlet sur l'inconnue
  140. ccli1 = 'MASQUE' ccli 'SUPERIEUR' vpetit ;
  141. ccli0 = '-' ccli1 1.D0 ;
  142. niccli = 'EXTRAIRE' ccli 'COMP' ;
  143. dcli = '-' ccli ('*' incovp ccli1 niccli niccli niccli) ;
  144. incovp = '+' incovp dcli ;
  145. *
  146. * Matrice : partie diagonale en vitesse et contrainte
  147. * on est obligé de reconstruire la transposée car KBBT
  148. * ne stocke que B et pas Bt
  149. *
  150. mkvit = 'EXTRAIRE' mat nivit nivit ;
  151. mkcontr = 'EXTRAIRE' mat nivit nipre ;
  152. mkcontrt = 'KOPS' 'TRANSPOS' mkcontr ;
  153. Kt = mkvit 'ET' mkcontr 'ET' mkcontrt ;
  154. *
  155. * Construction du résidu * -1
  156. *
  157. residu = '-' smbr ('*' Kt incovp) ;
  158. *
  159. * Calcul de la matrice de pression et de la matrice masse diagonalisée
  160. *
  161. * On lumpe la partie diagonale et on l'inverse
  162. chpo1 = 'MASQUE' un2 'SUPERIEUR' vpetit ;
  163. mkvitd = '*' mkvit chpo1 ;
  164. *
  165. * Choix 1 : si la matrice masse lumpée
  166. * a tous ses termes positifs, on en prend l'inverse
  167. * Avantage : simple, autorise le préconditionnement
  168. * car ne change pas à chaque itération
  169. * Inconvénient : ne marche pas en axi quadratique sur l'axe
  170. * Choix 2 : on prend la diagonale de la matrice totale
  171. * on la met à l'échelle pour avoir la même masse que mkvitd
  172. * Avantage : marche en axi, sur maillage déformés
  173. * Inconvénient : change à chaque itération quand le problème
  174. * est non linéaire (=> recalcul de la matrice de pression),
  175. * pas forcément positif.
  176. * Choix 3 : valeur absolue du précédent + petite tolérance
  177. * Avantage : marche peut-être
  178. * Inconvénient : change à chaque itération.
  179. *
  180. *
  181. mmi = 'MINIMUM' mkvitd ;
  182. 'SI' ('>' mmi 1.D-200) ;
  183. imkvitd1 = 'INVERSE' mkvitd ;
  184. 'SI' debug ; 'MESSAGE' 'Choix 1 masse lumpée' ; 'FINSI' ;
  185. 'SINON' ;
  186. mkdia = 'EXTR' mkvit 'DIAG' ;
  187. mmi = 'MINIMUM' mkdia ;
  188. 'SI' ('>' mmi 1.D-200) ;
  189. 'SI' debug ; 'MESSAGE' 'Choix 2 masse lumpée' ; 'FINSI' ;
  190. 'SINON' ;
  191. mkdia = 'ABS' mkdia ;
  192. mma = 'MAXIMUM' mkdia ;
  193. mkdia = '+' mkdia ('*' mma 1.D-8) ;
  194. 'SI' debug ; 'MESSAGE' 'Choix 3 masse lumpée' ; 'FINSI' ;
  195. 'FINSI' ;
  196. masvitd = 'MAXIMUM' ('RESULT' mkvitd) ;
  197. * 'MESSAGE' ('CHAINE' 'totmas=' masvitd) ;
  198. masdia = 'MAXIMUM' ('RESULT' mkdia) ;
  199. imkvitd1 = '*' ('INVERSE' mkdia) ('/' masdia masvitd) ;
  200. 'FINSI' ;
  201. * On surcharge par 0. l'inverse de la diagonale là où il y a des clims
  202. * de Dirichlet
  203. imkvitd = '-' imkvitd1 ('*' imkvitd1 ccli1 niccli niccli niccli) ;
  204. *
  205. * Attention ! On préconditionne éventuellement le calcul de la matrice
  206. * de pression
  207. *
  208. precok = FAUX ;
  209. *'MESSAGE' ('CHAINE' 'precok = ' precok) ;
  210. 'SI' ('NON' noprec) ;
  211. exis1 = 'EXISTE' rvgp 'mklapphi' ;
  212. exis2 = 'EXISTE' rvgp 'imkvitd' ;
  213. 'SI' (exis1 'ET' exis2) ;
  214. imkvd = rvgp . 'imkvitd' ;
  215. mkl = rvgp . 'mklapphi' ;
  216. * On vérifie l'égalité des mkvitd à un facteur constant près
  217. ech = 'MAXIMUM' imkvitd 'ABS' ;
  218. echp = 'MAXIMUM' imkvd 'ABS' ;
  219. alfa = '/' ech echp ;
  220. * 'MESSAGE' ('CHAINE' 'alfa = ' alfa) ;
  221. dimkv = '-' imkvitd ('*' imkvd alfa) ;
  222. ndimkv = 'EXTRAIRE' dimkv 'COMP' ;
  223. inimkv = 'INVERSE' ('+' imkvitd ccli1) ;
  224. dimkvs = '*' dimkv inimkv ndimkv ndimkv ndimkv ;
  225. mdimkvs = 'MAXIMUM' dimkvs 'ABS' ;
  226. egaimkv = 'EGA' mdimkvs 0. precrel ;
  227. *
  228. 'SI' debug ;
  229. 'MESSAGE' ('CHAINE' 'egaimkv = ' egaimkv
  230. ' a ' ('*' mdimkvs 100.) ' %' ) ;
  231. 'FINSI' ;
  232. precok = egaimkv ;
  233. 'FINSI' ;
  234. 'FINSI' ;
  235. *
  236. chdbg = 'CHAINE' 'Matrice de pression MP : ' ;
  237. 'SI' precok ;
  238. 'SI' ('NON' ('EGA' alfa 1.D0 precrel)) ;
  239. chdbg = 'CHAINE' chdbg 'MP(i) = MP(i-1) * ' alfa ;
  240. * mklapphi = '*' (rvgp . 'mklapphi') alfa ;
  241. * imkvitd = '*' (rvgp . 'imkvitd') alfa ;
  242. mklapphi = rvgp . 'mklapphi' ;
  243. imkvitd = rvgp . 'imkvitd' ;
  244. 'SINON' ;
  245. chdbg = 'CHAINE' chdbg 'MP(i) = MP(i-1)' ;
  246. mklapphi = rvgp . 'mklapphi' ;
  247. imkvitd = rvgp . 'imkvitd' ;
  248. 'FINSI' ;
  249. 'SINON' ;
  250. * Dans le cas général où mkcontr et mkcontrt sont différents,
  251. * il faudrait mettre :
  252. *mklapphi = 'KOPS' 'CMCT' mkcontr imchd (kops mkcontrt 'TRANSPOS') ;
  253. chdbg = 'CHAINE' chdbg 'calcul...' ;
  254. mklapphi = 'KOPS' 'CMCT' mkcontr imkvitd mkcontr ;
  255. alfa = 1.D0 ;
  256. 'FINSI' ;
  257. 'SI' debug ; 'MESSAGE' chdbg ; 'FINSI' ;
  258. 'SI' ('NON' noprec) ;
  259. rvgp . 'mklapphi' = mklapphi ;
  260. rvgp . 'imkvitd' = imkvitd ;
  261. 'FINSI' ;
  262. *'MESSAGE' 'Fin du calcul de la matrice de pression' ;
  263. *
  264. * Calcul éventuel d'une estimation de la pression (double projection)
  265. *
  266. *'MESSAGE' ('CHAINE' 'dblproj=' dblproj) ;
  267. 'SI' dblproj ;
  268. desqdm = 'EXCO' nivit residu nivit 'NOID' ;
  269. idesq = '*' desqdm imkvitd nivit nivit nivit ;
  270. * 1) Celui-ci se compense avec le 2)
  271. * donc pas nécessaire
  272. * 'SI' ('NEG' alfa 1.D0 1.D-13) ;
  273. * idesq = '*' idesq alfa ;
  274. * 'FINSI' ;
  275. didesq = '*' mkcontr idesq ;
  276. mtot = mklapphi ;
  277. *dbg dmtot = 'EXTRAIRE' mklapphi 'DIAG' ;
  278. *dbg HCHPO dmtot 'dmtot' ;
  279. ftot = didesq ;
  280. cltot = 'EXCO' nipre ccli0 nipre 'NOID' ;
  281. rvgm .'MATASS' =mtot ;
  282. rvgm . 'MAPREC' =mtot ;
  283. *dbg dphi1 tt = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvgm
  284. *dbg 'IMPR' 3 'LTIME' VRAI ;
  285. dphi1 = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvgm ;
  286. * 2) Se compense avec le 1)
  287. * donc pas nécessaire
  288. * 'SI' ('NEG' alfa 1.D0 1.D-13) ;
  289. * dphi1 = '/' dphi1 alfa ;
  290. * 'FINSI' ;
  291. * HCHPO dphi1 'dphi1' ;
  292. gdphi1 = '*' mkcontrt dphi1 ;
  293. residu = residu '-' gdphi1 ;
  294. 'FINSI' ;
  295. *
  296. * Calcul de l'incrément de vitesse dv* à partir de la qdm
  297. * On splitte les équations par direction
  298. *
  299. dvitstar = 0. ;
  300. ncomp = 'DIME' nivit ;
  301. 'REPETER' iicomp ncomp ;
  302. icomp = &iicomp ;
  303. licomp = 'LECT' icomp ;
  304. mcomp = 'EXTRAIRE' nivit licomp ;
  305. mtot = 'EXTRAIRE' mkvit mcomp mcomp ;
  306. ftot = 'EXCO' mcomp residu mcomp 'NOID' ;
  307. cltot = 'EXCO' mcomp ccli0 mcomp 'NOID' ;
  308. rvm . 'MATASS' =mtot ;
  309. rvm . 'MAPREC' =mtot ;
  310. *dbg ct tt = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvm
  311. *dbg 'IMPR' 3 'LTIME' VRAI ;
  312. ct = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvm ;
  313. dvitstar = '+' dvitstar ct ;
  314. 'FIN' iicomp ;
  315. * HCHPO dvitstar 'dvitstar' ;
  316. *
  317. * Calcul de l'incrément de pression dphi pour vérifier la conservation
  318. * de la masse
  319. *
  320. * Pas besoin de mettre dvitstar à 0 sur les clims Dirichlet car c'est
  321. * déjà le cas !
  322. smb1 = '*' mkcontr dvitstar ;
  323. mtot = mklapphi ;
  324. ftot = '-' smb1 ('EXCO' nipre residu nipre 'NOID') ;
  325. cltot = 'EXCO' nipre ccli0 nipre 'NOID' ;
  326. rvgm .'MATASS' =mtot ;
  327. rvgm . 'MAPREC' =mtot ;
  328. dphi = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvgm ;
  329. 'SI' ('NEG' alfa 1.D0 1.D-13) ;
  330. dphi = '/' dphi alfa ;
  331. 'FINSI' ;
  332. *HCHPO dphi 'dphi' ;
  333. *
  334. * Calcul de l'incrément final de vitesse dvit
  335. *
  336. gdphi = '*' mkcontrt dphi ;
  337. * On surcharge par 0. là où il y a des clims de Dirichlet
  338. * Pas besoin !!! vu que imkvitd est là !
  339. *!! gdphi = '-' gdphi ('*' gdphi ccli1 niccli niccli niccli) ;
  340. gdphi = '*' gdphi imkvitd nivit nivit nivit ;
  341. 'SI' ('NEG' alfa 1.D0 1.D-13) ;
  342. gdphi = '*' gdphi alfa ;
  343. 'FINSI' ;
  344. dvit = dvitstar '-' gdphi ;
  345. *HCHPO dvit 'dvit' ;
  346. *
  347. 'SI' dblproj ;
  348. dphi = dphi '+' dphi1 ;
  349. 'FINSI' ;
  350. incrfin = dvit '+' dphi ;
  351. solvp = incovp '+' incrfin ;
  352. *
  353. sol = '+' solvp sol ;
  354. *
  355. 'RESPRO' sol ;
  356. *
  357. * End of procedure file GRESP
  358. *
  359. 'FINPROC' ;
  360.  

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