Problème avec Freefem++

Yann
Yann
dans Logiciels pour les mathématiques
Bonjour,
J'ai un petit problème avec Freefem++.
J'essaie de programmer une équation d'advection-diffusion-réaction, avec coefficients de diffusion et de réaction variables avec le temps, et j'obtiens un code d'erreur à la ligne 79 que j'ai mise en rouge. Je ne vois pas où est le problème. Si vous avez une piste, je suis preneur.
Merci
Yann 
//Equation de diffusion avec variation du coefficient de diffusion D et d'advection par intervalles de temps

// saisie des données et déclaration des variables
int N=25;
real x0=3.25,y0=12.5,T=25,dt=T/N,dx=0.5,dy=0.5,eps=1/0.2,R0=4.30;
func u0x=exp((x0+dx-x)*eps)/(exp(-(x0+dx-x)*eps)+exp((x0+dx-x)*eps))-exp((x0-dx-x)*eps)/(exp(-(x0-dx-x)*eps+exp((x0-dx-x)*eps)));
func u0y=exp((y0+dy-y)*eps)/(exp(-(y0+dy-y)*eps)+exp((y0+dy-y)*eps))-exp((y0-dy-y)*eps)/(exp(-(y0-dy-y)*eps+exp((y0-dy-y)*eps)));
func u0=4*u0x*u0y;
real alphadv,muadv,alphdiff,mudiff;
//func advx0=0;
//func advy0=-0.21*exp(-(-x0)^2/2)/sqrt(2*pi);
//func diff0=0*exp(-(-x0-0.23)^2/2)/sqrt(2*pi);

//maillage et choix des éléments finis
border gamma1(t=0,25){x = t ; y = 0;}
border gamma2(t=0,25){x = 25 ; y =  t;}
border gamma3(t=0,25){x = 25-t ; y = 25;}
border gamma4(t=0,25){x = 0 ; y = 25-t;}

//on définit les problèmes variationnels : advection, diffusion
//problem advection(u,w)=
  //   int2d(Th)(u*w/dt +(dx(u)*vx+dy(u)*vy)*w +u*(dx(vx)+dy(vy))*w)
    // -int2d(Th)(uold*w/dt)
     //+on(gamma1,gamma2,gamma3,gamma4,u=0); 				//condition de Dirichlet homogène	

//problem diffusion(m,p)=
  //   int2d(Th)(m*p/dt +D*( dx(m)*dx(p)+dy(m)*dy(p)))
    // -int2d(Th)(mold*p/dt)
    // +on(gamma1,gamma2,gamma3,gamma4,m=0); 				//condition de Dirichlet homogène	
//ofstream ff("ADR.dat");

//construction du nuage de points
// (t,u(t,x)) où x=(0.5,0.5)
for(int t=1;t<=T;t=t+dt)
{	
	//uold=u;	
	if (t==1){alphadv=0.21;muadv=0;alphdiff=0;mudiff=0.23;R0=4.3;}
	else if (t==2){alphadv=0.21;muadv=0;alphdiff=0;mudiff=0.23;R0=1.57;}
	else if (t==3){alphadv=0.21;muadv=0;alphdiff=0;mudiff=0.23;R0=1.01;}
	else if (t==4){alphadv=1.06;muadv=3.24;alphdiff=0;mudiff=0.23;R0=0.57;}
	else if (t==5){alphadv=1.06;muadv=3.24;alphdiff=0;mudiff=0.23;R0=0.37;}
	else if (t==6){alphadv=1.06;muadv=3.24;alphdiff=0;mudiff=0.230;R0=0.31;}
	else if (t==7){alphadv=0.73;muadv=2.73;alphdiff=0.18;mudiff=3.49;R0=0.29;}
	else if (t==8){alphadv=0.73;muadv=2.73;alphdiff=0.18;mudiff=3.49;R0=0.15;}
	else if (t==9){alphadv=0.73;muadv=2.73;alphdiff=0.18;mudiff=3.49;R0=0.45;}
	else if (t==10){alphadv=0.73;muadv=2.79;alphdiff=0.62;mudiff=3.66;R0=0.39;}
	else if (t==11){alphadv=0.73;muadv=2.79;alphdiff=0.62;mudiff=3.66;R0=0.44;}
	else if (t==12){alphadv=0.73;muadv=2.79;alphdiff=0.62;mudiff=3.66;R0=0.34;}
	else if (t==13){alphadv=1.08;muadv=3.15;alphdiff=1.96;mudiff=3.94;R0=0.39;}
	else if (t==14){alphadv=1.08;muadv=3.15;alphdiff=1.96;mudiff=3.94;R0=0.35;}
	else if (t==15){alphadv=1.08;muadv=3.15;alphdiff=1.96;mudiff=3.94;R0=0.28;}
	else if (t==16){alphadv=1.12;muadv=4.25;alphdiff=1.95;mudiff=3.89;R0=0.28;}
	else if (t<=18){alphadv=1.12;muadv=4.25;alphdiff=1.95;mudiff=3.89;R0=0.27;}
	else if (t==19){alphadv=0.51;muadv=5.17;alphdiff=0.70;mudiff=4.83;R0=0.31;}
	else if (t==20){alphadv=0.51;muadv=5.17;alphdiff=0.70;mudiff=4.83;R0=0.30;}	
	else if (t==21){alphadv=0.51;muadv=5.17;alphdiff=0.70;mudiff=4.83;R0=0.31;}
	else if (t==22){alphadv=0.57;muadv=5.42;alphdiff=0.86;mudiff=4.11;R0=0.33;}
	else {alphadv=0.57;muadv=5.42;alphdiff=0.86;mudiff=4.11;R0=0.30;}
	func advx=0;
	func advy=-alphadv*exp(-((x-x0-muadv)^2+(y-y0-muadv)^2)/2)/sqrt(2*3.14159265);
	func diff=alphdiff*exp(-((x-x0-mudiff)^2+(y-y0-mudiff)^2)/2)/sqrt(2*3.14159265);

	mesh Th = buildmesh (gamma1(60)+gamma2(60)+gamma3(60)+gamma4(60));
	//fespace Uh(Th,P1);
	fespace Vh(Th,P2);
	//fespace Wh(Th,P2);
	//Uh n,q,R=R0,nold;
	Vh u=u0,w,uold,vx=advx,vy=advy,D=diff,R=R0;
	//Wh m=u0,p,mold,D=diff;

	uold=u;
	solve adr(u,w)=
     	[color=#FF0000]int2d(Th)(u*w/dt +(dx(u)*vx+dy(u)*vy)*w[/color]
		+u*(dx(vx)+dy(vy))*w +D*( dx(u)*dx(w)+dy(u)*dy(w))
		+R*u*w)
    	-int2d(Th)(uold*w/dt)
     	+on(gamma1,gamma2,gamma3,gamma4,u=0); 				//condition de Dirichlet homogène		
}

