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