Error using Non-Linear Poisson in Complex Fenicsx. Arity Mismatch

General
1

Hi,
I have a Non-Linear Poisson Solver in the complex Fenicsx type. I get an error when using the Non-Linear function. My code is

’
def rho_pointwise(self, phi):

    # Constants
    kB = 1.380649e-23  # [J/K]
    T = 300.0          # [K]
    KT = kB * T / 1.6e-19  # [eV] 

    kbT = fem.Constant(self.domain, PETSc.ScalarType(KT - 0j))
    qval = fem.Constant(self.domain, PETSc.ScalarType(1.6e-19 - 0j))

    # Take real parts for physical calculations
    phi_real = ufl.real(phi)
    chi_real = ufl.real(self.chi)
    Eg_real = ufl.real(self.Eg)
    Nc_real = ufl.real(self.Nc)
    Nv_real = ufl.real(self.Nv)
    phiref_real = ufl.real(self.phiref)
    Ef_real = ufl.real(self.Ef)
    kbT_real = ufl.real(kbT)
    
    # Compute band edges (now real-valued)
    Ec = -(phi_real - phiref_real + chi_real)  # conduction band edge [eV]
    Ev = Ec - Eg_real                          # valence band edge [eV]
    
    # Fermi level relative to bands
    xe = (Ef_real - Ec) / kbT_real
    xh = (Ev - Ef_real) / kbT_real
    
    def fermi_approx_smooth(x):

        exp_part = ufl.exp(ufl.min_value(x, 50.0))  # Clamped exponential
        power_part = 0.7523 * (ufl.sqrt((x**2 + 1e-6)))**1.5  # Always positive argument
        
        # Transition around x = 2.5
        transition_point = 2.5
        transition_width = 2.0
        
        # Sigmoid function for smooth blending: 0 at -inf, 1 at +inf
        blend = 0.5 * (1.0 + ufl.tanh((x - transition_point) / transition_width))
        
        return (1.0 - blend) * exp_part + blend * power_part
    
    # Carrier densities (all real-valued)
    nval = Nc_real * fermi_approx_smooth(xe)
    pval = Nv_real * fermi_approx_smooth(xh)
    
    # Charge density (real-valued expression)
    rho = ufl.real(qval) * (pval - nval)
    return rho



def nonlin_solve(self, gate_voltages, ohmic_gate_name, save_file=True, max_iterations=10, tol=1e-9):

    # Create function space and functions
    V = fem.functionspace(self.domain, ('Lagrange', 1))
    phi = fem.Function(V, dtype=np.complex128)  # Solution (complex-valued)
    # Define variational problem
    v = ufl.TestFunction(V)      # Test function
    u = ufl.TrialFunction(V)     # Trial function (for both linear and Newton increment)
    tdim = self.domain.topology.dim
    fdim = tdim - 1

    bcs = []
    
    for gate_name, voltage in gate_voltages.items():
        if gate_name not in self.device.facet_tags_map:
            continue
        
        tag = self.device.facet_tags_map[gate_name]
        dofs = fem.locate_dofs_topological(V, fdim, self.device.facet_tags.find(tag)) 
        if len(dofs) == 0:
            continue
        
        phi_v = voltage - (4.33 - 2.7)  
        bc = fem.dirichletbc(np.complex128(phi_v + 0j), dofs, V)
        bcs.append(bc)
        if self.domain.comm.rank == 0:
            print(f"{gate_name:<20} {phi_v:>10.3f} V {len(dofs):>8}")

    phi.x.array[:] = 0.0
        
    F = self.eps*ufl.inner(ufl.grad(phi), ufl.grad(v))*ufl.dx - self.rho_pointwise(phi)*v*ufl.dx
        
    J = ufl.derivative(F, phi, ufl.conj(v))

    problem = NonlinearProblem(F, phi, bcs=bcs, J=J, petsc_options_prefix="nonlinear_poisson_")

    solver = NewtonSolver(MPI.COMM_WORLD, problem)

’

