* fichier : equ_chaleur3D_VF2.dgibi ************************************************************************ * Section : Thermique Statique ************************************************************************ *********************************************************************** * 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' 'DIME' 3 'ELEM' TYPEL 'MODE' 'PLAN' ; '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' ; 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 ; * * 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' ; mt = 'VOLUME' sbas 'TRAN' nH pE ; * 'SI' ('NON' logreg) ; * * TET4 * f1 f2 f3 = 'FACE' mt ; 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 ; titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt) ' NBELEM=' ('NBEL' mt) titglob ; 'TRACER' mt 'TITRE' titgeo ; 'FINSI' ; * * Definition des cotés gauche et bas * cnt = 'ENVE' mt 'COULEUR' 'VERT' ; 'SI' (graph ET FAUX) ; 'TRACER' (mt 'ET' cnt) 'TITRE' 'Enveloppe' ; 'FINSI' ; * * Creation des modèles * * MODELS $mt = 'MODE' mt 'EULER' ; $cnt = 'MODE' cnt 'EULER' ; $f1 = 'MODE' f1 'EULER' ; * Tmt = 'DOMA' $mt 'VF' ; Tcnt = 'DOMA' $cnt 'VF' ; Tf1 = 'DOMA' $f1 'VF' ; * QUAF Mmt = Tmt . 'QUAF' ; Mcnt = Tcnt . 'QUAF' ; Mf1 = Tf1 . 'QUAF' ; 'ELIM' (Mmt 'ET' Mcnt 'ET' Mf1) 1.E-8 ; * CHAMPOINT vide et MATRIK vide chvid matvid = 'KOPS' 'MATRIK' ; * * Solution exacte bilinéaire (sur le centre du contour) * xx yy zz = 'COORDONNEE' ('DOMA' $cnt 'CENTRE') ; solex = (xx '*' alpha) '+' (yy '*' beta) '+' (zz '*' gamma) '+' (xx '*' yy '*' zz '*' eta) '+' delta ; * * Solution exacte aux centres * xxc yyc zzc = 'COORDONNEE' ('DOMA' $mt 'CENTRE') ; solexc = (xxc '*' alpha) '+' (yyc '*' beta) '+' (zzc '*' gamma) '+' (xxc '*' yyc '*' zzc '*' eta) '+' delta ; * * Gradient exact (aux faces) * xxf yyf zzf = 'COORDONNEE' ('DOMA' $mt 'FACE') ; gradtx = 'NOMC' 'UX' (alpha '+' (yyf '*' zzf '*' eta)) 'NATU' 'DISCRET' ; gradty = 'NOMC' 'UY' (beta '+' (xxf '*' zzf '*' eta)) 'NATU' 'DISCRET' ; gradtz = 'NOMC' 'UZ' (gamma '+' (xxf '*' yyf '*' eta)) 'NATU' '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. * * rho0 = 'MANUEL' 'CHPO' ('DOMA' $mt 'CENTRE') 1 'SCAL' 1.0 ; v0 = 'MANUEL' 'CHPO' ('DOMA' $mt 'CENTRE') 3 'UX' 0.0 'UY' 0.0 'UZ' 0.0 ; cdirv = 'MANUEL' 'POI1' ('POIN' 1 ('DOMA' $cnt 'CENTRE')) ; cvnv = 'DIFF' ('DOMA' $cnt 'CENTRE') cdirv ; v0li1 = 'MANUEL' 'CHPO' cdirv 3 'UX' 0.0 'UY' 0.0 'UZ' 0.0 ; v0li2 = 'MANUEL' 'CHPO' cvnv 9 'P1DX' 0.0 'P1DY' 0.0 'P1DZ' 0.0 'P2DX' 0.0 'P2DY' 0.0 'P2DZ' 0.0 'P3DX' 0.0 'P3DY' 0.0 'P3DZ' 0.0 ; taulim = 'MANUEL' 'CHPO' cvnv 6 'TXX' 0.0 'TYY' 0.0 'TXY' 0.0 'TZZ' 0.0 'TXZ' 0.0 'TYZ' 0.0 ; gradv0 mchamv = 'PENT' $mt 'FACE' 'DIAMAN2' ('MOTS' 'UX' 'UY' 'UZ') ('MOTS' 'P1DX' 'P1DY' 'P1DZ' 'P2DX' 'P2DY' 'P2DZ' 'P3DX' 'P3DY' 'P3DZ') v0 v0li1 v0li2 ; * t0 = 'MANU' 'CHPO' ('DOMA' $mt 'CENTRE') 1 'SCAL' 0.D0 ; * * We impose Dirichlet boundary conditions on cdir * cdir = ('DOMA' $f1 'CENTRE') ; * cvn = set of FACE centers where we impose von Neumann boundary * conditions cvn = ('DIFF' ('DOMA' $cnt 'CENTRE') cdir) 'COUL' 'TURQ' ; '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 = 'REDU' (gradtx 'ET' gradty 'ET' gradtz) cvn ; tlim2 = 'NOMC' ('MOTS' 'UX' 'UY' 'UZ') ('MOTS' 'P1DX' 'P1DY' 'P1DZ') gradtl ; gradt0 mchamt = 'PENT' $mt 'FACE' 'DIAMAN2' ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY' 'P1DZ') t0 tlim1 tlim2 ; * Jacobian * We use the 'NS' jacobian listinco = 'MOTS' 'RN' 'RUXN' 'RUYN' 'RUZN' 'RETN' ; qchal = -1.0 '*' kappa '*' (gradtx 'ET' gradty 'ET' gradtz) ; qlim = 'REDU' qchal cvn ; 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' 'RUZN') ('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' 'P1DZ') 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 ; 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' ; * * 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)) ** (1. '/' 3.) ; 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 ; lraff = 'LECT' 7 PAS 1 8 ; 'SINON' ; lraff = 'LECT' 4 PAS 1 5 ; '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' ;