# How do I solve "Expecting scalar." error when defining NonlinearProblem?

Hello everyone,

I am currently trying to implement a phase-field model for fracture using fenicsx. Currently I am stuck trying to define a dolfinx.fem.petsc.NonlinearProblem class for my problem. I’ll write a short summary of my problem here and add an MWE and the full error at the end of this post. I am defining the NonlinearProblem like so:

``````metadata = {"quadrature_degree": 4}

P = diff(pi, c_vec)
PId = (dot(P,c_vec_test)-dot(gamma_vec, u_vec_test))*dx

print('shape of integrand =',shape(dot(P,c_vec_test)-dot(gamma_vec, u_vec_test)))

problem = NonlinearProblem(PId,q,bcs=bcs)
``````

Refering to the NonlinearProblem “RuntimeError: Expecting scalar.” appears here consistently. I am not even quite sure, which input this refers to. The residual PId would be the obvious choice, which is why I implemented the print(…)-line to check whether, PId has the right dimensions. This does confirm at least the integrand to be scalar though, with the output:

``````shape of integrand = ()
``````

Am I missing something about my PId? Or does the error refer to something else?

Here is the MWE with the full error. Sadly not quite so short, but I did try to keep it compact!

Max

``````### imports

import numpy as np

from mpi4py import MPI
from dolfinx.fem import Constant, FunctionSpace, Function
from dolfinx import default_scalar_type
from dolfinx.mesh import create_box, CellType
from ufl import (Identity, TestFunction, TrialFunction, VectorElement, FiniteElement, MixedElement, grad, tr,variable, diff,
Measure, inner, shape, sqrt, ln, det, derivative, split, as_vector,dot)
from dolfinx.fem.petsc import NonlinearProblem

### definition of domain

voldim = 3

domain = create_box(MPI.COMM_WORLD,
[[0,0,0],[5,1,1]], [50,10,10],
CellType.tetrahedron,np.float64) #, ghost_mode.none)

### boundary conditions not needd to recreate error
bcs = []

### Constant parameters

lam = Constant(domain, np.float64(100))
# 2nd lamé
G = Constant(domain, np.float64(50))
# bulk modulus
blk = Constant(domain, np.float64(lam+2*G/3))
# length scale of diffusive crack modeling
ls = Constant(domain, np.float64(0.25))
# viscosity parameter to dampen evolution of phase-field
eta = Constant(domain, np.float64(0.1))
# time step, this is set in relation to viscority ety
tau = Constant(domain, np.float64(1))
# threshold for range with/without diffusive crack accumulation
gc = Constant(domain, np.float64(0.1))
# constant k for numerical stability
k = Constant(domain, np.float64(1E-5))
# Identity Matrix
eye = Identity(voldim)

### Elements and function spaces

u_VP2 = VectorElement("Lagrange", domain.ufl_cell(), 2)
beta_SP1 = FiniteElement("Lagrange", domain.ufl_cell(), 1)
dam_SP1 = FiniteElement("Lagrange", domain.ufl_cell(), 1)

mixEl = MixedElement(u_VP2, beta_SP1, dam_SP1)
MIXEL = FunctionSpace(domain, mixEl)

# old_ are values from previous time step in full code

oldq = Function(MIXEL)
oldu, oldbeta, oldd = split(oldq)

oldCGr = oldF.T * oldF
oldIc = tr(oldCGr)
oldJ = variable(det(oldF))

oldPsi_0 = (G / 2) * (oldIc - voldim) - G * ln(oldJ) + (lam / 2) * (ln(oldJ))**2
oldg = (1-oldd)**2
oldPsi = (oldg+k)*oldPsi_0

### values for current time step

q = Function(MIXEL)
u, beta, d = split(q)

dq = TestFunction(MIXEL)
du, dbeta, dd = split(dq)

### Neo-Hooke as in fenicsx example
# Right Cauchy-Green
CGr = variable(F.T * F)
# Invariants of deformation tensor
Ic = variable(tr(CGr))
J = variable(det(F))
# Deformation energy
Psi_0 = (G / 2) * (Ic - voldim) - G * ln(J) + (lam / 2) * (ln(J))**2

gamma = variable(1/(2*ls) * d**2 + ls/2 * abs(gradd)**2)
gamma_vec = variable(as_vector(np.hstack((gamma,[0]*(voldim+1)))))

u_vec = variable(as_vector(np.hstack((u,d,beta))))
u_vec_test = variable(as_vector(np.hstack((du, dd, dbeta))))

