Issue saving results in parrallel run with dolfinx

dolfinx
1

Hi all,

Recently, I have been trying to move from FEniCS to FEniCSx and I am trying to reproduce some of my older codes in FEniCSx. Previously, I used to save data with

filename5 = "Polarization_Beam.xdmf"
hdf5_file = XDMFFile(MPI.comm_world, filename5)
hdf5_file.write_checkpoint(polariz, "F", XDMFFile.Encoding.HDF5

and that worked for both serial and parallel simulation. Now, in FEniCSx, I am using

xdmf = XDMFFile(mesh.comm, "POT_Beam.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(pot)

which works fine when I run the code on one core, but when I run that with

 mpirun -np 6 python ScratchtestFenicsx3.py

I got weird results. I tried to find a solution for that from reading other similar posts but I couldn’t find a solution.

I am using paraview 5.10.0 to visualize the results.

Thanks,

2

Please provide a minimal reproducible example.

3

Thank you so much @dokken. Below is the code that I am running. Sorry, if it is not a minimal. This is a mixed element formulation, and I tried to make it a minimal reproducible code, but I couldn’t do that since I have to reduce the number of space and change the BC and the weak form to simplify that . If I do so, then there would be no issue.

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista

from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
                         locate_dofs_topological, assemble, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary, create_box, CellType
from petsc4py.PETSc import ScalarType
import meshio
from dolfinx.plot import create_vtk_mesh
from ufl import (FacetNormal, FiniteElement, TensorElement, Identity, TestFunctions, TrialFunctions, VectorElement,
                 div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, MixedElement, split)

# ......................................... Parameters
qe = 1.6021 * 1e-19
non_local = 1e-20
E = 1 * 152e9
nu = 0.33
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
e0 = 8.854 * 1e-12
perm = 11.8 * e0

f11 = 1 * 1e-6
f12 = 1 * 1e-6
f44 = 1 * 1e-6

f11_1 = 1 * 1e-6
f11_3 = 1 * 1e-6
#
f12_1 = 1 * 1e-6
f12_3 = 1 * 1e-6
f44_1 = 1 * 1e-6
f44_3 = 1 * 1e-6

mu11 = 1e-6
mu12 = 1e-6
mu44 = 1e-6

p0 = 1 * 1e23

mo_ho = 450 * 1e-4
dif_ho = 12 * 1e-4
# ......................................... Reading mesh file
L = 1000e-9
H = 100e-9
W = H
h = H / 2

mesh = create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
                  [20, 6, 6], cell_type=CellType.hexahedron)

# ......................................... Define Loading

Force = 1e-1  # M / (0.5 * L)
Area = H ** 2
tr_force = Force / Area
tr_S1 = Constant(mesh, ScalarType((0, 0, 0.0)))
# tr_S2 = fe.Constant((0.0, -tr_force, 0.0))
tr_S2 = Constant(mesh, ScalarType((0, 0, -tr_force)))
n = FacetNormal(mesh)
# ......................................... Define Elements
v = VectorElement("Lagrange", mesh.ufl_cell(), 2)
s = TensorElement("Lagrange", mesh.ufl_cell(), 1)
lm = TensorElement("Lagrange", mesh.ufl_cell(), 1)
p = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Ho = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
# ......................................... Define mixed element
MIXEL = MixedElement([v, s, lm, p, Ho])
Space = FunctionSpace(mesh, MIXEL)
# ......................................... Define Trail and Test Functions
(u, gradu, lamu, p, ho) = TrialFunctions(Space)
(del_u, del_gradu, del_lamu, del_p, del_ho) = TestFunctions(Space)
# ....................................... Mark boundary regions
tol = 0.5 * 10e-9


def Left(x, tol=0.1 * 10e-9):
    return np.isclose(x[0], 0, atol=tol)


def MidZ(x, tol=0.5 * 10e-9):
    return np.isclose(x[2], H / 2, atol=tol)


# ....................................... ds
boundaries = [(1, lambda x: np.isclose(x[0], 0, atol=tol)),
              (2, lambda x: np.isclose(x[0], L, atol=tol)),
              (3, lambda x: np.isclose(x[2], H / 2, atol=tol))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, 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 = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

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

# ....................................... BCs
S0, S0_to_Space = Space.sub(0).collapse()
# Clamp BC for left side
clamp = Function(S0)
clamp.x.set(0.0)
clamp_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Left)
clamp_dofs = locate_dofs_topological((Space.sub(0), S0), mesh.topology.dim - 1, clamp_facets)
bc0 = dirichletbc(clamp, clamp_dofs, Space.sub(0))
# .........Ground BC
P0, P0_to_Space = Space.sub(3).collapse()

ground = Function(P0)
ground.x.set(0.0)
ground_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, MidZ)
ground_dofs = locate_dofs_topological((Space.sub(3), P0), mesh.topology.dim - 1, ground_facets)
bc1 = dirichletbc(ground, ground_dofs, Space.sub(3))
# .........Hole BC
H0, H0_to_Space = Space.sub(4).collapse()

