Constructing discrete gradient operator from a group finite element formulation?

Question

Suppose we have the scalar advection equation with compressible velocity $\vec{v}$:
$$ \frac{du}{dt} + \nabla \cdot( \vec{v} u ) = 0 $$
in some domain $\Omega$.

Using the group finite element formulation with $\phi$ representing the basis functions such that

$$ u = \sum_j u_j \phi_j $$

and
$$ (\vec{v} u)_h = \sum_j u_j ( \vec{v}_j \cdot \nabla \phi_j ) $$

we can write the semi-discrete form of the advection equation as follows:

$$\sum_j \left( \int_{\Omega} \phi_i \phi_j dx \right) \frac{du_j}{dt} = - \sum_j \vec{v}_j \cdot \left( \int_{\Omega}\phi_i \nabla \phi_j dx \right) u_j $$

My question is this: is there a way to construct the discrete gradient defined by
$$ \vec{c}_{ij} = \left( \int_{\Omega}\phi_i \nabla \phi_j dx \right) $$

as a set of PETScMatrices (for each spatial variable) to be used for computing the discrete transport operator in the PDE above? This is to avoid reassembly of the gradient mainly since the velocity will be changing in time in my problem.

Comments

Hi, what is the function space for velocity $\vec{v}$ ?

---

$u$ lives in a CG1 space, while the velocity $\vec{v}$ can either live in CG1 or CG2 depending on the initial setup.

---

I've attempted to setup this operation as follows assuming $\vec{v}$ lives in a CG1 space.

For $\vec{c}_{ij}$ I've assembled each term separately, i.e.

$$ c^x_{ij} = \int_{\Omega} \phi \frac{\partial \phi}{\partial x} \; dx $$

$$ c^y_{ij} = \int_{\Omega} \phi \frac{\partial \phi}{\partial y} \; dx $$

the corresponding UFL code looks like:

  Q = FunctionSpace(mesh, "CG",1)
  q = TrialFunction(Q)
  u = TestFunction(Q)
  gradX_form = q * Dx(u,0) * dx
  gradY_form = q * Dx(u,1) * dx

With this, we can write the dot product $ \vec{c}_{ij} \cdot \vec{v}_j$ as:
$$ \vec{c}_{ij} \cdot \vec{v}_j = c^x_{ij} * v_{x,j} + c^y_{ij} * v_{y,j} $$
which I currently have setup as a loop over the scipy.coo_matrix forms of $c^x$ and $c^y$ with the corresponding $j$ indices for $v$.

However, comparing this construction to the assembly of the advection operator (assuming uvel is in a CG1 vector function space) directly, i.e.

adv =  -1.0 * div( uvel * u ) * q * dx

There is severe disagreement. It looks like the rows of the matrices assembled are not the same indexing as a CG1 function (both vector and scalar)? Is this true?

I should also comment that I have:

parameters['reorder_dofs_serial'] = False 

EDIT: This solution works. I just had to remember that the entries of this formulation is different than the standard Galerkin approach. Things seem to be doing what I expect now...

---

Method seems to be suffering from some consistency issues.

I've traced itback to assuming that:
$$ c^x_{ij} = \int_{\Omega} \phi \frac{\partial \phi}{\partial x} \; dx $$

can be written in UFL as:

gradX_form = q * Dx(u,0) * dx

Not really sure what to try next here...