If I take the Mixed Elements code from the Workshop and want to solve it with the Petsc Nonlinear problem Solver fem.petsc.NonlinearProblem ( I know it’s linear, just as an MWE), what is the right approach for doing that?
This MWE from that problem:
from mpi4py import MPI
import numpy as np
import basix.ufl
import dolfinx
import ufl
from petsc4py import PETSc
import warnings
import dolfinx.fem.petsc
M = 6
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, M, M)
el_u = basix.ufl.element(“Lagrange”, mesh.basix_cell(), 3, shape=(mesh.geometry.dim,))
el_p = basix.ufl.element(“Lagrange”, mesh.basix_cell(), 2)
el_mixed = basix.ufl.mixed_element([el_u, el_p])
W = dolfinx.fem.functionspace(mesh, el_mixed)
u, p = ufl.TrialFunctions(W)
v, q = ufl.TestFunctions(W)
wh = dolfinx.fem.Function(W)
uh, ph = ufl.split(wh)
ds = ufl.Measure(“ds”, domain=mesh)
g = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0, 0)))
f = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0, 0)))
F = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
F += ufl.inner(p, ufl.div(v)) * ufl.dx
F += ufl.inner(ufl.div(u), q) * ufl.dx
F -= ufl.inner(f, v) * ufl.dx
def boundary_marker(x):
return np.isclose(x[0], 0.0)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
left_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, boundary_marker)
def top_bottom_marker(x):
return np.isclose(x[1], 1.0) | np.isclose(x[1], 0.0)
tb_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, top_bottom_marker)
all_boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
remaining_facets = np.setdiff1d(all_boundary_facets, np.union1d(tb_facets, left_facets))
W0 = W.sub(0)
V, V_to_W0 = W0.collapse()
u_inlet_x = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(1))
u_inlet_y = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0))
dofs_inlet_x = dolfinx.fem.locate_dofs_topological(W.sub(0).sub(0), mesh.topology.dim - 1, left_facets)
dofs_inlet_y = dolfinx.fem.locate_dofs_topological(W.sub(0).sub(1), mesh.topology.dim - 1, left_facets)
bc_inlet_x = dolfinx.fem.dirichletbc(u_inlet_x, dofs_inlet_x, W.sub(0).sub(0))
bc_inlet_y = dolfinx.fem.dirichletbc(u_inlet_y, dofs_inlet_y, W.sub(0).sub(1))
u_wall = dolfinx.fem.Function(V)
u_wall.x.array[:] = 0
dofs_wall = dolfinx.fem.locate_dofs_topological((W0, V), mesh.topology.dim - 1, tb_facets)
bc_wall = dolfinx.fem.dirichletbc(u_wall, dofs_wall, W0)
bcs = [bc_inlet_x, bc_inlet_y, bc_wall]
Throws me an arity error when I try to naively squeeze it into the nonlinear problem
ptsc_opts_begin = {
“snes_max_it”:100,
“snes_linesearch_monitor”: None,
“snes_linesearch_type”: “None”,
“ksp_type”: “preonly”,
“pc_type”: “lu”,
“pc_factor_mat_solver_type”: “mumps”,
“mat_mumps_icntl_14”: 120,
“ksp_error_if_not_converged”: True,
“snes_atol”: 1e-12,
“snes_rtol”: 1e-10,
“snes_stol”: 1e-10,}
Resid = ufl.extract_blocks(F)
problem = dolfinx.fem.petsc.NonlinearProblem(Resid, [wh], bcs=bcs, petsc_options_prefix=“eulelastb”,
petsc_options=ptsc_opts_begin
)
Error Output ArityMismatch
---------------------------------------------------------------------------
ArityMismatch Traceback (most recent call last)
Cell In[4], line 2
1 Resid = ufl.extract_blocks(F)
----> 2 problem = dolfinx.fem.petsc.NonlinearProblem(Resid, [wh], bcs=bcs, petsc_options_prefix="eulelastb",
3 petsc_options=ptsc_opts_begin
4 )
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/petsc.py:1239, in NonlinearProblem.__init__(self, F, u, petsc_options_prefix, bcs, J, P, kind, petsc_options, form_compiler_options, jit_options, entity_maps)
1176 """
1177 Initialize solver for a nonlinear variational problem.
1178
(...) 1236 `EntityMap<dolfinx.mesh.EntityMap>` must be provided.
1237 """
1238 # Compile residual and Jacobian forms
-> 1239 self._F = _create_form(
1240 F,
1241 form_compiler_options=form_compiler_options,
1242 jit_options=jit_options,
1243 entity_maps=entity_maps,
1244 )
1246 if J is None:
1247 J = derivative_block(F, u)
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:449, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
446 else:
447 return form
--> 449 return _create_form(form)
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:445, in form.<locals>._create_form(form)
443 return _zero_form(form)
444 elif isinstance(form, collections.abc.Iterable):
--> 445 return list(map(lambda sub_form: _create_form(sub_form), form))
446 else:
447 return form
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:445, in form.<locals>._create_form.<locals>.<lambda>(sub_form)
443 return _zero_form(form)
444 elif isinstance(form, collections.abc.Iterable):
--> 445 return list(map(lambda sub_form: _create_form(sub_form), form))
446 else:
447 return form
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:445, in form.<locals>._create_form(form)
443 return _zero_form(form)
444 elif isinstance(form, collections.abc.Iterable):
--> 445 return list(map(lambda sub_form: _create_form(sub_form), form))
446 else:
447 return form
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:445, in form.<locals>._create_form.<locals>.<lambda>(sub_form)
443 return _zero_form(form)
444 elif isinstance(form, collections.abc.Iterable):
--> 445 return list(map(lambda sub_form: _create_form(sub_form), form))
446 else:
447 return form
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:441, in form.<locals>._create_form(form)
431 """Recursively convert ufl.Forms to dolfinx.fem.Form.
432
433 Args:
(...) 438 A ``dolfinx.fem.Form`` or a list of ``dolfinx.fem.Form``.
439 """
440 if isinstance(form, ufl.Form):
--> 441 return _form(form)
442 elif isinstance(form, ufl.ZeroBaseForm):
443 return _zero_form(form)
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:361, in form.<locals>._form(form)
359 if msh is None:
360 raise RuntimeError("Expecting to find a Mesh in the form.")
--> 361 ufcx_form, module, code = jit.ffcx_jit(
362 msh.comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
363 )
365 # For each argument in form extract its function space
366 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py:60, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
56 @functools.wraps(local_jit)
57 def mpi_jit(comm, *args, **kwargs):
58 # Just call JIT compiler when running in serial
59 if comm.size == 1:
---> 60 return local_jit(*args, **kwargs)
62 # Default status (0 == ok, 1 == fail)
63 status = 0
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py:215, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
213 # Switch on type and compile, returning cffi object
214 if isinstance(ufl_object, ufl.Form):
--> 215 r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
216 elif isinstance(ufl_object, tuple) and isinstance(ufl_object[0], ufl.core.expr.Expr):
217 r = ffcx.codegeneration.jit.compile_expressions([ufl_object], options=p_ffcx, **p_jit)
File /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:224, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
221 for name in form_names:
222 decl += form_template.format(name=name)
--> 224 impl = _compile_objects(
225 decl,
226 forms,
227 form_names,
228 module_name,
229 p,
230 cache_dir,
231 cffi_extra_compile_args,
232 cffi_verbose,
233 cffi_debug,
234 cffi_libraries,
235 visualise=visualise,
236 )
237 except Exception as e:
238 try:
239 # remove c file so that it will not timeout next time
File /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:349, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
345 libraries = _libraries + cffi_libraries if cffi_libraries is not None else _libraries
347 # JIT uses module_name as prefix, which is needed to make names of all struct/function
348 # unique across modules
--> 349 _, code_body = ffcx.compiler.compile_ufl_objects(
350 ufl_objects, prefix=module_name, options=options, visualise=visualise
351 )
353 # Raise error immediately prior to compilation if no support for C99
354 # _Complex. Doing this here allows FFCx to be used for complex codegen on
355 # Windows.
356 if sys.platform.startswith("win32"):
File /dolfinx-env/lib/python3.12/site-packages/ffcx/compiler.py:108, in compile_ufl_objects(ufl_objects, options, object_names, prefix, visualise)
106 # Stage 1: analysis
107 cpu_time = time()
--> 108 analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"]) # type: ignore
109 _print_timing(1, time() - cpu_time)
111 # Stage 2: intermediate representation
File /dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py:94, in analyze_ufl_objects(ufl_objects, scalar_type)
91 else:
92 raise TypeError("UFL objects not recognised.")
---> 94 form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
95 for data in form_data:
96 elements += data.unique_sub_elements # type: ignore
File /dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py:94, in <genexpr>(.0)
91 else:
92 raise TypeError("UFL objects not recognised.")
---> 94 form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
95 for data in form_data:
96 elements += data.unique_sub_elements # type: ignore
File /dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py:180, in _analyze_form(form, scalar_type)
177 complex_mode = np.issubdtype(scalar_type, np.complexfloating)
179 # Compute form metadata
--> 180 form_data: ufl.algorithms.formdata.FormData = ufl.algorithms.compute_form_data(
181 form,
182 do_apply_function_pullbacks=True,
183 do_apply_integral_scaling=True,
184 do_apply_geometry_lowering=True,
185 preserve_geometry_types=(ufl.classes.Jacobian,),
186 do_apply_restrictions=True,
187 do_append_everywhere_integrals=False, # do not add dx integrals to dx(i) in UFL
188 complex_mode=complex_mode,
189 )
191 # Determine unique quadrature degree and quadrature scheme
192 # per each integral data
193 for id, integral_data in enumerate(form_data.integral_data): # type: ignore
194 # Iterate through groups of integral data. There is one integral
195 # data for all integrals with same domain, itype, subdomain_id
(...) 199 # all integrals in this integral data group, i.e. must be the
200 # same for for the same (domain, itype, subdomain_id)
File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/compute_form_data.py:482, 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, do_replace_functions, coefficients_to_split, complex_mode, do_remove_component_tensors)
479 preprocessed_form = reconstruct_form_from_integral_data(self.integral_data)
481 # TODO: Test how fast this is
--> 482 check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)
484 # TODO: This member is used by unit tests, change the tests to
485 # remove this!
486 self.preprocessed_form = preprocessed_form
File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py:213, in check_form_arity(form, arguments, complex_mode)
211 """Check the arity of a form."""
212 for itg in form.integrals():
--> 213 check_integrand_arity(itg.integrand(), arguments, complex_mode)
File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py:194, in check_integrand_arity(expr, arguments, complex_mode)
192 arguments = tuple(sorted(set(arguments), key=lambda x: (x.number(), x.part())))
193 rules = ArityChecker(arguments)
--> 194 arg_tuples = map_expr_dag(rules, expr, compress=False)
195 args = tuple(a[0] for a in arg_tuples)
196 if args != arguments:
File /dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py:37, in map_expr_dag(function, expression, compress, vcache, rcache)
15 def map_expr_dag(function, expression, compress=True, vcache=None, rcache=None):
16 """Apply a function to each subexpression node in an expression DAG.
17
18 If the same function is called multiple times in a transformation
(...) 35 The result of the final function call
36 """
---> 37 (result,) = map_expr_dags(
38 function, [expression], compress=compress, vcache=vcache, rcache=rcache
39 )
40 return result
File /dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py:109, in map_expr_dags(function, expressions, compress, vcache, rcache)
107 r = handlers[v._ufl_typecode_](v)
108 else:
--> 109 r = handlers[v._ufl_typecode_](v, *(vcache[u] for u in v.ufl_operands))
111 # Optionally check if r is in rcache, a memory optimization
112 # to be able to keep representation of result compact
113 if compress:
File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py:58, in ArityChecker.sum(self, o, a, b)
56 """Apply to sum."""
57 if a != b:
---> 58 raise ArityMismatch(
59 f"Adding expressions with non-matching form arguments "
60 f"{tuple(map(_afmt, a))} vs {tuple(map(_afmt, b))}."
61 )
62 return a
ArityMismatch: Adding expressions with non-matching form arguments ('v_0',) vs ('v_0', 'v_1').
even though it solves linearly with
a, L = ufl.system(F)
a_compiled = dolfinx.fem.form(a)
L_compiled = dolfinx.fem.form(L)
A = dolfinx.fem.create_matrix(a_compiled)
b = dolfinx.fem.create_vector(L_compiled.function_spaces[0])
A_scipy = A.to_scipy()
bcs = [bc_inlet_x, bc_inlet_y, bc_wall]
dolfinx.fem.assemble_matrix(A, a_compiled, bcs=bcs)
dolfinx.fem.assemble_vector(b.array, L_compiled)
dolfinx.fem.apply_lifting(b.array, [a_compiled], [bcs])
b.scatter_reverse(dolfinx.la.InsertMode.add)
[bc.set(b.array) for bc in bcs]
import scipy.sparse
A_inv = scipy.sparse.linalg.splu(A_scipy)
wh = dolfinx.fem.Function(W)
wh.x.array[:] = A_inv.solve(b.array)
Is it related to the use of TrialFunctions, which are for linear problems? If I try:
F = ufl.replace(F, {u: uh, p:ph})
Resid = ufl.extract_blocks(F)
jac = ufl.derivative(F, uh, v) + ufl.derivative(F, ph, q)
J = ufl.extract_blocks(jac)
problem = dolfinx.fem.petsc.NonlinearProblem(Resid, wh, bcs=bcs, J=J,petsc_options_prefix=“eulelastb”,
petsc_options=ptsc_opts_begin
)
I get:
Attribute Error: 'FunctionSpace' object has no attribute '\_cpp_object'
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[7], line 4
2 jac = ufl.derivative(F, uh, v) + ufl.derivative(F, ph, q)
3 J = ufl.extract_blocks(jac)
----> 4 problem = dolfinx.fem.petsc.NonlinearProblem(Resid, wh, bcs=bcs, J=J,petsc_options_prefix="eulelastb",
5 petsc_options=ptsc_opts_begin
6 )
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/petsc.py:1239, in NonlinearProblem.__init__(self, F, u, petsc_options_prefix, bcs, J, P, kind, petsc_options, form_compiler_options, jit_options, entity_maps)
1176 """
1177 Initialize solver for a nonlinear variational problem.
1178
(...) 1236 `EntityMap<dolfinx.mesh.EntityMap>` must be provided.
1237 """
1238 # Compile residual and Jacobian forms
-> 1239 self._F = _create_form(
1240 F,
1241 form_compiler_options=form_compiler_options,
1242 jit_options=jit_options,
1243 entity_maps=entity_maps,
1244 )
1246 if J is None:
1247 J = derivative_block(F, u)
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:449, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
446 else:
447 return form
--> 449 return _create_form(form)
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:445, in form.<locals>._create_form(form)
443 return _zero_form(form)
444 elif isinstance(form, collections.abc.Iterable):
--> 445 return list(map(lambda sub_form: _create_form(sub_form), form))
446 else:
447 return form
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:445, in form.<locals>._create_form.<locals>.<lambda>(sub_form)
443 return _zero_form(form)
444 elif isinstance(form, collections.abc.Iterable):
--> 445 return list(map(lambda sub_form: _create_form(sub_form), form))
446 else:
447 return form
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:445, in form.<locals>._create_form(form)
443 return _zero_form(form)
444 elif isinstance(form, collections.abc.Iterable):
--> 445 return list(map(lambda sub_form: _create_form(sub_form), form))
446 else:
447 return form
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:445, in form.<locals>._create_form.<locals>.<lambda>(sub_form)
443 return _zero_form(form)
444 elif isinstance(form, collections.abc.Iterable):
--> 445 return list(map(lambda sub_form: _create_form(sub_form), form))
446 else:
447 return form
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:441, in form.<locals>._create_form(form)
431 """Recursively convert ufl.Forms to dolfinx.fem.Form.
432
433 Args:
(...) 438 A ``dolfinx.fem.Form`` or a list of ``dolfinx.fem.Form``.
439 """
440 if isinstance(form, ufl.Form):
--> 441 return _form(form)
442 elif isinstance(form, ufl.ZeroBaseForm):
443 return _zero_form(form)
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:366, in form.<locals>._form(form)
361 ufcx_form, module, code = jit.ffcx_jit(
362 msh.comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
363 )
365 # For each argument in form extract its function space
--> 366 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]
367 part = form_compiler_options.get("part", "full")
368 if part == "diagonal":
AttributeError: 'FunctionSpace' object has no attribute '_cpp_object'
Should one avoid the use of Mixed Elements in that case, perhaps follow instead Coupling PDEs instead? In that case do we set bcs on each function space separately?
