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
6

can you help in thin for non linear stokes equation?
Generalised Newtonian Fluids

So far, we have considered the Stokes equations with a constant viscosity parameter. Such fluids are called Newtonian. Many fluids, however, become “thinner” (= less viscous) when they are sheared. The viscosity then becomes a function of the unknown flow velocity, i.e. the Stokes problem becomes nonlinear and you will have to implement an iterative method such as Newton’s method.

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import gmsh
import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, io, la, mesh
from dolfinx.fem.petsc import LinearProblem
from ufl import Measure, TrialFunction, TrialFunctions, TestFunction, TestFunctions, SpatialCoordinate, div, dx, grad, inner
from pathlib import Path

def mesh_benchmark_2d(resolution=5):
    """Creates a mesh for the two-dimensional DFG benchmark for flow around a cyclinder and writes the mesh to the file "benchmark_2d.msh".

    The inlet is marked with facet tag 1.
    The top and bottom walls are marked with facet tag 2.
    The obstacle boundary is marked with facet tag 3.
    The outlet is marked with facet tag 4.

    Args:
        resolution: The target grid spacing will be the obstacle radius divided by this number.
    """

    gmsh.initialize()

    channel_length = 2.2
    channel_height = .41
    obstacle_centre = [.2, .2]
    obstacle_radius = .05

    channel = gmsh.model.occ.addRectangle(
        x=0,
        y=0,
        z=0,
        dx=channel_length,
        dy=channel_height
    )
    obstacle = gmsh.model.occ.addDisk(
        xc=obstacle_centre[0],
        yc=obstacle_centre[1],
        zc=0,
        rx=obstacle_radius,
        ry=obstacle_radius
    )

    gdim = 2 # geometric dimension of this model

    gmsh.model.occ.cut([(gdim, channel)], [(gdim, obstacle)])
    gmsh.model.occ.synchronize()

    # Domain
    domain_entities = [entity[1] for entity in gmsh.model.getEntities(dim=gdim)]
    gmsh.model.addPhysicalGroup(gdim, domain_entities, tag=1)

    # Pieces of the boundary
    boundary_entities = [entity[1] for entity in gmsh.model.getEntities(dim=gdim-1)]

    inlet_entities = []
    wall_entities = []
    obstacle_entities = []
    outlet_entities = []

    for boundary_part in boundary_entities:
        xc = gmsh.model.occ.getCenterOfMass(
            dim=gdim-1,
            tag=boundary_part
        )

        if np.isclose(xc[0], 0.):
            # This piece of the boundary is part of the inlet
            inlet_entities.append(boundary_part)
        elif np.isclose(xc[0], channel_length):
            # This piece of the boundary is part of the outlet
            outlet_entities.append(boundary_part)
        elif np.isclose(xc[1], 0.) or np.isclose(xc[1], channel_height):
            # This piece of the boundary is part of the bottom or top wall
            wall_entities.append(boundary_part)
        else:
            # This piece of the boundary is part of the obstacle
            obstacle_entities.append(boundary_part)

    # We label the inlet with tag 1
    gmsh.model.addPhysicalGroup(gdim-1, inlet_entities, tag=1)
    gmsh.model.setPhysicalName(gdim-1, 1, "Inlet")

    # We label the top and bottom walls with tag 2
    gmsh.model.addPhysicalGroup(gdim-1, wall_entities, tag=2)
    gmsh.model.setPhysicalName(gdim-1, 2, "Walls")

    # We label the obstacle with tag 3
    gmsh.model.addPhysicalGroup(gdim-1, obstacle_entities, tag=3)
    gmsh.model.setPhysicalName(gdim-1, 3, "Obstacle")

    # We label the outlet with tag 4
    gmsh.model.addPhysicalGroup(gdim-1, outlet_entities, tag=4)
    gmsh.model.setPhysicalName(gdim-1, 4, "Outlet")

    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", obstacle_radius/resolution)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", obstacle_radius/resolution)
    
    gmsh.model.occ.synchronize()
    gmsh.model.mesh.generate(gdim)
    gmsh.write("benchmark_2d.msh")

    gmsh.clear()
    gmsh.finalize()


