How to compute differential of a complex form?

[hawkspar]

Hello UFL users and developers,

This post is taken from the fenics forum with an adapted minimal working example.

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 q_1(x,r)e^(imtheta+(sigma+i omega)t) Going from the Navier-Stokes equations Nd_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)q_1=(sigma+i omega)Nq_1

The nonlinear operator M is well--known in the case of the Navier Stokes equations but in this decomposition it is complex (q_1 is complex as well). The core of the problem is computing dM/dq.

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 ufl import *
from ufl.classes import TestFunction, TrialFunction
from ufl.algorithms.compute_form_data import compute_form_data

# Physical parameters
m=-1
Re=200
cell = triangle
FE_vector=VectorElement("Lagrange",cell,2,3)
FE_scalar=FiniteElement("Lagrange",cell,1)
FE_TH = MixedElement([FE_vector,FE_scalar])
# Test and trial functions
test = TestFunction(FE_TH)
trial = TrialFunction(FE_TH)
# Initial conditions
q = Coefficient(FE_TH)

# 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]],
                      [v[1].dx(0), v[1].dx(1), m*1j*v[1]-v[2]],
                      [v[2].dx(0), v[2].dx(1), m*1j*v[2]+v[1]]])

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

# 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*dx
    #set_trace()
    #momentum (different test functions and IBP)
    M +=        dot(grad(u,m)*u,    test_u)   *dx # Convection
    M += 1/Re*inner(grad(u,m), grad(test_u,m))*dx # Diffusion
    M -=        dot(p,          div(test_u,m))*dx # Pressure
    return M

pert_form=compute_form_data(NS_form(m),complex_mode=True)

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

Traceback (most recent call last):
  File "MWE_ufl.py", line 46, 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.

[jorgensd]

The issue here is that you need to use inner whenever you multiply by a test function, i.e.

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

    #mass (variational formulation)
    M = inner(div(u,m),test_m)*dx
    #set_trace()
    #momentum (different test functions and IBP)
    M +=       inner(grad(u,m)*u,   test_u)   *dx # Convection
    M += 1/Re*inner(grad(u,m), grad(test_u,m))*dx # Diffusion
    M -=        inner(p,            div(test_u,m))*dx # Pressure
    return M

pert_form=compute_form_data(NS_form(m),complex_mode=True)