Multi-material issue on `fenicsx`

[RemDelaporteMathurin]

Run on the fenicsx branch

import dolfinx
import festim as F

dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)

volume_file = "mesh_0.05.xdmf"
facet_file = "mesh_0.05_facet.xdmf"

my_model = F.HydrogenTransportProblem()
my_model.mesh = F.MeshFromXDMF(volume_file=volume_file, facet_file=facet_file)

mat = F.Material(D_0=1, E_D=0.1)

vol1 = F.VolumeSubdomain(id=1, material=mat)
vol2 = F.VolumeSubdomain(id=2, material=mat)
vol3 = F.VolumeSubdomain(id=3, material=mat)

surface1 = F.SurfaceSubdomain(id=4)
surface2 = F.SurfaceSubdomain(id=5)

my_model.subdomains = [vol1, vol2, vol3, surface1, surface2]

mobile = F.Species(name="H", mobile=True)
trapped = F.Species(name="H_trapped", mobile=False)
empty = F.ImplicitSpecies(n=0.5, others=[trapped])

my_model.species = [mobile, trapped]

my_model.reactions = [
    F.Reaction(
        reactant=[mobile, empty],
        product=[trapped],
        k_0=1,
        E_k=mat.E_D,
        p_0=0.1,
        E_p=0.87,
        volume=vol,
    )
    for vol in my_model.volume_subdomains
]

my_model.temperature = 600

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=surface1, species=mobile, value=1),
    F.FixedConcentrationBC(subdomain=surface2, species=mobile, value=0),
]

my_model.settings = F.Settings(
    atol=1e-6,
    rtol=1e-6,
    transient=False,
)

my_model.initialise()

my_model.run()

[2024-10-29 13:58:28.903] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2024-10-29 13:58:28.903] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2024-10-29 13:58:28.919] [info] XDMF build map
[2024-10-29 13:58:29.009] [info] XDMF create meshtags
[2024-10-29 13:58:29.009] [info] Building MeshTags object from tagged entities (defined by vertices).
[2024-10-29 13:58:29.009] [info] Build list of mesh entity indices from the entity vertices.
[2024-10-29 13:58:30.733] [info] Column ghost size increased from 0 to 0
[2024-10-29 13:58:30.833] [info] Newton iteration 0: r (abs) = 22.65029703975096 (tol = 1e-06), r (rel) = inf (tol = 1e-06)
[2024-10-29 13:58:31.178] [info] PETSc Krylov solver starting to solve system.
[2024-10-29 13:58:53.497] [info] Newton iteration 1: r (abs) = -nan (tol = 1e-06), r (rel) = -nan (tol = 1e-06)
[2024-10-29 13:58:53.834] [info] PETSc Krylov solver starting to solve system.

[RemDelaporteMathurin]

mesh_0.05.zip

This is the mesh I use

[RemDelaporteMathurin]

If we comment out the Reaction objects then everything runs fine

[RemDelaporteMathurin]

The 1D equivalent works fine....

import numpy as np

import festim as F
import dolfinx

dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)

my_model = F.HydrogenTransportProblem()
interface_1 = 0.2
interface_2 = 0.8

vertices = np.concatenate(
    [
        np.linspace(0, interface_1, num=100),
        np.linspace(interface_1, interface_2, num=100),
        np.linspace(interface_2, 1, num=100),
    ]
)

my_model.mesh = F.Mesh1D(vertices)

mat = F.Material(D_0=1, E_D=0)

vol1 = F.VolumeSubdomain1D(id=1, material=mat, borders=[0, interface_1])
vol2 = F.VolumeSubdomain1D(id=2, material=mat, borders=[interface_1, interface_2])
vol3 = F.VolumeSubdomain1D(id=3, material=mat, borders=[interface_2, 1])
surface1 = F.SurfaceSubdomain1D(id=4, x=vertices[0])
surface2 = F.SurfaceSubdomain1D(id=5, x=vertices[-1])

my_model.subdomains = [vol1, vol2, vol3, surface1, surface2]

mobile = F.Species(name="H", mobile=True)
trapped_1 = F.Species(name="H_trapped_W1", mobile=False)
empty_1a = F.ImplicitSpecies(n=0.5, others=[trapped_1])

my_model.species = [mobile, trapped_1]

