# Partial derivative of a vector field

**URL:** <https://fenicsproject.discourse.group/t/partial-derivative-of-a-vector-field/10172>\
**Category:** General\
**Created:** [January 17, 2023, 9:29pm UTC](https://fenicsproject.discourse.group/t/partial-derivative-of-a-vector-field/10172 "2023-01-17T21:29:57Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![charithjeewantha](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/charithjeewantha/32/3948_2.png) [@charithjeewantha](https://fenicsproject.discourse.group/u/charithjeewantha)\
**Post date:** [January 17, 2023, 9:29pm UTC](https://fenicsproject.discourse.group/t/partial-derivative-of-a-vector-field/10172/1 "2023-01-17T21:29:57Z")

</div>

I’m trying to compute the partial derivatives for the velocity vectors I obtained by solving the Navier Stokes equations.

For a working example (and to avoid mistakes from my end), I will use the **exact** code mentioned in the [Test problem 1: Channel flow (Poiseuille flow)](https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code1.html) (Sorry it is a bit long but it’s just the exact code in the tutorial) :

```auto
import numpy as np
import gmsh
import sys
from dolfinx import plot
# Initialize gmsh:
from mpi4py import MPI
#from dolfinx import mesh
import pyvista
from dolfinx.io.gmshio import model_to_mesh
import matplotlib.pyplot as plt
from dolfinx import geometry
# from dolfinx import fem,mesh
import ufl
from petsc4py import PETSc
from prettytable import PrettyTable
import random
from dolfinx.mesh import create_unit_square

from ufl import (FacetNormal, FiniteElement, FacetArea,Identity,TestFunction, TrialFunction, VectorElement,
                 div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym)
from dolfinx.fem import Expression,Constant,Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,locate_dofs_topological
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
t = 0
T = 10
num_steps = 500
dt = T/num_steps

v_cg2 = VectorElement("CG", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

def walls(x):
    return np.logical_or(np.isclose(x[1],0), np.isclose(x[1],1))
wall_dofs = locate_dofs_geometrical(V, walls)
u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bc_noslip = dirichletbc(u_noslip, wall_dofs, V)

def inflow(x):
    return np.isclose(x[0], 0)
inflow_dofs = locate_dofs_geometrical(Q, inflow)
bc_inflow = dirichletbc(PETSc.ScalarType(8), inflow_dofs, Q)

def outflow(x):
    return np.isclose(x[0], 1)
outflow_dofs = locate_dofs_geometrical(Q, outflow)
bc_outflow = dirichletbc(PETSc.ScalarType(0), outflow_dofs, Q)
bcu = [bc_noslip]
bcp = [bc_inflow, bc_outflow]

u_n = Function(V)
u_n.name = "u_n"
U = 0.5 * (u_n + u)
n = FacetNormal(mesh)
f = Constant(mesh, PETSc.ScalarType((0,0)))
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(1))
rho = Constant(mesh, PETSc.ScalarType(1))

# Define strain-rate tensor
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(u.geometric_dimension())

# Define the variational problem for the first step
p_n = Function(Q)
p_n.name = "p_n"
F1 = rho*dot((u - u_n) / k, v)*dx
F1 += rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx
F1 += inner(sigma(U, p_n), epsilon(v))*dx
F1 += dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds
F1 -= dot(f, v)*dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))

A1 = assemble_matrix(a1, bcs=bcu)
A1.assemble()
b1 = create_vector(L1)

# Define variational problem for step 2
u_ = Function(V)
a2 = form(dot(nabla_grad(p), nabla_grad(q))*dx)
L2 = form(dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)

# Define variational problem for step 3
p_ = Function(Q)
a3 = form(dot(u, v)*dx)
L3 = form(dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)

# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")

# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.BCGS)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")

# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)

xdmf = XDMFFile(mesh.comm, "poiseuille.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(u_n, t)
xdmf.write_function(p_n, t)

for i in range(num_steps):
    # Update current time step
    t += dt

    # Step 1: Tentative veolcity step
    with b1.localForm() as loc_1:
        loc_1.set(0)
    assemble_vector(b1, L1)
    apply_lifting(b1, [a1], [bcu])
    b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b1, bcu)
    solver1.solve(b1, u_.vector)
    u_.x.scatter_forward()
    
    # Step 2: Pressure corrrection step
    with b2.localForm() as loc_2:
        loc_2.set(0)
    assemble_vector(b2, L2)
    apply_lifting(b2, [a2], [bcp])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b2, bcp)
    solver2.solve(b2, p_.vector)
    p_.x.scatter_forward()

    # Step 3: Velocity correction step
    with b3.localForm() as loc_3:
        loc_3.set(0)
    assemble_vector(b3, L3)
    b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    solver3.solve(b3, u_.vector)
    u_.x.scatter_forward()
    # Update variable with solution form this time step
    u_n.x.array[:] = u_.x.array[:]
    p_n.x.array[:] = p_.x.array[:]

    # Write solutions to file
    xdmf.write_function(u_n, t)
    xdmf.write_function(p_n, t)

# Close xmdf file
xdmf.close()

```

Now I have the velocity values saved in the variable `u_n`  
And my goal is to compute the partial derivative of this `u_n` with respect to `x` at a point such as `(0.5,0,0)`

By looking at [Differential operators documentation](https://fenics.readthedocs.io/projects/ufl/en/latest/manual/form_language.html#automatic-functional-differentiation:~:text=and%203x3%20matrices.-,Differential%20Operators,-Three%20different%20kinds),

I can understand if, \overline{u}=(u,v,w), then:

\nabla \overline{u}=\begin{pmatrix} u\_x &u\_y &u\_z \\ v\_x&v\_y &v\_z \\ w\_x&w\_y &w\_z \end{pmatrix}

Thus,  
u\_x = \begin{pmatrix} 1 &0 &0 \end{pmatrix}\begin{pmatrix} u\_x &u\_y &u\_z \\ v\_x&v\_y &v\_z \\ w\_x&w\_y &w\_z \end{pmatrix}\begin{pmatrix} 1\\ 0\\ 0 \end{pmatrix}

And it is this u\_x that I need at a given point.  
But I cannot figure out the exact FEniCSx syntax that I should use for this calculation.

I did go over a couple of forum posts such as: [Diff() is it a derivative or a partial derivative? - General - FEniCS Project](https://fenicsproject.discourse.group/t/diff-is-it-a-derivative-or-a-partial-derivative/9526) but in those cases, the variable is explicitly defined with respect to the other. But in my case it is not.

Your help is much appreciated.

---

<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:** [January 18, 2023, 7:21am UTC](https://fenicsproject.discourse.group/t/partial-derivative-of-a-vector-field/10172/2 "2023-01-18T07:21:19Z")

</div>

There are two ways of doing this:

1. Interpolate `ufl.grad(u_n)` into a suitable function space (`fem.TensorFunctionSpace(mesh, ("DG", 1))` for all components, or `fem.FunctionSpace(mesh, ("DG", 1)` for only one component), and call eval, as described in [Implementation — FEniCSx tutorial](https://jsdokken.com/dolfinx-tutorial/chapter1/membrane_code.html#making-curve-plots-throughout-the-domain)  
This would mean a setup like

```python
Z = dolfinx.fem.FunctionSpace(mesh, ("DG", 1))
expr = dolfinx.fem.Expression(ufl.grad(u_n)[0,0], Z.element.interpolation_points()))
z = dolfinx.fem.Function(Z)
z.interpolate(expr)
# Then call z.eval after determining what cell the point/points are in 

```

1. If you only want to do it for a single cell, I would suggest using the `expr` directly, but changing `Z.element.interpolation_points()` with the pulled back point on the reference element (i.e. use: [dolfinx/fem.cpp at 48ece0087588f5f274a36aec223abb129fa59e8e · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/48ece0087588f5f274a36aec223abb129fa59e8e/python/dolfinx/wrappers/fem.cpp#L1094))

---

<div class="post-metadata">

**Author:** ![charithjeewantha](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/charithjeewantha/32/3948_2.png) [@charithjeewantha](https://fenicsproject.discourse.group/u/charithjeewantha)\
**Post date:** [January 20, 2023, 7:11am UTC](https://fenicsproject.discourse.group/t/partial-derivative-of-a-vector-field/10172/3 "2023-01-20T07:11:45Z")

</div>

Thank you so much @dokken This was a huge help!