Facing problem with Taylor Hood element due to increase in nodes for larger mesh

when I solve this problem with Taylor hood element then it doesn’t work out but If I am considering P1-P1 element and add a stabilization term in the weak formulation then the code works well. I think, it happens due to increasing nodes. Can you suggest something in resolving this issue ?

Hope to hear from you soon. Hope so someone provide suggestions on the basis of experience.
Minimal example:

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

solve(F == 0, wh, bc, solver_parameters={"newton_solver":{"linear_solver":'mumps'},"newton_solver":{"relative_tolerance":1e-7}})
(uh, ph) = wh.split(True)

error:

   No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 2.035e+04 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-07)

UMFPACK V5.7.9 (Oct 20, 2019): ERROR: out of memory

Traceback (most recent call last):
  File "/home/ayush/a.py", line 46, in <module>
    solve(F == 0, wh, bc, solver_parameters={"newton_solver":{"linear_solver":'mumps'},"newton_solver":{"relative_tolerance":1e-7}})
  File "/home/ayush/.local/lib/python3.10/site-packages/fenics_adjoint/solving.py", line 51, in solve
    output = backend.solve(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 233, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 314, in _solve_varproblem
    solver.solve()
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 successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside ./dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0

Consider an iterative solver. See e.g:

1 Like