Exemples de jeux de données Cast3M

* fichier :  equ_chaleur2D_VF.dgibi
************************************************************************
************************************************************************
***********************************************************************
* NOM         : equ_chaleur2D_VF.dgibi
* ___
*
* DESCRIPTION : Solution stationnaire de l'équation de la chaleur  (2D)
* ___________  
*
*               GEOMETRIE : Un carré
*               ----------
*
*               y
*                 ^   y=1
*                 |------------
*                 |           |
*                 |           |
*                 |           |
*                 |x = 0      |x=1 
*                 |           |    
*                 |           |    
*                O-----------------------------> x
*                     y=0                                   
*
*
*               EQUATIONS :
*               ----------
*
*               - Equations :
*
*                  mu Laplacien(T) = 0 avec mu = 1 ici
*                  
*               - Conditions aux limites :
*
*                  conditions de von Neumann  sur tout le bord except en
*                  (0;0); 
*                  on prend la restriction de la solution exacte au bord
*
*               - Solution exacte :
*
*                  T(x,y)= alpha x + beta y + gamma xy + delta
*
*
* DISCRETISATION : une méthode de Volume Finis d'ordre 2 en espace. 
* ______________   
*                 
*                 
*                 
*
*       Le maillage est construit avec l'opérateur SURF, il est perturbé
*       par un bruit blanc gaussien.
*
* Opérateurs utilisés : PENT (methode 'DIAMANT')
*                       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      : A. Beccantini
*               Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
************************************************************************
* VERSION    : v1, 15/04/02, version initiale
* HISTORIQUE : v1, 15/04/02, création
* 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 ;
 logxmgr = FAUX ;
*
 'OPTION' 'DIME' 2 'ELEM' QUA4 'MODE' 'PLAN' ;
 'OPTION' 'ISOV' 'SULI' ;
 'OPTION' 'ECHO' 0 ;
 nbisov = 15 ;
 'SI' ('NON' interact) ;
   'OPTION' 'TRAC' 'PS' ;
   'OPTION' 'ISOV' 'LIGNE' ;
 'FINSI' ;  
*
** Erreur Linfini entre deux Champoints.
*
 DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ;
    err = MAXI (vitp1 - vit) 'ABS' ;
 FINPROC err ;
*
** Erreur Pseudo L2 entre deux Champoints.
*
*
* L2 = \sqrt{\frac{1}{vol} .
*            \sum_{i} err_i^2 vol_i}
*
 DEBPROC CALCERR2 vitp1*'CHPOINT' vit*'CHPOINT' vol*'CHPOINT' ;
   er2 = vitp1 '-' vit ;
   compv  = 'EXTRAIRE' er2 'COMP' ;
   er2 = 'PSCAL' er2 er2 compv compv ;
   suppv  = 'EXTRAIRE' vol 'MAIL' ;
   chpun = 'MANUEL' 'CHPO' suppv 1 'SCAL' 1 ;
   voltot = 'XTY' chpun ('MOTS' 'SCAL') vol ('MOTS' 'SCAL') ;
   error = 'XTY' er2 ('MOTS' 'SCAL') vol ('MOTS' 'SCAL') ;
   error = error '/' voltot ;
   error = error '**' 0.5 ;
 FINPROC error ;
*
* Procédure paramétrée (raffinement)
* renvoyant l'erreur en norme L2 sur la température.
* On calcule une solution de l'équation de Laplace
* (équations de la chaleur) ;
*
 'DEBPROC' CALCUL nraff*'ENTIER' ;
*nraff=2 ;
*
* titre global pour les dessins
*
 titglob = 'CHAINE' ' ; nraff=' nraff ;
*
* Paramètres physiques
*
* 
 cv=  1.0 ;
 mu = 0.0 ;
 kappa = 1.0 ;
 rho = 1.0 ;
*
* Conditions aux limites et solution exacte:
*
* solex = (alpha '*' x) '+' (beta '*' y) '+' (gamma '*' x '*' y) '+'
*         delta
*
 alpha = (** 2. 0.5) ;
 beta  = (** 3. 0.5) ;
 gamma = (** 5. 0.5)  ;
 delta = (** 1.5 0.5) ;
