Télécharger testalea.dgibi

Retour à la liste

* fichier :  testalea.dgibi
************************************************************************
************************************************************************
GRAPH = faux ;
 
********************* CAS TEST : testalea.dgibi ************************
*
* Cas-Test vérifiant le bon fonctionnement de l'opérateur 'ALEA'
* avec différentes options.
* Test de la bonne statistique du champ résultat.
*
************************************************************************
*
* Paramètres d'entrée :
* ---------------------
 
'OPTION' echo 1          ;
'OPTION' 'DIME' 2 'ELEM' 'QUA4' ;
'OPTION' 'TRACER' PSC          ;
 
* Vérifie la syntaxe et les différentes possibilités d'usage ?
verisynt = vrai ;
 
* Vérifie la corection statistique des résultats ?
veriresu = vrai ;
 
* Longueur de corrélation :
portee = .7 ;
 
* Moyenne et Ecart-type du Logarithme népérien de perméabilité :
VALYM  = -4.;
SIGMAY = 1. ;
VARIAN = sigmay * sigmay ;
 
 
*======================================================================
*==                                                                  ==
*==                          PROCEDURES                              ==
*==                                                                  ==
*======================================================================
 
 
  'DEBPROC' FCT xtab*'TABLE' p*'LISTREEL' ;
