Solve linear system in Mixed finite element method

General
1

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
from ufl import (dx,inner,dot,grad,div,sym,Identity,sin,pi,Measure,
                 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):
    return 2.0*mu*sym(grad(u))-p*I
f=-div(sigma(u_exact,p_exact))

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

a=inner(sigma(u,p),grad(v))*dx+div(u)*q*dx
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)

H1_error=np.sqrt(fem.assemble_scalar(fem.form(inner(grad(u-u_exact),grad(u-u_exact))*dx)))
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?

Thanks in advance!

2

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.8.0.0 documentation to make sure that you haven’t made any trivial errors.

6

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):
    return 2.0*mu*sym(grad(u)) - p*I
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)
res = (inner(sigma(u,p),grad(v)) + div(u)*q
       - 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)
grad_e = grad(u-u_exact)
print(sqrt(assemble(inner(grad_e,grad_e)*dx)))

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.

7

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)?

8

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.

13

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.

14

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