Segmentation violation while assembling matrices

dolfinx
, , , ,
1

Hello Everyone,
I am trying to solve eigenvalue problem J u=Lambda*u
where J is Jacobian of nonlinear convective diffusive reactive equations.
My code is as follows:

import basix.ufl
import dolfinx.fem
import dolfinx.fem.petsc
import numpy.typing
import petsc4py.PETSc
import ufl
import dolfinx.mesh
from mpi4py import MPI
from dolfinx.io import (VTXWriter, gmshio)
import numpy as np
gdim = 2
mesh_comm = MPI.COMM_WORLD
dtype0=np.float64
print("using dolfinx ",dolfinx.__version__," and \n data type",
      [dolfinx.default_scalar_type,petsc4py.PETSc.ScalarType])

#rectangular mesh
mesh=dolfinx.mesh.create_rectangle(mesh_comm, np.array([[0,-15.0],[15,15]]),np.array([65,129]),
                             cell_type=dolfinx.mesh.CellType.quadrilateral,dtype=dtype0)

# Function spaces
Vt_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)

W_element = basix.ufl.mixed_element([Vt_element,Vt_element,Vt_element])
W = dolfinx.fem.functionspace(mesh, W_element)

#functions
yg=dolfinx.fem.Function(W,dtype=dtype0)
vg=ufl.TestFunction(W)
yb=ufl.TrialFunction(W)

(yf,yx,yt)=ufl.split(yg)
(yfb,yxb,ytb) = ufl.split(yb)
(vf, vx, vt) = ufl.split(vg)

##paramters
LeF=petsc4py.PETSc.ScalarType(1.6)
LeX=petsc4py.PETSc.ScalarType(1.6)
Da=petsc4py.PETSc.ScalarType(100.0)
beta=petsc4py.PETSc.ScalarType(10.0)
gamma=petsc4py.PETSc.ScalarType(5.0)
phi=petsc4py.PETSc.ScalarType(1.0)
InvLeF=1.0/LeF
InvLeX=1.0/LeX

#governing equations in variational form
w=Da*(beta**3)*yf*yx*ufl.exp(beta*(yt-1)*(1+gamma)/(1+gamma*yt));
F=[(ufl.inner(yf.dx(0),vf)+InvLeF*ufl.inner(ufl.grad(yf),ufl.grad(vf))+ufl.inner(w,vf))*ufl.dx,
   (ufl.inner(yx.dx(0),vx)+InvLeX*ufl.inner(ufl.grad(yx),ufl.grad(vx))+phi*ufl.inner(w,vx))*ufl.dx,
   (ufl.inner(yt.dx(0),vt)+ufl.inner(ufl.grad(yt),ufl.grad(vt))-(1+phi)*ufl.inner(w,vt))*ufl.dx];

#Jacobian 
J = [[ufl.derivative(F[0], yf, yfb), ufl.derivative(F[0], yx, yxb),ufl.derivative(F[0], yt, ytb)],
[ufl.derivative(F[1], yf, yfb), ufl.derivative(F[1], yx, yxb),ufl.derivative(F[1], yt, ytb)],
     [ufl.derivative(F[2], yf, yfb), ufl.derivative(F[2], yx, yxb),ufl.derivative(F[2], yt, ytb)]]
rhs= [[-ufl.inner(yf, vf) * ufl.dx,None,None],[None,ufl.inner(yx, vx) * ufl.dx,None],
          [None,None,ufl.inner(yt, vt) * ufl.dx]]
##define x=0 inlet bc
def inlet_bc(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
    return np.isclose(x[0], 0)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, inlet_bc)

