Bilinear form for linearized steady NSE

General
1

Hi,

I’m trying to program a fixed point method for steady NSE. Unfortunately, I get the «classical»:

RuntimeError: Cannot create sparsity pattern. Form is not a bilinear form

Here my code:

import ufl
from dolfinx import fem, io, mesh, plot,  nls, log
from ufl import div, dx, grad, inner

from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc


from scipy.optimize import minimize
import scipy.linalg
import math
import logging
 
from logging.handlers import RotatingFileHandler

msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (1.0, 1.0)), n=(40, 40),
                            cell_type=mesh.CellType.triangle,)

nu = 0.001

# Function to mark x = 0, x = 1, y = 0 
def noslip_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)), np.isclose(x[1], 0.0))

# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)

# Lid velocity#
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))

P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V, Q = fem.FunctionSpace(msh, P2), fem.FunctionSpace(msh, P1)

# Create the function space
TH = P2 * P1
W = fem.FunctionSpace(msh, TH)
W0, _ = W.sub(0).collapse()

# No slip boundary condition
noslip = fem.Function(V)
facets = mesh.locate_entities_boundary(msh, 1, noslip_boundary)
dofs = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = fem.dirichletbc(noslip, dofs, W.sub(0))

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = fem.Function(W0)
lid_velocity.interpolate(lid_velocity_expression)
facets = mesh.locate_entities_boundary(msh, 1, lid)
dofs = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc1 = fem.dirichletbc(lid_velocity, dofs, W.sub(0))

# Since for this problem the pressure is only determined up to a constant, we pin the pressure at the point (0, 0)
zero = fem.Function(Q)
zero.x.set(0.0)
dofs = fem.locate_dofs_geometrical(
    (W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = fem.dirichletbc(zero, dofs, W.sub(1))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1, bc2]

# Define variational problem
w = fem.Function(W)

(u, p) = ufl.split(w)
(v, q) = ufl.TestFunctions(W)

log.set_log_level(log.LogLevel.INFO)

u_k = fem.Function(V)
def initU(x):
    k = [1.0, 0.0]
    return np.vstack((
        x[0] * k[0],
        x[1] * k[1]
    ))
u_k.interpolate(initU)  # previous (known) u

# Tentative velocity step
abilin = inner(grad(u)*u_k, v)*dx + nu*inner(grad(u), grad(v))*dx - div(v)*p*dx - q*div(u)*dx

Fext = fem.Function(V)
Fext.interpolate(initU) 

ZeroC = fem.Constant(msh,0.0)

Lrhs = inner(Fext,v) * dx + ZeroC*q*dx

problem_loc = fem.petsc.LinearProblem(abilin, Lrhs, bcs=bcs)
w_loc = problem_loc.solve()

Any ideas?
Thanks

2

If you want to solve a linear problem, you need to use TrialFunctions not a function in the definition of abilin.

3

Sorry for that and thank you…