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