Darcy et équation de chaleur
On discrétisé la vitesse v par l'élément de Raviart-Thomas (RT0), P0 pour la pression p et P1 pour la température T.
J’ai eu plusieurs problèmes.
Le premier est que la solution approchée n'est pas valable pour notre problème, c-à-d le code en n'a pas encore [été] validé et je ne sais pas est ce que si le code que j'utilise est valable ou non ?
Le deuxième est de calculer l'erreur, j'ai traité un code mais le compilateur n'accepte pas ce code et je ne sais pas où le problème ? $$
erreur=||u-u_h||_{L^2(\Omega)^d}+||\mathrm{div}(u-u_h)||_{L^2(\Omega)}+||p-p_h||_{L^2(\Omega)}+|T-T_h|_{H^1(\Omega)}.
$$ Voici le code que j'ai traité et des fichiers pdf joints
Yassine OUZROUR.
J’ai eu plusieurs problèmes.
Le premier est que la solution approchée n'est pas valable pour notre problème, c-à-d le code en n'a pas encore [été] validé et je ne sais pas est ce que si le code que j'utilise est valable ou non ?
Le deuxième est de calculer l'erreur, j'ai traité un code mais le compilateur n'accepte pas ce code et je ne sais pas où le problème ? $$
erreur=||u-u_h||_{L^2(\Omega)^d}+||\mathrm{div}(u-u_h)||_{L^2(\Omega)}+||p-p_h||_{L^2(\Omega)}+|T-T_h|_{H^1(\Omega)}.
$$ Voici le code que j'ai traité et des fichiers pdf joints
macro Gradient(u) [dx(u), dy(u)] // macro Divergence(ux, uy) (dx(ux) + dy(uy)) // macro Lablace(T) (dxx(T))+dyy(T)))// macro GradientNorme(T) (dx(T)*dx(T)+dy(T)*dy(T))// //donners int d=2; real gm=1; real alpha=3.; real bitha=5; //solution exacte func phi=exp(-2*((x-1)^2+(y-2)^2)); func Pexct=gm*cos(pi*x/3)*cos(pi*y/3); func Texct=gm*x^2*(x-3)^2*(y-3)^2*y^2; func UXexct=gm*2*bitha*(1-x)*exp(-bitha*((x-1)^2+(y-1)^2)); func UYexct=gm*2*bitha*(1-y)*exp(-bitha*((x-1)^2+(y-1)^2)); //les fonctions donners func fx =-(pi/3)*sin((pi/3)*x)*cos((pi/3)*y)+(Texct+1)*UXexct; func fy =-(pi/3)*cos((pi/3)*x)*sin((pi/3)*y)+(Texct+1)*UYexct; func g=-alpha*(2*(x-3)^2*(y-3)^2*y^2+4*x*(x-3)*(y-3)^2*y^2 +2*x^2*(y-3)^2*y^2+4*x*(x-3)*(y-3)^2*y^2 +2*(x-3)^2*(y-3)^2*x^2+4*(x-3)^2*y*(y-3)*x^2 +4*x^2*(y-3)*(x-3)^2*y+2*(y)^2*(x-3)^2*x^2 +UXexct*(2*x^2*(x-3)*(y-3)^2*y^2+2*x*(x-3)^2*(y-3)^2*y^2) +UYexct*(2*x^2*(x-3)^2*(y-3)*y^2+2*x^2*(x-3)^2*(y-3)^2*y)); ////////////////// func Uexct=[UXexct,UYexct]; //Mesh real L2erreur; real[int] L2error(2); //for(int n=0 ;n<2 ;n++){ //mesh Th=square(40*(n+1),40*(n+1)); mesh Th = square(100,100,[3*x,3*y]); //Fespace fespace Uh(Th,RT0);//Raviart-Thomas Uh [ux, uy]; Uh [vx, vy]; Uh [FEUXexct,FEUYexct]=Uexct; fespace Ph(Th, P0); Ph p; Ph q; Ph FEPexct=Pexct; fespace Zh(Th, P1); Zh T; Zh S; Zh FETexct=Texct; //Function func nu=T+1; func uT=Gradient(T)'*[ux,uy]; //func GN=dx(FETexct-T)*dx(FETexct-T)+dy(FETexct-T)*dy(FETexct-T); //Problèm problem Vh1([ux, uy,T,p],[vx, vy,S,q]) =int2d(Th)(nu*[ux,uy]'*[vx,vy]) +int2d(Th)(-p*(dx(vx)+dy(vy))) -int2d(Th)( fx*vx + fy*vy ) -int2d(Th)(q*(dx(ux)+dy(uy))) -int2d(Th)(alpha*[dx(T),dy(T)]'*[dx(S),dy(S)]) -int2d(Th)(uT*S) +int2d(Th)(g*S) +on(1,2,3,4,T=0) ; L2erreur= (int2d(Th)((FEPexct-p)^2))^(0.5) +int2d(Th)(GradientNorme(FETexct-T)^2)^(0.5) +(int2d(Th)((FEUXexct-ux)^2+(FEUYexct-uy)^2))^(0.5) +(int2d(Th)(Divergence((UXexct-ux),(UYexct-uy))^2))^0.5 ; Vh1; //if(L2erreur<0.9){ plot(Th,wait=true, value=true ,cmm="Maillage de domaine"); plot(p,wait=true, value=true , cmm="Pression"); plot(FEPexct,wait=true, value=true , cmm="Pression exacte"); plot([ux, uy],wait=true, value=true , cmm="Vitesse"); plot([FEUXexct,FEUYexct],wait=true, value=true , cmm="Vitesse exacte"); plot(T,wait=true, value=true ,cmm="cheleur"); plot(FETexct,wait=true, value=true ,cmm="cheleur exacte"); // } //cout << " L2error " << " = "<< L2erreur<<endl; //}Merci cordialement.
Yassine OUZROUR.
Connectez-vous ou Inscrivez-vous pour répondre.