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}
dx = Measure('dx', domain=domain, metadata=metadata)
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!
Thanks in advance!
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)
oldgradu = grad(oldu)
oldF = eye+oldgradu
oldCGr = oldF.T * oldF
oldIc = tr(oldCGr)
oldJ = variable(det(oldF))
oldgradd = grad(oldd)
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)
gradu = grad(u)
gradd = grad(d)
dq = TestFunction(MIXEL)
du, dbeta, dd = split(dq)
graddu = grad(du)
graddd = grad(dd)
### Neo-Hooke as in fenicsx example
# Deformation gradient
F = variable(eye + gradu)
F_test = variable(eye + graddu)
# 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
+ gc*ls/2 *(abs(gradd)**2 - abs(oldgradd)**2)
metadata = {"quadrature_degree": 4}
dx = Measure('dx', domain=domain, metadata=metadata)
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)
69 expression = replace_quadratureweight(expression)
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.