Accueil | Cast3M

* fichier :  tp3.dgibi
************************************************************************
************************************************************************
*
************************************************************************
* Cours MF307 - Tp3
* Diffusion d'un champ scalaire
* Solution stationnaire de l'équation de la chaleur
************************************************************************
*
* Pour la non-regression, on vérifie que dans le cas d'une solution
* bilinéaire sur un maillage régulier de quadrangles la solution
* calculée est exacte à la précision machine près.
*
'OPTI' 'ECHO' 1 ;
GRAPH = faux ;
*
*-----------------------------------------------------------------------
'DEBP' CALCUL ;
*-----------------------------------------------------------------------
* Cette procédure, utilisée comme un opérateur, calcule la différence
* entre deux pas de temps toutes les 5 itérations.
* L'evolution de cette différence (erreur absolue) au cours du temps
* est conservée en vue d'un post-traitement.
*-----------------------------------------------------------------------
'ARGU' RVX*'TABLE' ;
*
RV  = RVX . 'EQEX' ;
TN  = RV  . 'INCO' . 'TN' ;
*
'SI' ( 'EXIS' RVX 'COMPT') ;
  RVX . 'COMPT' = RVX . 'COMPT' + 1 ;
'SINO' ;
  RVX . 'COMPT'         = 1 ;
  RV  . 'INCO' . 'IT'   = 'PROG' 1  ;
  RV  . 'INCO' . 'ER'   = 'PROG' 0. ;
  RV  . 'INCO' . 'TBIS' = 'COPI' RV . 'INCO' . 'TN' ;
*  'FINS' ;
'FINS' ;
*
DD = RVX . 'COMPT' ;
NN = DD / 5 ;
LO = (DD-(5*NN)) 'EGA' 0 ;
'SI' LO ;
   ERR    = (RV . 'INCO' . 'TN') - (RV . 'INCO' . 'TBIS') ;
   RV . 'INCO' . 'TBIS' = 'COPI' (RV . 'INCO' . 'TN') ;
   ERRINF = 'MAXI' ERR 'ABS' ;
   ELI    = ('LOG' (ERRINF + 1.0E-50))/('LOG' 10.) ;
  'MESS' 'ITERATION ' (RVX . 'COMPT') '   LOG10 ERREUR ' ELI ;
   IT = 'PROG' RVX . 'COMPT' ;
   ER = 'PROG' ELI ;
   RV . 'INCO' . 'IT' = (RV . 'INCO' . 'IT') 'ET' IT ;
   RV . 'INCO' . 'ER' = (RV . 'INCO' . 'ER') 'ET' ER ;
 
'FINS' ;
*-----------------------------------------------------------------------
mat1 chp1 = 'KOPS' 'MATRIK' ; 
'FINP' mat1 chp1;
*-----------------------------------------------------------------------
*
*
*===================================
* OPTIONS FOURNIES PAR L'UTILISATEUR
*===================================
* (les choix par menus sont à réactiver en interactif)
*
*
*- Choix du type de problème
*
'SAUT' 15 'LIGN' ;
'MESS' 'Choix du problème à traiter :' ;
'MESS' '    1 -> Problème linéaire'  ;
'MESS' '    2 -> Problème bilinéaire'  ;
'MESS' '    3 -> Sinusoide amortie' ;
**** 'OBTE' ISOL0 ;
ISOL0 = 2 ;
TEST1 = (ISOL0 '>' 3) 'OU' (ISOL0 '<' 1) ;
'SI' TEST1 ;
     'MESS' 'Erreur lors de la saisie du choix.';
     'FIN';
'FINS' ;
*
*- Choix du type d'éléments
*
'SAUT' 15 'LIGN' ;
'MESS' 'Choix des elements maillant le domaine :' ;
'MESS' '    1 -> Triangles'  ;
'MESS' '    2 -> Rectangles' ;
'MESS' '    3 -> Quadrilateres quelconques' ;
**** 'OBTE' CHOIX ;
CHOIX = 2 ;
TEST1 = (CHOIX '>' 3) 'OU' (CHOIX '<' 1) ;
'SI' TEST1 ;
     'MESS' 'Erreur lors de la saisie du choix.';
     'FIN';
'FINS' ;
*
*- Choix du type d'interpolation
*
'SAUT' 15 'LIGN' ;
'MESS' 'Choix des degrés des éléments finis :' ;
'MESS' '    1 -> Degré 1'  ;
'MESS' '    2 -> Degré 2' ;
**** 'OBTE' IDEGR ;
IDEGR = 1 ;
TEST1 = (IDEGR '>' 2) 'OU' (IDEGR '<' 1) ;
'SI' TEST1 ;
     'MESS' 'Erreur lors de la saisie du choix.';
     'FIN';
