Télécharger equ_chaleur3D_VF2.dgibi
* fichier : equ_chaleur3D_VF2.dgibi ************************************************************************ ************************************************************************ *********************************************************************** * NOM : equ_chaleur3D_VF2.dgibi * ___ * * DESCRIPTION : Solution stationnaire de l'équation de la chaleur (3D) * ___________ * * GEOMETRIE : Un cube de côté 1 maillé avec 'VOLU'. * ---------- * * * EQUATIONS : * ---------- * * - Equations : * * mu Laplacien(T) = 0 avec mu = 1 ici * * - Conditions aux limites : * * conditions de Dirichlet sur une face * condition de von Neumann sur le reste de l'enveloppe. * * - Solution exacte : * * T(x,y)= alpha x + beta y gamma z + eta xyz + delta * * * DISCRETISATION : une méthode de Volume Finis d'ordre 2 en espace. * ______________ * * * * * Le maillage est construit avec l'opérateur VOLU, il est perturbé * par un bruit blanc gaussien. * * 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 ; logreg = FAUX ; TYPEL = 'CUB8' ; 'OPTION' 'ISOV' 'SULI' ; 'OPTION' 'ECHO' 0 ; nbisov = 15 ; 'SI' ('NON' interact) ; 'OPTION' 'TRAC' 'PS' ; 'OPTION' 'ISOV' 'LIGNE' ; 'FINSI' ; * ** Erreur Linfini entre deux Champoints. * DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ; 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 ; er2 = 'PSCAL' er2 er2 compv compv ; suppv = 'EXTRAIRE' vol 'MAIL' ; 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 ; * * Conditions aux limites et solution exacte: * * solex = (alpha '*' x) '+' (beta '*' y) '+' (gamma '*' z) '+' * (eta '*' x '*' y '*' z) '+' delta * alpha = (** 2. 0.5) ; beta = (** 3. 0.5) ; gamma = (** 5. 0.5) ; eta = (** 7. 0.5) ; delta = (** 1.5 0.5) ; * * Géométrie * pA = 0. 0. 0. ; pB = 1. 0. 0. ; pC = 1. 1. 0. ; pD = 0. 1. 0. ; pE = 0. 0. 1. ; pF = 1. 0. 1. ; pG = 1. 1. 1. ; pH = 0. 1. 1. ; * * Paramètres de la discrétisation de base * nAB = (nraff '+' 1) ; nBC = (nraff '+' 1) ; nCD = (nraff '+' 1) ; nDA = (nraff '+' 1) ; nH = (nraff '+' 1) ; * * Géométrie discrétisée * 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' ; * 'SI' ('NON' logreg) ; * * TET4 * f1 = 'CHANGER' f1 'TRI3' ; f2 = 'CHANGER' f2 'TRI3' ; f3 = 'CHANGER' f3 'TRI3' ; ss = (f1 'ET' f2 'ET' f3) ; mt = 'VOLUME' ss ; 'SINON' ; f1 = sbas ; 'FINSI' ; * * Eventuellement, on trace le résultat * 'SI' graph ; 'TRACER' mt 'TITRE' titgeo ; 'FINSI' ; * * Definition des cotés gauche et bas * 'SI' (graph ET FAUX) ; 'TRACER' (mt 'ET' cnt) 'TITRE' 'Enveloppe' ; 'FINSI' ; * * Creation des modèles * * MODELS * * QUAF Mmt = Tmt . 'QUAF' ; Mcnt = Tcnt . 'QUAF' ; Mf1 = Tf1 . 'QUAF' ; * CHAMPOINT vide et MATRIK vide * * Solution exacte bilinéaire (sur le centre du contour) * solex = (xx '*' alpha) '+' (yy '*' beta) '+' (zz '*' gamma) '+' (xx '*' yy '*' zz '*' eta) '+' delta ; * * Solution exacte aux centres * solexc = (xxc '*' alpha) '+' (yyc '*' beta) '+' (zzc '*' gamma) '+' (xxc '*' yyc '*' zzc '*' eta) '+' delta ; * * Gradient exact (aux faces) * 'DISCRET' ; 'DISCRET' ; 'DISCRET' ; * * 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. * * 'UZ' 0.0 ; 'P2DX' 0.0 'P2DY' 0.0 'P2DZ' 0.0 'P3DX' 0.0 'P3DY' 0.0 'P3DZ' 0.0 ; 'TZZ' 0.0 'TXZ' 0.0 'TYZ' 0.0 ; ('MOTS' 'P1DX' 'P1DY' 'P1DZ' 'P2DX' 'P2DY' 'P2DZ' 'P3DX' 'P3DY' 'P3DZ') v0 v0li1 v0li2 ; * * * We impose Dirichlet boundary conditions on cdir * * cvn = set of FACE centers where we impose von Neumann boundary * conditions 'SI' graph ; 'TRACER' (cnt et cvn) 'TITRE' 'Von Neumann B.C.' ; 'TRACER' (cnt et cdir) 'TITRE' 'Dirichlet B.C.' ; 'FINSI' ; gradtl ; * Jacobian * We use the 'NS' jacobian qchal = -1.0 '*' kappa '*' (gradtx 'ET' gradty 'ET' gradtz) ; $mt mu kappa cv rho0 v0 t0 gradv0 gradt0 mchamv mchamt 'VIMP' v0li1 'TAUI' taulim 'QIMP' qlim 'TIMP' tlim1 listinco ; * Identity matrix for momentum and density contribution matot = 'ET' mamat mati ; * rv . 'METHINV' . 'TYPINV' = 3 ; rv . 'METHINV' . 'PRECOND' = 3 ; rv . 'METHINV' . 'MATASS' = matot ; rv . 'METHINV' . 'MAPREC' = matot ; 'IMPR' 0 ; * We check that the residuum is 0 if computed with * gradt1 and t1 $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 ; 'SI' ('>' mres1 1.e-5) ; 'MESSAGE' 'La méthode de Newton na pas converge en un pas.' ; titt = 'CHAINE' 'Erreur sur le residu ' titglob ; 'TRACER' chm_err $mt 'TITRE' titt ; 'ERREUR' 5 ; 'FINSI' ; * * Résultats * 'SI' graph ; * * solutions exactes * tn = solexc ; 'TRACER' chm_tnex $mt nbisov 'TITRE' titt ; * * graphe de convergence de la méthode itérative * conver = (rv . 'METHINV' . 'CONVINV') ; lord = ('LOG' conver) '/' ('LOG' 10.D0) ; titev = 'CHAINE' 'Historique de convergence' titglob ; 'DESSIN' evtot 'TITR' titev 'LEGE' ; * * solutions calculées * tn = t1 ; 'TRACER' chm_tn $mt nbisov 'TITRE' titt ; * * erreur * titglob ; 'TRACER' ('ABS' (chm_tnex '-' chm_tn)) $mt nbisov 'TITRE' titt ; 'FINSI' ; * * Calcul des erreurs par rapport à la solution analytique * 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' '-------------------------------------------------' ; '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 ; 'SINON' ; '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 * 'FIN' iraff ; lh = lh '/' ('EXTRAIRE' lh 1) ; lh = ('LOG' lh) '/' ('LOG' 10.D0) ; lerl2t = ('LOG' lerl2t) '/' ('LOG' 10.D0) ; 'Err. tempér.' lerl2t ; * * Recherche des ordres de convergence * ordt = cpt . 1 ; * * 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' ; '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' ; 'REPETER' BL1 ndim ; 'MESSAGE' ('EXTRAIRE' lh &BL1) ' ' ('EXTRAIRE' lerl2t &BL1) ; 'FIN' BL1 ; lh1 = 'EXTRAIRE' tlipol 'ABSC' ; '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 - Tous droits réservés.
Mentions légales