Running `CR` in parallel: Error detected in HDF5

[bhaveshshrimali]

I am trying to run a hyperelasticity problem with periodic boundary conditions using the Crouziex-Raviart elements. The code seems to run fine when run in serial, namely

python3 <codename>.py

But the same when run in parallel with mpi, namely

mpirun -np 2 python3 <codename>.py

gives

HDF5-DIAG: Error detected in HDF5 (1.10.0-patch1) MPI-process 1:
  #000: ../../../src/H5Dio.c line 263 in H5Dwrite(): file selection+offset not within extent
    major: Dataspace
    minor: Out of range

I am using dolfin-2019.1.0 on a singularity image. Any ideas what could be going wrong?

[nate]

Please provide a minimal working example so we can debug your problem.

[bhaveshshrimali]

Thanks for your reply Nate! So I prepared a minimal code that still fails in parallel, albeit with a different error. For reference I follow a similar approach as mentioned in COMET for creating a periodic map between the dofs, and the hyperelasticity example in demos.

from dolfin import *
import numpy as np
import os
comm = MPI.comm_world

set_log_level(20)
workdir = os.path.join(os.getcwd(), "DiscourseTrial", "CR")
os.makedirs(workdir, exist_ok=True)
parameters["ghost_mode"] = "shared_facet"
parameters["form_compiler"]["quadrature_degree"] = 5
parameters["form_compiler"]["representation"]="uflacs"

ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

snes_solver_parameters = {"nonlinear_solver": "snes",
                          "snes_solver": {"linear_solver": "lu",   #lu or gmres or cg or bicgstab 'preconditioner: ilu, amg, jacobi'
                                        "preconditioner": "petsc_amg",						  
                                        "maximum_iterations": 50,
                                        "report": True,
                                        "line_search": 'basic',
                                        "error_on_nonconvergence": False,
                                        "relative_tolerance":1e-7,
                                        "absolute_tolerance":1e-8}}

msh = UnitCubeMesh(6, 6, 6) # tried with a periodic mesh 
class PeriodicBoundary(SubDomain):
    def __init__(self, vertices, tolerance=DOLFIN_EPS):
        """ vertices stores the coordinates of the 4 unit cell corners"""
        SubDomain.__init__(self, tolerance)
        self.tol = tolerance
        self.vv = vertices
        self.a1 = self.vv[1,:]-self.vv[0,:] # first vector generating periodicity
        self.a2 = self.vv[2,:]-self.vv[0,:] # second vector generating periodicity 
        self.a3 = self.vv[3,:]-self.vv[0,:] # third vector generating periodicity

        
    def inside(self, x, on_boundary):
        """
            return True if on left, bottom, or back faces
            and not on one of the top, front or right faces 
            or associate edges (and vertices) as defined below 
        """
        # faces
        left = near(x[0],self.vv[0,0]) 
        bottom = near(x[1],self.vv[0,1]) 
        back = near(x[2],self.vv[0,2])
        right = near(x[0],self.vv[1,0])
        top = near(x[1],self.vv[2,1])
        front = near(x[2],self.vv[3,2])

        # line-segments (bottom 4; top 4; vertical 4)
        bottom_front = bottom and front 
        bottom_right = bottom and right 
        
        top_left = top and left 
        top_back = top and back 

        left_front = left and front 
        right_back = right and back 

        return bool((left or back or bottom) and 
                    (not( (top_left) or (left_front) or (top_back) or (right_back) or (bottom_right) or (bottom_front))) and on_boundary)

        # return bool( bool(left and not((top_left) or (left_front))) or bool(back and not((top_back) or (right_back))) or 
        #         bool(bottom and not((bottom_right) or (bottom_front))) and on_boundary )

    def map(self, x, y):
        """ Mapping the right boundary to left and top to bottom"""
        
        # faces
        right = near(x[0],self.vv[1,0])
        top = near(x[1],self.vv[2,1])
        front = near(x[2],self.vv[3,2])


        # line-segments 
        top_right = top and right 
        top_front = top and front 
        right_front = right and front 
        point_6 = right and front and top 

        if point_6:
            y[0] = x[0] - (self.a1[0] + self.a2[0] + self.a3[0])
            y[1] = x[1] - (self.a1[1] + self.a2[1] + self.a3[1])
            y[2] = x[2] - (self.a1[2] + self.a2[2] + self.a3[2])
        elif top_right:
            y[0] = x[0] - (self.a1[0] + self.a2[0])
            y[1] = x[1] - (self.a1[1] + self.a2[1])
            y[2] = x[2] - (self.a1[2] + self.a2[2])
        elif top_front:
            y[0] = x[0] - (self.a2[0] + self.a3[0])
            y[1] = x[1] - (self.a2[1] + self.a3[1])
            y[2] = x[2] - (self.a2[2] + self.a3[2])
        elif right_front: 
            y[0] = x[0] - (self.a1[0] + self.a3[0])
            y[1] = x[1] - (self.a1[1] + self.a3[1])
            y[2] = x[2] - (self.a1[2] + self.a3[2])
        elif right:
            y[0] = x[0] - (self.a1[0])
            y[1] = x[1] - (self.a1[1])
            y[2] = x[2] - (self.a1[2])
        elif front:
            y[0] = x[0] - (self.a3[0])
            y[1] = x[1] - (self.a3[1])
            y[2] = x[2] - (self.a3[2])
        elif top:
            y[0] = x[0] - (self.a2[0])
            y[1] = x[1] - (self.a2[1])
            y[2] = x[2] - (self.a2[2])
        else: 
            y[0] = -1. 
            y[1] = -1. 
            y[2] = -1.


