Facing problem with solver()

Hello sir

I try to use LU+GMRES but it shows following error. Help me in resolving this error.

from fenics import *
from fenics_adjoint import *

nx = 32
mesh = UnitCubeMesh(nx, nx, nx)

element = mesh.ufl_cell()

H    = VectorElement("Lagrange", element, 2)   # velocity subspace
Q    = FiniteElement("Lagrange", element, 1)   # pressure subspace
HxQ  = FunctionSpace(mesh, MixedElement([H,Q]))     # Mixed subspace

wh = Function(HxQ)
(uh, ph) = split(wh)
(vh, qh) = TestFunctions(HxQ)


u_in = Expression(("(16*0.45*x[1]*x[2]*(0.41-x[1])*(0.41-x[2]))/(0.41*0.41*0.41*0.41)","0","0"), domain = mesh, degree=6)
p_exact = Expression("0",domain=mesh, degree=4)
f = Constant((0,0,0))

nu  = Constant(1)

# Strain:
def epsilon(v):
    return 0.5*(grad(v) + grad(v).T)
    

B_S = (2.0*nu*inner(epsilon(uh), epsilon(vh)) + div(vh)*ph - qh*div(uh))*dx + dot(grad(uh)*uh,vh)*dx
L_S = inner(f,vh)*dx

#u_n = project(u_exact, Z.sub(0).collapse(), solver_type="mumps", form_compiler_parameters=None)
tol=1E14
def boundary(x, on_boundary):
    return on_boundary and (near(x[1],0,tol) or near(x[1],1,tol) or near(x[2],0,tol) or near(x[2],1,tol))
def boundary_D(x, on_boundary):
    return on_boundary and (near(x[0],0,tol))
bc1 = DirichletBC(HxQ.sub(0), Constant((0,0,0)), boundary)
bc2 = DirichletBC(HxQ.sub(0), u_in, boundary_D)
bc  = [bc1, bc2]
F = B_S - L_S
# Assemble system.
a = lhs(F)
L = rhs(F)
A, b = assemble_system(a, L, bc)

MAX_ITERS = 10000
solver = PETScKrylovSolver('gmres','none')
solver.ksp().setGMRESRestart(MAX_ITERS)
solver.parameters["monitor_convergence"] = True
solver.parameters['relative_tolerance'] = 1e-6
solver.parameters['maximum_iterations'] = MAX_ITERS


solver.solve(A,wh.vector(),b)
(uh, ph) = wh.split(True)

error:

  A, b = assemble_system(a, L, bc)
  File "/home/ayush/.local/lib/python3.10/site-packages/fenics_adjoint/assembly.py", line 47, in assemble_system
    A, b = backend.assemble_system(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 423, in assemble_system
    assembler = cpp.fem.SystemAssembler(A_dolfin_form, b_dolfin_form, bcs)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to assemble system.
*** Reason:  expected a bilinear form for a.
*** Where:   This error was encountered inside SystemAssembler.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  ubuntu

When solving a linear problem, you need to use a TrialFunction for wh.

This is a non linear problem
look at the variational formulation:

B_S = (2.0*nu*inner(epsilon(uh), epsilon(vh)) + div(vh)*ph - qh*div(uh))*dx + dot(grad(uh)*uh,vh)*dx

if this is a non-linear problem, this does not make sense, as there is no lhs for a residual form.
You would need to form the jacobian.

For more info; see for instance
https://jsdokken.com/dolfinx-tutorial/chapter4/newton-solver.html

Using an iterative solver for the Navier-Stokes system is not trivial.

Consider for example the Stokes demo.

Sir, I want to use LU+GMRES. can you provide some insight in that direction?