Télécharger la plaquette Cast3M


fr en

Télécharger calctrac.procedur

Retour à la liste

Numérotation des lignes : 

  1. * CALCTRAC  PROCEDUR  GOUNAND   11/05/24    21:15:01     6976           
  2. 'DEBP'  CALCTRAC MoDARCY*'MMODEL' Difftot*'MCHAML' Cini*'CHPOINT'
  3.                  nomespec*'LISTMOTS' nbespece*'ENTIER' LMLump*'LOGIQUE'
  4.                  matrtr/'MATRIK' TABRES*'TABLE' Tbdartra*'TABLE'
  5.                  CHCLIM*'TABLE';                 
  6. * ATTENTION La matrice matrtr est optionnelle, L'ordre est important
  7. * et les types d'arguments qui se suivent aussi pour tester leur
  8. * présence
  9. *
  10. *  |-----------------------------------------------------------------|
  11. *  | Phrase d'appel (en GIBIANE)                                     |
  12. *  |-----------------------------------------------------------------|
  13. *  |                                                                 |
  14. *  |tcfin MatrTr = CALCTRAC MoDARCY Difftot Cini                     |
  15. *  |               nomespec nbespece  LMLump  (matrtr)               |
  16. *  |               TABRES Tbdartra CHCLIM;                           |
  17. *  |                                                                 |
  18. *  |-----------------------------------------------------------------|
  19. *  | Généralités : CALCTRAC calcule les traces de concentration      |
  20. *  |               associées à la donnée de concentrations initiales |
  21. *  |               Les Conditions limites de flux et de concentration|
  22. *  |               sont pris en compte.                              |
  23. *  |-----------------------------------------------------------------|
  24. *  |                                                                 |
  25. *  |-----------------------------------------------------------------|
  26. *  |                         ENTREES                                 |
  27. *  |-----------------------------------------------------------------|
  28. *  | MoDARCY  : modele Darcy.                                        |
  29. *  |                                                                 |
  30. *  | Difftot  : Coefficient de diffusion totale, integre decentrement|
  31. *  |                                                                 |
  32. *  | Cini     : Concentration initiale, CHPOINT centre.              |
  33. *  |            Composante 'H'.                                      |
  34. *  |                                                                 |
  35. *  | nomespec : liste des noms de composante des espèces dans Cini   |
  36. *  |                                                                 |
  37. *  | nbespece : nombre de composante de Cini, soit nombre d'especes  |
  38. *  |                                                                 |
  39. *  | LMLump   : Logique. Si vrai on effectue une condensation de     |
  40. *  |            masse de la matrice EFMH                             |
  41. *  |                                                                 |
  42. *  | TABRES   : table contenant les options de résolution pour KRES  |
  43. *  |                                                                 |
  44. *  | TbDarTra : Table Darcy transitoire utilisée par MHYB, SMTP ...  |
  45. *  |                                                                 |
  46. *  | CHCLIM   : table d'indice 'NEUMANN' et 'DIRICHLET' contenant les|
  47. *  |            Chpoint à n composantes contenant les conditions aux |
  48. *  |            limites de Neumann et Dirichlet par espece.          |
  49. *  |                                                                 |
  50. *  |-----------------------------------------------------------------|
  51. *  |                         ENTREES-SORTIES                         |
  52. *  |-----------------------------------------------------------------|
  53. *  |                                                                 |
  54. *  | MatrTr   : matrice globale sur les traces. MATRIK en entrée     |
  55. *  |            sort MATRIK si non modifiée, RIGIDITE sinon          |
  56. *  |                                                                 |
  57. *  |                                                                 |
  58. *  |-----------------------------------------------------------------|
  59. *  |                         SORTIES                                 |
  60. *  |-----------------------------------------------------------------|
  61. *  |                                                                 |
  62. *  |                                                                 |
  63. *  | Tcfin   : Trace de concentration aux faces (une composante par  |
  64. *  |           espece chimique)                                      |
  65. *  |                                                                 |
  66. *  |-----------------------------------------------------------------|
  67. *  |                  VARIABLES INTERNES                             |
  68. *  |-----------------------------------------------------------------|
  69. *  |                                                                 |
  70. *  | CCini    : concentration aux centres (une composante)           |
  71. *  |                                                                 |
  72. *  | LCALMATP : Logique, VRAI si on recalcule la matrice du systeme  |
  73. *  |            avec diffusion produite par MATP                     |
  74. *  |                                                                 |
  75. *  | TbDarTrb : table Darcy transitoire utilisée par MHYB, SMTP ...  |
  76. *  |                                                                 |
  77. *  | MassEFMH : matrice elementaire EFMH                             |
  78. *  |                                                                 |
  79. *  | Tccini   : Trace de concentration aux faces (une composante)    |
  80. *  |                                                                 |
  81. *  | SSMTr    : second membre sur les traces pour une espèce         |
  82. *  |                                                                 |
  83. *  | SMTr     : second membre sur les traces                         |
  84. *  |                                                                 |
  85. *  | TABSORT  : copie de TABRES                                      |
  86. *  |                                                                 |
  87. *  | dum1     : champoint vide                                       |
  88. *  |                                                                 |
  89. *  | matvide  : matrice rigidité vide                                |
  90. *  |                                                                 |
  91. *  | tccfin   : trace de concentration aux faces                     |
  92. *  |                                                                 |
  93. *  | DIRCLI   : logique valant VRAI si conditions aux                |
  94. *  |            limites de Dirichlet                                 | 
  95. *  |                                                                 |
  96. *  | CLDIRI   : Chpoint à n composantes contenant les conditions aux |
  97. *  |            limites de Dirichlet par espece.                     |
  98. *  |            il faudra en faire un nuage si supports géométriques |
  99. *  |            différents par espece. OPTIONNEL                     |
  100. *  |                                                                 |
  101. *  | FLUNEU   : LOGIQUE valant VRAI si conditions de Neumann         |
  102. *  |                                                                 |
  103. *  | CLFLUX   : Chpoint à n composantes contenant les flux imposés   |
  104. *  |            pour chaque espece chimique. nul si pas de flux      |
  105. *  |            OPTIONNEL                                            |
  106. *  |                                                                 |
  107. **********************************************************************
  108.  
  109.  
  110. *---------------------------------------------------------------------
  111. *---------- On récupere les conditions limites      ------------------
  112. *---------------------------------------------------------------------
  113.  
  114. * Flag sur conditions limites initialisés
  115.  FLUNEU = FAUX                                                         ;
  116.  FLUTOT = FAUX                                                         ;
  117.  FLUCLI = FAUX                                                         ;
  118.  DIRCLI = FAUX                                                         ;
  119.  FTOCLI = FAUX                                                         ;
  120.  
  121.  
  122. * Conditions flux total
  123. 'SI' ('EXISTE' CHCLIM 'FLUTOTAL')                                      ;
  124.    FTOCLI = VRAI                                                       ;
  125. 'FINSI'                                                                ;
  126.  
  127. * Neumann
  128. 'SI' ('EXISTE' CHCLIM 'NEUMANN')                                       ;
  129.    CLFLUX1 = CHCLIM . 'NEUMANN'                                        ;
  130.    FLUNEU = VRAI                                                       ;
  131. 'FINSI'                                                                ;
  132.  
  133. * Flux total
  134. 'SI' ('EXISTE' CHCLIM 'FLUTOTAL')                                     ;
  135.    CLFLUX2 = CHCLIM . 'FLUTOTAL'                                      ;
  136.    FLUTOT = VRAI                                                      ;
  137. 'FINSI'                                                               ;
  138.  
  139. * On fabrique le terme de flux complet
  140. 'SI' (FLUNEU 'ET' FLUTOT)                                             ;
  141.    CLFLUX = CLFLUX1 'ET' CLFLUX2                                      ;
  142.    FLUCLI = VRAI                                                      ;
  143. 'SINON'                                                               ;
  144.    'SI' (FLUNEU)                                                      ;
  145.       CLFLUX = CLFLUX1                                                ;
  146.       FLUCLI = VRAI                                                   ;
  147.    'FINSI'                                                            ;
  148.    'SI' (FLUTOT)                                                      ;
  149.       CLFLUX = CLFLUX2                                                ;
  150.       FLUCLI = VRAI                                                   ;       
  151.    'FINSI'                                                            ;
  152. 'FINSI'                                                               ; 
  153.  
  154.  
  155. * Dirichlet
  156. 'SI' ('EXISTE' CHCLIM 'DIRICHLET')                                     ;
  157.    CLDIRI = CHCLIM . 'DIRICHLET'                                       ;
  158.    DIRCLI = VRAI                                                       ;  
  159. 'FINSI'                                                                ;
  160.  
  161.  
  162. *---------------------------------------------------------------------
  163. *---------- Initialisations de tables, coefficients ------------------
  164. *---------------------------------------------------------------------
  165.  
  166.  
  167. * On regarde si une matrice est présente.
  168. 'SI' ('EXISTE' matrtr) ;
  169. *  On ne recalcule pas la matrice du probleme
  170.    LCALMATP = FAUX;
  171. 'SINON' ;
  172. *  On recalcule la matrice du probleme
  173.    LCALMATP = VRAI;
  174. 'FINSI' ;
  175.  
  176.  
  177.  
  178. *
  179. * On recopie la table de résolution TABRES dans TABSORT
  180. * Attention si une valeur contenue dans la table a le
  181. * debut d'un nom d'opérateur de castem, il y a probleme
  182. * d'ou démarrage apres l'indice soustype de valeur
  183. * METHINV, identique à opérateur METHode.
  184. *
  185.  
  186.  
  187. dumm = 'INDEX' TABRES;
  188. TABSORT = 'TABLE' METHINV;
  189. 'REPETER' bou1 (('DIME' dumm) '-' 1);
  190.    TABSORT . (mot dumm . (&bou1 '+' 1)) = TABRES . (mot dumm . (&bou1
  191.    '+' 1));
  192. 'FIN' bou1;
  193. * On impose un gradient conjugué pour la résolution (pour l'instant)
  194. TABSORT . 'TYPINV' = 2;
  195.  
  196.  
  197. *
  198. * on reconstruit une table darcy transitoire.
  199. *
  200.  
  201. TbDarTrb = 'TABLE' 'DARCY_TRANSITOIRE';
  202. * on prend un schéma Euler implicite car dt = 0.
  203. TbDarTrb . 'THETA' = 1.D0 ;
  204. * Pas non nul car test sur valeur négative pas terrible
  205. * dans SMTP
  206. TbDarTrb . 'PAS' = 1.D-90 ;
  207. * La valeur du terme suivant n'a pas d'influence car elle est 
  208. * multipliée par le pas de temps nul TbDarTrb . 'PAS' .
  209. TbDarTrb . 'SURF' = TbDarTra . 'SURF' ;
  210.  
  211.  
  212. *---------------------------------------------------------------------
  213. *-------------- Matrice masse inverse des EFMH -----------------------
  214. *---------------------------------------------------------------------
  215.  
  216.  
  217. * Calcul des  matrices de masse elementaires inverses
  218. * affectées par un changement de la diffusivité totale.
  219. * Cela ce porduit si la diffusivité change ou, en présence
  220. * de décentrement, si la diffusivité numérique change.
  221. 'SI' (LCALMATP);
  222.    'SI' (LMLump) ;
  223. *      masse lumping
  224.        MassEFMH = 'MHYB' MoDARCY  Difftot 'LUMP' ;
  225.    'SINON' ;
  226.        MassEFMH = 'MHYB' MoDARCY  Difftot ;
  227.    'FINSI' ;
  228. 'FINSI' ;
  229.  
  230.  
  231. *----------------------------------------------------------------------
  232. *---------- CALCUL DE LA MATRICE DE DIFFUSION DU PROBL ----------------
  233. *----------------------------------------------------------------------
  234.  
  235.  
  236. * Calcul de la matrice du probleme diffusion transitoire pour dt=0
  237.  
  238. 'SI' (LCALMATP) ;
  239.    MatrTr =  'MATP' MoDARCY MassEFMH  TbDarTrb ;
  240. 'FINSI' ;
  241.  
  242.  
  243.  
  244. *----------------------------------------------------------------------
  245. *--- CALCUL DE LA MATRICE DE CONVECTION ET DES SECONDS MEMBRES --------
  246. *----------------------------------------------------------------------
  247.  
  248.  
  249. * Nombres d'espèces à gérer est fixé et ne peut changer au cours d'une
  250. * simulation. Il sera d'ailleurs mis dans le modele
  251.  
  252. 'SI' (nbespece > 0) ;
  253.  
  254. *  pour un schéma en temps non euler implicite, il faut
  255. *  la trace à l'instant précédent pour le second membre
  256. *  pour la convection ou la diffusion. Elle est issue
  257. *  des calculs précédents en EFMH. C'est impératif.
  258. *  On ne prévoit pas de la recalculer car c'est très cher.
  259.  
  260. *
  261. *  CALCUL DES SECONDS MEMBRES
  262. *
  263.  
  264.  
  265.    'REPETER' bloc1 nbespece;
  266. *     On extrait la composante de Cini, Tcini et de la source
  267.       CCini = 'NOMC' 'H' ('EXCO' (extr &bloc1 nomespec) Cini);
  268. *     La trace de concentration initiale n'a pas de role et est mise à 0.
  269.       TCCINI = 'NOMC' 'TH' ('KCHT' Modarcy SCAL 'FACE' 0.D0);
  270. *     Conditions initiales
  271.       TbDarTrb . 'CHARGE' = CCini ;
  272.       TbDarTrb . 'TRACE'= TCCini ;
  273. *     prend en compte uniquement la concentration initiale
  274.       SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTrb ;   
  275. *     On reconstitue un champ de second membre
  276.       'SI' ((&bloc1) 'EGA' 1) ;
  277.          SSMTr = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) SSMTr;
  278.          SMTR =  SSMTr;
  279.       'SINON' ;
  280.          SSMTr = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) ('COPIER' SSMTr);
  281.          SMTr = SSMTr 'ET'  SMTr; 
  282.       'FINSI' ;
  283.    'FIN' bloc1;
  284.  
  285. 'FINSI' ;
  286.  
  287.  
  288. *---------------------------------------------------------------------
  289. *------ Conditions aux limites de Flux (mixtes et Neumann) -----------
  290. *---------------------------------------------------------------------
  291.  
  292. 'SI' (FLUCLI) ;  
  293.    SMTr = SMTR 'ET' CLFLUX ;
  294. 'FINSI' ;
  295.  
  296.  
  297. *---------------------------------------------------------------------
  298. *------ Conditions aux limites de Dirichlet                -----------
  299. *---------------------------------------------------------------------
  300.  
  301.  
  302.  
  303. *'SI' (EXISTE CLDIRI) ;
  304. * On fera peut etre une matrice bloque pour les méthodes directes.
  305. *'FINSI' ;
  306.  
  307.  
  308. *---------------------------------------------------------------------
  309. *------ Résolution                                         -----------
  310. *---------------------------------------------------------------------
  311.  
  312. 'SI' (LCALMATP)                                                       ;
  313. *  On recalcule les conditions aux limites flux total si existent
  314.    'SI' (FTOCLI)                                                      ;
  315.       'SI' ('EXISTE' QFACE)                                           ;
  316.          'SI' ('EGA' ('TYPE' CHCLIM . 'FLUTOTAL') 'CHPOINT')          ;
  317.             mayage = 'EXTRAIRE' CHCLIM . 'FLUTOTAL' maillage          ;
  318.             vites  = 'REDU' QFACE mayage                              ;
  319.             matcli = 'KOPS' 'MATDIAGO' ('NOMC' 'TH' vites) 'MATRIK'   ;
  320.          'FINSI'                                                      ;
  321.       'FINSI'                                                         ;
  322.    'FINSI'                                                            ;
  323.  
  324. *  Nouvelle matrice rigidité, on la transforme en Matrik
  325.    matrtr = 'KOPS' 'RIMA' matrtr                                      ;
  326.    matrtr = KOPS 'CHANINCO' matrtr
  327.            ('MOTS' 'TH' 'LX') ('MOTS' 'TH' 'LX')
  328.            ('MOTS' 'FLUX' 'FLX') ('MOTS' 'TH' 'LX')                   ;
  329. * On rajoute les conditions aux limites de flux totale (matrice U tC)
  330.    'SI' ('EGA' ('TYPE' matcli) 'MATRIK' )                             ;
  331.       matrtr = matrtr 'ET' matcli                                     ;
  332.    'FINSI'                                                            ;
  333. *  On initialise les indices pointants sur les matrices
  334. *  assembmées et préconditionnées
  335.    TABSORT . 'MATASS' = Matrtr ;
  336.    TABSORT . 'MAPREC' = Matrtr ;
  337. 'SINON' ;
  338.    TABSORT . 'MATASS' = Matrtr ;
  339.    TABSORT . 'MAPREC' = Matrtr ;
  340. 'FINSI' ;
  341.  
  342.  
  343.  
  344.  
  345. *---------------------------------------------------------------------
  346. *-------------- RESOLUTION EN TRACE DE CONCENTRATION -----------------
  347. *---------------------------------------------------------------------
  348.  
  349. * On fabrique un champoint vide et une matrice vide.
  350. dum1 matvide = 'KOPS' MATRIK;
  351.  
  352. * boucle sur les espèces.
  353. 'REPETER' bloc2 nbespece ;
  354. *  préparation solution initiale nulle car non connue
  355.    TABSORT . 'XINIT' = TCCINI;
  356. *  préparation second membre
  357.    SSmtr = 'NOMC' 'TH' ('EXCO' ('EXTRAIRE' &bloc2 nomespec) Smtr);
  358.  
  359. *  Solution en trace
  360. *  Si conditions de Dirichlet
  361.    'SI' (DIRCLI) ;
  362.       Tccfin = KRES matrtr 'TYPI' TABSORT
  363.             'SMBR' ('NOMC' 'TH' SSMTr)
  364.             'CLIM' ('NOMC' 'TH' ('EXCO' ('EXTRAIRE' &bloc2 nomespec)
  365.                      CLDIRI))
  366.             'IMPR' 0 ;
  367.    'SINON' ;
  368.       Tccfin = KRES matrtr 'TYPI' TABSORT
  369.             'SMBR' ('NOMC' 'TH' SSMTr)
  370.             'IMPR' 1 ;
  371.    'FINSI' ;   
  372.  
  373.  
  374. * On reconstitue les champoints à plusieurs composante
  375.  
  376.   'SI' (&bloc2 'EGA' 1);
  377.      Tccfin = 'NOMC' ('EXTRAIRE' &bloc2 nomespec) TCCfin;
  378.      Tcfin = tccfin;
  379.    'SINON' ;
  380.      Tccfin = 'NOMC' ('EXTRAIRE' &bloc2 nomespec) ('COPIER' TCCfin);
  381.      Tcfin = tcfin 'ET' tccfin;
  382.    'FINSI' ;
  383.  
  384. 'FIN' bloc2;   
  385.  
  386.  
  387.  
  388. 'FINP' tcfin MatrTr ;
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  

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