*
*  Géométrie
*
 pA = 0.   0. ;
 pB = 1.   0. ;
 pC = 1.   1. ;
 pD = 0.   1. ;
*
* Paramètres de la discrétisation de base
*
 'SI' complet ;
    nAB = 8  ;
    nBC = 7  ;
    nCD = 5  ;
    nDA = 9  ;
 'SINON' ;
    nAB = 3 ;
    nBC = 4 ;
    nCD = 5 ;
    nDA = 6 ;
 'FINSI' ;
*
* Rotation et translation aditionnelle + bruit blanc
* + raffinement
*
 vtran = ((** 2 0.5) (* -1 (** 3 0.5))) ;
 orig  = (0.D0 0.D0) ;
 arot  = 11.3D0 ;
 lesdens = 'PROG' (1. '/' nAB) (1. '/' nCD) (1. '/' nBC) (1. '/' nDA) ;
 ampbruit = (0.15 * ('MINIMUM' lesdens)) ;
*
* Géométrie discrétisée
*
 bas    = 'DROIT' nAB pA pB ;
 droite = 'DROIT' nBC pB pC ;
 haut   = 'DROIT' nCD pC pD ;
 gauche = 'DROIT' nDA pD pA ;
 pourtour = bas 'ET' droite 'ET' haut 'ET' gauche ;
 mt = 'SURFACE' pourtour 'PLAN' ;
*
* On rajoute le bruit
*
 pmt = 'CHANGER' mt 'POI1' ;
 pcnt= 'CHANGER' pourtour 'POI1' ;
 inmt= 'DIFF' pmt pcnt;
 brbl1 = 'BRUIT' 'BLAN' 'GAUS' 0.D0 ampbruit inmt ;
 brbl2 = 'BRUIT' 'BLAN' 'GAUS' 0.D0 ampbruit inmt ;
 brbl = 'ET' ('NOMC' 'UX' brbl1) ('NOMC' 'UY' brbl2) ;
 'DEPLACER' mt 'PLUS' brbl;
*
* On raffine nraff fois 
*
 'SI' (nraff > 0) ;
    'REPETER' bcl nraff ;
       mt  = 'CHANGER' mt 'QUADRATIQUE' ;
       $mt = 'DOMA' mt 'MACRO' ;
       mt  = ('DOMA '$mt 'MAILLAGE') ;
       'MENAGE' ;
    'FIN' bcl ;
 'FINSI' ;
*
* Eventuellement, on trace le résultat
*
 'SI' graph ;
     titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt)
                      ' NBELEM=' ('NBEL' mt) titglob ;
    'TRACER' mt 'NOEUD' 'TITRE' titgeo ;
 'FINSI' ;
*
* Definition des cotés gauche et bas
*
 cnt = 'CONTOUR' mt ;
 gau = cnt 'ELEM' 'COMPRIS' PD PA 'COULEUR' 'VERT' ;
 bas  = cnt 'ELEM' 'COMPRIS' PA PB 'COULEUR' 'TURQ' ;
 'SI' graph ;
    'TRACER' (gau 'ET' bas) 'TITRE' 'Coté bas et coté gauche' ;
 'FINSI' ;
*
* Creation des modèles
*
* MODELS 
 $mt   = 'MODE' mt 'EULER'  ;
 $bas  = 'MODE' bas 'EULER'  ;
 $gau  = 'MODE' gau 'EULER'  ;
 $cnt  = 'MODE' cnt 'EULER'  ;
* MODELS 
 Tmt   = 'DOMA' $mt 'VF'  ;
 Tbas  = 'DOMA' $bas 'VF'  ;
 Tgau  = 'DOMA' $gau 'VF'  ;
 Tcnt  = 'DOMA' $cnt 'VF'  ;
