How to create the advection matrix

variational formulation
1

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?

2

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.

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

will make the second addend a scalar too.

3

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:

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

or

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

where v is a function whose value is 1

4

Isn’t this question very similar to: Is it possible to create a matrix with linear form - #4 by dokken

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

5

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.

6

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).
1 Like
7

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:

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

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);
8

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.