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