# Stability problem related to SUPG

**URL:** <https://community.freefem.org/t/stability-problem-related-to-supg/3382>\
**Category:** General Discussion\
**Created:** [July 4, 2024, 2:54pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382 "2024-07-04T14:54:18Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Walid](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/walid/32/2477_2.png) [@Walid](https://community.freefem.org/u/Walid)\
**Post date:** [July 4, 2024, 2:54pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/1 "2024-07-04T14:54:18Z")

</div>

Hello All,

I hope you are fine.  
I am trying to solve 4 coupled PDE as below  
∇.∇V=−q(Np-Ne-Nn)/ϵ  
∂Ne/∂t +∇⋅(Ne.(-We)-De.∇Ne) = Se  
∂Np/∂t +∇⋅(Np.(Wp)-Dp.∇Np) = Sp  
∂Nn/∂t +∇⋅(Nn.(-Wn)-Dn.∇Nn) = Sn

I have tried to model the system but I faced instability I tried artificial diffusion and SUPG but also I face instability and divergence. I need your kind support please  
here is the code with SUPG trail.

```auto
load "iovtk";
verbosity = 1; 
//load "UMFPACK64";
				
real Uapp=-5500; // Applied voltage
real q=1.6022e-19; //[C] electron charge
real eps0=8.85418782e-12; // permittivity of free sace

real mup = 2.24e-4;//[m^2/(V*s)]
real mun = 2.16e-4;//[m^2/(V*s)]
real De = 0.18 ;//[m^2/s]
real Dp = 0.028e-4;//[m^2/s]
real Dn = 0.043e-4;//[m^2/s]
real beta=2e-13;//[m^3/s] 
real gamma= 0.01;// secondary emission coefficient
real w1 = -1.65e7; //[V/m]
real w2 = 3.5e5; //1/m
real w3 = 1500;// [1/m]
real w4	= -2.5e6;// [V/m]
real w5 = -60.6;// [m/s]
real w6	= 2.43e-4;// [m²/(V·s)]
real w7	= -2.7e-4; //[m²/(V·s)]
real w8 = 1.9163;// [m^2/(V*s)]

real dt = 5e-13; // Time step size
real T = 2e-6; // Final Time 
int numsteps = T / dt; // Number of time steps
real tt=0.;

real I;
real[int] It(numsteps);

int GND =1, needle =2, ax = 3, out1 = 4;// boundary numbering
real x0=0.03; // x-coordinate for domain 
real y0=(35e-6*((29)^2)+6e-3);// y-coordinate for domain 
real h=6e-3; // needle position  
real yy1=0.86; //needle split point
real yy2=2.3; //needle split point
real xn1=(2*35e-6*yy1); //x-coordinate for needle split point 
real yn1= (35e-6*(yy1^2)+6e-3); //y-coordinate for needle split point 
real xn2=(2*35e-6*yy2); //x-coordinate for needle split point 
real yn2= (35e-6*(yy2^2)+6e-3); //y-coordinate for needle split point 
real w = 2*35e-6*29;

border b1 (t=0. , yy1){x=2*35e-6*t; y=(35e-6*(t^2)+6e-3); label=needle;};
border b2 (t=yy1 , yy2){x=2*35e-6*t; y=(35e-6*(t^2)+6e-3); label=needle;};
border b3 (t=yy2, 29 ){x=2*35e-6*t; y=(35e-6*(t^2)+6e-3); label=needle;};

border b4 (t=5.8e-3, 6e-3 ){x=0; y=t; label=ax;}; 
border b5 (t=5.6e-3, 5.8e-3){x=0; y=t; label=ax;}; 
border b6 (t=0 , 5.6e-3){x=0; y=t; label=ax;}; 

border b7 (t=5.8e-3, yn1){x=xn1; y=t;};
border b8 (t=5.6e-3, yn2){x=xn2; y=t;};

border b9 (t=0.,x0){x=t; y=0.; label=GND;};
border b10(t=0, y0){x=x0; y=t; label=out1;};
border b11(t=w,x0){x=t; y=y0; label=out1;};

mesh Th;
Th = buildmesh (
				  b1 (-150) //needle
				+ b2 (-150) //needle
				+ b3 (-50) //needle
				+ b4 (-200)	//ax
				+ b5 (-200)	//ax
				+ b6 (-2000)	//ax
				+ b7 (150)	//mesh control edge
				+ b8 (300)	//mesh control edge
				+ b9 (50)	//GND
				+ b10(30)	//out1
				+ b11(-30) //out2
				);

//finite element space 
fespace VE(Th, P1); //
fespace VNe(Th, P1); //

// Define trial functions for electric potential and concentrations
VE u, v, E1, E2, Emag, E1L, E2L, EL, f, rohv;
VNe Ne, Neold, Se, S1, S2, S3, S4, We, Wex, Wey, alpha, eta, mue, tau1, dtf, Dstab, Np, Npold, Sp, Wp, Wpy, Wpx, tau2, Nn, Nnold, Sn, Wnx, Wn, Wny, tau3, v1, v2, v3;

macro grad(u) [dx(u), dy(u)] //

// initial condition for e, p, n distribution
func a = (x^2); func b = (y-h)^2; real c=2*((25e-6)^2);
real Nm = 1e16; //[m^-3] maximum density
Neold = (Nm)*exp(-(x^2/(c))-(y-0.006)^2/(c)); 
Npold = (Nm)*exp(-(x^2/(c))-(y-0.006)^2/(c)); 
Nnold = 0; 

problem Poisson(u,v) = int2d(Th)(grad(u)'*grad(v))
						-int2d(Th)(f*v)
						+on (1,u = 0)
						+on (2,u = Uapp)
						;
//matrix AP;
//real [int] bP(VE.ndof);						

//solve to get laplacian Electric field 		
	f=0;
//	AP = Poisson(VE,VE);
//	bP = Poisson(0 ,VE);	
//	u[] = AP^-1*bP ;
	Poisson;
	E1L = -dx(u);
	E2L = -dy(u);
	EL = sqrt(E1L*E1L + E2L*E2L);		
//	plot(EL, cmm="Laplacian electric field Mag",fill=true);//, hsv=colorhsv);
	
problem NeL(Ne, v1) =	int2d(Th)(Ne*v1/dt	
									+(De*(grad(Ne)'*grad(v1))
									+(-Wex * dx(Ne) + -Wey * dy(Ne))*v1))
						-int2d(Th)(Se*v1)
						-int2d(Th)(Neold*v1/dt)
						+int1d(Th,needle)(gamma*Npold*mup*Emag*v1)	
						+int2d(Th)((tau1*(-Wex*dx(v1)+-Wey*dy(v1)))*(Ne/dt-De*(dxx(Ne)+dyy(Ne))+-Wex*dx(Ne)+-Wey*dy(Ne))) //SUPG 
						-int2d(Th)((tau1*(-Wex*dx(v1)+-Wey*dy(v1)))*(Neold/dt)) //SUPG 				
						-int2d(Th)((tau1*(-Wex*dx(v1)+-Wey*dy(v1)))*(Se)) //SUPG 					
						;

problem NpL(Np, v2) =	int2d(Th)(Np*v2/dt
									+(Dp*(grad(Np)'*grad(v2))
									+(Wpx * dx(Np) + Wpy * dy(Np))*v2))		
						-int2d(Th)(Sp*v2)
						-int2d(Th)(Npold*v2/dt)
						+on (GND, Np = 0)	
						+int2d(Th)((tau2*( Wpx*dx(v2)+ Wpy*dy(v2)))*(Np/dt-Dp*(dxx(Np)+dyy(Np))+ Wpx*dx(Np)+ Wpy*dy(Np))) //SUPG 
						-int2d(Th)((tau2*( Wpx*dx(v2)+ Wpy*dy(v2)))*(Npold/dt)) //SUPG 
						-int2d(Th)((tau2*( Wpx*dx(v2)+ Wpy*dy(v2)))*(Sp)) //SUPG 
						;	
					
problem NnL(Nn, v3) =	int2d(Th)(Nn*v3/dt	
									+(Dn*(grad(Nn)'*grad(v3))
									+(-Wnx * dx(Nn) + -Wny * dy(Nn))*v3))
						-int2d(Th)(Sn*v3)
						-int2d(Th)(Nnold*v3/dt)
						+on (needle, Nn = 0)			
						+int2d(Th)((tau3*(-Wnx*dx(v3)+-Wny*dy(v3)))*(Nn/dt-Dn*(dxx(Nn)+dyy(Nn))+-Wnx*dx(Nn)+-Wny*dy(Nn))) //SUPG 
						-int2d(Th)((tau3*(-Wnx*dx(v3)+-Wny*dy(v3)))*(Nnold/dt)) //SUPG
						-int2d(Th)((tau3*(-Wnx*dx(v3)+-Wny*dy(v3)))*(Sn)) //SUPG 
						;	
						
//matrix AE, APo, An; 
//createMat(Th, AE, P1);
//createMat(Th, APo, P1);
//createMat(Th, An, P1);
//real [int] bE(VNe.ndof), bPo(VNp.ndof), bn(VNn.ndof);

real alphas=1;	

VNe hr1=hTriangle;
cout << hr1[].min << endl;
cout << hr1[].max << endl;

//int step=0; 
//while (tt<T){

for (int step = 0; step < numsteps; ++step){
	
	cout<< "time:" << tt << endl;
	cout<< "step #:" << step << endl;

	rohv = q*(Npold-Neold-Nnold);      
	f = -rohv/eps0; 	
	
//	AP=Poisson(VE,VE);
//	bP=Poisson(0,VE);
//	u[] = AP^-1*bP ;
//	cout << "Volt Max " << u[].max << endl ;
//	cout << "Volt Min " << u[].min << endl << endl ;
	Poisson;
	E1 = -dx(u); E2 = -dy(u); Emag = sqrt(E1*E1 + E2*E2);

	mue = w8/((Emag)^0.25);
	alpha = w2*exp(w1/Emag);
	eta =	w3*exp(w4/Emag);
	We = mue*Emag;
	Wp = mup*Emag;
	Wn = mun*Emag;
	Wex = E1*mue; 
	Wey = E2*mue;
	Wpx = E1*mup; 
	Wpy = E2*mup;
	Wnx = E1*mun;
	Wny = E2*mun;
	
	tau1= alphas*hr1/2/sqrt(Wex^2+Wey^2);
	tau2= alphas*hr1/2/sqrt(Wpx^2+Wpy^2);
	tau3= alphas*hr1/2/sqrt(Wnx^2+Wny^2);
	
//	dtf = hr/We;
//	dt = dtf[].min;
//	cout<< "Time step #:" << dt << endl;
//	cout<< "Calculated Time step #:" << dtf[].min << endl;

	S1	= Neold*alpha*abs(We); // [1/m^3.s] Ionizatoin 
	S2	= Neold*eta*abs(We); // [1/m^3.s] Attachment
	S3	=	Neold*Npold*beta; // [1/m^3.s] Recombination between electrons and positive ions
	S4	=	Nnold*Npold*beta; // [1/m^3.s] Recombination between positive ions and negative ions
  
	Se = S1-S2-S3;
	Sp = S1-S3-S4;
	Sn = S2-S4;
//	Dstab = De+(0.25*hr*We);
	
	NpL;
//	AE=NeL(VNe,VNe);
//	bE=NeR(0,VNe);
//	Ne[] = AE^-1*bE ;
//	cout << "Electron Max " << Ne[].max << endl ;
//	cout << "Electron Min " << Ne[].min << endl << endl ;

	NeL;
//	APo=NpL(VNp,VNp);
//	bPo=NeR(0,VNp);
//	Np[] = APo^-1*bPo ;
//	cout << "Positive Max " << Np[].max << endl ;
//	cout << "Positive Min " << Np[].min << endl << endl ;
	
	NnL;
//	An=NnL(VNn,VNn);
//	bn=NnR(0,VNn);
//	Nn[] = An^-1*bn ;
//	cout << "Negative Max " << Nn[].max << endl ;
//	cout << "Negative Min " << Nn[].min << endl << endl ;

	//Peclet Number
//	VNe hK = hTriangle ;
	VNe PeK = We*hr1/2./De ;
	cout << "Pe max= " << PeK[].max << endl ;
	cout << "Pe min= " << PeK[].min << endl ;
		
		if (step % 50 == 0)
		{
		savevtk("SUPG 3SP dt5e-13-Newmesh2/TNS "+step+".vtu", Th, Emag, Ne, Np, Nn, PeK);
		}
	
		for (int i=0; i<VNe.ndof; ++i){if(Ne[][i]<0.)Ne[][i]=0.;}
		for (int i=0; i<VNe.ndof; ++i){if(Np[][i]<0.)Np[][i]=0.;}
		for (int i=0; i<VNe.ndof; ++i){if(Nn[][i]<0.)Nn[][i]=0.;}
	
	I = int2d(Th)(2*pi*x*q*(mue*Ne + mup*Np + mun*Nn)*((E1L*E1)+(E2L*E2)));
	I = I/abs(Uapp);	
	cout << I << endl;
	It(step)=I;
//	Tt(step)=tt;
	
		{
		ofstream file("SUPG 3SP-dt-5e-13.txt");
		file << It << endl;
		}	
	
	Neold = Ne;
	Npold = Np;
	Nnold = Nn;
	tt += dt;
//	++step;
	}

```

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 9, 2024, 2:37pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/2 "2024-07-09T14:37:42Z")