def freeEnergyNG(u, params):
    I = Identity(3)             # Identity tensor
    F = I + grad(u)             # Deformation gradient
    C = F.T*F                   # Right Cauchy-Green tensor
    J = det(F)
    I1 = tr(C)

    mu1, mu2, alph1, alph2, mu_pr = params
    psiEq = 3**(1-alph1)/(2.*alph1) * mu1 * (I1**alph1 - 3**alph1) + \
            3**(1-alph2)/(2.*alph2) * mu2 * (I1**alph2 - 3**alph2) - \
                (mu1 + mu2) * ln(J) + mu_pr/2*(J-1)**2
    
    # psiEq = (mu1+mu2)/2.*(I1-3.) - (mu1+mu2)*ln(J) + mu_pr/2.*(J-1)**2 
    
    return psiEq

def assignLoading(i, fVal):
    Gamm_Voight = np.zeros(9)
    Gamm_Voight[0] = 1. * fVal
    Gamm_Voight[4] = Gamm_Voight[0]**(-0.5)
    Gamm_Voight[8] = Gamm_Voight[0]**(-0.5)

    return np.array([[Gamm_Voight[0]-1,    Gamm_Voight[1],       Gamm_Voight[2]],
                     [Gamm_Voight[3],      Gamm_Voight[4]-1 ,    Gamm_Voight[5]],
                     [Gamm_Voight[6],      Gamm_Voight[7],       Gamm_Voight[8]-1 ]])

deg = 1
Lx = Ly = Lz = 1
vertices = np.array([[0., 0,   0],   # 0: Origin
                     [Lx, 0,   0],   # 1: Right
                     [0,  Ly,  0],   # 2: Top
                     [0,  0,  Lz]])  # 3

kap_by_mu = Constant(1000)
rad = (0.2/(4./3*np.pi))**(1./3)
origin = CompiledSubDomain("near(x[0],0) && near(x[1],0) && near(x[2],0) && on_boundary")
rightPoint = CompiledSubDomain("near(x[0], Lx) && near(x[1],0) && near(x[2],0) && on_boundary", Lx = Lx)
topPoint = CompiledSubDomain("near(x[0],0) && near(x[1], Ly) && near(x[2],0) && on_boundary",   Ly = Ly)
frontPoint = CompiledSubDomain("near(x[0],0) && near(x[1],0) && near(x[2],Lz) && on_boundary",  Lz = Lz)

