Implémentation TDMA en 2D en C (Mécanique des fluides)

[invite596fa921]

Bonjour,

Je souhaiterais implémenter la méthode TDMA(algo de Thomas) en 2D, mais ya un truc qui cloche dans mon programme :

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#define Jmax 11
#define Imax 12
int main()
{
	int i,it,l,k,j;
	//****les constantes
	double nu=1e-3;
	double	dy=0.02;
	//le fichier de sortie(pour gnuplot)
	FILE *f;
	//les tableaux
	double a[Imax],b[Imax],c[Imax],d[Imax];
	double Ue[Imax];
	double x[Imax];
	double dx[Imax];
	//***************
	double P[Imax];
	double Q[Imax];
		// on crée les tableaux à deux dimensions U et V
		double U[Imax][Jmax];
		double V[Imax][Jmax];
//*****************************************************
	for(i=1;i<Imax;i++){
		x[i]=i+0.5;
	}
	for(i=1;i<=Imax;i++){
		dx[i]=x[i+1]-x[i]/2;
	}
	//
	for(i=1;i<Imax;i++)
	{
		Ue[i]=sin(x[i]);
	}
	//
	//Conditions aux limites***********************************************************
	//bord gauche
	for(i=1;i<=Imax;i++)
	{
		U[i][1]=0.000;
		V[i][1]=0.000;
	}
	//bord droite
	for(i=1;i<=Imax;i++)
	{
		U[i][Jmax]=Ue[i];
		V[i][Jmax]=0.0;
	}
	//bord extérieur(en haut du maillage)
	for(j=1;j<Jmax;j++)
	{
		U[1][j]=0.000; 
		V[1][j]=0.000;
	}
	//*******************************************************************************************
	//**************************Application : Méthode TDMA***************************************
	//*******************************************************************************************
	//*******sur les lignes Nord-Sud ************************************************************
	//*******************************************************************************************
	P[1]=0.0;
	for(i=2;i<Imax;i++)
	{	
		Q[1]=U[i][1];   		//tous nuls(Conditions aux limites gauche)
		for(j=2;j<Jmax;j++)
		{
			//calcul des coef. ligne par ligne
			a[j]=(V[i][j]/2*dy)-(nu/(dy*dy));
			b[j]=-(V[i][j]/2*dy)-(nu/(dy*dy));
			d[j]=(U[i][j]/dx[i])+(2*nu/(dy*dy));
			c[j]=(U[i][j]*U[i][j])/dx[i]+(Ue[i]*(Ue[i+1]-Ue[i])/dx[i]);
			//**********************************************************
			//on Calcule P(.,j) et Q(.,j) pour chaque ligne de j=2:Jmax
			//**********************************************************
			P[j]=-a[j]/(d[j]+b[j]*P[j-1]);
			Q[j]=(c[j]-b[j]*Q[j-1])/(d[j]+b[j]*P[j-1]);
			//**********************************************************
			//puis on calcule  U(.,j) pour chaque ligne de j=Jmax-1 à 1
			//**********************************************************
			do{
				U[i][j]=P[j]*U[i][j-1]+ Q[j];
				j=j-1;
			}while(j!=1);
		}
		//**************************************************************
		//une fois U(.,j) obtenu, on calcule V(.j) j=1:Jmax-1
		//**************************************************************
		for(j=1;j<Jmax;j++)
		{	V[i][j]=V[i][j-1]-(dy/dx[i])*(U[i][j-1]-U[i][j])+(dy/2.*dx[i])*(U[i][j-1]-U[i][j]);}
	}//*********************fin Nord-Sud*********************************************************
	//*******************************************************************************************
	//********************sur les colonnes Est-Ouest*********************************************
		P[1]=0.0;
	for(j=2;i<Jmax;j++)
	{	
		Q[1]=U[1][j];   		//tous nuls(Conditions aux limites gauche)
		for(i=2;i<Imax;i++)
		{
			//calcul des coef. colonne par colonne
			a[i]=(V[i][j]/2*dy)-(nu/(dy*dy));
			b[i]=-(V[i][j]/2*dy)-(nu/(dy*dy));
			d[i]=(U[i][j]/dx[j])+(2*nu/(dy*dy));
			c[i]=(U[i][j]*U[i][j])/dx[j]+(Ue[j]*(Ue[j+1]-Ue[j])/dx[j]);
			//**********************************************************
			//on Calcule P(i,.) et Q(i,.) pour chaque colonne de i=2:Imax-1
			//**********************************************************
			P[i]=-a[i]/(d[i]+b[i]*P[i-1]);
			Q[i]=(c[i]-b[i]*Q[i-1])/(d[i]+b[i]*P[i-1]);
			//**********************************************************
			//puis on calcule  U(i,.) pour chaque colonne de i=Imax-1:2
			//**********************************************************
			do{
				U[i][j]=P[i]*U[i-1][j]+ Q[i];
				printf("U[%d][%d]\n",i,j,U[i][j]);
				i=i-1;
			}while(i!=2);
		}
		//**************************************************************
		//une fois U(i,.) obtenu, on calcule V(i,.) i=1:Imax-1
		//**************************************************************
		for(i=1;i<Imax;i++)
		{	V[i][j]=V[i-1][j]-(dy/dx[j])*(U[i-1][j]-U[i][j])+(dy/2.*dx[j])*(U[i-1][j]-U[i][j]);
			printf("V[%d][%d]\n",i,j,V[i][j]);
		}
	}//*********************fin Est-Ouest********************************************************
	//************************************************************************************************************
	//**************************FIN TDMA-2D(U,V)******************************************************************
	//************************************************************************************************************
	//******************écriture des données pour gnuplot****************************************
	//**********************je veux simuler le profil de vitesse(i.e U/Ue en fonction x/c=0.5;1;1.5;2;2.5,...)*****
			f=fopen("meca-result.dat","w");
			for(i=1;i<=Imax;i++)
			{
				for(j=1;j<=Jmax;j++)
				{
					fprintf(f,"%lf %lf \n",U[i][j]/Ue[i],x[i]);
				}
			}
			fclose(f);
	return 1;
}