def dataset1(resolution=5):
    """Creates the DFG benchmark geometry and defines a set of parameters for a flow problem
    with parabolic inflow condition, no-slip conditions on the walls and a do-nothing condition
    at the outlet.

    Args:
        resolution: The target grid spacing will be the obstacle radius divided by this number.

    Returns:
        domain
        cell_tags
        facet_tags
        mu
        u_in
        u_wall
        f
        g
    """
    mesh_benchmark_2d(resolution)

    domain, cell_tags, facet_tags = io.gmshio.read_from_msh(
            filename="benchmark_2d.msh",
            comm=MPI.COMM_WORLD,
            gdim=2
    )

    # Viscosity
    mu = fem.Constant(domain, 1.)

    # Boundary conditions
    u_in = lambda x: (6. * x[1] * (.41 - x[1]) / (.41**2), 0. * x[0])
    u_wall = lambda x: (0.*x[0], 0.*x[0])

    # Right hand side
    f = fem.Constant(domain, ((0., 0.))) # body force density
    g = fem.Constant(domain, ((0., 0.))) # surface traction -µgrad(u).n + pn = g

    return domain, cell_tags, facet_tags, mu, u_in, u_wall, f, g


def dataset2(resolution=5):
    """Creates the DFG benchmark geometry and defines a set of parameters for the manufactured solution
    u = [1 + x² + 2y², -1 - 2xy]
    p = 3 + 4x + 5y

    Args:
        resolution: The target grid spacing will be the obstacle radius divided by this number.

    Returns:
        domain
        cell_tags
        facet_tags
        mu
        u_in
        u_wall
        f
        g
    """
    mesh_benchmark_2d(resolution)

    domain, cell_tags, facet_tags = io.gmshio.read_from_msh(
            filename="benchmark_2d.msh",
            comm=MPI.COMM_WORLD,
            gdim=2
    )

    # Viscosity
    mu = fem.Constant(domain, 1.)

    # Boundary conditions
    u_in = lambda x: (1. + x[0]**2. + 2. * x[1]**2., -1. - 2. * x[0] * x[1])
    u_wall = u_in

    # Right hand side
    f = fem.Constant(domain, ((-2., 5.))) # body force density
    x = SpatialCoordinate(domain)
    g = ufl.as_vector((7.4 + 5. * x[1], 2. * x[1])) # surface traction -µgrad(u).n + pn = g

    return domain, cell_tags, facet_tags, mu, u_in, u_wall, f, g


def exact_solution(domain):
    """ Returns the manufactured solution
    u = [1 + x² + 2y², -1 - 2xy]
    p = 3 + 4x + 5y
    as UFL expressions.
    """

    x = SpatialCoordinate(domain)
    ue = ufl.as_vector((1. + x[0]**2. + 2. * x[1]**2., -1. - 2. * x[0] * x[1]))
    pe = 3. + 4. * x[0] + 5. * x[1]

    return ue, pe