Wu = VectorElement("CR", msh.ufl_cell(), deg) # for displacement
Wlam = VectorElement("R", msh.ufl_cell(), 0)
V = FunctionSpace(msh, MixedElement([Wu, Wlam]), constrained_domain=PeriodicBoundary(vertices))
wper = Function(V)
uper, lamper = split(wper)

wper_  = TrialFunction(V)
uper_, lamper_ = split(wper)
dw_per = TestFunction(V)
du_per, dlam_per = split(dw_per)

Vplot = VectorFunctionSpace(msh, "CG", deg)
VQeplot = TensorFunctionSpace(msh, "DG", deg-1)  # at least one order lower

uplot = Function(Vplot, name="disp")

# *********************Non Gaussian Parameters ***********************
mu1 = Constant(1)
mu2 = Constant(1)
mu_pr = Constant(kap_by_mu*2)
alph1 = Constant(1) 
alph2 = Constant(1)
facMul = Constant(1000)
matParamsMatrix = [mu1, mu2, alph1, alph2, mu_pr]
matParamsParts = [mu1*facMul, mu2*facMul, alph1, alph2, mu_pr*facMul]
matParams = np.vstack((matParamsMatrix, matParamsParts))
qvals = Constant(20*(mu1 + mu2))
h = FacetArea(msh)
hAvg = 1./2*(h('+') + h('-'))

Hbar = Constant( ( (0,0,0), (0,0,0), (0,0,0) ))
x = SpatialCoordinate(msh)
bcOrigin = DirichletBC(V.sub(0), Constant((0., 0., 0.)), origin, method="pointwise")
bcAll = [bcOrigin]

amuv = derivative(freeEnergyNG(uper + dot(Hbar, x), matParams[0]), wper, dw_per) * dx
# amuv += inner(lamper_, du_per) * dx + inner(dlam_per, uper_) * dx
amuv += qvals/hAvg*inner(jump(uper), jump(du_per))*dS   # penalize the jump across inner facets
Jac = derivative(amuv, wper, wper_)

prblm  = NonlinearVariationalProblem(amuv, wper, bcAll, J=Jac)
solver = NonlinearVariationalSolver(prblm)
solver.parameters.update(snes_solver_parameters)

wfil = XDMFFile(comm, os.path.join(workdir, 'dispTrialSphereTension.xdmf'))
wfil.parameters["flush_output"] = True
wfil.parameters["functions_share_mesh"] = True
wfil.parameters["rewrite_function_mesh"] = False


Hbar.assign(Constant(assignLoading(0, 1.5)))
solver.solve()
uper, lamper = wper.split(True)
uplot.assign(project(uper + dot(Hbar, x) , Vplot))
wfil.write(uplot, 0)

gives the following error when run in parallel

RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to access vector of degrees of freedom.
*** Reason:  Cannot access a non-const vector from a subfunction.
*** Where:   This error was encountered inside Function.cpp.
*** Process: 1
***
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Edit: If requred I can furnish the entire code (with the mesh) that should reproduce the exact error that I posted in the first place. However I feel that the issue is with using “CR” elements on a periodic domain, and running the code in parallel

[bhaveshshrimali]

Any ideas what could be going wrong @nate?

[nate]

If you can produce a minimal (by which we mean as few lines as possible) example which reproduces your problem, you will be more likely to get people to help. See this, for example.

[bhaveshshrimali]

@nate

If you can produce a minimal (by which we mean as few lines as possible) example which reproduces your problem, you will be more likely to get people to help. See this , for example.

Alright, here you go (this is linear elasticity):

from dolfin import *
import numpy as np

