Solve linear system in Mixed finite element method

I want to solve a Stokes problem in a unit square but now I get the wrong results.My code is as follows

``````from dolfinx import mesh,fem,plot,io
from mpi4py import MPI
from petsc4py import PETSc
VectorElement,FiniteElement,MixedElement,as_vector,
TrialFunction,TestFunction,split,SpatialCoordinate)
import numpy as np

domain=mesh.create_unit_square(MPI.COMM_WORLD,64,64,cell_type=mesh.CellType.triangle)
mu=fem.Constant(domain,1.0)
k=2 #degree of function space
Ve=VectorElement("CG",domain.ufl_cell(),k)
Pe=FiniteElement("CG",domain.ufl_cell(),k-1)
We=MixedElement([Ve,Pe])
W=fem.FunctionSpace(domain,We)

#boundary conditons
def no_slip(x):
return np.logical_or(np.logical_or(np.isclose(x[0],0),np.isclose(x[0],1)),
np.logical_or(np.isclose(x[1],0),np.isclose(x[1],1)))
V,_=W.sub(0).collapse()
noslip_v=fem.Function(V)
noslip_v.x.set(0)
facets=mesh.locate_entities_boundary(domain,1,no_slip)
bc=fem.dirichletbc(noslip_v,fem.locate_dofs_topological((W.sub(0),V),1,facets),W.sub(0))

#Exact solution
x=SpatialCoordinate(domain)
u3=(sin(pi*x[0])*sin(pi*x[1]))**2
u_exact=as_vector([u3.dx(1),-u3.dx(0)])
p_exact=sin(2.0*pi*x[0])*sin(3.0*pi*x[1])
I=Identity(domain.geometry.dim)
def sigma(u,p):
f=-div(sigma(u_exact,p_exact))

#Numerical solution
w=TrialFunction(W)
u,p=split(w)
ww=TestFunction(W)
v,q=split(ww)

L=inner(f,v)*dx

solution=fem.Function(W)
problem=fem.petsc.LinearProblem(a,L,bcs=[bc],petsc_options={"ksp_type":"preonly","pc_type":"lu"})
solution=problem.solve()
u,p=split(solution)

print(H1_error)
``````

The result I get is H1_error = nan.So I doubt that there is something wrong when I use LinearProblem in this mixed FEM scheme ,maybe the ksp_type and pc_type involved here but I am not sure.
Also, I have some other questions:
1.when writing u3, I import ufl.sin,ufl.pi, but when I use
`u3=(np.sin(np.pi*x[0])*np.sin(np.pi*x[1]))**2` ,it sends

``````AttributeError: 'Product' object has no attribute 'sin'

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
u3=(np.sin(np.pi*x[0])*np.sin(np.pi*x[1]))**2
TypeError: loop of ufunc does not support argument 0 of type Product which has no callable sin method
``````

2.If I want to know the matrix A from the bilinear form a and vector b from linear form L after assembling the system, is there something to show the values of them ,maybe in a list form, or array-like form? And is there a way to import them to matlab?

You need to have a clear separation of `ufl-objects` and `numpy-objects`. Anything related to `dolfinx.fem.Function`, `dolfinx.fem.Constant`, `ufl.as_vector` or `ufl.sin/cos` works on data in the ufl form language.

numpy works on arrays. Thus it cannot work on `SpatialCoordinate` as that is a ufl object.

Please also have a look at Stokes equations using Taylor-Hood elements — DOLFINx 0.7.0.0 documentation to make sure that you haven’t made any trivial errors.

Thanks.The original code is from https://github.com/david-kamensky/mae-207-fea-for-coupled-problems/blob/master/incompressibility/taylor-hood.py
as follows

``````from dolfin import *
N = 64
k = 2
mu = Constant(1.0)
mesh = UnitSquareMesh(N,N)
x = SpatialCoordinate(mesh)
I = Identity(mesh.geometry().dim())
cell = mesh.ufl_cell()
Ve = VectorElement("CG",cell,k)
Qe = FiniteElement("CG",cell,k-1)
W = FunctionSpace(mesh,MixedElement([Ve,Qe]))

