Télécharger la plaquette Cast3M


fr en

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.  

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