Accueil | Cast3M

* fichier :  darcy4.dgibi
************************************************************************
************************************************************************
*
************************** CAS TEST : darcy4.dgibi ******************
*
 GRAPH = 'N'   ;
'SAUT' 'PAGE' ;
*
*---------------------------------------------------------------------
* TEST DARCY4
*               CALCUL DARCY ISOTROPE TRANSITOIRE
*                      AVEC TERME SOURCE
*    Résolution par une méthode d'éléments finis mixtes hybrides.
*
* Ce test permet de vérifier le bon fonctionnement des opérateurs
* utilisés afin de résoudre les équations de DARCY transitoire par
* une méthode d'éléments finis mixtes hybrides dans CASTEM2000.
*
* On effectue un calcul sur un domaine 2D maillé par des
* quadrangles réguliers. Le phénomène étudié est monodimensionnel.
* Les conditions aux limites sont de type trace de charge imposée
* à gauche du domaine, flux imposé à droite et flux nul en haut et
* en bas.
*
* La solution analytique en charge est un polynome de degré 2
* en x et en t, les variables x et t étant séparées.
* La vitesse est donc linéaire en x et indépendante du temps.
* Le terme source dépend de x et de t.
*
* Solution analytique :
*  H(x,t) = a * t**2 + b * t + c * x**2 + d * x + e
*  V(x,t) = - K * ( 2. * c * x + d )
*  S(x,t) = 2 * a * f * t + b * f - 2 * c * K
*
*  avec K  : Conductivité hydraulique .................    1.
*       a  : Coefficient ..............................   -1.
*       b  : Coefficient ..............................   20.
*       c  : Coefficient ..............................   -0.01
*       d  : Coefficient ..............................    2.
*       e  : Coefficient ..............................   10.
*       f  : Coefficient d'emmagasinnement ............    2.
*
* Remarque : Ne pas prendre t>40 car h devient négatif et on peut
*            planter le calcul d'erreur relative.
*
* Parametres de maillage :
*   L     : longueur du domaine .......................  100.
*   H     : hauteur du domaine ........................    1.
*   INUMX : Nombre d'éléments en x ....................   25
*   INUMY : Nombre d'éléments en y ....................    1 car 1D
*
*----------------------------------------------------------------
*
*------------------
* Options generales
*------------------
*
'OPTI' 'DIME' 2 'ELEM' 'QUA4' ;
'OPTI' 'ISOV' 'SURF'          ;
 'OPTI' 'ECHO' 0               ;
*
*
*=========
* MAILLAGE
*=========
*
*
*  d1------------d3
*  |              |      Post-traitement via les lignes b1b3 et c1c3
*  b1 c1------c3 b3  b1b3 : Extrémités de la droite joignant les faces
*  |              |  c1c3 : Extrémités de la droite joignant les centres
*  a1------------a3
*
*
*- Création des points supports des DROITES
*
L     = 100.D0                 ;
H     = 1.D0                   ;
X0    = 0.D0                   ;
X1    = L                      ;
Y0    = 0.D0                   ;
Y1    = Y0 + H                 ;
INUMX = 25                     ;
INUMY = 1                      ;
INUM1 = INUMX - 1              ;
DL    = L / INUMX              ;
DH    = H / INUMY              ;
Y01   = Y0 + Y1 * 0.5D0        ;
DX1   = X1 - X0 / INUMX / 2.D0 ;
XG    = X0 + DX1               ;
XD    = X1 - DX1               ;
*
A1 = X0  Y0  ;
A3 = X1  Y0  ;
D1 = X0  Y1  ;
D3 = X1  Y1  ;
B1 = X0  Y01 ;
B3 = X1  Y01 ;
C1 = XG  Y01 ;
C3 = XD  Y01 ;
*
*- Création des DROITES frontieres
*
DRBAS = A3 'DROI' INUMX A1 ;
DRGAU = A1 'DROI' INUMY D1 ;
DRHAU = D1 'DROI' INUMX D3 ;
DRDRO = D3 'DROI' INUMY A3 ;
DRMIF = B1 'DROI' INUMX B3 ;
DRMIC = C1 'DROI' INUM1 C3 ;
ELI0  = DH  / 10.D0        ;
*
*- Creation maillage GEOMETRIQUE
*
PTOT1 = 'DALL' DRBAS DRGAU DRHAU DRDRO ;
PTOT2 = 'ORIE' PTOT1                   ;
QFTOT = CHANGE PTOT2 QUAF ;
QFGAU = CHANGE DRGAU QUAF ;
QFDRO = CHANGE DRDRO QUAF ;
 ELIM ELI0 ( QFTOT ET QFGAU ET QFDRO ET DRMIF ET DRMIC) ;
*
 
*
*
*
*                                                    ----------------
*                                                    = MODELISATION =
*                                                    ----------------
MODHYB = MODE QFTOT 'DARCY' 'ISOTROPE' ;
MODGAU = MODE QFGAU 'DARCY' 'ISOTROPE' ;
MODDRO = MODE QFDRO 'DARCY' 'ISOTROPE' ;
CHYB1 = 'DOMA' MODHYB 'SURFACE'         ;
CHYB2 = 'DOMA' MODHYB 'NORMALE'         ;
CETOT  = DOMA MODHYB 'CENTRE' ;
CEDRO  = DOMA MODDRO 'CENTRE' ;
CEGAU  = DOMA MODGAU 'CENTRE' ;
*
*=====================================
* CONDITIONS INITIALES : CAS PERMANENT
*=====================================
*
*                                             -----------------------
*                                             = Solution analytique =
*                                             -----------------------
A  = -1.D0   ;
B  = 20.D0   ;
C  = -0.01D0 ;
D  =  2.D0   ;
E  = 10.D0   ;
F  =  2.D0   ;
VK =  1.D0   ;
*                                                    TRACE DE CHARGE
XX YY  = 'COOR' (DOMA MODHYB 'FACE'  )        ;
PANAF0 = C * XX * XX + ( D * XX ) + E                ;
PANAF0 = 'EXCO' PANAF0 'SCAL' 'TH'                   ;
PANAF0 = 'CHAN' 'ATTRIBUT' PANAF0 'NATURE' 'DISCRET' ;
*                                                             CHARGE
XXC YYC = 'COOR' CETOT                       ;
PANAC0  = C * XXC * XXC + ( D * XXC ) + E             ;
PANAC0  = 'EXCO' PANAC0 'SCAL' 'H'                    ;
PANAC0  = 'CHAN' 'ATTRIBUT' PANAC0 'NATURE' 'DISCRET' ;
*                                                            VITESSE
VXANA = -1.D0 * VK * ( 2.D0 * C * XXC + D )                         ;
VXANA = 'EXCO' VXANA 'SCAL' 'VX'                                    ;
VXANA = 'CHAN' VXANA 'ATTRIBUT' 'NATURE' 'DISCRET'                  ;
VYANA = 'MANU' 'CHPO' CETOT 1 'VY' 0.D0 'NATURE' 'DISCRET' ;
VANAC = VXANA 'ET' VYANA                                            ;
MOT1  = 'MOTS' 'VX' 'VY'                                            ;
*                                                    ----------------
*                                                    =  Tenseur K   =
*                                                    ----------------
MATI1 = 'MATE' MODHYB 'K' VK ;
*                                           -------------------------
*                                           = matrice Masse HYBride =
*                                           -------------------------
CND1A = 'MHYB' MODHYB  MATI1 ;
*                                           -------------------------
*                                           = MAtrice globale en TH =
*                                           -------------------------
HND1A = 'MATP' MODHYB  CND1A ;
*                                             -----------------------
*                                             = Second membre en TH =
*                                             -----------------------
AIRE   = DL * DH                                                      ;
SR0    = -2.D0 * VK * C * AIRE                                        ;
SRCVA0 = 'MANU' 'CHPO' CETOT 1 'SOUR' SR0 'NATURE' 'DISCRET' ;
SMSOUR = 'SMTP' MODHYB  CND1A SRCVA0                             ;
*                                                   -----------------
*                                                   =  TH imposée   =
*                                                   -----------------
BBGAU = 'BLOQ' CEGAU 'TH' ;
BBDRO = 'BLOQ' CEDRO 'TH' ;
EEGAU = 'DEPI' BBGAU PANAF0        ;
EEDRO = 'DEPI' BBDRO PANAF0        ;
*                                                   -----------------
*                                                   = Résolution TH =
*                                     *             -----------------
CCC1   = HND1A  'ET' BBDRO 'ET' BBGAU ;
FFF1   = SMSOUR 'ET' EEDRO 'ET' EEGAU ;
CHTER1 = 'RESO' CCC1 FFF1             ;
TH1    = 'EXCO' CHTER1 'TH' 'TH'      ;
*                                       -----------------------------
*                                       = Post-traitement Q, H et V =
*                                       -----------------------------
PCEN1  = 'HYBP' MODHYB  CND1A  TH1 SRCVA0       ;
QFACE1 = 'HDEB' MODHYB  CND1A  PCEN1 TH1   ;
VCENT1 = 'HVIT' MODHYB  QFACE1        ;
*
*
*
*===================================
* DONNEES SPECIFIQUES AU TRANSITOIRE
*===================================
*
*
*
TINI   = 0.D0                                  ;
TFINAL = 20.D0                                 ;
TSUP   = 1.5D0 * TFINAL                        ;
DELTAT = 1.0D0                                 ;
LICALC = 'PROG' 0.D0  'PAS' DELTAT TFINAL      ;
LISAUV = 'PROG' 0.D0  'PAS' 5.D0   TFINAL      ;
COEMM  = 'MANU' 'CHML' (DOMA MODHYB 'MAILLAGE') 'CK' F ;
*                                ------------------------------------
*                                = Nouvelles conditions aux limites =
*                                ------------------------------------
*
* Définition d'un listreel contenant les valeurs des temps de calcul
* Ce LISTREEL sert à définir les chargements variables en temps
*
TDEF   = PROG 0.D0 PAS 0.001D0 TSUP ;
*
* Valeur imposée à gauche = a*t*t + b*t + e
*
BLCGE  = BBGAU                                                ;
VALGA0 = 'MANU' 'CHPO' CEGAU 1 'TH' E                ;
VALGA0 = 'DEPI' BBGAU VALGA0                                  ;
VALI0  = 'CHAR' VALGA0
        ('EVOL' 'MANU' ('PROG' 0.D0 TSUP) ('PROG' 1.D0 1.D0)) ;