# Take the curl of a potential for solution
# to manufacture:
u3 = (sin(pi*x[0])*sin(pi*x[1]))**2
u_exact = as_vector([u3.dx(1),-u3.dx(0)])
p_exact = sin(2.0*pi*x[0])*sin(3.0*pi*x[1])
def sigma(u,p):
f = -div(sigma(u_exact,p_exact))

# u_exact is zero on boundaries.
bc = DirichletBC(W.sub(0),Constant((0,0)),
"on_boundary")

# Galerkin formulation:
u,p = TrialFunctions(W)
v,q = TestFunctions(W)
- dot(f,v))*dx
up = Function(W)
solve(lhs(res)==rhs(res),up,bc)

# Check H1 error; can verify optimal (k-th) order
# by modifying parameter N at top of script.
u,p = split(up)
``````

and I want to change this one to a Dolfinx code but I get the wrong results.So in his dolfin code , what is the default ksp and pc type in the solve function? I wonder maybe the problem is here but I am not sure.Please help.

Now I try to use
`problem=fem.petsc.LinearProblem(a,L,bcs=[bc],petsc_options={"ksp_type":"gmres","pc_type":"jacobi"})`

in my Dolfinx version and gets the following results

``````N=32 : H1_error=0.04002088634138398
N=64 : H1_error=0.010022885754297573
N=128 : H1_error=0.0025174172982043025
``````

and results from dolfin version

``````Solving linear variational problem.
N=32:H1_error=0.04002089708412184
Solving linear variational problem.
N=64:H1_error=0.010021651805481696
Solving linear variational problem.
N=128:H1_error=0.0025064461672576867
``````

So the two results seem like having little difference but my dolfinx code takes a much long time when solving N=128, so maybe I should change the linear algebra solver.
Also, when changing different ksp_type, the results are also different.So still the question about what the default solver in the dolfin version of the code for `solve(lhs(res)==rhs(res),up,bc)`?

Legacy dolfin defaults to `umfpack`, ref: Bitbucket

Thus setting:

``````solution = fem.Function(W)
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={
"ksp_type": "preonly", "pc_type": "lu",  "pc_factor_mat_solver_type": "umfpack"})
``````

Solving the 128x128 system takes `3.5 seconds` for me using docker (`dolfinx/dolfinx:nightly`), while the legacy code (using ` ghcr.io/scientificcomputing/fenics:2023-02-17`) you produced takes 4.18 seconds to run the solve command (the timing is a bit unfair, as solve(a==L) also creates the matrices), so a larger time than DOLFINx is expected. However, as you can see, the solve time is comparable.

Really helps!Thank you very much.
Sorry for another question: when I use N=256 with the same parameters as yours

but this time the solver fails and sends error

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

Traceback (most recent call last):
File "/home//fenics/Numerical_methods_for_NS/Taylor_Hood_Solver.py", line 56, in <module>
solution=problem.solve()
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 618, in solve
self._solver.solve(self._b, self._x)
File "PETSc/KSP.pyx", line 395, in petsc4py.PETSc.KSP.solve
petsc4py.PETSc.Error: error code 76
[0] KSPSolve() line 1085 in ./src/ksp/ksp/interface/itfunc.c
[0] KSPSolve_Private() line 850 in ./src/ksp/ksp/interface/itfunc.c
[0] KSPSetUp() line 406 in ./src/ksp/ksp/interface/itfunc.c
[0] PCSetUp() line 1017 in ./src/ksp/pc/interface/precon.c
[0] PCSetUp_LU() line 131 in ./src/ksp/pc/impls/factor/lu/lu.c
[0] MatLUFactorNumeric() line 3194 in ./src/mat/interface/matrix.c
[0] MatLUFactorNumeric_UMFPACK() line 192 in ./src/mat/impls/aij/seq/umfpack/umfpack.c
[0] Error in external library
[0] umfpack_UMF_numeric failed
``````

so is there a way to get out of this problem?My computer has 16GB memory, so maybe it’s too large to be solved.

I try to use
“pc-factor-mat-solver_type”:”mumps”
and this time it is ok.