* 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' ;
© Cast3M 2003 - All rights reserved.
Disclaimer