Hello everyone,
I wrote this simple code and everything works okay, but at the end when it comes to use matplotlib I get an error concerning PETSc
import matplotlib.pyplot as plt
import pyvista
import ufl
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, div, FacetNormal,VectorElement, inner
import numpy as np
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, io, plot
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
def V_exact(mode):
#return lambda x: mode.atan(mode.pi * x[1])
# return lambda x: (10000 -8000 * x[1]) * x[1]
return lambda x: -mode.cos(mode.pi*x[1])
V_numpy = V_exact(np) # which will be used for interpolation
V_ufl = V_exact(ufl) # which will be used for defining the source term
def solve_poisson(N=10, degree=1):
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
[N, N], mesh.CellType.triangle)
x = SpatialCoordinate(domain)
sigma = fem.Constant(domain, PETSc.ScalarType(500))
f = div(-sigma * grad(V_ufl(x)))
W = fem.FunctionSpace(domain, ("Lagrange", 1))
V = TrialFunction(W)
csi = TestFunction(W)
a = dot(sigma * grad(V), grad(csi)) * dx
L = f * csi * dx
# We start by identifying the facets contained in each boundary and create a custom integration measure ds
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], 1)),
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], 1))]
# We now loop through all the boundary conditions and create MeshTags identifying the facets for each boundary condition
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, 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 = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
# We can then inspect individual boundaries using the Threshold-filter in Paraview
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
# Dirichlet condition
facets = facet_tag.find(3)
dofs = fem.locate_dofs_topological(W, fdim, facets)
facets2 = facet_tag.find(4)
dofs2 = fem.locate_dofs_topological(W, fdim, facets2)
BCs = [fem.dirichletbc(PETSc.ScalarType(-1), dofs, W), fem.dirichletbc(PETSc.ScalarType(1), dofs2, W)]
default_problem = fem.petsc.LinearProblem(a, L, bcs=BCs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
return default_problem.solve(), V_ufl(x)
def error_L2_func(Vh, V_ex, degree_raise=3):
# Create higher order function space
degree = 1 #Vh.function_space.ufl_element().degree
family = Vh.function_space.ufl_element().family_name
mesh = Vh.function_space.mesh
Q = fem.functionspace(mesh, (family, degree + degree_raise))
# Interpolate approximate solution
V_W = fem.Function(Q)
V_W.interpolate(Vh)
# Interpolate exact solution, special handling if exact solution
# is a ufl expression or a python lambda function
V_ex_W = fem.Function(Q)
if isinstance(V_ex, ufl.core.expr.Expr):
u_expr = fem.Expression(V_ex, Q.element.interpolation_points)
V_ex_W.interpolate(u_expr)
else:
V_ex_W.interpolate(V_ex)
# Compute the error in the higher order function space
e_W = fem.Function(Q)
e_W.x.array[:] = V_W.x.array - V_ex_W.x.array
# Integrate the error
error = fem.form(ufl.inner(e_W, e_W) * ufl.dx)
error_local = fem.assemble_scalar(error)
error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
return np.sqrt(error_global)
N = [10, 20, 40, 80, 160]
error_L2 = []
error_H1 = []
h = []
mpi_rank = 5
for i in range(len(N)):
Vh, Vex = solve_poisson(N[i])
comm = Vh.function_space.mesh.comm
error_L2 += [error_L2_func(Vh, V_numpy)]
eh = Vh - Vex
error_H10 = fem.form(dot(grad(eh), grad(eh)) * dx)
E_H10 = np.sqrt(comm.allreduce(fem.assemble_scalar(error_H10), op=MPI.SUM))
error_H1 += [E_H10]
h += [1. / N[i]]
if comm.rank == 0:
mpi_rank = comm.rank
print(f"h: {h[i]:.2e} Error: {error_L2[i]:.2e}")
if mpi_rank == 0:
plt.figure(figsize=(10, 6))
plt.loglog(N, error_L2, label='$L^2$ error')
plt.loglog(N, error_H1, label='$H^1$ error')
plt.loglog(N, h)
h_square = [x**2 for x in h]
plt.loglog(N, h_square)
plt.xlabel('N')
plt.ylabel('Error')
plt.legend()
plt.grid(True)
plt.show()
the erro
[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 https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[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.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0
can someone tell me what is wrong? I look also at other posts related to matplotlib, but this problem is not related to plotting a mesh or something more involved.
Thanks in advance.
Michele.