$$$$ UPDAEFMH * UPDAEFMH PROCEDUR CHAT 12/08/07 21:15:14 7481 ********************************************************************** 'DEBP' UPDAEFMH MoDARCY*'MMODEL' Porosite*'CHPOINT' MateDiff*'CHPOINT' Diffdisp*'CHPOINT' ChPSour*'CHPOINT' Cini*'CHPOINT' Tcini*'CHPOINT' Deltat*'FLOTTANT' Qface/'CHPOINT' nomespec*'LISTMOTS' nbespece*'ENTIER' nbsource*'ENTIER' LMLump*'LOGIQUE' DECENTR*'LOGIQUE' massEFMH*'RIGIDITE' mattr/'MATRIK' mattm/'RIGIDITE' Tbdartra*'TABLE' TABMODI*'TABLE' CHCLIM*'TABLE'; * ATTENTION La vitesse est optionnelle, L'ordre est important * et les types d'arguments qui se suivent aussi pour tester leur * présence * * Attention il faudra transformer les vitesses en débits face * et sortir le champ * * |-----------------------------------------------------------------| * | Phrase d'appel (en GIBIANE) | * |-----------------------------------------------------------------| * | | * |SMTr MatrTr TbDarTra MassEFMH Difftot = UPDAEFMH MoDARCY Porosite| * | MateDiff difftot ChPSour cini tcini deltat | * | (Qface) nomespec nbespece nbsource LMLump | * | DECENTR massEFMH | mattr tbdartra TABMODI; | * | | mattm | * |-----------------------------------------------------------------| * | Généralités : MATTEFMH construit la matrice de discrétisation | * | du problème de transport convection-diffusion pour| * | le premier pas de tps d'un algorithme transitoire.| * | Le second membre et les Conditions limites de flux| * | sont pris en compte. | * | RESTE TCINI, DECENTR et TERME LIN | * |-----------------------------------------------------------------| * | | * |-----------------------------------------------------------------| * | ENTREES | * |-----------------------------------------------------------------| * | MoDARCY : modele Darcy. | * | | * | Porosite : champ par elements de composante 'CK' | * | | * | MateDiff : Tenseur de diffusion (type iso, ..) champ par | * | points de composante 'K' en isotrope, 'K11', 'K21', | * | 'K22' en anisotrope 2d et 'K11', 'K21', 'K22', 'K31'| * | 'K32', 'K33' en anisotrope 3d. Type 'CARACTERISTIQUE'| * | | * | Diffdisp : Tenseur de dispersion (type iso, ..) champ par | * | points de composante 'K' en isotrope, 'K11', 'K21', | * | 'K22' en anisotrope 2d et 'K11', 'K21', 'K22', 'K31'| * | 'K32', 'K33' en anisotrope 3d. Type 'CARACTERISTIQUE'| * | | * | ChPSour : Champ par points des sources volumiques par unité de | * | temps (support maillage centre). Composante ?????? | * | | * | Cini : Concentration initiale, CHPOINT centre. | * | Composante 'H'. | * | | * | Tcini : Trace de concentration aux faces (eventuellement à | * | plusieurs composantes (espèces) | * | | * | Deltat : Pas de temps | * | | * | Qface : vitesse aux faces, CHPO face de composantes Vx, Vy | * | en 2d et Vx, Vy, Vz en 3d. Il s'agit plus exatement | * | de (V.n)n, c'est à dire de la composante normale de | * | la vitesse aux faces. ???????? (je pressens que | * | castem va sortir des flux, cad intégrés sur surfaces)| * | | * | nomespec : liste des noms de composante des espèces dans Cini | * | | * | nbespece : nombre de composante de Cini, soit nombre d'especes | * | | * | nbsource : nombre de composantes du terme source qd X especes | * | | * | LMLump : Logique. Si vrai on effectue une condensation de | * | masse de la matrice EFMH | * | | * | DECENTR : Logique. Vrai veut dire schémas décentrés et faux | * | veut dire schéma convectif centré. | * | | * | MatTm : matrice globale sur les traces. MATRIK en entrée | * | sort MATRIK si non modifiée, RIGIDITE sinon | * | Soit on rentre cet argument soit le suivant Mattr | * | | * | MatTr : idem mais rigidité en entrée on ressort cette matrice| * | inchangée si les options MATMODI indiquent aucune | * | modif. Optionnel. On rentre Mattm si absent. | * | | * | TbDarTra : table Darcy transitoire utilisée par MHYB, SMTP ... | * | | * | TABMODI : table contenant des logiques indiquant la nécessité | * | ou non de reclalculer certains termes. | * | 'POROSITE' : VRAI si le coefficient devant D/DT | * | (porosité) est modifié depuis le dernier| * | appel | * | 'DELTAT' : VRAI si le pas de tps a changé | * | 'CONVECTI' : VRAI si la vitesse a changé | * | 'COEF_LIN' : VRAI si le coef en facteur de C a changé| * | 'DIFFUSI' : VRAI si les diffusivités ont changé | * | | * | CHCLIM : table d'indice 'NEUMANN' et 'DIRICHLET' contenant les| * | Chpoint à n composantes contenant les conditions aux | * | limites de Neumann et Dirichlet par espece. | * | | * | | * |-----------------------------------------------------------------| * | ENTREES-SORTIES | * |-----------------------------------------------------------------| * | | * | MassEFMH : matrice elementaire EFMH | * | | * | Remarque | * | -------- | * | On a toujours intéret à rentrer Mattm si on l'a et qu'il n'y a | * | pas de modification, afin de conserver les factorisations LU | * | dans solvEFMH en cas de résolution par méthode directe. | * | | * |-----------------------------------------------------------------| * | SORTIES | * |-----------------------------------------------------------------| * | | * | SMTr : second membre sur les traces | * | | * | Matrtr : Matrice du problème, rigidité si nouvelle, égale à la| * | matrice entrée (mattm ou mattr) sinon | * | | * |-----------------------------------------------------------------| * | VARIABLES INTERNES | * |-----------------------------------------------------------------| * | | * | LCALMATP : Logique, VRAI si on recalcule la matrice du systeme | * | avec diffusion produite par MATP | * | | * | LCALCONV : Logique, VRAI si on recalcule la matrice du systeme | * | avec convection produite par SMTP | * | | * | PCONV : Logique indiquant VRAI si présence de convection | * | | * | CoefDt : coeff devant dC/dt integre sur les elements | * | | * | Numdiff : diffusivité numérique due au décentrement | * | | * | CCini : concentration aux centres (une composante) | * | | * | Tccini : Trace de concentration aux faces (une composante) | * | | * | SSource : Source aux centre (une composante) | * | | * | SSMTr : second membre sur les traces pour une espèce | * | | * | MatConv : matrice globale sur les traces pour la convection | * | | * | Matrdiff : matrice globale sur les traces pour la diffusion | * | | * | toltheta : 1.D-4 seuil en dessous duquel on considère que la | * | valeur de theta du theta-schéma est nulle (schéma | * | explicite) ou au contraire euler-implicite si | * | theta > 0.9999 | * | | * | tetaconv : valeur du coef d'implicitation pour la convection | * | | * | MODDIFF : Logique, si VRAI on recalcule les termes liés à la | * | diffusion, en supposant qu'elle a changé, donnés par | * | Matediff. En implicite cela affecte la matrice | * | rigidité du problème. | * | | * | MODPOROS : Si VRAI on recalcule la contribution de la porosité | * | pour la matrice du problème et le coeficient devant | * | dC/dt (EFMH) intégré sur les mailles. la matrice | * | Identité des VF est recalculée. Logique. | * | | * | MODELTT : Logique. Si vrai on recalcule la matrice rigidité du | * | problème. Le pas de temps a changé. | * | | * | MODCONV : Logique. Si vrai on recalcule la matrice rigidité du | * | problème. La vitesse a changé. | * | | * | MODCOLIN : Logique. Si vrai on recalcule la matrice rigidité du | * | problème. Le coef linéaire devant C a changé. | * | | * | FLUNEU : LOGIQUE valant VRAI si conditions de Neumann | * | | * | CLFLUX : Chpoint à n composantes contenant les flux imposés | * | pour chaque espece chimique. nul si pas de flux | * | OPTIONNEL | * | | ********************************************************************** *--------------------------------------------------------------------- *---------- On récupere les conditions limites ------------------ *--------------------------------------------------------------------- FLUNEU = FAUX ; FLUTOT = FAUX ; FLUMIX = FAUX ; FLUCLI = FAUX ; * Neumann 'SI' ('EXISTE' CHCLIM 'NEUMANN') ; CLFLUX1 = CHCLIM . 'NEUMANN' ; FLUNEU = VRAI ; 'SINON' ; * on crée un champ vide CLFLUX1 dum = 'KOPS' MATRIK ; 'OUBLIER' dum ; 'FINSI' ; * Flux total 'SI' ('EXISTE' CHCLIM 'FLUTOTAL') ; CLFLUX2 = CHCLIM . 'FLUTOTAL' ; FLUTOT = VRAI ; 'SINON' ; * on crée un champ vide CLFLUX2 dum = 'KOPS' MATRIK ; 'OUBLIER' dum ; 'FINSI' ; * Flux mixte 'SI' ('EXISTE' CHCLIM 'FLUMIXTE') ; * comme on impose A Dgrad C + B C = flumix, on le traite sous * la forme D grad C + (B/A) C = flumix/A plus naturelle en EFMH car * D grad C est le flux diffusif COFA = -1.D0 * CHCLIM . 'FLUMIXTE' . 'COEFA' ; CLFLUX3 = CHCLIM . 'FLUMIXTE' . 'VAL' '/' COFA ; CLFLUX3 = CHAN 'ATTRIBUT' CLFLUX3 NATURE DISCRET ; FLUMIX = VRAI ; 'SINON' ; * on crée un champ vide CLFLUX3 dum = 'KOPS' MATRIK ; 'OUBLIER' dum ; 'FINSI' ; * On fabrique le terme de flux complet 'SI' (FLUNEU 'OU' FLUTOT 'OU' FLUMIX) ; CLFLUX = CLFLUX1 'ET' CLFLUX2 'ET' CLFLUX3 ; FLUCLI = VRAI ; 'FINSI' ; *--------------------------------------------------------------------- *---------- Initialisations de tables, coefficients ------------------ *--------------------------------------------------------------------- toltheta = 1.D-4 ; MODPOROS = TABMODI . 'POROSITE' ; MODCONV = TABMODI . 'CONVECTI' ; MODDELTT = TABMODI . 'DELTAT' ; MODCOLIN = TABMODI . 'COEF_LIN' ; MODDIFF = TABMODI . 'DIFFUSIV' ; * On ne recalcule pas a priori la matrice du probleme LCALMATP = FAUX ; * On ne recalculera pas a priori la matrice de convection du probleme LCALCONV = FAUX ; 'SI' ('EXISTE' QFACE) ; PCONV = VRAI ; 'SINON' ; PCONV = FAUX ; 'FINSI' ; 'SI' (PCONV) ; tetaconv = TbDarTra . 'THETA_CONVECTION' ; 'FINSI' ; * Si la porosité est modifiée (ou plus exactement le coef * devant dC/dt) 'SI' (MODPOROS) ; * Calcul du terme devant le dC/dt integré sur le volume * CHAMELEM 'SCAL'. CoefDt = ('NOMC' 'SCAL' Porosite) ; COefDt = Porosite * ('DOMA' modarcy VOLUME) ; CoefDt = 'KCHA' Modarcy Coefdt 'CHAM' ; TbDarTra . 'SURF' = CoefDt ; * on recalculera la matrice totale du probleme LCALMATP = VRAI ; 'FINSI' ; * * Si la convection est modifiée (vitesse modifiée, voir apparition * ou suppression d'un champ de vitesse) en cours de simulation. * 'SI' (PCONV 'ET' MODCONV) ; * TbDarTra . 'THETA_CONVECTION'= TetaConv ; * On recalculera la matrice de convection du probleme LCALCONV = VRAI ; 'FINSI' ; * * Si le pas de temps est modifié * 'SI' (MODDELTT) ; TbDarTra . 'PAS' = DeltaT ; * on recalculera la matrice totale du probleme LCALMATP = VRAI ; LCALCONV = VRAI ; 'FINSI' ; * * Si le coefficient de décroissance (ou linaire) en facteur de C est * modifié ????????????????, * 'SI' (MODCOLIN) ; * A priori sera géré dans le coeff de porosité 'FINSI' ; *--------------------------------------------------------------------- *--------------------- CALCUL DU DECENTREMENT ------------------------ *--------------------------------------------------------------------- compmat = 'EXTRAIRE' Modarcy 'MATERIAU' ; ANISO = ('DIME' compmat) > 1 ; * Seulement si présence de convection et option décentrée * Modifié seulement si la convection ou la diffusion changent (Péclet * de maille = u dx / diff) 'SI' (PCONV 'ET' MODCONV) ; * la convection change d'où recalcule de la matrice du problème * plus la matrice de convection. La matrice de diffusion numérique et * de dispersion change aussi. LCALCONV = VRAI ; LCALMATP = VRAI ; * les dispersivités changent d'où recalcule de matrice diffusion * plus la matrice du problème plus la matrice de convection * GBM IL RESTE A CODER LE RECALCUL DE LA DISPERSION (IDEM EN VF) * ET DU DECENTREMENT (fait en VF). 'FINSI' ; 'SI' (MODDIFF) ; * recalcul de la matrice complete du probleme + de la diffusion LCALMATP = VRAI ; zozo = DIFFANIS Matediff 'EFMH' ANISO ; Matediff = zozo ; 'FINSI' ; * Difftot calculé que si MODDIF ou dispersion change (a faire alors) SI (MODDIFF) ; Difftot = MateDiff '+' diffdisp ; difftot = 'KCHA' Modarcy difftot 'CHAM' ; difftot = CHANGER difftot TYPE CARACTERISTIQUES ; SI (ANISO) ; 'SI' ('EGA' ('VALEUR' 'DIME') 2) ; difftot = 'MATE' Modarcy 'DIRECTION' (1. 0. ) 'PARALLELE' 'K11' ('EXCO' 'K11' difftot) 'K21' ('EXCO' 'K21' difftot) 'K22' ('EXCO' 'K22' difftot) ; 'FINSI' ; 'SI' ('EGA' ('VALEUR' 'DIME') 3) ; difftot = 'MATE' Modarcy 'DIRECTION' (1. 0. 0.) (0. 1. 0.) 'PARALLELE' 'K11' ('EXCO' 'K11' difftot) 'K21' ('EXCO' 'K21' difftot) 'K22' ('EXCO' 'K22' difftot) 'K31' ('EXCO' 'K31' difftot) 'K32' ('EXCO' 'K32' difftot) 'K33' ('EXCO' 'K33' difftot) ; 'FINSI' ; SINON ; difftot = 'MATE' Modarcy 'K' difftot ; FINSI ; FINSI ; * 'SI' ('EGA' ('VALEUR' 'DIME') 2) ; * difftot = 'MATE' Modarcy 'DIRECTION' * (1. 0. ) 'PARALLELE' 'K11' ('EXCO' 'K11' difftot) * 'K21' ('EXCO' 'K21' difftot) * 'K22' ('EXCO' 'K22' difftot) ; * 'FINSI' ; * 'SI' ('EGA' ('VALEUR' 'DIME') 3) ; * difftot = 'MATE' Modarcy 'DIRECTION' * (1. 0. 0.) (0. 1. 0.) 'PARALLELE' * 'K11' ('EXCO' 'K11' difftot) * 'K21' ('EXCO' 'K21' difftot) * 'K22' ('EXCO' 'K22' difftot) * 'K31' ('EXCO' 'K31' difftot) * 'K32' ('EXCO' 'K32' difftot) * 'K33' ('EXCO' 'K33' difftot) ; * 'FINSI' ; *'FINSI' ; *--------------------------------------------------------------------- *-------------- Matrice masse inverse des EFMH ----------------------- *--------------------------------------------------------------------- * Calcul des matrices de masse elementaires inverses * affectées par un changement de la diffusivité totale. * Cela ce porduit si la diffusivité change ou, en présence * de décentrement, si la diffusivité numérique change. 'SI' (MODDIFF 'OU' (DECENTR 'ET' PCONV 'ET' MODCONV)) ; 'SI' (LMLump) ; * masse lumping MassEFMH = 'MHYB' MoDARCY Difftot 'LUMP' ; 'SINON' ; MassEFMH = 'MHYB' MoDARCY Difftot ; 'FINSI' ; LCALMATP = VRAI ; LCALCONV = VRAI ; 'FINSI' ; *---------------------------------------------------------------------- *---------- CALCUL DE LA MATRICE DE DIFFUSION DU PROBL ---------------- *---------------------------------------------------------------------- * Calcul de la matrice du probleme diffusion transitoire * Elle est modifiée si la matrice mass EFMH change, * ou si le pas de temps change, ou si le coefficient * en facteur de la dérivée temporelle change (porosité*retard), * ou si le coefficient du terme linéaire a(x,t) * C(x,t) change. * Pour des raisons pratiques on n'a pas stocké les matrices * de convection et de diffusion séparément. Si l'une est modifiée * on doit donc recalculer les deux matrices. Ce choix est discutable, * on gagne un facteur 3 en stockage (mat totale au lieu de mat totale * + diffusion + convection). Par ailleurs, le cout de l'assemblage + * resolution + créer une des deux matrice rend souvent négligeable * le recalcul inutile d'une des deux matrices précitées. On réduit * également le nombre d'entrées-sorties. 'SI' (LCALMATP 'OU' LCALCONV) ; Matrdiff = 'MATP' MoDARCY MassEFMH TbDarTra ; 'FINSI' ; *---------------------------------------------------------------------- *--- CALCUL DE LA MATRICE DE CONVECTION ET DES SECONDS MEMBRES -------- *---------------------------------------------------------------------- * Nombres d'espèces à gérer est fixé et ne peut changer au cours d'une * simulation. Il sera d'ailleurs mis dans le modele 'SI' (nbespece > 0) ; * pour un schéma en temps non euler implicite, il faut * la trace à l'instant précédent pour le second membre * pour la convection ou la diffusion. Elle est issue * des calculs précédents en EFMH. C'est impératif. * On ne prévoit pas de la recalculer car c'est très cher. 'REPETER' bloc1 nbespece ; * On extrait la composante de Cini, Tcini et de la source CCini = 'NOMC' 'H' ('EXCO' (extr &bloc1 nomespec) Cini) ; TCCini = 'NOMC' 'TH' ('EXCO' (extr &bloc1 nomespec) Tcini) ; 'SI' (nbsource > 1) ; SSource = 'NOMC' 'SOUR' ('EXCO' (extr &bloc1 nomespec) Chpsour) ; 'SINON' ; SSource = 'NOMC' 'SOUR' Chpsour ; 'FINSI' ; * Conditions initiales TbDarTra . 'CHARGE' = CCini ; TbDarTra . 'TRACE'= TCCini ; * Prise en compte du terme source et eventuellement * de la convection avec le schema centre 'SI' (PCONV) ; * convection 'SI' (TetaConv 'NEG' 0.0D0 toltheta) ; * schéma partiellement implicite, matrice MatConv * retirer pour voir Matconv à gauche si non nécessaire * ???????????? (avec test sur MODCCONV) MatConv SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTra SSource (NOMC 'FLUX' QFace) ; 'SINON' ; * schéma explicite, calcul du second membre uniquement SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTra SSource (NOMC 'FLUX' QFace) ; 'FINSI' ; 'SINON' ; * pas de convection, calcul du second membre restant SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTra SSource ; 'FINSI' ; * On reconstitue un champ de second membre 'SI' ((&bloc1) 'EGA' 1) ; SSMTr = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) SSMTr ; SMTR = SSMTr ; 'SINON' ; SSMTr = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) ('COPIER' SSMTr) ; SMTr = SSMTr 'ET' SMTr; 'FINSI' ; 'FIN' bloc1 ; * on stoque la matrice en assemblant la matrice de convection * du dernier calcul (elles sont toutes identiques). 'SI' (PCONV 'ET' (LCALMATP 'OU' LCALCONV)) ; 'SI' (TetaConv 'NEG' 0.0D0 toltheta) ; * calcul partiellement implicite pour la convection MatrTr = Matrdiff 'ET' MatConv ; 'DETRUIT' MatConv ; 'DETRUIT' Matrdiff ; 'SINON' ; MatrTr = Matrdiff ; 'FINSI' ; 'SINON' ; 'SI' (PCONV 'ET' (TetaConv 'NEG' 0.0D0 toltheta)) ; 'DETRUIT' MatConv ; 'FINSI' ; 'SI' (('NON' PCONV) 'ET' LCALMATP) ; MatrTr = Matrdiff ; 'FINSI' ; 'FINSI' ; 'FINSI' ; *--------------------------------------------------------------------- *------ Conditions aux limites de Flux (mixtes et Neumann) ----------- *--------------------------------------------------------------------- 'SI' (FLUCLI) ; SMTr = SMTR 'ET' CLFLUX ; 'FINSI' ; *--------------------------------------------------------------------- *------ Matrice sortante --------------------------------- ----------- *--------------------------------------------------------------------- * si pas de modif et matrice rigidité en entrée 'SI' (('EXISTE' mattr) 'ET' ('NON' (LCALMATP 'OU' LCALCONV))) ; Matrtr = mattr; 'FINSI' ; * si pas de modif et matrice matrik en entrée 'SI' (('EXISTE' mattm) 'ET' ('NON' (LCALMATP 'OU' LCALCONV))) ; Matrtr = mattm; 'FINSI' ; 'FINP' SMTr MatrTr TbDarTra MassEFMH ;