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)
• [x] 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

• [x] Black formatted
• [x] 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)

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 96.62162% with 40 lines in your changes missing coverage. Please review.

Project coverage is 97.69%. Comparing base (9eb3527) to head (7d0a6d4). Report is 264 commits behind head on fenicsx.

Files with missing lines	Patch %	Lines
src/festim/material.py	30.43%	16 Missing :warning:
src/festim/hydrogen_transport_problem.py	97.80%	12 Missing :warning:
src/festim/boundary_conditions/dirichlet_bc.py	95.45%	5 Missing :warning:
src/festim/mesh/mesh.py	91.37%	5 Missing :warning:
src/festim/exports/vtx.py	98.18%	1 Missing :warning:
src/festim/subdomain/interface.py	97.29%	1 Missing :warning:

Additional details and impacted files

@@             Coverage Diff             @@
##           fenicsx     #878      +/-   ##
===========================================
- Coverage    98.93%   97.69%   -1.24%     
===========================================
  Files           36       43       +7     
  Lines         1591     2042     +451     
===========================================
+ Hits          1574     1995     +421     
- Misses          17       47      +30     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[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

[jhdark]

Also, you should update the proposed user interface on the PR with the other classes that were introduced