VALGA1 = 'MANU' 'CHPO' CEGAU 1 'TH' B                ;
VALGA1 = 'DEPI' BBGAU VALGA1                                  ;
VALI1  = 'CHAR' VALGA1
        ('EVOL' 'MANU' ('PROG' 0.D0 TSUP) ('PROG' 0.D0 TSUP)) ;
VALGA2 = 'MANU' 'CHPO' CEGAU 1 'TH' A                ;
VALGA2 = 'DEPI' BBGAU VALGA2                                  ;
TDEF2  = TDEF * TDEF                                          ;
VALI2  = 'CHAR' VALGA2 ('EVOL' 'MANU' TDEF TDEF2)             ;
TPGAU  = VALI0 'ET' VALI1 'ET' VALI2                          ;
*
* Flux imposé à droite = -K * ( 2. * c * x + d ) * dh
*
FLVAL  = 2.D0 * C * X1 + D * (-1.D0) * VK * DH                 ;
FLDRO  = 'MANU' 'CHPO' CEDRO 1 'FLUX' FLVAL           ;
FLXDRO = 'CHAR' FLDRO
        ('EVOL' 'MANU' ('PROG' 0.D0 TSUP) ('PROG' 1.D0 1.D0) ) ;
*
* Définition de la source moyenne par élément sur un pas de temps
*      s(t) = (2 * a * f * tmoy + b * f - 2 * c * K) * aire
*      avec aire = aire des éléments
*           tmoy = t(n+1) - deltat/2  = t(n) + deltat/2
*
SRCE01 = B * F -  ( 2.D0 * VK * C )                                    ;
SRCE02 = (-1.D0) * F * A * DELTAT                                      ;
SR0    = SRCE01 + SRCE02 * AIRE                                        ;
SRCV10 = 'MANU' 'CHPO' CETOT 1 'SOUR' SR0 'NATURE' 'DISCRET'  ;
SRCVA0 = SRCV10                                                        ;
NTDEF  = 'DIME' TDEF                                                   ;
LIST1  = 'PROG' NTDEF * 1.D0                                           ;
CHSRC0 = 'CHAR' SRCVA0 ('EVOL' 'MANU' TDEF LIST1)                      ;
PRSR1  = 2.D0 * A * F * AIRE * TDEF                                    ;
SRCVA1 = 'MANU' 'CHPO' CETOT 1 'SOUR' 1.D0 'NATURE' 'DISCRET' ;
CHSRC1 = 'CHAR' SRCVA1 ('EVOL' 'MANU' TDEF PRSR1 )                     ;
CHSRCE = CHSRC0 'ET' CHSRC1                                            ;
*
*                                         ---------------------------
*                                         = Table DARCY_TRANSITOIRE =
*                                         ---------------------------
TRANS2                 = 'TABLE' ;
TRANS2. 'TEMPS'        = 'TABLE' ;
TRANS2. 'TRACE_CHARGE' = 'TABLE' ;
TRANS2. 'CHARGE'       = 'TABLE' ;
TRANS2. 'FLUX'         = 'TABLE' ;
*
TRANS2.'SOUSTYPE'          = 'DARCY_TRANSITOIRE' ;
TRANS2.'MODELE'            = MODHYB              ;
TRANS2.'CARACTERISTIQUES'  = MATI1               ;
TRANS2.'EMMAGASINEMENT'    = COEMM               ;
TRANS2.'DOMAINE'           = DOMA MODHYB 'TABLE'   ;
*
TRANS2. 'TEMPS' . 0        = TINI                ;
TRANS2. 'TRACE_CHARGE' . 0 = TH1                 ;
TRANS2. 'CHARGE' . 0       = PCEN1               ;
TRANS2. 'FLUX' . 0         = QFACE1              ;
*
TRANS2.'BLOCAGE'           = BLCGE               ;
TRANS2.'TRACE_IMPOSE'      = TPGAU               ;
TRANS2.'FLUX_IMPOSE'       = FLXDRO              ;
TRANS2.'SOURCE'            = CHSRCE              ;
*
TRANS2.'THETA'             = 0.5D0               ;
TRANS2.'TEMPS_CALCULES'    = LICALC              ;
TRANS2.'TEMPS_SAUVES'      = LISAUV              ;
*
*=======================
* Resolution transitoire
*=======================
*
DARCYTRA TRANS2 ;
*
*
*=================
* POST-TRAITEMENT
*=================
*                                           --------------------------
*                                           = Calcul ERREUR ET TRACE =
*                                           --------------------------
*--------------------------------------------------------------------
*              Dans chaque cas on trace (si GRAPH<>N)
*--------------------------------------------------------------------
*      Les traces de charge le long de DRMIF (ligne b1b3)
*      L'erreur relative et l'erreur absolue en Trace de charge
*      La charge le long de DRMIC (ligne c1c3)
*      L'erreur relative et l'erreur absolue sur la charge
*      La Vitesse aux centres des elements (calculé et analytique)
*      L'erreur relative et l'erreur absolue sur la vitesse au centre
*--------------------------------------------------------------------
*
*
VDVD  = 'PSCA' VANAC VANAC MOT1 MOT1     ;
VCANA = 'VECT' VANAC 0.01D0 VX VY 'ROUG' ;
*--------------------------------------------------
*
*                  VISURESU
* Calcul des solutions analytiques et des erreurs
* Visualisation sur DRMIF et DRMIC des courbes y=f(x) à t fixé
*
*    ISOR1 : Liste des indices de la table TRANS2.'TEMPS'
*    NTSOR : Nombre de pas de temps sauvegardés
*    IOK   : Booleen indiquant si le critère de précision
*            demandé est vérifiée (Gestion des erreurs)
*    IRESU : Compteur interne permettant de récuperer les
*            indices de table via l'index ISOR1.
*
ISOR1 = 'INDEX' ( TRANS2 . 'TEMPS') ;
NTSOR = 'DIME' ISOR1                ;
IOK   = FAUX                        ;
IRESU = 0                           ;
*--------------------------------------------------
'REPETER' VISURESU NTSOR ;
*
  IRESU = IRESU + 1                ;
  INDI1 = ISOR1 . IRESU            ;
  TTRA  = TRANS2 . 'TEMPS' . INDI1 ;