* QUAF
 Mmt    = Tmt   . 'QUAF' ;
 Mgau   = Tgau  . 'QUAF' ;
 Mbas   = Tbas  . 'QUAF' ;
 Mcnt   = Tcnt  . 'QUAF' ;
 'ELIM' (Mmt 'ET' Mgau 'ET' Mbas 'ET' Mcnt)  1.E-5 ;
*
* Solution exacte bilinéaire (sur le sommet)
*
 xx yy = 'COORDONNEE' ('DOMA' $mt 'SOMMET') ;
 solex = (xx '*' alpha) '+' (yy '*' beta) '+'
         (xx '*' yy '*' gamma) '+' delta ;
*
* Solution exacte aux centres
*
 xxc yyc = 'COORDONNEE' ('DOMA' $mt 'CENTRE') ;
 solexc = (xxc '*' alpha) '+' (yyc '*' beta) '+'
         (xxc '*' yyc '*' gamma) '+' delta ;
*
* Gradient exact (aux faces)
*
 xxf yyf = 'COORDONNEE' ('DOMA' $mt 'FACE') ;
 gradtx = 'NOMC' 'UX' (alpha '+' (yyf '*' gamma)) 'NATU' 'DISCRET' ;
 gradty = 'NOMC' 'UY' (beta  '+' (xxf '*' gamma)) 'NATU' 'DISCRET' ;
* 
* 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  donné sur gauche (CL de Dirichlet)
*        grad(T) \cdot next donné ailleurs
*
* 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.
* 
*
 rho0 = 'MANUEL' 'CHPO' ('DOMA' $mt 'CENTRE') 1 'SCAL' 1.0 ;
 v0   = 'MANUEL' 'CHPO' ('DOMA' $mt 'CENTRE') 2 'UX' 0.0 'UY' 0.0 ;
 gradv0 mchamv = 'PENT' $mt 'FACE' 'DIAMANT' v0 ;
*
 t0 = 'MANU' 'CHPO' ('DOMA' $mt 'CENTRE') 1 'SCAL' 0.D0 ;
*
* We want to impose the temperature value close to pA
* So we cannot impose von Neumann boundary conditions
* close to PA
*
 cdir = 'MANUEL' 'POI1' pa ;
 pdir1 = 'MANUEL' 'POI1' (('DOMA' $bas 'CENTRE') 'POIN' 'PROC' pA) ;
 pdir2 = 'MANUEL' 'POI1' (('DOMA' $gau 'CENTRE') 'POIN' 'PROC' pA) ;
* cvn = set of FACE centers where we impose von Neumann boundary
* conditions
 cvn = 'DIFF' ('DOMA' $cnt 'CENTRE') pdir1 ;
 cvn = 'DIFF' cvn pdir2 ;
 'SI' graph ;
   'TRACER' (mt et (cvn 'COULEUR' 'ROUG')) 'TITRE' 'Von Neumann B.C.' ;
 'FINSI' ;  
 tlim = 'REDU' solex cdir ;
 qchal = -1.0 '*' kappa  '*' (gradtx 'ET' gradty) ;
 qlim = 'REDU' qchal cvn ;
 gradt0 mchamt = 'PENT'  $mt
                       'FACE' 'DIAMANT' t0
                       'CLIM' tlim ;
* Jacobian
* We use the 'NS' jacobian                        
 listinco = 'MOTS' 'RN' 'RUXN' 'RUYN' 'RETN' ;
 jaco chpres dt = 'LAPN' 'VF' 'PROPCOST' 'RESI' 'IMPL'
     $mt mu kappa cv rho0 v0 t0 gradv0 gradt0 mchamv mchamt
     'QIMP' qlim 'TIMP' tlim  listinco ;
 mamat = 'KOPS' 'MULT' -1.0D0 jaco ;
* Identity matrix for momentum and density contribution 
 mati = 'KOPS' 'MATIDE' ('MOTS' 'RN' 'RUXN' 'RUYN')
      ('DOMA' $mt 'CENTRE') 'MATRIK' ;
 matot = 'ET' mamat mati ;