-----

[invite1c6b0acc]

Bonjour,

Si vraiment tu veux de l'aide, ce serait aussi bien de ne pas faire tout ton possible pour compliquer les choses ...
=> utilise de balises CODE : là, sans les indentations c'est carrément illisible.
=> "y a un truc qui cloche" : c'est vraiment l'indication la plus précise que tu puisses donner ? tu ne sais même pas si ton programme se compile ? s'il plante à l'exécution ? s'il fournit des résultats faux ? ou mal formatés ? s'il provoque des fuites radio-actives ? Parce que moi (et tous les autres lecteurs) j'en suis réduit à essayer de le deviner ...
=> tu sais qu'en langage C, on a le droit de donner aux variables un nom qui a un rapport avec leur signification ? Parce que a,b,c,d,i,j,k,l,x etc. c'est carrément cryptique ...

[invite1c6b0acc]

Bon, mon coup de colère passé, et comme mon compilateur n'a pas d'états d'âmes, il m'a suffit d'essayer de faire tourner ton programme sous debugger pour y trouver une boucle sans fin :

for(j=2;j<Jmax;j++) // j'incrémente j (un pas en avant)
   {
    //calcul des coef. ligne par ligne
    a[j]=(V[i][j]/2*dy)-(nu/(dy*dy));
    b[j]=-(V[i][j]/2*dy)-(nu/(dy*dy));
    d[j]=(U[i][j]/dx[i])+(2*nu/(dy*dy));
    c[j]=(U[i][j]*U[i][j])/dx[i]+(Ue[i]*(Ue[i+1]-Ue[i])/dx[i]);
    //****************************** ****************************
    //on Calcule P(.,j) et Q(.,j) pour chaque ligne de j=2:Jmax
    //****************************** ****************************
    P[j]=-a[j]/(d[j]+b[j]*P[j-1]);
    Q[j]=(c[j]-b[j]*Q[j-1])/(d[j]+b[j]*P[j-1]);
    //****************************** ****************************
    //puis on calcule U(.,j) pour chaque ligne de j=Jmax-1 à 1
    //****************************** ****************************
    do
       {
         U[i][j]=P[j]*U[i][j-1]+ Q[j];
         j=j-1; // je fais redescendre j à 1 (un pas en arrière)
       }while(j!=1);
   } // tant que j ne vaut pas 11 : On n'est pas arrivé ...

Tu as un debugger ?

[invite596fa921]

