Stability problem related to SUPG

Walid · July 22, 2024, 1:34pm

Hello Sir,

Unfortunately, this weak form

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

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.

fb77 · July 22, 2024, 3:35pm

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.

Walid · August 19, 2024, 6:01pm

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

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

Topic		Replies	Views
Problem related to SUPG General Discussion	7	270	July 3, 2024
Coupled high order PDE General Discussion	2	376	September 23, 2020
Third derivative General Discussion	18	450	June 20, 2024
Doubt in the variational formulation of LaplaceRT.edp General Discussion	5	245	May 24, 2023
Solving a BVP with an algebraic equation General Discussion	0	133	November 14, 2023

Stability problem related to SUPG

Related Topics