#degree of freedom boundary
(W0,W1,W2) = (W.sub(0),W.sub(2),W.sub(2))
YF, _ = W0.collapse()
bdofs_YF = dolfinx.fem.locate_dofs_topological((W0, YF), mesh.topology.dim - 1, boundary_facets)[0];
YX, _ = W1.collapse()
bdofs_YX = dolfinx.fem.locate_dofs_topological((W1, YX), mesh.topology.dim - 1, boundary_facets)[0];
YT, _ = W2.collapse()
bdofs_YT = dolfinx.fem.locate_dofs_topological((W2, YT), mesh.topology.dim - 1, boundary_facets)[0];

##generate inlet values with tanh profile
import pandas as pd
import numpy as np
from scipy import interpolate
y = np.linspace(-15, 15, 481)
dat = np.tanh(y)
yf_inlet = interpolate.interp1d(y, dat,fill_value="extrapolate")
yx_inlet = interpolate.interp1d(y, 1.0-dat,fill_value="extrapolate")  ##due to equal Lewis number case


#def func for interpolation in boundary
def yf_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the interpolated fuel mass fraction profile at the inlet."""
    values = np.zeros((x.shape[1]))
    values[:] = yf_inlet(x[1])
    return values
def yx_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the interpolated oxiidser mass fraction profile at the inlet."""
    values = np.zeros((x.shape[1]))
    values[:] = yx_inlet(x[1])
    return values
