Télécharger smithhutton_cvg.dgibi
* fichier : smithhutton_cvg.dgibi ************************************************************************ ************************************************************************ ************************************************************************ * NOM : SMITHHUTTON_CVG * DESCRIPTION : Une variante du cas-test de convection-diffusion d'un * scalaire dû à Smith et Hutton [1] * (voir aussi smithhutton.dgibi) * * On cherche à vérifier les ordres de convergence spatiale pour les * éléments quadratiques et linéaires avec une méthode de discrétisation * 'CENTREE' pour les termes de convection et à Péclet=20. * * Pour ce faire, on construit une solution de référence sur un maillage * fin et on calcule les ordres avec des maillages plus grossiers. * On compare les évolutions de la solution sur la frontière de sortie * avec une norme L2 discrète. * * Pour retrouver les ordres de convergence, il y a plusieurs subtilités * : * - la première est que le cas proposé par Smith et Hutton présente une * singularité au point (0,0) : en effet, on passe d'une condition à la * limite de Dirichlet (température imposée) à gauche, à une condition de * Neumann (flux diffusif nul) à droite, sur le bord rectiligne bas du * domaine. Il s'agit de la singularité Dirichlet-Neumann expliquée, par * exemple, dans mon cours sur le Laplacien p.26 [2]. Les ordres * théoriques de convergence sont de 1 lorsque la solution cherchée * présente ce type de singularité. * * Pour s'affranchir de cette singularité, on peut modifier le cas en ne * gardant que le demi-domaine x>0 et en imposant la condition d'entrée * en tangente hyperbolique sur le côté gauche (x=0) de ce demi-domaine. * On introduit la variable logique demi : si demi=faux, on a le cas * singulier et si demi=vrai, on a le cas régulier. * * - la deuxième subtilité est que lorsqu'on utilise IPOL pour calculer * l'erreur, on introduit également une erreur en O(h^2) car IPOL utilise * une interpolation linéaire par morceaux ! Pour remédier à ce problème, * on peut utiliser une interpolation par spline cubique en O(h^4) en * précisant le mot-clé 'SPLINE' à IPOL. On introduit la variable logique * spli pour ce faire dans le jeu de données. * * Les résultats de l'étude avec ces deux paramètres plus le type * d'éléments finis sont les suivants (complet = vrai) : * * EF demi spli || ordre p de minimum du log10 * || convergence de l'erreur * ------------------------------------------------------------------ * QUAF VRAI VRAI || 3.6 -7.5 * QUAF VRAI FAUX || 1.9 -4.7 * * LINE VRAI VRAI || 2.2 -4.7 * ------------------------------------------------------------------ * QUAF FAUX VRAI || 1.7 -3.0 * LINE FAUX VRAI || 1.3 -2.8 * ------------------------------------------------------------------ * * Les deux premières lignes illustrent le fait que si l'on n'utilise pas * l'option SPLINE, on massacre le calcul de l'erreur lorsque celle-ci * est petite. Les lignes 1 et 3 donnent des ordres de convergence * consistants avec la théorie pour le cas régulier : ordre 3 pour les * éléments quadratiques (ici 3.6) et 2 pour les éléments linéaires (ici * 2.2). * * Pour le cas singulier (les deux dernières lignes), on observe une * dégradation des ordres de convergence pour les deux types d'éléments * finis. On note que les erreurs sont également beaucoup plus grandes * que dans le cas régulier. L'ordre théorique est de 1, mais il * faudrait raffiner encore plus pour l'obtenir. * * Il serait intéressant d'investiguer le cas du décentrement SUPG et des * Péclets plus élevés pour voir si les résultats restent cohérents avec * ces explications. * * Pour le cas complet = faux avec des maillages moins fins (résultats * moins précis) : * * EF demi spli || ordre p de minimum du log10 * || convergence de l'erreur * ------------------------------------------------------------------ * QUAF VRAI VRAI || 6.0 -6.6 * QUAF VRAI FAUX || 2.2 -4.0 * * LINE VRAI VRAI || 2.7 -3.6 * ------------------------------------------------------------------ * QUAF FAUX VRAI || 2.9 -3.1 * LINE FAUX VRAI || 1.7 -2.4 * ------------------------------------------------------------------ * * Références : * [1] Smith, R. M. and Hutton, A. G., THE NUMERICAL TREATMENT OF * ADVECTION: A PERFORMANCE COMPARISON OF CURRENT METHODS , Numerical * Heat Transfer vol 5, n4, pp 439-461, 1982 * [2] S. Gounand, Introduction à la méthode des éléments finis en * mécanique des fluides incompressibles, 2012, Cours ENSTA B2-1 * * * LANGAGE : GIBIANE-CAST3M * AUTEURS : Stéphane GOUNAND (CEA/DEN/DM2S/STMF/LMSF) * Guillaume BARBA ROSSA (Ecole Polytechnique) * mél : stephane.gounand@cea.fr * ********************************************************************** * VERSION : v1, 22/10/2013, version initiale * HISTORIQUE : v1, 22/10/2013, création * HISTORIQUE : * HISTORIQUE : ************************************************************************ * 'SAUTER' 'PAGE' ; * interact= FAUX ; complet = FAUX ; graph = FAUX ; * 'SI' ('NON' interact) ; 'OPTION' 'TRAC' 'PSC' ; 'SINON' ; 'OPTION' 'TRAC' 'X' ; 'FINSI' ; ************************************************************************ * Procédure pour changer les maillages en QUAF (quadratique fluide) * faire les éliminations et construire les modèles Navier_Stokes. * * SYNTAXE * mail1 $mail1 mail2 $mail2 ... = MODIF 'LINE' mail1 mail2 ..... * mail1,2,... maillages avant et apres transformation en QUAF * $mail1,2 ... modeles associes aux maillages. * on boucle sur les arguments d'entree * on cree une table pour stocker les arguments toto = 'TABLE' ESCLAVE; tata = 'TABLE' ; 'REPETER' barg ; * On lit le premier argument 'ARGUMENT' titi/'MAILLAGE'; 'SI' ('EXISTE' titi) ; toto . &barg = 'CHANGER' titi 'QUAF' ; tata . &barg = 'MODELISER' (toto . &barg) 'NAVIER_STOKES' discr; 'SINON' ; nargu = &barg '-' 1; 'QUITTER' barg; 'FINSI' ; 'FIN' barg ; 'ELIMINATION' ('ET' toto) 1.D-6 ; 'REPETER' bsort nargu; * On sort les arguments * 'MESSAGE' 'coucou' &bsort; 'RESPRO' (tata . &bsort); * 'RESPRO' (tata . &bsort); 'FIN' bsort ; 'FINPROC' ; *************************************************************** * Procédure appelée dans EQEX pour tracer les champs en cours * de calcul. * 'DEBPROC' TRINCO; 'ARGUMENT' RVX * 'TABLE'; 'TRACER' 'NCLK' (RV.'INCO'.'CN') MT; 'RESPRO' CHVID MATVID; 'FINPROC'; *************************************************************** * Procédure effectuant la résolution du problème de * convection-diffusion et renvoyant l'erreur ou l'évolution * de référence (suivant la valeur de REF) * 'DEBPROC' CALCUL ; *'ARGUMENT' hh * 'FLOTTANT' REF * 'ENTIER'; 'ARGUMENT' NX*'ENTIER' REF*'ENTIER'; 'ARGUMENT' Pe*'FLOTTANT' ; 'ARGUMENT' demi*'LOGIQUE' ; 'ARGUMENT' spli*'LOGIQUE' ; NY =NX ; *MESSAGE (CHAINE 'NX = ' NX ' NY = ' NY); *MESSAGE (CHAINE 'MET = ' MET ' DISC = ' DISC); *MESSAGE (CHAINE 'Pe = ' Pe ' demi=' demi ' spli=' spli); *'h = ' hh ' L=1.; ML = '*' L -1. ; H=1.; * O=0. 0. ; A=0. H ; B=ML 0.; C=ML H ; D=L H ; E=L 0.; * AO= 'DROIT' NY A O; OB= 'DROIT' NX O B; BC= 'DROIT' NY B C; CA= 'DROIT' NX C A; AD= 'DROIT' NX A D; DE= 'DROIT' NY D E; OE= 'DROIT' NX O E; M1='DALLER' AO OB BC CA; M2='DALLER' AO OE DE AD; 'SI' demi ; MT = M2; 'SINON' ; MT = M1 'ET' M2; 'FINSI' ; *'TRACER' mt ; *'ELIMINATION' 1E-4 MT ; 'SI' demi ; MT $MT LIMP $LIMP LIN $LIN LOU $LOU = MODIF MT (AD 'ET' DE) AO OE DISC ; 'SINON' ; MT $MT LIMP $LIMP LIN $LIN LOU $LOU = MODIF MT (BC 'ET' CA 'ET' AD 'ET' DE) OB OE DISC ; 'FINSI' ; *LIST (DOMA $MT 'TABLE'); *TRACER ((DOMA $MT 'CENTRE') 'ET' (DOMA $MT 'MAILLAGE')) QUAL; X Y='COORDONNEE' MT; 'SI' demi ; CIN=1 '+' ('TANH' (10 * ((-2 * Y) '+' 1) )); CIMP=(0.0*X)+(1 '-' ('TANH' 10)); 'SINON' ; CIN=1 '+' ('TANH' (10 * ((2 * X) '+' 1) )); CIMP=(0.0*X)+(1 '-' ('TANH' 10)); 'FINSI' ; *EVC='EVOL' 'CHPO' CIN LIN; *DESSSIN EVC; *'TRACER' CIMP MT; V=UX0 + UY0; *CVEC='VECT' V 'UX' 'UY'; *TRACER CVEC MT; *description des equations 'SI' ('EGA' Pe 0. 1.D-6) ; iPe = 1. ; 'SINON' ; iPe = '/' 1. Pe ; 'FINSI' ; RV = 'EQEX' 'OPTI' 'EF' 'IMPL' * 'ITMA' 20 * 'NITER' 10 'ZONE' $MT 'OPER' 'LAPN' iPe 'INCO' 'CN' ; *RV = 'EQEX' RV * 'OPTI' 'EF' 'IMPL' 'CENTREE' * 'ZONE' $MT * 'OPER' 'DFDT' 1. 'CN' 1. 'INCO' 'CN' ; *'SI' (REF 'EGA' 1); * MET='CENTREE'; * **SINON; **toujours centree pour la reference * MET='CENTREE'; *FINSI; *MESSAGE (CHAINE 'methode : ' MET); 'SI' ('NEG' Pe 0. 1.D-6) ; RV = 'EQEX' RV 'OPTI' 'EF' 'IMPL' MET 'ZONE' $MT 'FINSI' ; *RV = 'EQEX' RV * 'ZONE' $MT * 'OPER' 'TRINCO' ; *description des inconnues, CI RV.'INCO' = 'TABLE' 'INCO'; RV.'INCO'.'VITESSE' = V; *conditions aux limites RV = 'EQEX' RV 'CLIM' 'CN' 'TIMP' (LIMP 'ET' LIN) (CIN 'ET' CIMP) ; rv . 'METHINV' . 'TYPINV' = 3 ; rv . 'METHINV' . 'IMPINV' = 1 ; rv . 'METHINV' . 'RESID' = 1.D-20 ; *resolution EXEC RV; CN = RV.'INCO'.'CN'; *tracer de la solution *'TRACER' CN MT; *post-traitement EVCO = EVCO 'COULEUR' 'ROUGE' ; 'SI' spli ; 'SINON' ; 'FINSI' ; *'LISTE' COUI; err=0; SI (REF 'EGA' 0); err = (err '/' NOUT) '**' 0.5 ; * 'LISTE' err ; 'MESSAGE' ('CHAINE' 'erreur =' err); * 'DESSIN' (EVCO 'ET' EVEX) 'GRILL' ; FINSI; 'RESPRO' err COUI; 'FINPROC'; *************************************************************** * Procédure appelant CALCUL et effectuant le calcul de l'ordre * de convergence. Elle renvoie également la valeur minimale de * l'erreur. * 'DEBPROC' ODCVG ; 'ARGUMENT' complet*'LOGIQUE' ; 'ARGUMENT' Pe*'FLOTTANT' ; 'ARGUMENT' demi*'LOGIQUE' ; 'ARGUMENT' spli*'LOGIQUE' ; * 'SI' complet ; 'SI' ('EGA' disc 'QUAF') ; 'SINON' ; 'FINSI' ; 'SINON' ; 'SI' ('EGA' disc 'QUAF') ; 'SINON' ; 'FINSI' ; 'FINSI' ; NXREF = 'EXTRAIRE' LNX 1 ; *MET = 'CENTREE' ; *DEMI = vrai ; SPLI = vrai ; * *XOUC = PROG 0. 'PAS' 0.1 1.0 ; href=2.0 * (1.0 '/' (FLOTTANT (NXREF))); err COUC=CALCUL nxref 1 MET DISC Pe DEMI spli ; EVEX = EVEX 'COULEUR' 'BLEU' ; * lh = '**' ('FLOTTANT' lnxr) -1. ; * *lh='PROG' 0.1 'PAS' -0.01 0.01; 'REPETER' bcl NLH; hh='EXTRAIRE' lnxr &bcl; err COUI=CALCUL hh 0 MET DISC Pe DEMI spli ; 'FIN' bcl; llh = '/' ('LOG' lh) ('LOG' 10.) ; llerr = '/' ('LOG' lerr ) ('LOG' 10.) ; *DESSIN EVERR 'GRILL'; * determination de la pente *X0='LOG' ('EXTRAIRE' lh NLH); *X0 = '/' X0 ('LOG' 10.) ; *X1='LOG' ('EXTRAIRE' lh ('-' NLH 1)); *X1 = '/' X1 ('LOG' 10.) ; *EVEI = 'EXTRAIRE' EVEI 'APRE' X0; ordp = TCN . 1 ; * tit = 'CHAINE' 'MET = ' MET ' DISC = ' DISC ' Pe = ' Pe ' demi=' demi ' spli=' spli ; 'SI' graph ; 'FINSI' ; 'MESSAGE' tit ; tr = 'CHAINE' 'ord.=' ordp ' min. log10 err.=' minr ; ; 'MESSAGE' tr ; 'RESPRO' ordp minr ; 'FINPROC' ; *************************************************************** * Jeu de donnés effectuant le test des ordres de convergences * pour diverses options. * Pe=20. ; tabtest = 'TABLE' ; tabtest . 1 = 'TABLE' ; tabtest . 1 . 'disc' = 'QUAF' ; tabtest . 1 . 'demi' = VRAI ; tabtest . 1 . 'spli' = VRAI ; tabtest . 1 . 'ordtheo' = 3 ; 'SI' complet ; tabtest . 1 . 'errminr' = -7.4 ; 'SINON' ; tabtest . 1 . 'errminr' = -6.5 ; 'FINSI' ; tabtest . 2 = 'TABLE' ; tabtest . 2 . 'disc' = 'QUAF' ; tabtest . 2 . 'demi' = VRAI ; tabtest . 2 . 'spli' = FAUX ; tabtest . 2 . 'ordtheo' = 2 ; 'SI' complet ; tabtest . 2 . 'errminr' = -4.6 ; 'SINON' ; tabtest . 2 . 'errminr' = -3.9 ; 'FINSI' ; tabtest . 3 = 'TABLE' ; tabtest . 3 . 'disc' = 'LINE' ; tabtest . 3 . 'demi' = VRAI ; tabtest . 3 . 'spli' = VRAI ; tabtest . 3 . 'ordtheo' = 2 ; 'SI' complet ; tabtest . 3 . 'errminr' = -4.6 ; 'SINON' ; tabtest . 3 . 'errminr' = -3.5 ; 'FINSI' ; tabtest . 4 = 'TABLE' ; tabtest . 4 . 'disc' = 'QUAF' ; tabtest . 4 . 'demi' = FAUX ; tabtest . 4 . 'spli' = VRAI ; tabtest . 4 . 'ordtheo' = 1 ; 'SI' complet ; tabtest . 4 . 'errminr' = -2.9 ; 'SINON' ; tabtest . 4 . 'errminr' = -2.9 ; 'FINSI' ; tabtest . 5 = 'TABLE' ; tabtest . 5 . 'disc' = 'LINE' ; tabtest . 5 . 'demi' = FAUX ; tabtest . 5 . 'spli' = VRAI ; tabtest . 5 . 'ordtheo' = 1 ; 'SI' complet ; tabtest . 5 . 'errminr' = -2.7 ; 'SINON' ; tabtest . 5 . 'errminr' = -2.4 ; 'FINSI' ; * Les calculs tt = tabtest . &itest ; disc = tt . 'disc' ; demi = tt . 'demi' ; spli = tt . 'spli' ; ordp minr = ODCVG complet Pe disc 'CENTREE' demi spli ; topt = 'CHAINE' 'disc=' disc ' demi=' demi ' spli=' spli ; tt . 'topt' = topt ; tt . 'ordp' = ordp ; tt . 'minr' = minr ; 'FIN' itest ; * Les tests lok = VRAI ; tt = tabtest . &itest ; ordtheo = tt . 'ordtheo' ; errminr = tt . 'errminr' ; topt = tt . 'topt' ; ordp = tt . 'ordp' ; minr = tt . 'minr' ; 'MESSAGE' ('CHAINE' 'Test ' &itest ' ' topt) ; 'MESSAGE' ('CHAINE' 'ordtheo =' ordtheo ' ordcalc =' ordp) ; 'MESSAGE' ('CHAINE' 'minrtheo=' errminr ' minrcalc=' minr) ; tok = (ordp '>' ('-' ordtheo 0.1)) 'ET' (minr '<' errminr) ; 'SI' ('NON' tok) ; 'MESSAGE' '!!!!! Pas bon !!!!!!!' ; 'FINSI' ; lok = lok 'ET' tok ; 'FIN' itest ; * 'SAUTER' 2 'LIGNE' ; 'SI' lok ; 'MESSAGE' 'Tout sest bien passe' ; 'SINON' ; 'MESSAGE' 'Il y a eu des erreurs' ; 'FINSI' ; 'SAUTER' 2 'LIGNE' ; 'SI' ('NON' lok) ; 'ERREUR' 5 ; 'FINSI' ; * * Fin des tests * ************************************************************************ 'SI' interact ; 'OPTION' 'DONN' 5 'ECHO' 1 ; 'FINSI' ; * 'FIN' ;
© Cast3M 2003 - Tous droits réservés.
Mentions légales