comm = MPI.comm_world 
comm_rank = MPI.rank(comm)
set_log_level(30)

parameters["ghost_mode"] = "shared_facet"

mesh = UnitSquareMesh(20,20)

def corner(x):
    return ((near(x[0],0.)) and (near(x[1],0.)))

class PeriodicBoundary(SubDomain):
    def inside(self, x, on_boundary):
        left=bool(near(x[0],0.))
        bottom=bool(near(x[1],0.))
        right=bool(near(x[0],1.))
        top=bool(near(x[1],1.))                
        return bool((left or bottom) and (not ((left and top) or (bottom and right))) and on_boundary)

    def map(self, x, y):
        if near(x[0], 1.) and near(x[1], 1.):
            y[0] = x[0] - 1.
            y[1] = x[1] - 1.
        elif near(x[0], 1.):
            y[0] = x[0] - 1.
            y[1] = x[1] 
        elif near(x[1], 1.):
            y[0] = x[0] 
            y[1] = x[1] - 1.
        else:
            y[0] = -1.
            y[1] = -1.

V = VectorFunctionSpace(mesh, "CR", 1, constrained_domain=PeriodicBoundary())
u = Function(V)
v = TrialFunction(V)
dv = TestFunction(V)

mu=1
lmbda=5

f = Expression(("4*(l + 2*m)*pi*pi*cos(2*pi*x[0])","4*(l + 2*m)*pi*pi*sin(2*pi*x[1])"), l=lmbda, m=mu, degree=5)

bc1 = DirichletBC(V, Expression(("1","0"),degree=1), corner, method="pointwise")

def eps(v):
    return sym(grad(v))
    
F = (lmbda*tr(eps(v))*tr(eps(dv))+2*mu*inner(eps(v),eps(dv)))*dx - inner(f, dv)*dx 

h=FacetArea(mesh)         
h_avg = (h('+') + h('-'))/2
F += (1.0/h_avg)*dot(jump(v),jump(dv))*dS 

solve(lhs(F)==rhs(F), u, bc1)

[nate]

Looks like the problem is related to the periodic boundary condition. I’m unfamiliar with CR1 elements. Are the DoFs located at the edge or facet midpoints?

I think the periodic boundary conditions in old-dolfin work by geometrically finding the matching DoF coordinates and tying the master/slave combinations. Perhaps this isn’t being done correctly due to the way CR1 elements work.

[bhaveshshrimali]

Thanks for your reply! . The dofs are located at Facet-midpoints (c.f. page 98 of the FEniCS book).

@bleyerj do you have any previous experience using CR for incompressible elasticity on periodic domains?

[bhaveshshrimali]

@nate

Perhaps this isn’t being done correctly due to the way CR1 elements work.

I am not sure if the issue is with CR. But in any case, why would this work in serial and not parallel?

[nate]

Probably due to the "shared_facet" mode. The periodic BC routine searches geometrically rather than topologically. So a DoF coordinate may exist on more than one process since extra cell and facet information is available on the process boundaries. Turn off "shared_facet" and your interior facet integrals, and you’ll see that the code works in parallel.

[bhaveshshrimali]

@nate

Turn off "shared_facet" and your interior facet integrals, and you’ll see that the code works in parallel.

I expect that to work, but not necessarily give a correct result. The interior facet integrals penalize the jump of the displacement field across the interior faces (which makes sense physically). I wonder what would happen if I remove that altogether.

[nate]

Naturally your FE approximation will not be correct if you turn off interior facet integrals. I’m just identifying the problem, which is an implementation issue. A fix would be to debug the constrained DofMap construction in FEniCS for the parallel case when you have ghost cells and facets.

[bleyerj]

NO sorry never tried it with periodic bcs nor in parallel

[bhaveshshrimali]

@nate

A fix would be to debug the constrained DofMap construction in FEniCS for the parallel case when you have ghost cells and facets.

I see. Thanks for you reply!