*                                           Erreurs en trace de charge
  PANAF1 = A*TTRA*TTRA + (B*TTRA) + (C*XX*XX) + (D*XX) + E ;
  PANAD1 = 'REDU' PANAF1 DRMIF                             ;
  CALCF1 = TRANS2 . 'TRACE_CHARGE' . INDI1                 ;
  CALCF1 = 'EXCO' CALCF1 'TH' 'SCAL'                       ;
  CALCD1 = 'REDU' CALCF1 DRMIF                             ;
  ETPAB1 = CALCD1 - PANAD1                                 ;
  ETPAB1 = 'ABS' ETPAB1                                    ;
  ETPRE1 = ETPAB1 / PANAD1                                 ;
  ETPRE1 = 'ABS' ETPRE1                                    ;
*                                                    Erreurs en charge
  K1     = C * AIRE * AIRE /12.D0 + E                          ;
  PANAC1 = A*TTRA*TTRA + (B*TTRA) + (C*XXC*XXC) + (D*XXC) + K1 ;
  CALC1  = TRANS2 . 'CHARGE' . INDI1                           ;
  CALC1  = 'EXCO' CALC1 'H' 'SCAL'                             ;
  EPAB1  = CALC1 - PANAC1                                      ;
  EPAB1  = 'ABS' EPAB1                                         ;
  EPRE1  = EPAB1 / PANAC1                                      ;
  EPRE1  = 'ABS' EPRE1                                         ;
