dolfinx

test.py

#!/usr/bin/env python

# https://fenicsproject.discourse.group/t/sudden-arity-mismatch-when-defining-non-linear-problem/8844/11

from mpi4py import MPI
from dolfinx import fem, mesh
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
import matplotlib.pyplot as plt


L = 4.0e-3  # depth of tray
t_step = 60
D = 6.5e-06  # diffusivity (m2/s)
k = 6e-2  # rate constant (1/s)
C_0 = 0.00133 * 1000  # atmospheric CO2 (mol/m3)
N_0 = 4.018 * 1000  # initial molar density of Ca(OH)2 (mol/m3)

domain = mesh.create_interval(MPI.COMM_WORLD, 100, [0, L])
P1 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
element = ufl.MixedElement([P1, P1])
V = fem.FunctionSpace(domain, element)

# test functions
v1, v2 = ufl.TestFunctions(V)

# fields to solve for
num_subs = V.num_sub_spaces
spaces = []
maps = []
for i in range(num_subs):
    space_i, map_i = V.sub(i).collapse()
    spaces.append(space_i)
    maps.append(map_i)

u = fem.Function(V)
u1 = fem.Function(spaces[0])
u2 = fem.Function(spaces[1])

# concentrations from previous timestep
u_n = fem.Function(V)
u_n1, u_n2 = ufl.split(u_n)

# apply timestep variable to each node in the mesh
dt = fem.Constant(domain, ScalarType(t_step))

# turn floats to Constant objects
D_ = fem.Constant(domain, ScalarType(D))
k_ = fem.Constant(domain, ScalarType(k))

# set initial conditions
u1.vector.array[:] = C_0 * np.ones(u1.vector.array.shape)
u2.vector.array[:] = N_0 * np.ones(u2.vector.array.shape)
u_n.vector.array[maps[0]] = u1.vector.array
u_n.vector.array[maps[1]] = u2.vector.array

# set initial guess
u.vector.array[maps[0]] = u1.vector.array
u.vector.array[maps[1]] = u2.vector.array

u1, u2 = ufl.split(u)

# function determining if a node is on the tray top
def on_top_boundary(x):
    return np.isclose(x[0], L)


# must use collapsed subspace to determine boundary DOFs
V0, submap = V.sub(0).collapse()

# determine boundary DOFs
boundary_dofs = fem.locate_dofs_geometrical((V.sub(0), V0), on_top_boundary)

# apply dirichlet BC to boundary DOFs
bc = fem.dirichletbc(ScalarType(C_0), boundary_dofs[0], V.sub(0))

u.x.array[bc.dof_indices()[0]] = bc.value.value

F = (
    (1 / dt) * (u1 - u_n1) * v1 * ufl.dx
    + (1 / dt) * (u2 - u_n2) * v2 * ufl.dx
    + D_ * ufl.inner(ufl.grad(u1), ufl.grad(v1)) * ufl.dx
    + k_ * u1 * u2 * v1 * ufl.dx
    + k_ * u1 * u2 * v2 * ufl.dx
)

problem = fem.petsc.NonlinearProblem(F, u, bcs=[bc])

test1.py

#!/usr/bin/env python

import dolfinx
from mpi4py import MPI
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
u_r = dolfinx.fem.Function(V, dtype=np.float64)
u_r.interpolate(lambda x: x[0])
u_c = dolfinx.fem.Function(V, dtype=np.complex128)
u_c.interpolate(lambda x: 0.5 * x[0] ** 2 + 1j * x[1] ** 2)
print(u_r.x.array)
print(u_c.x.array)

test2.py

#!/usr/bin/env python

from petsc4py import PETSc
import numpy as np

print(PETSc.ScalarType)
assert np.dtype(PETSc.ScalarType).kind == "c"

test3.py

#!/usr/bin/env python

import ufl
import dolfinx
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = dolfinx.fem.Constant(mesh, PETSc.ScalarType(-1 - 2j))
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = ufl.inner(f, v) * ufl.dx

L2 = f * ufl.conj(v) * ufl.dx
print(L)
print(L2)

u_c = dolfinx.fem.Function(V, dtype=np.complex128)
J = u_c**2 * ufl.dx
F = ufl.derivative(J, u_c, ufl.conj(v))
residual = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(F))
print(residual.array)