*
************************************************************************
* Cours MF307 - Tp3
* Diffusion d'un champ scalaire
* Solution stationnaire de l'équation de la chaleur
************************************************************************
*
'OPTI' 'ECHO' 0 ;
*
*-----------------------------------------------------------------------
'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' ;
*
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;
*-----------------------------------------------------------------------
*
*
* ==========
* ACQUISITION DES OPTIONS FOURNIES PAR L'UTILISATEUR
* ==========
*
*
*- 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 ;
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 ;
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 ;
'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 ;
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 authorisé 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' (NITER '<' 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 ;
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;
*
*- Modèle NAVIER_STOKES
*
MT = 'CHANGER' QUAF MT ;
'SI' ('EGA' IDEGR 1);
$MT = modele MT 'NAVIER_STOKES' 'LINE';
SINON ;
$MT = modele MT 'NAVIER_STOKES' 'QUAF';
'FINSI';
MT = 'DOMA' $mt maillage ;
CNT = 'CONT' mt ;
*
* ==========
* SOLUTION EXACTE
* ==========
*
*-> T(x,y) = 0.2 + sin(180*x) * e(-PI*y)
*
XX YY = 'COOR' MT ;
SINUS = 'SIN' ( XX * 180.) ;
EXPON = 'EXP' ( YY * (0. - PI)) ;
SOLEX = SINUS * EXPON + 0.2 ;
*
* ==========
* SOLUTION CASTEM
* ==========
*
* 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 ou EXIC)
* Les données numériques sont à ranger dans la table de
* soustype EQEX aux indices suivants :
* 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' ;
*
* Options pour le solveur (voir KRES)
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 ;
*
'TRAC' MT 'TITR' 'Maillage' ;
'TRAC' SOLEX MT CNT 'TITR' 'Solution exacte' ;
'TRAC' (RV . 'INCO' . 'TN') MT CNT 'TITR' 'Solution castem' ;
'TRAC' ER0 MT CNT 'TITR' 'Erreur absolue' ;
'TRAC' ER1 MT CNT 'TITR' 'Erreur relative' ;
* 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)';
'FINS' ;
* Tracé de la solution exacte en faux 3D
'OPTI' 'DIME' 3 'ISOV' 'SULI' ;
OEIL = -0.5 1.2 0.4 ;
MONTAGNE SOLEX MT 1.0 'Solution exacte' OEIL 3 ;
'OPTI' 'DIME' 2 ;
'FIN' ;