*
 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' ('EXCO' chpres 'RETN' 'RETN')
                 'IMPR' 0 ;
 t1 = t0 '+' ('EXCO' deltat 'RETN' 'SCAL') ;
 gradt1 mchamt = 'PENT'  $mt
                       'FACE' 'DIAMANT' t1
                       'CLIM' tlim ;
 listinco = 'MOTS' 'RN' 'RUXN' 'RUYN' 'RETN' ;
* We check that the residuum is 0 if computed with
* gradt1 and t1
 jacbid chpres1 dt = 'LAPN' 'VF' 'PROPCOST' 'RESI' 'IMPL'
    $mt mu kappa cv rho0 v0 t1 gradv0 gradt1 mchamv mchamt
   'QIMP' qlim 'TIMP' tlim  listinco ;
    mres1 = 'MAXIMUM' ('EXCO' chpres1 'RETN') 'ABS' ;
 'MESSAGE' ('CHAINE' 'Maxi. chpres1 =' 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_tnex = 'KCHA' $mt 'CHAM' tn ;
  titt = 'CHAINE' 'Température exacte' titglob ;
  'TRACER' chm_tnex $mt nbisov 'TITRE' titt ;
*
* graphe de convergence de la méthode itérative
*
    conver = (rv . 'METHINV' . 'CONVINV') ;
    dimcon = 'DIME' conver ;
    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' ;
*
* solutions calculées
*
   tn = t1 ;
   chm_tn = 'KCHA' $mt 'CHAM' tn ;
   titt = 'CHAINE' 'Température calculée' titglob ;
   'TRACER' chm_tn $mt nbisov 'TITRE' titt ;
*
* erreur
*
   titt = 'CHAINE' 'Abs (Température calculée - Température exacte)'
         titglob ;
   'TRACER' ('ABS' (chm_tnex '-' chm_tn)) $mt nbisov 'TITRE' titt ; 
 'FINSI' ;
*
* Calcul des erreurs par rapport à la solution analytique
*
 vol = 'DOMA' $mt 'VOLUME' ;
 echloc = (('MESURE' mt) '/' ('NBEL' mt)) ** 0.5 ;
 tn = t1 ;
 errlit = CALCERR  tn solexc ;
 errl2t = CALCERR2 tn solexc vol ;
 '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 3 ;
 'SINON' ;
   lraff  = 'LECT' 0 PAS 1 3 ;
 'FINSI' ;
*
*-----------------------------------------------------------*
* Boucles sur les différents paramètres du calcul           *
 precok = VRAI ;
* ordre théorique en espace de la méthode
 ordtth = 2 ;
* erreur L2 max pour la solution (raffinement 2, complet=FAUX)
 errtmax = 1.D-3 ;
*
*  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 = lh '/' ('EXTRAIRE' lh 1) ;
     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 = ordt '>' (ordtth '-' 0.5) ;
      '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 ;
    'SI' logxmgr ;
*     Sortie pour xmgr   
       'OPTION' ECHO 0 ;
       'OPTION' 'IMPR' 'file.txt' ;
       ndim = 'DIME' lh ;
       'REPETER' BL1 ndim ;
         'MESSAGE' ('EXTRAIRE' lh &BL1) '  '  ('EXTRAIRE' lerl2t &BL1) ;
       'FIN' BL1 ;
       lh1 = 'EXTRAIRE' tlipol 'ABSC' ;
       lerr = 'EXTRAIRE' tlipol 'ORDO' ;
       ndim = 'DIME' lh1 ;
       'OPTION' 'IMPR' 'file2.txt' ;
       'REPETER' BL1 ndim ;
          'MESSAGE' ('EXTRAIRE' lh1 &BL1) '  '  ('EXTRAIRE' lerr &BL1) ;
       'FIN' BL1 ;
       'OPTION' 'IMPR' 'poubelle.txt' ;
       'FIN' ; 
    'FINSI' ;
 'FINSI' ;
 
 'MESSAGE' 'Tout s_est bien passé' ;
 'FIN' ;