*                                                   Erreurs en vitesse
  EFL1  = TRANS2 . 'FLUX' . INDI1             ;
  EVI1  = 'HVIT' MODHYB  EFL1       ;
  VCAL1 = 'VECT' EVI1 0.01D0 'VX' 'VY' 'VERT' ;
  ERRV1 = EVI1 - VANAC                        ;
  EVAB1 = 'PSCA' ERRV1 ERRV1 MOT1 MOT1        ;
  EVRE1 = 'ABS' ( EVAB1 / VDVD )              ;
  EVAB1 = EVAB1 '**' 0.5D0                    ;
  EVRE1 = EVRE1 '**' 0.5D0                    ;
*                               Evolutions des solutions analytiques
  AV1 = 'EVOL' 'ROUG' 'CHPO' PANAD1 'SCAL' DRMIF ;
  AV3 = 'EVOL' 'ROUG' 'CHPO' PANAC1 'SCAL' DRMIC ;
*                                                             Tracés
  'SI' ('NEG' GRAPH 'N' ) ;
     LTITRE = 'TEXTE' 'darcy4 Th temps ' TTRA                   ;
     'TITR' LTITRE                                              ;
     EV1 = 'EVOL' 'VERT' 'CHPO' CALCD1 'SCAL' DRMIF             ;
     'DESS' (EV1 ET AV1)                                        ;
     LTITRE = 'TEXTE' 'darcy4 erreur relative Th temps ' TTRA   ;
     'TITR' LTITRE                                              ;
     EV2 = 'EVOL' 'VERT' 'CHPO' ETPRE1 'SCAL' DRMIF             ;
     'DESS' EV2                                                 ;
     LTITRE = 'TEXTE' 'darcy4 erreur absolue Th temps ' TTRA    ;
     'TITR' LTITRE                                              ;
     EV3 = 'EVOL' 'VERT' 'CHPO' ETPAB1 'SCAL' DRMIF             ;
     'DESS' EV3                                                 ;
     LTITRE = 'TEXTE' 'darcy4 h temps ' TTRA                    ;
     'TITR' LTITRE                                              ;
     EV4 = 'EVOL' 'VERT' 'CHPO' CALC1 'SCAL' DRMIC              ;
     'DESS' (EV4 ET AV3)                                        ;
     LTITRE = 'TEXTE' 'darcy4 erreur relative h temps' TTRA     ;
     'TITR' LTITRE                                              ;
     EV5 = 'EVOL' 'VERT' 'CHPO' EPRE1 'SCAL' DRMIC              ;
     'DESS' EV5                                                 ;
     LTITRE = 'TEXTE' 'darcy4 erreur absolue h temps' TTRA      ;
     'TITR' LTITRE                                              ;
     EV6 = 'EVOL' 'VERT' 'CHPO' EPAB1 'SCAL' DRMIC              ;
     'DESS' EV6                                                 ;
     LTITRE = 'TEXTE' 'darcy4 erreur relat V temps ' TTRA       ;
     'TITR' LTITRE                                              ;
     EV7 = 'EVOL' 'VERT' 'CHPO' EVRE1 'SCAL' DRMIC              ;
     'DESS' EV7                                                 ;
     LTITRE = 'TEXTE' 'darcy4 erreur absolue V temps ' TTRA     ;
     'TITR' LTITRE                                              ;
     EV8 = 'EVOL' 'VERT' 'CHPO' EVAB1 'SCAL' DRMIC              ;
     'DESS' EV8                                                 ;
  'FINSI' ;
