Differences between `dolfin` and `dolfinx`

bhaveshshrimali · July 28, 2024, 3:29pm

Hi everyone,

I was trying to solve a simple hyperelastic problem with prescribed facet displacements and was noticing discrepancies between the results from legacy dolfin and dolfinx.

I am running the following versions:

dolfin: version 2019.2.0.dev0
- commit 0d37780440ad11e15f8aeec1057001761fdd2ba7 if relevant,
dolfinx: version 0.8.0
- commit ef9b7c15e6bb063acba3a568678881bc2f6e3a20.

The physical problem is a unit cube with the left (x=0) face fixed and prescribed displacements in the x direction on the top/bottom (y=1, y=0) and front/back (z=1, z=0) faces. The prescribed displacement is u(x,y,z) = x \cdot u_0, where u_0 = 0.05.

Legacy dolfin code:

from dolfin import (
    UnitCubeMesh, NonlinearVariationalProblem, NonlinearVariationalSolver, MPI,
    parameters, set_log_level,
    FunctionSpace, Function, LocalSolver,
    grad, inner, det, Identity, VectorElement, FiniteElement, MixedElement, split,
    TrialFunction, TestFunction, derivative, Constant,
    dx, CompiledSubDomain, Expression,
    DirichletBC, XDMFFile
)

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["quadrature_degree"] = 5
parameters["form_compiler"]["representation"] = "uflacs"
parameters["linear_algebra_backend"] = "PETSc"
ffc_options = {"optimize": True, "eliminate_zeros": True, "precompute_basis_const": True,"precompute_ip_const": True}
set_log_level(10)
comm = MPI.comm_world

mesh = UnitCubeMesh(2,2,2)

Wu = VectorElement("CG", mesh.ufl_cell(), 2)
Wp = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(
    mesh, MixedElement([Wu, Wp])
)
u_p = Function(V)
v_q = TestFunction(V)
du_dp = TrialFunction(V)
_u, _p = split(u_p)
_F = Identity(3) + grad(_u)
I1 = inner(_F, _F)
J = det(_F)

mu = Constant(0.1)
energy = _p * (J - 1) + mu * (I1 - 3)
energy *= dx
residual = derivative(energy, u_p, v_q)
jacobian = derivative(residual, u_p, du_dp)

# Mark boundary subdomains
bottom_fc=CompiledSubDomain("near(x[1],0.0) && on_boundary")                                 # bottom face
top_fc=CompiledSubDomain("near(x[1],1) && on_boundary")                                      # top face
right_fc=CompiledSubDomain("near(x[0],1) && on_boundary")                                    # right face
left_fc=CompiledSubDomain("near(x[0],0.0) && on_boundary")                                   # left face
front_fc=CompiledSubDomain("near(x[2],1) && on_boundary")                                  # front face
back_fc=CompiledSubDomain("near(x[2],0.0) && on_boundary")                                   # back face

strtch = Constant(0.05)
applied_disp = Expression("x[0] * strtch", strtch=float(strtch.values()[0]), degree=1)
# applied_disp.compute_vertex_values()
bc_left = DirichletBC(V.sub(0), Constant((0., 0., 0.)), left_fc)
bc_top_x = DirichletBC(V.sub(0).sub(0), applied_disp, top_fc)
bc_bottom_x = DirichletBC(V.sub(0).sub(0), applied_disp, bottom_fc)
bc_back= DirichletBC(V.sub(0).sub(0), applied_disp, back_fc)
bc_front_x = DirichletBC(V.sub(0).sub(0), applied_disp, front_fc)

bcs = [
    bc_left, bc_top_x, bc_bottom_x, bc_back, bc_front_x
]
problem = NonlinearVariationalProblem(F=residual, u=u_p, bcs=bcs, J=jacobian, form_compiler_parameters=ffc_options)
solver = NonlinearVariationalSolver(problem)
solver.solve()

u, p = u_p.split(True)

with XDMFFile("dolfin_results.xdmf") as f:
    f.parameters["flush_output"] = True
    f.parameters["functions_share_mesh"]=True
    f.parameters["rewrite_function_mesh"] = False
    f.write(u, 0.)
    f.write(p, 0.)

gives plausible results for the displacement field (x component in the contour below):

whereas the corresponding dolfinx code

