1 : * DARCYSAT PROCEDUR LEPOTIER 08/12/10 21:15:20 6214 2 : * GBM REGLER LE PROBLEME DES TRACES - VF EFMH 3 : 4 : * DARCYSAT PROCEDUR MAUGIS 02/12/19 21:15:02 4527 5 : *********************************************************************** 6 : 'DEBPROC' DARCYSAT SATUR*'TABLE' ; 7 : 8 : NDIME = 'VALE' DIME ; 9 : * MESSAGE 10 : * Par défaut, on affiche beaucoup d'information 11 : 'SI' ('EXISTE' SATUR 'MESSAGES' ) ; 12 : DEBUG = SATUR.'MESSAGES' > 0 ; 13 : 'SINON' ; 14 : DEBUG = VRAI ; 15 : 'FINSI' ; 16 : 17 : * 18 : ******************************************************************** 19 : ********************************************************************* 20 : ********************************************************************* 21 : * LECTURES * 22 : ********************************************************************* 23 : ********************************************************************* 24 : ********************************************************************* 25 : 26 : * 27 : * CONDITIONS INITIALES OU DU DERNIER PAS CALCULE 28 : * 29 : TPSINI CHRG TRH FLU0 LAST1 = SATUTILS LICODINI SATUR ; 30 : * 31 : * DONNEES PHYSIQUES, GEOMETRIQUES ET MATERIELLES 32 : * 33 : MDMH MCHYB VHYB MAILSOM MAILCENT ABMC = SATUTILS LIMODELE SATUR ; 34 : * 35 : * LOI_SATURATION 36 : * 37 : LOIS = SATUTILS LILOISAT SATUR ; 38 : * 39 : * LOI_PERMEABILITE 40 : * 41 : LOIP = SATUTILS LILOIPER SATUR ; 42 : * 43 : * param physiques 44 : * 45 : COFEMMAG L_GRAV RHOWG DAXE CONVH = SATUTILS LIPHYSIK SATUR ; 46 : * 47 : * RECUPERATION DES DONNEES NUMERIQUES 48 : * 49 : NIVEAU TETA NPAS0 NITER0 ERR0 CFL0 ITMAXI OKPENAL 50 : cofpenal npenaldt cofdiv nmaxdt = SATUTILS LIPARANU SATUR ; 51 : 52 : NIVEAU = 'MOT' NIVEAU; 53 : 54 : * TEMPS_CALCULES 55 : isauv = LAST1 ; 56 : DCAL LTCALCUL DTAUTO ICAL tmin NPAS0 DELTAT TFINAL isor tpsor dsor 57 : = SATUTILS LITABTPS SATUR TPSINI debug NPAS0 ; 58 : * 59 : * initialisations pour TRANGEOL 60 : * 61 : CHCLIM GEOL1 TABMODI = SATUTILS YNYTYAL SATUR TETA ; 62 : * 63 : * integrale des chargements 64 : * 65 : startflu startmix startsou FLUIMP FLUMIX TERSOU = SATUTILS YNTGFLUX 66 : SATUR ; 67 : 68 : 69 : 70 : 71 : 72 : 73 : 74 : 75 : 76 : 77 : ******************************************************************* 78 : * PARAMETRES DIVERS PRECALCULES 79 : ******************************************************************* 80 : * Coordonnées à prendre en compte pour la gravité 81 : * (si AXE_G est nul, ZGRAV aussi, et donc pas d'effet de gravité). 82 : XAXE = 'COORDONNEE' 1 DAXE; 83 : YAXE = 'COORDONNEE' 2 DAXE; 84 : 'SI' (NDIME 'EGA' 2) ; 85 : XCO YCO = 'COOR' MAILSOM ; 86 : ZFF = 'KCHT' MDMH 'SCAL' 'SOMMET' 'COMP' 'SCAL' 87 : ((XCO * XAXE) + (YCO * YAXE)); 88 : XCO YCO = 'COOR' MAILCENT ; 89 : ZCC = 'KCHT' MDMH 'SCAL' 'CENTRE' 'COMP' 'SCAL' 90 : ((XCO * XAXE) + (YCO * YAXE)) ; 91 : 'SINON' ; 92 : ZAXE = 'COORDONNEE' 3 DAXE; 93 : XCO YCO ZCO = 'COOR' MAILSOM ; 94 : ZFF = 'KCHT' MDMH 'SCAL' 'SOMMET' 'COMP' 'SCAL' 95 : ((XCO * XAXE) + (YCO * YAXE) + (ZCO * ZAXE)); 96 : XCO YCO ZCO = 'COOR' MAILCENT ; 97 : ZCC = 'KCHT' MDMH 'SCAL' 'CENTRE' 'COMP' 'SCAL' 98 : ((XCO * XAXE) + (YCO * YAXE) + (ZCO * ZAXE)) ; 99 : 'FINSI' ; 100 : 101 : 102 : *--------------------------------------------------------------------* 103 : * RAPPEL DES DONNEES ENTREES DANS LA PROCEDURE * 104 : *--------------------------------------------------------------------* 105 : 106 : SATUTILS rekapitu SATUR L_GRAV rhowg DAXE convh cfl0 teta LAST1 debug ; 107 : 108 : 109 : *********************************************************************** 110 : *********************************************************************** 111 : *********************************************************************** 112 : * RESOLUTION * 113 : *********************************************************************** 114 : *********************************************************************** 115 : *********************************************************************** 116 : 117 : 118 : * 119 : *--------------------------------------------------------------------* 120 : * BOUCLE RESOLVANT LE SYSTEME POUR CHAQUE PAS DE TEMPS * 121 : *--------------------------------------------------------------------* 122 : * 123 : * Initialisations : 124 : * ================= 125 : * gbm a bugué sur PENAL à FAUX 126 : PENAL = OKPENAL ; 127 : TPS = TPSINI + DELTAT ; 128 : * 129 : *-- Paramètres physiques 130 : 131 : *- Pression tronquée à utiliser dans les lois de comportement : 132 : 133 : * tronkature pression et chgt 'unités 134 : PNS = SATUTILS TRONKP MDMH CHRG TRH L_GRAV NIVEAU RHOWG CONVH ; 135 : PNSC = SATUTILS TRONKP MDMH CHRG TRH L_GRAV 'CENTRE' RHOWG CONVH ; 136 : 137 : *- calcul teneur en eau, saturation et premier terme capacité capillaire 138 : * Attention, la teneur calculée utilise une porosité constante au cours 139 : * du pas de temps. GBM refait à l'identique dans boucle plus loin. 140 : * GBM - ATTENTION PEUT ETRE INUTIL ICI SAUF INITIALISATION FLUX 141 : 'SI' ( 'EXISTE' SATUR 'CONSERVATIF' ) ; 142 : CONSER = SATUR .'CONSERVATIF' ; 143 : 'FINSI' ; 144 : 145 : 'SI' ('EGA' CONSER VRAI) ; 146 : nbds SAT TENN CAPA PORO1 = SATUTILS KALSAT MDMH LOIS 'CENTRE' PNSC ; 147 : nbds SAT TEN CAPAD PORO1 = SATUTILS KALSAT MDMH LOIS NIVEAU PNS ; 148 : 'SINON' ; 149 : 150 : nbds SAT TEN CAPA PORO1 = SATUTILS KALSAT MDMH LOIS NIVEAU PNS ; 151 : 'FINSI'; 152 : 153 : 154 : 155 : *- pas de temps initial automatique - rendre calcul CFL optionnel 156 : VVOL = VHYB * PORO1 ; 157 : DELTAT DTI = SATUTILS YNYTDT MDMH FLU0 MCHYB VVOL CFL0 debug DTAUTO 158 : DELTAT ; 159 : 160 : *- Valeur des variables au pas de temps précédent 161 : * GBM PAS TOUTES UTILES - A REVOIR 162 : TRHANC = TRH ; 163 : CHRGANC = CHRG ; 164 : FLUANC = FLU0 ; 165 : PNSANC = PNS ; 166 : PNSCANC = PNSC ; 167 : TENANC = TEN ; 168 : SATANC = SAT ; 169 : TPSANC = TPSINI ; 170 : PERFANC = 0.D0 ; 171 : 172 : *- Valeur des variables à l'itération précédente 173 : * pour la trace de charge servant au calcul du résidu, 174 : * on met n'importe quoi qui ait la bonne structure. 175 : * Ici, TRH n'a pas de multiplicateur de lagrange, donc l'estimation 176 : * de FLRES sera foireuse au tout prsavreemier calcul du résidu. 177 : TRH2N = TRH ; 178 : CAPAN = 0.D0 ; 179 : PERFN = 0.D0 ; 180 : 181 : NOMESPL = 'H' ; 182 : 183 : * 184 : * initialisation du terme source intégral. 185 : * 186 : 187 : 'SI' ( 'EXISTE' SATUR 'SOURCE' ) ; 188 : TERSC2M1 = 'NOMC' (NOMESPL) ('TIRE' TERSOU TPSINI) ; 189 : 'FINSI' ; 190 : 191 : 'SI' ('EXISTE' SATUR 'FLUX_IMPOSE') ; 192 : FLUIMPM1 = 'TIRE' FLUIMP TPSINI ; 193 : 'FINSI' ; 194 : 195 : 'SI' ('EXISTE' SATUR 'FLUMIXTE') ; 196 : FLUMMPM1 = 'TIRE' FLUMIX TPSINI ; 197 : 'FINSI' ; 198 : 199 : FLU1 = FLU0 ; 200 : GEOL1 . 'CONCENTRATION' = CHRG ; 201 : 202 : 'SI' ('EGA' CONSER VRAI) ; 203 : TENNC = TENN ; 204 : 'FINSI' ; 205 : 206 : ************** NPX*************************************************** 207 : * 208 : * Sauvegarde 209 : * ========== 210 : 211 : MESS 'SAUVEGARDE INITIALE---------------------------------------' ; 212 : 213 : 214 : 215 : TMP TMP2 = SATUTILS SAVRESU SATUR TPSOR ISOR TPSANC TPS -1 DELTAT 216 : TRHANC TRH CHRGANC CHRG TENANC TEN PNSCANC PNSC SATANC SAT FLUANC 217 : FLU0 NIVEAU DSOR 0 MDMH ; 218 : 219 : ***************FIN NPX*********************************************** 220 : 221 : *======================================================== transitoire 222 : 223 : 'REPETER' TRANSI NPAS0 ; 224 : 225 : *==================================================================== 226 : * 227 : 228 : IPAS = ICAL + &TRANSI - 1 ; 229 : 230 : 231 : *----------------------------------------------------------------- 232 : 233 : 'REPETER' PENALDT NPENALDT ; 234 : 235 : CHRG = CHRGANC ; 236 : 237 : TPS = TPSANC + DELTAT ; 238 : 239 : 240 : *------------------------------------------------------------------ 241 : * Affichage information si debug vrai 242 : 'MESSAGE' 'deltat dti' deltat dti; 243 : SATUTILS AFFICH &PENALDT &TRANSI LAST1 TPS TPSANC DELTAT DTI 244 : debug ; 245 : 246 : ***************** INITIALISATION PAS DE TPS ************************** 247 : * 248 : *- Incorporation des CLs 249 : * 250 : 'SI' ('EXISTE' SATUR 'TRACE_IMPOSE') ; 251 : CHARIMPO = 'TIRE' SATUR . 'TRACE_IMPOSE' TPS ; 252 : CHCLIM . 'DIRICHLET' = 'NOMC' NOMESPL CHARIMPO ; 253 : 'FINSI' ; 254 : 255 : 256 : 'SI' ('EXISTE' SATUR 'FLUX_IMPOSE') ; 257 : 'SI' (tpsanc '<EG' startflu) ; 258 : FLUIMPO = 'TIRE' FLUIMP TPS ; 259 : 'SINON' ; 260 : FLUIMPO = 'COPIER' FLUIMPM1 ; 261 : 'FINSI' ; 262 : FLUIMPO = 'CHAN' 'ATTRIBUT' FLUIMPO 'NATURE' 'DISCRET' ; 263 : 'SI' (('EGA' (SATUR . 'TYPDISCRETISATION') 'VF') 264 : 'OU' (('EGA' (SATUR . 'TYPDISCRETISATION') 'EFMH') 265 : 'ET' ('EGA' teta 1.D0))) ; 266 : * On écrit les flux sous forme intégrale sauf en explicite et 267 : * kranck-Nickholson pour EFMH. On les différentie ici 268 : FLUMP = (FLUIMPO '-' FLUIMPM1) '/' DELTAT ; 269 : 'SINON' ; 270 : FLUMP = 'COPIER' FLUIMPO ; 271 : 'FINSI' ; 272 : CHCLIM . 'NEUMANN' = 'NOMC' NOMESPL FLUMP ; 273 : 'FINSI' ; 274 : 275 : 'SI' ('EXISTE' SATUR 'FLUMIXTE') ; 276 : 'SI' (tpsanc '<EG' startmix) ; 277 : FLUMMPO = 'TIRE' FLUMIX TPS ; 278 : 'SINON' ; 279 : FLUMMPO = 'COPIER' FLUMMPM1 ; 280 : 'FINSI' ; 281 : FLUMMPO = 'CHAN' 'ATTRIBUT' FLUMMPO 'NATURE' 'DISCRET' ; 282 : 'SI' (('EGA' (SATUR . 'TYPDISCRETISATION') 'VF') 283 : 'OU' (('EGA' (SATUR . 'TYPDISCRETISATION') 'EFMH') 284 : 'ET' ('EGA' teta 1.D0))) ; 285 : * On écrit les flux sous forme intégrale sauf en explicite et 286 : * kranck-Nickholson pour EFMH. On les différentie ici 287 : FLUMMP = ( FLUMMPO '-' FLUMMPM1) '/' DELTAT ; 288 : 'SINON' ; 289 : FLUMMP = 'COPIER' FLUMMPO ; 290 : 'FINSI' ; 291 : CHCLIM . 'FLUMIXTE' = TABLE ; 292 : CHCLIM . 'FLUMIXTE' . 'VAL' = 'NOMC' NOMESPL FLUMMP ; 293 : CHCLIM . 'FLUMIXTE' . 'COEFA' = SATUR . 'FLUMIXTE' . 'MIXCOFA' ; 294 : * GBM rappel il y aura des modifs à faire car HYDRAU 295 : CHCLIM . 'FLUMIXTE' . 'COEFB' = SATUR . 'FLUMIXTE' . 'MIXCOFB' ; 296 : 'FINSI' ; 297 : 298 : GEOL1 . 'CLIMITES' = CHCLIM ; 299 : 300 : * 301 : *- Calcul de la contribution des termes sources 302 : * 303 : TERSCE = 'MANU' 'CHPO' MAILCENT 1 'SOUR' 0. ; 304 : TERSCE = 'NOMC' NOMESPL TERSCE ; 305 : 306 : * Terme source propre : 307 : 'SI' ( 'EXISTE' SATUR 'SOURCE' ) ; 308 : 'SI' (tpsans <EG startsou) ; 309 : * la source varie encore et est non nulle 310 : TERSC2 = 'TIRE' TERSOU TPS ; 311 : 'SINON' ; 312 : * la source est nulle on garde l'intégrale précédente 313 : TERSC2 = TERSC2M1 ; 314 : 'FINSI' ; 315 : TERSC2 = 'NOMC' NOMESPL TERSC2 ; 316 : TERSCV = (TERSC2 '-' TERSC2M1) '/' DELTAT ; 317 : TERSCE = TERSCV + TERSCE ; 318 : 'FINSI' ; 319 : 320 : * Le terme source physique est tersce. 321 : 322 : GEOL1 . 'DELTAT' = DELTAT ; 323 : GEOL1 . 'METHINV' = SATUR . 'METHINV' ; 324 : 325 : *----------------------------------------------------- non linéaire 326 : 327 : 'REPETER' LINEAR itmaxi ; 328 : 'SI' ('EGA' CONSER VRAI) ; 329 : GEOL1 . 'CONCENTRATION' = CHRG ; 330 : 'FINSI' ; 331 : 332 : *__________________________________________________________________ 333 : 334 : * Préparation relaxation - GBM tester efficacité 335 : * ====================== 336 : 'SI' ('EGA' &LINEAR 1) ; 337 : * pas de relaxation à la première itération 338 : cofrelax = 1.D0 ; 339 : 'SINON' ; 340 : 'SI' (&LINEAR '<EG' (itmaxi/2)) ; 341 : * on relaxe à partir du deuxième pas de temps 342 : cofrelax = SATUR.'SOUS_RELAXATION' ; 343 : 'SINON' ; 344 : * on relaxe plus si on dépasse itmaxi/2 345 : * GBM GBM GBMGBM 0.5 normalement 346 : cofrelax = 0.5D0 ; 347 : 'FINSI' ; 348 : 'FINSI' ; 349 : 350 : * Calcul des nouveaux paramètres 351 : * ============================== 352 : 353 : *-- Pression tronquée à utiliser dans les lois de comportement: 354 : * P1 pression pascal, PNS pression tronquée 355 : PNS = SATUTILS TRONKP MDMH CHRG TRH L_GRAV NIVEAU RHOWG CONVH ; 356 : PNSC = SATUTILS TRONKP MDMH CHRG TRH L_GRAV 'CENTRE' RHOWG CONVH ; 357 : 358 : 359 : * GBM reprendre calcul sat, ten et capa 360 : * mettre un flag pour calculer CAPA ou pas 361 : * retirer PNSANC !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 362 : * TENANC 363 : 364 : 'SI' ('EGA' CONSER VRAI) ; 365 : nbds SATC TENC CAPA = SATUTILS KALSOT MDMH LOIS 'CENTRE' PNSC ; 366 : dum1 SAT TEN CAPAD = SATUTILS KALSOT MDMH LOIS NIVEAU PNS ; 367 : 'SINON' ; 368 : nbds SATC TEN CAPA = SATUTILS KALSOT MDMH LOIS 'CENTRE' PNSC ; 369 : dum1 SAT TEN CAPA = SATUTILS KALSOT MDMH LOIS NIVEAU PNS ; 370 : 'FINSI' ; 371 : 372 : 373 : 374 : *-- Calcul perméabilité 375 : * PERF = 0. ; 376 : 'REPETER' BB ('DIME' LOIP.'INDEX') ; 377 : NOMT = 'CHAINE' LOIP.'INDEX'.&BB ; 378 : LOI1 = LOIP . NOMT ; 379 : MOD1 = LOI1.'MODELE' ; 380 : NOMKR = 'MOT' ('TEXTE' ('CHAINE' LOI1.'NOM_PROCEDURE')) ; 381 : 382 : 'SI' (('EGA' LOI1.'SOUSTYPE' ('CHAINE' PUISSANCE) ) 'OU' 383 : ('EGA' LOI1.'SOUSTYPE' ('CHAINE' PERSONNELLE)) 'OU' 384 : ('EGA' LOI1.'SOUSTYPE' ('CHAINE' MUALEM)) 'OU' 385 : ('EGA' LOI1.'SOUSTYPE' ('CHAINE' BURDINE)) 'OU' 386 : ('EGA' LOI1.'SOUSTYPE' ('CHAINE' BROOKS_COREY)) 'OU' 387 : ('EGA' LOI1.'SOUSTYPE' ('CHAINE' MUALEM_BURDINE))) ; 388 : * réduction de la saturation au sous-domaine (local) 389 : SATL = 'REDU' SAT ('DOMA' MOD1 NIVEAU) ; 390 : PERFL = ('TEXTE' NOMKR) LOI1 SATL ; 391 : PERFL = RECENTRE PERFL MOD1 NIVEAU ; 392 : 'SINON' ; 393 : PNL = 'REDU' PNS ('DOMA' MOD1 NIVEAU) ; 394 : PERFL = ('TEXTE' NOMKR) LOI1 PNL ; 395 : PERFL = RECENTRE PERFL MOD1 NIVEAU ; 396 : 'FINSI' ; 397 : 398 : * 'SI' ('EGA' NIVEAU 'FACE') ; 399 : * LOI1.'ABMC' = 'ABS' ('DOMA' MOD1 'ORIENTAT') ; 400 : * 'FINSI' ; 401 : 402 : 403 : * concaténation perméabilité 404 : LCOMP = 'EXTR' PERFL 'COMP' ; 405 : NCOMP = 'DIME' LCOMP ; 406 : 'REPETER' CONSPER NCOMP ; 407 : J = &CONSPER ; 408 : COMPI = 'EXTR' LCOMP J ; 409 : PERFINI = 'MANU' 'CHPO' ('DOMA' MDMH 'CENTRE') 1 COMPI 0. 410 : 'NATURE' 'DISCRET' ; 411 : PERFI = 'EXCO' PERFL COMPI COMPI ; 412 : PERFI = 'KCHT' MDMH 'SCAL' 'CENTRE' 'COMP' COMPI PERFINI PERFI; 413 : * 414 : 'SI' (&BB 'EGA' 1) ; 415 : 'SI' ('EGA' J 1) ; 416 : PERFT = PERFI ; 417 : 'SINON' ; 418 : PERFT = PERFT 'ET' PERFI ; 419 : 'FINSI' ; 420 : 'SINON' ; 421 : 'SI' ('EGA' J 1) ; 422 : PERFT = PERFT 'ET' PERFI ; 423 : 'SINON' ; 424 : PERFT = PERFT 'ET' PERFI ; 425 : 'FINSI' ; 426 : 'FINSI' ; 427 : * 428 : 'FIN' CONSPER ; 429 : * 430 : 'OUBLIER' satl ; 431 : 'OUBLIER' perl; 432 : 433 : 'FIN' BB ; 434 : 435 : *-- coefficient d'emmagasinement si milieu saturé 436 : * GBM ATTENTION ICI BIZARRE 437 : * GBM MET A 0 CAR DISCONTINUITE DU TERME DEVANT DT, 0 AVANT 438 : * SAT ET POSITIF APRES. DOIT ETRE MIS DANS LES LOIS DE SAT 439 : * ET EVOLUER CONTINUEMENT. 440 : EMMAG = ('MASQUE' ('NOMC' SATC 'SCAL') 'EGSUP' 0.99 ) * COFEMMAG; 441 : 442 : * EMMAG = NOMC SCAL (SATC * COFEMMAG) ; 443 : EMMAG = 'KCHT' MDMH 'SCAL' 'CENTRE' 'COMP' 'SCAL' EMMAG ; 444 : * 'MESSAGE' 'maxi min emmag' ('MINIMUM' emmag) ('MAXIMUM' emmag); 445 : 446 : 447 : * Construction système matriciel 448 : * ============================== 449 : 450 : *-- Penalisation - 451 : * GBM ERREUR sur NOMC 'SOUR' 452 : 'SI' PENAL ; 453 : TERSC1 = TERSCE + ('NOMC' NOMESPL 454 : ((CHRG - CHRGANC) * VHYB / DELTAT * cofpenal)) ; 455 : 'SINON' ; 456 : TERSC1 = TERSCE ; 457 : 'FINSI' ; 458 : 459 : 460 : 'SI' ('EGA' CONSER VRAI) ; 461 : TERSC1 = TERSC1 - ('NOMC' NOMESPL 462 : ((TENC - TENNC) * VHYB / DELTAT )) ; 463 : 'FINSI' ; 464 : GEOL1 . 'SOURCE' = 'NOMC' NOMESPL TERSC1 ; 465 : 466 : * relaxation 467 : 468 : 469 : * 'MESSAGE' 'coefrelaaaaaaaaaaaaaaaa' cofrelax; 470 : 'SI' ('NEG' cofrelax 1.D0 1.D-14) ; 471 : * on relaxe CAPA et PERF 472 : CAPAR = (CAPA '*' cofrelax) 473 : '+' (CAPAN '*' (1. - cofrelax)) ; 474 : PERFR = ( PERFT '*' cofrelax) 475 : '+' (PERFN '*' (1. - cofrelax)) ; 476 : 'SINON' ; 477 : CAPAR = CAPA ; 478 : PERFR = PERFT ; 479 : 'FINSI' ; 480 : 481 : 482 : 483 : 484 : 485 : *-- coef devant la derivée temps 486 : COFDT = CAPAR '+' EMMAG ; 487 : 488 : * permeabilité 489 : * LAM = ('NOMC' PERFR 'K' NATURE DIFFUS) ; 490 : 491 : * GBM - penalisation un peu particuliere ! ecrire ce que cela 492 : * veut dire 493 : 'SI' PENAL ; 494 : COFDT = COFDT '+' COFPENAL ; 495 : 'FINSI' ; 496 : 497 : * initialisation TRANSGEOL - gbm bien verif les tabmodi 498 : GEOL1 . 'POROSITE' = COFDT ; 499 : TABMODI . 'POROSITE' = VRAI ; 500 : GEOL1 . 'DIFFUSIVITE' = PERFR ; 501 : TABMODI . 'DIFFUSIV' = VRAI ; 502 : 503 : * Résolution : 504 : * ============ 505 : 506 : * GBM appel trangeol - tester peut etre aussi premier penal 507 : 'SI' ( (&TRANSI 'EGA' 1)) ; 508 : GEOLPF1 GEOLPF2 = TRANGEOL MDMH GEOL1 ; 509 : 'SINON' ; 510 : GEOLPF1 GEOLPF2 = TRANGEOL MDMH GEOL1 GEOL2 ; 511 : 'FINSI' ; 512 : CHRG = GEOLPF1 . 'CONCENTRATION' ; 513 : FLU1 = GEOLPF1 . 'FLUXDIFF' ; 514 : 'SI' ('EGA' SATUR . 'TYPDISCRETISATION' EFMH) ; 515 : TRH2 = GEOLPF2 . 'TRACE_CONC' ; 516 : 'SINON' ; 517 : * GBM §§§§§§§§§§§§§§§§§§§§§§§ 518 : TRH2 = 0.D0 * FLU1; 519 : 'FINSI' ; 520 : 521 : * GBM - a clarifier 522 : TRH = TRH2 ; 523 : 524 : * FIN DE BOUCLE : 525 : * ================= 526 : 527 : 528 : 529 : * test de convergence sur le résidu 530 : * ================================= 531 : * A-t-on convergé à l'itération précédente (avec les nouveaux paramètres) ? 532 : * flux correspondant au résidu système résolu sans relaxation 533 : 534 : ** 535 : ** GBM reflechir comment améliorer la conservation 536 : 537 : 'SI' (&LINEAR 'EGA' 1) ; 538 : OKCONV = FAUX ; 539 : maxer = 1.D0; 540 : 'SINON' ; 541 : 542 : RES1 = 'MAXIMUM' ('RESULT' ('ABS' (TEN - TENN))) ; 543 : RES2 = 544 : ((0.1D0 * ('MAXIMUM' ('RESULT' ('ABS' (TEN))))) 545 : + ( ('MAXIMUM' ('RESULT' ('ABS' (TEN '-' TENANC)))))); 546 : 547 : 548 : 549 : 550 : RES3 = 'MAXIMUM' ('RESULT' ('ABS' (PERFR - PERFN))) ; 551 : RES4 = ('MAXIMUM' ('RESULT' ('ABS' ((perfr '-' perfn))))) 552 : '/' ('MAXIMUM' ('RESULT' PERFR)) ; 553 : RES3 = RES3 '/' ( 554 : 'MAXIMUM' ('RESULT' ('ABS' (PERFR - PERFANC))) + 1.D-100) ; 555 : 556 : 557 : 558 : 559 : * 'MESSAGE' 'minmax res1 res2' res1 res2 ('MAXIMUM' 560 : * (TEN - TENANC)) ; 561 : MAXER = RES1 '/' RES2 ; 562 : MAXER = MAXER '+' RES3 ; 563 : 564 : * on prend en compte la taille du pas de tps comparé à CFL 565 : 'SI' (&TRANSI 'EGA' 1) ; 566 : * Pour le premier pas de temps les flux ne sont pas connus 567 : * donc le pas de temps CFL est mal évalué. On ne l'inclue 568 : * pas dans le crtiere de convergence 569 : OKCONV = MAXER < (ERR0) ; 570 : 'SINON' ; 571 : OKCONV = MAXER < (ERR0 * DTI '/' DELTAT) ; 572 : * 'MESSAGE' 'error4 err3' res4 res3 ; 573 : OKCONV = OKCONV 'ET' (RES4 < ( ERR0)) ; 574 : 'FINSI' ; 575 : 576 : 'SI' (DTAUTO) ; 577 : 'SI' (&LINEAR 'EGA' 2) ; 578 : * on impose de boucler au moins 2 fois en pas de tps auto 579 : OKCONV = FAUX ; 580 : 'FINSI' ; 581 : 'FINSI' ; 582 : 'FINSI' ; 583 : 584 : 585 : * sorties textes 586 : * ============== 587 : * 588 : 589 : 'SI' debug ; 590 : * nb d'espaces avant le nb d'itérations, pour faire joli 591 : * moins joli si &LINEAR > 9999. 592 : IL = &LINEAR - 1 ; 593 : TXT = '| ' ; 594 : 'SI' (IL < 1000) ; txt = 'CHAINE' txt ' ' ; 'FINSI' ; 595 : 'SI' (IL < 100) ; txt = 'CHAINE' txt ' ' ; 'FINSI' ; 596 : 'SI' (IL < 10) ; txt = 'CHAINE' txt ' ' ; 'FINSI' ; 597 : * nb d'espaces avant le nb d'éléments désaturés, pour faire joli 598 : * moins joli si nbds > 9999. 599 : esp = ' ' ; 600 : 'SI' (nbds < 1000) ; esp = 'CHAINE' esp ' ' ; 'FINSI' ; 601 : 'SI' (nbds < 100) ; esp = 'CHAINE' esp ' ' ; 'FINSI' ; 602 : 'SI' (nbds < 10) ; esp = 'CHAINE' esp ' ' ; 'FINSI' ; 603 : 604 : * affichage des paramètres pour chaque itération 605 : * GBM change les maxi affichés - pas maxi TP1 mais maxi CHRG 606 : 'MESS' ('CHAINE' TXT IL ' |' MAXER ' | ' ('MINI' CHRG) 607 : ' | ' ('MAXI' CHRG) ' |' ('MAXI' CAPA) 608 : ' |' ('MINI' PERFR) ' | ' esp nbds ' | ') ; 609 : 'SI' OKCONV ; 610 : 'MESS' ('CHAINE' ' -------------------------------' 611 : '-------------------------------' 612 : '-------------------------------') ; 613 : 'FINSI' ; 614 : 'FINSI' ; 615 : 616 : * tests sortie 617 : * ============ 618 : 619 : 620 : *--- sortie si convergence 621 : 'SI' OKCONV ; 622 : * GBM - NON ON NE PEUT PAS SORTIR SI ON n'A PAS FAIT DE CALCUL 623 : * Le test précédent n'a pas de sens au premier indice. 624 : 'QUITTER' LINEAR ; 625 : 'FINSI' ; 626 : 627 : 628 : * petit message pour prévenir qu'on a changé la relaxation 629 : 'SI' (DEBUG 'ET' ('EGA' &LINEAR (itmaxi/2))) ; 630 : mess 'On accroît la sous-relaxation à 0.5' ; 631 : 'FINSI' ; 632 : 633 : 634 : *-- Sauvegarde des valeurs du pas de l'itération précédente 635 : TRH2N = TRH2 ; 636 : CAPAN = CAPAR ; 637 : PERFN = PERFR ; 638 : CHRGN = CHRG ; 639 : TENN = TEN ; 640 : 641 : 642 : *_____________________________________________________ non linéaire 643 : 644 : 'FIN' LINEAR ; 645 : *__________________________________________________________________ 646 : 647 : 'SI' ('EGA' CONSER VRAI) ; 648 : TENNC = TENC ; 649 : 'FINSI' ; 650 : 651 : * GBM - gerer les tabmodi à faux !!!!!!!!!!!!!!!!!!!!!!!!! 652 : 653 : 'SI' ('NON' OKCONV) ; 654 : * message, adaptation du pas de tps ou penalisation 655 : DELTAT CFL0 = SATUTILS TESARRET &PENALDT NPENALDT LAST1 &LINEAR 656 : DELTAT maxer OKPENAL 657 : cofpenal CFL0 COFDIV debug 658 : SATUR tps ; 659 : 'SINON' ; 660 : * on a convergé 661 : * PENAL = FAUX ; 662 : 'QUITTER' PENALDT ; 663 : 'FINSI' ; 664 : 665 : 666 : *--------------------------------------------------------- artifice 667 : 668 : 'FIN' PENALDT ; 669 : 670 : *------------------------------------------------------------------ 671 : 672 : 673 : FLU0 = FLU1 ; 674 : LAST1 = LAST1 + 1 ; 675 : 676 : * Sauvegarde 677 : * ========== 678 : 679 : ISOR ISAUV = SATUTILS SAVRESU SATUR TPSOR ISOR TPSANC TPS ISAUV DELTAT 680 : TRHANC TRH CHRGANC CHRG TENANC TEN PNSCANC PNSC SATANC SAT FLUANC 681 : FLU0 NIVEAU DSOR &LINEAR MDMH ; 682 : 683 : 684 : 'MESSAGE' ' ' ; 685 : 'MESSAGE' ' ' ; 686 : 'MESSAGE' 'RESIDU EN TEMPS ' 687 : ((tfinal '-' tpsini) * 688 : ('MAXIMUM' (('RESULT' ('ABS' (TEN '-' TENANC))) '/' 689 : (deltat * ('RESULT' ('ABS' TEN)))))) ; 690 : 'MESSAGE' ' ' ; 691 : 'MESSAGE' ' ' ; 692 : 693 : 694 : 695 : * Temps caractéristique et nouveau pas de temps : 696 : * =============================================== 697 : 698 : *-- pas de temps idéal lié à la cinétique d'hydratation : 699 : * nb de face avec des flux pas trop faibles (> 1.D-9 x maximum) 700 : 701 : 702 : 703 : *-- Nouveau pas de temps : 704 : 'SI' ('EXISTE' SATUR 'TEMPS_CALCULES') ; 705 : 'SI' (ICAL < DCAL) ; 706 : ICAL = ICAL + 1 ; 707 : DELTAT = ('EXTR' LTCALCUL ICAL) - TPS ; 708 : * GBM verifier que DTAUTO = FAUX dans ce cas sinon écrasé apres. 709 : 'FINSI' ; 710 : 'SINON' ; 711 : * quand on a donné SATUR.'DT_INITIAL', le pas de temps n'est 712 : * imposé qu'au départ. Après, il est automatique. 713 : DTAUTO = VRAI ; 714 : 'FINSI' ; 715 : 'SI' DTAUTO ; 716 : * pas de temps automatique 717 : * on modifie le CFL de façon à approcher nb d'itérations visé. 718 : 'MESSAGE' 'avant' CFL0; 719 : CFL0 = ((((Niter0 '/' 1.)) / &LINEAR) ** 0.5) ; 720 : 721 : 'MESSAGE' 'apres' CFL0; 722 : SATUR.'CFL' = CFL0 ; 723 : 'FINSI' ; 724 : 725 : VVOL = VHYB * PORO1 ; 726 : 'MESS' 'mini max capa' ('MINIMUM' (CAPA)) 727 : ('MAXIMUM' (CAPA)) ; 728 : *'MESS' 'mini max perm' ('MINIMUM' (PERFR)) 729 : * ('MAXIMUM' (PERFR)) ; 730 : DELTAT DTI = SATUTILS YNYTDT MDMH FLU0 MCHYB VVOL CFL0 debug 731 : DTAUTO DELTAT TERSCE ; 732 : 733 : *-- sortie en cas de temps limite dépassé 734 : 'SI' ( TPS '>EG' TFINAL ) ; 735 : 'MENAGE' ; 736 : 'QUITTER' TRANSI ; 737 : 'FINSI' ; 738 : 739 : * Préparation pas de temps suivant : 740 : * ================================== 741 : 742 : * Sauvegarde des valeurs du pas de temps précédent 743 : TRHANC = TRH ; 744 : CHRGANC = CHRG ; 745 : FLUANC = FLU0 ; 746 : PNSANC = PNS ; 747 : PNSCANC = PNSC ; 748 : TENANC = TEN ; 749 : SATANC = SAT ; 750 : TPSANC = TPS ; 751 : PERFANC = PERFR ; 752 : * GBM REGARDE 753 : 'SI' ('EXISTE' SATUR 'FLUX_IMPOSE') ; 754 : FLUIMPM1 = FLUIMPO ; 755 : 'FINSI' ; 756 : 'SI' ('EXISTE' SATUR 'FLUMIXTE') ; 757 : FLUMMPM1 = FLUMMPO ; 758 : 'FINSI' ; 759 : 'SI' ('EXISTE' SATUR 'SOURCE') ; 760 : TERSC2M1 = TERSC2 ; 761 : 'FINSI' ; 762 : 763 : GEOL1 = 'TABLE' GEOLPF1 ; 764 : GEOL2 = 'TABLE' GEOLPF2 ; 765 : 766 : 767 : 768 : *================================================= boucle transitoire 769 : 770 : 'FIN' TRANSI ; 771 : 772 : *==================================================================== 773 : * 774 : 'FINPROC' ; 775 : 776 : 777 :
© Cast3M 2003 - All rights reserved.
Disclaimer