1 : * INITEFMH PROCEDUR GOUNAND 05/02/16 21:15:51 5029 2 : ********************************************************************** 3 : 'DEBP' INITEFMH MoDARCY*'MMODEL' Porosite*'CHPOINT' MateDiff*'CHPOINT' 4 : ChPSour*'CHPOINT' DeltaT*'FLOTTANT' Cini*'CHPOINT' 5 : TetaDiff*'FLOTTANT' TetaConv*'FLOTTANT' 6 : TetaLin*'FLOTTANT' Qface/'CHPOINT' QELEM/'CHPOINT' 7 : DISPL/'CHPOINT' DISPT/'CHPOINT' 8 : LMLump*'LOGIQUE' DECENTR*'LOGIQUE' CHCLIM*'TABLE' 9 : OPTRESOL/'TABLE'; 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 MatrTr CoefDt TbDarTra MassEFMH nomespec | 22 : * | nbespece Difftot Tcini TABRES TABMODI= INITEFMH MoDARCY Porosite| 23 : * | MateDiff ChPSour DeltaT Cini TetaDiff | 24 : * | TetaConv TetaLin fluimp dircli (QFACE) | 25 : * | LMLump DECENTR CHCLIM optresol ; | 26 : * | | 27 : * |-----------------------------------------------------------------| 28 : * | Généralités : MATTEFMH construit la matrice de discrétisation | 29 : * | du problème de transport convection-diffusion pour| 30 : * | le premier pas de tps d'un algorithme transitoire.| 31 : * | Le second membre et les Conditions limites de flux| 32 : * | sont pris en compte. | 33 : * | RESTE TCINI, DECENTR et TERME LIN | 34 : * |-----------------------------------------------------------------| 35 : * | | 36 : * |-----------------------------------------------------------------| 37 : * | ENTREES | 38 : * |-----------------------------------------------------------------| 39 : * | MoDARCY : modele Darcy. | 40 : * | | 41 : * | Porosite : champoint de composante 'CK' | 42 : * | | 43 : * | MateDiff : Tenseur de diffusion (type iso, ..) champoint | 44 : * | de composante 'K' en isotrope, 'K11', 'K21', | 45 : * | 'K22' en anisotrope 2d et 'K11', 'K21', 'K22', 'K31'| 46 : * | 'K32', 'K33' en anisotrope 3d. Type 'CARACTERISTIQUE'| 47 : * | | 48 : * | DeltaT : Pas de temps. | 49 : * | | 50 : * | ChPSour : Champ par points des sources volumiques par unité de | 51 : * | temps (support maillage centre). Composante ?????? | 52 : * | | 53 : * | Cini : Concentration initiale, CHPOINT centre. | 54 : * | Composante 'H'. | 55 : * | | 56 : * | Qface : vitesse aux faces, CHPO face de composantes Vx, Vy | 57 : * | en 2d et Vx, Vy, Vz en 3d. Il s'agit plus exatement | 58 : * | de (V.n)n, c'est à dire de la composante normale de | 59 : * | la vitesse aux faces. ???????? (je pressens que | 60 : * | castem va sortir des flux, cad intégrés sur surfaces)| 61 : * | | 62 : * | TetaDiff : Valeur de theta pour theta-schéma en temps, opérateur| 63 : * | de diffusion. Entre 0 et 1. 0 = explicite, 1 = euler | 64 : * | implicite. | 65 : * | | 66 : * | TetaConv : Valeur de theta pour theta-schéma en temps, opérateur| 67 : * | de convection. Entre 0 et 1. 0 = explicite, 1 = Euler| 68 : * | implicite. | 69 : * | | 70 : * | TetaLin : valeur de theta pour theta-schéma en temps, opérateur| 71 : * | linéaire du type coef * C, où C est l'inconnue. | 72 : * | Entre 0 et 1. 0 = explicite, 1 = euler implicite. | 73 : * | ??????????? A voir car peut etre identique à Tetadiff| 74 : * | | 75 : * | LMLump : Logique. Si vrai on effectue une condensation de | 76 : * | masse de la matrice EFMH | 77 : * | | 78 : * | DECENTR : Logique. Vrai veut dire schémas décentrés et faux | 79 : * | veut dire schéma convectif centré. | 80 : * | | 81 : * | CHCLIM : table d'indice 'NEUMANN' et 'DIRICHLET' contenant les| 82 : * | Chpoint à n composantes contenant les conditions aux | 83 : * | limites de Neumann et Dirichlet par espece. | 84 : * | L'indice 'FLUXTOT' contient les conditions limites | 85 : * | de flux total et 'FLUMIXTE' concerne une condition | 86 : * | de flux mixte : 'FLUMIXTE' . 'VAL' contient le champ | 87 : * | à n composantes indiquant le flux, 'FLUMIXTE' . 'A' | 88 : * | et 'FLUMIXTE' . 'B' les coef (champoints SCAL) tels | 89 : * | que A D grad (C) + B (C) = VAL | 90 : * | | 91 : * | OPTRESOL : Table dont l'entree est optionnelle définissant | 92 : * | les options de résolution pour 'KRES'. | 93 : * | | 94 : * |-----------------------------------------------------------------| 95 : * | SORTIES | 96 : * |-----------------------------------------------------------------| 97 : * | | 98 : * | | 99 : * | MassEFMH : matrice elementaire EFMH | 100 : * | | 101 : * | MatrTr : matrice globale sur les traces | 102 : * | | 103 : * | SMTr : second membre sur les traces | 104 : * | | 105 : * | TbDarTra : table Darcy transitoire utilisée par MHYB, SMTP ... | 106 : * | | 107 : * | nomespec : liste des noms de composante des espèces dans Cini | 108 : * | | 109 : * | nbespece : nombre de composante de Cini, soit nombre d'especes | 110 : * | | 111 : * | nbsource : nombre de composantes du terme source qd X especes | 112 : * | | 113 : * | Diffdisp : Dipersivité, tenseur chpoint K11 K22 K33 K21 K31 K32 | 114 : * | | 115 : * | TABRES : Table complète définissant les options de résolution | 116 : * | pour 'KRES'. | 117 : * | | 118 : * | Tcini : Trace de concentration aux faces (eventuellement à | 119 : * | plusieurs composantes (espèces) | 120 : * | | 121 : * | TABMODI : table contenant des logiques indiquant la nécessité | 122 : * | ou non de reclalculer certains termes. | 123 : * | 'POROSITE' : VRAI si le coefficient devant D/DT | 124 : * | (porosité) est modifié depuis le dernier| 125 : * | appel | 126 : * | 'DELTAT' : VRAI si le pas de tps a changé | 127 : * | 'CONVECTI' : VRAI si la vitesse a changé | 128 : * | 'COEF_LIN' : VRAI si le coef en facteur de C a changé| 129 : * | 'DIFFUSI' : VRAI si les diffusivités ont changé | 130 : * | | 131 : * |-----------------------------------------------------------------| 132 : * | VARIABLES INTERNES | 133 : * |-----------------------------------------------------------------| 134 : * | | 135 : * | CoefDt : coeff devant dC/dt integre sur les elements | 136 : * | | 137 : * | PCONV : Logique indiquant VRAI si présence de convection | 138 : * | | 139 : * | toltheta : 1.D-4 seuil en dessous duquel on considère que la | 140 : * | valeur de theta du theta-schéma est nulle (schéma | 141 : * | explicite) ou au contraire euler-implicite si | 142 : * | theta > 0.9999 | 143 : * | | 144 : * | Tccini : Trace de concentration aux faces (une composante) | 145 : * | | 146 : * | SSource : Source aux centre (une composante) | 147 : * | | 148 : * | CCini : concentration aux centres (une composante) | 149 : * | | 150 : * | lstcps : liste des noms de composante des espèces dans Chpsour| 151 : * | | 152 : * | SSMTr : second membre sur les traces pour une espèce | 153 : * | | 154 : * | MatConv : matrice globale sur les traces pour la convection | 155 : * | | 156 : * | Numdiff : diffusivité numérique due au décentrement | 157 : * | | 158 : * | FLUNEU : LOGIQUE valant VRAI si conditions de Neumann | 159 : * | | 160 : * | CLFLUX : Chpoint à n composantes contenant les flux imposés | 161 : * | pour chaque espece chimique. nul si pas de flux | 162 : * | OPTIONNEL | 163 : * | | 164 : ********************************************************************** 165 : 166 : *--------------------------------------------------------------------- 167 : *---------- On récupere les conditions limites ------------------ 168 : *--------------------------------------------------------------------- 169 : 170 : 171 : FLUNEU = FAUX ; 172 : FLUTOT = FAUX ; 173 : FLUMIX = FAUX ; 174 : FLUCLI = FAUX ; 175 : 176 : * Neumann 177 : 'SI' ('EXISTE' CHCLIM 'NEUMANN') ; 178 : CLFLUX1 = CHCLIM . 'NEUMANN' ; 179 : FLUNEU = VRAI ; 180 : 'SINON' ; 181 : * on crée un champ vide 182 : CLFLUX1 dum = 'KOPS' MATRIK ; 183 : 'OUBLIER' dum ; 184 : 'FINSI' ; 185 : 186 : * Flux total 187 : 'SI' ('EXISTE' CHCLIM 'FLUTOTAL') ; 188 : CLFLUX2 = CHCLIM . 'FLUTOTAL' ; 189 : FLUTOT = VRAI ; 190 : 'SINON' ; 191 : * on crée un champ vide 192 : CLFLUX2 dum = 'KOPS' MATRIK ; 193 : 'OUBLIER' dum ; 194 : 'FINSI' ; 195 : 196 : * Flux mixte 197 : 'SI' ('EXISTE' CHCLIM 'FLUMIXTE') ; 198 : * comme on impose A Dgrad C + B C = flumix, on le traite sous 199 : * la forme D grad C + (B/A) C = flumix/A plus naturelle en EFMH car 200 : * D grad C est le flux diffusif 201 : COFA = -1.D0 * CHCLIM . 'FLUMIXTE' . 'COEFA' ; 202 : CLFLUX3 = CHCLIM . 'FLUMIXTE' . 'VAL' '/' COFA ; 203 : CLFLUX3 = CHAN 'ATTRIBUT' CLFLUX3 NATURE DISCRET ; 204 : FLUMIX = VRAI ; 205 : 'SINON' ; 206 : * on crée un champ vide 207 : CLFLUX3 dum = 'KOPS' MATRIK ; 208 : 'OUBLIER' dum ; 209 : 'FINSI' ; 210 : 211 : 212 : * On fabrique le terme de flux complet 213 : 'SI' (FLUNEU 'OU' FLUTOT 'OU' FLUMIX) ; 214 : CLFLUX = CLFLUX1 'ET' CLFLUX2 'ET' CLFLUX3 ; 215 : FLUCLI = VRAI ; 216 : 'FINSI' ; 217 : 218 : 219 : *--------------------------------------------------------------------- 220 : *---------- Initialisations de tables, coefficients ------------------ 221 : *--------------------------------------------------------------------- 222 : 223 : * Table de logiques indiquant des modifications. Initialisation 224 : TABMODI = TABLE; 225 : TABMODI . 'POROSITE' = FAUX; 226 : TABMODI . 'CONVECTI' = FAUX; 227 : TABMODI . 'DELTAT' = FAUX; 228 : TABMODI . 'COEF_LIN' = FAUX; 229 : TABMODI . 'DIFFUSIV' = FAUX; 230 : 231 : 232 : 'SI' ('EXISTE' QFACE) ; 233 : PCONV = VRAI; 234 : 'SI' ('EXISTE' DISPL) ; 235 : DISPERSI = VRAI ; 236 : 'SINON' ; 237 : DISPERSI = FAUX; 238 : 'FINSI' ; 239 : 'SINON' ; 240 : PCONV = FAUX; 241 : DISPERSI = FAUX; 242 : 'FINSI' ; 243 : 244 : 245 : * tolerance sur theta du theta schéma de discrétisation en temps. 246 : * il faudrait remmettre les theta à 0 ou 1 si nécessaire dans 247 : * procédure amont. 248 : toltheta = 1.D-4; 249 : 250 : 'SI' (TetaConv 'EGA' 0.D0 toltheta) ; 251 : TetaConv = 0.D0; 252 : 'FINSI' ; 253 : 'SI' (TetaConv 'EGA' 1.D0 toltheta) ; 254 : TetaConv = 1.D0; 255 : 'FINSI' ; 256 : 'SI' (TetaDiff 'EGA' 0.D0 toltheta) ; 257 : TetaDiff = 0.D0; 258 : 'FINSI' ; 259 : 'SI' (TetaDiff 'EGA' 1.D0 toltheta) ; 260 : TetaDiff = 1.D0; 261 : 'FINSI' ; 262 : 'SI' (TetaLin 'EGA' 0.D0 toltheta) ; 263 : TetaLin = 0.D0; 264 : 'FINSI' ; 265 : 'SI' (TetaLin 'EGA' 1.D0 toltheta) ; 266 : TetaLin = 1.D0; 267 : 'FINSI' ; 268 : 269 : 270 : * Calcul du terme devant le dC/dt integré sur le volume 271 : * CHAMELEM 'SCAL'. Voir pour VF ????? 272 : CoefDt = ('NOMC' 'SCAL' Porosite); 273 : COefDt = Porosite * ('DOMA' modarcy VOLUME); 274 : CoefDt = 'KCHA' Modarcy CoefDt 'CHAM'; 275 : 276 : * Creation de la table TbDarTra utilisée par les operateurs MATP, 277 : * SMTP, HYBP et HDEBI, dans le cadre des EFMH. invariante si 278 : * le pas de temps ne bouge pas. 279 : TbDarTra = 'TABLE' 'DARCY_TRANSITOIRE'; 280 : TbDarTra . 'THETA' = TetaDiff ; 281 : * essayer si marche avec, sans indice !!!! 282 : 'SI' (PCONV) ; 283 : TbDarTra . 'THETA_CONVECTION'= TetaConv ; 284 : 'FINSI' ; 285 : TbDarTra . 'PAS' = DeltaT ; 286 : TbDarTra . 'SURF' = CoefDt ; 287 : 288 : *--------------------------------------------------------------------- 289 : *----------------------- CREATION TABLE POUR RESOLUTION -------------- 290 : *--------------------------------------------------------------------- 291 : 292 : 293 : **************** OPTIONS PAR DEFAUT ************************** 294 : 295 : * création de la table de résolution pour l'opérateur KRES 296 : * On crée la table de résolution avec les options par défaut 297 : * On y remplacera les valeurs définit par l'utilisateur ensuite. 298 : TABRES = 'TABLE' 'METHINV' ; 299 : * type d'inversion 300 : 'SI' (PCONV) ; 301 : * option BCGSTAB par défaut pour une matrice non symétrique 302 : METHRES = 3; 303 : 'SINON' ; 304 : * option gradient conjugué par défaut pour une matrice symétrique 305 : METHRES = 2; 306 : 'FINSI' ; 307 : TABRES . 'TYPINV' = METHRES ; 308 : 309 : * niveau d'impression. 310 : TABRES . 'IMPINV' = 0 ; 311 : 312 : * Type de renumérotation. Option SLOANE par défaut 313 : *TABRES . 'TYRENU' = 'SLOANE' ; 314 : TABRES . 'TYRENU' = 'SLOA' ; 315 : 316 : * La gestion des multiplicateurs sera modifiée 317 : * par la suite. Pas d'option pour l'instant 318 : TABRES . 'PCMLAG' = 'APR2' ; 319 : TABRES . 'OUBMAT' = 0 ; 320 : TABRES . 'SCALING' = 0 ; 321 : 322 : *INDICES SPÉCIFIQUES POUR UNE MÉTHODE ITÉRATIVE 323 : * Nombre maxi d'itérations 324 : TABRES . 'NITMAX' = 1500 ; 325 : * résidu pour la convergence de la méthode 326 : TABRES . 'RESID' = 1.D-15 ; 327 : * valeur minimale du pivot de la méthode 328 : TABRES . 'BCGSBTOL' = 1.D-120 ; 329 : * preconditionnement ILU(0) 330 : TABRES . 'PRECOND' = 3 ; 331 : *relaxation pour MILU0 332 : TABRES . 'MILURELX' = 1.D0 ; 333 : *GMRESTART 334 : TABRES . 'GMRESTRT' = 100 ; 335 : *ILUTLFIL 336 : TABRES . 'ILUTLFIL' = 2; 337 : * drop tolerence pour ILUT2 338 : TABRES . 'ILUTDTOL' = 0.D0; 339 : TABRES . 'ILUTPTOL' = 0.01D0; 340 : TABRES . 'ILUTALPH' = 0.D0; 341 : 342 : ************** OPTIONS UTILISATEUR ************************** 343 : 344 : 345 : 'SI' ('EGA' ('TYPE' OPTRESOL) 'TABLE') ; 346 : * L'utilisateur a défini des options pour la méthode 347 : * de résolution. 348 : 349 : * Type d'inversion 350 : 'SI' ('EXISTE' OPTRESOL 'TYPINV') ; 351 : TABRES . 'TYPINV' = OPTRESOL . 'TYPINV'; 352 : 'FINSI' ; 353 : 354 : * Niveau d'impression 355 : 'SI' ('EXISTE' OPTRESOL 'IMPINV') ; 356 : TABRES . 'IMPINV' = OPTRESOL . 'IMPINV'; 357 : 'FINSI' ; 358 : 359 : * Type de renumérotation. 360 : 'SI' ('EXISTE' OPTRESOL 'TYRENU') ; 361 : TABRES . 'TYRENU' = OPTRESOL . 'TYRENU'; 362 : 'FINSI' ; 363 : 364 : * Indices spécifiques aux méthodes itératives 365 : 'SI' ((TABRES . 'TYPINV') > 1); 366 : * Nombre maxi d'iterations 367 : 'SI' ('EXISTE' OPTRESOL 'NITMAX') ; 368 : TABRES . 'NITMAX' = OPTRESOL . 'NITMAX'; 369 : 'FINSI' ; 370 : * Valeur du résidu de la méthode 371 : 'SI' ('EXISTE' OPTRESOL 'RESID') ; 372 : TABRES . 'RESID' = OPTRESOL . 'RESID'; 373 : 'FINSI' ; 374 : * valeur minimal du pivot de la méthode 375 : 'SI' ('EXISTE' OPTRESOL 'BCGSBTOL') ; 376 : TABRES . 'BCGSBTOL' = OPTRESOL . 'BCGSBTOL'; 377 : 'FINSI' ; 378 : * precond par diagonale 379 : 'SI' ('EXISTE' OPTRESOL 'PRECOND') ; 380 : TABRES . 'PRECOND' = OPTRESOL . 'PRECOND'; 381 : 'FINSI' ; 382 : * precon ILUT2 383 : 'SI' ('EXISTE' OPTRESOL 'ILUTLFIL') ; 384 : TABRES . 'ILUTLFIL' = OPTRESOL . 'ILUTLFIL' ; 385 : 'FINSI' ; 386 : 'FINSI' ; 387 : * Pour GMRES 388 : 'SI' ((TABRES . 'TYPINV') EGA 5); 389 : 'SI' ('EXISTE' OPTRESOL 'GMRESTRT') ; 390 : TABRES . 'GMRESTRT' = OPTRESOL . 'GMRESTRT'; 391 : 'SINON' ; 392 : TABRES . 'GMRESTRT' = 50; 393 : 'FINSI' ; 394 : 'FINSI' ; 395 : 'FINSI' ; 396 : 397 : 398 : SI (('EGA' TABRES . 'PRECOND' 8) 'OU' ('EGA' TABRES . 'PRECOND' 7)); 399 : TABRES . 'ILUTDTOL' = 0.1D-2; 400 : FINSI ; 401 : 402 : *--------------------------------------------------------------------- 403 : *--------------------- CALCUL DE LA DISPERSIVITE---------------------- 404 : *--------------------------------------------------------------------- 405 : 406 : * On test si le modèle est anisotrope ou non (plus d'une composante si 407 : * anisotrope) 408 : compmat = 'EXTRAIRE' Modarcy 'MATERIAU' ; 409 : ANISO = ('DIME' compmat) > 1 ; 410 : *'SI' ('NON' ANISO) ; 411 : * 'SI' (DISPERSI 'OU' DECENTR) ; 412 : * 'MESSAGE' 'ERREUR - le modèle est déclaré isotrope en présence 413 : * de décentrement ou de dispersivité' ; 414 : * ERREUR 5 ; 415 : * 'SINON' ; 416 : * difftot = 'NOMC' 'K' Matediff ; 417 : * 'FINSI' ; 418 : *'FINSI' ; 419 : 420 : * Remise de la diffusion aux bonnes composantes aniso car + général 421 : 422 : * la ligne suivante ne marche pas, je dois introduire zozo ???? 423 : *Matediff = DIFFANIS Matediff 'EFMH' ; 424 : zozo = DIFFANIS Matediff 'EFMH' ANISO ; 425 : Matediff = zozo; 426 : 427 : * Seulement si présence de convection et présence de 428 : * dispersivité et disp_l et disp_t 429 : 430 : 'SI' (DISPERSI 'ET' PCONV) ; 431 : DISPL = 'NOMC' 'SCAL' DISPL ; 432 : DISPT = 'NOMC' 'SCAL' DISPT ; 433 : * on calcul la dispersivité 434 : diffdisp = CALCDISP QELEM DISPL DISPT ; 435 : 'SINON' ; 436 : diffdisp = 0.D0 * Matediff ; 437 : 'FINSI' ; 438 : 439 : Difftot = Matediff + diffdisp ; 440 : 441 : *--------------------------------------------------------------------- 442 : *--------------------- CALCUL DU DECENTREMENT ------------------------ 443 : *--------------------------------------------------------------------- 444 : 445 : * Seulement si présence de convection et option décentrée 446 : * Deuxième matrice masse ???? 447 : 448 : 449 : 'SI' (DECENTR 'ET' PCONV) ; 450 : 451 : 'MESSAGE' 'on utilise un schéma décentré pour la convection' ; 452 : 453 : * projection sur axe des x de la normale 454 : normax = 'DOMA' modarcy 'NORMALE' ; 455 : normax = ('EXCO' 'UX' normax) 'ABS' ; 456 : normax = 'NOMC' 'FLUX' normax ; 457 : normax = ('DOMA' modarcy SURFACE) * normax ; 458 : dxmail = 'DIVU' modarcy normax ('ABS' ('DOMA' modarcy ORIENTATION)) ; 459 : * x 2 car surface volume / (surf haut + surf bas) 460 : dxmail = 2.D0 * ('DOMA' modarcy 'VOLUME') '/' (dxmail) ; 461 : 462 : * projection sur axe des x de la normale 463 : normay = 'DOMA' modarcy 'NORMALE' ; 464 : normay = ('EXCO' 'UY' normay) 'ABS' ; 465 : normay = 'NOMC' 'FLUX' normay ; 466 : normay = ('DOMA' modarcy SURFACE) * normay ; 467 : dymail = 'DIVU' modarcy normay ('ABS' ('DOMA' modarcy ORIENTATION)) ; 468 : dymail = 2.D0 * ('DOMA' modarcy 'VOLUME') '/' (dymail) ; 469 : 470 : * porjection sur axe des x de la normale 471 : 'SI' ('EGA' ('VALEUR' 'DIME') 3) ; 472 : normaz = 'DOMA' modarcy 'NORMALE' ; 473 : normaz = ('EXCO' 'UZ' normaz) 'ABS' ; 474 : normaz = 'NOMC' 'FLUX' normaz ; 475 : normaz = ('DOMA' modarcy SURFACE) * normaz ; 476 : dzmail = 'DIVU' modarcy normaz ('ABS' ('DOMA' modarcy 477 : ORIENTATION)) ; 478 : dzmail = 2.D0 * ('DOMA' modarcy 'VOLUME') '/' (dzmail) ; 479 : 'FINSI' ; 480 : 481 : * Vitesses extraites par composantes et comme chpo centre 482 : V1 = 'EXCO' Qelem 'VX' 'SCAL' ; 483 : V2 = EXCO Qelem 'VY' 'SCAL' ; 484 : 'SI' ('EGA' ('VALEUR' 'DIME') 3) ; 485 : V3 = EXCO Qelem 'VZ' 'SCAL' ; 486 : 'FINSI' ; 487 : 488 : Pe = OPTRESOL . 'PECLET' ; 489 : *---on incorpore la diffusion numerique suivant X si PECLET > 2 ; 490 : dum = 'ABS' (dxmail * V1 '/' Pe) ; 491 : D11P = exco 'K11' difftot 'SCAL' ; 492 : m1 = masq D11P SUPERIEUR dum ; 493 : D11P = m1*D11P + ((1 - m1) * dum) ; 494 : D11P = 'NOMC' 'K11' D11P ; 495 : 496 : *---on incorpore la diffusion numerique suivant Y si PECLET > 2 ; 497 : dum = 'ABS' (dymail * V2 '/' Pe) ; 498 : D22P = exco 'K22' difftot 'SCAL' ; 499 : m1 = masq D22P SUPERIEUR dum ; 500 : D22P = m1*D22P + ((1 - m1) * dum) ; 501 : D22P = 'NOMC' 'K22' D22P ; 502 : 503 : *---on incorpore la diffusion numerique suivant Z si PECLET > 2 ; 504 : 'SI' ('EGA' ('VALEUR' 'DIME') 3) ; 505 : dum = 'ABS' (dzmail * V3 '/' Pe) ; 506 : D33P = exco 'K33' difftot 'SCAL' ; 507 : m1 = masq D33P SUPERIEUR dum ; 508 : D33P = m1*D33P + ((1 - m1) * dum) ; 509 : D33P = 'NOMC' 'K33' D33P ; 510 : 511 : diffdisp = D11P '+' D22P '+' D33P 512 : '+' ('EXCO' 'K21' difftot 'K21') 513 : '+' ('EXCO' 'K31' difftot 'K31') 514 : '+' ('EXCO' 'K32' difftot 'K32') 515 : - Matediff ; 516 : 'SINON' ; 517 : diffdisp = D11P '+' D22P '+' ('EXCO' 'K21' difftot 'K21') 518 : - Matediff ; 519 : 'FINSI' ; 520 : 521 : 'FINSI' ; 522 : 523 : * diffusion totale (dispers + effective + numérique) 524 : 525 : difftot = Matediff '+' diffdisp ; 526 : 527 : * changement en chamelem pour EFMH 528 : 529 : difftot = 'KCHA' Modarcy difftot 'CHAM' ; 530 : difftot = CHANGER difftot TYPE CARACTERISTIQUES ; 531 : 532 : SI (ANISO) ; 533 : 'SI' ('EGA' ('VALEUR' 'DIME') 2) ; 534 : difftot = 'MATE' Modarcy 'DIRECTION' 535 : (1. 0. ) 'PARALLELE' 'K11' ('EXCO' 'K11' difftot) 536 : 'K21' ('EXCO' 'K21' difftot) 537 : 'K22' ('EXCO' 'K22' difftot) ; 538 : 'FINSI' ; 539 : 'SI' ('EGA' ('VALEUR' 'DIME') 3) ; 540 : difftot = 'MATE' Modarcy 'DIRECTION' 541 : (1. 0. 0.) (0. 1. 0.) 'PARALLELE' 542 : 'K11' ('EXCO' 'K11' difftot) 543 : 'K21' ('EXCO' 'K21' difftot) 544 : 'K22' ('EXCO' 'K22' difftot) 545 : 'K31' ('EXCO' 'K31' difftot) 546 : 'K32' ('EXCO' 'K32' difftot) 547 : 'K33' ('EXCO' 'K33' difftot) ; 548 : 'FINSI' ; 549 : SINON ; 550 : difftot = 'MATE' Modarcy 'K' difftot ; 551 : FINSI ; 552 : 553 : 554 : *--------------------------------------------------------------------- 555 : *-------------- Matrice masse inverse des EFMH ----------------------- 556 : *--------------------------------------------------------------------- 557 : 558 : 559 : * Calcul des matrices de masse elementaires inverses 560 : 'SI' (LMLump) ; 561 : * masse lumping 562 : MassEFMH = 'MHYB' MoDARCY Difftot 'LUMP' ; 563 : 'SINON' ; 564 : MassEFMH = 'MHYB' MoDARCY Difftot ; 565 : 'FINSI' ; 566 : 567 : 568 : 569 : *---------------------------------------------------------------------- 570 : *---------- CALCUL DE LA MATRICE DE DIFFUSION DU PROBL ---------------- 571 : *---------------------------------------------------------------------- 572 : 573 : 574 : * Calcul de la matrice du probleme diffusion transitoire 575 : MatrTr = 'MATP' MoDARCY MassEFMH TbDarTra ; 576 : 577 : 578 : 579 : *---------------------------------------------------------------------- 580 : *--- CALCUL DE LA MATRICE DE CONVECTION ET DES SECONDS MEMBRES -------- 581 : *---------------------------------------------------------------------- 582 : 583 : 584 : ************** Nombres d'espèces à gérer ****************************** 585 : 586 : * évaluation du nombre d'espèces. 587 : nomespec = 'EXTRAIRE' Cini 'COMP'; 588 : nbespece = 'DIME' nomespec; 589 : * nb de termes sources, commun (1), ou egal au nombre 590 : * d'especes (nbcompi) 591 : lstcps = 'EXTRAIRE' Chpsour 'COMP'; 592 : nbsource = 'DIME' lstcps; 593 : 'SI' ((nbsource 'NEG' 1) 'ET' (nbsource 'NEG' nbespece)) ; 594 : 'MESSAGE' 'La source doit avoir le meme nombre de composantes'; 595 : 'MESSAGE' 'que les espèces ou 1 seule composante'; 596 : ERREUR 5; 597 : 'FINSI' ; 598 : 599 : 600 : * La matrice de convection ne dépend pas de Tcini mais 601 : * est calculée en meme temps que le calcul du second 602 : * membre. On effectue donc un traitement particulier 603 : * dans un cas multiespece pour gagner en temps de calcul. 604 : 605 : 606 : ********************* Une ou plusieurs espèces ************************ 607 : 608 : 'SI' (nbespece > 0) ; 609 : 610 : * pour un schéma en temps non euler implicite, il faut 611 : * la trace à l'instant précédent pour le second membre 612 : * pour la convection ou la diffusion 613 : 'SI' ((TetaDiff 'NEG' 1.D0 toltheta) 614 : 'OU' (TetaConv 'NEG' 1.D0 toltheta)) ; 615 : * Calcul de Tcini 616 : tcini dummy = CALCTRAC MoDARCY Difftot Cini 617 : nomespec nbespece LMLump 618 : TABRES Tbdartra CHCLIM; 619 : 'OUBLIER' dummy; 620 : 'SINON' ; 621 : * La trace de charge n'est pas réellement utilisée car multipliée 622 : * par 0. On la met à 0 pour simplifier les calculs. 623 : TCCINI = 'NOMC' 'TH' ('KCHT' Modarcy SCAL 'FACE' 0.D0); 624 : 'REPETER' bloctc nbespece; 625 : 'SI' (&bloctc 'EGA' 1); 626 : tcini = ('NOMC' ('EXTRAIRE' &bloctc nomespec) 627 : TCCINI); 628 : 'SINON'; 629 : tcini = tcini 'ET' ('NOMC' ('EXTRAIRE' 630 : &bloctc nomespec) ('COPIER' TCCINI)); 631 : 'FINSI' ; 632 : 'FIN' bloctc; 633 : 'FINSI'; 634 : 635 : 'REPETER' bloc1 nbespece; 636 : * On extrait la composante de Cini, Tcini et de la source 637 : CCini = 'NOMC' 'H' ('EXCO' (extr &bloc1 nomespec) Cini); 638 : TCCini = 'NOMC' 'TH' ('EXCO' (extr &bloc1 nomespec) Tcini); 639 : 'SI' (nbsource > 1) ; 640 : SSource = 'NOMC' 'SOUR' ('EXCO' (extr &bloc1 nomespec) 641 : Chpsour); 642 : 'SINON' ; 643 : SSource = 'NOMC' 'SOUR' Chpsour; 644 : 'FINSI' ; 645 : * Conditions initiales 646 : TbDarTra . 'CHARGE' = CCini ; 647 : TbDarTra . 'TRACE'= TCCini ; 648 : * Prise en compte du terme source et eventuellement 649 : * de la convection avec le schema centre 650 : 'SI' (PCONV); 651 : * convection 652 : 'SI' (TetaConv 'NEG' 0.0D0 toltheta); 653 : * schéma partiellement implicite, matrice MatConv 654 : MatConv SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTra SSource 655 : (NOMC 'FLUX' QFace); 656 : 'SINON' ; 657 : * schéma explicite, calcul du second membre uniquement 658 : SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTra SSource 659 : (NOMC 'FLUX' QFace) ; 660 : 'FINSI' ; 661 : 'SINON'; 662 : * pas de convection, calcul du second membre restant 663 : SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTra SSource; 664 : 'FINSI'; 665 : * On reconstitue un champ de second membre 666 : 'SI' ((&bloc1) 'EGA' 1) ; 667 : SSMTr = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) SSMTr; 668 : SMTR = SSMTr; 669 : 'SINON' ; 670 : SSMTr = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) ('COPIER' SSMTr); 671 : SMTr = SSMTr 'ET' SMTr; 672 : 'FINSI' ; 673 : 'FIN' bloc1; 674 : * on stoque la matrice en assemblant la matrice de convection 675 : * du dernier calcul (elles sont toutes identiques). On mettra 676 : * un jour une option pour les sortir si besoin seulement. 677 : 'SI' (PCONV 'ET' (TetaConv 'NEG' 0.0D0 toltheta)); 678 : MatrTr = MatrTr 'ET' MatConv ; 679 : 'DETRUIT' MatConv; 680 : 'MENAGE' ; 681 : 'FINSI' ; 682 : 'FINSI' ; 683 : 684 : 685 : 686 : *--------------------------------------------------------------------- 687 : *------ Conditions aux limites de Flux (mixtes et Neumann) ----------- 688 : *--------------------------------------------------------------------- 689 : 690 : 691 : 'SI' (FLUCLI) ; 692 : SMTr = SMTR 'ET' CLFLUX ; 693 : 'FINSI' ; 694 : 695 : 'MENAGE' ; 696 : 697 : 'FINP' SMTr MatrTr TbDarTra MassEFMH nomespec 698 : nbespece nbsource Diffdisp Tcini TABRES TABMODI; 699 : 700 : 701 : 702 :
© Cast3M 2003 - All rights reserved.
Disclaimer