hground = Function(H0)
hground.x.set(0.0)
hground_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, MidZ)
hground_dofs = locate_dofs_topological((Space.sub(3), P0), mesh.topology.dim - 1, hground_facets)
bc2 = dirichletbc(hground, hground_dofs, Space.sub(4))

# ........................
bc = [bc0, bc1, bc2]

# ......................................... Define tensor
delta = Identity(len(u))
i, j, k, l, m, n = indices(6)


# Strain Tensor
def StrainT(u):
    return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])


# Stress Tensor
def Elas(la, mu):
    return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l],
                     (i, j, k, l))  # Elasticity tensor


def StressT(u):
    return as_tensor((Elas(la, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])


# Relaxed Strain Tensor
def RStrain(u):
    return as_tensor((u[j].dx(i)), [i, j])


# Strain Gradient Tensor using relaxed strain
def RStrianGrad(gradu):
    return as_tensor((1. / 2. * (gradu[k, j].dx(i) + gradu[j, k].dx(i))), [i, j, k])


# Higher Order Stress
def SGElas(la, non_local, mu):
    return as_tensor((la * (non_local ** 2) * delta[i, j] * delta[l, m] * delta[k, n] +
                      2 * mu * (non_local ** 2) * delta[i, l] * delta[j, m] * delta[k, n]),
                     (i, j, k, l, m, n))  # Strain Gradient Elasticity Tensor


def HStress(gradu):
    return as_tensor(
        (SGElas(la, non_local, mu)[i, j, k, l, m, n] * RStrianGrad(gradu)[l, m, n]),
        (i, j, k))


# Flexoelectricity Tensor
fl1 = np.zeros((3, 3, 3, 3))
fl1[0, 0, 0, 0] = 1
fl1[1, 1, 1, 1] = 1
fl1[2, 2, 2, 2] = 1

flexo1 = Constant(mesh, PETSc.ScalarType(fl1))
flexo = as_tensor(((mu11 - mu12 - 2 * mu44) * flexo1[i, j, k, l] + mu12 * delta[i, j] * delta[
    k, l] + mu12 * delta[i, k] * delta[j, l] + mu12 * delta[i, l] * delta[j, k]), (i, j, k, l))


def HStressElec(p):
    return as_tensor((-1 * flexo[l, i, j, k] * (-1 * nabla_grad(p))[l]), (i, j, k))


def ElecDis(gradu):
    return as_tensor(perm * delta[i, j] * (-1 * nabla_grad(p))[j] + (flexo[i, j, k, l] * RStrianGrad(gradu)[j, k, l]),
                     i)


def Polarization(gradu):
    return as_tensor((flexo[i, j, k, l] * RStrianGrad(gradu)[j, k, l]), i)


def Jp(ho):
    return as_tensor(-1 * qe * p0 * mo_ho * nabla_grad(p)[i] - qe * dif_ho * nabla_grad(ho)[i], i)


# ......................................... Dimensionless Form
Fl = 1  # 10 ** (-6)
Fs = 1  # E
NONDIM = Constant(mesh, PETSc.ScalarType(1.0))
NONDIM2 = Constant(mesh, PETSc.ScalarType(1.0))
# ......................................... Define Variational Formulation
F1 = NONDIM * StressT(u)[j, k] * StrainT(del_u)[j, k] * dx - NONDIM * lamu[j, k] * del_u[k].dx(j) * dx
F2 = NONDIM2 * HStress(gradu)[i, j, k] * RStrianGrad(del_gradu)[i, j, k] * dx + NONDIM * lamu[j, k] * \
     del_gradu[
         j, k] * dx

F3 = HStressElec(p)[i, j, k] * RStrianGrad(del_gradu)[i, j, k] * dx - ElecDis(gradu)[i] * (
        -1 * nabla_grad(del_p)[i]) * dx

F4 = p.dx(i).dx(i) * ds(3)  # Electrod BC

# j=0 and k=0
C00 = NONDIM * gradu[0, 0] * del_lamu[0, 0] * dx - NONDIM * u[0].dx(0) * del_lamu[0, 0] * dx
# j=0 and k=1
C01 = NONDIM * gradu[0, 1] * del_lamu[0, 1] * dx - NONDIM * u[1].dx(0) * del_lamu[0, 1] * dx
# j=0 and k=2
C02 = NONDIM * gradu[0, 2] * del_lamu[0, 2] * dx - NONDIM * u[2].dx(0) * del_lamu[0, 2] * dx
# j=1 and k=0
C10 = NONDIM * gradu[1, 0] * del_lamu[1, 0] * dx - NONDIM * u[0].dx(1) * del_lamu[1, 0] * dx
# j=1 and k=1
C11 = NONDIM * gradu[1, 1] * del_lamu[1, 1] * dx - NONDIM * u[1].dx(1) * del_lamu[1, 1] * dx
# j=1 and k=2
C12 = NONDIM * gradu[1, 2] * del_lamu[1, 2] * dx - NONDIM * u[2].dx(1) * del_lamu[1, 2] * dx
# j=2 and k=0
C20 = NONDIM * gradu[2, 0] * del_lamu[2, 0] * dx - NONDIM * u[0].dx(2) * del_lamu[2, 0] * dx
# j=2 and k=1
C21 = NONDIM * gradu[2, 1] * del_lamu[2, 1] * dx - NONDIM * u[1].dx(2) * del_lamu[2, 1] * dx
# j=2 and k=2
C22 = NONDIM * gradu[2, 2] * del_lamu[2, 2] * dx - NONDIM * u[2].dx(2) * del_lamu[2, 2] * dx

FDD = -qe * ho * del_p * dx - (1 / qe) * (Jp(ho)[i] * del_ho.dx(i)) * dx

VarForm = F1 + F2 + F3 + C00 + C01 + C02 + C10 + C11 + C12 + C20 + C21 + C22 + FDD  # + F4

# ......................................... Preparing the form
a = VarForm
L = tr_S1[i] * del_u[i] * ds(1) + tr_S2[i] * del_u[i] * ds(2)
##########
Result = Function(Space)
problem = LinearProblem(a, L, bcs=bc,
                        petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
Result = problem.solve()
uf = Result.sub(0).collapse()
pot = Result.sub(3).collapse()
hole = Result.sub(4).collapse()
# ......................................... Export Results
xdmf = XDMFFile(mesh.comm, "POT_Beam.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(pot)
4

Thank you @dokken , Here is a minimal example with the issue I already posted. I have a simple mesh and I want to find the area of one surface on that. If I run it with one core or 6 core, the code will print different value.

from mpi4py import MPI

proc = MPI.COMM_WORLD.rank
from petsc4py import PETSc
import numpy as np
import pyvista

from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
                         locate_dofs_topological, assemble, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary, create_box, CellType
from petsc4py.PETSc import ScalarType
import meshio
from dolfinx.plot import create_vtk_mesh
from ufl import (FacetNormal, FiniteElement, TensorElement, Identity, TestFunctions, TrialFunctions, VectorElement,
                 div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, MixedElement, split)

# ......................................... Reading mesh file
L = 1000e-9
H = 100e-9
h = H / 2


mesh = create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, H, H])],
                  [40, 10, 20], cell_type=CellType.hexahedron)

# ....................................... Mark boundary regions
tol = 1.0 * 1e-11

def Left(x, tol=1.0 * 1e-11):
    return np.isclose(x[0], 0, atol=tol)


def MidZ(x, tol=1.0 * 1e-11):
    return np.isclose(x[2], H / 2, atol=tol)


# ....................................... ds
boundaries = [(1, lambda x: np.isclose(x[0], 0, atol=tol)),
              (2, lambda x: np.isclose(x[0], L, atol=tol)),
              (3, lambda x: np.isclose(x[2], H / 2, atol=tol))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, 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 = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

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

#...............................................................
one = Constant(mesh, ScalarType(1.))
f = form(one * ds(2))
area = assemble_scalar(f)
print(area)
5

Assemble scalar does no parallel communication, so you have to sum these up with a mpi allreduce, as shown at
https://jsdokken.com/dolfinx-tutorial/chapter2/heat_code.html?highlight=allreduce#verifying-the-numerical-solution