Télécharger updaVF.procedur

Retour à la liste

Numérotation des lignes :

  1. * UPDAVF PROCEDUR GOUNAND 11/05/24 21:16:24 6976
  2. **********************************************************************
  3. 'DEBP' UPDAVF MoDARCY*'MMODEL' Porosite*'CHPOINT' Matediff*'CHPOINT'
  4. Diffdisp*'CHPOINT' ChPSour*'CHPOINT' Deltat*'FLOTTANT'
  5. Cini*'CHPOINT' TetaDiff*'FLOTTANT' TetaConv*'FLOTTANT'
  6. Qface/'CHPOINT' nomespec*'LISTMOTS' nbespece*'ENTIER'
  7. nbsource*'ENTIER' mattr*'MATRIK' JACO*'MATRIK'
  8. MPOR2*'MATRIK' MCHAMT*'MCHAML' MCHAMT1*'MCHAML'
  9. TABMODI*'TABLE' CHCLIM*'TABLE' LCONV*'LOGIQUE';
  10. * ATTENTION La vitesse est optionnelle, L'ordre est important
  11. * et les types d'arguments qui se suivent aussi pour tester leur
  12. * présence
  13. *
  14. * Attention il faudra transformer les vitesses en débits face
  15. * et sortir le champ
  16. *
  17. * |-----------------------------------------------------------------|
  18. * | Phrase d'appel (en GIBIANE) |
  19. * |-----------------------------------------------------------------|
  20. * | |
  21. * SMTr Mattt Difftot Mctot Mdiff Nouvmat = UPDAVF
  22. * MoDARCY Porosite Matediff Diffdisp ChPSour
  23. * DeltaT Cini TetaDiff TetaConv
  24. * Qface nomespec nbespece nbsource
  25. * Matot Jaco Mctot Mdiff Mpor TABMODI CHCLIM ;
  26. * |-----------------------------------------------------------------|
  27. * | Généralités : MATTVF construit la matrice de discrétisation |
  28. * | du problème de transport convection-diffusion pour|
  29. * | le premier pas de tps d'un algorithme transitoire.|
  30. * | Le second membre et les Conditions limites de flux|
  31. * | sont pris en compte. |
  32. * | RESTE TCINI, DECENTR et TERME LIN |
  33. * |-----------------------------------------------------------------|
  34. * | |
  35. * |-----------------------------------------------------------------|
  36. * | ENTREES |
  37. * |-----------------------------------------------------------------|
  38. * | MoDARCY : modele Darcy. |
  39. * | |
  40. * | Porosite : champ par elements de composante 'CK' |
  41. * | |
  42. * | MateDiff : Tenseur de diffusion (type iso, ..) champ par |
  43. * | points de composante 'K' en isotrope, 'K11', 'K21',|
  44. * | 'K22' en anisotrope 2d et 'K11', 'K21', 'K22', 'K31'|
  45. * | 'K32', 'K33' en anisotrope 3d. Type 'CARACTERISTIQUE'|
  46. * | |
  47. * | Diffdisp : Tenseur de dispersion (type iso, ..) champ par |
  48. * | points de composante 'K' en isotrope, 'K11', 'K21',|
  49. * | 'K22' en anisotrope 2d et 'K11', 'K21', 'K22', 'K31'|
  50. * | 'K32', 'K33' en anisotrope 3d. Type 'CARACTERISTIQUE'|
  51. * | |
  52. * | ChPSour : Champ par points des sources volumiques par unité de |
  53. * | temps (support maillage centre). Composante ?????? |
  54. * | |
  55. * | Cini : Concentration initiale, CHPOINT centre. |
  56. * | Composante 'H'. |
  57. * | |
  58. * | Deltat : Pas de temps |
  59. * | |
  60. * | Qface : vitesse aux faces, CHPO face de composantes Vx, Vy |
  61. * | en 2d et Vx, Vy, Vz en 3d. Il s'agit plus exatement |
  62. * | de (V.n)n, c'est à dire de la composante normale de |
  63. * | la vitesse aux faces. ???????? (je pressens que |
  64. * | castem va sortir des flux, cad intégrés sur surfaces)|
  65. * | |
  66. * | nomespec : liste des noms de composante des espèces dans Cini |
  67. * | |
  68. * | nbespece : nombre de composante de Cini, soit nombre d'especes |
  69. * | |
  70. * | nbsource : nombre de composantes du terme source qd X especes |
  71. * | |
  72. * | Matot : matrice globale de discretisation en VF |
  73. * | |
  74. * | Jaco : matrice globale de discretisation en VF pour le problème
  75. * | stationnaire
  76. *
  77. * | Mpor : matrice globale de discretisation en VF pour le problème
  78. * | stationnaire
  79. * | |
  80. * | Mchamt : Coef permettant de calculer le flux total
  81. * | |
  82. * | Mchamt1 : Coef permettant de calculer le flux diffusif
  83. * | |
  84. * | |
  85. * | TABMODI : table contenant des logiques indiquant la nécessité |
  86. * | ou non de reclalculer certains termes. |
  87. * | 'POROSITE' : VRAI si le coefficient devant D/DT |
  88. * | (porosité) est modifié depuis le dernier|
  89. * | appel |
  90. * | 'DELTAT' : VRAI si le pas de tps a changé |
  91. * | 'CONVECTI' : VRAI si la vitesse a changé |
  92. * | 'COEF_LIN' : VRAI si le coef en facteur de C a changé|
  93. * | 'DIFFUSI' : VRAI si les diffusivités ont changé |
  94. * | |
  95. * | CHCLIM : table d'indice 'NEUMANN' et 'DIRICHLET' contenant les|
  96. * | Chpoint à n composantes contenant les conditions aux |
  97. * | limites de Neumann et Dirichlet par espece. |
  98. * | |
  99. * | |
  100. * |-----------------------------------------------------------------|
  101. * | ENTREES-SORTIES |
  102. * |-----------------------------------------------------------------|
  103. * | |
  104. * | Difftot : Coefficient de diffusion totale, integre decentrement|
  105. * | |
  106. * | |
  107. * |-----------------------------------------------------------------|
  108. * | SORTIES |
  109. * |-----------------------------------------------------------------|
  110. * | |
  111. * | |
  112. * | RESI : second membre |
  113. * | |
  114. * | Matot : matrice globale de discretisation en VF |
  115. * | |
  116. * | Difftot : Coefficient de diffusion totale, integre decentrement|
  117. *
  118. * | Mpor : matrice globale de discretisation en VF pour le problème
  119. * | stationnaire
  120. * | |
  121. * | Mchamt : Coef permettant de calculer le flux total
  122. * | |
  123. * | Mchamt1 : Coef permettant de calculer le flux diffusif
  124. * | |
  125. * | |
  126. * |-----------------------------------------------------------------|
  127. * | VARIABLES INTERNES |
  128. * |-----------------------------------------------------------------|
  129. * | |
  130. * | |
  131. * | PCONV : Logique indiquant VRAI si présence de convection |
  132. * | |
  133. * | CoefDt : coeff devant dC/dt integre sur les elements |
  134. * | |
  135. * | CCini : concentration aux centres (une composante) |
  136. * | |
  137. * | SSource : Source aux centre (une composante) |
  138. * | |
  139. * | SSMTr : second membre sur les traces pour une espèce |
  140. * | |
  141. * | toltheta : 1.D-4 seuil en dessous duquel on considère que la |
  142. * | valeur de theta du theta-schéma est nulle (schéma |
  143. * | explicite) ou au contraire euler-implicite si |
  144. * | theta > 0.9999 |
  145. * | |
  146. * | tetaconv : valeur du coef d'implicitation pour la convection |
  147. * | |
  148. * | MODDIFF : Logique, si VRAI on recalcule les termes liés à la |
  149. * | diffusion, en supposant qu'elle a changé, donnés par |
  150. * | Matediff. En implicite cela affecte la matrice |
  151. * | rigidité du problème. |
  152. * | |
  153. * | MODPOROS : Si VRAI on recalcule la contribution de la porosité |
  154. * | pour la matrice du problème et le coeficient devant |
  155. * | dC/dt (VF) intégré sur les mailles. la matrice |
  156. * | Identité des VF est recalculée. Logique. |
  157. * | |
  158. * | MODELTT : Logique. Si vrai on recalcule la matrice rigidité du |
  159. * | problème. Le pas de temps a changé. |
  160. * | |
  161. * | MODCONV : Logique. Si vrai on recalcule la matrice rigidité du |
  162. * | problème. La vitesse a changé. |
  163. * | |
  164. * | MODCOLIN : Logique. Si vrai on recalcule la matrice rigidité du |
  165. * | problème. Le coef linéaire devant C a changé. |
  166. * | |
  167. * | FLUNEU : LOGIQUE valant VRAI si conditions de Neumann |
  168. * | |
  169. * | CLFLUX : Chpoint à n composantes contenant les flux imposés |
  170. * | pour chaque espece chimique. nul si pas de flux |
  171. * | OPTIONNEL |
  172. * | |
  173. **********************************************************************
  174.  
  175.  
  176. *---------------------------------------------------------------------
  177. *---------- On récupere les conditions limites ------------------
  178. *---------------------------------------------------------------------
  179.  
  180. FLUNEU = FAUX;
  181. DIRCLI = FAUX;
  182. FLUTOT = FAUX;
  183. FLUMIX = FAUX;
  184. CLFLUX MA = 'KOPS' 'MATRIK';
  185. CLDIRI MA = 'KOPS' 'MATRIK';
  186. CLFLUT MA = 'KOPS' 'MATRIK';
  187. EXFLU = CLFLUX;
  188. EXDIR = CLDIRI;
  189. EXFLUT = CLFLUT;
  190. * Neumann
  191. 'SI' ('EXISTE' CHCLIM 'NEUMANN') ;
  192. CLFLUX = CHCLIM . 'NEUMANN';
  193. FLUNEU = VRAI;
  194. 'FINSI';
  195.  
  196. 'SI' ('EXISTE' CHCLIM 'DIRICHLET') ;
  197. CLDIRI = CHCLIM . 'DIRICHLET';
  198. DIRCLI = VRAI;
  199. 'FINSI';
  200.  
  201. 'SI' ('EXISTE' CHCLIM 'FLUTOTAL') ;
  202. CLFLT = CHCLIM . 'FLUTOTAL';
  203. FLUTOT = VRAI;
  204. 'FINSI';
  205.  
  206. * Flux mixte
  207. 'SI' ('EXISTE' CHCLIM 'FLUMIXTE') ;
  208. * comme on impose A Dgrad C + B C = flumix, on le traite sous
  209. * la forme D grad C + (B/A) C = flumix/A plus naturelle en EFMH car
  210. * D grad C est le flux diffusif
  211. COFA = -1.D0 * CHCLIM . 'FLUMIXTE' . 'COEFA' ;
  212. CLFLUX3 = CHCLIM . 'FLUMIXTE' . 'VAL' '/' COFA ;
  213. CLFLUX3 = CHAN 'ATTRIBUT' CLFLUX3 NATURE DISCRET ;
  214. FLUMIX = VRAI ;
  215. mayage = 'EXTRAIRE' CHCLIM . 'FLUMIXTE' . 'VAL' maillage ;
  216. cofb = (doma modarcy SURFACE) *CHCLIM . 'FLUMIXTE' . 'COEFB' ;
  217. cofb = CHCLIM . 'FLUMIXTE' . 'COEFB' ;
  218. cofa = 'REDU' (CHCLIM . 'FLUMIXTE' . 'COEFA') mayage ;
  219. cofb = 'REDU' cofb mayage ;
  220. coefm = (1.D0 * cofb) '/' cofa ;
  221. CLFLUX3 = 'REDU' CLFLUX3 mayage ;
  222. 'OUBLIER' cofa ;
  223. 'OUBLIER' cofb ;
  224. 'FINSI' ;
  225.  
  226.  
  227. *---------------------------------------------------------------------
  228. *---------- Initialisations de tables, coefficients ------------------
  229. *---------------------------------------------------------------------
  230.  
  231. toltheta = 1.D-4 ;
  232.  
  233. MODPOROS = TABMODI . 'POROSITE' ;
  234. MODCONV = TABMODI . 'CONVECTI' ;
  235. MODDELTT = TABMODI . 'DELTAT' ;
  236. MODCOLIN = TABMODI . 'COEF_LIN' ;
  237. MODDIFF = TABMODI . 'DIFFUSIV' ;
  238.  
  239.  
  240.  
  241.  
  242. * Si la porosité est modifiée (ou plus exactement le coef
  243. * devant dC/dt)
  244.  
  245.  
  246. 'SI' (MODPOROS) ;
  247. MPOR2 = 'KOPS' 'MATDIAGO' ('NOMC' 'RETN' Porosite) 'MATRIK' ;
  248. 'FINSI' ;
  249.  
  250. *
  251. * Si le pas de temps est modifié
  252. *
  253.  
  254.  
  255. *
  256. * Si le coefficient de décroissance (ou linaire) en facteur de C est
  257. * modifié ????????????????,
  258. *
  259.  
  260. 'SI' (MODCOLIN) ;
  261. * A priori sera géré dans le coeff de porosité
  262. 'FINSI' ;
  263.  
  264.  
  265.  
  266.  
  267.  
  268. *----------------------------------------------------------------------
  269. *---------- CALCUL DE LA MATRICE DE DIFFUSION DU PROBL ----------------
  270. *----------------------------------------------------------------------
  271.  
  272.  
  273. * Calcul de la matrice du probleme diffusion transitoire
  274. * Elle est modifiée si la diffusion ,
  275. * ou si le pas de temps change, ou si le coefficient
  276. * en facteur de la dérivée temporelle change (porosité*retard),
  277. * ou si le coefficient du terme linéaire a(x,t) * C(x,t) change.
  278.  
  279. * Pour des raisons pratiques on n'a pas stocké les matrices
  280. * de convection et de diffusion séparément. Si l'une est modifiée
  281. * on doit donc recalculer les deux matrices. Ce choix est discutable,
  282. * on gagne un facteur 3 en stockage (mat totale au lieu de mat totale
  283. * + diffusion + convection). Par ailleurs, le cout de l'assemblage +
  284. * resolution + créer une des deux matrice rend souvent négligeable
  285. * le recalcul inutile d'une des deux matrices précitées. On réduit
  286. * également le nombre d'entrées-sorties.
  287.  
  288.  
  289. 'SI' (MODDIFF) ;
  290. * La diffusivité effective a changé, on remmet la matrice au propre
  291. * c'est à dire avec des noms de composantes isotropes K11 K22 etc ...
  292. zozo = DIFFANIS Matediff ;
  293. Matediff = zozo ;
  294. 'FINSI' ;
  295.  
  296. Difftot = MateDiff '+' Diffdisp ;
  297.  
  298.  
  299. *----------------------------------------------------------------------
  300. *--- CALCUL DE LA MATRICE DE CONVECTION ET DES SECONDS MEMBRES --------
  301. *----------------------------------------------------------------------
  302.  
  303.  
  304. * Nombres d'espèces à gérer est fixé et ne peut changer au cours d'une
  305. * simulation. Il sera d'ailleurs mis dans le modele
  306.  
  307. 'SI' (nbespece > 0) ;
  308.  
  309.  
  310. DTI = 1.D0/DeltaT ;
  311. 'REPETER' bloc1 nbespece ;
  312.  
  313. espc = ('EXTRAIRE' &bloc1 nomespec) ;
  314.  
  315. * On extrait la composante de Cini, Tcini et de la source
  316. CCini = 'NOMC' 'H' ('EXCO' espc Cini) ;
  317. 'SI' (nbsource > 1) ;
  318. SSource = 'NOMC' 'SOUR' ('EXCO' espc
  319. Chpsour) ;
  320. 'SINON' ;
  321. SSource = 'NOMC' 'SOUR' Chpsour ;
  322. 'FINSI' ;
  323. * Conditions initiales
  324. * Prise en compte du terme source et eventuellement
  325. * de la convection avec le schema centre
  326. HHS = 'NOMC' 'SCAL' CCini ;
  327. * Prise en compte des conditions aux limites
  328. 'SI' (FLUNEU) ;
  329. EXFLU =
  330. ('NOMC' 'FLUX' ('EXCO' espc
  331. CLFLUX)) ;
  332. * MESS 'EXFLU ='; LIST EXFLU;
  333. 'FINSI' ;
  334. 'SI' (DIRCLI) ;
  335. EXDIR =
  336. ('NOMC' 'SCAL' ('EXCO' espc
  337. CLDIRI)) ;
  338. * MESS 'EXDIR ='; LIST EXDIR;
  339. 'FINSI' ;
  340. 'SI' (FLUTOT) ;
  341. EXFLUT =
  342. ('NOMC' 'FLUX' ('EXCO' espc
  343. CLFLT )) ;
  344. QLIM = NOMC(EXFLUT) 'FLUX' ;
  345. SUP = 'EXTR' EXFLUT 'MAIL' ;
  346. XPAR1 = 1.D0 + (0.0*CLFLT) ;
  347. XPAR1 = 'NOMC' XPAR1 'PAR1' ;
  348. 'SI' LCONV ;
  349. USCNR = 'REDU' QFACE SUP ;
  350. MUSCN = USCNR*(-1.D0) ;
  351. 'SINON' ;
  352. MUSCN = (0.D0*CLFLT) ;
  353. 'FINSI' ;
  354.  
  355. XPAR2 = 'NOMC' MUSCN 'PAR2' ;
  356. EXFLUT = XPAR1 + XPAR2 + QLIM ;
  357. EXFLUT = CHAN 'ATTRIBUT' EXFLUT NATURE DISCRET ;
  358. * MESS 'EXFLUT ='; LIST EXFLUT;
  359. 'FINSI' ;
  360.  
  361. 'SI' (FLUMIX) ;
  362. EXFLUM =
  363. ('NOMC' 'FLUX' ('EXCO' espc
  364. CLFLUX3 )) ;
  365. QLIM = NOMC(EXFLUM) 'FLUX' ;
  366. SUP = 'EXTR' EXFLUM 'MAIL' ;
  367. XPAR1 = 1.D0 + (0.0*CLFLUX3) ;
  368. * MESS 'XPAR1 =';
  369. XPAR1 = 'NOMC' XPAR1 'PAR1' ;
  370. MUSCN = coefm;
  371. XPAR2 = 'NOMC' MUSCN 'PAR2' ;
  372. EXFLUM = XPAR1 + XPAR2 + QLIM ;
  373. EXFLUM = CHAN 'ATTRIBUT' EXFLUM NATURE DISCRET ;
  374. EXFLUT = EXFLUT 'ET' EXFLUM ;
  375. 'FINSI' ;
  376. * MESS 'UPDAVF EXFLUT ='; LIST EXFLUT;
  377. * MESS 'UPDAVF EXFLU ='; LIST EXFLU;
  378.  
  379. * 'MESS' 'MODDIFF'; LIST MODDIFF;
  380. * 'MESS' 'MODCONV'; LIST MODCONV;
  381. 'SI' (MODDIFF 'OU' MODCONV) ;
  382. * On recalcule deux fois les coefficients. En général, cela doit etre
  383. * négligeable en temps de calcul;
  384. GRADT0 MCHAMT1 = 'PENT' MoDARCY 'FACE' 'MPFA'
  385. HHS 'DISPDIF' difftot 'TIMP' EXDIR 'QIMP' EXFLU
  386. 'MIXT' EXFLUT ;
  387. MCHAMT = MCHAMT1;
  388. MCHAJA = TetaDiff*MCHAMT1;
  389.  
  390. 'SI' LCONV ;
  391. GRADT0 MCHAMT = 'PENT' MoDARCY 'FACE' 'MPFA'
  392. HHS 'DISPDIF' difftot 'TIMP' EXDIR 'QIMP' EXFLU
  393. 'MIXT' EXFLUT 'UPWICENT' QFace ;
  394. * ATTENTION PEUT ETRE FAUX SI TetaConv different de TetaDiff
  395. MCHCONV = MCHAMT + ((-1.D0)*MCHAMT1) ;
  396. MCHAJA = (TetaConv*MCHCONV) + (TetaDiff*MCHAMT1) ;
  397. 'DETRUIT' MCHCONV ;
  398. 'FINSI' ;
  399.  
  400.  
  401. * On ne calcule la jacobienne que pour la première espèce
  402. 'SI' ((&bloc1) 'EGA' 1) ;
  403. JACO CHPRES DT = 'LAPN' 'VF' 'CLAUDEIS' 'IMPL'
  404. MoDARCY HHS GRADT0 MCHAJA 'QIMP' EXFLU 'MIXT' EXFLUT
  405. 'TIMP' EXDIR ;
  406. 'SINON' ;
  407. JACOB CHPRES DT = 'LAPN' 'VF' 'CLAUDEIS' 'EXPL'
  408. MoDARCY HHS GRADT0 ;
  409. 'FINSI' ;
  410. * SINON;
  411. * 'DETRUIT' MCHAJA ;
  412. 'SI' ( ((&bloc1) 'EGA' 1) ) ;
  413. MATOT = 'KOPS' 'MULT' -1.000D0 JACO ;
  414. MATI = '*' MPOR2 DTI ;
  415. MATOT = MATOT 'ET' MATI ;
  416. 'FINSI' ;
  417. 'SINON';
  418. 'SI' ( ((&bloc1) 'EGA' 1) 'ET' (MODPOROS 'OU' MODDELTT)) ;
  419. MATOT = 'KOPS' 'MULT' -1.000D0 JACO ;
  420. MATI = '*' MPOR2 DTI ;
  421. MATOT = MATOT 'ET' MATI ;
  422. 'FINSI';
  423. * 'SINON';
  424. * On reconstitue un champ de second membre
  425. GRADT0 = 'PENT' MoDARCY 'FACE' 'MPFA'
  426. HHS 'DISPDIF' difftot
  427. 'TIMP' EXDIR 'QIMP' EXFLU 'MIXT' EXFLUT 'UPWICENT' QFace
  428. 'GRADGEO' MCHAMT ;
  429. JACO CHPRES DT = 'LAPN' 'VF' 'CLAUDEIS' 'EXPL'
  430. MoDARCY HHS GRADT0 ;
  431.  
  432. 'FINSI' ;
  433.  
  434. 'SI' ((&bloc1) 'EGA' 1) ;
  435. SSMTr = 'NOMC' espc CHPRES ;
  436. RESI = SSMTr ;
  437. 'SINON' ;
  438. CHPRES = 'NOMC' espc
  439. ('COPIER' CHPRES) ;
  440. RESI = CHPRES 'ET' RESI;
  441. 'FINSI' ;
  442. 'FIN' bloc1 ;
  443.  
  444. 'FINSI' ;
  445.  
  446.  
  447. *---------------------------------------------------------------------
  448. *------ Conditions aux limites de Flux (mixtes et Neumann) -----------
  449. *---------------------------------------------------------------------
  450.  
  451.  
  452. *---------------------------------------------------------------------
  453. *------ Matrice sortante --------------------------------- -----------
  454. *---------------------------------------------------------------------
  455.  
  456.  
  457. * si pas de modif et matrice matrik en entrée
  458. NOUVMAT = VRAI;
  459. 'SI' (('EXISTE' mattr) 'ET' ('NON' (MODDIFF 'OU' MODCONV 'OU'
  460. MODPOROS 'OU' MODDELTT))) ;
  461. Matot = mattr ;
  462. NOUVMAT = FAUX ;
  463. 'FINSI' ;
  464.  
  465.  
  466. 'FINP' RESI Matot Difftot MCHAMT MCHAMT1 NOUVMAT ;
  467.  
  468.  
  469.  
  470.  
  471.  
  472.  
  473.  
  474.  
  475.  
  476.  

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