Cahn Hilliard solver problem dolfinx

dolfinx
1

Hello

I’m trying to write a script in line with this example for dolfinx.

I get an error when I invoke dolfinx.fem.assemble_matrix :

raise ArityMismatch(“Failure to conjugate test function in complex Form”)
ufl.algorithms.check_arities.ArityMismatch: Failure to conjugate test function in complex Form

Any ideas on what can be done?

2

It’s difficult to debug without a MWE. But have you not correctly used inner to signify an inner product in your FE space?

3

I have not used inner. This variational form used to work fine in dolfin. Has this been changed in dolfinx?

MWE below (without mesh):

from ufl import *
from dolfinx import *
import matplotlib.pyplot as plt
from dolfinx.plotting import plot
from dolfinx.mesh import locate_entities_boundary
from dolfinx.fem import locate_dofs_geometrical, locate_dofs_topological, assemble_matrix
from dolfinx.fem.assemble import set_bc, assemble_vector, apply_lifting
import numpy as np
import time, sys

class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L):
        super().__init__()
        self.L = Form(L)
        self.a = Form(a)
        self._F = None
        self._J = None

    def form(self, x):
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

    def F(self, x):
        if self._F is None:
            self._F = assemble_vector(self.L)
        else:
            with self._F.localForm() as f_local:
                f_local.set(0.0)
            self._F = assemble_vector(self._F, self.L)
        self._F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        return self._F

    def J(self, x):
        if self._J is None:
            self._J = assemble_matrix(self.a)
        else:
            self._J.zeroEntries()
            self._J = assemble_matrix(self._J, self.a)
        self._J.assemble()
        return self._J

# Load mesh
mesh = getDolfinx2DMesh('attempt2')

# Build function space on mesh
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
V = FunctionSpace(mesh, TH)
V0 = V.sub(0).collapse()

# Label boundaries 
def inflow(x):
    return np.isclose(x[0],-110.0)
def outflow(x):
    return np.isclose(x[0],170.0)
def walls(x):
    return np.logical_or(np.isclose(x[0],-60.0),np.isclose(x[0],60.0))
def cylinder(x):
    return np.isclose(x[0]*x[0] + x[1]*x[1],0.25)

# Define boundary conditions
def inflow_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))
def walls_expression(x):
    return np.stack((np.zeros(x.shape[1])))
def cylinder_expression(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1])))
def outflow_expression(x):
    return np.stack((np.zeros(x.shape[1])))

inflow_velocity = Function(V0)
facets = locate_entities_boundary(mesh,1,inflow)
dofs = locate_dofs_topological(V.sub(0),1,facets)
bc1 = DirichletBC(inflow_velocity,dofs)

walls_velocity = Function(V0)
facets = locate_entities_boundary(mesh,1,walls)
dofs = locate_dofs_topological(V.sub(0).sub(1),1,facets)
bc2 = DirichletBC(walls_velocity,dofs)

cylinder_velocity = Function(V0)
facets = locate_entities_boundary(mesh,1,cylinder)
dofs = locate_dofs_topological(V.sub(0),1,facets)
bc3 = DirichletBC(cylinder_velocity,dofs)


outflow_pressure = Function(V0)
facets = locate_entities_boundary(mesh,1,outflow)
dofs = locate_dofs_topological(V.sub(1),1,facets)
bc4 = DirichletBC(outflow_pressure,dofs)

bcs = [bc1,bc2,bc3,bc4]

# Declare test and trial functions
up = Function(V)
(u,p) = up.split()
(v,q) = TestFunctions(V)

Re = Constant(mesh,50.)
i,j = indices(2)

# Variational form
a = q*div(u)*dx - p*div(v)*dx + (1.0/Re)*(Dx(u[i],j)*Dx(v[i],j))*dx \
    + dot(dot(u,nabla_grad(u)),v)*dx 
b1 = 0

J = derivative(a,up)
problem = CahnHilliardEquation(J, a)
solver = NewtonSolver(MPI.COMM_WORLD)
solver.convergence_criterion = "residual"
solver.linear_solver = "mumps"
solver.rtol = 1e-6
solver.solve()

I’d appreciate any inputs!