# Stability problem related to SUPG

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.

load "iovtk";
verbosity = 1;

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;

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

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


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…

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.”

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

	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 (8.0 KB)

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

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.

problem NeL(Ne, v1) =	int2d(Th)(Ne*v1/dt
+Wex*Ne*dx(v1) + Wey*Ne*dy(v1)))
problem NpL(Np, v2) =	int2d(Th)(Np*v2/dt
-Wpx*Np*dx(v2) - Wpy*Np*dy(v2)))
problem NnL(Nn, v3) =	int2d(Th)(Nn*v3/dt
+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

I really appreciate your kind response

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

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


or should it be

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


This was my question to you also. You wrote the equations to solve

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?

here is the document

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

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?

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

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

Here is the code
convdiffcoupled.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.

@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

//fb77 problem syntax
problem NeL(Ne, v1) =	int2d(Th)(Ne*v1/dt
+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
-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
+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
+(-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
+(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
+(-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

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

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

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


In your formulation you miss the terms of the type dx(Wex)*Ne*v1
Thus the right way to write your formulation is

problem NeL(Ne, v1) =	int2d(Th)(Ne*v1/dt
//+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
//-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
//+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 (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.

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

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

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


Or

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


I hope I am not annoying with too much replying

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

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}]|.

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 =

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


∫∇(-We⋅Ne).v1

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


∫DeΔNe)v1 =

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


the boundary conditions

-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

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


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


why it is negative ?

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.

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

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


and for Np the term to be

 =   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 ?