The error that i get is
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
self._F = _create_form(
raise ArityMismatch(
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments (‘v_0’,) vs (‘conj(v_0)’,).
Exception ignored in: <function NonlinearProblem.del at 0x146e13735800>
Traceback (most recent call last):
lambda obj: obj is not None, (self._snes, self._A, self._b, self._x, self._P_mat)
^^^^^^^^^^
AttributeError: ‘NonlinearProblem’ object has no attribute ‘_snes’

Thank you

2

I don’t have time to check your code, but your formulation looks a bit strange:

    F = self.eps*ufl.inner(ufl.grad(phi), ufl.grad(v))*ufl.dx - self.rho_pointwise(phi)*v*ufl.dx
        
    J = ufl.derivative(F, phi, ufl.conj(v))

Bear in mind that

(a, b) = \int_\Omega a \cdot \overline{b} \; \mathrm{d}\mathbf{x}

So I would be expecting something like

    F = ufl.inner(self.eps*ufl.grad(phi), ufl.grad(v))*ufl.dx - ufl.inner(self.rho_pointwise(phi), v)*ufl.dx
    J = ufl.derivative(F, phi, v)
3

Sorry, yes it is
F = ufl.inner(self.eps*ufl.grad(phi), ufl.grad(v))*ufl.dx - ufl.inner(self.rho_pointwise(phi), v)*ufl.dx
J = ufl.derivative(F, phi, v)

But even with these changes I get the same error.

1 Like
4

Hi,

So it looks like the earlier issue was caused by importing NewtonSolver from the wrong module. Switching to

from dolfinx.fem.petsc import NonlinearProblem

resolved the ArityMismatch error.

However, I’m now running into a problem when I use any ufl.conditional (or ufl.lt/le/min/max ) with complex-valued UFL expressions. I tried converting the expression to its real part for the comparison, but the generated FFCx code still attempts to evaluate the comparison on a complex value and fails at compile time.

I’ve attached my updated code below. The error I see is:

libffcx_forms_...c:880:32: error: invalid operands to binary < (have ‘complex double’ and ‘double’)
  bool sv_... = sv_... < 5.0;

Code snippet:

from petsc4py import PETSc
from dolfinx.io import XDMFFile
from dolfinx import fem, io, nls, log
from dolfinx.fem.petsc import NonlinearProblem

V = fem.functionspace(domain, ('Lagrange', 1))
phi = fem.Function(V, dtype=np.complex128)
dphi = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

def rho_test(phi):
    kB = 1.380649e-23
    T  = 10.0
    KT = kB*T/1.6e-19
    kbT = KT
    q   = 1.6e-19

    Ec = -(phi - phiref + chi)
    Ev = Ec - Eg

    xe = (Ef - Ec) / kbT
    xh = (Ev - Ef) / kbT

    def cond_real(x):
        x_real = ufl.real(x)
        return ufl.conditional(ufl.lt(x_real, 5.0),
                               ufl.exp(x_real),
                               0.7523 * (x_real**1.0))

    Xe_real = cond_real(xe)
    Xp_real = cond_real(xh)

    n = Nc * Xe_real
    p = Nv * Xp_real
    return q * (p - n)

F = ufl.inner(eps*ufl.grad(phi), ufl.grad(v))*ufl.dx - rho_test(phi)*ufl.conj(v)*ufl.dx
J = ufl.derivative(F, phi, dphi)

petsc_options = {
    "snes_type": "newtonls",
    "snes_linesearch_type": "none",
    "snes_atol": 1e-6,
    "snes_rtol": 1e-5,
    "snes_monitor": None,
    "ksp_error_if_not_converged": True,
    "ksp_type": "gmres",
    "ksp_rtol": 1e-8,
    "ksp_monitor": None,
    "pc_type": "hypre",
    "pc_hypre_type": "boomeramg",
    "pc_hypre_boomeramg_max_iter": 1,
    "pc_hypre_boomeramg_cycle_type": "v",
}

problem = NonlinearProblem(F, phi, J=J, bcs=bcs,
                           petsc_options=petsc_options,
                           petsc_options_prefix="nonlinear_poisson_")
phi = problem.solve()

Could you please advise whether there is a supported approach for implementing conditionals in complex mode or an alternative strategy to avoid this compilation error or if I am missing something entirely?

Thank you,
Kelvin

1 Like