Problème de Darcy couplée avec la chaleur

yassine azrour

j'ai essayé de programmer le code Freefem de la première formulation variationnelle qui est le suivante. $$

(V_{h,1})\left\{\begin{aligned}
\text { Chercher }\left(\mathbf{u}_{h}, p_{h}, T_{h}\right) \in \mathcal{W}_{h, 1} \times M_{h, 1} \times X_{h} \text { tel que : } \quad \quad \quad \quad \quad \quad \\
\qquad \forall \mathbf{v}_{h} \in \mathcal{W}_{h, 1}, \quad \int_{\Omega} \nu\left(T_{h}(\mathbf{x})\right) \mathbf{u}_{h}(\mathbf{x}) \cdot \mathbf{v}_{h}(\mathbf{x}) d \mathbf{x}-\int_{\Omega} p_{h}(\mathbf{x}) \operatorname{div} \mathbf{v}_{h}(\mathbf{x}) d \mathbf{x}= \\
\int_{\Omega} \mathbf{f}(\mathbf{x}) \cdot \mathbf{v}_{h}(\mathbf{x}) d \mathbf{x}, \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \\
\forall \mathbf{q}_{h} \in M_{h, 1}\text{, } \quad \int_{\Omega} q_{h}(\mathbf{x}) \operatorname{div} \mathbf{u}_{h}(\mathbf{x}) d \mathbf{x}=0,
\quad \quad \quad \quad \quad \quad \quad \quad \quad \\
\forall S_{h} \in X_{h}, \alpha \int_{\Omega} \nabla T_{h}(\mathbf{x}) \cdot \nabla S_{h}(\mathbf{x}) d \mathbf{x}+\int_{\Omega}\left(\mathbf{u}_{h} \cdot \nabla T_{h}\right)(\mathbf{x}) S_{h}(\mathbf{x}) d \mathbf{x}= \\
\int_{\Omega} g(\mathbf{x}) S_{h}(\mathbf{x}) d \mathbf{x}. \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad
\end{aligned}
\right.

$$ 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 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}+||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 ces joints des fichiers pdf
Merci cordialement.
Yassine OUZROUR.

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;
//}