def solve_monolithic_lu(domain, cell_tags, facet_tags, mu, inflow_velocity, wall_velocity, f, g):
    """Solves the Stokes problem by applying LU decompoition to the complete system of equations

    Args:
        domain
        cell_tags
        facet_tags
        mu
        inflow_velocity
        wall_velocity
        f
        g
    
    Returns:
        uh, ph: Numerical solution
    """

    # Add information about boundary sections to the ds measure
    ds = Measure("ds", domain=domain, subdomain_data=facet_tags)

    # Finite elements and function spaces for velocity and pressure
    V_el = element("P", domain.topology.cell_name(), 2, shape=(2, ))
    Q_el = element("P", domain.topology.cell_name(), 1)
    VQ_el = mixed_element([V_el, Q_el])
    VQ = fem.functionspace(domain, VQ_el)

    (u, p) = TrialFunctions(VQ)
    (v, q) = TestFunctions(VQ)

    V, _ = VQ.sub(0).collapse()

    u_in = fem.Function(V)
    u_in.interpolate(inflow_velocity)

    u_wall = fem.Function(V)
    u_wall.interpolate(wall_velocity)

    inlet_dofs = fem.locate_dofs_topological((VQ.sub(0), V), 1, facet_tags.find(1))
    bc_in = fem.dirichletbc(u_in, inlet_dofs, VQ.sub(0))

    wall_dofs = fem.locate_dofs_topological((VQ.sub(0), V), 1, np.concat([facet_tags.find(2), facet_tags.find(3)]))
    bc_wall = fem.dirichletbc(u_wall, wall_dofs, VQ.sub(0))
    bcs = [bc_in, bc_wall]

    # Weak formulation of the Stokes problem
    a = mu * inner(grad(u), grad(v)) * dx - p * div(v) * dx - q * div(u) * dx
    L = inner(f, v) * dx - inner(g, v) * ds(4)

    # "Brute-force" solver: LU decomposition of the full matrix
    lu_solver = {
        "ksp_type": "preonly",
        "pc_type": "lu"
    }

    problem = LinearProblem(a, L, bcs=bcs, petsc_options=lu_solver)
    uph = problem.solve()

    # Split solution into velocity and pressure
    uh, ph = uph.sub(0).collapse(), uph.sub(1).collapse()
    uh.name, ph.name = "Velocity", "Pressure"

    return uh, ph


def solve_iterative_schur(domain, cell_tags, facet_tags, mu, inflow_velocity, wall_velocity, f, g):
    """Solves the Stokes problem by a Krylov subspace method and Schur-complement based block preconditioning.

    Implementation based in the FEniCSx demo https://docs.fenicsproject.org/dolfinx/main/python/demos/demo_stokes.html
    
    Args:
        domain
        cell_tags
        facet_tags
        mu
        inflow_velocity
        wall_velocity
        f
        g
    
    Returns:
        uh, ph: Numerical solution
    """

    # Add information about boundary sections to the ds measure
    ds = Measure("ds", domain=domain, subdomain_data=facet_tags)

    # Finite elements and function spaces for velocity and pressure
    V_el = element("P", domain.topology.cell_name(), 2, shape=(2, ))
    Q_el = element("P", domain.topology.cell_name(), 1)

    V = fem.functionspace(domain, ("P", 2, (2, )))
    Q = fem.functionspace(domain, ("P", 1))

    u, p = TrialFunction(V), TrialFunction(Q)
    v, q = TestFunction(V), TestFunction(Q)

    # Boundary conditions
    u_in = fem.Function(V)
    u_in.interpolate(inflow_velocity)

    u_wall = fem.Function(V)
    u_wall.interpolate(wall_velocity)

    inlet_dofs = fem.locate_dofs_topological(V, 1, facet_tags.find(1))
    bc_in = fem.dirichletbc(u_in, inlet_dofs)

    wall_dofs = fem.locate_dofs_topological(V, 1, np.concat([facet_tags.find(2), facet_tags.find(3)]))
    bc_wall = fem.dirichletbc(u_wall, wall_dofs)
    bcs = [bc_in, bc_wall]
    
    # Weak formulation of the Stokes problem
    a = fem.form([[mu * inner(grad(u), grad(v)) * dx, -p * div(v) * dx], [- q * div(u) * dx, None]])
    L = fem.form([inner(f, v) * dx - inner(g, v) * ds(4), fem.Constant(domain, 0.) * q * dx])
    
    # Block-diagonal preconditioner
    a_p11 = fem.form(p * q * dx)
    a_p = [[a[0][0], None], a_p11]

    # Assemble nested matrix operators
    A = fem.petsc.assemble_matrix_nest(a, bcs=bcs)
    A.assemble()

    # Create a nested matrix P to use as the preconditioner. The
    # top-left block of P is shared with the top-left block of A. The
    # bottom-right diagonal entry is assembled from the form a_p11:
    P11 = fem.petsc.assemble_matrix(a_p11, [])
    P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]])
    P.assemble()

    A00 = A.getNestSubMatrix(0, 0)
    A00.setOption(PETSc.Mat.Option.SPD, True)

    P00, P11 = P.getNestSubMatrix(0, 0), P.getNestSubMatrix(1, 1)
    P00.setOption(PETSc.Mat.Option.SPD, True)
    P11.setOption(PETSc.Mat.Option.SPD, True)

    # Assemble right-hand side vector
    b = fem.petsc.assemble_vector_nest(L)

    # Modify ('lift') the RHS for Dirichlet boundary conditions
    fem.petsc.apply_lifting_nest(b, a, bcs=bcs)

    # Sum contributions for vector entries that are shared across
    # parallel processes
    for b_sub in b.getNestSubVecs():
        b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

    # Set Dirichlet boundary condition values in the RHS vector
    bcs0 = fem.bcs_by_block(fem.extract_function_spaces(L), bcs)
    fem.petsc.set_bc_nest(b, bcs0)

    # Create a MINRES Krylov solver and a block-diagonal preconditioner
    # using PETSc's additive fieldsplit preconditioner
    ksp = PETSc.KSP().create(domain.comm)
    ksp.setOperators(A, P)
    ksp.setType("minres")
    ksp.setTolerances(rtol=1e-9)
    ksp.getPC().setType("fieldsplit")
    ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)

    # Define the matrix blocks in the preconditioner with the velocity
    # and pressure matrix index sets
    nested_IS = P.getNestISs()
    ksp.getPC().setFieldSplitIS(("u", nested_IS[0][0]), ("p", nested_IS[0][1]))

    # Set the preconditioners for each block. For the top-left
    # Laplace-type operator we use one algebraic multigrid iteration. For the
    # lower-right block we use a Jacobi (diagonal) preconditioner.
    ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()
    
    ksp_u.setType("preonly")
    ksp_u.getPC().setType("hypre")
    ksp_u.getPC().setHYPREType("boomeramg")

    ksp_p.setType("preonly")
    ksp_p.getPC().setType("jacobi")

    # Create finite element Functions for the velocity (on the space V)
    # and for the pressure (on the space Q). The vectors for uh and ph are
    # combined to form a nested vector and the system is solved.
    uh, ph = fem.Function(V), fem.Function(Q)
    uh.name, ph.name = "Velocity", "Pressure"
    x = PETSc.Vec().createNest([la.create_petsc_vector_wrap(uh.x), la.create_petsc_vector_wrap(ph.x)])
    ksp.solve(b, x)

    # Update ghost values
    uh.x.scatter_forward()
    ph.x.scatter_forward()

    return uh, ph


