Hi everyone,
I have built and tested my code using UnitSquareMesh( ) available in Fenics library and it is working properly. However, I need to work with more complex domains. So, I have chosen gmsh software for this purpose. I have saved the mesh in format .msh and convert it to .xdmf file or converting it directly from .geo to .xdmf .
When I use the mesh input with *.xdmf files, two error message shows up according to the tests that I have made:
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
Or this other error message :
HDF5-DIAG: Error detected in HDF5 (1.10.4) MPI-process 0:
#000: ../../../src/H5F.c line 444 in H5Fcreate(): unable to create file
major: File accessibilty
minor: Unable to open file
#001: ../../../src/H5Fint.c line 1364 in H5F__create(): unable to open file
major: File accessibilty
minor: Unable to open file
#002: ../../../src/H5Fint.c line 1615 in H5F_open(): unable to lock the file
major: File accessibilty
minor: Unable to open file
#003: ../../../src/H5FD.c line 1640 in H5FD_lock(): driver lock request failed
major: Virtual File Layer
minor: Can't update object
#004: ../../../src/H5FDsec2.c line 941 in H5FD_sec2_lock(): unable to lock file, errno = 11, error message = 'Resource temporarily unavailable'
major: File accessibilty
minor: Bad file ID accessed
HDF5-DIAG: Error detected in HDF5 (1.10.4) MPI-process 0:
#000: ../../../src/H5F.c line 444 in H5Fcreate(): unable to create file
major: File accessibilty
minor: Unable to open file
#001: ../../../src/H5Fint.c line 1364 in H5F__create(): unable to open file
major: File accessibilty
minor: Unable to open file
#002: ../../../src/H5Fint.c line 1615 in H5F_open(): unable to lock the file
major: File accessibilty
minor: Unable to open file
#003: ../../../src/H5FD.c line 1640 in H5FD_lock(): driver lock request failed
major: Virtual File Layer
minor: Can't update object
#004: ../../../src/H5FDsec2.c line 941 in H5FD_sec2_lock(): unable to lock file, errno = 11, error message = 'Resource temporarily unavailable'
major: File accessibilty
minor: Bad file ID accessed
HDF5-DIAG: Error detected in HDF5 (1.10.4) MPI-process 0:
#000: ../../../src/H5F.c line 444 in H5Fcreate(): unable to create file
major: File accessibilty
minor: Unable to open file
#001: ../../../src/H5Fint.c line 1364 in H5F__create(): unable to open file
major: File accessibilty
minor: Unable to open file
#002: ../../../src/H5Fint.c line 1615 in H5F_open(): unable to lock the file
major: File accessibilty
minor: Unable to open file
#003: ../../../src/H5FD.c line 1640 in H5FD_lock(): driver lock request failed
major: Virtual File Layer
minor: Can't update object
#004: ../../../src/H5FDsec2.c line 941 in H5FD_sec2_lock(): unable to lock file, errno = 11, error message = 'Resource temporarily unavailable'
major: File accessibilty
minor: Bad file ID accessed
I am proving the code with UnitSquareMesh( ) code (working) just for checking if necessary, the code with xdmf input (not working), and the code used to convert from .geo file to .xdmf.
Any suggestions would be very helpful.
xdmf Mesh:
def lmbdainv(s, mu_rel, no, nw):
return 1.0 / ((1.0 / mu_rel) * s ** nw + (1.0 - s) ** no)
# Fractional flow function
def F(s, mu_rel, no, nw):
return s ** nw / (s ** nw + mu_rel * (1.0 - s) ** no)
def tensor_jump(v, n):
return ufl.outer(v, n)("+") + ufl.outer(v, n)("-")
class PiecewiseConstant(UserExpression):
def __init__(self, values, markers, **kwargs):
self._values = values
self._markers = markers
super().__init__(**kwargs)
def eval_cell(self, values, x, cell):
values[0] = self._values[self._markers[cell.index]]
def value_shape(self):
return tuple()
### Relative viscosity of water w.r.t. crude oil
Kinv = Constant(1 / (1e-6))
# Time step
dt = Constant(0.01)
mu = 1
mu_b = 1
phi = Constant(0.4)
pin = 6000
pout = 5000
sbar = Constant(1)
t = 0.0
T = 1000 * float(dt)
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("domains.xdmf") as infile:
infile.read(mvc)
Markers = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc1 = MeshValueCollection("size_t", mesh, 1)
with XDMFFile("boundaries.xdmf") as infile:
infile.read(mvc1)
boundaries = cpp.mesh.MeshFunctionSizet(mesh, mvc1)
order = 1
BDM = FiniteElement("BDM", mesh.ufl_cell(), order)
DG = FiniteElement("DG", mesh.ufl_cell(), order - 1)
# DG1 = FiniteElement("DG", mesh.ufl_cell(), order)
Element = MixedElement([BDM, DG, DG])
mixed_space = FunctionSpace(mesh, Element)
n = FacetNormal(mesh)
# Function spaces and functions
V = TestFunction(mixed_space)
dU = TrialFunction(mixed_space)
U = Function(mixed_space)
U0 = Function(mixed_space)
v, q, r = split(V)
u, p, s = split(U)
u0, p0, s0 = split(U0)
s_mid = 0.5 * (s0 + s)
# ============= DEFINITION OF SPATIALLY-VARYING PARAMETERS =====================
marker_inner = 1
marker_outer = 0
no_outer = 2
nw_outer = 2
no_inner = 1
nw_inner = 1
no = {marker_inner: no_inner, marker_outer: no_outer}
nw = {marker_inner: nw_inner, marker_outer: nw_outer}
VVV = FunctionSpace(mesh, "DG", 0)
noo = PiecewiseConstant(no, Markers)
noo_proj = project(noo, VVV)
nww = PiecewiseConstant(nw, Markers)
nww_proj = project(nww, VVV)
# =========== END DEFINITION OF SPATIALLY-VARYING PARAMETERS ===================
ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dx = Measure("dx", domain=mesh, subdomain_data=Markers)
bc2 = DirichletBC(mixed_space.sub(0), Constant((0.0, 0.0)), boundaries, 2)
bc4 = DirichletBC(mixed_space.sub(0), Constant((0.0, 0.0)), boundaries, 4)
bcs = [bc2, bc4] # velocity BC
alpha = 35
h = CellDiameter(mesh)
h2 = ufl.Min(h("+"), h("-"))
stab = (
mu * (alpha / h2) * inner(tensor_jump(u, n), tensor_jump(v, n)) * dS
- mu * inner(avg(grad(u)), tensor_jump(v, n)) * dS
- mu * inner(avg(grad(v)), tensor_jump(u, n)) * dS
)
L1 = (
mu_b * inner(grad(u), grad(v)) * dx(1)
+ inner(v, lmbdainv(s_mid, mu, no_outer, nw_outer) * Kinv * u) * dx(0)
- div(v) * p * dx
+ pin * dot(v, n) * ds(1)
+ pout * dot(v, n) * ds(3)
+ stab
)
L2 = q * div(u) * dx(0) + q * div(u) * dx(1)
un = 0.5 * (inner(u0, n) + abs(inner(u0, n)))
un_h = 0.5 * (inner(u0, n) - abs(inner(u0, n)))
stabilisation = (
dt("+")
* inner(
jump(r),
un("+") * F(s_mid, mu, noo_proj, nww_proj)("+")
- un("-") * F(s_mid, mu, noo_proj, nww_proj)("-"),
)
* dS
)
L3 = (
phi * r * (s - s0) * dx(0)
+ r * (s - s0) * dx(1)
- dt * inner(grad(r), F(s_mid, mu, no_outer, nw_outer) * u) * dx(0)
- dt * inner(grad(r), F(s_mid, mu, no_outer, nw_outer) * u) * dx(1)
+ dt * r * F(s_mid, mu, no_outer, nw_outer) * un * ds
+ stabilisation
+ dt * r * un_h * sbar * ds(1)
)
# Total L
L = L1 + L2 + L3
J = derivative(L, U, dU)
problem = NonlinearVariationalProblem(L, U, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)
solver.parameters["newton_solver"]["relative_tolerance"] = 1e-6
solver.parameters["newton_solver"]["absolute_tolerance"] = 1e-6
solver.parameters["newton_solver"]["maximum_iterations"] = 10
u_file = XDMFFile("dir2" + "/velocity.xdmf")
p_file = XDMFFile("dir2" + "/pressure.xdmf")
s_file = XDMFFile("dir2" + "/saturation.xdmf")
while t < T:
t += float(dt)
U0.assign(U)
solver.solve()
u, p, s = U.split()
p_file.write(p, t)
s_file.write(s, t)
u_file.write(u, t)
UnitSquareMesh( ) code:
# Total mobility
def lmbdainv(s, mu_rel, no, nw):
return 1.0 / ((1.0 / mu_rel) * s ** nw + (1.0 - s) ** no)
# Fractional flow function
def F(s, mu_rel, no, nw):
return s ** nw / (s ** nw + mu_rel * (1.0 - s) ** no)
def tensor_jump(v, n):
return ufl.outer(v, n)("+") + ufl.outer(v, n)("-")
class PiecewiseConstant(UserExpression):
def __init__(self, values, markers, **kwargs):
self._values = values
self._markers = markers
super().__init__(**kwargs)
def eval_cell(self, values, x, cell):
values[0] = self._values[self._markers[cell.index]]
def value_shape(self):
return tuple()
class Obstacle(SubDomain):
def inside(self, x, on_boundary):
return between(x[1], (0.3, 0.7)) and between(x[0], (0.1, 0.5))
# Relative viscosity of water w.r.t. crude oil
Kinv = Constant(1 / (1e-6))
# Time step
dt = Constant(0.01)
mu = 1
mu_b = 1
phi = Constant(0.4)
pin = 6000
pout = 5000
sbar = Constant(1)
t = 0.0
T = 100 * float(dt)
mesh = UnitSquareMesh(90, 90)
n = FacetNormal(mesh)
order = 1
BDM = FiniteElement("BDM", mesh.ufl_cell(), order)
DG = FiniteElement("DG", mesh.ufl_cell(), order - 1)
Element = MixedElement([BDM, DG, DG])
mixed_space = FunctionSpace(mesh, Element)
# Function spaces and functions
V = TestFunction(mixed_space)
dU = TrialFunction(mixed_space)
U = Function(mixed_space)
U0 = Function(mixed_space)
v, q, r = split(V)
u, p, s = split(U)
u0, p0, s0 = split(U0)
s_mid = 0.5 * (s0 + s)
# ============= DEFINITION OF SPATIALLY-VARYING PARAMETERS =====================
marker_inner = 1
marker_outer = 0
no_outer = 2
nw_outer = 2
no_inner = 1
nw_inner = 1
obstacle = Obstacle()
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(marker_outer)
obstacle.mark(domains, marker_inner)
no = {marker_inner: no_inner, marker_outer: no_outer}
nw = {marker_inner: nw_inner, marker_outer: nw_outer}
VVV = FunctionSpace(mesh, "DG", 0)
noo = PiecewiseConstant(no, domains)
noo_proj = project(noo, VVV)
nww = PiecewiseConstant(nw, domains)
nww_proj = project(nww, VVV)
# =========== END DEFINITION OF SPATIALLY-VARYING PARAMETERS ===================
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
left = AutoSubDomain(lambda x: near(x[0], 0.0))
right = AutoSubDomain(lambda x: near(x[0], 1.0))
bottom = AutoSubDomain(lambda x: near(x[1], 0.0))
top = AutoSubDomain(lambda x: near(x[1], 1.0))
# Define boundary markers
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dx = Measure("dx", domain=mesh, subdomain_data=domains)
bc2 = DirichletBC(mixed_space.sub(0), Constant((0.0, 0.0)), boundaries, 2)
# # bc3 = DirichletBC(VQ.sub(0), Constant((0.0, 0.0)), boundaries, 3)
bc4 = DirichletBC(mixed_space.sub(0), Constant((0.0, 0.0)), boundaries, 4)
bcs = [bc2, bc4] # velocity BC
alpha = 35
h = CellDiameter(mesh)
h2 = ufl.Min(h("+"), h("-"))
stab = (
mu * (alpha / h2) * inner(tensor_jump(u, n), tensor_jump(v, n)) * dS
- mu * inner(avg(grad(u)), tensor_jump(v, n)) * dS
- mu * inner(avg(grad(v)), tensor_jump(u, n)) * dS
)
L1 = (
mu_b * inner(grad(u), grad(v)) * dx(1)
+ inner(v, lmbdainv(s_mid, mu, no_outer, nw_outer) * Kinv * u) * dx(0)
- div(v) * p * dx
+ pin * dot(v, n) * ds(1)
+ pout * dot(v, n) * ds(3)
+ stab
)
L2 = q * div(u) * dx(0) + q * div(u) * dx(1)
un = 0.5 * (inner(u0, n) + abs(inner(u0, n)))
un_h = 0.5 * (inner(u0, n) - abs(inner(u0, n)))
stabilisation = (
dt("+")
* inner(
jump(r),
un("+") * F(s_mid, mu, noo_proj, nww_proj)("+")
- un("-") * F(s_mid, mu, noo_proj, nww_proj)("-"),
)
* dS
)
L3 = (
phi * r * (s - s0) * dx(0)
+ r * (s - s0) * dx(1)
- dt * inner(grad(r), F(s_mid, mu, no_outer, nw_outer) * u) * dx(0)
- dt * inner(grad(r), F(s_mid, mu, no_outer, nw_outer) * u) * dx(1)
+ dt * r * F(s_mid, mu, no_outer, nw_outer) * un * ds
+ stabilisation
+ dt * r * un_h * sbar * ds(1)
)
# Total L
L = L1 + L2 + L3
J = derivative(L, U, dU)
problem = NonlinearVariationalProblem(L, U, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)
u_file = XDMFFile("dir2" + "/velocity.xdmf")
p_file = XDMFFile("dir2" + "/pressure.xdmf")
s_file = XDMFFile("dir2" + "/saturation.xdmf")
while t < T:
t += float(dt)
U0.assign(U)
solver.solve()
u, p, s = U.split()
p_file.write(p, t)
s_file.write(s, t)
u_file.write(u, t)
Convert from .geo to .xdmf
import os
from dolfin import *
import meshio
file_in = "holesuCT"
file_out = "holesuCT"
os.system("gmsh " + file_in + ".geo -2 -o " + file_out + ".msh")
os.system("meshio-convert " + file_out + ".msh " + file_out + ".xdmf")
msh = meshio.read("holesuCT.msh")
tri_data = msh.cell_data_dict["gmsh:physical"]["triangle"]
line_data = msh.cell_data_dict["gmsh:physical"]["line"]
meshio.write(
"mesh.xdmf",
meshio.Mesh(
points=msh.points[:, :2], cells={"triangle": msh.cells_dict["triangle"]}
),
)
meshio.write(
"boundaries.xdmf",
meshio.Mesh(
points=msh.points,
cells={"line": msh.cells_dict["line"]},
cell_data={"bnd_marker": [msh.cell_data_dict["gmsh:physical"]["line"]]},
),
)
meshio.write(
"domains.xdmf",
meshio.Mesh(
points=msh.points,
cells={"triangle": msh.cells_dict["triangle"]},
cell_data={"dom_marker": [msh.get_cell_data("gmsh:physical", "triangle")]},
),
)
