* fichier : warrickVF.dgibi ************************************************************************ * Section : Fluides Darcy ************************************************************************ ******************************************************************** * 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;