Interface concentration discontinuities (fenicsx)

[RemDelaporteMathurin]

Proposed changes

A massive thank you to @jorgensd who was instrumental in getting this functionality in FESTIM 🥇

This PR is a first (big) step towards interface discontinuities in the fenicsx branch (#719).

(requires #764)

This new approach to conservation of chemical potential relies on the mixed domain approach implemented in dolfinx 0.8.0 and based off a demo available here

I created a new problem class HTransportProblemDiscontinuous - which is not yet complete - to avoid breaking existing functionality. We will complete this class by integrating transient, post-processing, ICs, etc. Eventually this whole class will be merged with HydrogenTransportProblem.

In the current interface, users need to provide correspondence between surfaces and volumes as well as listing interfaces and which volumes they are connected to:

list_of_interfaces = {5: [left_domain, right_domain]}
my_model.subdomains = [left_domain, right_domain, left_surface, right_surface]

I believe this can be streamlined by FESTIM under the hood but I'll keep that for a future PR.

Proposed user interface:

import festim as F

import dolfinx
import numpy as np
import festim as F

my_model = F.HTransportProblemDiscontinuous()

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)

material_left = F.Material(D_0=2.0, E_D=0.1)
material_mid = F.Material(D_0=2.0, E_D=0.1)
material_right = F.Material(D_0=2.0, E_D=0.1)

material_left.K_S_0 = 2.0
material_left.E_K_S = 0
material_mid.K_S_0 = 4.0
material_mid.E_K_S = 0
material_right.K_S_0 = 6.0
material_right.E_K_S = 0

left_domain = F.VolumeSubdomain1D(3, borders=[0, interface_1], material=material_left)
middle_domain = F.VolumeSubdomain1D(4, borders=[interface_1, interface_2], material=material_mid)
right_domain = F.VolumeSubdomain1D(5, borders=[interface_2, 1], material=material_right)

left_surface = F.SurfaceSubdomain1D(id=1, x=vertices[0])
right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1])

# the ids here are arbitrary in 1D, you can put anything as long as it's not the same as the surfaces
my_model.interfaces = {6: [left_domain, middle_domain], 7: [middle_domain, right_domain]}
my_model.subdomains = [left_domain, middle_domain, right_domain, left_surface, right_surface]
my_model.surface_to_volume = {right_surface: right_domain, left_surface: left_domain}

H = F.Species("H", mobile=True)
trapped_H = F.Species("H_trapped", mobile=False)
empty_trap = F.ImplicitSpecies(n=0.5, others=[trapped_H])

my_model.species = [H, trapped_H]

for species in [H, trapped_H]:
    species.subdomains = [left_domain, middle_domain, right_domain]


my_model.reactions = [
    F.Reaction(
        reactant=[H, empty_trap],
        product=[trapped_H],
        k_0=2,
        E_k=0,
        p_0=0.1,
        E_p=0,
        volume=domain,
    )
    for domain in [left_domain, middle_domain, right_domain]
]

my_model.boundary_conditions = [
    F.DirichletBC(left_surface, value=0.05, species=H),
    F.DirichletBC(right_surface, value=0.2, species=H),
]


my_model.temperature = lambda x: 300 + 100 * x[0]

my_model.settings = F.Settings(atol=None, rtol=None, transient=False)


my_model.initialise()
my_model.run()


for subdomain in my_model.volume_subdomains:
    u_sub_0 = subdomain.u.sub(0).collapse()
    u_sub_0.name = "u_sub_0"

    u_sub_1 = subdomain.u.sub(1).collapse()
    u_sub_1.name = "u_sub_1"
    bp = dolfinx.io.VTXWriter(
        my_model.mesh.mesh.comm, f"u_{subdomain.id}.bp", [u_sub_0, u_sub_1], engine="BP4"
    )
    bp.write(0)
    bp.close()

Types of changes

What types of changes does your code introduce to FESTIM?

• Bugfix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected)
• Code refactoring
• Documentation Update (if none of the other choices apply)
• New tests

Checklist

• Black formatted
• Unit tests pass locally with my changes
• I have added tests that prove my fix is effective or that my feature works
• I have added necessary documentation (if appropriate)

[RemDelaporteMathurin]

We have noticed crashed in 1D simulations with coarse meshes. This is a fix for it (needs testing and review) and the working example is here (works with this PR).

[RemDelaporteMathurin]

We had a discussion with @jorgensd @jhdark and @JonathanGSDUFOUR and noticed weird things happening with the performance so we need to investigate. Otherwise I'll carry on writing tests for this