*
*- Gestion ERREUR
*
  MAXTP1 = 'MAXI' ETPRE1 ;
  MAXP1  = 'MAXI' EPRE1  ;
  MAXVC1 = 'MAXI' EVRE1  ;
  LOGRES = IRESU 'EGA' 1 ;
  'SI' (LOGRES) ;
     'SAUT' 'PAGE'                                                    ;
     'SAUT' 2 'LIGNE'                                                 ;
     'MESS' '                       ERREURS RELATIVES                ';
     'SAUT' 1 'LIGNE'                                                 ;
     'MESS' '  temps              TH                H              Vc';
     'SAUT' 1 'LIGNE'                                                 ;
  'FINSI' ;
  'MESS' ttra  ' ' maxtp1  ' ' maxp1 ' ' maxvc1 ;
  EPS0  = 1.E-04                   ;
  LOG1  = MAXTP1 > EPS0            ;
  LOG2  = MAXP1  > EPS0            ;
  LOG3  = MAXVC1 > EPS0            ;
  L0    = LOG1 'OU' LOG2 'OU' LOG3 ;
  'SI' ( L0 ) ;
     IOK = VRAI ;
  'FINSI'     ;
'FIN' VISURESU ;
*
*********************************
*-     En cas de tracé uniquement
*********************************
*
'SI' ('NEG' GRAPH 'N' ) ;
*--------------------------------------------------
*
*                      CONSTAB
*             (préparation pour VISPTS)
* Transformation des champoints TH=f(t) en listreels
*
*    TABEVOL : Table contenant TH=f(t) sous forme de LISTREEL
*    IRESU   : Compteur interne permettant de récuperer les
*              indices de table via l'index ISOR1.
*
TABEVOL = 'TABLE' ;
IRESU   = 0       ;
*
*--------------------------------------------------
'REPETER' CONSTAB NTSOR ;
    IRESU           = IRESU + 1                              ;
    INDI1           = ISOR1 . IRESU                          ;
    TTRA            = TRANS2 . 'TEMPS' . INDI1               ;
    ETP1            = TRANS2 . 'TRACE_CHARGE' . INDI1        ;
    ETP1            = 'EXCO' ETP1 'TH' 'SCAL'                ;
    EV1             = 'EVOL' 'VERT' 'CHPO' ETP1 'SCAL' DRMIF ;
    ORD1            = 'EXTR' EV1 'ORDO'                      ;
    TABEVOL . IRESU = ORD1                                   ;
    'SI' ( IRESU 'EGA' 1)                 ;
       LITRAC = 'PROG' TTRA               ;
    'SINO'                                ;
       LITRAC = LITRAC 'ET' ('PROG' TTRA) ;
    'FINSI'                               ;