Réponses

  • remarque
    remarque
    Sans la version de FreeFem, le message d'erreur et un fichier un peu plus minimal qui engendre cette même erreur, difficile de se faire une idée.
  • Yann
    Yann
    Bonne Remarque ;-)

    La version de Freefem est 14.3
    Je ne peux joindre le programme car fichier non accepté.
    Et le soucis est justement qu'avec un programme plus simple (dans les conditions) mais "identique", je n'avais pas de message d'erreur.
  • Yann
    Yann
    Exemple : problème d'advection-diffusion... (avec un calcul d'erreur, mais bon...) 
    //Equation de diffusion avec coefficient de diffusion D constant

// saisie des données et déclaration des variables
verbosity=0;
int N=50;
real D0=0.01, T=10., dt=T/N; 	// D0 coefficient de diffusion 
func vx0=0.;
func vy0=-0.01; 		// Direction de croissance verticale (advection)
real[int] abscisse(N);
real[int] L2error(N);

//maillage et choix des éléments finis : carré [-1;1]*[-1;1] & éléments P1
border gamma1(t=-1,1){x = t ; y = -1;}
border gamma2(t=-1,1){x = 1 ; y =  t;}
border gamma3(t=-1,1){x = -t ; y = 1;}
border gamma4(t=-1,1){x = -1 ; y = -t;}

//on définit les problème variationnels : advection, diffusion	
			//condition de Dirichlet homogène	
ofstream ff("ADR.dat");

for(real t=0;t<T;t=t+dt)
{
	func nexact=(4/(4+0.01*t))*exp(-(x^2+(y+0.01*t)^2)/(16+0.04*t));
	func n0=exp((-x^2-y^2)/16);
	mesh Th = buildmesh (gamma1(60)+gamma2(60)+gamma3(60)+gamma4(60));
	fespace Vh(Th,P1);
	fespace Wh(Th,P2);
	Vh n=n0,w,nold,vx=vx0,vy=vy0;
	Wh m,p,mold,D=D0;

	nold=n;
	solve  advection(n,w)=
     int2d(Th)(n*w/dt +(dx(n)*vx+dy(n)*vy)*w)
     -int2d(Th)(nold*w/dt)
     +on(gamma1,gamma2,gamma3,gamma4,n=1); 	
	
	mold=n;
	solve diffusion(m,p)=
     int2d(Th)(m*p/dt +D*(dx(m)*dx(p)+dy(m)*dy(p)))
     -int2d(Th)(mold*p/dt)
     +on(gamma1,gamma2,gamma3,gamma4,m=1); 	;	

	plot(n,value=1,fill=1);
	abscisse[t]=t;
	L2error[t]= sqrt(int2d(Th)((n-nexact)^2));
	cout <<"Erreur L2 "<<t <<" = "<<L2error[t]<<endl;
}
  • remarque
    remarque
    Ok, mais quel est le message d'erreur (on n'a pas forcément envie de faire tourner FreeFem++ rien que pour le voir) ?
  • Yann
    Yann
    Le voici : 
    Compile error : 
	line number :71, )
error Compile error : 
	line number :71, )
 code = 1 mpirank: 0
FreeFem++ returned error 1
  • remarque
    remarque
    Ah oui, ce n'est pas très parlant... :-D C'est laquelle la ligne 71 ? (En comptant à l'œil nu, il n'y en a qu'une cinquantaine dans le code que tu as posté).
  • Yann
    Yann
    C'est bon, j'ai trouvé ! Ce sont les variables dx et dy (pas en espace) qui ne conviennent pas au début du programme.
    En effet, si u est une fonction définie dans un espace variationnel Vh, dx(u) représente la dérivée de u par rapport à x.
    Donc nom réservé. Con, mais chronophage pour trouver l'erreur !

    Merci à remarque de m'avoir répondu.

    Yann
  • remarque
    remarque
    Ah oui, effectivement. C'est un peu dommage que le message d'erreur ne soit pas plus explicite. Je me demande s'il n'y a pas un réglage de verbosité qui permettrait d'en savoir plus...
Connectez-vous ou Inscrivez-vous pour répondre.

Bonjour!

Pour participer au forum, cliquer sur l'un des boutons :

Connexion S'inscrire

Catégories

In this Discussion

Qui est en ligne 33

APfAmadouBibixCidrolinGravityMath CossThierry PomaVassilliaabcdbiguine_equationgai requin +22 Invités