# Guidance on how to implement non-Newtonian viscosity in Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark)

I’m am trying to make it such that viscosity varies with shear rate according to Power Law in Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark). A similar question was discussed a while ago for Fenics, but I am not familiar enough with FEA / Fenics / Fenicsx to use that example to piece together an effective solution: https://fenicsproject.org/qa/11271/nonlinear-non-newtonian-navier-stokes/.

After banging my head against the problem for a while and getting nowhere, I am hoping to get some guidance on how to implement this change. Here is one example of some unsuccessful code I’ve used:

# Define strain-rate tensor
def epsilon(u):

# Defining shear-rate tensor
def gamma(u):
return pow(2*inner(epsilon(u),epsilon(u)),0.5)

# Define viscosity using power law
def mufunc(u):
return K*pow(gamma(u),(N-1))

# # Define stress tensor
# def sigma(u, p):
#     return 2*mufunc(u)*epsilon(u) - p*Identity(len(u))

class InletVelocity():
def __init__(self, t, dp, vel, inrad):
self.t = t
self.dp = dp
self.vel = vel

def __call__(self, x):
values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = vel
return values

# Define the variational problem for the first step
f = Constant(mesh, PETSc.ScalarType((0, 0)))
F1 = rho / k * dot(u - u_n, v) * dx
F1 += inner(dot(1.5 * u_n - 0.5 * u_n1, 0.5 * nabla_grad(u + u_n)), v) * dx
F1 += 0.5 * mufunc(u) * inner(grad(u + u_n), grad(v)) * dx - dot(p_, div(v)) * dx
F1 += dot(f, v) * dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
A1 = create_matrix(a1)
b1 = create_vector(L1)

Running this modification results in the following error message:

Traceback (most recent call last):
File "/home/kyle/Desktop/testenv/fenicstest/syringe_flow_sim_2d_tenth.py", line 297, in <module>
a1 = form(lhs(F1))
^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/formoperators.py", line 71, in lhs
return compute_form_lhs(form)
^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 363, in compute_form_lhs
return compute_form_with_arity(form, 2)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 331, in compute_form_with_arity
return map_integrands(_transform, form)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 25, in map_integrands
mapped_integrals = [map_integrands(function, itg, only_integral_type)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 33, in map_integrands
return itg.reconstruct(function(itg.integrand()))
^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 327, in _transform
e, provides = pe.visit(e)
^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 109, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 109, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 109, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 113, in visit
r = h(o)
^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 62, in expr
raise ValueError(f"Found Argument in {ufl_err_str(x)}, this is an invalid expression.")
ValueError: Found Argument in <Power id=140541038225104>, this is an invalid expression

Since you have non-Newtonian viscosity, your problem is no longer linear and cannot use lhs/rhs. You must define a non linear residual form and use a NonlinearProblem, see for instance A nonlinear Poisson equation — FEniCSx tutorial

2 Likes

Thank you bleyerj for the response. I attempted to piece together a solution using the variational form in Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial and the non-linear solver approach in Implementation — FEniCSx tutorial, but I’m running into issues (both code and error message below).

import os
import numpy as np
import matplotlib.pyplot as plt
import gmsh
import pyvista
import tqdm.autonotebook
from mpi4py import MPI
from petsc4py import PETSc
from PIL import Image, ImageDraw, ImageFont
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction,
TrialFunction, VectorElement, as_vector, div, dot, ds, dx, inner,

import dolfinx.log as log
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.plot import vtk_mesh
from dolfinx.fem import (Constant, Function, FunctionSpace,
assemble_scalar, dirichletbc, form, locate_dofs_topological,
set_bc, locate_dofs_geometrical, petsc)

from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, NonlinearProblem,
create_vector, create_matrix, set_bc)