</div>

You have many different orders of magnitude in your data. I suspect this is an issue. I would try to use quadruple precision, but I don’t know if it is possible in FreeFem…

---

<div class="post-metadata">

**Author:** ![Walid](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/walid/32/2477_2.png) [@Walid](https://community.freefem.org/u/Walid)\
**Post date:** [July 9, 2024, 3:30pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/3 "2024-07-09T15:30:17Z")

</div>

Thank you for you attention and reply.

I do not know what is " quadruple precision " is, but the model runs stable to a certain point then diverge. this is a common problem in different FEM software like COMSOL for this type of simulation and usually there are different approaches to stabilize the problem like SUPG or GLS but I really do not know why SUPG applying here did not add any difference.  
and please what do you mean by “You have many different orders of magnitude in your data.”

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 9, 2024, 4:39pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/4 "2024-07-09T16:39:17Z")

</div>

The different orders of magnitude means that you add and subtract quantities so that the difference becomes much smaller than each term separately and you need to evaluate the difference accurately.  
This happens in

```auto
	Se = S1-S2-S3;
	Sp = S1-S3-S4;
	Sn = S2-S4;

```

because one should have Se-Sp+Sn=0 but computationally this is not the case because of roundoff errors.  
Quadruple precision means that real numbers are computed with 32 digits instead of 16 for standard double precision.

