Arity Mismatch error in mixed space with FENICS

Errors
1

Dear. Researchers
I am trying to solve coupled time dependent differential equations in Fenics.
After running the program I have received the following error. I would be really appreciated if someone assist me to resolve this issue.
Adding expressions with non-matching form arguments (‘v_0’,) vs ()

Here is my program:

import random
from dolfin import *
import ufl

dt = 5.0e-06 # time step
theta = 0.5 # time stepping family, e.g. theta=1 → backward Euler, theta=0.5 → Crank-Nicolson

mesh = IntervalMesh(500,0,100)
P1 = FiniteElement(“Lagrange”, interval, 3)
ME = FunctionSpace(mesh,MixedElement(P1, P1, P1, P1, P1))
mesh_file = File(“Mesh1D.xml”)
mesh_file << mesh

Define functions

u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
du = TrialFunction(ME)

class K(UserExpression):

  def eval(self, values, x):
      left_x = 0.0
      right_x = 100    
      # int_width = 0.05
      tol = 1e-14

      if (x[0] <= 30 + tol):
          values[0] = 0.0
          values[1] = 0.0
          values[2] = 0.0
          values[3] = 0.6
          values[4] = 0.0
      
      elif(x[0] >= 30 + tol):
          values[0] = 1.0
          values[1] = 0.0
          values[2] = 0.4
          values[3] = 0.0
          values[4] = 0.0
           
      
  def value_shape(self):
      return(5, )

order_value = K(mesh)

u0 = interpolate(order_value, ME)
u.assign(u0)

u = interpolate(order_value, ME)

(v1 ,v2 ,v3 , v4, v5) = TestFunctions(ME)

(u1, u2, u3, u4, u5) = split(u)

(r0, s0, uu0, v0, w0) = split(u0)

define swtiching function

u1 = variable(u1)
hfunction1 = (u13) (6(u12) -15*u1 + 10)
hder1 = diff(hfunction1,u1)

hfunction2 = ((1-u1)**3)(6((1-u1)*2) -15(1-u1) + 10)

wg = Constant(3)
kappa = Constant(1.5)
Lp = Constant(2)

g = wg*(u1**2)*(1-u1)**2

gder = diff(g,u1)

#formation the weak form of Eq number 8

#free energy of phase a

u3 = variable(u3)
f_a = 100*(u3-0.4)**2 - 10

fder_a = diff(f_a, u3)

#define free energy of phase b
u4 = variable(u4)
f_b = 100*(u4-0.6)**2
fder_b = diff(f_b, u4)

Form_eta = (u1v1dx - r0v1dx + dtLp(hder1*(f_a - f_b) +
fder_a* hder1 (u4 - u3))v1dx + dtLp*(wggder +
0.5kappa*dot(grad(u1), grad(v1)))*dx)

Form_cons = u5v2dx - fder_av2dx

Form_mu = (u2v3dx - s0v3dx + (dtdot(grad(u5), grad(v3))(10*
(1 - hfunction1) + hfunction1))*dx)

Form_xia = fder_av4dx - fder_bv4dx
Form_xib = u2v5dx - hfunction1u3v5dx - (1-hfunction1)u4v5dx

Form = Form_eta + Form_cons + Form_mu + Form_xia + Form_xib

Gain = derivative(Form, u, du)

Step in time

t = 0.0
T = 50*dt
while (t < T):
t += dt
solve(Form == 0, u)
u0.assign(u)

2

You need to properly format your code with 3x` encapsulation,
i.e.

```python
#add code

```
3

#!/usr/bin/env python3

-- coding: utf-8 --

“”"
Created on Fri Feb 21 14:58:17 2025

@author: mohammad
“”"


import random
from dolfin import *
import ufl 


dt     = 5.0e-06  # time step
theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson


mesh = IntervalMesh(500,0,100)
P1 = FiniteElement("Lagrange", interval, 3)
ME = FunctionSpace(mesh,MixedElement(P1, P1, P1, P1, P1))
mesh_file = File("Mesh1D.xml")
mesh_file << mesh


# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step
du    = TrialFunction(ME)


class K(UserExpression):
    
    
      
      def eval(self, values, x):
          left_x = 0.0
          right_x = 100    
          # int_width = 0.05
          tol = 1e-14

          if (x[0] <= 30 + tol):
              values[0] = 0.0
              values[1] = 0.0
              values[2] = 0.0
              values[3] = 0.6
              values[4] = 0.0
          
          elif(x[0] >= 30 + tol):
              values[0] = 1.0
              values[1] = 0.0
              values[2] = 0.4
              values[3] = 0.0
              values[4] = 0.0
               
          
      def value_shape(self):
          return(5, )
   
        
order_value = K(mesh)



u = interpolate(order_value, ME)
u0.assign(u)

(v1 ,v2 ,v3 , v4, v5) = TestFunctions(ME)



(u1, u2, u3, u4, u5) = split(u)


(r1, s1, uu1, vv1, w1) = split(u0)

(du1, du2, du3, du4, du5) = split(du)


# define swtiching function

u1 = variable(u1)
hfunction1 = (u1**3) *(6*(u1**2) -15*u1 + 10)
hder1 = diff(hfunction1,u1)


hfunction2 = ((1-u1)**3)*(6*((1-u1)**2) -15*(1-u1) + 10)


wg = Constant(3)
kappa = Constant(1.5)
Lp = Constant(2)


g = wg*(u1**2)*(1-u1)**2

gder = diff(g,u1)



#formation the weak form of Eq number 8

#free energy of phase a

u3 = variable(u3)
f_a = 100*(u3-0.4)**2 - 10

fder_a = diff(f_a, u3)




#define free energy of phase b
u4 = variable(u4)
f_b = 100*(u4-0.6)**2
fder_b = diff(f_b, u4)


Form_eta = (u1*v1*dx - r1*v1*dx + dt*Lp*(hder1*(f_a - f_b) +
                                         fder_a* hder1 *(u4 - u3))*v1*dx +
            dt*Lp*(wg*gder + 
            0.5*kappa*dot(grad(u1), grad(v1)))*dx)
                                        
  
                                        
                                        
Form_cons = u5*v2*dx - fder_a*v2*dx


Form_mu = (u2*v3*dx - s1*v3*dx + (dt*dot(grad(u5), grad(v3))*(10*
(1 - hfunction1) + hfunction1))*dx)


Form_xia = fder_a*v4*dx - fder_b*v4*dx
Form_xib = u2*v5*dx - hfunction1*u3*v5*dx - (1-hfunction1)*u4*v5*dx

Form = Form_eta + Form_cons + Form_mu + Form_xia + Form_xib


Gain = derivative(Form, u, du)



# Step in time
t = 0.0
T = 50*dt
while (t < T):
    t += dt
    solve(Form == 0, u, J = Gain)
    u0.assign(u)
4

Finally I managed to debug the code. The error stemmed from the existence of terms with arity zero in the weak form.