Procédures | Cast3M

* UPDAVF    PROCEDUR  GOUNAND   11/05/24    21:16:24     6976           
**********************************************************************
'DEBP'  UPDAVF MoDARCY*'MMODEL' Porosite*'CHPOINT' Matediff*'CHPOINT'
                 Diffdisp*'CHPOINT' ChPSour*'CHPOINT' Deltat*'FLOTTANT'
                 Cini*'CHPOINT' TetaDiff*'FLOTTANT' TetaConv*'FLOTTANT' 
                 Qface/'CHPOINT' nomespec*'LISTMOTS' nbespece*'ENTIER' 
                 nbsource*'ENTIER' mattr*'MATRIK' JACO*'MATRIK' 
                 MPOR2*'MATRIK'  MCHAMT*'MCHAML'  MCHAMT1*'MCHAML' 
                 TABMODI*'TABLE' CHCLIM*'TABLE' LCONV*'LOGIQUE';
* 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 Mattt Difftot Mctot Mdiff Nouvmat = UPDAVF
*                 MoDARCY Porosite Matediff Diffdisp ChPSour
*                 DeltaT Cini TetaDiff TetaConv
*                 Qface nomespec nbespece nbsource
*                 Matot Jaco Mctot Mdiff Mpor TABMODI CHCLIM           ;
*  |-----------------------------------------------------------------|
*  | Généralités : MATTVF 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'.                                      |
*  |                                                                 |
*  | 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   |
*  |                                                                 |
*  | Matot    : matrice globale de discretisation en VF              |
*  |                                                                 |
*  | Jaco     : matrice globale de discretisation en VF pour le problème
*  |          stationnaire
*
*  | Mpor     : matrice globale de discretisation en VF pour le problème
*  |          stationnaire
*  |                                                                 |
*  | Mchamt   : Coef permettant de calculer le flux total
*  |                                                                 |
*  | Mchamt1  : Coef permettant de calculer le flux diffusif
*  |                                                                 |
*  |                                                                 |
*  | 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                         |
*  |-----------------------------------------------------------------|
*  |                                                                 |
*  | Difftot  : Coefficient de diffusion totale, integre decentrement|
*  |                                                                 |
*  |                                                                 |
*  |-----------------------------------------------------------------|
*  |                         SORTIES                                 |
*  |-----------------------------------------------------------------|
*  |                                                                 |
*  |                                                                 |
*  | RESI     : second membre                                        |
*  |                                                                 |
*  | Matot    : matrice globale de discretisation en VF              |
*  |                                                                 |
*  | Difftot  : Coefficient de diffusion totale, integre decentrement|
*
*  | Mpor     : matrice globale de discretisation en VF pour le problème
*  |          stationnaire
*  |                                                                 |
*  | Mchamt   : Coef permettant de calculer le flux total
*  |                                                                 |
*  | Mchamt1  : Coef permettant de calculer le flux diffusif
*  |                                                                 |
*  |                                                                 |
*  |-----------------------------------------------------------------|
*  |                  VARIABLES INTERNES                             |
*  |-----------------------------------------------------------------|
*  |                                                                 |
*  |                                                                 |
*  | PCONV    : Logique indiquant VRAI si présence de convection     |
*  |                                                                 |
*  | CoefDt   : coeff devant dC/dt integre sur les elements          |
*  |                                                                 |
*  | CCini    : concentration aux centres (une composante)           |
*  |                                                                 |
*  | SSource  : Source aux centre (une composante)                   |
*  |                                                                 |
*  | SSMTr    : second membre sur les traces pour une espèce         |
*  |                                                                 |
*  | 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 (VF) 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;
DIRCLI = FAUX;
FLUTOT = FAUX;
FLUMIX = FAUX;
CLFLUX MA = 'KOPS' 'MATRIK';
CLDIRI MA = 'KOPS' 'MATRIK';
CLFLUT MA = 'KOPS' 'MATRIK';
EXFLU = CLFLUX;
EXDIR = CLDIRI;
EXFLUT = CLFLUT;
* Neumann
'SI' ('EXISTE' CHCLIM 'NEUMANN') ;
   CLFLUX = CHCLIM . 'NEUMANN';
   FLUNEU = VRAI;
'FINSI';
 
'SI' ('EXISTE' CHCLIM 'DIRICHLET') ;
   CLDIRI = CHCLIM . 'DIRICHLET';
   DIRCLI = VRAI;
'FINSI';
 
'SI' ('EXISTE' CHCLIM 'FLUTOTAL') ;
   CLFLT = CHCLIM . 'FLUTOTAL';
   FLUTOT = VRAI;