import numpy as np
from mpi4py import MPI
from dolfinx.fem import (
    Constant,
    Function,
    functionspace,
    dirichletbc,
    locate_dofs_topological,
    petsc,
)
from dolfinx.mesh import (
    CellType,
    locate_entities_boundary,
    create_box,
)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile
from ufl import (
    Measure,
    grad,
    inner,
    Identity,
    det,
    TrialFunction,
    TestFunction,
    derivative,
    split,
)
from dolfinx.log import set_log_level, LogLevel

set_log_level(LogLevel.INFO)
comm = MPI.COMM_WORLD
mesh = create_box(
    comm,
    [np.array([0.0, 0.0, 0.0]), np.array([1.0, 1.0, 1.0])],
    n=[2, 2, 2],
    cell_type=CellType.tetrahedron,
)
ord = 2
elemu = element("P", mesh.basix_cell(), ord, shape=(3,))
elem_u_project = element("P", mesh.basix_cell(), ord - 1, shape=(3,))
elemp = element("P", mesh.basix_cell(), ord - 1)
elem = mixed_element([elemu, elemp])

V = functionspace(mesh, elem)
Vu = functionspace(mesh, elem_u_project)
Vp = functionspace(mesh, elemp)
u_p = Function(V)
v_q = TestFunction(V)
du_dp = TrialFunction(V)
u_h = Function(Vu, name="u")

u, p = split(u_p)

mu = Constant(mesh, default_scalar_type(1.0))
lmbda = Constant(mesh, default_scalar_type(1.0e3))

F = Identity(len(u)) + grad(u)
J = det(F)
I1 = inner(F, F)

dx = Measure("dx", domain=mesh, metadata={"quadrature_degree": 5})
energy = p * (J - 1) + mu * (I1 - 3)
energy *= dx
residual_form = derivative(energy, u_p, v_q)
jacobian_form = derivative(residual_form, u_p, du_dp)


def left(x):
    return np.isclose(x[0], 0.0)


def right(x):
    return np.isclose(x[0], 1)


def front(x):
    return np.isclose(x[2], 1)


def back(x):
    return np.isclose(x[2], 0.0)


def bottom(x):
    return np.isclose(x[1], 0.0)


def top(x):
    return np.isclose(x[1], 1)


# facets
left_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, left)
right_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, right)
front_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, front)
back_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, back)
bottom_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, bottom)
top_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, top)

boundary_dofs_bottom_x = locate_dofs_topological(
    (V.sub(0).sub(0), V), mesh.topology.dim - 1, bottom_boundary_facets
)
boundary_dofs_top_x = locate_dofs_topological(
    (V.sub(0).sub(0), V), mesh.topology.dim - 1, top_boundary_facets
)
boundary_dofs_front_x = locate_dofs_topological(
    (V.sub(0).sub(0), V), mesh.topology.dim - 1, front_boundary_facets
)
boundary_dofs_back_x = locate_dofs_topological(
    (V.sub(0).sub(0), V), mesh.topology.dim - 1, back_boundary_facets
)
boundary_dofs_left = locate_dofs_topological(
    (V.sub(0), V), mesh.topology.dim - 1, left_boundary_facets
)

_V, _ = V.sub(0).collapse()
_V0, _ = V.sub(0).sub(0).collapse()
u_zero = Function(_V)
u_zero.x.array[:] = 0.0
fixed_bc = dirichletbc(u_zero, boundary_dofs_left, V.sub(0))

u_bdry = Function(_V0, name="u_bdry")
disp_x = Constant(mesh, default_scalar_type(0.05))


def disp_x_facets(x):
    return x[0] * disp_x.value


u_bdry.interpolate(disp_x_facets)

# u_bdry.sub(0).sub(0).interpolate(disp_x_facets)

top_facet_bc = dirichletbc(u_bdry, boundary_dofs_top_x, V.sub(0).sub(0))
bottom_facet_bc = dirichletbc(u_bdry, boundary_dofs_bottom_x, V.sub(0).sub(0))
front_facet_bc = dirichletbc(u_bdry, boundary_dofs_front_x, V.sub(0).sub(0))
back_facet_bc = dirichletbc(u_bdry, boundary_dofs_back_x, V.sub(0).sub(0))
right_bc = dirichletbc(
    Constant(mesh, default_scalar_type(0.01)),
    locate_dofs_topological(V.sub(0).sub(0), 2, right_boundary_facets),
    V.sub(0).sub(0),
)

