How to compute differential of a complex form?

[hawkspar]

Hello Fenics users and developers,

I’m trying to use a fenics code to reproduce the results of this paper as part of a validation process.

This is a fluid mechanics paper that involves a specific decomposition : q(x,r,\theta,t)=q_0(x,r)+\epsilon\hat{q}_1(x,r)e^{im\theta+(\sigma+i\omega)t}

Going from the Navier-Stokes equations N\partial_tq+M(q)=0 and obtaining a base flow q_0 from M(q_0)=0, this allows for computation of the perturbations as a generalised eigenvalue problem dM/dq(q_0)\hat{q}_1=(\sigma+i\omega)N\hat{q}_1

The nonlinear operator M is well–known in the case of the Navier Stokes equations but in this decomposition it is complex (\hat{q}_1 is complex as well).

I’ve read here that the UFL can handle complex values. The code below is an attempt to leverage automatic differentiation to obtain dM/dq(q_0).

from dolfin import *
from pdb import set_trace

# Trivial geometry
mesh = UnitSquareMesh(2,2)
r = SpatialCoordinate(mesh)[1]
# Physical parameters
m=-1
Re=200
# Taylor Hodd elements ; stable element pair
FE_vector=VectorElement("Lagrange",mesh.ufl_cell(),2,3)
FE_scalar=FiniteElement("Lagrange",mesh.ufl_cell(),1)
V=FunctionSpace(mesh,MixedElement([FE_vector,FE_scalar]))
# Test and trial functions
test = TestFunction(V)
trial = TrialFunction(V)
# Initial conditions
q = interpolate(Constant([1,0,0,0]), V)

# Gradient with x[0] is x, x[1] is r, x[2] is theta
def grad(v,m):
	return as_tensor([[v[0].dx(0), v[0].dx(1), m*1j*v[0]/r],
				 	  [v[1].dx(0), v[1].dx(1), m*1j*v[1]/r-v[2]/r],
					  [v[2].dx(0), v[2].dx(1), m*1j*v[2]/r+v[1]/r]])
def grad_re(v):
	return as_tensor([[v[0].dx(0), v[0].dx(1),  0],
				 	  [v[1].dx(0), v[1].dx(1), -v[2]/r],
					  [v[2].dx(0), v[2].dx(1),  v[1]/r]])

def div(v,m): return v[0].dx(0) + (r*v[1]).dx(1)/r + m*1j*v[2]/r

# Navier Stokes variational form
def NS_form(m):
    u,p=split(q)
    test_u,test_m=split(test)

    #mass (variational formulation)
    M = div(u,m)*test_m*r*dx
    #set_trace()
    #momentum (different test functions and IBP)
    M += 	    dot(grad(u,m)*u,	test_u)   *r*dx # Convection
    M += 1/Re*inner(grad(u,m), grad(test_u,m))*r*dx # Diffusion
    M -= 	    dot(p, 			div(test_u,m))*r*dx # Pressure
    return M
    
# Compute baseflow
solve(NS_form(0) == 0, q, [])

from ufl.algorithms.compute_form_data import compute_form_data

pert_form=compute_form_data(NS_form(m),complex_mode=True)
parameters['linear_algebra_backend'] = 'Eigen'
L = as_backend_type(assemble(derivative(pert_form,q,trial))).sparray()

With my installation, dolfin 2019.1.0 installed using docker I obtain the resulting traceback :

Traceback (most recent call last):
  File "MWE.py", line 51, in <module>
    pert_form=compute_form_data(NS_form(m),complex_mode=True)
  File "/usr/local/lib/python3.6/dist-packages/ufl/algorithms/compute_form_data.py", line 418, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)  # Currently testing how fast this is
  File "/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py", line 177, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments, complex_mode)
  File "/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py", line 159, in check_integrand_arity
    arg_tuples = map_expr_dag(rules, expr, compress=False)
  File "/usr/local/lib/python3.6/dist-packages/ufl/corealg/map_dag.py", line 37, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress)
  File "/usr/local/lib/python3.6/dist-packages/ufl/corealg/map_dag.py", line 86, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py", line 48, in sum
    raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments ('v_0',) vs ('conj(v_0)',).

I would be grateful on any insight as to why this error arises, how to go around it, or more generally how to handle complex values in UFL (I’ve found few examples) and how to solve my problem.

Cheers,

[dokken]

Dolfin does not handle complex values. It is handles in dolfinx .

However i seem to recall that there is an issue with using derivative on a complex valued form. However i couldnt find where I read this.

It might be worthwhile creating a minimal example (using ufl syntax only) and submit an issue to Issues · FEniCS/ufl · GitHub

[hawkspar]

Thank you for your quick response.

I have tried to shift my whole code from dolfin to dolfinx, but there are tons of subtle changes everywhere (file format change, importing meshes work differently, boundary conditions are different too) that doused my enthusiasm.