'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                                                      ;
   mayage = 'EXTRAIRE' CHCLIM . 'FLUMIXTE' . 'VAL' maillage           ;
   cofb   = (doma modarcy SURFACE) *CHCLIM . 'FLUMIXTE' . 'COEFB'    ;
   cofb   = CHCLIM . 'FLUMIXTE' . 'COEFB'    ;
   cofa   = 'REDU' (CHCLIM . 'FLUMIXTE' . 'COEFA')  mayage            ;
   cofb   = 'REDU' cofb  mayage                                       ;
   coefm  = (1.D0 * cofb) '/' cofa                                   ;
   CLFLUX3   = 'REDU' CLFLUX3  mayage                                 ;
   'OUBLIER' cofa                                                     ;
   'OUBLIER' cofb                                                     ;
'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'                                       ;
 
 
 
 
* Si la porosité est modifiée (ou plus exactement le coef
* devant dC/dt)
 
 
'SI' (MODPOROS)                                                       ;
MPOR2 = 'KOPS' 'MATDIAGO' ('NOMC' 'RETN' Porosite) 'MATRIK'           ; 
'FINSI'                                                               ;
 
*
* Si le pas de temps est modifié
*
 
 
*
* 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 DE LA MATRICE DE DIFFUSION DU PROBL ----------------
*----------------------------------------------------------------------
 
 
* Calcul de la matrice du probleme diffusion transitoire
* Elle est modifiée si la diffusion  ,
* 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' (MODDIFF)                                                        ;
*   La diffusivité effective a changé, on remmet la matrice au propre
*   c'est à dire avec des noms de composantes isotropes K11 K22 etc ...
    zozo = DIFFANIS Matediff                                          ;
    Matediff = zozo                                                   ;
'FINSI'                                                               ;
 
Difftot = MateDiff '+' Diffdisp                                       ;
 
 
*----------------------------------------------------------------------
*--- 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)                                                   ;
 
 
   DTI = 1.D0/DeltaT                                                  ;
   'REPETER' bloc1 nbespece ;
 
      espc =  ('EXTRAIRE' &bloc1 nomespec) ;
 
*     On extrait la composante de Cini, Tcini et de la source
      CCini = 'NOMC' 'H' ('EXCO' espc Cini)         ;
      'SI' (nbsource > 1)                                             ;
         SSource = 'NOMC' 'SOUR' ('EXCO' espc
                               Chpsour)                               ;
      'SINON'                                                         ;
         SSource = 'NOMC' 'SOUR' Chpsour                              ;
      'FINSI'                                                         ;    
*     Conditions initiales
*     Prise en compte du terme source et eventuellement 
*     de la convection avec le schema centre 
       HHS  = 'NOMC' 'SCAL'  CCini                                    ;
* Prise en compte des conditions aux limites 
       'SI' (FLUNEU)                                                  ;
             EXFLU =  
              ('NOMC' 'FLUX' ('EXCO' espc
              CLFLUX))                                                ;
*              MESS 'EXFLU ='; LIST EXFLU;
       'FINSI'                                                        ;
       'SI' (DIRCLI) ;
             EXDIR =  
             ('NOMC' 'SCAL' ('EXCO' espc
             CLDIRI))                                                 ;
*             MESS 'EXDIR ='; LIST EXDIR;
       'FINSI'                                                        ;
        'SI' (FLUTOT) ;
               EXFLUT =
               ('NOMC' 'FLUX' ('EXCO' espc
               CLFLT ))                                               ;
               QLIM = NOMC(EXFLUT) 'FLUX'                             ;
               SUP = 'EXTR' EXFLUT  'MAIL'                            ;
               XPAR1 = 1.D0 + (0.0*CLFLT)                             ;
               XPAR1  = 'NOMC' XPAR1 'PAR1'                           ;
               'SI' LCONV                                             ;
               USCNR = 'REDU' QFACE SUP                               ;
               MUSCN = USCNR*(-1.D0)                                  ;
               'SINON'                                                ;
               MUSCN =  (0.D0*CLFLT)                                  ;
               'FINSI'                                                ;
 
               XPAR2  = 'NOMC' MUSCN 'PAR2'                           ;
               EXFLUT = XPAR1 + XPAR2 + QLIM                          ;
             EXFLUT = CHAN 'ATTRIBUT' EXFLUT NATURE DISCRET           ;
*               MESS 'EXFLUT ='; LIST EXFLUT;
         'FINSI'                                                      ;
 
            'SI' (FLUMIX) ;
               EXFLUM =
               ('NOMC' 'FLUX' ('EXCO' espc
               CLFLUX3 ))                                             ;
               QLIM = NOMC(EXFLUM) 'FLUX'                             ;
               SUP = 'EXTR' EXFLUM  'MAIL'                            ;
               XPAR1 = 1.D0 + (0.0*CLFLUX3)                           ;
*              MESS 'XPAR1 =';
               XPAR1  = 'NOMC' XPAR1 'PAR1'                           ;
               MUSCN = coefm;
               XPAR2  = 'NOMC' MUSCN 'PAR2'                           ;
               EXFLUM = XPAR1 + XPAR2 + QLIM                          ;
              EXFLUM = CHAN 'ATTRIBUT' EXFLUM NATURE DISCRET          ;
               EXFLUT = EXFLUT 'ET'  EXFLUM                            ;
            'FINSI'                                                   ;
