Accueil | Cast3M

* fichier : warrickVF.dgibi
************************************************************************
************************************************************************
********************************************************************
*
GRAPH = FAUX ;
*
******************** CAS TEST : warrick.dgibi **********************
*
* Test de fonctionnement de DARCYSAT en 2D avec effet de gravité
* en regime permanent.
*
*-------------------------------------------------------------------
*
*     A. GENTY, G. BERNARD-MICHEL - DM2S/SFME/MTMS - 02/2006
*
* ==================================================================
* Infiltration d'eau depuis une source ponctuelle placee dans un
* milieu 2D non sature limite par des surfaces inferieures et
* superieure a pression d'eau (h<0) imposees.
*
*                   h = h1 (charge H = h1 + z0)
*            z=z0   ---------------------------
*                   I                         I
*                   I                         I
*            Q1  -> I                         I    
*         (z=z0-d)  I                         I   
*                   I hi = (z/z0)(h1-h0) + h0 I
*                   I  (charge Hi = hi + z)   I
*         flux nul  I                         I flux nul
*                   I                         I
*                   I                         I
*                   I                         I
*                   I                         I
*                   I                         I
*                   I                         I
*             z=0   ---------------------------
*                x=0  h = h0 (charge H = h0)   x=x0
*
* Ce probleme est du meme type que le Probleme V9 decrit dans la 
* notice de validation v 2.50 de Porflow.
*
* La solution analytique (en pression) du probleme est obtenue en 
* suivant la demarche de Warrick et Lomen, 1977, 
* Soil Sci. Soc. Am. J., Vol 41, pp 849-851.
* ==================================================================
*
* - conditions initiales : 
*   il s'agit de trouver la solution d'un regime permanent, les 
*   conditions initiales sont donc a priori sans influence. On 
*   choisira cependant, dans l'objectif d'accelerer la convergence
*   du calcul, la pression initiale dans le domaine sous la forme
*   hi = (z/z0)(h1-h0) + h0 (charge initiale Hi = hi + z).
*
* - conditions aux limites :
*   > la pression sur la frontiere inferieure du domaine est imposee
*     a h = h0 (charge H = h + z = h0)
*   > la pression sur la frontiere superieure du domaine est imposee
*     a h = h1 (charge H = h + z = h1 + z0)
*   > Un flux entrant Q1 est imposé sur une surface d'extention 2*r0
*     centree a la cote z = z0 - d sur la frontiere gauche du domaine
*   > les autres faces exterieures du domaine sont imposees a flux nul
*
* ===================================================================
*    Les options de modélisation declarées dans la table transmise à
*    la procédure DARCYSAT sont les suivantes :
*
* -  les effets gravitationnels sont pris en compte : 
*                   FORCE_GRAVITE = 1.
*
* -  Le calcul est réalisé avec l'option de calcul de pas de temps
*    automatique avec dt_initial = 10000 s sur une duree totale
*    de 100000 s (regime permanent).
*
* -  Les options numeriques testees sont les suivantes:
*       EFMH  CENTRE    TYPINV = 3  PRECOND = 5
*       EFMH  DECENTRE  TYPINV = 3  PRECOND = 5
*       VF    CENTRE    TYPINV = 3  PRECOND = 5
*       VF    DECENTRE  TYPINV = 3  PRECOND = 5
*
* -  La précision de convergence demandée est de 1e-03 m
*
* ===================================================================
*                          RESULTATS ATTENDUS
*
* Le permanent de pression sera compare a la solution analytique
* par l'intermediaire des erreurs L_infini et L_2.
* Une erreur inferieure a 1.5 % est attendue.
*
* ===================================================================
****************** CAS TEST : warrick.dgibi **********************
*
'OPTION' 'ECHO' 1 ;
'SAUTER' 'PAGE';
*
'TITRE' 'Infiltration source ponctuelle : warrick'   ;
OPTI DIME 2 ELEM QUA4   ;
OPTI ISOV SULI ;
*OPTI TRAC PSC ;
*
********************************************************************
*  Fonction qui calcule la perméabilité en fonction de la
*  saturation (loi exponentielle)
************************************************************************
'DEBPROC' KR_EXPO1 LOI*'TABLE' SAT/'CHPOINT' TEST/'MOT' ;
 
   'SI' ('NON' ('EXISTE' TEST)) ;
      TEST = 'MOT' 'OK' ;
   'FINSI' ;
 
