I have been solving the following problem using FEniCSx for over a month now without any issues. All of a sudden, I am getting ArityMismatch: Adding expressions with non-matching form arguments ('v_0',) vs ('conj(v_0)',) when I define the non-linear problem. I do not believe I have changed anything leading up to the problem definition. Here is the code:
L = 4.e-3
t_step = 60
D = 6.5e-06
k = 6e-2
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])
When I get to this step, it returns:
ArityMismatch Traceback (most recent call last)
Input In [8], in <cell line: 7>()
1 F = (1/dt)*(u1 - u_n1)*v1*ufl.dx + \
2 (1/dt)*(u2 - u_n2)*v2*ufl.dx + \
3 D_*ufl.inner(ufl.grad(u1), ufl.grad(v1))*ufl.dx + \
4 k_*u1*u2*v1*ufl.dx + \
5 k_*u1*u2*v2*ufl.dx
----> 7 problem = fem.petsc.NonlinearProblem(F, u, bcs=[bc])
8 solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/petsc.py:678, in NonlinearProblem.__init__(self, F, u, bcs, J, form_compiler_params, jit_params)
656 def __init__(self, F: ufl.form.Form, u: _Function, bcs: typing.List[DirichletBCMetaClass] = [],
657 J: ufl.form.Form = None, form_compiler_params={}, jit_params={}):
658 """Initialize solver for solving a non-linear problem using Newton's method, :math:`dF/du(u) du = -F(u)`.
659
660 Args:
(...)
676
677 """
--> 678 self._L = _create_form(F, form_compiler_params=form_compiler_params,
679 jit_params=jit_params)
681 # Create the Jacobian matrix, dF/du
682 if J is None:
File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py:168, in form(form, dtype, form_compiler_params, jit_params)
165 return list(map(lambda sub_form: _create_form(sub_form), form))
166 return form
--> 168 return _create_form(form)
File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py:163, in form.<locals>._create_form(form)
160 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
161 return form argument"""
162 if isinstance(form, ufl.Form):
--> 163 return _form(form)
164 elif isinstance(form, collections.abc.Iterable):
165 return list(map(lambda sub_form: _create_form(sub_form), form))
File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py:137, in form.<locals>._form(form)
134 if mesh is None:
135 raise RuntimeError("Expecting to find a Mesh in the form.")
--> 137 ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
138 form_compiler_params=form_compiler_params,
139 jit_params=jit_params)
141 # For each argument in form extract its function space
142 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]
File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/jit.py:56, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
51 @functools.wraps(local_jit)
52 def mpi_jit(comm, *args, **kwargs):
53
54 # Just call JIT compiler when running in serial
55 if comm.size == 1:
---> 56 return local_jit(*args, **kwargs)
58 # Default status (0 == ok, 1 == fail)
59 status = 0
File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/jit.py:204, in ffcx_jit(ufl_object, form_compiler_params, jit_params)
202 # Switch on type and compile, returning cffi object
203 if isinstance(ufl_object, ufl.Form):
--> 204 r = ffcx.codegeneration.jit.compile_forms([ufl_object], parameters=p_ffcx, **p_jit)
205 elif isinstance(ufl_object, ufl.FiniteElementBase):
206 r = ffcx.codegeneration.jit.compile_elements([ufl_object], parameters=p_ffcx, **p_jit)
File /usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py:168, in compile_forms(forms, parameters, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
165 for name in form_names:
166 decl += form_template.format(name=name)
--> 168 impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
169 cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
170 except Exception:
171 # remove c file so that it will not timeout next time
172 c_filename = cache_dir.joinpath(module_name + ".c")
File /usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py:232, in _compile_objects(decl, ufl_objects, object_names, module_name, parameters, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
228 import ffcx.compiler
230 # JIT uses module_name as prefix, which is needed to make names of all struct/function
231 # unique across modules
--> 232 _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, parameters=parameters)
234 ffibuilder = cffi.FFI()
235 ffibuilder.set_source(module_name, code_body, include_dirs=[ffcx.codegeneration.get_include_path()],
236 extra_compile_args=cffi_extra_compile_args, libraries=cffi_libraries)
File /usr/local/lib/python3.10/dist-packages/ffcx/compiler.py:97, in compile_ufl_objects(ufl_objects, object_names, prefix, parameters, visualise)
95 # Stage 1: analysis
96 cpu_time = time()
---> 97 analysis = analyze_ufl_objects(ufl_objects, parameters)
98 _print_timing(1, time() - cpu_time)
100 # Stage 2: intermediate representation
File /usr/local/lib/python3.10/dist-packages/ffcx/analysis.py:85, in analyze_ufl_objects(ufl_objects, parameters)
82 else:
83 raise TypeError("UFL objects not recognised.")
---> 85 form_data = tuple(_analyze_form(form, parameters) for form in forms)
86 for data in form_data:
87 elements += [convert_element(e) for e in data.unique_sub_elements]
File /usr/local/lib/python3.10/dist-packages/ffcx/analysis.py:85, in <genexpr>(.0)
82 else:
83 raise TypeError("UFL objects not recognised.")
---> 85 form_data = tuple(_analyze_form(form, parameters) for form in forms)
86 for data in form_data:
87 elements += [convert_element(e) for e in data.unique_sub_elements]
File /usr/local/lib/python3.10/dist-packages/ffcx/analysis.py:168, in _analyze_form(form, parameters)
165 complex_mode = "_Complex" in parameters["scalar_type"]
167 # Compute form metadata
--> 168 form_data = ufl.algorithms.compute_form_data(
169 form,
170 do_apply_function_pullbacks=True,
171 do_apply_integral_scaling=True,
172 do_apply_geometry_lowering=True,
173 preserve_geometry_types=(ufl.classes.Jacobian,),
174 do_apply_restrictions=True,
175 do_append_everywhere_integrals=False, # do not add dx integrals to dx(i) in UFL
176 complex_mode=complex_mode)
178 # Determine unique quadrature degree, quadrature scheme and
179 # precision per each integral data
180 for id, integral_data in enumerate(form_data.integral_data):
181 # Iterate through groups of integral data. There is one integral
182 # data for all integrals with same domain, itype, subdomain_id
(...)
188
189 # Extract precision
File /usr/local/lib/python3.10/dist-packages/ufl/algorithms/compute_form_data.py:407, in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
403 # TODO: This is a very expensive check... Replace with something
404 # faster!
405 preprocessed_form = reconstruct_form_from_integral_data(self.integral_data)
--> 407 check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode) # Currently testing how fast this is
409 # TODO: This member is used by unit tests, change the tests to
410 # remove this!
411 self.preprocessed_form = preprocessed_form
File /usr/local/lib/python3.10/dist-packages/ufl/algorithms/check_arities.py:177, in check_form_arity(form, arguments, complex_mode)
175 def check_form_arity(form, arguments, complex_mode=False):
176 for itg in form.integrals():
--> 177 check_integrand_arity(itg.integrand(), arguments, complex_mode)
File /usr/local/lib/python3.10/dist-packages/ufl/algorithms/check_arities.py:159, in check_integrand_arity(expr, arguments, complex_mode)
156 arguments = tuple(sorted(set(arguments),
157 key=lambda x: (x.number(), x.part())))
158 rules = ArityChecker(arguments)
--> 159 arg_tuples = map_expr_dag(rules, expr, compress=False)
160 args = tuple(a[0] for a in arg_tuples)
161 if args != arguments:
File /usr/local/lib/python3.10/dist-packages/ufl/corealg/map_dag.py:36, in map_expr_dag(function, expression, compress, vcache, rcache)
17 def map_expr_dag(function, expression,
18 compress=True,
19 vcache=None,
20 rcache=None):
21 """Apply a function to each subexpression node in an expression DAG.
22
23 If *compress* is ``True`` (default) the output object from
(...)
34 Return the result of the final function call.
35 """
---> 36 result, = map_expr_dags(function, [expression], compress=compress,
37 vcache=vcache,
38 rcache=rcache)
39 return result
File /usr/local/lib/python3.10/dist-packages/ufl/corealg/map_dag.py:99, in map_expr_dags(function, expressions, compress, vcache, rcache)
97 r = handlers[v._ufl_typecode_](v)
98 else:
---> 99 r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
101 # Optionally check if r is in rcache, a memory optimization
102 # to be able to keep representation of result compact
103 if compress:
File /usr/local/lib/python3.10/dist-packages/ufl/algorithms/check_arities.py:48, in ArityChecker.sum(self, o, a, b)
46 def sum(self, o, a, b):
47 if a != b:
---> 48 raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
49 return a
ArityMismatch: Adding expressions with non-matching form arguments ('v_0',) vs ('conj(v_0)',).
Any idea what might be causing this? And why it seemingly started happening out of the blue?