Nevertheless I think that this is not the main problem here. I think that the main problem is to keep the quantities Ne, Np, Nn positive. The way you do, to cut them as  
`if(Ne[][i]<0.)Ne[][i]=0.;`  
is not good because it kills the consistency of your model and the conservativity properties.  
You should avoid that and rather look for a model which intrinsically keeps these quantities positive.  
I have tried to set boundary conditions Ne,Np,Nn vanishing on the whole boundary, and it seems to work robustly. Thus you have to check and select appropriate physically relevant boundary conditions. Here is the code with vanishing Dirichlet BC  
[convdiffcoupled.edp](https://community.freefem.org/uploads/short-url/bEdJwIp7rHZr56HYo1904spyX4d.edp) (8.0 KB)

---

<div class="post-metadata">

**Author:** ![Walid](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/walid/32/2477_2.png) [@Walid](https://community.freefem.org/u/Walid)\
**Post date:** [July 9, 2024, 5:23pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/5 "2024-07-09T17:23:35Z")

</div>

I agree with you regarding the line  
`if(Ne[][i]<0.)Ne[][i]=0.;`

but I really can not get or understand this line ?  
`real intNbil=int2d(Th)(Neold-Npold+Nnold);` why this signs ?

and there is a missing boundary condition could you help me to form in in FF++ it is applicable on boundaries number b10, b11 named out1

![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/e/ec233a60af75aa3a5f85e7671612ba0be38d9cd6.png)

and for the positivity of the concentrations there is an approach named Logarithmic Variational Formulation based using u = exp(U) instead of u such  
exp(U) ∂U/∂t + ∇(exp(U) · v − D exp(U) ∇U) = s.  
The main issue is that the variational problem is not linear anymore. and I can not understand hot to apply or solve this problem.

---

<div class="post-metadata">

**Author:** ![Walid](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/walid/32/2477_2.png) [@Walid](https://community.freefem.org/u/Walid)\
**Post date:** [July 9, 2024, 7:37pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/6 "2024-07-09T19:37:15Z")

</div>

Another question please I reviewed the code you had attached

```auto
problem NeL(Ne, v1) =	int2d(Th)(Ne*v1/dt	
									+(De*(grad(Ne)'*grad(v1))
									+Wex*Ne*dx(v1) + Wey*Ne*dy(v1)))
problem NpL(Np, v2) =	int2d(Th)(Np*v2/dt
									+(Dp*(grad(Np)'*grad(v2))
									-Wpx*Np*dx(v2) - Wpy*Np*dy(v2)))
problem NnL(Nn, v3) =	int2d(Th)(Nn*v3/dt	
									+(Dn*(grad(Nn)'*grad(v3))
									+Wnx*Nn*dx(v3) + Wny*Nn*dy(v3)))

```

you had changed the sign of the convection term in Ne and Nn to positive and Np to negative could you please explain the reason while the original equation is as below

![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/0/00e6934f53ee17c8777d55d9dfc808088d1d01ec.png)

I really appreciate your kind response

---

<div class="post-metadata">

**Author:** ![Walid](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/walid/32/2477_2.png) [@Walid](https://community.freefem.org/u/Walid)\
**Post date:** [July 10, 2024, 2:32am UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/7 "2024-07-10T02:32:28Z")

</div>

@fb77  
In the code you had attached there is differences in the weak form definition I am really confused about this line is correct

```auto
Wex*Ne*dx(v1) + Wey*Ne*dy(v1)))

```

or should it be

```auto
(Wex * dx(Ne) + Wey * dy(Ne))*v1))

```

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 10, 2024, 10:33am UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/8 "2024-07-10T10:33:50Z")

</div>

This was my question to you also. You wrote the equations to solve  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/0/00e6934f53ee17c8777d55d9dfc808088d1d01ec.png)  
The code lines I took correspond to that, according to (neglecting boundary terms)  
-\int v\_1\nabla\cdot (N\_e W\_e)=\int N\_e W\_e\cdot \nabla v\_1  
Hence we can use either  
`+Wex*Ne*dx(v1) + Wey*Ne*dy(v1)`  
or  
`-(Wex*dx(Ne)+Ne*dx(Wex) + Wey*dy(Ne)+Ne*dy(Wey))*v1`  
In any case we must choose appropriately the boundary conditions.

Do you have a document you refer to, describing all the boundary conditions, or do you have to set them up by yourself?

---

<div class="post-metadata">

**Author:** ![Walid](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/walid/32/2477_2.png) [@Walid](https://community.freefem.org/u/Walid)\
**Post date:** [July 10, 2024, 2:45pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/9 "2024-07-10T14:45:52Z")

</div>

[dordizadeh2015.pdf](https://community.freefem.org/uploads/short-url/ezClevTfc8VkdUWObTqkR5sVFoW.pdf) (2.8 MB)

here is the document

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 10, 2024, 2:58pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/10 "2024-07-10T14:58:47Z")

</div>

Ok I can set these boundary conditions of Table 1, and give it back to you.  
I take  
boundary#1=needle  
boundary#2=GND  
boundary#3=out1  
what about “ax”?

It seems to me that the potentiel V is represented as the unknown `u`, thus the equation  
-\Delta V=\frac{e(n\_p-n\_e-n\_n)}{\varepsilon\_0}\equiv f  
gives for any test function v  
\int \nabla V\cdot\nabla v-\int f v=\mbox{bdry terms}  
which is your problem Poisson, but your definition of f has the wrong sign, it should be  
`f = rohv/eps0`  
Do you agree?

---

<div class="post-metadata">

**Author:** ![Walid](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/walid/32/2477_2.png) [@Walid](https://community.freefem.org/u/Walid)\
**Post date:** [July 10, 2024, 3:33pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/11 "2024-07-10T15:33:32Z")

</div>

boundary ax is the axis of symmetry has no effect on the physics we used it to make the domain smaller and has no conditions

for Poisson equation  
V is unknown and representer by u correct I should modify the sign you are correct  
it would be as below

```auto
problem Poisson(u,v) = int2d(Th)(grad(u)'*grad(v))
						-int2d(Th)(f*v)
						+on (1,u = 0)
						+on (2,u = Uapp)
						;
	rohv = q*(Npold-Neold-Nnold);      
	f = rohv/eps0; 	

```

I am waiting you  
I would like to thank you very much for your kind attention and reply

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 10, 2024, 4:11pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/12 "2024-07-10T16:11:20Z")

</div>

Here is the code  
[convdiffcoupled.edp](https://community.freefem.org/uploads/short-url/aRVqqvxGeYN5HB0pqPoWEFDdbdd.edp) (8.3 KB)  
I removed the term S4 in the definition of Sp, because it is not in eq (6) of the paper. You can put it back if you want.  
Still we have to take a small enough timestep. Indeed if we want to strictly keep positive densities, dt should be restricted as  
dt\*D/h^2\leq 1/2  
where h is the mesh size, and D the diffusion coefficient.  
Your mesh has a fine part with h=4e-7, and for electrons D=0.18, which gives dt\< 4.4e-13  
I used dt=1e-11, it seems to work but negative values are present, they could lead to instability. It seems that instability occurs as soon as \int N\<0 for N=N\_e or N\_p or N\_n  
For dt=1e-11 it indeed fails at step #:259 time:2.59e-09 unfortunately.

---

<div class="post-metadata">

**Author:** ![Walid](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/walid/32/2477_2.png) [@Walid](https://community.freefem.org/u/Walid)\
**Post date:** [July 10, 2024, 4:50pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/13 "2024-07-10T16:50:49Z")

</div>

@fb77  
first of all I really want to thank you very much for your reply and attention.

I modified dt to be 1e-13 and run the software and will wait and see what will happen

but please I want to understand the differences between

```auto
//fb77 problem syntax 
problem NeL(Ne, v1) =	int2d(Th)(Ne*v1/dt	
									+(De*(grad(Ne)'*grad(v1))
									+Wex*Ne*dx(v1) + Wey*Ne*dy(v1)))
						-int2d(Th)(Se*v1)
						-int2d(Th)(Neold*v1/dt)
                       	-int1d(Th,needle)(gamma*Npold*mup*Emag*v1)
                        +int1d(Th,out1)(1.e30*Ne*v1*(E1*N.x+E2*N.y > 0. ? 1.:0.))					
						;
problem NpL(Np, v2) =	int2d(Th)(Np*v2/dt
									+(Dp*(grad(Np)'*grad(v2))
									-Wpx*Np*dx(v2) - Wpy*Np*dy(v2)))		
						-int2d(Th)(Sp*v2)
						-int2d(Th)(Npold*v2/dt)
  						+on (GND, Np = 0)	
                                                +int1d(Th,out1)(1.e30*Np*v2*(E1*N.x+E2*N.y < 0. ? 1.:0.))
						;	
					
problem NnL(Nn, v3) =	int2d(Th)(Nn*v3/dt	
									+(Dn*(grad(Nn)'*grad(v3))
									+Wnx*Nn*dx(v3) + Wny*Nn*dy(v3)))
						-int2d(Th)(Sn*v3)
						-int2d(Th)(Nnold*v3/dt)
						+on (needle, Nn = 0)			
                                                +int1d(Th,out1)(1.e30*Nn*v3*(E1*N.x+E2*N.y > 0. ? 1.:0.))
						;

// walid syntax 	
problem NeL(Ne, v1) =	int2d(Th)(Ne*v1/dt	
									+(De*(grad(Ne)'*grad(v1))
									+(-Wex * dx(Ne) + -Wey * dy(Ne))*v1))
						-int2d(Th)(Se*v1)
						-int2d(Th)(Neold*v1/dt)
						+int1d(Th,needle)(gamma*Npold*mup*Emag*v1)	
					;

problem NpL(Np, v2) =	int2d(Th)(Np*v2/dt
									+(Dp*(grad(Np)'*grad(v2))
									+(Wpx * dx(Np) + Wpy * dy(Np))*v2))		
						-int2d(Th)(Sp*v2)
						-int2d(Th)(Npold*v2/dt)
						+on (GND, Np = 0)	
						+int1d(Th,out1)(1.e30*Np*v2*(E1*N.x+E2*N.y < 0. ? 1.:0.))
						;	
					
problem NnL(Nn, v3) =	int2d(Th)(Nn*v3/dt	
									+(Dn*(grad(Nn)'*grad(v3))
									+(-Wnx * dx(Nn) + -Wny * dy(Nn))*v3))
						-int2d(Th)(Sn*v3)
						-int2d(Th)(Nnold*v3/dt)
						+on (needle, Nn = 0)
						+int1d(Th,out1)(1.e30*Nn*v3*(E1*N.x+E2*N.y > 0. ? 1.:0.))						
					;

```

these are two different ways to write the weak form and boundary conditions the first one is your script and the second is mine are they the same or one of them is correct or has a stable performance over the other?  
the part was writen by me is inspired from this line in the documentation of FF++ Pure Convection : The Rotating Hill

```auto
problem Adual(cc, w) = int2d(Th)( (cc/dt+(v1*dx(cc)+v2*dy(cc)))*w )

```

**I really want to understand the differences.**

the second question please if I need to apply streamline stability  
would it be add this term for each equation

```auto
+int2d(Th)(tau1*We*(dx(Ne)+dy(Ne))*We*(dx(v1)+dy(v1)))

tau1= hr1/2/sqrt(Wex^2+Wey^2);

```

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 10, 2024, 5:06pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/14 "2024-07-10T17:06:12Z")

</div>

This is my pleasure to help you.  
In your formulation you miss the terms of the type` dx(Wex)*Ne*v1`  
Thus the right way to write your formulation is

```auto
problem NeL(Ne, v1) =	int2d(Th)(Ne*v1/dt	
									+De*(grad(Ne)'*grad(v1))
									//+Wex*Ne*dx(v1) + Wey*Ne*dy(v1)
                                                                        -(Wex*dx(Ne)+dx(Wex)*Ne + Wey*dy(Ne)+dy(Wey)*Ne)*v1
                                                                        )
						-int2d(Th)(Se*v1)
						-int2d(Th)(Neold*v1/dt)
						-int1d(Th,needle)(gamma*Npold*mup*Emag*v1)
                                                +int1d(Th,out1)(1.e30*Ne*v1*(E1*N.x+E2*N.y > 0. ? 1.:0.))						
						;

problem NpL(Np, v2) =	int2d(Th)(Np*v2/dt
									+Dp*(grad(Np)'*grad(v2))
									//-Wpx*Np*dx(v2) - Wpy*Np*dy(v2)
                                                                        +(Wpx*dx(Np)+dx(Wpx)*Np + Wpy*dy(Np)+dy(Wpy)*Np)*v2
                                                                        )		
						-int2d(Th)(Sp*v2)
						-int2d(Th)(Npold*v2/dt)
						+on (GND, Np = 0)	
                                                +int1d(Th,out1)(1.e30*Np*v2*(E1*N.x+E2*N.y < 0. ? 1.:0.))
						;	
					
problem NnL(Nn, v3) =	int2d(Th)(Nn*v3/dt	
									+Dn*(grad(Nn)'*grad(v3))
									//+Wnx*Nn*dx(v3) + Wny*Nn*dy(v3)
                                                                        -(Wnx*dx(Nn)+dx(Wnx)*Nn + Wny*dy(Nn)+dy(Wny)*Nn)*v3
                                                                        )
						-int2d(Th)(Sn*v3)
						-int2d(Th)(Nnold*v3/dt)
						+on (needle, Nn = 0)			
                                                +int1d(Th,out1)(1.e30*Nn*v3*(E1*N.x+E2*N.y > 0. ? 1.:0.))
						;	

```

The boundary conditions are slightly different between this version and my previous version. This version corresponds strictly to the BC of the paper. My previous version corresponds to replacing all the Neumann conditions \nabla n\cdot\nu=0 (\nu the normal) by (Wn+D\nabla n)\cdot \nu=0 (this is for the n\_e equation, and similarly for the n\_p, n\_n equations). I think that (Wn+D\nabla n)\cdot \nu=0 is more natural, but I observe that with the condition \nabla n\cdot\nu=0 the computation is more stable (which is good!). Thus probably it is better to use this last version, full code:  
[convdiffcoupled.edp](https://community.freefem.org/uploads/short-url/81mXXaYYM6dqSDGozw06E2Lgvth.edp) (8.9 KB)  
that seems to work for dt=1e-10

About SUPG or streamline stability I don’t know.

The two formulations are related by  
-\int\_\Omega v\_1\nabla\cdot (N\_e \vec W\_e)=\int\_\Omega N\_e \vec W\_e\cdot \nabla v\_1-\int\_{\partial\Omega}v\_1N\_e\vec W\_e\cdot\vec{n}  
This formula is valid for any scalar functions v\_1 and N\_e, and any vector function \vec W\_e. In particular it is valid for P1 functions. This explains the difference between the two formulations: the difference is only by the boundary term.

---

<div class="post-metadata">

**Author:** ![Walid](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/walid/32/2477_2.png) [@Walid](https://community.freefem.org/u/Walid)\
**Post date:** [July 10, 2024, 5:22pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/15 "2024-07-10T17:22:55Z")

</div>

Thank you very much  
but I think there is something I need to clarify  
the We is the drift velocity and Wex, Wey are the x derivative and y derivative of drift velocity

```auto
We = mue * Emag;
Wex = dx(We) = mue * E1; where E1= - dx(u) 
Wey = dy(We) = mue * E2; where E2= - dy(u)

```

So the convective term in the weak form ∇(NeWe).v1  
will be

```auto
(We*dx(Ne)+dx(We)*Ne + We*dy(Ne)+dy(We)*Ne)*v1)

```

Or

```auto
(We*dx(Ne)+Wex*Ne + We*dy(Ne)+Wey*Ne)*v1)

```

I hope I am not annoying with too much replying  
thank you for your patience and reply

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 10, 2024, 5:34pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/16 "2024-07-10T17:34:59Z")

</div>

What is misleading is that mue is not constant. According to the paper eq (1), one has

```auto
Wex=mue*E1;
Wey=mue*E2;

```

But Wex is not equal to dx(We).  
The convective term is v\_1\nabla\cdot(N\_e[W\_{ex},W\_{ey}]) hence  
`(Wex*dx(Ne)+dx(Wex)*Ne + Wey*dy(Ne)+dy(Wey)*Ne)*v1`

One has W\_e=|[W\_{ex},W\_{ey}]|.

---

<div class="post-metadata">

**Author:** ![Walid](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/walid/32/2477_2.png) [@Walid](https://community.freefem.org/u/Walid)\
**Post date:** [July 10, 2024, 6:15pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/17 "2024-07-10T18:15:04Z")

</div>

Thank you I will modify the variational formulation  
but this is to wrap up all the post and ensure that I understand  
The strong form is  
∂Ne/∂t +∇(-We⋅Ne)-DeΔNe=Se -----------(1)  
∂Np/∂t +∇(Wp⋅Np)-DpΔNp=Sp -----------(2)  
∂Nn/∂t +∇(Wn⋅Nn)-DnΔNn=Sn -----------(2)

to obtain the weak form multiply by test function v1, v2, v3 and integrate over the domain  
∫(∂Ne/∂t +∇(-We⋅Ne)-DeΔNe)v1=Se.v1 -----------(1)

∫v1.∂Ne/∂t =

```auto
  int2d(Th)(Ne*v1/dt) - int2d(Th)(Neold*v1/dt)

```

∫∇(-We⋅Ne).v1

```auto
 = int2d(Th)(-(Wex*dx(Ne)+dx(Wex)*Ne + Wey*dy(Ne)+dy(Wey)*Ne)*v1)

```

∫DeΔNe)v1 =

```auto
int2d(Th)(De*(grad(Ne)'*grad(v1))

```

the boundary conditions

```auto
-int1d(Th,needle)(gamma*Npold*mup*Emag*v1)
 +int1d(Th,out1)(1.e30*Ne*v1*(E1*N.x+E2*N.y > 0. ? 1.:0.))	

```

and for Np the convection term sign to be changed from - to + correct ?  
so the three problems will be as below

```auto
problem NeL(Ne, v1) =	int2d(Th)(Ne*v1/dt	
									+De*(grad(Ne)'*grad(v1))
									//+Wex*Ne*dx(v1) + Wey*Ne*dy(v1)
                                                                        -(Wex*dx(Ne)+dx(Wex)*Ne + Wey*dy(Ne)+dy(Wey)*Ne)*v1)
						-int2d(Th)(Se*v1)
						-int2d(Th)(Neold*v1/dt)
						-int1d(Th,needle)(gamma*Npold*mup*Emag*v1)
                        +int1d(Th,out1)(1.e30*Ne*v1*(E1*N.x+E2*N.y > 0. ? 1.:0.))						
						;

problem NpL(Np, v2) =	int2d(Th)(Np*v2/dt
									+Dp*(grad(Np)'*grad(v2))
									+(Wpx*dx(Np)+dx(Wpx)*Np + Wpy*dy(Np)+dy(Wpy)*Np)*v2)		
						-int2d(Th)(Sp*v2)
						-int2d(Th)(Npold*v2/dt)
						+on (GND, Np = 0)	
                        +int1d(Th,out1)(1.e30*Np*v2*(E1*N.x+E2*N.y < 0. ? 1.:0.))
						;	
					
problem NnL(Nn, v3) =	int2d(Th)(Nn*v3/dt	
									+Dn*(grad(Nn)'*grad(v3))
									-(Wnx*dx(Nn)+dx(Wnx)*Nn + Wny*dy(Nn)+dy(Wny)*Nn)*v3 )
						-int2d(Th)(Sn*v3)
						-int2d(Th)(Nnold*v3/dt)
						+on (needle, Nn = 0)			
                                             +int1d(Th,out1)(1.e30*Nn*v3*(E1*N.x+E2*N.y > 0. ? 1.:0.))
						;	

```

---

<div class="post-metadata">

**Author:** ![Walid](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/walid/32/2477_2.png) [@Walid](https://community.freefem.org/u/Walid)\
**Post date:** [July 10, 2024, 6:48pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/18 "2024-07-10T18:48:47Z")

</div>

I am a little confused about this term

```auto
- int1d(Th,needle)(gamma*Npold*mup*Emag*v1)

```

why it is negative ?

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 10, 2024, 8:13pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/19 "2024-07-10T20:13:21Z")

</div>

I agree with your formulas.  
About the minus sign, it corresponds to this boundary term contributing to increase Ne. This is what I understood from the paper. If its effect should rather be to decrease Ne, the sign must be reversed.

---

<div class="post-metadata">

**Author:** ![Walid](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/walid/32/2477_2.png) [@Walid](https://community.freefem.org/u/Walid)\
**Post date:** [July 10, 2024, 8:25pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/20 "2024-07-10T20:25:11Z")

</div>

it is to increase yes.  
for Ne and Nn the sign would be like this

```auto
 = int2d(Th)(-(Wex*dx(Ne)+dx(Wex)*Ne + Wey*dy(Ne)+dy(Wey)*Ne)*v1)

```

and for Np the term to be

```auto
 = int2d(Th)((Wpx*dx(Np)+dx(Wpx)*Np + Wpy*dy(Np)+dy(Wpy)*Np)*v1)

```

and I will calculate only Wex, Wey, Wpx, Wpy, Wnx, Wny  
then calculate  
We = sqrt (Wex^2+Wey^2)  
Wp = sqrt (Wpx^2+Wpy^2)  
Wn = sqrt (Wnx^2+Wny^2)  
right ?

[Next page](https://community.freefem.org/t/stability-problem-related-to-supg/3382.md?page=2)