def yt_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the zero temperature profile at the inlet."""
    values = np.zeros((x.shape[1]))
    return values

# Boundary conditions
#fuel
yf_in=dolfinx.fem.Function(YF,dtype=dtype0)
yf_in.interpolate(yf_in_eval)
bcf=dolfinx.fem.dirichletbc(yf_in, bdofs_YF)
#oxy
yx_in=dolfinx.fem.Function(YX,dtype=dtype0)
yx_in.interpolate(yx_in_eval)
bcx=dolfinx.fem.dirichletbc(yx_in, bdofs_YX)

#temp
yt_in=dolfinx.fem.Function(YT,dtype=dtype0)
yt_in.interpolate(yt_in_eval)
bct=dolfinx.fem.dirichletbc(yt_in, bdofs_YT)
#combined bcs
bc=[bcf,bcx,bct]

# Assemble lhs and rhs matrices
print("assembling A matrices")
A = dolfinx.fem.petsc.assemble_matrix_block(dolfinx.fem.form(J), bcs=bc)
 #   A.assemble()
print("assembled lhs matrix successfully ")
print(" Assembling rhs matrix...")
B = dolfinx.fem.petsc.assemble_matrix_block(dolfinx.fem.form(rhs), bcs=bc)
B.assemble()
print("assembled rhs matrix successfully ")

when I run this, I get error at assembly of rhs matrix
The full log is:
using dolfinx 0.8.0 and
data type [<class ‘numpy.float64’>, <class ‘numpy.float64’>]
assembling A matrices
assembled lhs matrix successfully
Assembling rhs matrix…
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see FAQ — PETSc 3.22.1 documentation
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.

MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 59.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
Any help is appreciated. Thanks

2

Looks like you’re accidentally using the wrong variables in your RHS formulation. Consider:

rhs= [[-ufl.inner(yfb, vf) * ufl.dx,None,None],
      [None,ufl.inner(yxb, vx) * ufl.dx,None],
      [None,None,ufl.inner(ytb, vt) * ufl.dx]]

Where you’ll notice I’ve replaced the Functions from

yg=dolfinx.fem.Function(W,dtype=dtype0)
(yf,yx,yt)=ufl.split(yg)

with the TrialFunctions from

yb=ufl.TrialFunction(W)
(yfb,yxb,ytb) = ufl.split(yb)
3

Thanks for pointing out. I am able to fix this error.
However, I want to replace the unknown functions (yf,yx,yt) defined by

yg=dolfinx.fem.Function(W,dtype=dtype0)
(yf,yx,yt)=ufl.split(yg)

which appear in functional F and its Jacobian J
as

#governing equations in variational form
w=Da*(beta**3)*yf*yx*ufl.exp(beta*(yt-1)*(1+gamma)/(1+gamma*yt));
F=[(ufl.inner(yf.dx(0),vf)+InvLeF*ufl.inner(ufl.grad(yf),ufl.grad(vf))+ufl.inner(w,vf))*ufl.dx,
   (ufl.inner(yx.dx(0),vx)+InvLeX*ufl.inner(ufl.grad(yx),ufl.grad(vx))+phi*ufl.inner(w,vx))*ufl.dx,
   (ufl.inner(yt.dx(0),vt)+ufl.inner(ufl.grad(yt),ufl.grad(vt))-(1+phi)*ufl.inner(w,vt))*ufl.dx];

#Jacobian 
J = [[ufl.derivative(F[0], yf, yfb), ufl.derivative(F[0], yx, yxb),ufl.derivative(F[0], yt, ytb)],
[ufl.derivative(F[1], yf, yfb), ufl.derivative(F[1], yx, yxb),ufl.derivative(F[1], yt, ytb)],
     [ufl.derivative(F[2], yf, yfb), ufl.derivative(F[2], yx, yxb),ufl.derivative(F[2], yt, ytb)]]

. So I want to replace these unknown function by known function yfbase,yxbase and ytbase (these are base flows at which the Jacobian is taken for stability analysis)
these are defined in python using def.
I tried with
lhs=ufl.replace(J,{yf:yfbase,yx:yxbase,yt:ytbase})
but I am getting error: ValueError: Not an UFL type: <class ‘list’>
How do I achieve this ?
Any help appreciated. Thank You

4

If you remove the «unknowns» from the jacobian (ie replace a trial function with a function), it is no longer a matrix, but a vector. This can’t be sent into a NewtonSolver.

I can’t quite wrap my head around how your Newton method should look like.

5

This is for eigen value calculation and not Newton solver. Thanks

6

Still, for an eigenvalue calculation, you need two matrices, not a matrix and a vector, which is what you would get if you replace the test-function/trial-function by a known coefficient.

7

Thanks for pointing out. I have used known function, loaded from a file.
I have put the complete code for your reference. It is working ok. Still, I would appreciate if you would check if it is doing what it intends to do. i,e. calculating the eigenvalues for Linear stability analysis.
the data file can be found in
BF2D_Le1.6_phi1_z0.5_da100.dat

import basix.ufl
import dolfinx.fem
import dolfinx.fem.petsc
import numpy.typing
import petsc4py.PETSc
import scipy.special
import slepc4py.SLEPc
import ufl
import dolfinx.mesh
from mpi4py import MPI
from dolfinx.io import (VTXWriter, gmshio)
import numpy as np
gdim = 2
mesh_comm = MPI.COMM_WORLD
dtype0=np.float64
print("using dolfinx ",dolfinx.__version__," and \n data type",
      [dolfinx.default_scalar_type,petsc4py.PETSc.ScalarType])

#rectangular mesh
mesh=dolfinx.mesh.create_rectangle(mesh_comm, np.array([[0,-15.0],[15,15]]),np.array([16,32]),
                             cell_type=dolfinx.mesh.CellType.quadrilateral,dtype=dtype0)


# Function spaces
Vt_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)

W_element = basix.ufl.mixed_element([Vt_element,Vt_element,Vt_element])
W = dolfinx.fem.functionspace(mesh, W_element)

#functions
yb=ufl.TrialFunction(W)
(yfb,yxb,ytb) = ufl.split(yb)
yg=dolfinx.fem.Function(W,dtype=dtype0)


##base flow
##load base flows
import pandas as pd
import numpy as np
from scipy.interpolate import RegularGridInterpolator

# load full data file
df = pd.read_csv('./BF2D_Le1.6_phi1_z0.5_da100.dat', sep='\s+',header=None)
##define x and y array
x = np.linspace(0, 15, 241)
y = np.linspace(-15, 15, 481)
# Convert the DataFrame to a NumPy array (matrix)
matrix = df.values
#extract values
temp=matrix[:,2].reshape(481,241)
fuel=matrix[:,3].reshape(481,241)
oxy=matrix[:,4].reshape(481,241)
Yt_base = RegularGridInterpolator((x,y), np.transpose(temp),bounds_error=False,fill_value=None)
Yf_base = RegularGridInterpolator((x, y), np.transpose(fuel),bounds_error=False,fill_value=None)
Yx_base = RegularGridInterpolator((x, y), np.transpose(oxy),bounds_error=False,fill_value=None)


##form numpy functoions
#def func for interpolation of base flow
def yf_base_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the interpolated fuel mass fraction profile at the inlet."""
    values = np.zeros((x.shape[1]))
    values[:] = Yf_base(np.transpose([x[0],x[1]]))
    return values
