Download equ_chaleur3D_VFconv.dgibi

Back to the list

* fichier : equ_chaleur3D_VFconv.dgibi
************************************************************************
* NOM         : equ_chaleur3D_VFconv.dgibi         
* ___
*
* DESCRIPTION : Solution stationnaire de l'équation de la chaleur  (3D)
* ___________  
*
*               GEOMETRIE : Un cube de côté 1 maillé avec 'VOLU'.
*
*               EQUATIONS :
*               ----------
*
*               EQUATIONS :
*               ----------
*
*               - Equations :
*
*                   div K GRAD T  - UT = 0
*
*                           |          
*                  avec K = |x+1  z    y|
*                           |z    y+1  x|
*                           |y    x    z+1|
*
*
*                  avec U = |1 + x + y + z|
*                           |1 + x + y + z|
*                           |1 + x + y + z|
*
*               - Conditions aux limites :
*
*                  conditions de Dirichlet, conditions de flux,
*                  conditions mixtes calculées à partir de la
*                  restriction de la solution exacte sur le bord
*                  exacte
*
*               - Solution exacte :
*
*                  T(x,y)= e(x+y+z)
*
*
*
*
* DISCRETISATION : une méthode de Volume Finis d'ordre 2 en espace. 
*
*
*
*       Le maillage est construit avec l'opérateur VOLU.
*
* Opérateurs utilisés :  PENT (option VF implicite)
*                        LAPN (option VF implicite)
*
*
* RESOLUTION :     - Solveur BiCGStab
* __________       - Préconditionneur ILU(0)
*
*
* TESTS EFFECTUES : Vérification de l'ordre 2 en espace de la méthode
* _______________   (utilisation d'une norme pseudo-L2) et de la 
*                   précision absolue sur le maillage le plus fin.
*
*
*
*
*
* LANGAGE     : GIBIANE-CASTEM 2000
* AUTEUR      : Christophe Le Potier (CEA/DEN/DM2S/SFME/LSET)
*               mél : clepotier@cea.fr
************************************************************************
* VERSION    : v1, 23/01/2010, version initiale
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
interact=  FAUX ;
complet =  FAUX ;
graph   = FAUX ;
*
'OPTION' 'DIME' 3 'ELEM' CUB8   ;
'OPTION' 'ISOV' 'SULI' ;
'OPTION' 'ECHO' 1 ;
nbisov = 15 ;
'SI' ('NON' interact) ;
  'OPTION' 'TRAC' 'PSC' ;
*  'OPTION' 'ISOV' 'LIGNE' ;
'FINSI' ;  
*******************************************************************
** Erreur Linfini entre deux Champoints.
*
*******************************************************************
DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ;
  AUX = (vitp1 - vit) ;
  err = MAXI (AUX) 'ABS' ;
  err = err / (MAXI(ABS(vit)));
FINPROC err ;
*
*
*******************************************************************       
** Erreur Pseudo L2 entre deux Champoints.
*
*******************************************************************       
DEBPROC CALCERR2 vitp1*'CHPOINT' vit*'CHPOINT' ;
  er0 = 'KOPS' vitp1 '-' vit ;
  er2 = 'KOPS' er0 '*' er0 ;
  er3 = 'KOPS' vit '*' vit ;
  err = 0.D0 ;
  erl = 0.D0;
  suppv  = 'EXTRAIRE' er2 'MAIL' ;
  compv  = 'EXTRAIRE' er2 'COMP' ;
  nbcomp = 'DIME' compv ;
  nptd   = 'NBNO' suppv ;
  'REPETER' numcomp nbcomp ;
  lacomp = 'EXTRAIRE' compv &numcomp ;
    'REPETER' numpo nptd ;
      lepoi = suppv 'POIN' &numpo ;
      err = err '+' ('EXTRAIRE' er2 lacomp lepoi) ;
      erl = erl '+' ('EXTRAIRE' er3 lacomp lepoi) ;
    'FIN' numpo ;
  'FIN' numcomp ;
  err = err '/' erl ;
  err = err ** 0.5D0 ;  
FINPROC err ;
*
*******************************************************************       
* Procédure paramétrée (raffinement)
* renvoyant l'erreur en norme L2 sur la température.
* On calcule une solution exacte de l'équation de Laplace
* (équation de la chaleur) ;
*
*nraff=1 ;
'DEBPROC' CALCUL nraff*'ENTIER' ;
*
* titre global pour les dessins
*
titglob = 'CHAINE' ' ; nraff=' nraff ;
*
*
*  Géométrie
*
pA = 0.   0.  0. ;
pB = 1.0   0.  0. ;
pC = 1.0   1.0  0. ;
pD = 0.    1.0  0. ;

pE = 0.   0.  1. ;
pF = 1.   0.  1. ;
pG = 1.   1.  1. ;
pH = 0.   1.  1. ;
*
* Paramètres de la discrétisation de base
*
'SI' complet ;
  nAB = 6 ;
  nBC = 6 ;
  nCD = 6 ;
  nDA = 6 ;
  nH  = 6 ;
'SINON' ;
  nAB = 2 '*' (nraff '+' 1) ;
  nBC = 2 '*' (nraff '+' 1) ;
  nCD = 2 '*' (nraff '+' 1) ;
  nDA = 2 '*' (nraff '+' 1) ;
  nH  = 2 '*' (nraff '+' 1) ;
'FINSI' ;  
*
*******************************************************************       
* Géométrie discrétisée
*******************************************************************       
*
'OPTION' 'ELEM' 'TRI3'  ;
bas    = 'DROIT' nAB pA pB ;
droite = 'DROIT' nBC pB pC ;
haut   = 'DROIT' nCD pC pD ;
gauche = 'DROIT' nDA pD pA ;
spourt = bas 'ET' droite 'ET' haut 'ET' gauche ;
sbas   = 'SURFACE' spourt 'PLAN' ;

'OPTION' 'DIME' 3 'ELEM' TET4   ;
mt = 'VOLUME' sbas 'TRAN' nH  pE ;
elim mt 0.001;
f1 f2 f3  = 'FACE' mt  ;
ss = (f1 'ET' f2 'ET' f3) ;
'OPTION' 'ELEM' 'TET4' ;
mt = 'VOLUME' ss ;
pourtour = 'ENVELOPPE' mt ;
*
* Eventuellement, on trace le résultat
*
*   list mt;
'SI' graph ;
   titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt)
                      ' NBELEM=' ('NBEL' mt) titglob ;
   'TRACER' mt 'NOEUD' 'TITRE' titgeo ;