'FINS' ;
*
*- Choix du nombre d'éléments
*
'SAUT' 15 'LIGN' ;
'MESS' 'Nombre d elements par cote :';
**** 'OBTE' NX ;
NX = 10 ;
'SI' (NX '<' 1) ;
     'MESS' 'Erreur lors de la saisie du choix.';
     'FIN';
'FINS' ;
NY = NX ;
*
*- Choix de la méthode
*
'SAUT' 15 'LIGN' ;
'MESS' 'Type de methode de resolution :' ;
'MESS' '    1 -> Instationnaire explicite'  ;
'MESS' '    2 -> Instationnaire implicite' ;
'MESS' '    3 -> Stationnaire' ;
**** 'OBTE' IMETH ;
IMETH = 3 ;
TEST1 = (IMETH '>' 3) 'OU' (IMETH '<' 1) ;
'SI' TEST1 ;
     'MESS' 'Erreur lors de la saisie du choix.';
     'FIN';
'FINS' ;
'SI' (('EGA' IMETH 1) 'ET' ('EGA' IDEGR 2));
     'MESS' 'Seul degré 1 autorisé avec méthode explicite.';
     'FIN';
'FINS' ;
*
*- Nombre de pas de temps
*
'SI' ('NEG' IMETH 3) ;
     'SAUT' 15 'LIGN' ;
     'MESS' 'Nombre de pas de temps :';
     'OBTE' NITER ;
     'SI' (NX '<' 1) ;
         'MESS' 'Erreur lors de la saisie du choix.';
         'FIN';
     'FINS' ;
'FINS' ;
*
*- Choix de la méthode d'inversion
*
'SAUT' 15 'LIGN' ;
'MESS' 'Type de resolution matricielle :' ;
'MESS' '    1 -> Méthode directe (Crout)'  ;
'MESS' '    2 -> Gradient conjugué' ;
'MESS' '    3 -> Bi gradient conjugué' ;
**** 'OBTE' ININV ;
ININV = 1 ;
TEST1 = (ININV '>' 3) 'OU' (ININV '<' 1) ;
'SI' TEST1 ;
     'MESS' 'Erreur lors de la saisie du choix.';
     'FIN';
'FINS' ;
'SAUT' 15 'LIGN' ;
*
*
*=========
* MAILLAGE
*=========
*
*
*- Initialisation des options
*
'SI' (CHOIX 'EGA' 1) ;
  'OPTI' 'DIME' 2 'ELEM' 'TRI3' ;
'SINO' ;
  'OPTI' 'DIME' 2 'ELEM' 'QUA4' ;
FINSI;
*
*- Definition des sommets de la plaque
*
LX = 1.D0 ;
LY = 1.D0 ;
XMIN = 0.D0 ;
YMIN = 0.D0 ;
XMAX = XMIN + LX ;
YMAX = YMIN + LY ;
*
A1 = XMIN YMIN ;
A2 = XMAX YMIN ;
A3 = XMAX YMAX ;
A4 = XMIN YMAX ;
*
*- Definition des aretes
*
'SI' (CHOIX 'EGA' 3) ;
  R  = 10. ** (1. / (NX - 1));
  DA = 1. * ((R - 1.) / ((R ** NX) - 1));
  DB = 10. * DA ;
  FBAS = 'DROI' A1 A2 'DINI' DA 'DFIN' DB ;
  FDRO = 'DROI' A2 A3 'DINI' DB 'DFIN' DA ;
  FHAU = 'DROI' A3 A4 'DINI' DA 'DFIN' DB ;
  FGAU = 'DROI' A4 A1 'DINI' DB 'DFIN' DA ;
'SINO' ;
  FBAS = A1 'DROI' NX A2 ;
  FDRO = A2 'DROI' NY A3 ;
  FHAU = A3 'DROI' NX A4 ;
  FGAU = A4 'DROI' NY A1 ;
'FINS' ;
*
*- Maillage de la plaque
*
SI (CHOIX EGA 1);
  CNT = FBAS ET FDRO ET FHAU ET FGAU ;
  MT = CNT 'SURF' 'PLAN' ;
SINON;
  MT = 'DALLER' FBAS FDRO FHAU FGAU 'PLAN' ;
