Error with PETSc

Karlosjusus · September 25, 2023, 3:01pm

hi everyone.
I’m trying to solve a linear system with PETSc, but I get the following error


 File "/home/c.ruiz/Fenicx/Venturi/Codes/Venturi_adjoint_laminar_nobl.py", line 356, in <module>
    wa = problem.solve()
  File "/opt/ohpc/pub/apps/miniconda/miniconda3/envs/fenicsx-env-0.5.2/lib/python3.10/site-packages/dolfinx/fem/petsc.py", line 618, in solve
    self._solver.solve(self._b, self._x)
  File "PETSc/KSP.pyx", line 403, in petsc4py.PETSc.KSP.solve
petsc4py.PETSc.Error: error code 378
[0] KSPSolve() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/ksp/ksp/interface/itfunc.c:1078
[0] KSPSolve_Private() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/ksp/ksp/interface/itfunc.c:843
[0] KSPSetUp() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/ksp/ksp/interface/itfunc.c:407
[0] PCSetUp() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/ksp/pc/interface/precon.c:990
[0] PCSetUp_LU() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/ksp/pc/impls/factor/lu/lu.c:93
[0] MatLUFactorSymbolic() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/mat/interface/matrix.c:2966
[0] MatLUFactorSymbolic_SeqAIJ() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/mat/impls/aij/seq/aijfact.c:378
[0] PetscMallocA() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/sys/memory/mal.c:401
[0] PetscMallocAlign() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/sys/memory/mal.c:49
[0] MatLUFactorSymbolic_SeqAIJ

The code is

import dolfinx
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import math
import ufl
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, MixedElement)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, 
                               create_vector, create_matrix, set_bc)

from dolfinx.fem import (Constant, Function, FunctionSpace, 
                         assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)


with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "Mesh_nobl/venturi_nobl.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")
    cf = xdmf.read_meshtags(domain, name = 'Grid')

    
    
tdim = domain.topology.dim
tree = dolfinx.geometry.BoundingBoxTree(domain, tdim)
num_entities_local = domain.topology.index_map(tdim).size_local + domain.topology.index_map(tdim).num_ghosts
entities = np.arange(num_entities_local, dtype = np.int32)
mid_tree = dolfinx.geometry.create_midpoint_tree(domain, tdim,entities)
mesh_geom = domain.geometry.x
tol = 1e-8

def closest_point(mesh, point,tol):
    points = point
    while True:
        try:
            entity = dolfinx.geometry.compute_closest_entity(tree, mid_tree, mesh, points)
            break
        except RuntimeError:
            print(points)
            for j in range(len(mesh.geometry.x)):
                if (np.abs(mesh_geom[j][0]-points[0])< tol and np.abs(mesh_geom[j][1]-points[1])< tol and np.abs(mesh_geom[j][2]-points[2])< tol):
                     return points, j
    
    geom_dof = dolfinx.cpp.mesh.entities_to_geometry(mesh,tdim,[entity], False)
    #print('geom',geom_dof)
    mesh_nodes = mesh_geom[geom_dof][0]
    #print('mesh_nodes', mesh_nodes)
    dis = dolfinx.geometry.compute_distance_gjk(points, mesh_nodes)
    index = -100
    for i in range(len(mesh_nodes)):
        #print(mesh_nodes[i])
        if np.abs(mesh_nodes[i][0]-points[0]+ dis[0]) < tol and np.abs(mesh_nodes[i][1]-points[1]+ dis[1]) < tol and np.abs(mesh_nodes[i][2]-points[2]+ dis[2]) < tol :
            index = i
            
    if (index == -100):
        return point, index
    else:
        return points - dis , geom_dof[0][index]


vx_np = np.loadtxt('Data_nobl/data_vp.txt')
vx_aux = vx_np[:,1:4]

b = []
indices = []
ind = int(0)
for i in range(len(vx_aux)):
    a,j = closest_point(domain,vx_aux[i],tol)
    indices.append(int(j))





########################################################################################################

    # Definition of Spaces. (dim and polinomial type)

v_cg2 = VectorElement("CG", domain.ufl_cell(), 2)
v_cg1 = VectorElement("CG", domain.ufl_cell(), 1)
s_cg1 = FiniteElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg1)
Q = dolfinx.fem.FunctionSpace(domain, s_cg1)
Z = MixedElement([v_cg2, s_cg1])
W = dolfinx.fem.FunctionSpace(domain, Z)


# In[4]:


Lx  =0.945

def upstream(x): 
    return np.isclose(x[0],0)


def downstream(x):
    return np.isclose(x[0],Lx)

def walls(x):
    a = np.greater(x[0],0)
    b = np.less(x[0],Lx)
    return np.logical_and(np.greater(x[2]**2+x[1]**2, 0.00000001), np.logical_and(a,b))



fdim = domain.topology.dim -1
up_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, upstream)
down_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, downstream)
wall_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, walls)



V1, _ = W.sub(0).collapse()
Q1, _ = W.sub(1).collapse()

uuw_condition = Function(V1)
pd_condition = Function(Q1)




