I am trying to use Carreau-Yasuda (CY) model in the simulation. I have successfully applied generalized cross model but when it comes to CY model i am getting error.
Here is the error detail
Traceback (most recent call last):
File “/home/adnan/simulation/non_newton/ttst.py”, line 18, in
shear_rate = project(shear_rate_expr, Ww)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/projection.py”, line 132, in project
A, b = assemble_system(a, L, bcs=bcs,
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py”, line 422, in assemble_system
b_dolfin_form = _create_dolfin_form(b_form, form_compiler_parameters)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py”, line 60, in _create_dolfin_form
return Form(form,
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py”, line 43, in init
ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py”, line 50, in mpi_jit
return local_jit(*args, **kwargs)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py”, line 100, in ffc_jit
return ffc.jit(ufl_form, parameters=p)
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 130, in jit_build
module, signature = dijitso.jit(jitable=ufl_object,
File “/usr/lib/python3/dist-packages/dijitso/jit.py”, line 165, in jit
header, source, dependencies = generate(jitable, name, signature, params[“generator”])
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 65, in jit_generate
code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
File “/usr/lib/python3/dist-packages/ffc/compiler.py”, line 142, in compile_form
return compile_ufl_objects(forms, “form”, object_names,
File “/usr/lib/python3/dist-packages/ffc/compiler.py”, line 190, in compile_ufl_objects
ir = compute_ir(analysis, prefix, parameters, jit)
File “/usr/lib/python3/dist-packages/ffc/representation.py”, line 182, in compute_ir
irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
File “/usr/lib/python3/dist-packages/ffc/representation.py”, line 182, in
irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
File “/usr/lib/python3/dist-packages/ffc/representation.py”, line 455, in _compute_integral_ir
ir = r.compute_integral_ir(itg_data,
File “/usr/lib/python3/dist-packages/ffc/uflacs/uflacsrepresentation.py”, line 43, in compute_integral_ir
ir = initialize_integral_ir(“uflacs”, itg_data, form_data, form_id)
File “/usr/lib/python3/dist-packages/ffc/representationutils.py”, line 203, in initialize_integral_ir
“needs_oriented”: needs_oriented_jacobian(form_data),
File “/usr/lib/python3/dist-packages/ffc/representationutils.py”, line 156, in needs_oriented_jacobian
element = create_element(ufl_element)
File “/usr/lib/python3/dist-packages/ffc/fiatinterface.py”, line 100, in create_element
element = _create_fiat_element(ufl_element)
File “/usr/lib/python3/dist-packages/ffc/fiatinterface.py”, line 198, in _create_fiat_element
element = ElementClass(fiat_cell, degree)
File “/usr/lib/python3/dist-packages/FIAT/lagrange.py”, line 43, in init
poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
File “/usr/lib/python3/dist-packages/FIAT/polynomial_set.py”, line 168, in init
dv = expansion_set.tabulate_derivatives(degree, pts)
File “/usr/lib/python3/dist-packages/FIAT/expansions.py”, line 277, in tabulate_derivatives
data = _tabulate_dpts(self._tabulate, 2, n, order, numpy.array(pts))
File “/usr/lib/python3/dist-packages/FIAT/expansions.py”, line 74, in _tabulate_dpts
symbolic_tab = tabulator(n, X)
File “/usr/lib/python3/dist-packages/FIAT/expansions.py”, line 227, in _tabulate
ref_pts = [sum(self.A[i][j] * pts[j] for j in range(m2)) + self.b[i]
File “/usr/lib/python3/dist-packages/FIAT/expansions.py”, line 227, in
ref_pts = [sum(self.A[i][j] * pts[j] for j in range(m2)) + self.b[i]
File “/usr/lib/python3/dist-packages/sympy/core/decorators.py”, line 261, in _func
other = sympify(other, strict=True)
File “/usr/lib/python3/dist-packages/sympy/core/sympify.py”, line 365, in sympify
return convertersuperclass
File “/usr/lib/python3/dist-packages/sympy/core/numbers.py”, line 1087, in new
num = _convert_numpy_types(num)
File “/usr/lib/python3/dist-packages/sympy/core/sympify.py”, line 86, in _convert_numpy_types
return Float(a, precision=prec)
File “/usr/lib/python3/dist-packages/sympy/core/numbers.py”, line 1143, in new
mpf = mlib.from_str(num, precision, rnd)
File “/usr/lib/python3/dist-packages/mpmath/libmp/libmpf.py”, line 1331, in from_str
man, exp = str_to_man_exp(x, base=10)
File “/usr/lib/python3/dist-packages/mpmath/libmp/libmpf.py”, line 1294, in str_to_man_exp
float(x)
ValueError: could not convert string to float: ‘np.float64(-1.0)’
And here is sample code
from dolfin import *
# Create a mesh
mesh = UnitSquareMesh(10, 10)
# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2) # Velocity field (vector)
Ww = FunctionSpace(mesh, "CG", 1) # Scalar function space for viscosity
# Define a Function (not TrialFunction) for velocity
u = Function(V) # Solution variable
# Compute shear rate (magnitude of strain rate tensor)
D = sym(grad(u)) # Symmetric gradient (strain rate tensor)
shear_rate_expr = sqrt(inner(D, D) + 1e-12) # Frobenius norm (stabilized)
# Project shear rate onto function space
shear_rate = project(shear_rate_expr, Ww)
# Define Carreau-Yasuda model
mu_0 = Constant(0.16) # Zero-shear viscosity
mu_inf = Constant(0.0036) # Infinite-shear viscosity
k = Constant(0.5) # Consistency index
n = Constant(0.4) # Power-law index
# Define mu_pp using projected shear_rate
mu_pp_expr = mu_inf + (mu_0 - mu_inf) * (1 + (1 / (k * shear_rate**2))**((n - 1) / 2))
# Project viscosity onto function space
mu_pp = project(mu_pp_expr, Ww)
# Solve weak formulation (if needed)
mu = Function(Ww)
v = TestFunction(Ww)
a = inner(mu, v) * dx
L = inner(mu_pp, v) * dx
solve(a == L, mu)