Bonjour,
mes excuses pour la présentation, je suis nouveau et j'ai pas pris le temps de jeter un coup d'oeil sur le réglement(j'ai eu tort!).

en fait, c'est un système d'équations en 2D avec deux composantes(U,V),,
et t'as pas tort, j'ai du mal copier-coller(la boucle for sur j dois être fermée avant l'entrée du boucle dou--while,,,

i pour parcourir l'axe horizontal, j pour l'axe vertical

le système à résoudre :

$$\left\{ 
\begin{array}{cc}\label{m6}
\hspace{0.5cm}U_i^{j+1}=P_iU_{i+1}^{j+1}+Q_i\quad i=2:Imax-1,\dots,j=2:Jmax-1\vspace{0.3cm}\\
V_i^{j+1}=V_{i-1}^{j+1}-e_i(U_{i-1}^{j+1}-U_i^{j+1})+d_j\quad\quad i=1:Imax,\dots,j=1:Jmax
\end{array}\right.$$
Avec 
$$\left\{ 
\begin{array}{cc}\label{m6}
P_i=-A_i\left(D_i+B_iP_{i-1}\right)^{-1}\vspace{0.3cm}\nonumber\\
Q_i=\left(C_i-B_iQ_{i-1}\right)\left(D_i+B_iP_{i-1}\right)^{-1}\nonumber\vspace{0.3cm}\\
B_i=-\frac{V_i^j}{2\Delta y}-\frac{\nu}{(\Delta y)^2}\nonumber\vspace{0.3cm}\\
D_i=\frac{U_i^j}{\Delta x}+\frac{2\nu}{(\Delta y)^2}\nonumber\vspace{0.3cm}\\
A_i=\frac{V_i^j}{2\Delta y}-\frac{\nu}{(\Delta y)^2}\nonumber\vspace{0.3cm}\\
C_i=\frac{(U_i^j)^2}{\Delta x}+U_e^j\frac{(U_e^{j+1}-U_e^j)}{\Delta x}\nonumber\vspace{0.3cm}\\
d_j=\frac{\Delta y}{2\Delta x}(U_{i-1}^j+U_i^j)\nonumber\vspace{0.3cm}\\
e_i=\frac{\Delta y}{\Delta x}\nonumber\vspace{0.3cm}\\
P_1=Q_1=0\quad\quad P_{Imax}=0\quad Q_{Imax}=U_e\quad U_e<0.99) \nonumber\vspace{0.3cm}\\
V_1^{j+1}=0,\quad V_{Imax}^{j+1}=0\nonumber\vspace{0.3cm}\\
\end{array}\right.$$
Pour l'instant,j'ai choisi $x=0.5,1,1.5,\dots$, et $dx=x_{i+1}-xi/2$, $Ue=sin(x)$ $dy=0.02$ $U_1=0$ et $U_N=Ue$, $nu=1e-3$

le principe de l'algo :

-Calculer (P_2,Q_2),(P_3,Q_3),...(P_{N-1},Q_{N-1})
_puis U_{N-1},U_{N-2},...,U_2
_puis V_1,...,V_{N-1}(à partir de U)

*********************************************************************************************
rectification:(sur la ligne Nord-Sud)
//*******sur les lignes Nord-Sud ************************************************************
	//*******************************************************************************************
	//*******sur les lignes Nord-Sud ************************************************************
	//*******************************************************************************************

	P[1]=0.0;
	for(i=2;i<Imax;i++)
	{	
		Q[1]=U[i][1];   		//tous nuls(Conditions aux limites gauche)
		for(j=2;j<Jmax;j++)
		{
			//calcul des coef. ligne par ligne
			a[j]=(V[i][j]/2*dy)-(nu/(dy*dy));
			b[j]=-(V[i][j]/2*dy)-(nu/(dy*dy));
			d[j]=(U[i][j]/dx[i])+(2*nu/(dy*dy));
			c[j]=(U[i][j]*U[i][j])/dx[i]+(Ue[i]*(Ue[i+1]-Ue[i])/dx[i]);
			//
			//**********************************************************
			//on Calcule P(.,j) et Q(.,j) pour chaque ligne de j=2:Jmax
			//**********************************************************
			P[j]=-a[j]/(d[j]+b[j]*P[j-1]);
			Q[j]=(c[j]-b[j]*Q[j-1])/(d[j]+b[j]*P[j-1]);
			//**********************************************************
			//puis on calcule  U(.,j) pour chaque ligne de j=Jmax-1 à 1
			//**********************************************************
			k=Jmax-1 ;
			do{
				U[i][k]=P[j]*U[i][k+1]+ Q[j];
				k=k-1;
			}while(k!=2);
		//**************************************************************
		//une fois U(.,j) obtenu, on calcule V(.j) j=1:Jmax-1
		//**************************************************************
		for(l=1;l<Jmax;l++)
		{	V[i][l]=V[i][l-1]-(dy/dx[i])*(U[i][l-1]-U[i][l])+(dy/2.*dx[i])*(U[i][l-1]-U[i][l]);}
	}
}//*********************fin Nord-Sud************
	//*******************************************************************************************

[A voir en vidéo sur Futura]

[invite596fa921]

ci-joint le sujet en pdf ...