from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import (XDMFFile, VTXWriter, distribute_entity_data, gmshio)
from dolfinx.mesh import create_mesh, meshtags_from_entities
from ufl import (FacetNormal, FiniteElement, MixedElement, Identity, Measure, TestFunction,
TestFunctions, TrialFunction, VectorElement, as_vector, div, dot, ds,

## Creating mesh for a syringe
comm = MPI.COMM_WORLD
model_rank = 0

# Creating gmsh mesh
gdim = 2
gmsh.initialize()

# Define mesh sizes
L = 1
H = 1
res = 0.05

gmsh.option.setNumber("Mesh.MeshSizeMax", 30)

# Order points
ptlst = [
(0,0),
(0,H),
(L,H),
(L,0),
]

# Defining gmsh points from precompiled point list
pcount = 1
lcount = 1
lnlst = []
for pt in range(len(ptlst)):
if pt == 0:
pcount += 1
else:
lnlst.append(lcount)
pcount += 1
lcount += 1

lnlst.append(lcount)
pcount += 1
lcount += 1

# Creating surface

# Synchronizing mesh parts
gmsh.model.occ.synchronize()

# Defining nodes where fluid will flow
surfaces = gmsh.model.getEntities(dim=gdim)
assert (len(surfaces) == 1)
fm = surfaces[0][1]
gmsh.model.setPhysicalName(surfaces[0][0], fm, "Fluid")

# Finding boundary, inflow, and outflow nodes
inlet_marker, outlet_marker, wall_marker, idk_marker = fm+1,fm+2,fm+3,fm+4
inflow, outflow, walls = [], [], []
boundaries = gmsh.model.getBoundary(surfaces,oriented=False)
for boundary in boundaries:
center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
if np.allclose(center_of_mass, [0, H/2, 0]):
inflow.append(boundary[1])
elif np.allclose(center_of_mass, [L, H/2, 0]):
outflow.append(boundary[1])
else:
walls.append(boundary[1])

gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")

# Generating the mesh
gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh_instance = gmsh.model.mesh.generate(gdim)

# Creating mesh model with gmsh mesh
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model,comm=comm,rank=comm.rank,gdim=gdim)
ft.name = "Facet Markers"
ct.name = "Cell Markers"

with XDMFFile(mesh.comm, "results/mesh.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)

# Finalize gmsh
gmsh.finalize()

# The following resources were used to build this script:
# https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html
# https://fenicsproject.org/qa/11271/nonlinear-non-newtonian-navier-stokes/
# https://fenicsproject.discourse.group/t/ipcs-solver-with-time-dependent-pressure-bc/10616/5
# https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code1.html

## Defining problem specific parameters
t = 0                       # Start time
T = 1                      # Final time
dt = 1/1000                 # Time step size
num_steps = int(T/dt)       # Number of steps

# Physical properties of water (for testing)
# Units are kg, mm, s, mPa, (mm/s),
RHO = 1     # kg/(mm**3)

# Keeping viscosity constant for initial testing, but would like to make it vary with shear rate according to Power Law
# Defining visocity for a shear-thinning fluid using the Power Law
K_VALUE = 1
N_VALUE = 0.5
K = Constant(mesh, PETSc.ScalarType(K_VALUE))  # Rheological constant
N = Constant(mesh, PETSc.ScalarType(N_VALUE))  # Power-law index

n = FacetNormal(mesh)
f = Constant(mesh, PETSc.ScalarType((0, 0)))
k = Constant(mesh, PETSc.ScalarType(dt))
rho = Constant(mesh, PETSc.ScalarType(RHO))     # Density
mu = Constant(mesh, PETSc.ScalarType(K_VALUE))  # Setting initial viscosity
vel = 0.1 # Inlet velocity is constant for testing

# Defining Taylor-Hood function space
v_cg2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

u = Function(V)
v = TestFunction(V)
u_ = Function(V)
u_.name = "u"
u_s = Function(V)
u_n = Function(V)
u_n1 = Function(V)
p = TrialFunction(Q)
q = TestFunction(Q)
p_ = Function(Q)
p_.name = "p"
phi = Function(Q)
fdim = mesh.topology.dim-1

# Define strain-rate tensor
def epsilon(u):

# Compute the magnitude of the strain rate tensor (shear rate)
def gamma(u):
return pow(2*inner(epsilon(u),epsilon(u)),0.5)

# Define viscosity using power law
def mufunc(u):
return K*pow(gamma(u),(N-1))

class InletVelocity():
def __init__(self, t, vel):
self.t = t
self.vel = vel

def __call__(self, x):
values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = vel
return values

# Define boundary conditions
# Walls
wall_dofs = locate_dofs_topological(V, fdim, ft.find(wall_marker))
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_walls = dirichletbc(u_nonslip, wall_dofs, V)

# Inlet
u_inlet = Function(V)
inlet_velocity = InletVelocity(t,vel)
u_inlet.interpolate(inlet_velocity)
inflow_dofs = locate_dofs_topological(V, fdim, ft.find(inlet_marker))
bcu_inflow = dirichletbc(u_inlet, inflow_dofs)
bcu = [bcu_inflow, bcu_walls]

# Outlet
p_outlet = PETSc.ScalarType(0)
outflow_dofs = locate_dofs_topological(Q, fdim, ft.find(outlet_marker))
bcp_outlet = dirichletbc(p_outlet, outflow_dofs, Q)
bcp = [bcp_outlet]

# Define the variational problem for the first step
f = Constant(mesh, PETSc.ScalarType((0, 0)))
F1 = rho / k * dot(u - u_n, v) * dx
F1 += inner(dot(1.5 * u_n - 0.5 * u_n1, 0.5 * nabla_grad(u + u_n)), v) * dx
F1 += 0.5 * mufunc(u_n) * inner(grad(u + u_n), grad(v)) * dx - dot(p_, div(v)) * dx
F1 += dot(f,v) * dx
problem = NonlinearProblem(F1,u,bcs=[bcu,bcp])

# Define variational problem for the second step
L2 = form(-rho / k * dot(div(u_s), q) * dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)

# Define variational problem for the third step
a3 = form(rho * dot(u, v) * dx)
L3 = form(rho * dot(u_s, v) * dx - k * dot(nabla_grad(phi), v) * dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)

# Solver for step 1
solver1 = NewtonSolver(comm, problem)
solver1.convergence_criterion = "incremental"
solver1.rtol = 1e-6
solver1.report = True
ksp = solver1.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.MINRES)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")

# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)

## Solving
from pathlib import Path
folder = Path("results")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, "dfg2D-3-u.bp", [u_], engine="BP4")
vtx_p = VTXWriter(mesh.comm, "dfg2D-3-p.bp", [p_], engine="BP4")
vtx_u.write(t)
vtx_p.write(t)
progress = tqdm.autonotebook.tqdm(desc="Solving PDE", total=num_steps)

for i in range(num_steps):
progress.update(1)
# Update current time step
t += dt
# Update inlet velocity
inlet_velocity.t = t
inlet_velocity.vel = vel
u_inlet.interpolate(inlet_velocity)

# Step 1: Tentative velocity step
log.set_log_level(log.LogLevel.INFO)
n, converged = solver1.solve(u_s.vector)
assert (converged)
print(f"Number of interations: {n:d}")

# Step 2: Pressure corrrection step
with b2.localForm() as loc:
loc.set(0)
assemble_vector(b2, L2)
apply_lifting(b2, [a2], [bcp])
set_bc(b2, bcp)
solver2.solve(b2, phi.vector)
phi.x.scatter_forward()

p_.vector.axpy(1, phi.vector)
p_.x.scatter_forward()

# Step 3: Velocity correction step
with b3.localForm() as loc:
loc.set(0)
assemble_vector(b3, L3)
solver3.solve(b3, u_.vector)
u_.x.scatter_forward()

# Write solutions to file
vtx_u.write(t)
vtx_p.write(t)

# Update variable with solution from this time step
with u_.vector.localForm() as loc_, u_n.vector.localForm() as loc_n, u_n1.vector.localForm() as loc_n1:
loc_n.copy(loc_n1)
loc_.copy(loc_n)