F_vec = [F[0,0],F[1,1],F[2,2],F[1,2]+F[2,1],F[0,2]+F[2,0],F[1,0]+F[0,1]]
F_vec_test = [F_test[0,0],F_test[1,1],F_test[2,2],F_test[1,2]+F_test[2,1],F_test[0,2]+F_test[2,0],F_test[1,0]+F_test[0,1]]
c_vec = variable(as_vector(np.hstack((F_vec, d, gradd, beta))))
c_vec_test = variable(as_vector(np.hstack((F_vec_test, dd, graddd, dbeta))))

f = beta-gc/ls*oldd
pi = Psi_0 - oldPsi_0 + beta*(d-oldd)
- tau/(2*eta)*f*2

P = diff(pi, c_vec)
PId = (dot(P,c_vec_test)-dot(gamma_vec, u_vec_test))*dx
#- inner(tN,u_vec_test)*ds(2) #This is commented, because the problem will be driven by displacement driven, not traction driven
print('shape of integrand =',shape(dot(P,c_vec_test)-dot(gamma_vec, u_vec_test)))

problem = NonlinearProblem(PId,q,bcs=bcs)

``````
``````shape of integrand = ()

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[19], line 106
103 #- inner(tN,u_vec_test)*ds(2) #This is commented, because the problem will be driven by displacement driven, not traction driven
104 print('shape of integrand =',shape(dot(P,c_vec_test)-dot(gamma_vec, u_vec_test)))
--> 106 problem = NonlinearProblem(PId,q,bcs=bcs)

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/dolfinx/fem/petsc.py:712, in NonlinearProblem.__init__(self, F, u, bcs, J, form_compiler_options, jit_options)
689 def __init__(self, F: ufl.form.Form, u: _Function, bcs: typing.List[DirichletBC] = [],
690              J: ufl.form.Form = None, form_compiler_options: typing.Optional[dict] = None,
691              jit_options: typing.Optional[dict] = None):
692     """Initialize solver for solving a non-linear problem using Newton's method, :math:`(dF/du)(u) du = -F(u)`.
693
694     Args:
(...)
710
711     """
--> 712     self._L = _create_form(F, form_compiler_options=form_compiler_options,
713                            jit_options=jit_options)
715     # Create the Jacobian matrix, dF/du
716     if J is None:

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/dolfinx/fem/forms.py:188, in form(form, dtype, form_compiler_options, jit_options)
185         return list(map(lambda sub_form: _create_form(sub_form), form))
186     return form
--> 188 return _create_form(form)

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/dolfinx/fem/forms.py:183, in form.<locals>._create_form(form)
180 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
181 return form argument"""
182 if isinstance(form, ufl.Form):
--> 183     return _form(form)
184 elif isinstance(form, collections.abc.Iterable):
185     return list(map(lambda sub_form: _create_form(sub_form), form))

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/dolfinx/fem/forms.py:141, in form.<locals>._form(form)
139 if mesh is None:
140     raise RuntimeError("Expecting to find a Mesh in the form.")
--> 141 ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
142                                        form_compiler_options=form_compiler_options,
143                                        jit_options=jit_options)
145 # For each argument in form extract its function space
146 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]

File ~/miniforge3/envs/py23/lib/python3.11/site-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 ~/miniforge3/envs/py23/lib/python3.11/site-packages/dolfinx/jit.py:204, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
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], options=p_ffcx, **p_jit)
205 elif isinstance(ufl_object, ufl.FiniteElementBase):
206     r = ffcx.codegeneration.jit.compile_elements([ufl_object], options=p_ffcx, **p_jit)

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/ffcx/codegeneration/jit.py:199, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
197     except Exception:
198         pass
--> 199     raise e
201 obj, module = _load_objects(cache_dir, module_name, form_names)
202 return obj, module, (decl, impl)

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/ffcx/codegeneration/jit.py:190, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
187     for name in form_names:
188         decl += form_template.format(name=name)
--> 190     impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
191                             cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
192 except Exception as e:
193     try:
194         # remove c file so that it will not timeout next time

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/ffcx/codegeneration/jit.py:260, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
256 import ffcx.compiler
258 # JIT uses module_name as prefix, which is needed to make names of all struct/function
259 # unique across modules
--> 260 _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
262 ffibuilder = cffi.FFI()
263 ffibuilder.set_source(module_name, code_body, include_dirs=[ffcx.codegeneration.get_include_path()],
264                       extra_compile_args=cffi_extra_compile_args, libraries=cffi_libraries)

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/ffcx/compiler.py:102, in compile_ufl_objects(ufl_objects, object_names, prefix, options, visualise)
100 # Stage 2: intermediate representation
101 cpu_time = time()
--> 102 ir = compute_ir(analysis, object_names, prefix, options, visualise)
103 _print_timing(2, time() - cpu_time)
105 # Stage 3: code generation

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/ffcx/ir/representation.py:197, in compute_ir(analysis, object_names, prefix, options, visualise)
191 ir_elements = [_compute_element_ir(e, analysis.element_numbers, finite_element_names)
192                for e in analysis.unique_elements]
194 ir_dofmaps = [_compute_dofmap_ir(e, analysis.element_numbers, dofmap_names)
195               for e in analysis.unique_elements]
--> 197 irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
198                             options, visualise)
199        for (i, fd) in enumerate(analysis.form_data)]
200 ir_integrals = list(itertools.chain(*irs))
202 ir_forms = [_compute_form_ir(fd, i, prefix, form_names, integral_names, analysis.element_numbers,
203                              finite_element_names, dofmap_names, object_names)
204             for (i, fd) in enumerate(analysis.form_data)]

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/ffcx/ir/representation.py:197, in <listcomp>(.0)
191 ir_elements = [_compute_element_ir(e, analysis.element_numbers, finite_element_names)
192                for e in analysis.unique_elements]
194 ir_dofmaps = [_compute_dofmap_ir(e, analysis.element_numbers, dofmap_names)
195               for e in analysis.unique_elements]
--> 197 irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
198                             options, visualise)
199        for (i, fd) in enumerate(analysis.form_data)]
200 ir_integrals = list(itertools.chain(*irs))
202 ir_forms = [_compute_form_ir(fd, i, prefix, form_names, integral_names, analysis.element_numbers,
203                              finite_element_names, dofmap_names, object_names)
204             for (i, fd) in enumerate(analysis.form_data)]

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/ffcx/ir/representation.py:492, in _compute_integral_ir(form_data, form_index, element_numbers, integral_names, finite_element_names, options, visualise)
489 integrands = {rule: integral.integrand() for rule, integral in sorted_integrals.items()}
491 # Build more specific intermediate representation
--> 492 integral_ir = compute_integral_ir(itg_data.domain.ufl_cell(), itg_data.integral_type,
493                                   ir["entitytype"], integrands, ir["tensor_shape"],
494                                   options, visualise)
496 ir.update(integral_ir)
498 # Fetch name

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/ffcx/ir/integral.py:72, in compute_integral_ir(cell, integral_type, entitytype, integrands, argument_shape, p, visualise)
71 # Build initial scalar list-based graph representation
---> 72 S = build_scalar_graph(expression)
74 # Build terminal_data from V here before factorization. Then we
75 # can use it to derive table properties for all modified
76 # terminals, and then use that to rebuild the scalar graph more
77 # efficiently before argument factorization. We can build
78 # terminal_data again after factorization if that's necessary.
80 initial_terminals = {i: analyse_modified_terminal(v['expression'])
81                      for i, v in S.nodes.items()
82                      if is_modified_terminal(v['expression'])}

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/ffcx/ir/analysis/graph.py:80, in build_scalar_graph(expression)
77 G = build_graph_vertices([expression], skip_terminal_modifiers=False)
79 # Build more fine grained computational graph of scalar subexpressions
---> 80 scalar_expressions = rebuild_with_scalar_subexpressions(G)
82 # Build new list representation of graph where all
83 # vertices of V represent single scalar operations
84 G = build_graph_vertices(scalar_expressions, skip_terminal_modifiers=True)

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/ffcx/ir/analysis/graph.py:180, in rebuild_with_scalar_subexpressions(G)
177 wops = [tuple(W[k] for k in so) for so in sops]
179 # Reconstruct scalar subexpressions of v
--> 180 ws = reconstruct(expr, wops)
182 # Store all scalar subexpressions for v symbols
183 if len(vs) != len(ws):

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/ffcx/ir/analysis/reconstruct.py:161, in reconstruct(o, *args)
159 f = _reconstruct_call_lookup.get(type(o), False)
160 if f:
--> 161     return f(o, *args)
162 else:
163     # Look for parent class types instead
164     for k in _reconstruct_call_lookup.keys():

File ~/miniforge3/envs/py23/lib/python3.11/site-packages/ffcx/ir/analysis/reconstruct.py:12, in handle_scalar_nary(o, ops)
10 def handle_scalar_nary(o, ops):
11     if o.ufl_shape != ():
---> 12         raise RuntimeError("Expecting scalar.")
13     sops = [op[0] for op in ops]
14     return [o._ufl_expr_reconstruct_(*sops)]

RuntimeError: Expecting scalar.
``````

I guess this is the offending lines.
Here you are mixing in numpy and ufl, a dangerous combinatation. Could you elaborate on what you are trying to do when defining `c_vec` and `c_vec_test`?