*- Vérification des paramètres
   'SI' ('NEG' TEST 'NOTEST') ;
      NOMT = 'MOT' LOI.'NOMZONE' ;
 
      'SI' ('NON' ('EXISTE' LOI 'ALPHA' )) ;
         'ERREUR' 'Il manque le coefficient ALPHA dans'
                  ' la loi de perméabilité de ' NOMT ;
         'QUITTER' KR_PRO ;
      'FINSI' ;
      'SI' ('NON' ('EXISTE' LOI 'PERMSAT' )) ;
         'ERREUR' 'Il manque le coefficient PERMSAT dans'
                  ' la loi de perméabilité de ' NOMT;
         'QUITTER' KR_PRO ;
      'FINSI' ;
   'FINSI' ;
*
    se1 = 'KOPS' SAT '**' -1. ;
    se2 = 'KOPS' se1 '-' 1. ;
    se3 = 'KOPS' se2 '*' -1. ;
    K11 = 'KOPS' se3 '*' LOI.'ALPHA' ;
    K12 = EXP(K11) ;
    K2 = 'KOPS' K12 '*' LOI.'PERMSAT' ;
 
    'DETRUIT' K11 ;
    'DETRUIT' K12 ;
    'DETRUIT' se1 ;
    'DETRUIT' se2 ;
    'DETRUIT' se3 ;
 