Even if it is not possible to directly handle complex values directly in dolfin - and I for one am rather confused about where dolfin ends and where UFL starts, since UFL seems to have a direct python API too - don’t you think there should be a way around this to address my problem in dolfin ?

I understand there could be an issue with derivating a complex form automatically, though naively I would have thought as long as the complex form is differentiable, symbolic derivative on such a complex form really works the same as differentiating two real valued functions, namely the real and imaginary part of the operator.

I know that is not the case in the above example, but I have tried to split real and imaginary part and differentiate both separately. In fact I have tried a number of things before posting here :

1. Split real and imaginary part of the operator
   a. Pretending to forget its argument is complex and differentiate : gives results that are different from reference paper
   b. Double the size of q to include imaginary parts then try to differentiate at q_0 which is real : only I can’t seem to force differentiation along the required direction
2. Double the size of everything ; double the size of q and that of the form, handle only reals and recompose the actual complex matrix after a bit of as_backend_type(assemble(2n_form)).sparray() magic : can’t seem to get fenics to produce the 2n form
3. Convert everything to dolfinx : even importing the .msh mesh is a nightmare

I’m getting pretty desperate and I feel this whole thing shouldn’t be this hard. Am I missing anything ?

[dokken]

@hawkspar

I have tried to shift my whole code from dolfin to dolfinx , but there are tons of subtle changes everywhere (file format change, importing meshes work differently, boundary conditions are different too) that doused my enthusiasm.

I have made a tutorial that shows how to use dolfinx: The FEniCSx tutorial — FEniCSx tutorial

@hawkspar

Even if it is not possible to directly handle complex values directly in dolfin - and I for one am rather confused about where dolfin ends and where UFL starts, since UFL seems to have a direct python API too - don’t you think there should be a way around this to address my problem in dolfin ?

Dolfin uses UFL to generate integration kernels. UFL supports complex variables in forms (on the form a+j*b), but the kernel generator in dolfin does not. Therefore, to use actual complex numbers, you need to change to dolfin-x.

One thing i missed in your previous code is that you use the dot operator, when using complex numbers, you cannot use the dot operation when multiplying with a test function; you need to use the inner function.

I’ve added a minimal code example that uses dolfinx, in complex mode:

docker run -ti -e DISPLAY=$DISPLAY -v /tmp/.X11-unix:/tmp/.X11-unix -v $(pwd):/root/shared -w /root/shared --rm dolfinx/dolfinx
source dolfinx-complex-mode

import ufl
import dolfinx
from mpi4py import MPI


cell = ufl.triangle
coordinate_element = ufl.VectorElement("Lagrange", cell, 1)
element = ufl.FiniteElement("Lagrange", cell, 1)
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 10, 10)
V = dolfinx.FunctionSpace(mesh, element)

u = dolfinx.Function(V)
u.x.array[:] = 1
v = ufl.TestFunction(V)
a = u*ufl.inner(u, v)*ufl.dx


da = ufl.derivative(a, u)
A = dolfinx.fem.assemble_matrix(da)
A.assemble()
import numpy as np
np.set_printoptions(precision=5)

print(A[:, :])

@hawkspar

Split real and imaginary part of the operator
a. Pretending to forget its argument is complex and differentiate : gives results that are different from reference paper
b. Double the size of q qq to include imaginary parts then try to differentiate at q_0 q0q_0 which is real : only I can’t seem to force differentiation along the required direction

There are plenty of posts on how to split a complex problem into two real problems.

@hawkspar

Convert everything to dolfinx : even importing the .msh mesh is a nightmare

This is for instance covered in
https://jorgensd.github.io/dolfinx-tutorial/chapter3/subdomains.html#convert-msh-files-to-xdmf-using-meshio

[hawkspar]

@dokken

Dolfin uses UFL to generate integration kernels. UFL supports complex variables in forms (on the form a+j*b), but the kernel generator in dolfin does not. Therefore, to use actual complex numbers, you need to change to dolfin-x.

Thank you for the clarification.

@dokken

One thing i missed in your previous code is that you use the dot operator, when using complex numbers, you cannot use the dot operation when multiplying with a test function; you need to use the inner function.

That’s what the UFL developers told me, and indeed it fixes the original traceback, but fails at the assemble stage - no doubt because of the kernel problem.

@dokken

There are plenty of posts on how to split a complex problem into two real problems.

I’m sorry if my question is redundant, I posted it here after combing through documentation and Q&A for some time. I had seen and implemented the snippet on .msh meshes in dolfin-x - it mostly convinced me transition would not be painless.

I guess it’s now clear there is no easy way around this problem without going to dolfin-x. I have to admit it is very frustrating to have a differentiated UFL complex form and stuck at changing all the rest of the code to get the assembly part right.

@dokken

I have made a tutorial that shows how to use dolfinx: The FEniCSx tutorial — FEniCSx tutorial

I had not seen this specific tutorial, it will serve nicely as a starting point. Thanks for your time.