def yx_base_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the interpolated fuel mass fraction profile at the inlet."""
    values = np.zeros((x.shape[1]))
    values[:] = Yx_base(np.transpose([x[0],x[1]]))
    return values
def yt_base_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the interpolated fuel mass fraction profile at the inlet."""
    values = np.zeros((x.shape[1]))
    values[:] = Yt_base(np.transpose([x[0],x[1]]))
    return values


##interpolate the base_flow onto the yg function
yg.sub(0).interpolate(yf_base_eval)
yg.sub(1).interpolate(yx_base_eval)
yg.sub(2).interpolate(yt_base_eval)

(yf,yx,yt)=ufl.split(yg)
vg=ufl.TestFunction(W)
(vf, vx, vt) = ufl.split(vg)

##paramters
LeF=petsc4py.PETSc.ScalarType(1.6)
LeX=petsc4py.PETSc.ScalarType(1.6)
Da=petsc4py.PETSc.ScalarType(100.0)
beta=petsc4py.PETSc.ScalarType(10.0)
gamma=petsc4py.PETSc.ScalarType(5.0)
phi=petsc4py.PETSc.ScalarType(1.0)
InvLeF=1.0/LeF
InvLeX=1.0/LeX

#governing equations in variational form
w=Da*(beta**3)*yf*yx*ufl.exp(beta*(yt-1)*(1+gamma)/(1+gamma*yt));
F=(ufl.inner(yf.dx(0),vf)+InvLeF*ufl.inner(ufl.grad(yf),ufl.grad(vf))+ufl.inner(w,vf))*ufl.dx\
+(ufl.inner(yx.dx(0),vx)+InvLeX*ufl.inner(ufl.grad(yx),ufl.grad(vx))+phi*ufl.inner(w,vx))*ufl.dx\
+(ufl.inner(yt.dx(0),vt)+ufl.inner(ufl.grad(yt),ufl.grad(vt))-(1+phi)*ufl.inner(w,vt))*ufl.dx;

#Jacobian 
J = ufl.derivative(F, yg, yb)

rhs= ufl.inner(yfb, vf) * ufl.dx+ufl.inner(yxb, vx) * ufl.dx+ufl.inner(ytb, vt) * ufl.dx
#rhs= [ufl.inner(yf, vf) * ufl.dx,ufl.inner(yx, vx) * ufl.dx,ufl.inner(yt, vt) * ufl.dx]
##define x=0 inlet bc
def inlet_bc(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
    return np.isclose(x[0], 0)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, inlet_bc)

#degree of freedom boundary
(W0,W1,W2) = (W.sub(0),W.sub(2),W.sub(2))
YF, _ = W0.collapse()
bdofs_YF = dolfinx.fem.locate_dofs_topological((W0, YF), mesh.topology.dim - 1, boundary_facets)[0];
YX, _ = W1.collapse()
bdofs_YX = dolfinx.fem.locate_dofs_topological((W1, YX), mesh.topology.dim - 1, boundary_facets)[0];
YT, _ = W2.collapse()
bdofs_YT = dolfinx.fem.locate_dofs_topological((W2, YT), mesh.topology.dim - 1, boundary_facets)[0];