'FIN' CONSTAB ;
*--------------------------------------------------
*
*                      VISPTS
* Visualisation des courbes TH=f(t) pour un x fixé
*
*    NTTRA : Nombre de valeurs en temps à considérer
*    LISTX : Listreel des abscisses
*    NXTOT : Nombre de points sur la ligne DRMIF
*    IX    : Position dans LISTX du premier point à dépouiller
*    PASX  : Incrément pour le pas suivant à dépouiller
*
*
NTTRA = 'DIME' LITRAC     ;
LISTX = 'EXTR' EV1 'ABSC' ;
NXTOT = 'NBNO' DRMIF      ;
IX    = 5                 ;
PASX  = 10                ;
*
*--------------------------------------------------
'REPETER' VISTPS ;
*                                             Solution analytique
    VALX   = 'EXTR' LISTX IX                                 ;
    PPROG0 = 'PROG' NTTRA * (C*VALX*VALX + (D*VALX) + E )    ;
    PPROG1 = A * LITRAC * LITRAC + (B * LITRAC)              ;
    PEVOL1 = PPROG0 + PPROG1                                 ;
    EVOL1  = 'EVOL' 'ROUG' 'MANU' 'Temps' LITRAC 'Th' PEVOL1 ;
*                                               Solution calculée
    ICRB  = 0 ;
    'REPETER' CONSTR NTSOR                   ;
        ICRB  = ICRB + 1                     ;
        VAL1  = 'EXTR' TABEVOL.ICRB IX       ;
        'SI' ( ICRB 'EGA' 1 )                ;
           PROG1 = 'PROG' VAL1               ;
        'SINO'                               ;
           PROG1 = PROG1 'ET' ('PROG' VAL1 ) ;
        'FINSI'                              ;
    'FIN' CONSTR                             ;
    EVOL2  = 'EVOL' 'VERT' 'MANU' 'Temps' LITRAC 'Th' PROG1 ;
*                                                          Tracés
    LTITRE = 'TEXTE' 'darcy4 Th en x=' VALX                  ;
    'TITR' LTITRE                                            ;
    'DESS' ( EVOL1 ET EVOL2 )                                ;
    LTITRE = 'TEXTE' 'erreur relative Th x=' VALX            ;
    'TITR' LTITRE                                            ;
    ERRTP1 = PEVOL1 - PROG1 / PEVOL1                         ;
    EVOL3  = 'EVOL' 'VERT' 'MANU' 'Temps' LITRAC 'Th' ERRTP1 ;
    'DESS' EVOL3                                             ;
*                                         Test de sortie de boucle
    'DETR' PROG1 ;
    IX = IX + PASX      ;
    'SI' (IX > NXTOT )  ;
       'QUITTER' VISTPS ;
    'FINSI'             ;
'FIN' VISTPS            ;
*
'FINSI' ;
*
****************************
*            FIN des tracés
****************************
*
'SAUT' 4 'LIGNE' ;
*
'SI' ( IOK )   ;
   'ERRE' 5    ;
'SINO'         ;
   'ERRE' 0    ;
'FINSI'        ;
*
'FIN' ;