bcs = [fixed_bc, front_facet_bc, back_facet_bc, bottom_facet_bc, top_facet_bc]

problem = petsc.NonlinearProblem(F=residual_form, u=u_p, bcs=bcs, J=jacobian_form)
solver = NewtonSolver(comm, problem=problem)

solver.solve(u=u_p)
u, p = u_p.split()

u_h.interpolate(u)


import os

trial_dir = os.path.join(os.getcwd(), "dolfinx_results")
os.makedirs(trial_dir, exist_ok=True)
with XDMFFile(comm, os.path.join(trial_dir, "dolfinx_edgeDisp.xdmf"), "w") as f:
    f.write_mesh(mesh)
    f.write_function(u_h)
    # f.write_function(p)

crashes with

[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

Stepping through a debugger, it seems to be happening in apply_lifting inside the solve call, which could suggest a problem with the bcs but its not obvious to me what that is.

Would anyone know what could be going wrong?
Thanks

bhaveshshrimali · July 28, 2024, 8:53pm

I think it is due to passing the wrong function spaces in locate_dofs_topological. Although it is still not intuitive to me as to what function spaces should be passed…

Here is a slightly modified working code:

import numpy as np
from mpi4py import MPI
from dolfinx.fem import (
    Constant,
    Function,
    functionspace,
    dirichletbc,
    locate_dofs_topological,
    petsc,
    locate_dofs_geometrical,
)
from dolfinx.mesh import (
    CellType,
    locate_entities_boundary,
    create_box,
)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile
from ufl import (
    Measure,
    grad,
    inner,
    Identity,
    det,
    TrialFunction,
    TestFunction,
    derivative,
    split,
)
from dolfinx.log import set_log_level, LogLevel
from petsc4py import PETSc

opts = PETSc.Options()


set_log_level(LogLevel.INFO)
comm = MPI.COMM_WORLD
Lx, Ly, Lz = 1., 1., 1.
mesh = create_box(
    comm,
    [np.array([0.0, 0.0, 0.0]), np.array([Lx, Ly, Lz])],
    n=[3, 10, 5],
    cell_type=CellType.tetrahedron,
)
ord = 2
elemu = element("P", mesh.basix_cell(), ord, shape=(3,))
elem_u_project = element("P", mesh.basix_cell(), ord - 1, shape=(3,))
elemp = element("P", mesh.basix_cell(), ord - 1)
elem = mixed_element([elemu, elemp])

V = functionspace(mesh, elem)
Vu = functionspace(mesh, elem_u_project)
Vp = functionspace(mesh, elemp)
u_p = Function(V)
v_q = TestFunction(V)
du_dp = TrialFunction(V)
u_h = Function(Vu, name="u")

u, p = split(u_p)

mu = Constant(mesh, default_scalar_type(1.0))
lmbda = Constant(mesh, default_scalar_type(1.0e3))

F = Identity(len(u)) + grad(u)
J = det(F)
I1 = inner(F, F)

dx = Measure("dx", domain=mesh, metadata={"quadrature_degree": 5})
energy = p * (J - 1) + mu * (I1 - 3)
energy *= dx
residual_form = derivative(energy, u_p, v_q)
jacobian_form = derivative(residual_form, u_p, du_dp)


def left(x):
    return np.isclose(x[0], 0.0)


def right(x):
    return np.isclose(x[0], Lx)


def front(x):
    return np.isclose(x[2], Lz)


def back(x):
    return np.isclose(x[2], 0.0)


def bottom(x):
    return np.isclose(x[1], 0.0)


def top(x):
    return np.isclose(x[1], Ly)


# facets
left_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, left)
right_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, right)
front_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, front)
back_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, back)
bottom_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, bottom)
top_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, top)

_V, _ = V.sub(0).collapse()
_V0, _ = V.sub(0).sub(0).collapse()


boundary_dofs_bottom_x = locate_dofs_topological(
    (V.sub(0).sub(0), _V0), mesh.topology.dim - 1, bottom_boundary_facets
)
boundary_dofs_top_x = locate_dofs_topological(
    (V.sub(0).sub(0), _V0), mesh.topology.dim - 1, top_boundary_facets
)
boundary_dofs_front_x = locate_dofs_topological(
    (V.sub(0).sub(0), _V0), mesh.topology.dim - 1, front_boundary_facets
)
boundary_dofs_back_x = locate_dofs_topological(
    (V.sub(0).sub(0), _V0), mesh.topology.dim - 1, back_boundary_facets
)

