Download equ_chaleur2D_VFcyl.dgibi
* fichier : equ_chaleur2D_VFcyl.dgibi
***********************************************************************
* NOM : equ_chaleur2D_VFcyl.dgibi
* ___
*
* DESCRIPTION : Solution stationnaire de l'équation de la chaleur (3D)
* dans un doamine cylindrique
* ___________
*
* GEOMETRIE : 1 < r < 2
* ----------
*
*
* EQUATIONS :
* ----------
*
* - Equations :
*
* mu Laplacien(T) = 0 avec mu = 1 ici
*
* - Conditions aux limites :
*
* T = 0 @ r = 1
* DT/dr = 1 @ r = 2
*
* - Solution exacte :
*
* T(r)= 2 ln (r)
*
*
* DISCRETISATION : une méthode de Volume Finis d'ordre 2 en espace.
* ______________
*
*
* Opérateurs utilisés : PENT ('DIAMAN2')
* 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, 25/02/03, 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 ;
TYPEL = 'QUA4' ;
'OPTION' 'DIME' 2 'ELEM' TYPEL 'MODE' 'PLAN' ;
'OPTION' 'ISOV' 'SULI' ;
'OPTION' 'ECHO' 1 ;
nbisov = 15 ;
'SI' ('NON' interact) ;
'OPTION' 'TRAC' 'PS' ;
'OPTION' 'ISOV' 'LIGNE' ;
'OPTION' 'ECHO' 0 ;
'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 ;
*
* Géométrie
*
p00 = 0. 0. ;
p01 = 0. 0. ;
pA = 1. 0. ;
pB = 2. 0. ;
*
* Paramètres de la discrétisation de base
*
nAB = (nraff '+' 1) ;
*
* Géométrie discrétisée
*
angle = (180 '/' PI) * (1. '/' nab) ;
xPA1 = ('COORDONNEE' 1 PA) '*' ('COS' ANGLE) ;
yPA1 = ('COORDONNEE' 1 PA) '*' ('SIN' ANGLE) ;
PA1 = xpa1 ypa1;
xPB1 = ('COORDONNEE' 1 PB) '*' ('COS' ANGLE) ;
yPB1 = ('COORDONNEE' 1 PB) '*' ('SIN' ANGLE) ;
PB1 = xpb1 ypb1 ;
facint = PA1 'DROIT' 1 PA ;
facext = PB1 'DROIT' 1 PB ;
mt = facint 'REGL' nab facext ;
cnt = 'CONT' mt ;
faccon = 'DIFF' cnt (facint 'ET' facext) ;
*
* Eventuellement, on trace le résultat
*
'SI' graph ;
'TRACER' faccon ;
titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt)
' NBELEM=' ('NBEL' mt) titglob ;
'TRACER' mt 'TITRE' titgeo ;
'FINSI' ;
*
*
* Creation des modèles
*
* MODELS
MDNS = 'EULER' ;
$mt = 'MODE' mt 'EULER' ;
$cnt = 'MODE' cnt 'EULER' ;
$facint = 'MODE' facint 'EULER' ;
$facext = 'MODE' facext 'EULER' ;
$faccon = 'MODE' faccon 'EULER' ;
*
Tmt = 'DOMA' $mt 'VF' ;
Tcnt = 'DOMA' $cnt 'VF' ;
Tfacint = 'DOMA' $facint 'VF' ;
Tfacext = 'DOMA' $facext 'VF' ;
Tfaccon = 'DOMA' $faccon 'VF' ;
* QUAF
Mmt = Tmt . 'QUAF' ;
Mcnt = Tcnt . 'QUAF' ;
Mfacint = Tfacint . 'QUAF' ;
Mfacext = Tfacext . 'QUAF' ;
Mfaccon = Tfaccon . 'QUAF' ;
'ELIM' (Mmt 'ET' Mcnt 'ET' Mfacint 'ET' Mfacext 'ET' Mfaccon) 1.E-8 ;
* CHAMPOINT vide et MATRIK vide
chvid matvid = 'KOPS' 'MATRIK' ;
*
* Solution exacte bilinéaire (sur le centre du contour)
*
xx yy = 'COORDONNEE' ('DOMA' $cnt 'CENTRE') ;
*
rr2 = (xx '*' xx) '+' (yy '*' yy) ;
rr = rr2 '**' 0.5 ;
solex = 2 '*' ('LOG' rr) ;
*
* Solution exacte aux centres
*
xx yy = 'COORDONNEE' ('DOMA' $mt 'CENTRE') ;
rr2 = (xx '*' xx) '+' (yy '*' yy) ;
rr = rr2 '**' 0.5 ;
solexc = 2 '*' ('LOG' rr) ;
*
'SI' FAUX ;
*
* Gradient exact (aux faces)
*
* Ca ne marche pas (probleme 2)
* En effet j'ai 1.0002 à la place de 1
* Probleme de quadr. du cercle.
*
xx yy = 'COORDONNEE' ('DOMA' $mt 'FACE') ;
rr2 = (xx '*' xx) '+' (yy '*' yy) ;
gradtx = 'NOMC' 'UX' (2 '*' (xx '/' rr2)) 'NATU'
'DISCRET' ;
gradty = 'NOMC' 'UY' (2 '*' (yy '/' rr2)) 'NATU'
'DISCRET' ;
gradt = gradtx 'ET' gradty ;
* Test
nn = 'DOMA' $mt 'XXNORMAF' ;
aa = 'REDU' nn ('DOMA' $faccon 'CENTRE') ;
aa = 'MAXIMUM' ('PSCAL' aa gradt ('MOTS' 'UX' 'UY')
('MOTS' 'UX' 'UY')) 'ABS' ;
'SI' (aa > 1.0D-11) ;
'LISTE' aa ;
'MESSAGE' 'Probleme 1' ;
'ERREUR' 21 ;
'FINSI' ;
aa = 'REDU' nn ('DOMA' $facext 'CENTRE') ;
aa = 'MAXIMUM' ('PSCAL' aa gradt ('MOTS' 'UX' 'UY')
('MOTS' 'UX' 'UY')) 'ABS' ;
'SI' ((aa '-' 1.0) > 1.0D-11) ;
'LISTE' aa ;
'MESSAGE' 'Probleme 2' ;
'ERREUR' 21 ;
'FINSI' ;
gradtcon = 'REDU' gradt ('DOMA' $faccon 'CENTRE') ;
gradtext = 'REDU' gradt ('DOMA' $facext 'CENTRE') ;
'SINON' ;
*
* Gradient exact sur le bord
*
nn = 'DOMA' $mt 'XXNORMAF' ;
gradtcon = ('REDU' nn ('DOMA' $faccon 'CENTRE')) '*' 0.0 ;
gradtext = ('REDU' nn ('DOMA' $facext 'CENTRE')) '*' 1.0 ;
'FINSI' ;
*
* Mise en place du calcul numérique
*
* 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 ;
cdirv = 'MANUEL' 'POI1' ('POIN' 1 ('DOMA' $cnt 'CENTRE')) ;
cvnv = 'DIFF' ('DOMA' $cnt 'CENTRE') cdirv ;
v0li1 = 'MANUEL' 'CHPO' cdirv 2 'UX' 0.0 'UY' 0.0 ;
v0li2 = 'MANUEL' 'CHPO' cvnv 4 'P1DX' 0.0 'P1DY' 0.0
'P2DX' 0.0 'P2DY' 0.0 ;
taulim = 'MANUEL' 'CHPO' cvnv 3 'TXX' 0.0 'TYY' 0.0 'TXY' 0.0 ;
gradv0 mchamv = 'PENT' $mt 'FACE' 'DIAMAN2' ('MOTS' 'UX' 'UY')
('MOTS' 'P1DX' 'P1DY' 'P2DX' 'P2DY') v0 v0li1 v0li2 ;
*
t0 = 'MANU' 'CHPO' ('DOMA' $mt 'CENTRE') 1 'SCAL' 0.D0 ;
*
* We impose Dirichlet boundary conditions on cdir
*
cdir = ('DOMA' $facint 'CENTRE') ;
* cvn = set of FACE centers where we impose von Neumann boundary
* conditions
cvn = (('DOMA' $facext 'CENTRE') 'ET' ('DOMA' $faccon 'CENTRE')) ;
'SI' graph ;
'TRACER' (cnt et cvn) 'TITRE' 'Von Neumann B.C.' ;
'TRACER' (cnt et cdir) 'TITRE' 'Dirichlet B.C.' ;
'FINSI' ;
tlim1 = 'REDU' solex cdir ;
gradtl = gradtext 'ET' gradtcon ;
tlim2 = 'NOMC' ('MOTS' 'UX' 'UY') ('MOTS' 'P1DX' 'P1DY')
gradtl ;
gradt0 mchamt = 'PENT' $mt 'FACE' 'DIAMAN2' ('MOTS' 'SCAL')
('MOTS' 'P1DX' 'P1DY') t0 tlim1 tlim2 ;
* Jacobian
* We use the 'NS' jacobian
listinco = 'MOTS' 'RN' 'RUXN' 'RUYN' 'RETN' ;
qlim = -1.0 '*' kappa '*' gradtl ;
jaco chpres dt = 'LAPN' 'VF' 'PROPCOST' 'RESI' 'IMPL'
$mt mu kappa cv rho0 v0 t0 gradv0 gradt0 mchamv mchamt
'VIMP' v0li1 'TAUI' taulim 'QIMP' qlim 'TIMP' tlim1 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' 'DIAMAN2' ('MOTS' 'SCAL')
('MOTS' 'P1DX' 'P1DY') t1 tlim1 tlim2 ;
* We check that the residuum is 0 if computed with
* gradt1 and t1
jaco chpres1 dt1 = 'LAPN' 'VF' 'PROPCOST' 'RESI' 'IMPL'
$mt mu kappa cv rho0 v0 t1 gradv0 gradt1 mchamv mchamt
'VIMP' v0li1 'TAUI' taulim 'QIMP' qlim 'TIMP' tlim1 listinco ;
* jacbid chpres1 dt = 'LAPN' 'VF' 'PROPCOST' 'RESI' 'IMPL'
* $mt mu kappa cv rho0 v0 t1 gradv0 gradt1 mchamv mchamt
* 'QIMP' qlim 'TIMP' tlim1 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.' ;
chm_err = 'KCHA' $mt 'CHAM' (chpres1 '-' chpres) ;
titt = 'CHAINE' 'Erreur sur le residu ' titglob ;
'TRACER' chm_err $mt 'TITRE' titt ;
'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 ;
'SI' ('EGA' dimcon 1) ;
labs = 'PROG' 1. ;
'SINON' ;
labs = 'PROG' 1.D0 'PAS' 1.D0 dimcon ;
'FINSI' ;
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' 7 PAS 1 8 ;
'SINON' ;
lraff = 'LECT' 25 PAS 1 28 ;
'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-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 = 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.9 ) ;
'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' ;
© Cast3M 2003 - All rights reserved.
Disclaimer