def errornorm(uh, ue, norm="L2"):
    """Computes the L² norm or the H¹ norm of the error between the numerical and the exact solution

    Args:
        uh: Numerical solution
        ue: Exact solution
        norm: "L2" (default) or "H1"

    Returns:
        error: Value of the error
    """
    if norm=="H1":
        errorform = fem.form(inner(uh - ue, uh - ue) * dx + inner(grad(uh - ue), grad(uh - ue)) * dx)
    else:
        errorform = fem.form(inner(uh - ue, uh - ue) * dx)
    
    error2_subdomain = fem.assemble_scalar(errorform)
    error2_total = MPI.COMM_WORLD.allreduce(error2_subdomain, op=MPI.SUM)

    return np.sqrt(error2_total)


def save_solution(uh):
    """Exports the numerical solution in VTX format

    Args:
        uh: Numerical solution
    """

    domain = uh.function_space.mesh

    results_folder = Path("results")
    results_folder.mkdir(exist_ok=True, parents=True)

    # We replace any whitespace characters in the name of uh with an underscore
    filename = "_".join(uh.name.split())

    with io.VTXWriter(domain.comm, results_folder / ("stokes_" + filename + ".bp"), [uh]) as vtx:
        vtx.write(t=0.)


if __name__ == "__main__":
    parameters = dataset1(5)

    uh, ph = solve_monolithic_lu(*parameters)
    uh, ph = solve_iterative_schur(*parameters)

    save_solution(uh)
    save_solution(ph)
