# 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:** 3\
**Page:** 2

<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 22, 2024, 1:34pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/21 "2024-07-22T13:34:23Z")

</div>

Hello Sir,

Unfortunately, this weak form

```auto
problem NeL(Ne, v1) =	int2d(Th)(Ne*v1/dt	
									+(De*(grad(Ne)'*grad(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.))							
						;

```

did not give a correct result.

Do you thing calculating We then using the below weak form could be correct or this is not correct

```auto
We = mue * Emag;

problem NeL(Ne, v1) =	int2d(Th)(Ne*v1/dt	
									+(De*(grad(Ne)'*grad(v1))
									+(-dx(We)* dx(Ne) + -dy(We) * dy(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.))							
						;

```

Thanks 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 22, 2024, 3:35pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/22 "2024-07-22T15:35:22Z")

</div>

Hello Walid,  
according to the paper you sent, the equation for n\_e is  
\partial\_t n\_e+\nabla\cdot(-\mu\_e\vec E n\_e -D\_e\nabla n\_e)=S\_e.  
The advection term is  
-\nabla\cdot(\mu\_e\vec E n\_e)=-\bigl(\partial\_x(\mu\_e E\_x n\_e)+\partial\_y(\mu\_e E\_y n\_e)\bigr)  
=-\bigl(\partial\_x(\mu\_e E\_x)\times n\_e+\mu\_eE\_x\partial\_x n\_e +\partial\_y(\mu\_e E\_y)\times n\_e+\mu\_eE\_y\partial\_y n\_e\bigr).  
Multiplied by `v1` it well corresponds to the formula  
`-(Wex*dx(Ne)+dx(Wex)*Ne + Wey*dy(Ne)+dy(Wey)*Ne)*v1`  
with W\_x=\mu\_eE\_x, W\_y=\mu\_eE\_y.  
But \partial\_x(W\_x n\_e)\not=(\partial\_x W\_x)\times(\partial\_x n\_e), thus your proposed formulation is not correct.  
The reason for getting not good results must be somewhere else.  
Reading the code I don’t see any mistake remaining.  
How do you evaluate your numerical results as incorrect?  
The system you try to solve may blow up in finite time because of the quadratic terms. Therefore problems may come either from a numerical instability, or from the nature of the system.

---

