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 :
© Cast3M 2003 - All rights reserved.
Disclaimer