*              MESS 'UPDAVF EXFLUT ='; LIST EXFLUT;
*              MESS 'UPDAVF EXFLU ='; LIST EXFLU;
 
*        'MESS' 'MODDIFF'; LIST MODDIFF; 
*        'MESS' 'MODCONV'; LIST MODCONV; 
          'SI' (MODDIFF 'OU' MODCONV)                                 ;
* On recalcule deux fois les coefficients. En général, cela doit etre
* négligeable en temps de calcul; 
            GRADT0 MCHAMT1 = 'PENT'  MoDARCY 'FACE' 'MPFA'
            HHS 'DISPDIF' difftot 'TIMP' EXDIR 'QIMP' EXFLU            
            'MIXT' EXFLUT                                             ;
            MCHAMT = MCHAMT1;
            MCHAJA = TetaDiff*MCHAMT1;
 
            'SI' LCONV                                                ;
            GRADT0 MCHAMT = 'PENT'  MoDARCY 'FACE' 'MPFA'
            HHS 'DISPDIF' difftot   'TIMP' EXDIR 'QIMP' EXFLU 
            'MIXT' EXFLUT 'UPWICENT' QFace                            ;
* ATTENTION PEUT ETRE  FAUX SI TetaConv different de TetaDiff
             MCHCONV = MCHAMT + ((-1.D0)*MCHAMT1)                     ;
             MCHAJA = (TetaConv*MCHCONV) + (TetaDiff*MCHAMT1)         ;
             'DETRUIT' MCHCONV                                        ;
            'FINSI'                                                   ;
 
 
* On ne calcule la jacobienne que pour la première espèce 
         'SI' ((&bloc1) 'EGA' 1)                                      ;
         JACO CHPRES DT = 'LAPN' 'VF' 'CLAUDEIS'  'IMPL'
          MoDARCY HHS GRADT0 MCHAJA 'QIMP' EXFLU 'MIXT' EXFLUT 
          'TIMP' EXDIR                                                ; 
          'SINON'                                                     ;
          JACOB CHPRES DT = 'LAPN' 'VF' 'CLAUDEIS'  'EXPL'
          MoDARCY HHS GRADT0                                          ;
         'FINSI'                                                      ;         
*       SINON;
*         'DETRUIT' MCHAJA                                            ; 
      'SI' ( ((&bloc1) 'EGA' 1) )        ;
              MATOT = 'KOPS' 'MULT' -1.000D0 JACO                     ;     
              MATI = '*'  MPOR2         DTI                           ;
              MATOT = MATOT 'ET' MATI                                 ; 
      'FINSI'                                                         ;
      'SINON'; 
      'SI' ( ((&bloc1) 'EGA' 1) 'ET' (MODPOROS 'OU' MODDELTT))        ;
              MATOT = 'KOPS' 'MULT' -1.000D0 JACO                     ;     
              MATI = '*' MPOR2   DTI                                  ;
              MATOT = MATOT 'ET' MATI                                 ; 
      'FINSI';
*       'SINON';
*     On reconstitue un champ de second membre
            GRADT0  = 'PENT'  MoDARCY 'FACE' 'MPFA'
            HHS 'DISPDIF' difftot
            'TIMP' EXDIR 'QIMP' EXFLU 'MIXT' EXFLUT 'UPWICENT' QFace 
            'GRADGEO' MCHAMT                                           ;
          JACO CHPRES DT = 'LAPN' 'VF' 'CLAUDEIS'  'EXPL'
          MoDARCY HHS GRADT0                                           ;
 
     'FINSI'                                                           ;
 
      'SI' ((&bloc1) 'EGA' 1)                                          ;
         SSMTr = 'NOMC' espc CHPRES            ;
         RESI =  SSMTr                                                 ;
      'SINON'                                                          ;
         CHPRES = 'NOMC' espc 
         ('COPIER' CHPRES) ;
         RESI = CHPRES 'ET'  RESI; 
      'FINSI'                                                          ;
   'FIN' bloc1                                                         ;
 
'FINSI'                                                                ;
 
 
*---------------------------------------------------------------------
*------ Conditions aux limites de Flux (mixtes et Neumann) -----------
*---------------------------------------------------------------------
 
 
*---------------------------------------------------------------------
*------ Matrice sortante --------------------------------- -----------
*---------------------------------------------------------------------
 
 
* si pas de modif et matrice matrik en entrée
NOUVMAT = VRAI;
'SI' (('EXISTE' mattr) 'ET' ('NON' (MODDIFF 'OU' MODCONV 'OU' 
                                    MODPOROS 'OU' MODDELTT)))         ;
   Matot = mattr                                                      ;
   NOUVMAT = FAUX                                                     ;
'FINSI'                                                               ;
 
 
'FINP' RESI Matot Difftot MCHAMT MCHAMT1 NOUVMAT                      ;