1 : * INITVF PROCEDUR GOUNAND 11/05/24 21:15:27 6976 2 : ********************************************************************** 3 : 'DEBP' INITVF 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 : CHCLIM*'TABLE' OPTRESOL/'TABLE' LCONV*'LOGIQUE' ; 9 : * ATTENTION La vitesse est optionnelle, L'ordre est important 10 : * et les types d'arguments qui se suivent aussi pour tester leur 11 : * présence 12 : * 13 : * Attention il faudra transformer les vitesses en débits face 14 : * et sortir le champ 15 : * 16 : * |-----------------------------------------------------------------| 17 : * | Phrase d'appel (en GIBIANE) | 18 : * |-----------------------------------------------------------------| 19 : * | | 20 : * | RESI Matot jaco Mpor2 Mchamt mchamt1 difftot nomespec nbespece | 21 : * | nbsource TABRES TABMODI NOUVMAT = INITVF MoDARCY Porosite | 22 : * | MateDiff ChPSour DeltaT Cini TetaDiff | 23 : * | TetaConv TetaLin (QFACE) QELEM | 24 : * | DISPL DISPT CHCLIM optresol ; | 25 : * | | 26 : * |-----------------------------------------------------------------| 27 : * | Généralités : INITVF construit la matrice de discrétisation | 28 : * | du problème de transport convection-diffusion pour| 29 : * | le premier pas de tps d'un algorithme transitoire.| 30 : * | Le second membre et les Conditions limites de flux| 31 : * | sont pris en compte. | 32 : * | RESTE TCINI, DECENTR et TERME LIN | 33 : * |-----------------------------------------------------------------| 34 : * | | 35 : * |-----------------------------------------------------------------| 36 : * | ENTREES | 37 : * |-----------------------------------------------------------------| 38 : * | MoDARCY : modele Darcy. | 39 : * | | 40 : * | Porosite : champ par elements de composante 'CK' | 41 : * | | 42 : * | MateDiff : Tenseur de diffusion (type iso, ..) champ par | 43 : * | elements de composante 'K' en isotrope, 'K11', 'K21',| 44 : * | 'K22' en anisotrope 2d et 'K11', 'K21', 'K22', 'K31'| 45 : * | 'K32', 'K33' en anisotrope 3d. Type 'CARACTERISTIQ 46 : * | | 47 : * | ChPSour : Champ par points des sources volumiques par unité de | 48 : * | temps (support maillage centre). Composante ?????? | 49 : * | | 50 : * | DeltaT : Pas de temps. | 51 : * | | 52 : * | Cini : Concentration initiale, CHPOINT centre. | 53 : * | | 54 : * | TetaDiff : Valeur de theta pour theta-schéma en temps, opérateur| 55 : * | de diffusion. Entre 0 et 1. 0 = explicite, 1 = euler | 56 : * | implicite. | 57 : * | | 58 : * | TetaConv : Valeur de theta pour theta-schéma en temps, opérateur| 59 : * | de convection. Entre 0 et 1. 0 = explicite, 1 = Euler| 60 : * | implicite. | 61 : * | | 62 : * | TetaLin : valeur de theta pour theta-schéma en temps, opérateur| 63 : * | linéaire du type coef * C, où C est l'inconnue. | 64 : * | Entre 0 et 1. 0 = explicite, 1 = euler implicite. | 65 : * | ??????????? A voir car peut etre identique à Tetadiff| 66 : * | | 67 : * | | 68 : * | Qface : vitesse aux faces, CHPO face de composantes Vx, Vy | 69 : * | en 2d et Vx, Vy, Vz en 3d. Il s'agit plus exatement | 70 : * | de (V.n)n, c'est à dire de la composante normale de | 71 : * | la vitesse aux faces. ???????? (je pressens que | 72 : * | castem va sortir des flux, cad intégrés sur surfaces)| 73 : * | | 74 : * | CHCLIM : table d'indice 'NEUMANN' et 'DIRICHLET' contenant les| 75 : * | Chpoint à n composantes contenant les conditions aux | 76 : * | limites de Neumann et Dirichlet par espece. | 77 : * | | 78 : * | OPTRESOL : Table dont l'entree est optionnelle définissant | 79 : * | les options de résolution pour 'KRES'. | 80 : * | | 81 : * |-----------------------------------------------------------------| 82 : * | SORTIES | 83 : * |-----------------------------------------------------------------| 84 : * | | 85 : * | | 86 : * | RESI : second membre | 87 : * | | 88 : * | Matot : matrice globale de discretisation en VF | 89 : * | | 90 : * | Jaco : matrice globale de discretisation en VF pour le problème 91 : * | stationnaire | 92 : * | | 93 : * | Mpor : matrice globale de discretisation en VF pour le problème 94 : * | stationnaire | 95 : * | | 96 : * | Mchamt : Coef permettant de calculer le flux total | 97 : * | | 98 : * | Mchamt1 : Coef permettant de calculer le flux diffusif | 99 : * | | 100 : * | Difftot : Coefficient de diffusion totale, integre decentrement| 101 : * | | 102 : * | Diffdisp : Coefficient de dispersivité | 103 : * | | 104 : * | nomespc : liste des noms de composante des espèces dans Cini | 105 : * | | 106 : * | nbespece : nombre de composante de Cini, soit nombre d'especes | 107 : * | | 108 : * | nbsource : nombre de composantes du terme source qd X especes | 109 : * | | 110 : * | TABRES : Table complète définissant les options de résolution | 111 : * | pour 'KRES'. | 112 : * | | 113 : * | TABMODI : table contenant des logiques indiquant la nécessité | 114 : * | ou non de reclalculer certains termes. | 115 : * | 'POROSITE' : VRAI si le coefficient devant D/DT | 116 : * | (porosité) est modifié depuis le dernier| 117 : * | appel | 118 : * | 'DELTAT' : VRAI si le pas de tps a changé | 119 : * | 'CONVECTI' : VRAI si la vitesse a changé | 120 : * | 'COEF_LIN' : VRAI si le coef en facteur de C a changé| 121 : * | 'DIFFUSI' : VRAI si les diffusivités ont changé | 122 : * | | 123 : * | NOUVMAT : Logique affecté à VRAI lorsque que Matot vient 124 : * | d'etre calculée 125 : * |-----------------------------------------------------------------| 126 : * | VARIABLES INTERNES | 127 : * |-----------------------------------------------------------------| 128 : * | | 129 : * | PCONV : Logique indiquant VRAI si présence de convection | 130 : * | | 131 : * | toltheta : 1.D-4 seuil en dessous duquel on considère que la | 132 : * | valeur de theta du theta-schéma est nulle (schéma | 133 : * | explicite) ou au contraire euler-implicite si | 134 : * | theta > 0.9999 | 135 : * | | 136 : * | SSource : Source aux centre (une composante) | 137 : * | | 138 : * | CCini : concentration aux centres (une composante) | 139 : * | | 140 : * | lstcps : liste des noms de composante des espèces dans Chpsour| 141 : * | | 142 : * | SSMTr : second membre sur les traces pour une espèce | 143 : * | | 144 : * | FLUNEU : LOGIQUE valant VRAI si conditions de Neumann | 145 : * | | 146 : * | CLFLUX : Chpoint à n composantes contenant les flux imposés | 147 : * | pour chaque espece chimique. nul si pas de flux | 148 : * | OPTIONNEL | 149 : * | | 150 : ********************************************************************** 151 : 152 : *--------------------------------------------------------------------- 153 : *---------- On récupere les conditions limites ------------------ 154 : *--------------------------------------------------------------------- 155 : 156 : *MESS 'MODARCY'; 157 : *LIST MODARCY; 158 : FLUNEU = FAUX; 159 : DIRCLI = FAUX; 160 : FLUTOT = FAUX; 161 : FLUMIX = FAUX; 162 : CLFLUX MA = 'KOPS' 'MATRIK'; 163 : CLDIRI MA = 'KOPS' 'MATRIK'; 164 : CLFLUT MA = 'KOPS' 'MATRIK'; 165 : CLFLUX3 MA = 'KOPS' 'MATRIK'; 166 : EXFLU = CLFLUX; 167 : EXDIR = CLDIRI; 168 : EXFLUT = CLFLUT; 169 : EXFLUX3 = CLFLUX3; 170 : 171 : * Neumann 172 : 'SI' ('EXISTE' CHCLIM 'NEUMANN') ; 173 : CLFLUX = CHCLIM . 'NEUMANN'; 174 : FLUNEU = VRAI; 175 : 'FINSI'; 176 : 177 : 'SI' ('EXISTE' CHCLIM 'DIRICHLET') ; 178 : CLDIRI = CHCLIM . 'DIRICHLET'; 179 : DIRCLI = VRAI; 180 : 'FINSI'; 181 : 182 : 'SI' ('EXISTE' CHCLIM 'FLUTOTAL') ; 183 : CLFLT = CHCLIM . 'FLUTOTAL'; 184 : FLUTOT = VRAI; 185 : 'FINSI'; 186 : 187 : * Flux mixte 188 : 'SI' ('EXISTE' CHCLIM 'FLUMIXTE') ; 189 : * comme on impose A Dgrad C + B C = flumix, on le traite sous 190 : * la forme D grad C + (B/A) C = flumix/A plus naturelle en EFMH car 191 : * D grad C est le flux diffusif 192 : COFA = -1.D0 * CHCLIM . 'FLUMIXTE' . 'COEFA' ; 193 : CLFLUX3 = CHCLIM . 'FLUMIXTE' . 'VAL' '/' COFA ; 194 : CLFLUX3 = CHAN 'ATTRIBUT' CLFLUX3 NATURE DISCRET ; 195 : FLUMIX = VRAI ; 196 : mayage = 'EXTRAIRE' CHCLIM . 'FLUMIXTE' . 'VAL' maillage ; 197 : cofb = (doma modarcy SURFACE) *CHCLIM . 'FLUMIXTE' . 'COEFB' ; 198 : cofb = CHCLIM . 'FLUMIXTE' . 'COEFB' ; 199 : cofa = 'REDU' (CHCLIM . 'FLUMIXTE' . 'COEFA') mayage ; 200 : cofb = 'REDU' cofb mayage ; 201 : coefm = (1.D0 * cofb) '/' cofa ; 202 : CLFLUX3 = 'REDU' CLFLUX3 mayage ; 203 : 'OUBLIER' cofa ; 204 : 'OUBLIER' cofb ; 205 : 'FINSI' ; 206 : *MESS 'DIRICHLET'; 207 : *LIST CLDIRI ; 208 : *--------------------------------------------------------------------- 209 : *---------- Initialisations de tables, coefficients ------------------ 210 : *--------------------------------------------------------------------- 211 : 212 : * Table de logiques indiquant des modifications. Initialisation 213 : TABMODI = TABLE; 214 : TABMODI . 'POROSITE' = FAUX; 215 : TABMODI . 'CONVECTI' = FAUX; 216 : TABMODI . 'DELTAT' = FAUX; 217 : TABMODI . 'COEF_LIN' = FAUX; 218 : TABMODI . 'DIFFUSIV' = FAUX; 219 : 220 : 221 : 'SI' ('EXISTE' QFACE) ; 222 : PCONV = VRAI; 223 : 'SI' ('EXISTE' DISPL) ; 224 : DISPERSI = VRAI ; 225 : 'SINON' ; 226 : DISPERSI = FAUX; 227 : 'FINSI' ; 228 : 'SINON' ; 229 : PCONV = FAUX; 230 : DISPERSI = FAUX; 231 : 'FINSI' ; 232 : 233 : 234 : 235 : * tolerance sur theta du theta schéma de discrétisation en temps. 236 : * il faudrait remmettre les theta à 0 ou 1 si nécessaire dans 237 : * procédure amont. 238 : toltheta = 1.D-4; 239 : 240 : 'SI' (TetaConv 'EGA' 0.D0 toltheta) ; 241 : TetaConv = 0.D0; 242 : 'FINSI' ; 243 : 'SI' (TetaConv 'EGA' 1.D0 toltheta) ; 244 : TetaConv = 1.D0; 245 : 'FINSI' ; 246 : 'SI' (TetaDiff 'EGA' 0.D0 toltheta) ; 247 : TetaDiff = 0.D0; 248 : 'FINSI' ; 249 : 'SI' (TetaDiff 'EGA' 1.D0 toltheta) ; 250 : TetaDiff = 1.D0; 251 : 'FINSI' ; 252 : 'SI' (TetaLin 'EGA' 0.D0 toltheta) ; 253 : TetaLin = 0.D0; 254 : 'FINSI' ; 255 : 'SI' (TetaLin 'EGA' 1.D0 toltheta) ; 256 : TetaLin = 1.D0; 257 : 'FINSI' ; 258 : 259 : 260 : * Calcul du terme devant le dC/dt integré sur le volume 261 : *MPOR2 = 'KOPS' 'MATDIAGO' ('MOTS' 'RETN') Porosite ; 262 : MPOR2 = 'KOPS' 'MATDIAGO' ('NOMC' 'RETN' Porosite) 'MATRIK' ; 263 : 264 : 265 : 266 : *--------------------------------------------------------------------- 267 : *----------------------- CREATION TABLE POUR RESOLUTION -------------- 268 : *--------------------------------------------------------------------- 269 : 270 : 271 : **************** OPTIONS PAR DEFAUT ************************** 272 : 273 : * création de la table de résolution pour l'opérateur KRES 274 : * On crée la table de résolution avec les options par défaut 275 : * On y remplacera les valeurs définit par l'utilisateur ensuite. 276 : TABRES = 'TABLE' 'METHINV' ; 277 : 278 : * option BCGSTAB par défaut pour une matrice non symétrique 279 : METHRES = 3; 280 : TABRES . 'TYPINV' = METHRES ; 281 : 282 : * niveau d'impression. 283 : TABRES . 'IMPINV' = 0 ; 284 : 285 : * Type de renumérotation. Option SLOANE par défaut 286 : TABRES . 'TYRENU' = 'SLOANE' ; 287 : 288 : * La gestion des multiplicateurs sera modifiée 289 : * par la suite. Pas d'option pour l'instant 290 : TABRES . 'PCMLAG' = 'APR2' ; 291 : TABRES . 'OUBMAT' = 0 ; 292 : TABRES . 'SCALING' = 0 ; 293 : 294 : *INDICES SPÉCIFIQUES POUR UNE MÉTHODE ITÉRATIVE 295 : * Nombre maxi d'itérations 296 : TABRES . 'NITMAX' = 1500 ; 297 : * résidu pour la convergence de la méthode 298 : TABRES . 'RESID' = 1.D-15 ; 299 : * valeur minimale du pivot de la méthode 300 : TABRES . 'BCGSBTOL' = 1.D-120 ; 301 : * preconditionnement ILU(0) 302 : TABRES . 'PRECOND' = 3 ; 303 : *relaxation pour MILU0 304 : TABRES . 'MILURELX' = 1.D0 ; 305 : *GMRESTART 306 : TABRES . 'GMRESTRT' = 100 ; 307 : *ILUTLFIL 308 : TABRES . 'ILUTLFIL' = 2; 309 : * drop tolerence pour ILUT2 310 : TABRES . 'ILUTDTOL' = 0.D0; 311 : TABRES . 'ILUTPTOL' = 0.01D0; 312 : TABRES . 'ILUTALPH' = 0.D0 ; 313 : 314 : ************** OPTIONS UTILISATEUR ************************** 315 : 316 : 317 : 'SI' ('EGA' ('TYPE' OPTRESOL) 'TABLE') ; 318 : * L'utilisateur a défini des options pour la méthode 319 : * de résolution. 320 : 321 : * Type d'inversion 322 : 'SI' ('EXISTE' OPTRESOL 'TYPINV') ; 323 : TABRES . 'TYPINV' = OPTRESOL . 'TYPINV'; 324 : 'FINSI' ; 325 : 326 : * Niveau d'impression 327 : 'SI' ('EXISTE' OPTRESOL 'IMPINV') ; 328 : TABRES . 'IMPINV' = OPTRESOL . 'IMPINV'; 329 : 'FINSI' ; 330 : 331 : * Type de renumérotation. 332 : 'SI' ('EXISTE' OPTRESOL 'TYRENU') ; 333 : TABRES . 'TYRENU' = OPTRESOL . 'TYRENU'; 334 : 'FINSI' ; 335 : 336 : * Indices spécifiques aux méthodes itératives 337 : 'SI' ((TABRES . 'TYPINV') > 1); 338 : * Nombre maxi d'iterations 339 : 'SI' ('EXISTE' OPTRESOL 'NITMAX') ; 340 : TABRES . 'NITMAX' = OPTRESOL . 'NITMAX'; 341 : 'FINSI' ; 342 : * Valeur du résidu de la méthode 343 : 'SI' ('EXISTE' OPTRESOL 'RESID') ; 344 : TABRES . 'RESID' = OPTRESOL . 'RESID'; 345 : 'FINSI' ; 346 : * valeur minimal du pivot de la méthode 347 : 'SI' ('EXISTE' OPTRESOL 'BCGSBTOL') ; 348 : TABRES . 'BCGSBTOL' = OPTRESOL . 'BCGSBTOL'; 349 : 'FINSI' ; 350 : * precond par diagonale 351 : 'SI' ('EXISTE' OPTRESOL 'PRECOND') ; 352 : TABRES . 'PRECOND' = OPTRESOL . 'PRECOND'; 353 : 'FINSI' ; 354 : * precon ILUT2 355 : 'SI' ('EXISTE' OPTRESOL 'ILUTLFIL') ; 356 : TABRES . 'ILUTLFIL' = OPTRESOL . 'ILUTLFIL' ; 357 : 'FINSI' ; 358 : 'FINSI' ; 359 : * Pour GMRES 360 : 'SI' ((TABRES . 'TYPINV') EGA 5); 361 : 'SI' ('EXISTE' OPTRESOL 'GMRESTRT') ; 362 : TABRES . 'GMRESTRT' = OPTRESOL . 'GMRESTRT'; 363 : 'SINON' ; 364 : TABRES . 'GMRESTRT' = 50; 365 : 'FINSI' ; 366 : 'FINSI' ; 367 : 'FINSI' ; 368 : 369 : 370 : 371 : SI (('EGA' TABRES . 'PRECOND' 8) 'OU' ('EGA' TABRES . 'PRECOND' 7)); 372 : TABRES . 'ILUTDTOL' = 0.1D-2; 373 : FINSI ; 374 : 375 : *--------------------------------------------------------------------- 376 : *--------------------- CALCUL DE LA DISPERSIVITE---------------------- 377 : *--------------------------------------------------------------------- 378 : 379 : * On remmet la matrice de diffusivité avec les bons noms de composantes 380 : * anisotropes 381 : 382 : zozo = DIFFANIS Matediff ; 383 : Matediff = zozo ; 384 : 385 : * Seulement si présence de convection et présence de 386 : * dispersivité et disp_l et disp_t 387 : 388 : 'SI' (DISPERSI 'ET' PCONV) ; 389 : Diffdisp = CALCDISP QELEM DISPL DISPT ; 390 : 'SINON' ; 391 : Diffdisp = 0.D0 * Matediff ; 392 : 'FINSI' ; 393 : 394 : Difftot = Matediff '+' Diffdisp ; 395 : 396 : 397 : *---------------------------------------------------------------------- 398 : *--- CALCUL DE LA MATRICE DE CONVECTION ET DES SECONDS MEMBRES -------- 399 : *---------------------------------------------------------------------- 400 : 401 : 402 : ************** Nombres d'espèces à gérer ****************************** 403 : 404 : ** évaluation du nombre d'espèces. 405 : nomespec = 'EXTRAIRE' Cini 'COMP'; 406 : nbespece = 'DIME' nomespec; 407 : * nb de termes sources, commun (1), ou egal au nombre 408 : * d'especes (nbcompi) 409 : lstcps = 'EXTRAIRE' Chpsour 'COMP'; 410 : nbsource = 'DIME' lstcps; 411 : 'SI' ((nbsource 'NEG' 1) 'ET' (nbsource 'NEG' nbespece)) ; 412 : 'MESSAGE' 'La source doit avoir le meme nombre de composantes'; 413 : 'MESSAGE' 'que les espèces ou 1 seule composante'; 414 : ERREUR 5; 415 : 'FINSI' ; 416 : 417 : 418 : ********************* Une ou plusieurs espèces ************************ 419 : 420 : 'SI' (nbespece > 0) ; 421 : 422 : 423 : DTI = 1.D0/DeltaT; 424 : 'REPETER' bloc1 nbespece; 425 : * On extrait la composante de Cini, Tcini et de la source 426 : CCini = 'NOMC' 'H' ('EXCO' (extr &bloc1 nomespec) Cini); 427 : 'SI' (nbsource > 1) ; 428 : SSource = 'NOMC' 'SOUR' ('EXCO' (extr &bloc1 nomespec) 429 : Chpsour); 430 : 'SINON' ; 431 : SSource = 'NOMC' 'SOUR' Chpsour; 432 : 'FINSI' ; 433 : * Conditions initiales 434 : HHS = NOMC (CCini) 'SCAL' ; 435 : 436 : 437 : 438 : * Extraction des conditions aux limites 439 : 'SI' (FLUNEU) ; 440 : EXFLU = 441 : ('NOMC' 'FLUX' ('EXCO' ('EXTRAIRE' &bloc1 nomespec) 442 : CLFLUX)) ; 443 : * MESS 'EXFLU ='; LIST EXFLU; 444 : 'FINSI' ; 445 : 'SI' (DIRCLI) ; 446 : EXDIR = 447 : ('NOMC' 'SCAL' ('EXCO' ('EXTRAIRE' &bloc1 nomespec) 448 : CLDIRI)) ; 449 : * MESS 'EXDIR ='; LIST EXDIR; 450 : 'FINSI' ; 451 : 452 : 'SI' (FLUTOT) ; 453 : EXFLUT = 454 : ('NOMC' 'FLUX' ('EXCO' ('EXTRAIRE' &bloc1 nomespec) 455 : CLFLT )) ; 456 : * MESS 'QLIM ='; 457 : QLIM = NOMC(EXFLUT) 'FLUX' ; 458 : SUP = 'EXTR' EXFLUT 'MAIL' ; 459 : XPAR1 = 1.D0 + (0.0*CLFLT) ; 460 : * MESS 'XPAR1 ='; 461 : XPAR1 = 'NOMC' XPAR1 'PAR1' ; 462 : 'SI' LCONV ; 463 : USCNR = 'REDU' QFACE SUP ; 464 : MUSCN = USCNR*(-1.D0) ; 465 : * MESS 'XPAR2 ='; 466 : 'SINON' ; 467 : MUSCN = (0.D0*CLFLT) ; 468 : 'FINSI' ; 469 : XPAR2 = 'NOMC' MUSCN 'PAR2' ; 470 : EXFLUT = XPAR1 + XPAR2 + QLIM ; 471 : EXFLUT = CHAN 'ATTRIBUT' EXFLUT NATURE DISCRET ; 472 : * MESS 'EXFLUT ='; LIST EXFLUT; 473 : 'FINSI' ; 474 : 475 : 'SI' (FLUMIX) ; 476 : EXFLUM = 477 : ('NOMC' 'FLUX' ('EXCO' ('EXTRAIRE' &bloc1 nomespec) 478 : CLFLUX3 )) ; 479 : QLIM = NOMC(EXFLUM) 'FLUX' ; 480 : SUP = 'EXTR' EXFLUM 'MAIL' ; 481 : XPAR1 = 1.D0 + (0.0*CLFLUX3) ; 482 : * MESS 'XPAR1 ='; 483 : XPAR1 = 'NOMC' XPAR1 'PAR1' ; 484 : MUSCN = coefm; 485 : XPAR2 = 'NOMC' MUSCN 'PAR2' ; 486 : EXFLUM = XPAR1 + XPAR2 + QLIM ; 487 : EXFLUM = CHAN 'ATTRIBUT' EXFLUM NATURE DISCRET ; 488 : EXFLUT = EXFLUT 'ET' EXFLUM ; 489 : 'FINSI' ; 490 : 491 : * calcul des coefs mchamt1 et mchamt qui permettent d'assembler la 492 : * matrice et de calculer les flux diffusifs et convectifs 493 : * On recalcule deux fois les coefficients. En général, cela doit etre 494 : * négligeable en temps de calcul; 495 : 496 : GRADT0 MCHAMT1 = 'PENT' MoDARCY 'FACE' 'MPFA' 497 : HHS 'DISPDIF' difftot 'TIMP' EXDIR 'QIMP' EXFLU 498 : 'MIXT' EXFLUT ; 499 : MCHAMT = MCHAMT1; 500 : MCHAJA = TetaDiff*MCHAMT1 ; 501 : 'SI' LCONV ; 502 : GRADT0 MCHAMT = 'PENT' MoDARCY 'FACE' 'MPFA' 503 : HHS 'DISPDIF' difftot 'TIMP' EXDIR 'QIMP' EXFLU 504 : 'MIXT' EXFLUT 'UPWICENT' QFace ; 505 : * ATTENTION PEUT ETRE FAUX SI TetaConv different de TetaDiff 506 : MCHCONV = MCHAMT + ((-1.D0)*MCHAMT1) ; 507 : MCHAJA = (TetaConv*MCHCONV) + (TetaDiff*MCHAMT1) ; 508 : * VERRUE 509 : MCHAJA = tetadiff*MCHAMT; 510 : 'DETRUIT' MCHCONV ; 511 : 'FINSI' ; 512 : 513 : 514 : * On ne calcule la jacobienne que pour la première espèce 515 : 'SI' ((&bloc1) 'EGA' 1) ; 516 : JACO CHPRES DT = 'LAPN' 'VF' 'CLAUDEIS' 'IMPL' 517 : MoDARCY HHS GRADT0 MCHAJA 'QIMP' EXFLU 'MIXT' EXFLUT 518 : 'TIMP' EXDIR ; 519 : MATOT = 'KOPS' 'MULT' -1.000D0 JACO ; 520 : MATI = '*' MPOR2 DTI ; 521 : MATOT = MATOT 'ET' MATI ; 522 : 'SINON' ; 523 : JACOB CHPRES DT = 'LAPN' 'VF' 'CLAUDEIS' 'EXPL' 524 : MoDARCY HHS GRADT0 ; 525 : 'FINSI'; 526 : * On reconstitue un champ de second membre 527 : 'SI' ((&bloc1) 'EGA' 1) ; 528 : SSMTr = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) CHPRES ; 529 : RESI = SSMTr ; 530 : 'SINON' ; 531 : CHPRES = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) 532 : ('COPIER' CHPRES) ; 533 : RESI = CHPRES 'ET' RESI ; 534 : 'FINSI' ; 535 : 'FIN' bloc1; 536 : 'FINSI' ; 537 : 538 : 'MENAGE' ; 539 : 540 : NOUVMAT = VRAI; 541 : 'FINP' RESI MATOT JACO MPOR2 MCHAMT MCHAMT1 diffdisp difftot nomespec 542 : nbespece nbsource TABRES TABMODI NOUVMAT; 543 : 544 : 545 : 546 : 547 : 548 :
© Cast3M 2003 - All rights reserved.
Disclaimer