boundary_dofs_left = locate_dofs_topological(
    (V.sub(0), _V), mesh.topology.dim - 1, left_boundary_facets
)

u_zero = Function(_V)
u_zero.x.array[:] = 0.0
fixed_bc = dirichletbc(u_zero, boundary_dofs_left, V.sub(0))

u_bdry = Function(_V0, name="u_bdry")
disp_x = Constant(mesh, default_scalar_type(0.05))


def disp_x_facets(x):
    return x[0] * disp_x.value


u_bdry.interpolate(disp_x_facets)

# u_bdry.sub(0).sub(0).interpolate(disp_x_facets)

top_facet_bc = dirichletbc(u_bdry, boundary_dofs_top_x, V.sub(0).sub(0))
bottom_facet_bc = dirichletbc(u_bdry, boundary_dofs_bottom_x, V.sub(0).sub(0))
front_facet_bc = dirichletbc(u_bdry, boundary_dofs_front_x, V.sub(0).sub(0))
back_facet_bc = dirichletbc(u_bdry, boundary_dofs_back_x, V.sub(0).sub(0))
right_bc = dirichletbc(
    Constant(mesh, default_scalar_type(0.01)),
    locate_dofs_topological(V.sub(0).sub(0), 2, right_boundary_facets),
    V.sub(0).sub(0),
)

bcs = [fixed_bc, front_facet_bc, back_facet_bc, bottom_facet_bc, top_facet_bc]

problem = petsc.NonlinearProblem(F=residual_form, u=u_p, bcs=bcs, J=jacobian_form)
solver = NewtonSolver(comm, problem=problem)

solver.solve(u=u_p)
u, p = u_p.split()

u_h.interpolate(u)

import os

trial_dir = os.path.join(os.getcwd(), "trialDofs")
os.makedirs(trial_dir, exist_ok=True)
with XDMFFile(comm, os.path.join(trial_dir, "mixedFormulationDofs_Box.xdmf"), "w") as f:
    f.write_mesh(mesh)
    f.write_function(u_h)

A helpful place to look for me was the tests. Any docs elaborating what do the two function spaces in the V argument represent in locate_dofs_topological would be helpful!

github.com

FEniCS/dolfinx/blob/5d15d263f20a9d8ada051dd2ece0f35894d38d4a/python/test/unit/fem/test_bcs.py#L34-L67


      
          def test_locate_dofs_geometrical():
              """Test that locate_dofs_geometrical, when passed two function
              spaces, returns the correct degrees of freedom in each space"""
              mesh = create_unit_square(MPI.COMM_WORLD, 4, 8)
              p0, p1 = 1, 2
              P0 = element("Lagrange", mesh.basix_cell(), p0, dtype=default_real_type)
              P1 = element("Lagrange", mesh.basix_cell(), p1, dtype=default_real_type)
          
              W = functionspace(mesh, mixed_element([P0, P1]))
              V = W.sub(0).collapse()[0]
          
              with pytest.raises(RuntimeError):
                  locate_dofs_geometrical(W, lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
          
              dofs = locate_dofs_geometrical((W.sub(0), V), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
          
              # Collect dofs (global indices) from all processes
              dofs0_global = W.sub(0).dofmap.index_map.local_to_global(dofs[0])
              dofs1_global = V.dofmap.index_map.local_to_global(dofs[1])
              all_dofs0 = set(np.concatenate(MPI.COMM_WORLD.allgather(dofs0_global)))

This file has been truncated. show original

francesco-ballarin · July 29, 2024, 7:31am

dolfinx is an open source project. If you have ideas on how to improve the documentation feel free to open a pull request.

Topic		Replies	Views
Difference in DolfinX & Dolfin behaviour for elasticity dolfinx	0	502	April 14, 2023
Results differ between dolfin and dolfinx dolfinx	5	443	December 2, 2023
Unequal results for linear and nonlinear solver with dolfinx_mpc dolfinx	5	953	February 1, 2023
Hyper elasticity tutorial in 2D and 3D comparison General	10	374	January 30, 2024
Inconsistent results in parallel computation dolfinx	3	474	July 5, 2021

Differences between `dolfin` and `dolfinx`

Related topics