# Weak form definition

Dear All, i am asking for your kind support for the below problem

I am trying to write the weak form of the below PDE
∂n/∂t +∇⋅(W.n+D∇n)=S
where:
W drift velocity
n unknow concentration
D the diffusion coefficient
S the source term

when I use problem statement to define the weak form the results are correct but unstable.
the problem statement is as below
problem HydrodynamicNe(Ne, v1) = int2d(Th)(Nev1/dt+(De(grad(Ne)'grad(v1)) +(-Wex * dx(Ne) + -Wey * dy(Ne))v1))
-int2d(Th)(Neold
v1/dt)
+int1d(Th,needle) gamma
NpoldmupEmagv1)
-int2d(Th)(Se
v1)
;
and to solve it
HydrodynamicNe;

when I use varf statement I got wrong results
the varf statement is as below
varf NeL(Ne, v1) = int2d(Th)(Nev1/dt
+(De
+(-Wex * dx(Ne) + -Wey * dy(Ne))v1))
-int2d(Th)(Neold
v1/dt)
+int1d(Th,needle)(gamma
NpoldmupEmagv1)
;
varf NeR(Ne, v1) = int2d(Th)(Se
v1);

and to solve it

matrix AE;
real [int] bE(VNe.ndof);
AE=NeL(VNe,VNe);
bE=NeR(0,VNe);
Ne [ ] = AE^-1*bE ;

I need your support to figure out the problem

The two last lines of varf NeL do not contain the unknown Ne. Therefore you should put them in the varf NeR (with a minus sign), they are part of the right-hand side.

1 Like

Thank you very much for your aatention and reply. hence, the weak formulation presentation would be

the varf statement is as below
varf NeL(Ne, v1) = int2d(Th)(Nev1/dt
+(-Wex * dx(Ne) + -Wey * dy(Ne))v1))
+int1d(Th,needle)(gamma
Npold
mup
Emagv1)
;
varf NeR(Ne, v1) = int2d(Th)(Se * v1)
+int2d(Th)(Neold
v1/dt);

OR
varf NeL(Ne, v1) = int2d(Th)(Nev1/dt
+(-Wex * dx(Ne) + -Wey * dy(Ne))v1))
;
varf NeR(Ne, v1) = int2d(Th)(Se * v1)
+int2d(Th)(Neold
v1/dt)
-int1d(Th,needle)(gamma
Npold
mup* Emag*v1);

also I would like to confirm that this part
int1d(Th,needle)(gammaNpoldmup* Emag*v1);
is NBC.

Please put your code lines between three backqotes (a backquote is this char ), three at the beginning and three at the end, so that it is easier to read it.

I am sorry I modified it

the varf statement is as below

varf NeL(Ne, v1) = int2d(Th)(Ne*v1/dt
+(-Wex * dx(Ne) + -Wey * dy(Ne))* v1))
+int1d(Th,needle)(gamma* Npold* mup*  Emag* v1);

varf NeR(Ne, v1) = int2d(Th)(Se * v1)
+int2d(Th)(Neold* v1/dt);


OR

varf NeL(Ne, v1) = int2d(Th)(Ne* v1/dt
+(-Wex * dx(Ne) + -Wey * dy(Ne))* v1));

varf NeR(Ne, v1) = int2d(Th)(Se * v1)
+int2d(Th)(Neold* v1/dt)
-int1d(Th,needle)(gamma* Npold* mup* Emag*v1);


also I would like to confirm that this part

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


is NBC.

Your varf definitions correspond now to the problem HydrodynamicNe you have defined. However for being consistent with the conservative drift term and for applying Neumann BC the NeL part should be

varf NeL(Ne, v1) = int2d(Th)(Ne* v1/dt
-Wex * Ne*dx(v1) - Wey * Ne*dy(v1)   );

and the Ner part should be
varf NeR(Ne, v1) = int2d(Th)(Se * v1)
+int2d(Th)(Neold* v1/dt);


thus no boundary term (I don’t know what is your gamma* Npold* mup* Emag, it does not appear in your equation statement).
Note also that I consider your sign in the diffusion term to be wrong in your equation statement, it should be ∂n/∂t +∇⋅(W.n-D∇n)=S

Yes sir, you are right regarding the De sign I mistype it

gamma and mup are real numbers
Npold, Emag is a finite element variables are being calculated from other PDEs during the software running.

So the sign of the Neumann BC will be positive in the right hand side varf definition, correct?

I have another small question
loading PETSc will enhance the solution stability ? or I must apply one of the stabilization techniques like (artificial diffusion or SUPG)

The lines I sent were to set homogeneous Neumann BC, which is
(Wn-D\nabla n)\cdot N=0 on the whole boundary.
If you include the boundary integral as

varf NeR(Ne, v1) = int2d(Th)(Se * v1)
+int2d(Th)(Neold* v1/dt)
-int1d(Th,needle)(gamma* Npold* mup* Emag*v1);


it means that (Wn-D\nabla n)\cdot N=\gamma N_p \mu_p E on “needle” (and =0 on the rest of the boundary).

PETSc will not change the stability issue.
If you have an instability, since your system is a bit complex (several coupled unknowns), it is difficult to say a priori. Either you have an analysis that ensures you that a good choice is known and should work, or you can just try several recipes of stabilization or slightly change the way you discretize equations (for example solve both equations simultaneously in a single variational formulation instead of one after the other).

Is there is an example to show how this could be done ?

An example is as follows. Consider the coupled problem
\partial_t u-\Delta u-v=0,
\partial_t v+\partial_x u-\Delta v=0,
with homogeneous Dirichlet BC for both.
A resolution of one equation after the other is

solve probu(u,ut)=
int2d(u*ut/dt+dx(u)*dx(ut)+dy(u)*dy(ut))
-int2d(uold*ut/dt)
-int2d(vold*ut)
+on(1,u=0.)
;
solve probv(v,vt)=
int2d(v*vt/dt+dx(v)*dx(vt)+dy(v)*dy(vt))
-int2d(vold*vt/dt)
-int2d(u*dx(vt))
+on(1,v=0.)
;


A simultaneous resolution is

solve probuv(u,v,ut,vt)=
int2d(u*ut/dt+dx(u)*dx(ut)+dy(u)*dy(ut))
-int2d(uold*ut/dt)
-int2d(v*ut)
+int2d(v*vt/dt+dx(v)*dx(vt)+dy(v)*dy(vt))
-int2d(vold*vt/dt)
-int2d(u*dx(vt))
+on(1,u=0.,v=0.)
;


The difference is here only in the term -int2d(v*ut) which is now implicit, instead of explicit as -int2d(vold*ut)` previously.
The simultaneous resolution is in general more stable, but more CPU consuming than the resolution of one equation after the other.