FINSI;
*
*- MODELE NAVIER_STOKES
*
MT = 'CHANGER' QUAF mt;
'SI' ('EGA' IDEGR 1); 
   $MT = 'MODE' MT 'NAVIER_STOKES' 'LINE';
SINON ;
   $MT = 'MODE' MT 'NAVIER_STOKES' 'QUAF';
'FINSI';
MT  = 'DOMA' $MT 'MAILLAGE' ;
CNT = 'CONTOUR' MT;
*
*
*================
* SOLUTION EXACTE
*================
*
*
XX YY = 'COOR' MT ;
'SI' (ISOL0 '<' 3) ;
*
*-> T(x,y) = 11xy + 7x + 5y + 3 + NORMALISATION
*
   'SI' (ISOL0 'EGA' 1) ;
      AXY = 0. ;
   'SINO' ;
      AXY = 11. ;
   'FINSI' ;
   SOLEX = AXY*XX*YY + (7.*XX) + (5.*YY) + 3. ;
   SOLEX = SOLEX / (MAXI SOLEX) ;
'SINO' ;
*
*-> T(x,y) = 0.2 + sin(180*x) * e(-PI*y)
*
   SINUS = 'SIN' ( XX * 180.) ;
   EXPON = 'EXP' ( YY * (0. - PI)) ;
   SOLEX = SINUS * EXPON + 0.2 ;
'FINSI' ;
*
*
*==================
* SOLUTION CALCULEE
*==================
*
*
* La structure du bloc de chaque méthode est identique :
*    -> description de la modélisation physique et numérique (EQEX)
*    -> conditions aux limites (EQEX, option CLIM)
*    -> initialisations des champs (table INCO)
*    -> initialisations des données numériques (table RV)
*    -> resolution (EXEC)
* Les données numériques suivantes permettent de piloter le calcul :
*    NITER : Nombre d'iterations internes (cas non-lineaire)
*    OMEGA : Coefficient de relaxation (cas non-linéaire)
*    IMPR  : Flag d'impression (1 pour oui, 0 pour non)
*    ITMA  : Nombre de pas de temps
*    EPS   : Critere de convergence pour les iterations internes
*    DT    : Pas de temps
*
*------------------------------------------------
* Méthode de résolution instationnaire explicite
* Utilisation de l'opérateur TSCAL
*------------------------------------------------
*
'SI' (IMETH 'EGA' 1) ;
   RV = 'EQEX' $MT 'ITMA' NITER 'ALFA' 0.9
              'ZONE' $MT 'OPER' CALCUL
              'OPTI' 'EFM1' 'CENTREE' 'EXPL'
              'ZONE' $MT 'OPER' 'TSCA' 1. 'UN' 0.0 'INCO' 'TN'
              'OPTI' 'EFM1' 'CENTREE' 'EXPL'
              'ZONE' $MT 'OPER' 'DFDT' 1. 'TN' 'DELTAT' 'INCO' 'TN'
              ;
 
   RV = 'EQEX' RV 'CLIM' CNT 'TN' 'TIMP' ('REDU' SOLEX CNT) ;
   RV . 'INCO'        = 'TABLE' 'INCO' ;
   RV . 'INCO' . 'TN' = 'KCHT' $MT 'SCAL' 'SOMMET' 0. ;
   RV . 'INCO' . 'UN' = 'KCHT' $MT 'VECT' 'SOMMET' (0. 0.) ;
'FINS' ;
*
*------------------------------------------------
* Méthode de résolution instationnaire implicite
* Utilisation des opérateurs DFDT et LAPN
*------------------------------------------------
*
'SI' (IMETH 'EGA' 2) ;
   RV = 'EQEX' $MT 'OPTI' 'EF' 'IMPL' 'CENTREE'
*        'ZONE' $MT 'OPER' 'DFDT' 1. 'TNM' 'DT' (0. 0.) 0. 'INCO' 'TN'
        'ZONE' $MT 'OPER' 'DFDT' 1. 'TNM' 'DT' 'INCO' 'TN'
        'ZONE' $MT 'OPER' 'LAPN' 1.            'INCO' 'TN'
        'ZONE' $MT 'OPER' CALCUL
   ;
   RV = 'EQEX' RV 'CLIM' CNT 'TN' 'TIMP' ('REDU' SOLEX CNT) ;
   RV . 'INCO'         = 'TABLE' 'INCO' ;
   RV . 'INCO' . 'TN'  = 'KCHT' $MT 'SCAL' 'SOMMET' 0. ;
   RV . 'INCO' . 'TNM' = 'KCHT' $MT 'SCAL' 'SOMMET' 0. ;
   RV . 'NITER'        = 1      ;
   RV . 'OMEGA'        = 1.D0   ;
   RV . 'IMPR'         = 1      ;
   RV . 'ITMA'         = NITER  ;
   RV . 'EPS'          = 1.D-15 ;
   RV . 'DT'           = 1.D-1  ;
   RV . 'INCO' . 'DT'  = RV . 'DT' ;