'FINSI' ;
*******************************************************************       
*
* Modèle
*
*******************************************************************       
Mmt   = 'CHAN' mt   'QUAF' ;
'ELIM' Mmt 1.E-5 ;
MDNS  = 'NAVIER_STOKES' ;
$mt   = 'MODE' Mmt MDNS 'LINE' ;
mt    = 'DOMA' $mt 'MAILLAGE' ;
MESS 'NOMBRE DE MAILLES' (NBEL (MT)); 

*   list ('DOMA' $mt 'FACEP') ;
mtm = ('DOMA' $mt 'MAILLAGE') ;
mtc = 'COULEUR' ('DOMA' $mt 'CENTRE') 'JAUN' ;
mtf = 'COULEUR' ('DOMA' $mt 'FACE') 'BLEU' ;
*cmt = 'DOMA' $mt 'ENVELOPP' ;
cmt = pourtour;
Mcmt   = 'CHAN' cmt   'QUAF' ;
'ELIM' (Mmt 'ET' Mcmt) 1.E-5 ;
$cmt   = 'MODE' Mcmt MDNS  'QUAF';

* DEFINITION DES DIFFERENTS BORD DE LA FRONTIERE
bas1 = cmt 'ELEM' 'APPUYE' 'STRICTEMENT' f3 ;
gau1 = cmt 'ELEM' 'APPUYE' 'STRICTEMENT' (f1) ;
   titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt);
