* fichier : warrickEFMH.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. '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
*
*- 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' = 'EFMH' ;
* 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;
© Cast3M 2003 - All rights reserved.
Disclaimer