Darcy et équation de chaleur — Les-mathematiques.net The most powerful custom community solution in the world

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
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.105658
Connectez-vous ou Inscrivez-vous pour répondre.
Success message!