'FINPROC' K2 ;
*-------------------------------------------------------------------
*-------------------------------------------------------------------
*--------------------- Création du maillage 2D ---------------------
*
*- Densités de mailles
NDX1 = 20 ;
NDY1T = 6 ;
NDY12S = 8 ;
NDY1B = 35 ;
*
NDY1 = NDY1T + NDY12S + NDY1B ;
*
*- Cotes du domaine rectangulaire (x0 et z0)
XCA1 = 0.61 ;
YCA1 = 1.22 ;
*
*- Profondeur du point d'injection (d)
YDD1 = 0.15 ;
*- demi surface d'injection 
R0 = 0.05 ;
*
*- Création des points du domaine
A1 = 0. 0.     ;
A2 = XCA1 0.   ;
A3 = XCA1 YCA1 ;
A4 = 0. YCA1   ;
*
*- Point d'injection
YCAZ1 = YCA1 - YDD1 ;
A5T = 0. (YCAZ1 + R0) ;
A5B = 0. (YCAZ1 - R0) ;
*
* densites
DMINI1 = (2.D0 * R0) / NDY12S ;
DMOY1 = (YCA1 - YDD1 - R0) / NDY1B ;
DMOY2 = (YDD1 - R0) / NDY1T ;
ALP1 = DMINI1 / DMOY1 ;
ALP2 = DMINI1 / DMOY2 ;
DMAX1 = DMOY1 / ALP1 ;
DMAX2 = DMOY2 / ALP2 ;
*
*- Création des lignes
LIN1 = DROIT NDX1 A1 A2 ;
LIN2 = DROIT NDY1 A2 A3 ;
LIN3 = DROIT NDX1 A3 A4 ;
LIN4T = DROIT (-1 * NDY1T) A4 A5T 'DINI' DMAX2 'DFIN' DMINI1 ;
LIN412S = DROIT NDY12S A5T A5B ;
* Attention avec l'operateur droit qui gere comme il peut lorsque
* N (négatif), DINI et DFINI sont donnés et non nécessairement compatibles
DMAX1 = DMAX1 / 2.D0 ;
LIN4B = DROIT (-1 * NDY1B) A5B A1 'DINI' DMINI1 'DFIN' DMAX1 ;
LIN4 = LIN4T ET LIN412S ET LIN4B ;
*
*- Face pour imposition du flux entrant
LSOURCE = COUL VERT LIN412S ;
*
*- Face droite
LDROIT = LIN2 ;
*
*- Face gauche
LGAU = LIN4T ET LIN4B ;
*
*- Création du domaine
MASSIF0 = DALLER LIN1 LIN2 LIN3 LIN4 PLAN   ; 
LHAUT = LIN3 COUL BLEU ;
LBAS = LIN1 COUL ROUG ;
*
*-Tracé du maillage avec faces sup et inf pour conditions limites
* et surface d'entree du debit
ELIM 0.00001 (MASSIF0 ET LBAS ET LHAUT ET LGAU ET LSOURCE) ;
SI GRAPH;
TRAC CACH (MASSIF0 ET LBAS ET LHAUT ET LGAU ET LSOURCE);
FINSI;
*
*
***************************************************************************
***************************************************************************
*
*- Création des maillages contenant tous les points
QFTOT  = 'CHANGER' MASSIF0 'QUAF' ;
QFSORT = 'CHANGER' LBAS    'QUAF' ;
QFSORF = 'CHANGER' LHAUT    'QUAF' ;
QFENT  = 'CHANGER' LSOURCE   'QUAF' ;
*
'ELIMINATION' 0.00001 (QFTOT 'ET' QFSORT 'ET' QFENT 'ET' QFSORF) ;
*
*-------------------------------------------------------------------
*----------------------- Création du Modèle ------------------------
*
MODHYB  = 'MODELE' QFTOT  'DARCY' 'ISOTROPE'        ;
MODSORT = 'MODELE' QFSORT 'DARCY' 'ISOTROPE'        ;
MODSORF = 'MODELE' QFSORF 'DARCY' 'ISOTROPE'        ;
MODENT  = 'MODELE' QFENT  'DARCY' 'ISOTROPE'        ;
CESORT  = 'DOMA'  MODSORT 'CENTRE' ;
CESORF  = 'DOMA'  MODSORF 'CENTRE' ;
CEENT   = 'DOMA'  MODENT  'CENTRE' ;
HYCEN   = 'DOMA'  MODHYB  'CENTRE' ;
HYFAC   = 'DOMA'  MODHYB  'FACE'   ;
HYSOM   = 'DOMA'  MODHYB  'SOMMET' ;
*
ZCC = COOR 2 HYCEN ;
ZFF = COOR 2 HYFAC ;
*-------------------------------------------------------------------
*--------------- Conditions aux limites (blocages) ------------------
*
*- surfaces inferieure et superieure
* valeur imposee de pression P = HN0 et HN1
HN0 = -0.387 ;
HN1 = HN0 ;
*HN1 = -0.500 ;
*
* on impose la charge H = P + z
ESORT = 'MANU' 'CHPO' CESORT 1 'TH' HN0 'NATURE' 'DISCRET';
ESORF = 'MANU' 'CHPO' CESORF 1 'TH' (HN1 + YCA1) 'NATURE' 'DISCRET';
ESORT = ESORT + ESORF;
*
*--------------------------------------------------------------------
*- Chargement des conditions à la limite
*
LICALC = 'PROG' 0.D0  1.D8 ;
LIST1 = 'PROG' 2 * 1.D0   ;
VALI0 = 'CHAR'  'THIM' (ESORT)  ('EVOL' 'MANU' LICALC LIST1) ;
*
*- Flux impose
Q1 = -5.25E-08 ;
QCA1 = Q1 / NDY12S ; 
*EENT = 'MANU' 'CHPO' CEENT 1 'FLUX' Q1 'NATURE' 'DISCRET';
EENT = 'MANU' 'CHPO' CEENT 1 'FLUX' QCA1 'NATURE' 'DISCRET';
VALI1 = 'CHAR' 'FLUX' (EENT)  ('EVOL' 'MANU' LICALC LIST1) ;
*
*--------------------------------------------------------------------
*----------------- initialisation des inconnues --------------------- 
*
*- charge initiale (en metres d'eau) dans le domaine (H = P + z)
* Pression h
PM = ZCC '*' ((HN1 - HN0)/YCA1) '+' HN0 ;
* Charge H
HM = PM '+' ZCC;
*
H0   =  nomc (HM) 'H' 'NATURE' 'DISCRET';
*
*- Flux
QFACE0 = MANU CHPO HYFAC 1 'FLUX' 0.D0 'NATURE' 'DISCRET' ;
*
*--------------------------------------------------------------------
*----------------------------- Calcul -------------------------------
*
*                                         ---------------------------
*                                         = Table DARCY_TRANSITOIRE =
*                                         ---------------------------
*- initialisation table
SATUR                     = 'TABLE' ;
SATUR. 'TEMPS'            = 'TABLE' ;
SATUR. 'TRACE_CHARGE'     = 'TABLE' ;
SATUR. 'CHARGE'           = 'TABLE' ;
SATUR. 'FLUX'             = 'TABLE' ;
*
*- données géométriques
SATUR. 'SOUSTYPE'         = 'DARCY_TRANSATUR' ;
SATUR. 'MODELE'           = MODHYB            ;
*
*- instant initial
SATUR. 'TEMPS' . 0        = 0.                ;
SATUR. 'CHARGE' . 0       = H0                ;
SATUR. 'FLUX' . 0         = QFACE0            ;
*
*- flux impose 
SATUR. 'FLUX_IMPOSE'      = VALI1            ;
***********************************************************
*                SEULS GROS CHGTS GBM
***********************************************************
* GBM - optionnel calculée si absente
*SATUR. 'TRACE_CHARGE' . 0 = TH0               ;
*
*- conditions aux limites et chargements
* GBM - MODIFIE
*SATUR. 'BLOCAGES_DARCY'   = BSORT  ;
*SATUR . 'CHARGEMENT'       = VALI0  ;
SATUR . 'TRACE_IMPOSE' = VALI0 ;
*
* GBM - NOUVEAU
SATUR . 'LUMP'   = FAUX;
SATUR . 'TYPDISCRETISATION'   = 'VF' ;
*SATUR . 'DECENTR'  = FAUX;
* gravite 
SATUR . 'FORCE_GRAVITE' = 1.D0 ;
************************************************************
*- données physiques
************************************************************
* loi de perméabilité (Cf procedure en debut de fichier)
LoiP                      = 'TABLE' 'PERSONNELLE';
LoiP. 'ALPHA'             = 12.58 ;
LoiP. 'PERMSAT'           = 1.12E-5 ;
LoiP. 'NOM_PROCEDURE'     = 'MOT' 'KR_EXPO1' ;
SATUR.'LOI_PERMEABILITE'  = 'TABLE' LoiP ;
*
*****************************************************************
* loi de succion simple
LoiS                      = 'TABLE' 'VAN_GENUCHTEN';
LoiS. 'PORO'              = 0.1 ;
LoiS. 'TERESIDU'          = 0.001 ;
LoiS. 'NEXP'              = 1.0 ;
LoiS. 'MEXP'              = 1.0  ;
LoiS. 'BHETA'             = 1.0 ;
SATUR.'LOI_SATURATION'    = 'TABLE' LoiS ;
*
*****************************************************************
*- données numériques
*
TF1 = 1.D5 ;
SATUR. 'TEMPS_FINAL'     = TF1 ;
SATUR. 'DT_INITIAL'      = 1.D4 ;
SATUR. 'ITMAX'           = 30;
SATUR. 'HOMOGENEISATION' = 'CHAINE' 'CENTRE' ;
SATUR. 'RESIDU_MAX'      = 1.D-3 ;
SATUR. 'NITER'           = 40 ;
SATUR. 'MESSAGES'        = 2 ;
*
*****************************************************************
*------------------------- Exécution procédure ----------------------
*
TABRES = table METHINV;
TABRES . 'TYPINV' = 3;
TABRES . 'PRECOND' = 5;
*
SATUR . 'METHINV' = TABRES;
SATUR . 'METHINV' . RESID = 1.D-15;
*
*********************************************************************
*
DARCYSAT SATUR ;
*
*********************************************************************
*
* Date de trace  
T11 = TF1 ;
*
*-Légende des tracés
TABI1 = TABLE ;
TABI1.'TITRE'= TABLE ;
TABI1.1='MARQ LOSA' ;
TABI1.'TITRE' . 1  = MOT 'permanent' ;
*
*-------------------------------------------------------------------
*-------------------------- Extration des charges ---------------------
T11 = SATUR . TEMPS . (('DIME' SATUR . TEMPS) '-' 1);
CALCTH1 = PECHE SATUR CHARGE T11 ;
*
*-------------------------------------------------------------------
*------------------- Calcul et trace de la pression ----------------
CPRE1 = CALCTH1 '-' ('NOMC' 'H' ZCC) ;
CALOO1 = KCHA MODHYB CPRE1 'CHAM';
savcal = 'NOMC' SCAL CPRE1 ;
* Trace sur elements
'SI' GRAPH  ;
TRAC MODHYB CALOO1 ;
'FINSI' ;
* Trace sur noeuds
CALOO2 = ELNO MODHYB CPRE1 ;
CONT1 = CONTOUR MASSIF0 ;
'SI' GRAPH  ;
TRAC CALOO2 (MASSIF0) CONT1 ;
'FINSI' ;
*
*---------------------------------------------------------------------------
* =========================================================================
* Solution analytique pour
*      une pression h(m) imposee a h0 sur la frontiere inferieure
*      une pression h(m) imposee a h1 sur la frontiere superieure
*      flux nul impose sur la frontiere droite
*      flux impose a Q1 sur la partie d-r0 < z < d+r0 de la frontiere gauche
*                  a flux nul sur le reste de la frontiere gauche
* =========================================================================
* Attention, la solution analytique est calculee en pression dans un repere
* ou l'axe z est dirige vers le bas avec 
*    Xnew = (x*alpha)/2 et Znew = ((z0-z)*alpha)/2
*---------------------------------------------------------------------------
*
* recuperation des coordonnees X et Z aux centres
XCCZ1 = COOR 1 HYCEN ;
ZCCZ1 = COOR 2 HYCEN ;
*
* constantes utiles
KZER1 = LoiP. 'PERMSAT' ;
ALPHA1 = ABS((LoiP. 'ALPHA')) ;
ALPHA2 = ALPHA1 / 2.D0 ;
HZER1 = HN0 ;
HZUN1 = HN1 ;
DELTA1 = ALPHA2 * R0 ;
ZZER1 = YCA1 * ALPHA2 ;
XZER1 = XCA1 * ALPHA2 ;
D111 = YDD1 * ALPHA2 ;
QI1 = Q1 / (2.D0 * DELTA1) ;
*
* X et Z addimentionnels (et reorientation de z)
XCC1 = XCCZ1 * ALPHA2 ;
ZCC1 = (YCA1 - ZCCZ1) * ALPHA2 ;
*
*----------------------------------------------------------------------------
* calcul de phi1(Z) 
* (solution en Znew sans source pour pressions imposee en haut et en bas)
* phi1(Z) = (K0.exp(alpha.h1))/alpha * 
*           ( 1 + (exp(alpha(h0-h1)) -1)(exp(2Z) -1)/(exp(2/Z0) - 1))
*----------------------------------------------------------------------------
NUM1 = KZER1 * (EXP(ALPHA1*HZUN1)) ;
CTE1 = NUM1 / ALPHA1 ;
LIST CTE1 ;
*
DD3 = 2.D0 * ZZER1 ;
NUM2 = (EXP((ALPHA1*(HZER1 - HZUN1))) - 1.D0) / (EXP(DD3) - 1.D0) ;
*
DD1 = KOPS ZCC1 '*' 2.D0 ;
DD2 = EXP(DD1) ;
DD4 = KOPS DD2 '-' 1.D0 ;
DD5 = KOPS DD4 '*' NUM2 ;
DD6 = KOPS DD5 '+' 1.D0 ;
DD7 = KOPS DD6 '*' CTE1 ;
*
* trace de phi1
*PHI1E = KCHA MODHYB DD7 'CHAM';
*TRAC MODHYB PHI1E ;
*PHI1N = ELNO MODHYB DD7 ;
*TRAC MASSIF0 PHI1N ;
*
* calcul de h1 associe a phi1
CC1 = ALPHA1 / KZER1 ;
VPH1 = DD7 * CC1 ;
VPH2 = LOG(VPH1) ;
VCH1 = VPH2 / ALPHA1 ;
*
* trace de pression h1
*H1EL = KCHA MODHYB VCH1 'CHAM';
*TRAC MODHYB H1EL ;
*H1NO = ELNO MODHYB VCH1 ;
*TRAC MASSIF0 H1NO ;
*
*----------------------------------------------------------------------------
*
*----------------------------------------------------------------------------
* calcul de phi2(X,Z)
* phi2(X,Z) = Exp(Z-D) * 
*             Somme_n (Cn * sin(Mu_n.Z) * (Exp(L_n.(X-2X0)) + Exp(-L_n.X)))
* avec
*    Mu_n = n.Pi/Z0
*    L_n  = Sqrt(1 + Mu_n**2)
*    Cn   = 2.Q1/(Z0.L_n.((n.Pi/Z0)**2 + 1)) * ((Exp(-2.L_n.X0) - 1)**(-1)) * 
*   (Exp(R0) * (n.Pi/Z0 * Cos(n.Pi/Z0 * (D-R0)) + Sin(n.Pi/Z0 * (D-R0))) - 
*    Exp(-R0) * (n.Pi/Z0 * Cos(n.Pi/Z0 * (D+R0)) + Sin(n.Pi/Z0 * (D+R0))))
*
* Remarque: l'indice n=0 n'est pas calcule dans la somme car alors Mu_n=0
*           et sin(Mu_n.Z)=0
*----------------------------------------------------------------------------
* constantes trigo pour calcul
PI1 = PI ;
CSD1 = 180.D0/PI1 ;
* 
CSTE1 = PI1 / ZZER1 ;
D2XZER1 = 2.D0 * XZER1 ;
*
NBOUC1 = 1000 ;
NNN1 = 1.D0 ;
SOM1 = 0.D0 ;
*---------------------------------------------
REPETER BOU1 NBOUC1 ;
*
* boucle de calcul des mu_n, L_n, sin(Mu_n.Z),
* Exp(L_n.(X-2X0)) + Exp(-L_n.X) et Cn a 
* partir de n=1.
* Calcul de la somme.
*
**********************
* Mu_n 
  MU1 = NNN1 * CSTE1 ;
*******************************
* L_n 
  MU2 = MU1*MU1 ;
  LM2 = 1.D0 + MU2 ;
  LM1 = LM2 ** 0.5D0 ;
*******************************
* calcul de SI1 = sin(Mu_n Z)
  MU1Z = KOPS ZCC1 '*' MU1 ;
* en degres
  MU2Z = MU1Z * CSD1 ;
  SI1 = SIN(MU2Z) ;
***************************************
* calcul des exponentielles en L_n X
  MLM1 = -1.D0 * LM1 ;
* expo - lamda X
  TEX1 = KOPS XCC1 '*' MLM1 ;
  NP1 = EXP(TEX1) ;
* expo lambda (X -2X0)
  NTEX1 = KOPS XCC1 '-' D2XZER1 ;
  TEX2 = KOPS NTEX1 '*' LM1 ;
  NP2 = EXP(TEX2) ;
  TOTEX2 = NP1 + NP2;
***************************************
* calcul de Cn
* calcul terme multiplicatif devant Cn
* on retire 2.Q1/Z0 de la somme
  CM1 = (LM1 * LM2) * (EXP((D2XZER1 * MLM1)) - 1.D0) ;
* calcul terme second membre
  PHA1 = MU1 * (D111 + DELTA1) ;
  PHA2 = MU1 * (D111 - DELTA1) ;
* en degres
  PHA1_D = PHA1 * CSD1 ;
  PHA2_D = PHA2 * CSD1 ;
* premier terme trigo
  PT1 = (MU1 * (COS(PHA1_D))) + (SIN(PHA1_D)) ;
* second terme trigo
  PT2 = (MU1 * (COS(PHA2_D))) + (SIN(PHA2_D)) ;
* terme trigo avec exponentielle delta
  TT1 = (PT2 * (EXP(DELTA1))) - (PT1 * (EXP((-1.D0 * DELTA1)))) ;
* total second membre
  CN1 = TT1 / CM1 ;
  MESS 'Indice N =' NNN1 'Mu_n =' MU1 'Lambda_n =' LM1 'Cn =' CN1;
***************************************
* calcul produit total et somme
  COEF1 = CN1 * TOTEX2 * SI1 ;
  SOM1 = SOM1 + COEF1 ;
  NNN1 = NNN1 + 1. ;
FIN BOU1 ;
***************************************
* fin de calcul du terme somme
**************************************************************
* multiplication par Exp(Z - D)
VEX1 = ZCC1 - D111 ;
VEX2 = EXP(VEX1) ;
EXT1 = (2.D0 * QI1) / ZZER1 ;
* Calcul final de phi2
VEX3 = VEX2 * SOM1 * EXT1;
**************************************************************
* fin de calcul de phi2
* *********************
* Trace phi2
*PHI2E = KCHA MODHYB VEX3 'CHAM';
*TRAC MODHYB PHI2E ;
*PHI2N = ELNO MODHYB VEX3 ;
*TRAC MASSIF0 PHI2N ;
*
******************************
* calcul de phi = phi1 + phi2
PHIT1 = DD7 + VEX3 ;
* Trace phi
*PHIT1E = KCHA MODHYB PHIT1 'CHAM';
*TRAC MODHYB PHIT1E ;
*PHIT1N = ELNO MODHYB PHIT1 ;
*TRAC MASSIF0 PHIT1N ;
*
* calcul de h associe a phi
CC1 = ALPHA1 / KZER1 ;
VPH1 = PHIT1 * CC1 ;
VPH2 = LOG(VPH1) ;
VCH1 = VPH2 / ALPHA1 ;
* Trace h sur elements
HFI1E = KCHA MODHYB VCH1 'CHAM';
'SI' GRAPH  ;
TRAC MODHYB HFI1E ;
'FINSI' ;
* Trace h sur noeuds
HFI1N = ELNO MODHYB VCH1 ;
'SI' GRAPH  ;
TRAC HFI1N (MASSIF0) CONT1 ;
'FINSI' ;
*
**************************************************************
* Calcul de l'erreur relative par element (en pression)
error = ((savcal '-' vch1) '/' vch1) 'ABS' ;
error = 'KCHA' modhyb error CHAM;
'SI' GRAPH  ;
'TRACER' modhyb error;
'FINSI' ;
*
* Erreur Linfini
NUME1 = ABS(savcal '-' vch1) ;
DENU1 = ABS(vch1) ;
DENU2 = MAXI (DENU1) ;
ERR1 = NUME1 / DENU2 ;
ER1C = KCHA MODHYB ERR1 'CHAM';
'SI' GRAPH  ;
TRAC MODHYB ER1C ; 
'FINSI' ;
*
ERINF1 = MAXI (ERR1) ;
MESS 'ERREUR L INFINI' ERINF1;
*
**********************
* Calcul de l'erreur L2 par element (en pression)
ERRL21 = NUME1 ** 2.D0 ;
ERRL22 = DENU1 ** 2.D0 ;
ERRL2 = (ERRL21 '/' ERRL22) ** 0.5D0 ;
ERRL3 = 'KCHA' MODHYB ERRL2 CHAM ;
'SI' GRAPH  ;
TRAC MODHYB ERRL3 ;
'FINSI' ;
*
* Calcul erreur L2
SURF1 = 'DOMA'  MODHYB  'VOLUME' ;
ERNUM1 = ERRL21 *  SURF1  ;
ERDEN1 = ERRL22 *  SURF1  ;
INT1 = RESU ERNUM1 ;
INT2 = RESU ERDEN1 ;
ERL2 = (MAXI (INT1 / INT2)) ** 0.5 ;
MESS 'ERREUR L2' ERL2 ;
*
**************************************************************
* Test de non regression
ERMAX1 = 0.015 ;
SI (ERL2 > ERMAX1) ;
   mess 'ERREUR ANORMALE' ERL2 '>' ERMAX1 ;
  'ERREUR' 5;
SINON ;
  'ERREUR' 0;
FINSI ;
*
*-------------------------------------------------------------------
*
FIN;