<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:** [August 19, 2024, 6:01pm UTC](https://community.freefem.org/t/stability-problem-related-to-supg/3382/23 "2024-08-19T18:01:43Z")

</div>

This system is already solved by me using COMSOL Multiphysics so I compare the results with the one I got from COMSOL  
and also the system divergence however I am using SUPG (as I understand)  
and the documentation of COMSOL mentions that the use SUPG and the final code is as below if you can review it one more time I really appreciate your attention

```auto
load "iovtk";
verbosity = 1; 
				
real Uapp=-4000; // 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 = 1e-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;

VNe ADe, ADp, ADn;

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

macro ffSaveData4(u1, u2, u3, u4, filename){
  {
  ofstream file( filename+".txt");
  int datalen=u1[].n;
  if ((u2[].n!=datalen) | (u3[].n!=datalen) | (u4[].n!=datalen)){
    cout << "error: arguments must have same size" << endl;
    exit(1);
  }
  file.precision(1);
  for (int j=0; j<datalen; j++){
    file << u1[][j] << " " << u2[][j] << " " << u3[][j] << " " << u4[][j] << endl;
  }
  }
} //EOM
macro ffLoadData(u1, u2, u3, u4, Vh, filename){
  {
  
  ifstream file(filename + ".txt");

  for (int j=0; j<Vh.ndof; j++){
    file >> u1[][j] >> u2[][j] >> u3[][j] >> u4[][j];
  }
  }
} //EOM

macro ffSaveVh(Th, Vh, filename){
  {
  ofstream file(+ filename + ".txt");
  file.precision(10);
  for (int i=0; i<Th.nt; i++){
    for (int j=0; j<Vh.ndofK; j++){
      file << Vh(i,j) << "\n";
    }
  }
  }
} //EOM

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

//ffLoadData (Neold, Nnold, Npold, Emag, VE, "Concentrations45250");	

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)+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.))

						+int2d(Th)((tau1*(Wex*dx(v1)+Wey*dy(v1)))*(Ne/dt-De*(dxx(Ne)+dyy(Ne))+-(Wex*dx(Ne)+dx(Wex)*Ne + Wey*dy(Ne)+dy(Wey)*Ne)))
						-int2d(Th)((tau1*(Wex*dx(v1)+Wey*dy(v1)))*(Neold/dt)) //SUPG 				
						-int2d(Th)((tau1*(Wex*dx(v1)+Wey*dy(v1)))*(Se)) //SUPG 							
						//+int2d(Th)(tau1*We*(dx(Ne)+dy(Ne))*We*(dx(v1)+dy(v1)))
						;

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.))
						
						+int2d(Th)((tau2*(Wpx*dx(v2)+Wpy*dy(v2)))*(Np/dt-Dp*(dxx(Np)+dyy(Np))+ (Wpx*dx(Np)+dx(Wpx)*Np + Wpy*dy(Np)+dy(Wpy)*Np)))
						-int2d(Th)((tau2*(Wpx*dx(v2)+Wpy*dy(v2)))*(Npold/dt)) //SUPG 				
						-int2d(Th)((tau2*(Wpx*dx(v2)+Wpy*dy(v2)))*(Sp)) //SUPG 
						//+int2d(Th)(tau2*Wp*(dx(Np)+dy(Np))*Wp*(dx(v2)+dy(v2)))
						;	
					
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.))						
						
						+int2d(Th)((tau3*(Wnx*dx(v3)+Wny*dy(v3)))*(Nn/dt-Dn*(dxx(Nn)+dyy(Nn))+-(Wnx*dx(Nn)+dx(Wnx)*Nn + Wny*dy(Nn)+dy(Wny)*Nn)))
						-int2d(Th)((tau3*(Wnx*dx(v3)+Wny*dy(v3)))*(Nnold/dt)) //SUPG 				
						-int2d(Th)((tau3*(Wnx*dx(v3)+Wny*dy(v3)))*(Sn)) //SUPG 
						//+int2d(Th)(tau3*Wn*(dx(Nn)+dy(Nn))*Wn*(dx(v3)+dy(v3)))
						;	
						
//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;

        real intNe=int2d(Th)(Neold);
        real intNp=int2d(Th)(Npold);
        real intNn=int2d(Th)(Nnold);
        real intNbil=int2d(Th)(Neold-Npold+Nnold);
        cout << " intN = " << intNe << " " << intNp << " " << intNn << " " << intNbil << 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);

	Wex = E1*mue; 
	Wey = E2*mue;
	Wpx = E1*mup; 
	Wpy = E2*mup;
	Wnx = E1*mun;
	Wny = E2*mun;

	We = sqrt(Wex^2+Wey^2);
	Wp = sqrt(Wpx^2+Wpy^2);
	Wn = sqrt(Wnx^2+Wny^2);
	
	tau1= hr1/2/We;
	tau2= hr1/2/Wp;
	tau3= hr1/2/Wn;

	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;
	
//	ADe = De+We*hr1; 
//	ADp = Dp+Wp*hr1;
//	ADn = Dn+Wn*hr1;
	
	NeL;
//	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 ;

	NpL;
//	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 PeK1 = We*hr1/2./De ;
//	cout << "Pe max= " << PeK1[].max << endl ;
//	cout << "Pe min= " << PeK1[].min << endl ;

//	VNe PeK2 = Wp*hr1/2./Dp ;
//	cout << "Pep max= " << PeK2[].max << endl ;
//	cout << "Pep min= " << PeK2[].min << endl ;		

//	VNe PeK3 = Wn*hr1/2./Dn ;
//	cout << "Pen max= " << PeK3[].max << endl ;
//	cout << "Pen min= " << PeK3[].min << endl ;
	
		if (step % 50 == 0)
		{
		savevtk("results/TNS "+step+".vtu", Th, Emag, Ne, Np, Nn);
		ffSaveData4(Neold, Nnold, Npold, Emag, "Concentrations"+step);
		}
	
// 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*Neold + mup*Npold + mun*Nnold)*((E1L*E1)+(E2L*E2)));
	I = I/abs(Uapp);	
	cout << I << endl;
	It(step)=I;
//	Tt(step)=tt;
	
		{
		ofstream file("AD 3SP-dt-1e-13.txt");
		file << It << endl;
		}	
	
	Neold = Ne;
	Npold = Np;
	Nnold = Nn;
	tt += dt;
//	++step;
	}

```

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