7

Please format your code appropriately, and please try to minimize the amount of code. Have you looked at Implementation — FEniCSx tutorial for inspiration?

8

yes but I found it bit harder for the stokes problem.
I actually have to change the viscosity parameter but not understanding how to implement it in the code.

9 
 # Viscosity
    mu = fem.Constant(domain, 1.)

    # Boundary conditions
    u_in = lambda x: (1. + x[0]**2. + 2. * x[1]**2., -1. - 2. * x[0] * x[1])
    u_wall = u_in

    # Right hand side
    f = fem.Constant(domain, ((-2., 5.))) # body force density
    x = SpatialCoordinate(domain)
    g = ufl.as_vector((7.4 + 5. * x[1], 2. * x[1])) # surface traction -µgrad(u).n + pn = g

    return domain, cell_tags, facet_tags, mu, u_in, u_wall, f, g

hoe to make change in viscosity here?? according to statement" The viscosity then becomes a function of the unknown flow velocity, i.e. the Stokes problem becomes nonlinear and you will have to implement an iterative method such as Newton’s method.

10

please reply ??? How to change viscosity here in form of coding

11

It takes some time to get your question answered by the development team members. So, be patient.

As mu is a function, so you should not define it as a fem.Constant but

mu = fem.Funtion(Q)
12

What explicit relation is mu(u)?

13

just I have to change this in the code ore other parameters too?

14

can please someone help me to correct this code according to statement

15

Your question is fairly vague. You have not stated what relationship there is between u and my, ie what is the explicit formula that ties them together.

16

mu = fem.Constant(domain, 1.)
Mu is constant in the code and now we have to make it function of unknown parameter for generalised Newtonian fluids.
I have mentioned above the statement of my problem.its for 2D non linear navier stokes.
We have to use carreau fluid here,but exactly how I don’t know.

17

Define the Carreau viscosity model

def carreau_viscosity(u):
    mu_0=1.0
    mu_inf=0.1 
    lam=1.0
    n=0.5
    D = sym(grad(u))  # Rate-of-strain tensor
    shear_rate = sqrt(2 * inner(D, D))  # Magnitude of shear rate
    return mu_inf + (mu_0 - mu_inf) * (1 + (lam * shear_rate)**2)**((n - 1) / 2)
mu_expr = carreau_viscosity(u)

# Use Carreau viscosity in the weak formulation
mu = carreau_viscosity(u)
a = inner(2 *mu_expr * sym(grad(u)), sym(grad(v))) * dx - div(v) * p * dx - q * div(u) * dx
L = inner(f, v) * dx - inner(g, v) * ds(4)


div_u = div(u)
div_p = div(p)


# Solver configuration using LU decomposition
lu_solver = {
    "ksp_type": "preonly",
    "pc_type": "lu"
}

# Solve the problem
problem = NonlinearProblem(F, u, bcs=bcs)
uph = problem.solve()

# Extract and name velocity and pressure solutions
uh, ph = uph.sub(0).collapse(), uph.sub(1).collapse()
uh.name = "Velocity"
ph.name = "Pressure"

return uh, ph

type or paste code here
```Traceback (most recent call last):
  File "/home/pool/smuzsheh/Desktop/MAT877-Project/stokes.py", line 477, in <module>
    uh, ph = solve_monolithic_lu(*parameters)
             ~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^
  File "/home/pool/smuzsheh/Desktop/MAT877-Project/stokes.py", line 295, in solve_monolithic_lu
    problem = NonlinearProblem(F, u, bcs=bcs)
                               ^
NameError: name 'F' is not defined. Did you mean: 'f'? I am getting this error for define F how to define F in the code``
18

please help me I am using newtonian fluid model carreau model

19

Being pushy will not increase your chances of getting an answer. People do their best to provide help when they can…
Regarding your problem, the error message is explicit. F is not defined. The redidual F should be nonlinear in terms of the unknown Function. Nonlinear problems are covered here Implementation — FEniCSx tutorial

2 Likes