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

General
1

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):
    return grad(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
        self.inrad = inrad

    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
2

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
3

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, 
                 lhs, grad, nabla_grad, rhs, sym, sqrt, derivative, SpatialCoordinate)

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.graph import adjacencylist
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, 
                 dx, inner, lhs, grad, nabla_grad, rhs, sym, split)

## 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:
        gmsh.model.occ.addPoint(ptlst[pt][0],ptlst[pt][1],0,res,pcount)
        pcount += 1
    else:
        gmsh.model.occ.addPoint(ptlst[pt][0],ptlst[pt][1],0,res,pcount)
        gmsh.model.occ.addLine(pcount-1,pcount,lcount)
        lnlst.append(lcount)
        pcount += 1
        lcount += 1

# Adding the last line
gmsh.model.occ.addPoint(ptlst[pt][0],ptlst[pt][1],0,res,pcount)
gmsh.model.occ.addLine(pcount,1,lcount)
lnlst.append(lcount)
pcount += 1
lcount += 1

# Creating surface
gmsh.model.occ.addCurveLoop(lnlst,lnlst[-1]+1) # base surface
gmsh.model.occ.addPlaneSurface([lnlst[-1]+1],lnlst[-1]+2)
    
# 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.addPhysicalGroup(surfaces[0][0], [fm], fm)
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.addPhysicalGroup(1, walls, wall_marker)
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
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):
    return sym(grad(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
a2 = form(dot(grad(p), grad(q)) * dx)
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])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    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)
    b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    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.

4

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.

5

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

1 Like