# How to create the advection matrix

**URL:** <https://fenicsproject.discourse.group/t/how-to-create-the-advection-matrix/16262>\
**Category:** variational formulation\
**Tags:** dolfinx\
**Created:** [November 10, 2024, 11:16pm UTC](https://fenicsproject.discourse.group/t/how-to-create-the-advection-matrix/16262 "2024-11-10T23:16:12Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![WjjSteve](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wjjsteve/32/4389_2.png) [@WjjSteve](https://fenicsproject.discourse.group/u/WjjSteve)\
**Post date:** [November 10, 2024, 11:16pm UTC](https://fenicsproject.discourse.group/t/how-to-create-the-advection-matrix/16262/1 "2024-11-10T23:16:12Z")

</div>

Hello,  
I’m trying to solve a advection equation  
_\partial\_t f+\nabla f=0_  
after apply FEM, we obtain the weak form  
_(\partial\_t f,\phi\_i)+(\nabla f,\phi\_i)=0_  
Usually, the second term is assembled as a vector in the linear form.

But I want to use this **advection matrix** , whose element is computed as  
_a\_{ij}=\int\_{\Omega}\nabla\phi\_j\phi\_i\ {\rm d}x_

so I create the bilinear form (I work in 2D)

> f = TrialFunction(V)  
> w = TestFunction(V)  
> C = f.dx(0) \* w \* dx + f.dx(1) \* w \* dx

Then I assemble the matrix as I do for the Mass matrix, but the obtained matrix **is not correct** ,  
I used a coarse mesh and compute the matrix using my pen, and I think the elements of the matrix are not the correct number  
Does anyone have any idea about this question?

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [November 11, 2024, 5:18am UTC](https://fenicsproject.discourse.group/t/how-to-create-the-advection-matrix/16262/2 "2024-11-11T05:18:36Z")

</div>

> [@WjjSteve](#):
>
> Hello,  
> I’m trying to solve a advection equation  
> _\partial\_t f+\nabla f=0∂tf+∇f=0\partial\_t f+\nabla f=0_

This equation makes no sense. If f is a scalar, \partial\_t f is a scalar, while \nabla f is a vector, and you definitely cannot sum them up.

> [@WjjSteve](#):
>
> after apply FEM, we obtain the weak form  
> _(\partial\_t f,\phi\_i)+(\nabla f,\phi\_i)=0(∂tf,ϕi)+(∇f,ϕi)=0(\partial\_t f,\phi\_i)+(\nabla f,\phi\_i)=0_  
> Usually, the second term is assembled as a vector in the linear form.

Again, the first addend will be the integral of a scalar, while the second addend the integral of a vector, but

> [@WjjSteve](#):
>
> ```auto
> C = f.dx(0) * w * dx + f.dx(1) * w * dx
> 
> ```

will make the second addend a scalar too.

---

<div class="post-metadata">

**Author:** ![WjjSteve](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wjjsteve/32/4389_2.png) [@WjjSteve](https://fenicsproject.discourse.group/u/WjjSteve)\
**Post date:** [November 11, 2024, 8:34am UTC](https://fenicsproject.discourse.group/t/how-to-create-the-advection-matrix/16262/3 "2024-11-11T08:34:15Z")

</div>

Sorry I forgot the advection field  
_\partial\_t f+\pmb{\beta}\cdot\nabla f= 0_  
if we consider the simpliest case where _\pmb{\beta}=(1\quad1)^\mathsf{T}_

can I write my form like this:

```auto
C = f.dx(0) * w * dx + f.dx(1) * w * dx

```

or

```auto
C =v * f.dx(0) * w * dx + v * f.dx(1) * w * dx

```

where v is a function whose value is 1

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 11, 2024, 8:43am UTC](https://fenicsproject.discourse.group/t/how-to-create-the-advection-matrix/16262/4 "2024-11-11T08:43:19Z")

</div>

Isn’t this question very similar to: [Is it possible to create a matrix with linear form - #4 by dokken](https://fenicsproject.discourse.group/t/is-it-possible-to-create-a-matrix-with-linear-form/16243/4)

I.e. use trial functions for `f` instead of a `Function`. Then you can post-multiply by the `f` coefficients.

---

<div class="post-metadata">

**Author:** ![WjjSteve](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wjjsteve/32/4389_2.png) [@WjjSteve](https://fenicsproject.discourse.group/u/WjjSteve)\
**Post date:** [November 11, 2024, 8:45am UTC](https://fenicsproject.discourse.group/t/how-to-create-the-advection-matrix/16262/5 "2024-11-11T08:45:01Z")

</div>

Yes this is the same question actually,  
now I have obtained this matrix, without any coefficient, but I think the numbers are not correct, so I wonder if that I made any mistake.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 11, 2024, 8:50am UTC](https://fenicsproject.discourse.group/t/how-to-create-the-advection-matrix/16262/6 "2024-11-11T08:50:10Z")

</div>

A simple way to verify the numbers are:

- Create a standard “vector” assembly, where you choose a given `f`. Assemble that vector
- Compare that to the matrix approach (post multiplication).

---

<div class="post-metadata">

**Author:** ![WjjSteve](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wjjsteve/32/4389_2.png) [@WjjSteve](https://fenicsproject.discourse.group/u/WjjSteve)\
**Post date:** [November 11, 2024, 9:33am UTC](https://fenicsproject.discourse.group/t/how-to-create-the-advection-matrix/16262/7 "2024-11-11T09:33:52Z")

</div>

Ok, now I think this is a problem of my code about assemble the matrix,

I checked it in a short python code, and the obtained matrix is just as expected, but this doesn’t work in CPP, maybe there is some bug here:

```auto
C = f.dx(0) * w * dx + f.dx(1) * w * dx

```

```auto
auto C = std::make_shared<Matrix>(fem::petsc::create_matrix(*c_blinear), false);
MatZeroEntries(C->mat());
dolfinx::fem::assemble_matrix(Matrix::set_block_fn(C->mat(), ADD_VALUES), *c_blinear, {});
MatAssemblyBegin(C->mat(), MAT_FLUSH_ASSEMBLY);
MatAssemblyEnd(C->mat(), MAT_FLUSH_ASSEMBLY);
fem::set_diagonal<T>(Matrix::set_fn(C->mat(), INSERT_VALUES), *V, {});
MatAssemblyBegin(C->mat(), MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(C->mat(), MAT_FINAL_ASSEMBLY);  

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 11, 2024, 10:00am UTC](https://fenicsproject.discourse.group/t/how-to-create-the-advection-matrix/16262/8 "2024-11-11T10:00:00Z")

</div>

> [@WjjSteve](#):
>
> ```auto
> MatAssemblyBegin(C->mat(), MAT_FLUSH_ASSEMBLY);
> MatAssemblyEnd(C->mat(), MAT_FLUSH_ASSEMBLY);
> fem::set_diagonal<T>(Matrix::set_fn(C->mat(), INSERT_VALUES), *V, {});
> 
> ```

Since you don’t have bcs, this is not needed.

I would start with a single cell problem (Unit square with single quadrilateral cell). Then you can verify that a single cell contribution in your assembly is the same.

Note that no-one can really help you debug this as long as you are only sharing snippets of code.