*   pour la procédure AJUSTE
*   Calcul de la fonction objectif
*   tbfonc . 'F' . i   (fonction non-linéaire précédée d'un coef. inconnu n°i)
*   tbfonc . 'G'       (fonction non-linéaire non précédée d'un coef. inconnu)
*   y = q^i f_i(x,p) + g(x,p)
    x = xtab . 1 ;
    n = 'DIME' x ;
    b = 'EXTR' p 1;
    tbfonc     = 'TABLE';
    tbfonc.'F' = 'TABLE';
    tbfonc.'F' . 1 = 'EXP' (-1. / b * x);
  'FINPROC' tbfonc;
 
  'DEBPROC' DERI xtab*'TABLE' p*'LISTREEL';
*   pour la procédure AJUSTE
*   Calcul de la dérivée de la fonction objectif par rapport aux
*   paramètres non linéaires
*   tbderi . 'F' . i . k  (dérivée par rapport au param k de la fonction f_i)
*   tbderi . 'G' . k      (dérivée par rapport au param k de la fonction g)
*   dy_k / dp_k= q^i (df_i(x,p) / dp_k) + (dg(x,p) / dp_k)
    x = xtab . 1;
    n = 'DIME' x;
    b = 'EXTR' p 1;
    tbderi = 'TABLE';
    tbderi . 'F' = 'TABLE';
    tbderi . 'F' . 1  = 'TABLE';
    tbderi . 'F' . 1 . 1 =  x / b / b * ('EXP' (-1. * x / b)) ;
  'FINPROC' tbderi;
 
*======================================================================
*==                                                                  ==
*==                      PROGRAMME PRINCIPAL                         ==
*==                                                                  ==
*======================================================================
 
* ===========================================
* Test de la syntaxe et des différents usages
* ===========================================
 
'SI' verisynt ;
 
* Création maillage (grossier) :
* ------------------------------
* Nombre de mailles pour 1 lambda :
  Np = 1.;
 
* Dimensions :
  lx = 2. * portee;
  ly = 2. * portee;
  lz = 2. * portee;
 
* Discrétisation :
* le +10d-10 sert à compenser l'imprécision machine pour ne pas avoir
* des valeurs de nx et ny trop petites de 1
  nx = 'ENTIER' (lx / portee * Np + 0.0001) ;
  ny = 'ENTIER' (ly / portee * Np + 0.0001) ;
  nz = 'ENTIER' (lz / portee * Np + 0.0001) ;
 
* Points :
  p1 = 0. 0. ;
  p2 = lx 0. ;
  p3 = lx ly ;
  p4 = 0. ly ;
 
* Droites :
  p1p2   = 'DROIT' p1 p2 nx ;
  p2p3   = 'DROIT' p2 p3 ny ;
  p3p4   = 'DROIT' p3 p4 nx ;
  p4p1   = 'DROIT' p4 p1 ny ;
 
* Surface :
  carre  = 'SURFACE' (p1p2 'ET' p2p3 'ET' p3p4 'ET' p4p1) ;
 
* domaine
  DOMHYB = 'CHANGER' CARRE  'QUAF';
 
* modèle :
  MODHYB = 'MODELE' DOMHYB 'DARCY' 'ISOTROPE'        ;
 
* points centres :
  HYSOM = 'DOMA' MODHYB 'SOMMET';
  HYCEN = 'DOMA' MODHYB 'CENTRE';
 
* -- Début tests syntaxiques --
* Syntaxe de base
  YY  = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'GRAVITE'
               'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM
               'LAMBDA1' portee (0. 1.) ;
 
  YY = 'ALEA' 'BANDES_TOURNANTES' modhyb 'GRAVITE'
              'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
 
* Syntaxes alternatives
 
*- en 2D
* champ aléatoire dans une direction 1D
  YY  = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'GRAVITE'
               'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
 
  'SI' GRAPH ;
    'TITRE' 'Champ aleatoire unidirectionnel 2D, bandes tournantes';
    'TRACER' MODHYB YY ;
  'FINSI' ;
 
* champ aléatoire dans une direction 1D en spécifiant la direction
  YY  = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'GRAVITE'
               'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM
               'LAMBDA1' portee (0. 1.) ;
 
* champ point centre
  YY  = 'ALEA' 'BANDES_TOURNANTES' HYCEN
               'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
 
* champ point sommet (on fournit en test le maillage de 'QUA4')
  YY  = 'ALEA' 'BANDES_TOURNANTES' CARRE
               'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
 
* champ par élément sommet
  YY  = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'NOEUD'
               'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee  ;
 
* champ par élément centre
  YY  = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'GRAVITE'
               'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
 
* sans le mot-clef 'MOYENNE'
  YY  = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'GRAVITE'
               'EXPO' 'SIGMA' sigmay 'LAMBDA' portee  ;
 
* avec un modèle mécanique
  modmec = 'MODELE' CARRE 'MECANIQUE' ;
  YY  = 'ALEA' 'BANDES_TOURNANTES' modmec 'MASSE'
               'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
 
  YY  = 'ALEA' 'BANDES_TOURNANTES' modmec 'STRESSES'
               'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
 
* avec un double modèle dont l'un seulement convient
  mod2 = modhyb 'ET' modmec ;
  YY  = 'ALEA' 'BANDES_TOURNANTES' mod2 'RIGIDITE'
               'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee  ;
 
* anisotrope penché
  portee1 = 3*portee ;
  portee2 = .2*portee ;
  vec1    = 2. 1. ;
  vec2    = -1. 2.;
  YY  = 'ALEA' 'BANDES_TOURNANTES' CARRE
               'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM
               'LAMBDA1' portee1 vec1 'LAMBDA2' portee2 vec2 ;
 
*- en 3D
  'SI' vrai ;
    'OPTION' 'DIME' 3 'ELEM' 'CUB8' ;
    cube     = 'VOLUME' carre nz 'TRANSLATION' (0 0 lz)  ;
    domhyb3D = 'CHANGER' cube 'QUAF' ;
    modhyb3D = 'MODELE' DOMHYB3D 'DARCY' 'ISOTROPE' ;
 
*   plan
    YY  = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'GRAVITE'
                 'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee;
    'SI' GRAPH;
      'TITRE' 'Champ aleatoire plan en 3D, bandes tournantes';
      'TRACER' MODHYB YY ;
    'FINSI' ;
 
*   isotrope
    YY  = 'ALEA' 'BANDES_TOURNANTES' MODHYB3D 'GRAVITE'
                 'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee;
    'SI' GRAPH;
      'TITRE' 'Champ aleatoire isotrope 3D, bandes tournantes';
      'TRACER' MODHYB YY ;
    'FINSI' ;
 
*   unidirectionnel
    YY  = 'ALEA' 'BANDES_TOURNANTES' MODHYB3D 'GRAVITE'
                 'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM
                 'LAMBDA1' portee (0 0 1) ;
    'SI' GRAPH;
      'TITRE' 'Champ aleatoire monodirectionnel 3D, bandes tournantes';
      'TRACER' MODHYB YY ;
    'FINSI' ;
 
*   bidirectionnel
    YY  = 'ALEA' 'BANDES_TOURNANTES' MODHYB3D 'GRAVITE'
                 'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM
                 'LAMBDA1' portee 'LAMBDA2' portee ;
*    'SI' GRAPH;
      'TITRE' 'Champ aleatoire bidirectionnel 3D, bandes tournantes';
      'TRACER' MODHYB YY ;
*    'FINSI' ;
 
*   tridirectionnel stratifié
    YY  = 'ALEA' 'BANDES_TOURNANTES' MODHYB3D 'GRAVITE'
                 'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM
                'LAMBDA1' portee 'LAMBDA2' portee 'LAMBDA3' (.1*portee);
    'SI' GRAPH;
      'TITRE' 'Champ aleatoire stratifie 3D, bandes tournantes';
      'TRACER' MODHYB YY ;
    'FINSI' ;
 
    'OPTION' 'DIME' 2 ;
  'FINSI';
 
'FINSI' ;
 
* =====================================
* Test de la qualité de la génération :
* =====================================
* en 2D
 
'SI' veriresu ;
 
*-- on utilise un maillage plus fin et plus vaste :
  portee1 = portee ;
  portee2 = portee ;
  vec1    = (1. 0.) ;
  vec2    = (0. 1.) ;
  lx      = 20. * portee1;
  ly      = 20. * portee2;
  Np      = 5.;
  nx      = 'ENTIER' (lx / portee * Np + 0.0001) ;
  ny      = 'ENTIER' (ly / portee * Np + 0.0001) ;
 
  p1      = 0. 0. ;
  p2      = lx 0. ;
  p3      = lx ly ;
  p4      = 0. ly ;
  p1p2    = 'DROIT' p1 p2 nx ;
  p2p3    = 'DROIT' p2 p3 ny ;
  p3p4    = 'DROIT' p3 p4 nx ;
  p4p1    = 'DROIT' p4 p1 ny ;
  carre   = 'SURFACE' (p1p2 'ET' p2p3 'ET' p3p4 'ET' p4p1) ;
  DOMHYB  = 'CHANGER' CARRE  'QUAF';
  MODHYB  = 'MODELE' DOMHYB 'DARCY' 'ISOTROPE'    ;
 
* génération :
  YY  = 'ALEA' 'BANDES_TOURNANTES' CARRE
               'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
 
* Tracé :
  'SI' GRAPH;
    'TITRE' ('CHAINE' 'Champ aleatoire fin 2D, bandes tournantes, '
             ((lx / portee) @ARR 0) 'x' ((ly / portee) @ARR 0) ', l/'
             (Np @ARR 0) ) ;
    'TRACER' CARRE YY ('CONTOUR' CARRE) ;
  'FINSI' ;
 
*-- Calcul de la moyenne
  moyc = ('MINIMUM' ('RESULT' yy)) / ('NBNO' ('EXTR' yy 'MAIL')) ;
  mess 'moyenne : ' (moyc @ARR 3) ', objectif : ' (valym @ARR 3) ;
  err1 = 'ABS' ((moyc - valym) / ('MAXIMUM' ('PROG' moyc valym) 'ABS'));
 
*-- Calcul de l'écart-type
  sigc = ( ('MINIMUM' ('RESULT' ((yy - moyc) * (yy - moyc))))
         / ('NBEL' ('EXTR' yy 'MAIL')) ) ** .5 ;
  mess 'écart-type : ' (sigc @ARR 3) ', objectif : ' (sigmay @ARR 3) ;
  err2 = 'ABS' ((sigc - sigmay) / sigmay) ;
 
******************************************************
'SI' faux ;
* l'opérateur STAT n'est pas encore évolué.
 
*-- Calcul de la fonction de corrélation
*   ====================================
* on opère une transformation de l'espace réel qui ramène au cas
* isotrope de longueur de corrélation 1.
  'DEPLACER' carre 'AFFI' (1./portee1) (0. 0.) vec1 ;
  'DEPLACER' carre 'AFFI' (1./portee2) (0. 0.) vec2 ;
  portee0 = 1. ;
 
  lidist = 'PROG' 0. 'PAS' (portee0 / Np) lx  ;
  YYM    = 'MANU' 'CHPO' CARRE 1 'SCAL' moyc ;
  'SI' (portee1 'EGA' portee2) ;
    ev0    = 'STAT' 'COVA' YY YYM lidist 'REGU' ;
  'SINON' ;
    ev0    = 'STAT' 'COVA' YY YYM lidist ;
  'FINSI' ;
 
  evcov  = 'COULEUR'  ('EXTR' ev0 'COURBE' 1)  'ROUGE' ;
  evnb   = 'COULEUR' (('EXTR' ev0 'COURBE' 2)) 'VERT' ;
  evnb   = evnb * (varian / ('MAXIMUM' ('EXTR' evnb 'ORDO'))) ;
 
  liabs  = 'EXTR' evcov 'ABSC' ;
  licov  = 'EXTR' evcov 'ORDO'  ;
* on impose la valeur de la variance calculée avant (rétroactif)
  'REMPLACER' licov 1 (sigc * sigc) ;
 
  liordth = varian * ('EXP' (-1. * liabs / portee0)) ;
  evcovth = 'EVOL' 'DEFA' 'MANU' 'h' liabs 'Cyy' liordth ;
 
*-- Calcul de la longueur de corrélation
*   ====================================
* On cale la fonction evcov sur
*      y = varian * exp (- l * x)
* sur les 10 premières longueurs :
  litest = 'MASQUE' liabs 'EGINFE' (10. * portee0) ;
  ntest  = 'ENTIER' ('RESULT' litest) ;
  liind  = 'LECT' 1 'PAS' 1 ntest ;
  liabs2 = 'EXTR' liabs liind ;
  licov2 = 'EXTR' licov liind ;
 
* on part d'un jeu de couple yi = f(xi)
 
* utilisation de la procedure AJUSTE :
* on cherche a identifier les parametres a et b de la fonction
*                 y = a * (exp (- x / b))
*
* on part d'un jeu de couple yj = f(xj)
* p = ensemble des paramètres k mis-en-jeu dans des expressions non linéaires
* q = ensemble des paramètres i mis-en-jeu dans des expressions linéaires
*
 
* Table des couples (x_j,y_j) à atteindre :
* NOTA : l'indice 'F' suivant aurait dû s'appeller plutôt 'Y'.
  tab1 = 'TABLE' ;
  tab1.'X' = 'TABLE' ;
  tab1.'X' . 1 = liabs ;
  tab1.'F' = licov;
 
* nombre de paramètres non linéaires
  tab1.'K' = 1;
* nombre de paramètres linéaires (nb de jeux de valeurs à caler)
  tab1.'L' = 1;
 
* espace de discrétisation pour les paramètres
  pmin = 'PROG' 0.001;
  pmax = 'PROG' 10.;
  prc  = 'PROG' 0.0001;
 
  tab1.'MXTER'     = 150;
  tab1.'PMIN'      = pmin ;
  tab1.'PMAX'      = pmax ;
  tab1.'PRECISION' = prc ;
 
* Calcul de l'ajustement de la variance et de la longueur de
* corrélation :
  q p = AJUSTE tab1;
 
  varical = 'EXTR' q 1 ;
  lambcal = 'EXTR' p 1 ;
 
  ycal  = varical * ('EXP' (-1. / lambcal * liabs));
  evcal = 'EVOL' 'ROUG' 'MANU'  'X' liabs 'Y1' ycal;
 
  portcal = 1. / lambcal ;
  err3    = 'ABS' ((portcal - portee0) / portee0) ;
 
*-- tracé de la corrélation spatiale
  'SI' GRAPH;
    tabdess   = 'TABLE'  ;
    tabdess.1 = 'MOT' ' ';
    tabdess.2 = 'MOT' 'TIRR';
    tabdess.3 = 'MOT' 'TIRC';
    tabdess.4 = 'MOT' 'TIRL';
    tabdess.titre = table ;
    tabdess.titre.1 = 'CHAINE' 'Bantou';
    tabdess.titre.2 = 'CHAINE' 'Calage' ;
    tabdess.titre.3 = 'CHAINE' 'Theorique' ;
    tabdess.titre.4 = 'CHAINE' 'Nb couples' ;
    titdess  = 'Correlation du champ aleatoire, bandes tournantes' ;
    'DESSIN' (evcov 'ET' evcal 'ET' evcovth 'ET' evnb)
             'MIMA' 'LEGE' 'TITRE' titdess tabdess 'AXES'
             'XBOR' 0. ('MINIMUM' ('PROG' lx (10.*portee0))) ;
  'FINSI' ;
 
'SINON' ;
  err3 = 0;
'FINSI' ;
******************************************************
 
*-- Les valeurs calculées doivent être conformes aux valeurs théoriques
*   à 5% pour la moyenne et l'écart-type, 20% pour la longueur de corrélation.
  mess 'erreurs sur la moyenne, écart-type et covariance (%) = ';
  mess (100 * err1) (100 * err2) (100 * err3) ;
 
  'SI' ((err1 > 0.05) 'OU' (err2 > 0.05) 'OU' (err3 > 0.20));
     'ERREUR' 5;
  'SINON';
     'ERREUR' 0;
  'FINSI';
 
'FINSI' ;  'COMM' 'finsi veriresu' ;
 
'TEMPS' 'SGAC' 'IMPR' ;
 
'FIN' ;
 
 
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales