How to represent directional derivative of 1D elements in 3D space?

variational formulation
1

Hello FEniCS Team, I try to solve a convection-diffusion-reaction NonlinearProblem which involves a network consisting of many 1D line elements in 3D space. Because there are so many line elements in 3D space, I need to compute the directional derivative rather than ‘u.dx(0)’. But there are some problems related to ‘ufl.sqrt’. Any comments/suggestions would be greatly appreciated. I use “DOLFINx version: 0.6.0.0 based on GIT commit: 72c70a59b8bed7cc5caa4c68a84182c9c4be87cf of GitHub - FEniCS/dolfinx: Next generation FEniCS problem solving environment”. Thanks a lot!

Here is a Minimal Reproducible Example:
####################################

import numpy as np
import dolfinx
from ufl import (TrialFunction, TestFunction, dx, ds, dS, dot, inner, jump, div, grad,
                 lhs, rhs, as_vector, FacetNormal, Measure, derivative)
from ufl.operators import sqrt
from dolfinx.fem import *
from dolfinx.mesh import create_interval, exterior_facet_indices, locate_entities, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import NonlinearProblem, assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.nls.petsc import NewtonSolver

# Compute directional derivative
def dudx(u):
    directional_derivative_u = sqrt(dot(grad(u), grad(u)))
    return directional_derivative_u

# Create mesh     
L = 1.0; dt = 0.1
domain = create_interval(MPI.COMM_WORLD, 10, (0.0, L))  # In fact, I have many 1D lines in 3D space!
V = FunctionSpace(domain, ('CG', 1))
du = TrialFunction(V)
v = TestFunction(V)
u = Function(V)
u0 = Function(V)

# BCs
dofs_L = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0.0))
dofs_R = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], L))
uL, uR = Function(V), Function(V)
uL.vector.array = 0.0
uR.vector.array = 0.0
bcL, bcR = dirichletbc(uL, dofs_L), dirichletbc(uR, dofs_R)
bcs = [bcL, bcR]


F = v*(u-u0)/dt*dx + 0.01*v*dudx(u)*dx + 0.2*dudx(u)*dudx(v)*dx + 0.001*u/(1.0+u)*v*dx
a, L = lhs(F), rhs(F)
Jac = derivative(F, u, du)

# initial conditions
formFb = form(rhs(F))
init_cond = create_vector(formFb)
with init_cond.localForm() as loc:
    loc.set(ScalarType(50.0))
u0.interpolate(init_cond)

problem = NonlinearProblem(F, u, bcs=bcs, J=Jac)
solver = NewtonSolver(domain.comm, problem)

solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

num_steps = 10
for n in range(num_steps):
    num_its, converged = solver.solve(u)
    print(f"--> Time step {n}, Number of iterations {num_its}")

####################################

Here is part of the error message:
…
Found Argument in , this is an invalid expression.
ERROR:UFL:Found Argument in , this is an invalid expression.

2

Will you need to integrate have a 1D integral in the 3D domain to write your weak formulation, or only compute a directional derivative but still integrate it in the 3D domain?

3

Thank you, Francesco! I was trying to compute the directional derivative which later as a SCALAR enters the weak formulation. If it is challenging to directly use directional_derivative_u = sqrt(dot(grad(u), grad(u))), I also wonder if it is possible to extract the coordinates (or its two vertices/nodes) of each line cell to compute the normalized directional vector of each 1D cell in 3D space, and then compute the ‘directional derivative’ as follows:

direction_vector = cell_node1 - cell_node2
magnitude = np.sqrt(np.dot(direction_vector, direction_vector))
direction_vector /= magnitude if magnitude != 0 else direction_vector
direction_vector_ufl = as_vector(direction_vector)
directional_derivative = dot(grad(u), direction_vector_ufl)

Would the following be the correct way to compute the direction vector of each 1D cell in 3D space?

x = ufl.SpatialCoordinate(domain)
xMag = sqrt(x[0]**2 + x[1]**2 + x[2]**2)
direction_vector_ufl = as_vector([x[0]/xMag, x[1]/xMag, x[2]/xMag])

BTW, how to extract the values of the cells’ direction vectors so that I can double-check the implementation? I tried the following but got errors:

print(f'---> dirVector: {direction_vector_ufl}')
print(f'---> dirVector.vector.array: {direction_vector_ufl.vector.array}')

------------ Here is the error message:
—> dirVector: [x[0] / sqrt(x[0] ** 2 + x[1] ** 2 + x[2] ** 2), x[1] / sqrt(x[0] ** 2 + x[1] ** 2 + x[2] ** 2), x[2] / sqrt(x[0] ** 2 + x[1] ** 2 + x[2] ** 2)]
Traceback (most recent call last):
File “/home/whf/FEniCSx/mixedElement_TODO/mixedElem_1v_TODO.py”, line 427, in
print(f’—> dirVector.vector.array: {direction_vector_ufl.vector.array}')
AttributeError: ‘ListTensor’ object has no attribute ‘vector’

4

I tried to emphasize that the directional derivative is a scalar but just now realized that this might not clearly answer your question. To clarify, I only compute a directional derivative but still integrate it in the 3D domain. Thank you!