#def func for interpolation in boundary
def yf_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the interpolated fuel mass fraction profile at the inlet."""
    values = np.zeros((x.shape[1]))
    #values[:] = yf_inlet(x[1])
    return values
def yx_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the interpolated oxiidser mass fraction profile at the inlet."""
    values = np.zeros((x.shape[1]))
    #values[:] = yx_inlet(x[1])
    return values
def yt_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the zero temperature profile at the inlet."""
    values = np.zeros((x.shape[1]))
    return values

# Boundary conditions
#fuel
yf_in=dolfinx.fem.Function(YF,dtype=dtype0)
yf_in.interpolate(yf_in_eval)
bcf=dolfinx.fem.dirichletbc(yf_in, bdofs_YF)
#oxy
yx_in=dolfinx.fem.Function(YX,dtype=dtype0)
yx_in.interpolate(yx_in_eval)
bcx=dolfinx.fem.dirichletbc(yx_in, bdofs_YX)

#temp
yt_in=dolfinx.fem.Function(YT,dtype=dtype0)
yt_in.interpolate(yt_in_eval)
bct=dolfinx.fem.dirichletbc(yt_in, bdofs_YT)
#combined bcs
bc=[bcf,bcx,bct]

# Assemble lhs and rhs matrices
print("assembling A matrices")
A = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(J), bcs=bc)
A.assemble()
print("assembled lhs matrix successfully ")
print(" Assembling rhs matrix...")
B = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(rhs), bcs=bc)
B.assemble()
print("assembled rhs matrix successfully ")


##eigen problem
eps = slepc4py.SLEPc.EPS().create(mesh.comm)
eps.setOperators(A, B)
eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GNHEP)
eps.setDimensions(10, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)
eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.SMALLEST_REAL)
eps.getST().getKSP().setType("preonly")
eps.getST().getKSP().getPC().setType("lu")
eps.getST().getKSP().getPC().setFactorSolverType("mumps")
eps.solve()
assert eps.getConverged() >= 1


import multiphenicsx.fem
import multiphenicsx.fem.petsc
Print = petsc4py.PETSc.Sys.Print
its = eps.getIterationNumber()
Print("Number of iterations of the method: %i" % its)
sol_type = eps.getType()
Print("Solution method: %s" % sol_type)
nev, ncv, mpd = eps.getDimensions()
Print("Number of requested eigenvalues: %i" % nev)
tol, maxit = eps.getTolerances()
Print("Stopping condition: tol=%.4g, maxit=%d" % (tol, maxit))
nconv = eps.getConverged()
Print("Number of converged eigenpairs: %d" % nconv)

vr = dolfinx.cpp.fem.petsc.create_vector_block(
    [(space.dofmap.index_map, restriction_.dofmap.index_map_bs) for space in [YF,YX,YT]])
vi = dolfinx.cpp.fem.petsc.create_vector_block(
    [(space.dofmap.index_map, restriction_.dofmap.index_map_bs) for space in [YF,YX,YT]])
vrfun_global=[];
eig_global=[]
for i in range(nconv):
    print(i)
    eigv = eps.getEigenpair(i, vr, vi)
    er, ei = eigv.real, eigv.imag
    eig_global.append([er,ei])
    error = eps.computeError(i)
    if i != 0.0:
        Print(" %9f%+9f j  %12g" % (er, ei, error))
    else:
        Print(" %12f       %12g" % (r, error))
    Print("")
       # Transform eigenvector into an eigenfunction that can be plotted
    vr.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
    vr_fun = dolfinx.fem.Function(YF)
    with vr_fun.x.petsc_vec.localForm() as vr_fun_local, \
            multiphenicsx.fem.petsc.VecSubVectorWrapper(vr, YF.dofmap) as vr_wrapper:
        vr_fun_local[:] = vr_wrapper
    vrfun_global.append(vr_fun)