vtx_u.close()
vtx_p.close()
Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 20%] Meshing curve 2 (Line)
Info    : [ 40%] Meshing curve 3 (Line)
Info    : [ 60%] Meshing curve 4 (Line)
Info    : [ 80%] Meshing curve 5 (Line)
Info    : Done meshing 1D (Wall 0.000327134s, CPU 0.000131s)
Info    : Meshing 2D...
Info    : Meshing surface 6 (Plane, Frontal-Delaunay for Quads)
Info    : Done meshing 2D (Wall 0.0151339s, CPU 0.015126s)
Info    : 461 nodes 905 elements
x86_64-conda-linux-gnu-cc: error: unrecognized command-line option '-march'
Traceback (most recent call last):
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/_distutils/unixccompiler.py", line 185, in _compile
self.spawn(compiler_so + cc_args + [src, '-o', obj] + extra_postargs)
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/_distutils/ccompiler.py", line 1041, in spawn
spawn(cmd, dry_run=self.dry_run, **kwargs)
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/_distutils/spawn.py", line 70, in spawn
raise DistutilsExecError(
distutils.errors.DistutilsExecError: command '/home/kyle/anaconda3/envs/fenicsx-env/bin/x86_64-conda-linux-gnu-cc' failed with exit code 1

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/cffi/ffiplatform.py", line 48, in _build
dist.run_command('build_ext')
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/dist.py", line 963, in run_command
super().run_command(command)
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/_distutils/dist.py", line 988, in run_command
cmd_obj.run()
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/command/build_ext.py", line 88, in run
_build_ext.run(self)
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/_distutils/command/build_ext.py", line 345, in run
self.build_extensions()
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/_distutils/command/build_ext.py", line 467, in build_extensions
self._build_extensions_serial()
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/_distutils/command/build_ext.py", line 493, in _build_extensions_serial
self.build_extension(ext)
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/command/build_ext.py", line 249, in build_extension
_build_ext.build_extension(self, ext)
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/_distutils/command/build_ext.py", line 548, in build_extension
objects = self.compiler.compile(
^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/_distutils/ccompiler.py", line 600, in compile
self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts)
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/setuptools/_distutils/unixccompiler.py", line 187, in _compile
raise CompileError(msg)
distutils.errors.CompileError: command '/home/kyle/anaconda3/envs/fenicsx-env/bin/x86_64-conda-linux-gnu-cc' failed with exit code 1

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/home/kyle/Desktop/testenv/fenicstest/nonnew_error_demo.py", line 224, in <module>
problem = NonlinearProblem(F1,u,bcs=[bcu,bcp])
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/petsc.py", line 712, in __init__
self._L = _create_form(F, form_compiler_options=form_compiler_options,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 188, in form
return _create_form(form)
^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 183, in _create_form
return _form(form)
^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 141, in _form
ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/jit.py", line 56, in mpi_jit
return local_jit(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/jit.py", line 204, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 199, in compile_forms
raise e
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 190, in compile_forms
impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 284, in _compile_objects
ffibuilder.compile(tmpdir=cache_dir, verbose=True, debug=cffi_debug)
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/cffi/api.py", line 725, in compile
return recompile(self, module_name, source, tmpdir=tmpdir,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/cffi/recompiler.py", line 1564, in recompile
outputfilename = ffiplatform.compile('.', ext,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/cffi/ffiplatform.py", line 20, in compile
outputfilename = _build(tmpdir, ext, compiler_verbose, debug)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/kyle/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/cffi/ffiplatform.py", line 54, in _build
raise VerificationError('%s: %s' % (e.__class__.__name__, e))
cffi.VerificationError: CompileError: command '/home/kyle/anaconda3/envs/fenicsx-env/bin/x86_64-conda-linux-gnu-cc' failed with exit code 1

Basically just changed u from TrialFunction(V) to Function(V), updated F1 variational equation with mufunc(u_n), called on NonLinearProblem() to solve for F1, and updated the downstream solver scripts to reflect what was in Implementation — FEniCSx tutorial. Any suggestions on what can be done here to avoide the error message? Many thanks ahead of time for any and all insight.

I ran the code in a Docker container and received more clarification on the error:

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 20%] Meshing curve 2 (Line)
Info    : [ 40%] Meshing curve 3 (Line)
Info    : [ 60%] Meshing curve 4 (Line)
Info    : [ 80%] Meshing curve 5 (Line)
Info    : Done meshing 1D (Wall 0.000349457s, CPU 0.00032s)
Info    : Meshing 2D...
Info    : Meshing surface 6 (Plane, Frontal-Delaunay for Quads)
Info    : Done meshing 2D (Wall 0.0220554s, CPU 0.018219s)
Info    : 461 nodes 905 elements
Traceback (most recent call last):
File "/root/fenicstestslim/nonnew_error_demo.py", line 236, in <module>
A3 = assemble_matrix(a3)
File "/usr/lib/python3.10/functools.py", line 889, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py", line 350, in assemble_matrix
A = _cpp.fem.petsc.create_matrix(a._cpp_object)
RuntimeError: Cannot create sparsity pattern. Form is not a bilinear.

It seems that solving this variational setup with variable mu isn’t as simple as making u a function, as the velocity correction step calls on u.

Since you are trying to create a bilinear form here, you need to use a trial-function, not a function in a3.

1 Like