Télécharger solvEFMH.procedur

Retour à la liste

Numérotation des lignes :

  1. * SOLVEFMH PROCEDUR GOUNAND 11/05/24 21:16:18 6976
  2. **********************************************************************
  3. 'DEBP' SOLVEFMH MoDARCY*'MMODEL' ChPSour*'CHPOINT'
  4. MassEFMH*'RIGIDITE' MatTR/'RIGIDITE' matTM/'MATRIK'
  5. SMTR*'CHPOINT' Tcini*'CHPOINT' cini*'CHPOINT'
  6. QFACE/'CHPOINT' nomespec*'LISTMOTS'
  7. nbespece*'ENTIER' nbsource*'ENTIER' TABRES*'TABLE'
  8. tbdartra*'TABLE' CHCLIM*'TABLE';
  9. * |-----------------------------------------------------------------|
  10. * | Phrase d'appel (en GIBIANE) |
  11. * |-----------------------------------------------------------------|
  12. * | |
  13. * | mattr TABRES tfin tcfin = SOLVEFMH MoDARCY ChPSour |
  14. * | MassEFMH (MatTR ou mattM) SMTR Tcini cini |
  15. * | nomespec nbespece nbsource TABRES tbdartra |
  16. * | CHCLIM ; |
  17. * | |
  18. * | |
  19. * |-----------------------------------------------------------------|
  20. * | Généralités : MATTEFMH construit la matrice de discrétisation |
  21. * | du problème de transport convection-diffusion pour|
  22. * | le premier pas de tps d'un algorithme transitoire.|
  23. * | Le second membre et les Conditions limites de flux|
  24. * | sont pris en compte. |
  25. * | RESTE TCINI, DECENTR et TERME LIN |
  26. * |-----------------------------------------------------------------|
  27. * | |
  28. * |-----------------------------------------------------------------|
  29. * | ENTREES |
  30. * |-----------------------------------------------------------------|
  31. * | MoDARCY : modele Darcy. |
  32. * | |
  33. * | Deltat : pas de temps utilisé pour calculer la concentration |
  34. * | |
  35. * | ChPSour : Champ par points des sources volumiques par unité de |
  36. * | temps (support maillage centre). Composante associées|
  37. * | aux especes |
  38. * | |
  39. * | MassEFMH : matrice elementaire EFMH |
  40. * | |
  41. * | MatTr : matrice globale sur les traces - rigidité |
  42. * | argument optionnel on donne alors martTM |
  43. * | |
  44. * | MatTM : matrice globale sur les traces - Matrik |
  45. * | argument optionnel on donne alors matrTR |
  46. * | |
  47. * | SMTr : second membre sur les traces |
  48. * | |
  49. * | Tcini : Trace de concentration aux faces (eventuellement à |
  50. * | plusieurs composantes (espèces) - sert a initialiser |
  51. * | le XINIT de KRES et à calculer la valeur de la |
  52. * | concentration au centre |
  53. * | |
  54. * | nomespec : liste des noms de composante des espèces dans Cini |
  55. * | |
  56. * | nbespece : nombre de composante de Cini, soit nombre d'especes |
  57. * | |
  58. * | nbsource : nombre de composantes du terme source qd X especes |
  59. * | |
  60. * | TABRES : Table complète définissant les options de résolution |
  61. * | pour 'KRES'. |
  62. * | |
  63. * | TbDarTra : table Darcy transitoire utilisée par MHYB, SMTP ... |
  64. * | |
  65. * | CHCLIM : table d'indice 'NEUMANN' et 'DIRICHLET' contenant les|
  66. * | Chpoint à n composantes contenant les conditions aux |
  67. * | limites de Neumann et Dirichlet par espece. |
  68. * | |
  69. * |-----------------------------------------------------------------|
  70. * | SORTIES |
  71. * |-----------------------------------------------------------------|
  72. * | |
  73. * | Tcfin : Trace de concentration aux faces (eventuellement à |
  74. * | plusieurs composantes (espèces) - état final apres |
  75. * | résolution |
  76. * | |
  77. * | cfin : concentration apres calcul pour toutes les especes |
  78. * | |
  79. * | TABSORT : Table complète définissant les options de résolution |
  80. * | pour 'KRES'. |
  81. * | |
  82. * | Matk : matrice globale sur les traces pour la convection |
  83. * | en format matrik. Elle diffère de la matrice entree |
  84. * | si cette derniere est une rigidité car traduite en |
  85. * | Matrik. Elle contient également les préconditionnemen|
  86. * | crée par l'opérateur de résolution KRES |
  87. * | |
  88. * |-----------------------------------------------------------------|
  89. * | VARIABLES INTERNES |
  90. * |-----------------------------------------------------------------|
  91. * | |
  92. * | Tccfin : Trace de concentration aux faces (une composante) |
  93. * | |
  94. * | CChpsour : Source aux centre (une composante) |
  95. * | |
  96. * | CCfin : concentration aux centres (une composante) |
  97. * | |
  98. * | SSMTr : second membre sur les traces pour une espèce |
  99. * | |
  100. * | Nouvmatr : Logique, si VRAI on transforme la matrice rigidité |
  101. * | en matrik |
  102. * | |
  103. * | DIRCLI : logique valant VRAI si conditions aux |
  104. * | limites de Dirichlet |
  105. * | |
  106. * | CLDIRI : Chpoint à n composantes contenant les conditions aux |
  107. * | limites de Dirichlet par espece. |
  108. * | il faudra en faire un nuage si supports géométriques |
  109. * | différents par espece. OPTIONNEL |
  110. * | |
  111. **********************************************************************
  112.  
  113.  
  114.  
  115.  
  116. *---------------------------------------------------------------------
  117. *---------- On récupere les conditions limites ------------------
  118. *---------------------------------------------------------------------
  119.  
  120. * Flag sur conditions limites initialisés
  121. DIRCLI = FAUX ;
  122. FTOCLI = FAUX ;
  123. MIXCLI = FAUX ;
  124.  
  125. * Dirichlet
  126. 'SI' ('EXISTE' CHCLIM 'DIRICHLET') ;
  127. CLDIRI = CHCLIM . 'DIRICHLET' ;
  128. DIRCLI = VRAI ;
  129. 'FINSI' ;
  130.  
  131. * Conditions flux total
  132. 'SI' ('EXISTE' CHCLIM 'FLUTOTAL') ;
  133. FTOCLI = VRAI ;
  134. 'FINSI' ;
  135.  
  136. * Conditions flux mixte
  137. 'SI' ('EXISTE' CHCLIM 'FLUMIXTE') ;
  138. MIXCLI = VRAI ;
  139. 'FINSI' ;
  140.  
  141.  
  142. *
  143. * On recopie la table de résolution TABRES dans TABSORT
  144. * Attention si une valeur contenue dans la table a le
  145. * debut d'un nom d'opérateur de castem, il y a probleme
  146. * d'ou démarrage apres l'indice soustype de valeur
  147. * METHINV, identique à opérateur METHode.
  148. *
  149.  
  150. dumm = 'INDEX' TABRES ;
  151. TABSORT = 'TABLE' METHINV ;
  152. 'REPETER' bou1 (('DIME' dumm) '-' 1) ;
  153. TABSORT . (mot dumm . (&bou1 '+' 1)) = TABRES . (mot dumm . (&bou1
  154. '+' 1)) ;
  155. 'FIN' bou1 ;
  156.  
  157. * idem avec tbdartra pour la preserver
  158.  
  159. dumm = 'INDEX' TbDartra ;
  160. TbDartrb = 'TABLE' 'DARCY_TRANSITOIRE' ;
  161. 'REPETER' bou1 (('DIME' dumm) '-' 1) ;
  162. TbDartrb . (mot dumm . (&bou1 '+' 1)) = TbDartra . (mot dumm . (&bou1
  163. '+' 1)) ;
  164. 'FIN' bou1 ;
  165.  
  166.  
  167.  
  168.  
  169. *---------------------------------------------------------------------
  170. *-------------- TRANSFORMATION MATRIK RIGIDITE -----------------------
  171. *---------------------------------------------------------------------
  172.  
  173. 'SI' ('EXISTE' mattr) ;
  174. * matrice rigidité donc vient d'etre calculée
  175. * il faudra la traduire en matrik.
  176. NOUVMATR = VRAI ;
  177. 'SINON' ;
  178. * matrice issue d'un précédent calcul
  179. * on la stocke dans mattr
  180. NOUVMATR = FAUX ;
  181. mattr = mattm ;
  182. 'FINSI' ;
  183.  
  184. 'SI' (NOUVMATR) ;
  185. * On recalcule les conditions aux limites flux total si existent
  186. 'SI' (FTOCLI) ;
  187. 'SI' ('EXISTE' QFACE) ;
  188. 'SI' ('EGA' ('TYPE' CHCLIM . 'FLUTOTAL') 'CHPOINT') ;
  189. mayage = 'EXTRAIRE' CHCLIM . 'FLUTOTAL' maillage ;
  190. vites = 'REDU' QFACE mayage ;
  191. matcli = 'KOPS' 'MATDIAGO' ('NOMC' 'TH' vites) 'MATRIK' ;
  192. 'FINSI' ;
  193. 'FINSI' ;
  194. 'FINSI' ;
  195. * On recalcule les conditions aux limites flux mixte si existent
  196. 'SI' (MIXCLI) ;
  197. mayage = 'EXTRAIRE' CHCLIM . 'FLUMIXTE' . 'VAL' maillage ;
  198. cofb = (doma modarcy SURFACE) * CHCLIM . 'FLUMIXTE' . 'COEFB' ;
  199. cofa = 'REDU' (CHCLIM . 'FLUMIXTE' . 'COEFA') mayage ;
  200. cofb = 'REDU' cofb mayage ;
  201. coefm = (-1.D0 * cofb) '/' cofa ;
  202. matmix = 'KOPS' 'MATDIAGO' ('NOMC' 'TH' coefm) 'MATRIK' ;
  203. 'OUBLIER' cofa ;
  204. 'OUBLIER' cofb ;
  205. 'OUBLIER' coefm ;
  206. 'FINSI' ;
  207.  
  208. * Nouvelle matrice rigidité, on la transforme en Matrik
  209. mattr = 'KOPS' 'RIMA' mattr ;
  210. mattr = KOPS 'CHANINCO' mattr
  211. ('MOTS' 'TH' 'LX') ('MOTS' 'TH' 'LX')
  212. ('MOTS' 'FLUX' 'FLX') ('MOTS' 'TH' 'LX') ;
  213. * On rajoute les conditions aux limites de flux totale (matrice U tC)
  214. 'SI' ('EGA' ('TYPE' matcli) 'MATRIK') ;
  215. mattr = mattr 'ET' matcli ;
  216. 'FINSI' ;
  217. * On rajoute les conditions aux limites de flux mixte
  218. 'SI' ('EGA' ('TYPE' matmix) 'MATRIK') ;
  219. mattr = mattr 'ET' matmix ;
  220. 'FINSI' ;
  221. * On initialise les indices pointants sur les matrices
  222. * assembmées et préconditionnées
  223. TABSORT . 'MATASS' = Mattr ;
  224. TABSORT . 'MAPREC' = Mattr ;
  225. 'FINSI' ;
  226.  
  227.  
  228. *---------------------------------------------------------------------
  229. *-------------- RESOLUTION EN TRACE DE CONCENTRATION -----------------
  230. *---------------------------------------------------------------------
  231.  
  232. * boucle sur les espèces.
  233. 'REPETER' bloc1 nbespece ;
  234. * préparation solution initiale - ie trace initiale
  235. TABSORT . 'XINIT' = 'NOMC' 'TH' ('EXCO' ('EXTRAIRE' &bloc1 nomespec)
  236. Tcini) ;
  237. tccini = TABSORT . 'XINIT' ;
  238. * préparation second membre
  239. SSmtr = 'NOMC' 'TH' ('EXCO' ('EXTRAIRE' &bloc1 nomespec) Smtr) ;
  240.  
  241.  
  242. * Solution en trace
  243. * Si conditions de Dirichlet
  244. 'SI' (DIRCLI) ;
  245. Tccfin = KRES mattr 'TYPI' TABSORT
  246. 'SMBR' SSMTr
  247. 'CLIM' ('NOMC' 'TH' ('EXCO' ('EXTRAIRE' &bloc1 nomespec)
  248. CLDIRI))
  249. 'IMPR' 0 ;
  250. 'SINON' ;
  251. Tccfin = KRES mattr 'TYPI' TABSORT
  252. 'SMBR' SSMTr
  253. 'IMPR' 0 ;
  254. 'FINSI' ;
  255. TABSORT . 'MATASS' = mattr ;
  256. TABSORT . 'MAPREC' = mattr ;
  257.  
  258. * préparation du terme source
  259. 'SI' (nbsource 'EGA' 1) ;
  260. CChpsour = 'NOMC' 'SOUR' Chpsour ;
  261. 'SINON' ;
  262. CChpsour = 'NOMC' 'SOUR' ('EXCO' ('EXTRAIRE' &bloc1 nomespec)
  263. Chpsour) ;
  264. 'FINSI' ;
  265.  
  266. TbDartrb . 'CHARGE' = 'NOMC' 'H' ('EXCO' ('EXTRAIRE' &bloc1 nomespec)
  267. cini) ;
  268. TbDartrb . 'TRACE' = TABSORT . 'XINIT' ;
  269. *
  270. 'SI' ('EXISTE' QFACE) ;
  271. TETAC = TbDartra . 'THETA_CONVECTION' ;
  272. * On ajoute au terme source total la contribution de la convection
  273. CONC1 = ((1.D0 - TETAC) * tccini) + (TETAC * tccfin) ;
  274. CONC1 = (-1.D0) * ('NOMC' 'SCAL' CONC1) * QFACE ;
  275. CONC1 = 'DIVU' MoDarcy CONC1 ('DOMA' MoDarcy 'ORIENTAT') ;
  276. CONC1 = 'NOMC' 'SOUR' CONC1 ;
  277. CChpsour = CChpsour + CONC1 ;
  278. oubli CONC1 ;
  279. 'FINSI' ;
  280.  
  281. * on reconstruit la concentration au centre
  282. CCfin = 'HYBP' MoDARCY MassEFMH Tccfin TbDartrb CChpsour ;
  283. * On reconstruit le flux diffusif aux faces
  284. CCFLU = 'HDEB' MoDARCY MassEFMH CCfin Tccfin ;
  285. 'SI' ('EXISTE' QFACE) ;
  286. * On construit le flux convectif aux faces
  287. * CCFLUCO = ((1.D0 - TETAC) * tccini) + (TETAC * tccfin) ;
  288. * VERRUE
  289. CCFLUCO = tccfin ;
  290. CCFLUCO = ('NOMC' 'SCAL' CCFLUCO) * QFACE ;
  291. 'SINON' ;
  292. * flux convectif nul
  293. CCFLUCO = 0.D0 * CCFLU ;
  294. 'FINSI' ;
  295.  
  296. * On reconstitue les champoints à plusieurs composante
  297.  
  298. 'SI' (&bloc1 'EGA' 1) ;
  299. CCfin = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) CCfin ;
  300. Tccfin = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) TCCfin ;
  301. CCflu = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) CCFLU ;
  302. CCFLUCO = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) CCFLUCO
  303. NATU DISCRET ;
  304. Tcfin = tccfin ;
  305. cfin = ccfin ;
  306. cflu = ccflu ;
  307. cfluco = ccfluco ;
  308. 'SINON' ;
  309. * CCfin = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) ('COPIER' CCfin) ;
  310. * Tccfin = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) ('COPIER' TCCfin) ;
  311. * ccflu = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) ('COPIER' CCflu) ;
  312. * CCFLUCO = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) ('COPIER' CCFLUCO)
  313. * NATU DISCRET ;
  314. CCfin = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) CCfin ;
  315. Tccfin = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) TCCfin ;
  316. ccflu = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) CCflu ;
  317. CCFLUCO = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) CCFLUCO
  318. NATU DISCRET ;
  319. Tcfin = tcfin 'ET' tccfin ;
  320. cfin = cfin 'ET' ccfin ;
  321. cflu = cflu 'ET' ccflu ;
  322. cfluco = cfluco 'ET' ccfluco ;
  323. 'FINSI' ;
  324.  
  325. 'FIN' bloc1 ;
  326.  
  327. 'FINP' mattr TABSORT cfin tcfin cflu cfluco ;
  328.  
  329.  
  330.  
  331.  
  332.  

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