'FINS' ;
*
*------------------------------------------------
* Méthode de résolution stationnaire (implicite)
* Utilisation de l'opérateur LAPN
*------------------------------------------------
*
'SI' (IMETH 'EGA' 3) ;
   RV = 'EQEX' $MT 'OPTI' 'EF' 'IMPL'
        'ZONE' $MT 'OPER' 'LAPN' 1.D0 'INCO' 'TN' ;
   RV = 'EQEX' RV 'CLIM' CNT 'TN' 'TIMP' ('REDU' SOLEX CNT) ;
   RV . 'INCO'        = 'TABLE' 'INCO' ;
   RV . 'INCO' .'TN'  = 'KCHT' $MT 'SCAL' 'SOMMET' 0.D0 ;
   RV . 'NITER'       = 2       ;
   RV . 'OMEGA'       = 1.D0    ;
   RV . 'IMPR'        = 1       ;
   RV . 'ITMA'        = 1       ;
   RV . 'EPS'         = 1.D-15  ;
   RV . 'DT'          = 1.D3    ;
   RV . 'INCO' . 'DT' = RV . 'DT' ;
'FINS' ;
*
*
rv. 'METHINV' . 'TYPINV'  = ININV ;
rv. 'METHINV' . 'IMPINV'  = 2 ;
rv. 'METHINV' . 'NITMAX'  = 1500 ;
rv. 'METHINV' . 'PRECOND' = 3 ;
rv. 'METHINV' . 'RESID'   = 1.e-14 ;
rv. 'METHINV' . 'FCPRECT' = 1 ;
rv. 'METHINV' . 'FCPRECI' = 1 ;
rv. 'METHINV' . 'XINIT'   = RV.INCO.'TN';
rv. 'METHINV' . 'TYRENU' = 'SLOANE' ;
*
EXEC rv;
*
*
*==================================
* POST TRAITEMENT ET VISUALISATIONS
*==================================
*
* ER0 : Erreur absolue
* ER1 : Erreur relative
* ER2 : Erreur absolue au carrée
* ERR : Erreur L2 discrète
*
ER0 = 'ABS' (RV.INCO.'TN' - SOLEX) ;
ER1 = ER0 '/' SOLEX ;
ER2 = ER0 '*' ER0 ;
ERR = ('MAXI' ('RESU' ER2)) / ('NBEL' mt) ** 0.5 ;
*
*------------------------> Début des tracés
'SI' GRAPH ;
  'TRAC' MT 'TITR' 'Maillage' 'NCLK' ;
  'TRAC' SOLEX MT CNT 'TITR' 'Solution exacte' 'NCLK' ;
  'TRAC' (RV . 'INCO' . 'TN') MT CNT 'TITR' 'Solution castem' 'NCLK' ;
  'TRAC' ER0 MT CNT 'TITR' 'Erreur absolue' 'NCLK' ;
  'TRAC' ER1 MT CNT 'TITR' 'Erreur relative' 'NCLK' ;
*                        Evolution de l'erreur au cours des itérations
  'SI' ('NEG' IMETH 3) ;
     EVO1 = 'EVOL' 'MANU' 'n' RV.INCO.'IT' 'Log Epsi' RV.INCO.'ER';
     'DESS' EVO1 'TITRE' 'Histoire de l erreur absolue'
                 'TITX' 'Iterations' 'TITY' 'Log(Erreur)';
  'FINSI' ;
*                               Tracé de la solution exacte en faux 3D
  'OPTI' 'ISOV' 'SULI' ;
  'SI' ('<' ISOL0 3) ;
     OEIL = 'PROG' -0.5 -1.2 2. ;
  'SINON' ;
     OEIL = 'PROG' -0.5 1.2 0.4 ;
  'FINSI' ;
  MONTAGNE  SOLEX MT 1.0 OEIL 'TITRE' 'Solution exacte' 'CACHE' ;
'FINSI' ;
*------------------------> Fin des tracés
*
*
*=======================
* TEST DE NON-REGRESSION
*=======================
*
TEST1 = ERR < 1E-15 ;
si test1 ;
   erre 0 ;
sinon ;
   erre 5 ;
finsi ;
'FIN' ;