Seg fault when using mesh from gmshio.read_from_msh

General
1

Hello,

I’m trying to use a .msh file generated from a .geo file, reading it in using gmshio.read_from_msh. Everything works until it is time to solve the problem (Stokes flow), when the program crashes giving ‘Segmentation fault: 11’ or ‘Bus error:10’. Output below:

Info    : Reading 'basic_box_pillar_from_2d.msh'...
Info    : 33 entities
Info    : 12397 nodes
Info    : 14688 elements
Info    : Done reading 'basic_box_pillar_from_2d.msh'
Setting up variational problem
Solving Stokes flow
Bus error: 10

I have seen a similar issue in another post:

And I guess I’m looking to see if my error can be replicated in other folks’ installations of fenicsx.

The gmsh (4.13) and fenicsx (v0.8) are shown below:

GMSH (.geo)

SetFactory("OpenCASCADE");

// Element parameters
lc_coarse = 0.1;
lc_med = 0.08;
lc_fine = 0.04;

// Geometry
L = 6.0;
W = 1.0;
H = 0.25;

// Number of layers in extrusion
n_layers = Round(H/lc_fine);

x_min = -2.0;
x_max = x_min + L;

y_min = -0.5*W;
y_max = 0.5*W;

// Build geometry
Rectangle(1) = {x_min, y_min,  0, L, W, 0};
Circle(5) = {0, 0, 0, 0.25, 0, 2*Pi};

Curve Loop(2) = {3, 4, 1, 2};
Curve Loop(3) = {5};

// Form plane with hole for cylinder
Plane Surface(2) = {2, 3};

Recursive Delete {
  Curve{3}; 
}
Recursive Delete {
  Surface{1}; 
}

// Set element size over entire domain
Field[1] = Box;
Field[1].VIn = lc_med;
Field[1].VOut = lc_coarse;
Field[1].XMax = x_max*1.1;
Field[1].XMin = x_min*1.1;
Field[1].YMax = y_max*1.1;
Field[1].YMin = y_min*1.1;

// Set refined element size near cylinder
Field[2] = Box;
Field[2].VIn = lc_fine;
Field[2].VOut = lc_med;
Field[2].XMax = 1;
Field[2].XMin = -0.5;
Field[2].YMax = y_max*1.1;
Field[2].YMin = y_min*1.1;

// Set background mesh
Field[3] = Min;
Field[3].FieldsList = {1, 2};
Background Field = 3;

Mesh.Algorithm = 8;
Recombine Surface{2};

Extrude {0, 0, H} {
  Surface{2}; Layers {n_layers}; Recombine;
}

Physical Volume("fluid") = {1};
Physical Surface("inlet") = {4};
Physical Surface("outlet") = {6};
Physical Surface("walls") = {2,3,5,7,8};

FENICSX CODE (using an msh file from the above called ‘basic_box_pillar_from_2d.msh’

#!/usr/bin/env python

import os

os.environ['OMP_NUM_THREADS'] = '1'

import numpy as np
from basix.ufl import element, mixed_element
import ufl
import dolfinx
from dolfinx import fem, io, mesh, plot, log
from dolfinx.fem import Constant, Function, form, Expression
from dolfinx.fem import functionspace
from dolfinx.fem.petsc import LinearProblem 
from dolfinx.io import XDMFFile, gmshio
from ufl import dx, grad, inner, TestFunctions, TrialFunctions, dot, div

from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType  

comm = MPI.COMM_WORLD

# Read in mesh
msh, _, ft = gmshio.read_from_msh('basic_box_pillar_from_2d.msh', comm, 0, 3)
ft.name = 'Facet markers'

if comm.rank == 0:
    print('Setting up variational problem', flush=True)

# P1-P1 elements with GLS/Grad-Div stabilization
P1_v = element('Lagrange', msh.basix_cell(), 1, shape=(msh.geometry.dim,))
P1_p = element('Lagrange', msh.basix_cell(), 1)

MX = mixed_element([P1_v, P1_p])
W = functionspace(msh, MX)

# No slip condition on walls (meshtag 4)
W0,_ = W.sub(0).collapse()
noslip = Function(W0)
wall_dofs = fem.locate_dofs_topological((W.sub(0), W0), msh.topology.dim-1, ft.find(4))
bc_walls = fem.dirichletbc(noslip, dofs=wall_dofs, V=W.sub(0))

# Use inlet of all 1's for now
def inlet_velocity_expr(x):
    zero_val = np.zeros(x.shape[1], dtype=np.float64)
    one_val = np.ones(x.shape[1], dtype=np.float64)
    return (one_val, zero_val, zero_val)

# Inlet velocity of 1's for now (meshtag 2)
inlet_velocity = Function(W0)
inlet_velocity.interpolate(inlet_velocity_expr)
inlet_dofs = fem.locate_dofs_topological((W.sub(0), W0), msh.topology.dim-1, ft.find(2))
bc_inlet = fem.dirichletbc(inlet_velocity, dofs=inlet_dofs, V=W.sub(0))

# Outlet pressure condition of 0 (meshtag 3)
W0,_ = W.sub(1).collapse()
outflow_pressure = Function(W0)
outlet_dofs = fem.locate_dofs_topological((W.sub(0), W0), msh.topology.dim-1, ft.find(3))
bc_outlet = fem.dirichletbc(outflow_pressure, dofs=outlet_dofs, V=W.sub(1))

bcs = [bc_walls, bc_inlet, bc_outlet]

# STOKES FLOW
W0 = W.sub(0)
Q, _ = W0.collapse()
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = Function(Q)

# GLS Stabilization parameter
h = ufl.CellDiameter(msh)
Beta = 0.1 # Seems ok for 3D
mu_T = Beta*h*h # Stabilization coefficient

# Stokes flow residual
a = inner(grad(u), grad(v)) * dx
a -= inner(p, div(v)) * dx
a += inner(div(u), q) * dx
L = inner(f,v) * dx

# Stabilization terms
a += mu_T*inner(grad(p), grad(q)) * dx # GLS
L -= mu_T * inner(f, grad(q)) * dx # GLS
a += 1.0 * inner(div(u), div(v)) * dx # Grad-Div

problem = LinearProblem(a, L, bcs = bcs)

u_stokes = problem.solve()
2

Please note that you should massively simplify the reading in code if your segfault comes from reading in the mesh. Then we do not need all the other code after

and you can certainly reduce the number of imports

3

Hi Dokken, in fact, as I said in the original post, everything is fine in the code until it actually comes time to solve the problem. The segfault happens with u_stokes = problem.solve()

4

I’m sorry, I read your message to quickly.
I’ll get back to this later once I’m at a computer.

One thing you could try in the meantime is to check if it runs with docker.

5

I used the following docker command to create a docker terminal and test the code, and found the same result:

docker run -ti -v ${PWD}:/home/shared dolfinx/dolfinx:v0.8.0
6

Here you are mixing different sub spaces, alternating between sub(0) and sub(1)

7

You’ve done it again, Dokken! What a foolish mistake, thank you for spotting it.