my_model.reactions = [
    F.Reaction(
        reactant=[mobile, empty_1a],
        product=[trapped_1],
        k_0=1,
        E_k=mat.E_D,
        p_0=0.1,
        E_p=0.87,
        volume=vol,
    )
    for vol in my_model.volume_subdomains
]

my_model.temperature = 600

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=surface1, species=mobile, value=2),
    F.FixedConcentrationBC(subdomain=surface2, species=mobile, value=0),
]

my_model.settings = F.Settings(
    atol=1e-6,
    rtol=1e-6,
    transient=False,
)

my_model.initialise()

my_model.run()

[RemDelaporteMathurin]

I tried to make a 2D MWE with a dolfinx mesh but couldn't reproduce the issue with this:

import dolfinx
import festim as F
import numpy as np
import mpi4py.MPI as MPI

def make_mesh(N=10):

    mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N)

    def top_half(x):
        return x[1] > 0.5 - 1e-10

    def bottom_half(x):
        return x[1] < 0.5 + 1e-10

    vdim = mesh.topology.dim
    fdim = mesh.topology.dim - 1

    num_cells = mesh.topology.index_map(vdim).size_local
    mesh_cell_indices = np.arange(num_cells, dtype=np.int32)
    tags_volumes = np.full(num_cells, 0, dtype=np.int32)

    entities_top = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, top_half)
    entities_bot = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, bottom_half)

    tags_volumes[entities_bot] = 1
    tags_volumes[entities_top] = 2
    ct = dolfinx.mesh.meshtags(
        mesh, dim=mesh.topology.dim, entities=mesh_cell_indices, values=tags_volumes
    )

    # find all facets in domain and mark them as 0

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

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

    entities_top = dolfinx.mesh.locate_entities(mesh, fdim, top_surface)
    entities_bot = dolfinx.mesh.locate_entities(mesh, fdim, bottom_surface)

    # concatenate entities in one array and mark them as 1 and 2
    entities = np.concatenate([entities_bot, entities_top])
    tags_facets = np.concatenate(
        [np.full_like(entities_bot, 4), np.full_like(entities_top, 5)]
    )

    ft = dolfinx.mesh.meshtags(mesh, dim=fdim, entities=entities, values=tags_facets)

    return mesh, ct, ft

dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)

my_model = F.HydrogenTransportProblem()

mesh, ct, ft = make_mesh(N=10)

# with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "w") as xdmf:
#     xdmf.write_mesh(mesh)
#     xdmf.write_meshtags(ct, mesh.geometry)

# with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "ft.xdmf", "w") as xdmf:
#     xdmf.write_mesh(mesh)
#     xdmf.write_meshtags(ft, mesh.geometry)

my_model.mesh = F.Mesh(mesh=mesh)
my_model.volume_meshtags = ct
my_model.facet_meshtags = ft

mat = F.Material(D_0=1, E_D=0.1)

vol1 = F.VolumeSubdomain(id=1, material=mat)
vol2 = F.VolumeSubdomain(id=2, material=mat)

surface1 = F.SurfaceSubdomain(id=4)
surface2 = F.SurfaceSubdomain(id=5)

my_model.subdomains = [vol1, vol2, surface1, surface2]

mobile = F.Species(name="H", mobile=True)
trapped = F.Species(name="H_trapped", mobile=False)
empty = F.ImplicitSpecies(n=0.5, others=[trapped])

my_model.species = [mobile, trapped]

my_model.reactions = [
    F.Reaction(
        reactant=[mobile, empty],
        product=[trapped],
        k_0=1,
        E_k=mat.E_D,
        p_0=0.1,
        E_p=0.87,
        volume=vol,
    )
    for vol in my_model.volume_subdomains
]

my_model.temperature = 600

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=surface1, species=mobile, value=1),
    F.FixedConcentrationBC(subdomain=surface2, species=mobile, value=0),
]

my_model.settings = F.Settings(
    atol=1e-6,
    rtol=1e-6,
    transient=False,
)

my_model.initialise()

my_model.run()

[jorgensd]

Setting the following in ProblemBase.create_solver resolves it:

        ksp = self.solver.krylov_solver
        ksp.setType("preonly")
        ksp.getPC().setType("lu")
        ksp.getPC().setFactorSolverType("mumps")
        ksp.setErrorIfNotConverged(True)