'SI' graph ;
trac bas1 'TITRE' 'bas1';
trac gau1 'TITRE' 'gau1';
'FINSI' ;
basgau = ( bas1 'ET' gau1);
'SI' graph ;
trac basgau 'TITRE' 'basgau';
'FINSI' ;
cvn = 'DIFF' cmt basgau ;
'SI' graph ;
trac cvn 'TITRE' 'cvn';
'FINSI' ;
basgau = ( bas1 'ET' gau1);
cvn = cmt 'ELEM' 'APPUYE' 'STRICTEMENT' f1 ;
bas1 = cmt 'ELEM' 'APPUYE' 'STRICTEMENT' f2 ;
basgau = ( bas1 'ET' cvn);
gau1 = diff cmt basgau  ;
Mbas  = 'CHAN' bas1  'QUAF' ;
Mgau  = 'CHAN' gau1  'QUAF' ;
Mcvn  = 'CHAN' cvn   'QUAF' ;
Mf1  = 'CHAN' f1   'QUAF' ;
'ELIM' (Mmt 'ET' Mbas 'ET' Mgau 'ET' Mcvn 'ET' Mf1) 1.E-5 ;
$bas  = 'MODE' Mbas MDNS 'QUAF' ;
$gau  = 'MODE' Mgau MDNS 'QUAF' ;
$cvn  = 'MODE' Mcvn MDNS 'QUAF' ;
$f1  = 'MODE' Mf1 MDNS  'QUAF';
cvn   = 'DOMA' $cvn 'CENTRE' ;
cvn1   = 'DOMA' $cvn 'MAILLAGE' ;
bas   = 'DOMA' $bas 'MAILLAGE' ;
gau   = 'DOMA' $gau 'CENTRE' ;
*

'SI' graph ;
pdir1ma = ('DOMA' $bas 'MAILLAGE');
pdir2ma = ('DOMA' $gau 'MAILLAGE');
cvnma = ('DOMA' $cvn 'MAILLAGE');
mai =  (cvnma 'COULEUR' 'VERT') et (pdir1ma 'COULEUR' ROUGE)
      et (pdir2ma 'COULEUR' BLEU);
 
 'TRACER'  (mai) 'TITRE' 'CONDITIONS AUX LIMITES'  ;
'FINSI' ;

*******************************************************************       
* Solution exacte aux centres
*******************************************************************       
xxc yyc zzc = 'COORDONNEE' ('DOMA' $mt 'CENTRE') ;
solexc = (exp ((+1.D0)*(xxc + yyc + zzc))) ; 
solexc = 'KCHT' $mt 'SCAL' 'CENTRE' solexc ;


*******************************************************************       
* terme source 
*******************************************************************       
sour = 0.D0 + (0.0*xxc);
sour = 'KCHT' $mt 'SCAL' 'CENTRE' sour ;

*******************************************************************
* calcul de la convection
*******************************************************************
CHYB2 = 'DOMA' $mt 'NORMALE'         ;
xxg yyg zzg = 'COORDONNEE' ('DOMA' $mt 'FACE') ;      
ux = (1.D0 + xxg + yyg + zzg);
uy = (1.D0 + xxg + yyg + zzg) ;
uz = (1.D0 + xxg + yyg + zzg) ;
ux = 'NOMC' 'UX' ux 'NATU' 'DISCRET' ;
uy = 'NOMC' 'UY' uy 'NATU' 'DISCRET' ;
uz = 'NOMC' 'UZ' uZ 'NATU' 'DISCRET' ;
qchal = (ux + uy + uz);
MOT1  = 'MOTS' 'UX' 'UY' 'UZ'                     ;
MOT2  = 'MOTS' 'UX' 'UY' 'UZ'                  ;
FLU1  = 'PSCA' CHYB2 qchal MOT1 MOT2         ;
FLUC = NOMC(FLU1) 'FLUX';
*******************************************************************
*   CONDITIONS DE NEUMANN
*******************************************************************
xxf yyf zzf = 'COORDONNEE' ('DOMA' $gau 'CENTRE') ;
gradtx = (+1.D0)*(1.D0 + (xxf + yyf + zzf)) * 
                 (exp ((+1.D0)*(xxf + yyf + zzf)))  ;
gradty = (+1.D0)*(1.D0 + (xxf+ yyf + zzf )) * 
                 (exp ((+1.D0)*(xxf + yyf + zzf))) ;
gradtz = (+1.D0)*(1.D0 + (xxf+ yyf + zzf )) * 
                  (exp ((+1.D0)*(xxf + yyf + zzf)))  ;
gradtx = 'NOMC' 'UX' gradtx 'NATU' 'DISCRET' ;
gradty = 'NOMC' 'UY' gradty 'NATU' 'DISCRET' ;
gradtz = 'NOMC' 'UZ' gradtz 'NATU' 'DISCRET' ;
CHYB2 = 'DOMA' $mt 'NORMALE'         ;
qchal = -1.D0 * (gradtx + gradty + gradtz) ;
qlim2 = qchal;
CHYB2R = 'REDU' CHYB2 gau    ;
MOT1  = 'MOTS' 'UX' 'UY' 'UZ'                 ;
MOT2  = 'MOTS' 'UX' 'UY' 'UZ'                ;
FLU1  = 'PSCA' CHYB2R qlim2 MOT1 MOT2         ;
qlim2 = NOMC(FLU1) 'FLUX';
*******************************************************************
* CONDITION DE DIRICHLET
*******************************************************************
xxlim yylim zzlim = 'COORDONNEE' ('DOMA' $bas  'CENTRE') ;
solim = (exp ((+1.D0)*(xxlim + yylim + zzlim))) ; 
solim = 'KCHT' $bas 'SCAL' 'CENTRE' solim ;

