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)(Neoldv1/dt)
+int1d(Th,needle) gammaNpoldmupEmagv1)
-int2d(Th)(Sev1)
; 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(grad(Ne)'grad(v1))
+(-Wex * dx(Ne) + -Wey * dy(Ne))v1))
-int2d(Th)(Neoldv1/dt)
+int1d(Th,needle)(gammaNpoldmupEmagv1)
;
varf NeR(Ne, v1) = int2d(Th)(Sev1);
and to solve it
matrix AE;
real [int] bE(VNe.ndof);
AE=NeL(VNe,VNe);
bE=NeR(0,VNe);
Ne [ ] = AE^-1*bE ;
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.
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.
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
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)
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).
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
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.
but I want to make it more clear; the main goal is solve 4 PDE (1 Poisson and 3 Hydrodynamic equation) The code you offered could be also modified to solve the 4 together right ?
Finally, I would like to express my gratitude and sincere thanksز
In principle you could solve the 4 together if they are linear in the vector of all unknowns (which is the case for my example).
If the 4 together is not possible because of some nonlinearities or because it is computationally too demanding, then to find a “good” scheme (stable and not too much expensive in cpu) is difficult in general, it depends on the structure of your system and needs some elaborate numerical analysis…