# Weak form definition

**URL:** <https://community.freefem.org/t/weak-form-definition/3325>\
**Category:** General Discussion\
**Created:** [June 19, 2024, 7:37pm UTC](https://community.freefem.org/t/weak-form-definition/3325 "2024-06-19T19:37:45Z")\
**Posts on this page:** 12\
**Page:** 1

<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:** [June 19, 2024, 7:37pm UTC](https://community.freefem.org/t/weak-form-definition/3325/1 "2024-06-19T19:37:45Z")

</div>

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)(Ne_v1/dt+(De_(grad(Ne)'_grad(v1)) +(-Wex \* dx(Ne) + -Wey \* dy(Ne))v1))  
-int2d(Th)(Neoldv1/dt)  
+int1d(Th,needle) gamma_Npold_mup_Emag_v1)  
-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)(Ne_v1/dt   
+(De_(grad(Ne)'_grad(v1))  
+(-Wex \* dx(Ne) + -Wey \* dy(Ne))v1))  
-int2d(Th)(Neoldv1/dt)  
+int1d(Th,needle)(gamma_Npold_mup_Emag_v1)   
;  
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

---

<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:** [June 20, 2024, 11:02am UTC](https://community.freefem.org/t/weak-form-definition/3325/2 "2024-06-20T11:02:03Z")

</div>

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.

---

<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:** [June 20, 2024, 2:59pm UTC](https://community.freefem.org/t/weak-form-definition/3325/3 "2024-06-20T14:59:29Z")

</div>

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)(Ne_v1/dt  
+(De.(grad(Ne)'grad(v1))  
+(-Wex \* dx(Ne) + -Wey \* dy(Ne))v1))  
+int1d(Th,needle)(gammaNpoldmup_ 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  
+(De.(grad(Ne)'grad(v1))  
+(-Wex \* dx(Ne) + -Wey \* dy(Ne))v1))  
;  
varf NeR(Ne, v1) = int2d(Th)(Se \* v1)  
+int2d(Th)(Neoldv1/dt)  
-int1d(Th,needle)(gammaNpold_mup\* Emag\*v1);

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

---

<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:** [June 20, 2024, 3:28pm UTC](https://community.freefem.org/t/weak-form-definition/3325/4 "2024-06-20T15:28:37Z")

</div>

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.

---

<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:** [June 20, 2024, 3:53pm UTC](https://community.freefem.org/t/weak-form-definition/3325/5 "2024-06-20T15:53:44Z")

</div>

I am sorry I modified it

the varf statement is as below

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

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

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

```

is NBC.

---

<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:** [June 20, 2024, 5:22pm UTC](https://community.freefem.org/t/weak-form-definition/3325/6 "2024-06-20T17:22:43Z")

</div>

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

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

```

```auto
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**

---

<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:** [June 20, 2024, 5:40pm UTC](https://community.freefem.org/t/weak-form-definition/3325/7 "2024-06-20T17:40:00Z")

</div>

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)

---

<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:** [June 21, 2024, 9:12am UTC](https://community.freefem.org/t/weak-form-definition/3325/8 "2024-06-21T09:12:57Z")

</div>

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

```auto
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).

---

<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:** [June 21, 2024, 11:42pm UTC](https://community.freefem.org/t/weak-form-definition/3325/9 "2024-06-21T23:42:38Z")

</div>

> [@fb77](#):
>
> (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 ?

---

<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:** [June 23, 2024, 12:29pm UTC](https://community.freefem.org/t/weak-form-definition/3325/10 "2024-06-23T12:29:31Z")

</div>

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

```auto
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

```auto
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.

---

<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:** [June 23, 2024, 3:17pm UTC](https://community.freefem.org/t/weak-form-definition/3325/11 "2024-06-23T15:17:42Z")

</div>

Thank you for your attention and reply.

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ز

---

<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:** [June 23, 2024, 5:28pm UTC](https://community.freefem.org/t/weak-form-definition/3325/12 "2024-06-23T17:28:44Z")

</div>

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…