*******************************************************************
* CONDITIONS MIXTES
*******************************************************************
* DEFINITION DES PARAMETRES POUR LES CONDITIONS MIXTES
xxm yym zzm = 'COORDONNEE' ('DOMA' $cvn 'CENTRE') ;
gradmx = (+1.D0)*(1.D0 + (xxm + yym + zzm)) *
                 (exp ((+1.D0)*(xxm + yym + zzm)))  ;
gradmx = 'NOMC' 'UX' gradmx 'NATU' 'DISCRET' ;
gradmy = (+1.D0)*(1.D0 + (xxm + yym + zzm)) *
                 (exp ((+1.D0)*(xxm + yym + zzm)))  ;
gradmy = 'NOMC' 'UY' gradmy 'NATU' 'DISCRET' ;
gradmz = (+1.D0)*(1.D0 + (xxm + yym + zzm)) *
                 (exp ((+1.D0)*(xxm + yym + zzm)))  ;
gradmz = 'NOMC' 'UZ' gradmz 'NATU' 'DISCRET' ;
AUX = 1.D0 + (0.0D0*xxm);
XPAR1  = 'NOMC' 'PAR1' AUX 'NATU' 'DISCRET' ;
XPAR1 = 'REDU' XPAR1 cvn ;
qchal =  (-1.0) * (gradmx + gradmy + gradmz) ;
qaux =  'REDU' qchal cvn  ;
qlim = (1.D0)* (qaux);
CHYB2R = 'REDU' CHYB2 cvn    ;
MOT1  = 'MOTS' 'UX' 'UY' 'UZ'                     ;
MOT2  = 'MOTS' 'UX' 'UY' 'UZ'                     ;
FLU1  = 'PSCA' CHYB2R qlim MOT1 MOT2         ;
qlim = NOMC(FLU1) 'FLUX';
MOT1  = 'MOTS' 'UX' 'UY' 'UZ'                     ;
MOT2  = 'MOTS' 'UX' 'UY' 'UZ'                     ;
FLU1  = 'PSCA' CHYB2R qaux MOT1 MOT2         ;
XPAR2  = 'NOMC' FLU1 'PAR2';
XPAR2= (0.0D0*XPAR2);
* ON SIMULE UNE CONDITON DE NEUMANN VIA UNE CONDITION MIXTE  
*XPAR2 doit ETRE POSITIF
VAL = XPAR1 + XPAR2 + QLIM;

*******************************************************************
*
* Calcul du tenseur
*
*******************************************************************
K11 = 'KOPS' xxc  '+' 1.D0 ;
K22 = 'KOPS' yyc  '+' 1.D0 ;
K33 = 'KOPS' zzc  '+' 1.D0 ;
K21 =  zzc  ;
K31 =  (yyc)  ;
K32 =  (xxc)  ;
 
K11 = NOMC(K11) 'K11' 'NATU' 'DISCRET';
K22 = NOMC(K22) 'K22' 'NATU' 'DISCRET';
K33 = NOMC(K33) 'K33' 'NATU' 'DISCRET';
K21 = NOMC(K21) 'K21' 'NATU' 'DISCRET';
K31 = NOMC(K31) 'K31' 'NATU' 'DISCRET';
K32 = NOMC(K32) 'K32' 'NATU' 'DISCRET';

PERM = (K11 'ET' K22) ;
PERM = PERM 'ET' K33;
PERM = PERM 'ET' K21;
PERM = PERM 'ET' K31;
PERM = PERM 'ET' K32;

