1 : * UPDAVF PROCEDUR GOUNAND 11/05/24 21:16:24 6976 2 : ********************************************************************** 3 : 'DEBP' UPDAVF MoDARCY*'MMODEL' Porosite*'CHPOINT' Matediff*'CHPOINT' 4 : Diffdisp*'CHPOINT' ChPSour*'CHPOINT' Deltat*'FLOTTANT' 5 : Cini*'CHPOINT' TetaDiff*'FLOTTANT' TetaConv*'FLOTTANT' 6 : Qface/'CHPOINT' nomespec*'LISTMOTS' nbespece*'ENTIER' 7 : nbsource*'ENTIER' mattr*'MATRIK' JACO*'MATRIK' 8 : MPOR2*'MATRIK' MCHAMT*'MCHAML' MCHAMT1*'MCHAML' 9 : TABMODI*'TABLE' CHCLIM*'TABLE' LCONV*'LOGIQUE'; 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 Mattt Difftot Mctot Mdiff Nouvmat = UPDAVF 22 : * MoDARCY Porosite Matediff Diffdisp ChPSour 23 : * DeltaT Cini TetaDiff TetaConv 24 : * Qface nomespec nbespece nbsource 25 : * Matot Jaco Mctot Mdiff Mpor TABMODI CHCLIM ; 26 : * |-----------------------------------------------------------------| 27 : * | Généralités : MATTVF 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 : * | points de composante 'K' en isotrope, 'K11', 'K21',| 44 : * | 'K22' en anisotrope 2d et 'K11', 'K21', 'K22', 'K31'| 45 : * | 'K32', 'K33' en anisotrope 3d. Type 'CARACTERISTIQUE'| 46 : * | | 47 : * | Diffdisp : Tenseur de dispersion (type iso, ..) champ par | 48 : * | points de composante 'K' en isotrope, 'K11', 'K21',| 49 : * | 'K22' en anisotrope 2d et 'K11', 'K21', 'K22', 'K31'| 50 : * | 'K32', 'K33' en anisotrope 3d. Type 'CARACTERISTIQUE'| 51 : * | | 52 : * | ChPSour : Champ par points des sources volumiques par unité de | 53 : * | temps (support maillage centre). Composante ?????? | 54 : * | | 55 : * | Cini : Concentration initiale, CHPOINT centre. | 56 : * | Composante 'H'. | 57 : * | | 58 : * | Deltat : Pas de temps | 59 : * | | 60 : * | Qface : vitesse aux faces, CHPO face de composantes Vx, Vy | 61 : * | en 2d et Vx, Vy, Vz en 3d. Il s'agit plus exatement | 62 : * | de (V.n)n, c'est à dire de la composante normale de | 63 : * | la vitesse aux faces. ???????? (je pressens que | 64 : * | castem va sortir des flux, cad intégrés sur surfaces)| 65 : * | | 66 : * | nomespec : liste des noms de composante des espèces dans Cini | 67 : * | | 68 : * | nbespece : nombre de composante de Cini, soit nombre d'especes | 69 : * | | 70 : * | nbsource : nombre de composantes du terme source qd X especes | 71 : * | | 72 : * | Matot : matrice globale de discretisation en VF | 73 : * | | 74 : * | Jaco : matrice globale de discretisation en VF pour le problème 75 : * | stationnaire 76 : * 77 : * | Mpor : matrice globale de discretisation en VF pour le problème 78 : * | stationnaire 79 : * | | 80 : * | Mchamt : Coef permettant de calculer le flux total 81 : * | | 82 : * | Mchamt1 : Coef permettant de calculer le flux diffusif 83 : * | | 84 : * | | 85 : * | TABMODI : table contenant des logiques indiquant la nécessité | 86 : * | ou non de reclalculer certains termes. | 87 : * | 'POROSITE' : VRAI si le coefficient devant D/DT | 88 : * | (porosité) est modifié depuis le dernier| 89 : * | appel | 90 : * | 'DELTAT' : VRAI si le pas de tps a changé | 91 : * | 'CONVECTI' : VRAI si la vitesse a changé | 92 : * | 'COEF_LIN' : VRAI si le coef en facteur de C a changé| 93 : * | 'DIFFUSI' : VRAI si les diffusivités ont changé | 94 : * | | 95 : * | CHCLIM : table d'indice 'NEUMANN' et 'DIRICHLET' contenant les| 96 : * | Chpoint à n composantes contenant les conditions aux | 97 : * | limites de Neumann et Dirichlet par espece. | 98 : * | | 99 : * | | 100 : * |-----------------------------------------------------------------| 101 : * | ENTREES-SORTIES | 102 : * |-----------------------------------------------------------------| 103 : * | | 104 : * | Difftot : Coefficient de diffusion totale, integre decentrement| 105 : * | | 106 : * | | 107 : * |-----------------------------------------------------------------| 108 : * | SORTIES | 109 : * |-----------------------------------------------------------------| 110 : * | | 111 : * | | 112 : * | RESI : second membre | 113 : * | | 114 : * | Matot : matrice globale de discretisation en VF | 115 : * | | 116 : * | Difftot : Coefficient de diffusion totale, integre decentrement| 117 : * 118 : * | Mpor : matrice globale de discretisation en VF pour le problème 119 : * | stationnaire 120 : * | | 121 : * | Mchamt : Coef permettant de calculer le flux total 122 : * | | 123 : * | Mchamt1 : Coef permettant de calculer le flux diffusif 124 : * | | 125 : * | | 126 : * |-----------------------------------------------------------------| 127 : * | VARIABLES INTERNES | 128 : * |-----------------------------------------------------------------| 129 : * | | 130 : * | | 131 : * | PCONV : Logique indiquant VRAI si présence de convection | 132 : * | | 133 : * | CoefDt : coeff devant dC/dt integre sur les elements | 134 : * | | 135 : * | CCini : concentration aux centres (une composante) | 136 : * | | 137 : * | SSource : Source aux centre (une composante) | 138 : * | | 139 : * | SSMTr : second membre sur les traces pour une espèce | 140 : * | | 141 : * | toltheta : 1.D-4 seuil en dessous duquel on considère que la | 142 : * | valeur de theta du theta-schéma est nulle (schéma | 143 : * | explicite) ou au contraire euler-implicite si | 144 : * | theta > 0.9999 | 145 : * | | 146 : * | tetaconv : valeur du coef d'implicitation pour la convection | 147 : * | | 148 : * | MODDIFF : Logique, si VRAI on recalcule les termes liés à la | 149 : * | diffusion, en supposant qu'elle a changé, donnés par | 150 : * | Matediff. En implicite cela affecte la matrice | 151 : * | rigidité du problème. | 152 : * | | 153 : * | MODPOROS : Si VRAI on recalcule la contribution de la porosité | 154 : * | pour la matrice du problème et le coeficient devant | 155 : * | dC/dt (VF) intégré sur les mailles. la matrice | 156 : * | Identité des VF est recalculée. Logique. | 157 : * | | 158 : * | MODELTT : Logique. Si vrai on recalcule la matrice rigidité du | 159 : * | problème. Le pas de temps a changé. | 160 : * | | 161 : * | MODCONV : Logique. Si vrai on recalcule la matrice rigidité du | 162 : * | problème. La vitesse a changé. | 163 : * | | 164 : * | MODCOLIN : Logique. Si vrai on recalcule la matrice rigidité du | 165 : * | problème. Le coef linéaire devant C a changé. | 166 : * | | 167 : * | FLUNEU : LOGIQUE valant VRAI si conditions de Neumann | 168 : * | | 169 : * | CLFLUX : Chpoint à n composantes contenant les flux imposés | 170 : * | pour chaque espece chimique. nul si pas de flux | 171 : * | OPTIONNEL | 172 : * | | 173 : ********************************************************************** 174 : 175 : 176 : *--------------------------------------------------------------------- 177 : *---------- On récupere les conditions limites ------------------ 178 : *--------------------------------------------------------------------- 179 : 180 : FLUNEU = FAUX; 181 : DIRCLI = FAUX; 182 : FLUTOT = FAUX; 183 : FLUMIX = FAUX; 184 : CLFLUX MA = 'KOPS' 'MATRIK'; 185 : CLDIRI MA = 'KOPS' 'MATRIK'; 186 : CLFLUT MA = 'KOPS' 'MATRIK'; 187 : EXFLU = CLFLUX; 188 : EXDIR = CLDIRI; 189 : EXFLUT = CLFLUT; 190 : * Neumann 191 : 'SI' ('EXISTE' CHCLIM 'NEUMANN') ; 192 : CLFLUX = CHCLIM . 'NEUMANN'; 193 : FLUNEU = VRAI; 194 : 'FINSI'; 195 : 196 : 'SI' ('EXISTE' CHCLIM 'DIRICHLET') ; 197 : CLDIRI = CHCLIM . 'DIRICHLET'; 198 : DIRCLI = VRAI; 199 : 'FINSI'; 200 : 201 : 'SI' ('EXISTE' CHCLIM 'FLUTOTAL') ; 202 : CLFLT = CHCLIM . 'FLUTOTAL'; 203 : FLUTOT = VRAI; 204 : 'FINSI'; 205 : 206 : * Flux mixte 207 : 'SI' ('EXISTE' CHCLIM 'FLUMIXTE') ; 208 : * comme on impose A Dgrad C + B C = flumix, on le traite sous 209 : * la forme D grad C + (B/A) C = flumix/A plus naturelle en EFMH car 210 : * D grad C est le flux diffusif 211 : COFA = -1.D0 * CHCLIM . 'FLUMIXTE' . 'COEFA' ; 212 : CLFLUX3 = CHCLIM . 'FLUMIXTE' . 'VAL' '/' COFA ; 213 : CLFLUX3 = CHAN 'ATTRIBUT' CLFLUX3 NATURE DISCRET ; 214 : FLUMIX = VRAI ; 215 : mayage = 'EXTRAIRE' CHCLIM . 'FLUMIXTE' . 'VAL' maillage ; 216 : cofb = (doma modarcy SURFACE) *CHCLIM . 'FLUMIXTE' . 'COEFB' ; 217 : cofb = CHCLIM . 'FLUMIXTE' . 'COEFB' ; 218 : cofa = 'REDU' (CHCLIM . 'FLUMIXTE' . 'COEFA') mayage ; 219 : cofb = 'REDU' cofb mayage ; 220 : coefm = (1.D0 * cofb) '/' cofa ; 221 : CLFLUX3 = 'REDU' CLFLUX3 mayage ; 222 : 'OUBLIER' cofa ; 223 : 'OUBLIER' cofb ; 224 : 'FINSI' ; 225 : 226 : 227 : *--------------------------------------------------------------------- 228 : *---------- Initialisations de tables, coefficients ------------------ 229 : *--------------------------------------------------------------------- 230 : 231 : toltheta = 1.D-4 ; 232 : 233 : MODPOROS = TABMODI . 'POROSITE' ; 234 : MODCONV = TABMODI . 'CONVECTI' ; 235 : MODDELTT = TABMODI . 'DELTAT' ; 236 : MODCOLIN = TABMODI . 'COEF_LIN' ; 237 : MODDIFF = TABMODI . 'DIFFUSIV' ; 238 : 239 : 240 : 241 : 242 : * Si la porosité est modifiée (ou plus exactement le coef 243 : * devant dC/dt) 244 : 245 : 246 : 'SI' (MODPOROS) ; 247 : MPOR2 = 'KOPS' 'MATDIAGO' ('NOMC' 'RETN' Porosite) 'MATRIK' ; 248 : 'FINSI' ; 249 : 250 : * 251 : * Si le pas de temps est modifié 252 : * 253 : 254 : 255 : * 256 : * Si le coefficient de décroissance (ou linaire) en facteur de C est 257 : * modifié ????????????????, 258 : * 259 : 260 : 'SI' (MODCOLIN) ; 261 : * A priori sera géré dans le coeff de porosité 262 : 'FINSI' ; 263 : 264 : 265 : 266 : 267 : 268 : *---------------------------------------------------------------------- 269 : *---------- CALCUL DE LA MATRICE DE DIFFUSION DU PROBL ---------------- 270 : *---------------------------------------------------------------------- 271 : 272 : 273 : * Calcul de la matrice du probleme diffusion transitoire 274 : * Elle est modifiée si la diffusion , 275 : * ou si le pas de temps change, ou si le coefficient 276 : * en facteur de la dérivée temporelle change (porosité*retard), 277 : * ou si le coefficient du terme linéaire a(x,t) * C(x,t) change. 278 : 279 : * Pour des raisons pratiques on n'a pas stocké les matrices 280 : * de convection et de diffusion séparément. Si l'une est modifiée 281 : * on doit donc recalculer les deux matrices. Ce choix est discutable, 282 : * on gagne un facteur 3 en stockage (mat totale au lieu de mat totale 283 : * + diffusion + convection). Par ailleurs, le cout de l'assemblage + 284 : * resolution + créer une des deux matrice rend souvent négligeable 285 : * le recalcul inutile d'une des deux matrices précitées. On réduit 286 : * également le nombre d'entrées-sorties. 287 : 288 : 289 : 'SI' (MODDIFF) ; 290 : * La diffusivité effective a changé, on remmet la matrice au propre 291 : * c'est à dire avec des noms de composantes isotropes K11 K22 etc ... 292 : zozo = DIFFANIS Matediff ; 293 : Matediff = zozo ; 294 : 'FINSI' ; 295 : 296 : Difftot = MateDiff '+' Diffdisp ; 297 : 298 : 299 : *---------------------------------------------------------------------- 300 : *--- CALCUL DE LA MATRICE DE CONVECTION ET DES SECONDS MEMBRES -------- 301 : *---------------------------------------------------------------------- 302 : 303 : 304 : * Nombres d'espèces à gérer est fixé et ne peut changer au cours d'une 305 : * simulation. Il sera d'ailleurs mis dans le modele 306 : 307 : 'SI' (nbespece > 0) ; 308 : 309 : 310 : DTI = 1.D0/DeltaT ; 311 : 'REPETER' bloc1 nbespece ; 312 : 313 : espc = ('EXTRAIRE' &bloc1 nomespec) ; 314 : 315 : * On extrait la composante de Cini, Tcini et de la source 316 : CCini = 'NOMC' 'H' ('EXCO' espc Cini) ; 317 : 'SI' (nbsource > 1) ; 318 : SSource = 'NOMC' 'SOUR' ('EXCO' espc 319 : Chpsour) ; 320 : 'SINON' ; 321 : SSource = 'NOMC' 'SOUR' Chpsour ; 322 : 'FINSI' ; 323 : * Conditions initiales 324 : * Prise en compte du terme source et eventuellement 325 : * de la convection avec le schema centre 326 : HHS = 'NOMC' 'SCAL' CCini ; 327 : * Prise en compte des conditions aux limites 328 : 'SI' (FLUNEU) ; 329 : EXFLU = 330 : ('NOMC' 'FLUX' ('EXCO' espc 331 : CLFLUX)) ; 332 : * MESS 'EXFLU ='; LIST EXFLU; 333 : 'FINSI' ; 334 : 'SI' (DIRCLI) ; 335 : EXDIR = 336 : ('NOMC' 'SCAL' ('EXCO' espc 337 : CLDIRI)) ; 338 : * MESS 'EXDIR ='; LIST EXDIR; 339 : 'FINSI' ; 340 : 'SI' (FLUTOT) ; 341 : EXFLUT = 342 : ('NOMC' 'FLUX' ('EXCO' espc 343 : CLFLT )) ; 344 : QLIM = NOMC(EXFLUT) 'FLUX' ; 345 : SUP = 'EXTR' EXFLUT 'MAIL' ; 346 : XPAR1 = 1.D0 + (0.0*CLFLT) ; 347 : XPAR1 = 'NOMC' XPAR1 'PAR1' ; 348 : 'SI' LCONV ; 349 : USCNR = 'REDU' QFACE SUP ; 350 : MUSCN = USCNR*(-1.D0) ; 351 : 'SINON' ; 352 : MUSCN = (0.D0*CLFLT) ; 353 : 'FINSI' ; 354 : 355 : XPAR2 = 'NOMC' MUSCN 'PAR2' ; 356 : EXFLUT = XPAR1 + XPAR2 + QLIM ; 357 : EXFLUT = CHAN 'ATTRIBUT' EXFLUT NATURE DISCRET ; 358 : * MESS 'EXFLUT ='; LIST EXFLUT; 359 : 'FINSI' ; 360 : 361 : 'SI' (FLUMIX) ; 362 : EXFLUM = 363 : ('NOMC' 'FLUX' ('EXCO' espc 364 : CLFLUX3 )) ; 365 : QLIM = NOMC(EXFLUM) 'FLUX' ; 366 : SUP = 'EXTR' EXFLUM 'MAIL' ; 367 : XPAR1 = 1.D0 + (0.0*CLFLUX3) ; 368 : * MESS 'XPAR1 ='; 369 : XPAR1 = 'NOMC' XPAR1 'PAR1' ; 370 : MUSCN = coefm; 371 : XPAR2 = 'NOMC' MUSCN 'PAR2' ; 372 : EXFLUM = XPAR1 + XPAR2 + QLIM ; 373 : EXFLUM = CHAN 'ATTRIBUT' EXFLUM NATURE DISCRET ; 374 : EXFLUT = EXFLUT 'ET' EXFLUM ; 375 : 'FINSI' ; 376 : * MESS 'UPDAVF EXFLUT ='; LIST EXFLUT; 377 : * MESS 'UPDAVF EXFLU ='; LIST EXFLU; 378 : 379 : * 'MESS' 'MODDIFF'; LIST MODDIFF; 380 : * 'MESS' 'MODCONV'; LIST MODCONV; 381 : 'SI' (MODDIFF 'OU' MODCONV) ; 382 : * On recalcule deux fois les coefficients. En général, cela doit etre 383 : * négligeable en temps de calcul; 384 : GRADT0 MCHAMT1 = 'PENT' MoDARCY 'FACE' 'MPFA' 385 : HHS 'DISPDIF' difftot 'TIMP' EXDIR 'QIMP' EXFLU 386 : 'MIXT' EXFLUT ; 387 : MCHAMT = MCHAMT1; 388 : MCHAJA = TetaDiff*MCHAMT1; 389 : 390 : 'SI' LCONV ; 391 : GRADT0 MCHAMT = 'PENT' MoDARCY 'FACE' 'MPFA' 392 : HHS 'DISPDIF' difftot 'TIMP' EXDIR 'QIMP' EXFLU 393 : 'MIXT' EXFLUT 'UPWICENT' QFace ; 394 : * ATTENTION PEUT ETRE FAUX SI TetaConv different de TetaDiff 395 : MCHCONV = MCHAMT + ((-1.D0)*MCHAMT1) ; 396 : MCHAJA = (TetaConv*MCHCONV) + (TetaDiff*MCHAMT1) ; 397 : 'DETRUIT' MCHCONV ; 398 : 'FINSI' ; 399 : 400 : 401 : * On ne calcule la jacobienne que pour la première espèce 402 : 'SI' ((&bloc1) 'EGA' 1) ; 403 : JACO CHPRES DT = 'LAPN' 'VF' 'CLAUDEIS' 'IMPL' 404 : MoDARCY HHS GRADT0 MCHAJA 'QIMP' EXFLU 'MIXT' EXFLUT 405 : 'TIMP' EXDIR ; 406 : 'SINON' ; 407 : JACOB CHPRES DT = 'LAPN' 'VF' 'CLAUDEIS' 'EXPL' 408 : MoDARCY HHS GRADT0 ; 409 : 'FINSI' ; 410 : * SINON; 411 : * 'DETRUIT' MCHAJA ; 412 : 'SI' ( ((&bloc1) 'EGA' 1) ) ; 413 : MATOT = 'KOPS' 'MULT' -1.000D0 JACO ; 414 : MATI = '*' MPOR2 DTI ; 415 : MATOT = MATOT 'ET' MATI ; 416 : 'FINSI' ; 417 : 'SINON'; 418 : 'SI' ( ((&bloc1) 'EGA' 1) 'ET' (MODPOROS 'OU' MODDELTT)) ; 419 : MATOT = 'KOPS' 'MULT' -1.000D0 JACO ; 420 : MATI = '*' MPOR2 DTI ; 421 : MATOT = MATOT 'ET' MATI ; 422 : 'FINSI'; 423 : * 'SINON'; 424 : * On reconstitue un champ de second membre 425 : GRADT0 = 'PENT' MoDARCY 'FACE' 'MPFA' 426 : HHS 'DISPDIF' difftot 427 : 'TIMP' EXDIR 'QIMP' EXFLU 'MIXT' EXFLUT 'UPWICENT' QFace 428 : 'GRADGEO' MCHAMT ; 429 : JACO CHPRES DT = 'LAPN' 'VF' 'CLAUDEIS' 'EXPL' 430 : MoDARCY HHS GRADT0 ; 431 : 432 : 'FINSI' ; 433 : 434 : 'SI' ((&bloc1) 'EGA' 1) ; 435 : SSMTr = 'NOMC' espc CHPRES ; 436 : RESI = SSMTr ; 437 : 'SINON' ; 438 : CHPRES = 'NOMC' espc 439 : ('COPIER' CHPRES) ; 440 : RESI = CHPRES 'ET' RESI; 441 : 'FINSI' ; 442 : 'FIN' bloc1 ; 443 : 444 : 'FINSI' ; 445 : 446 : 447 : *--------------------------------------------------------------------- 448 : *------ Conditions aux limites de Flux (mixtes et Neumann) ----------- 449 : *--------------------------------------------------------------------- 450 : 451 : 452 : *--------------------------------------------------------------------- 453 : *------ Matrice sortante --------------------------------- ----------- 454 : *--------------------------------------------------------------------- 455 : 456 : 457 : * si pas de modif et matrice matrik en entrée 458 : NOUVMAT = VRAI; 459 : 'SI' (('EXISTE' mattr) 'ET' ('NON' (MODDIFF 'OU' MODCONV 'OU' 460 : MODPOROS 'OU' MODDELTT))) ; 461 : Matot = mattr ; 462 : NOUVMAT = FAUX ; 463 : 'FINSI' ; 464 : 465 : 466 : 'FINP' RESI Matot Difftot MCHAMT MCHAMT1 NOUVMAT ; 467 : 468 : 469 : 470 : 471 : 472 : 473 : 474 : 475 :
© Cast3M 2003 - All rights reserved.
Disclaimer