up_f = dolfinx.fem.locate_dofs_topological((W.sub(0),V1), fdim, up_facets)
w_f = dolfinx.fem.locate_dofs_topological((W.sub(0),V1), fdim, wall_facets)
d_f = dolfinx.fem.locate_dofs_topological((W.sub(1),Q1), fdim, down_facets)

bc_u = dirichletbc(uuw_condition, up_f,W.sub(0))
bc_w = dirichletbc(uuw_condition, w_f,W.sub(0))
bc_p = dirichletbc(pd_condition, d_f,W.sub(1))

bc = [bc_u, bc_w, bc_p]




boundaries = [(1,up_facets),(2,down_facets),(3,wall_facets)]

facet_indices, facet_markers = [], []
for (marker, face ) in boundaries:
    facet_indices.append(face)
    facet_markers.append(np.full_like(face, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])




ds = Measure("ds", domain=domain, subdomain_data=facet_tag)



    # Parameters
mu = Constant(domain, PETSc.ScalarType(0.0303))

    # Normal and operators
    
n = FacetNormal(domain)

def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))




dataname = 'Data_nobl/data_vp.txt'
data_np = np.loadtxt(dataname)

u = Function(V1)
p = Function(Q1)

u_aux = Function(V)

for i in range(len(indices)):  
    u_aux.x.array[3*indices[i]]=data_np[i,6]
    u_aux.x.array[3*indices[i]+1]=data_np[i,5]
    u_aux.x.array[3*indices[i]+2]=data_np[i,4]
    p.x.array[indices[i]]=data_np[i,7]

u.interpolate(u_aux)




print(len(p.x.array))    
print(len(indices))



wa = TrialFunction(W)
ua, pa = ufl.split(wa)
wtest = TestFunction(W)
va,qa = ufl.split(wtest)
f = Function(Q1)

x_upstream = 0.1
x_downstream = 0.8
x_throat = 0.3

   

x_up, ind_up =closest_point(domain,[x_upstream,0,0],0.001)

x_throat, ind_throat =closest_point(domain,[x_throat,0,0],0.001)

x_down, ind_down =closest_point(domain,[x_downstream,0,0],0.01)




delta_up = Function(Q1)
delta_up.name = "d_u"

delta_down = Function(Q1)
delta_down.name = "d_d"

delta_throat = Function(Q1)
delta_throat.name = "d_t"


delta_up.x.array[ind_up] = 1
delta_throat.x.array[ind_throat] = 1
delta_down.x.array[ind_down] = 1


# In[9]:


p_up = p.x.array[ind_up]
p_throat = p.x.array[ind_throat]
p_dn = p.x.array[ind_down]
umd = (p_up-p_dn)**2

di_ua = Function(Q1)

di_ua.x.array[:] = 1/(umd)*((p_throat - p_dn)*delta_up.x.array[:] + (p_dn-p_up)*delta_throat.x.array[:] + (p_up - p_throat)*delta_down.x.array[:] )


    
a1 = dot(dot(ua,nabla_grad(u).T),va)*dx - dot(dot(nabla_grad(ua).T,u),va)*dx + inner(sigma(ua,-pa),epsilon(va))*dx + dot(div(ua),qa)*dx + dot(dot(u,n)*ua,va)*ds - dot(pa*n,va)*ds - mu*dot(dot(nabla_grad(ua),n),va)*ds
a2 = dot(div(ua),qa)*dx 
#a2 = dot(div(ua)-f,qa)*dx 
a3 = dot(dot(u,n)*ua,va)*ds - dot(pa*n,va)*ds - mu*dot(dot(nabla_grad(ua),n),va)*ds
a0 = a1+a2+a3
L = dot(di_ua,qa)*dx

problem = dolfinx.fem.petsc.LinearProblem(a0, L, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
wa = problem.solve()

ua,pa = wa.split()
with dolfinx.io.XDMFFile(domain.comm, "Plot/ua.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(ua)
with dolfinx.io.XDMFFile(domain.comm, "Plot/pa.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(pa)

The files are here:

Any suggestions?
Thank you!

nate · September 25, 2023, 3:18pm

I don’t have the bandwidth to try and run your example, but are you sure you’re not running out of memory? You could also specify which LU solver you’re employing as UMFPACK, for example, has a very small memory limit as a default.

Karlosjusus · September 26, 2023, 9:32am

I am working with my university’s supercomputer, memory should not be a problem, how can I check what kind of LU I am applying?

dokken · September 26, 2023, 10:44am

Add "pc_factor_mat_solver_type": "mumps" to petsc_options to use something else than the default (“umfpack”).

Karlosjusus · October 4, 2023, 2:04pm

Thank you @dokken and @nate!

Topic		Replies	Views
Error on the cluster : PETScLUSolver.cpp Errors	3	568	November 23, 2021
Error: Unable to solve linear system using PETSc Krylov solver Errors	1	499	August 10, 2022
Problem with PETSc solver. error code 77 Linear Algebra	13	1653	August 8, 2021
Issue with PETSc SNES and LU Solver General	5	79	July 11, 2024
PETSc error when running code on a fine mesh Linear Algebra	0	441	October 7, 2022

Error with PETSc

Related topics