* Mise en place du calcul numérique
*
* équation de Laplace
*
*
* on utilise une méthode de Newton pour résoudre :
* - \Delta T = 0     (\Delta opérateur laplacien)
* - avec T_{\partial \Omega} donné (CL de Dirichlet)
* T_0 : estimation initiale de la solution
* On a T_1 = T_0 - {\Delta'}^{-1} (\Delta T_0)
*
* L'opérateur 'LAPN' 'VF' nous donne la matrice \Delta' et le
* résidu \Delta T_0.
* On n'inverse bien évidemment pas \Delta' mais on résoud le système:
* \Delta' IncT = \Delta T_0
* => IncT = {\Delta'}^{-1} (\Delta T_0)
*
* La méthode de Newton doit converger en un pas (on vérifie que le
* résidu (\Delta T_1) est nul après le premier pas.
*
*
t0 = solexc;
t0 = 'KCHT' $mt 'SCAL' 'CENTRE' 0.0 ;
menage;
gradt0 mchamt = 'PENT'  $mt
                       'FACE' 'MPFA' t0
                        'DISPDIF' PERM 'TIMP' 
                         solim 'QIMP' qlim2 'MIXT' val
                         'CENTERED' FLUC;
qaux =  'REDU' gradt0 cvn  ;
*list qaux;
jaco chpres dt = 'LAPN' 'VF' 'CLAUDEIS' 'IMPL'
$mt t0 gradt0 mchamt 'QIMP' qlim2 'MIXT' val 'TIMP' solim ;
mamat = 'KOPS' 'MULT' -1.0D0 jaco ;
matot =  mamat ;


AUX = ('EXCO' chpres 'RETN' 'RETN') - (NOMC(sour) 'RETN'); 
rv = 'EQEX' ;
rv . 'METHINV' . 'TYPINV'   = 3 ;
rv . 'METHINV' . 'PRECOND'  = 3 ;
rv . 'METHINV' . 'MATASS'   = matot ;
rv . 'METHINV' . 'MAPREC'   = matot ;
deltat = 'KRES' matot 'TYPI' (rv . 'METHINV')
                      'SMBR' AUX   
                      'IMPR' 0 ;
t1 = t0 '+' ('EXCO' deltat 'RETN' 'SCAL') ;
*******************************************************************
* Verification de la convergence
*******************************************************************
gradt1 mchamt    = 'PENT'  $mt  'FACE' 'MPFA' t1
                    'DISPDIF' PERM 'TIMP' solim 
                     'QIMP' qlim2 'MIXT' val 
                     'CENTERED' FLUC ;
jacbid chpres1 dt = 'LAPN' 'VF' 'CLAUDEISl' 'IMPL'
$mt t0 gradt1 mchamt 'QIMP' qlim2 'MIXT' val 'TIMP' solim ;
AUX = ('EXCO' chpres1 'RETN' 'RETN') - (NOMC(sour) 'RETN');
res = 'KRES' matot 'TYPI' (rv . 'METHINV')
                      'SMBR'AUX
                      'IMPR' 0 ;
mres1 = 'MAXIMUM' res  'ABS' ;
'MESSAGE' ('CHAINE' 'Maxi. mres1 =' mres1) ;
'SI' ('>' mres1 1.e-5) ;
    'MESSAGE' 'La méthode de Newton na pas converge en un pas.'
    'ERREUR' 5 ;
'FINSI' ;
*
* Résultats
*
'SI' graph ;
*
* solutions exactes
*
  tn = solexc ;
  chm_tn = 'KCHA' $mt 'CHAM' tn ;
  titt = 'CHAINE' 'Température exacte' titglob ;
  'TRACER' chm_tn $mt 'TITRE' titt ;
*
* graphe de convergence de la méthode itérative
*
    conver = (rv . 'METHINV' . 'CONVINV') ;
    dimcon = 'DIME' conver ;
    'SI' ('>' dimcon 1) ;
        labs = 'PROG' 1.D0 'PAS' 1.D0 dimcon ;
        lord = ('LOG' conver) '/' ('LOG' 10.D0) ;
        evtot = 'EVOL' 'MANU' 'ITER' labs 'RESID' lord ;
        titev = 'CHAINE' 'Historique de convergence' titglob ;
        'DESSIN' evtot 'TITR' titev 'LEGE' ;
    'FINSI' ;
*
* solutions calculées
*
  tn = t1 ;
  chm_tn = 'KCHA' $mt 'CHAM' tn ;
  titt = 'CHAINE' 'Température calculée' titglob ;
  'TRACER' chm_tn $mt  'TITRE' titt ;
*
* Calcul des erreurs par rapport à la solution analytique
*
* erreur
*
  erl = ABS(solexc '-' tn)/(maxi solexc);
  erlcham = 'KCHA' $mt 'CHAM' erl ;
  titt = 'CHAINE' 'erreur relative' ;
  'TRACER' (erlcham) $mt  'TITRE' titt ;
'FINSI' ;
echloc = '**' (('MESURE' mt) '/' ('NBEL' mt))
              (1.D0 '/' ('VALEUR' 'DIME')) ;
tn = t1 ;
errlit = CALCERR  tn solexc ;
errl2t = CALCERR2 tn solexc ;
'MESSAGE' '-------------------------------------------------' ;
'MESSAGE' ('CHAINE' 'Erreur en norme LINF : ' errlit) ;
'MESSAGE' ('CHAINE' 'Erreur en norme L2   : ' errl2t) ;
'MESSAGE' '-------------------------------------------------' ;
*'OPTION' 'DONN' 5 ;
'FINPROC' echloc errl2t ;
* Fin de la procédure de calcul                             *
*___________________________________________________________*
*-----------------------------------------------------------
*-----------------------------------------------------------
* Paramètres du calcul
*
*
* lraff : nb raffinement du maillage (à chaque fois, on découpe
*         les éléments en quatre).
'SI' complet ;
  lraff  = 'LECT' 0 PAS 1 2 ;
'SINON' ;
  lraff  = 'LECT' 1 PAS 1 4 ;
'FINSI' ;
*
*-----------------------------------------------------------*
* Boucles sur les différents paramètres du calcul           *
ordok = VRAI ;
precok = VRAI ;
* ordre théorique en espace de la méthode
ordtth = 2 ;
* erreur L2 max pour la solution (raffinement 1, complet=FAUX)
errtmax = 1.D-2 ;
*
*  Boucle sur les raffinements
*
    lh = 'PROG' ;
    lerl2t = 'PROG' ;
    'REPETER' iraff ('DIME' lraff)  ;
      raff = 'EXTRAIRE' lraff &iraff ;
      h errt = CALCUL raff ;
      lh     = 'ET' lh     ('PROG' h ) ;
      lerl2t = 'ET' lerl2t ('PROG' errt) ;
    'FIN' iraff ;
    lh     = ('LOG' lh)     '/' ('LOG' 10.D0) ;
    lerl2t = ('LOG' lerl2t) '/' ('LOG' 10.D0) ;
    tl2 = 'EVOL' 'MANU' 'Long. carac.' lh
                        'Err. tempér.' lerl2t ;
*
* Recherche des ordres de convergence
*
    cpt tlipol = @POMI tl2 1 'IDEM' ;
    ordt = cpt . 1 ;
    'MESSAGE' ('CHAINE' 'ordre température=' ordt) ;
*
* Tracés graphiques
*  
    'SI' graph ;
      tableg = 'TABLE' ;
      tableg . 'TITRE' = 'TABLE' ;
      tableg . 1 = 'MARQ CROI NOLI' ;
      tableg . 'TITRE' . 1 = 'Erreur L2 calculée' ;
      tableg . 2 = 'TIRR' ;                   
      tableg . 'TITRE' . 2 = 'Erreur L2 interpolée' ;
      titdest= 'CHAINE' 'Err L2 Tempér., ordre=' ordt ;
      'DESSIN' (tl2 'ET' tlipol) 'TITRE' titdest tableg ;
    'FINSI' ;
*
* Tests des ordres de convergence
* Tests des erreurs L2 sur le maillage le plus fin dans le
* cas complet=faux
*
     ordok = 'ET' ordok ('EGA' ('ENTIER' ('+' ordt 0.5)) ordtth) ;
     'SI' ('EGA' complet FAUX) ;
       precok = 'ET' precok ('<' errt errtmax) ;
     'FINSI' ;
     'SI' ('NON' ordok) ;
       'MESSAGE' 'On n_obtient pas un ordre de convergence correct' ;
       'ERREUR' 5 ;
     'FINSI' ;
     'SI' ('NON' precok) ;
       'MESSAGE' 'On n_obtient pas une précision correcte' ;
       'ERREUR' 5 ;
     'FINSI' ;
'SI' interact ;
   'OPTION' 'ECHO' 1 ;
   'OPTION' 'DONN' 5 ;
'FINSI' ;
'MESSAGE' 'Tout s_est bien passé' ;
'FIN' ;

 

 

 

© Cast3M 2003 - All rights reserved.
Disclaimer