Download equ_chaleur3D_VFSYM.dgibi
* fichier : equ_chaleur3D_VFSYM.dgibi
************************************************************************
* NOM : equ_chaleur3D_VF.dgibi
* ___
*
* DESCRIPTION : Solution stationnaire de l'équation de la chaleur (3D)
* ___________
*
* GEOMETRIE : Un cube de côté 1 maillé avec 'VOLU'.
*
* EQUATIONS :
* ----------
*
* EQUATIONS :
* ----------
*
* - Equations :
*
* div K GRAD T = F
*
* |
* avec K = |+1 0 0|
* |0 1 0|
* |0 0 1|
*
*
* - Conditions aux limites :
*
* conditions de Dirichlet
*
* - Solution exacte :
*
* T(x,y)= sin pix sin piy sin pi z
*
*
*
*
* DISCRETISATION : une méthode de Volume Finis d'ordre 2 en espace, avec la
* methode VFSYM .
* SCHEMA PROPOSE PAR Christophe Le Potier
* Références : {C. Le Potier}
* \emph{Schema volumes finis pour des operateurs de diffusion
* fortement anisotropes sur des maillages non structures},
* C. R. Acad. Sci. Ser. I \textbf{340}, 2005, pp. 921--926.
* interfaces
*
*
*
* Le maillage est construit avec l'opérateur VOLU.
*
* Opérateurs utilisés : PENT (option VF implicite)
* 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 : Christophe Le Potier(CEA/DEN/DM2S/SFME/LSET)
* mél : clepotier@cea.fr
************************************************************************
* VERSION : v1, 24/03/2010, version initiale
* 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 ;
*
'OPTION' 'DIME' 3 'ELEM' CUB8 ;
'OPTION' 'ISOV' 'SULI' ;
'OPTION' 'ECHO' 1 ;
nbisov = 15 ;
'SI' ('NON' interact) ;
'OPTION' 'TRAC' 'PSC' ;
* 'OPTION' 'ISOV' 'LIGNE' ;
'FINSI' ;
*******************************************************************
** Erreur Linfini entre deux Champoints.
*
*******************************************************************
DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ;
AUX = (vitp1 - vit) ;
err = MAXI (AUX) 'ABS' ;
err = err / (MAXI(ABS(vit)));
FINPROC err ;
*
*
*******************************************************************
** Erreur Pseudo L2 entre deux Champoints.
*
*******************************************************************
DEBPROC CALCERR2 vitp1*'CHPOINT' vit*'CHPOINT' ;
er0 = 'KOPS' vitp1 '-' vit ;
er2 = 'KOPS' er0 '*' er0 ;
er3 = 'KOPS' vit '*' vit ;
err = 0.D0 ;
erl = 0.D0;
suppv = 'EXTRAIRE' er2 'MAIL' ;
compv = 'EXTRAIRE' er2 'COMP' ;
nbcomp = 'DIME' compv ;
nptd = 'NBNO' suppv ;
'REPETER' numcomp nbcomp ;
lacomp = 'EXTRAIRE' compv &numcomp ;
'REPETER' numpo nptd ;
lepoi = suppv 'POIN' &numpo ;
err = err '+' ('EXTRAIRE' er2 lacomp lepoi) ;
erl = erl '+' ('EXTRAIRE' er3 lacomp lepoi) ;
'FIN' numpo ;
'FIN' numcomp ;
err = err '/' erl ;
err = err ** 0.5D0 ;
FINPROC err ;
*
*******************************************************************
* Procédure paramétrée (raffinement)
* renvoyant l'erreur en norme L2 sur la température.
* On calcule une solution exacte de l'équation de Laplace
* (équation de la chaleur) ;
*
'DEBPROC' CALCUL nraff*'ENTIER' ;
*
* titre global pour les dessins
*
titglob = 'CHAINE' ' ; nraff=' nraff ;
*
* Géométrie
*
pA = 0. 0. 0. ;
pB = 1.0 0. 0. ;
pC = 1.0 1.0 0. ;
pD = 0. 1.0 0. ;
pE = 0. 0. 1. ;
*
* Paramètres de la discrétisation de base
*
'SI' complet ;
nAB = 1 ;
nBC = 1 ;
nCD = 1 ;
nDA = 1 ;
nH = 1 ;
'SINON' ;
nAB = 2 '*' (nraff '+' 1) ;
nBC = 2 '*' (nraff '+' 1) ;
nCD = 2 '*' (nraff '+' 1) ;
nDA = 2 '*' (nraff '+' 1) ;
nH = 2 '*' (nraff '+' 1) ;
'FINSI' ;
*
*******************************************************************
* Géométrie discrétisée
*******************************************************************
*
*'OPTION' 'ELEM' 'TRI3' ;
bas = 'DROIT' nAB pA pB ;
droite = 'DROIT' nBC pB pC ;
haut = 'DROIT' nCD pC pD ;
gauche = 'DROIT' nDA pD pA ;
spourt = bas 'ET' droite 'ET' haut 'ET' gauche ;
sbas = 'SURFACE' spourt 'PLAN' ;
*'OPTION' 'DIME' 3 'ELEM' TET4 ;
mt = 'VOLUME' sbas 'TRAN' nH pE ;
elim mt 0.001;
f1 f2 f3 = 'FACE' mt ;
f1 = 'CHANGER' f1 'TRI3' ;
f2 = 'CHANGER' f2 'TRI3' ;
f3 = 'CHANGER' f3 'TRI3' ;
ss = (f1 'ET' f2 'ET' f3) ;
*'OPTION' 'ELEM' 'TET4' ;
mt = 'VOLUME' ss ;
*'OPTION' 'ELEM' 'TET4' ;
* VERRUE
pourtour = 'ENVELOPPE' mt ;
*
* Eventuellement, on trace le résultat
*
* list mt;
'SI' graph ;
titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt)
' NBELEM=' ('NBEL' mt) titglob ;
'TRACER' mt 'NOEUD' 'TITRE' titgeo ;
'FINSI' ;
*******************************************************************
*
* Modèle
*
*******************************************************************
L = 10.D0;
mt = DEPLA (mt) AFFI (L) (0. 0. 0.) (0. 0. 1.);
Mmt = 'CHAN' mt 'QUAF' ;
'ELIM' Mmt 1.E-5 ;
MDNS = 'DARCY' ;
$mt = 'MODE' Mmt MDNS 'ANISOTROPE' ;
mt = 'DOMA' $mt 'MAILLAGE' ;
MESS 'NOMBRE DE MAILLES' (NBEL (MT));
modmt = 'MODELISER' ('DOMA' $mt 'MAILLAGE') 'THERMIQUE' ;
* list ('DOMA' $mt 'FACEP') ;
mtm = ('DOMA' $mt 'MAILLAGE') ;
mtc = 'COULEUR' ('DOMA' $mt 'CENTRE') 'JAUN' ;
mtf = 'COULEUR' ('DOMA' $mt 'FACE') 'BLEU' ;
cmt = pourtour;
Mcmt = 'CHAN' cmt 'QUAF' ;
'ELIM' (Mmt 'ET' Mcmt) 1.E-5 ;
$cmt = 'MODE' Mcmt MDNS 'ANISOTROPE' ;
*******************************************************************
* Solution exacte aux centres
*******************************************************************
re = 180.;
ep = 1.D0;
xxc yyc zzc = 'COORDONNEE' ('DOMA' $mt 'CENTRE') ;
solexc = (sin(re*xxc))*(sin(re*yyc))*(sin(re*zzc/L));
solexc = 'KCHT' $mt 'SCAL' 'CENTRE' solexc ;
*******************************************************************
* terme source
*******************************************************************
sour = 0.0D0 -(((2.D0 + (ep*(L**-2))) *(pi*pi))*solexc);
sour = 'KCHT' $mt 'SCAL' 'CENTRE' sour ;
*******************************************************************
* CONDITION DE DIRICHLET
*******************************************************************
xxlim yylim zzlim = 'COORDONNEE' ('DOMA' $cmt 'CENTRE') ;
solim = (sin(re*xxlim))*(sin(re*yylim))*(sin(re*zzlim/L));
solim = 'KCHT' $cmt 'SCAL' 'CENTRE' solim ;
*******************************************************************
*
* Calcul du tenseur
*
*******************************************************************
K11 = 'KOPS' xxc '+' 1.D0 ;
K22 = 'KOPS' yyc '+' 1.D0 ;
K33 = 'KOPS' zzc '+' 1.D0 ;
K11 = 1.D0 + (0.0*xxc);
K22 = 1.D0 + (0.0*yyc);
K33 = 1.D0 + (0.0*zzc);
K21 = (0.0*zzc);
K31 = (0.0*zzc);
K32 = (0.0*zzc);
K11 = NOMC(K11) 'K11' 'NATU' 'DISCRET';
K22 = NOMC(K22) 'K22' 'NATU' 'DISCRET';
K33 = NOMC(K33) 'K33' 'NATU' 'DISCRET';
K21 = NOMC(K21) 'K21' 'NATU' 'DISCRET';
K31 = NOMC(K31) 'K31' 'NATU' 'DISCRET';
K32 = NOMC(K32) 'K32' 'NATU' 'DISCRET';
PERM = (K11 'ET' K22) ;
PERM = PERM 'ET' K33;
PERM = PERM 'ET' K21;
PERM = PERM 'ET' K31;
PERM = PERM 'ET' K32;
* 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_{\partial \Omega} donné (CL de Dirichlet)
* 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.
*
*
t0 = 'KCHT' $mt 'SCAL' 'CENTRE' 0.0 ;
gradt0 mchamt = 'PENT' $mt
'FACE' 'VFSYM' t0
'DISPDIF' PERM 'TIMP' solim
;
jaco chpres dt = 'LAPN' 'VF' 'CLAUDEIS' 'IMPL'
$mt t0 gradt0 mchamt 'TIMP' solim ;
mamat = 'KOPS' 'MULT' -1.0D0 jaco ;
matot = mamat ;
AUX = ('EXCO' chpres 'RETN' 'RETN') - (NOMC(sour) 'RETN');
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' AUX
'IMPR' 0 ;
t1 = t0 '+' ('EXCO' deltat 'RETN' 'SCAL') ;
*******************************************************************
* Verification de la convergence
*******************************************************************
gradt1 = 'PENT' $mt 'FACE' 'VFSYM' t1
'DISPDIF' PERM 'TIMP' solim
'GRADGEO' mchamt;
jacbid chpres1 dt = 'LAPN' 'VF' 'CLAUDEISl' 'IMPL'
$mt t0 gradt1 mchamt 'TIMP' solim ;
AUX = ('EXCO' chpres1 'RETN' 'RETN') - (NOMC(sour) 'RETN');
res = 'KRES' matot 'TYPI' (rv . 'METHINV')
'SMBR'AUX
'IMPR' 0 ;
mres1 = 'MAXIMUM' res 'ABS' ;
'MESSAGE' ('CHAINE' 'Maxi. mres1 =' 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_tn = 'KCHA' $mt 'CHAM' tn ;
titt = 'CHAINE' 'Température exacte' titglob ;
'TRACER' chm_tn modmt 'TITRE' titt ;
*
* graphe de convergence de la méthode itérative
*
conver = (rv . 'METHINV' . 'CONVINV') ;
dimcon = 'DIME' conver ;
'SI' ('>' dimcon 1) ;
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' ;
'FINSI' ;
*
* solutions calculées
*
tn = t1 ;
chm_tn = 'KCHA' $mt 'CHAM' tn ;
titt = 'CHAINE' 'Température calculée' titglob ;
'TRACER' chm_tn modmt 'TITRE' titt ;
* erreur
*
erl = ABS(solexc '-' tn)/(maxi solexc);
erlcham = 'KCHA' $mt 'CHAM' erl ;
titt = 'CHAINE' 'erreur relative' ;
'TRACER' (erlcham) modmt 'TITRE' titt ;
'FINSI' ;
*
* Calcul des erreurs par rapport à la solution analytique
*
echloc = '**' (('MESURE' mt) '/' ('NBEL' mt))
(1.D0 '/' ('VALEUR' 'DIME')) ;
tn = t1 ;
errlit = CALCERR tn solexc ;
errl2t = CALCERR2 tn solexc ;
'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 2 ;
'SINON' ;
lraff = 'LECT' 0 PAS 1 5 ;
'FINSI' ;
*
*-----------------------------------------------------------*
* Boucles sur les différents paramètres du calcul *
ordok = VRAI ;
precok = VRAI ;
* ordre théorique en espace de la méthode
ordtth = 2 ;
* erreur L2 max pour la solution (raffinement 1, complet=FAUX)
errtmax = 1.6605D-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 = ('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 = 'ET' ordok ('EGA' ('ENTIER' ('+' ordt 0.5)) ordtth) ;
'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 ;
'FINSI' ;
'MESSAGE' 'Tout s_est bien passé' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer