Interface condition diagnostic

[jhdark]

Description

Fix for issue #1180

Summary

This PR adds an additional export which evaluates the residual of the interface condition, which should confirm whether the penalty term used is high enough to enforce the condition properly.

Here is a small MWE which plots the max value of the residual against the penalty term used in both Sievert-Sieverts and Sieverts-Henry cases:

from mpi4py import MPI

import dolfinx
import matplotlib.pyplot as plt
import numpy as np

import festim as F


def sim(penalty, mat_r_law):
    my_model = F.HydrogenTransportProblemDiscontinuous()

    msh = dolfinx.mesh.create_unit_square(
        MPI.COMM_WORLD,
        10,
        10,
        cell_type=dolfinx.mesh.CellType.triangle,
        diagonal=dolfinx.mesh.DiagonalType.crossed,
    )

    my_mesh = F.Mesh(msh)
    my_model.mesh = my_mesh

    mat_L = F.Material(D_0=1e-03, E_D=0.0, K_S_0=1, E_K_S=0.0, solubility_law="sievert")
    mat_R = F.Material(D_0=1e-03, E_D=0.0, K_S_0=2, E_K_S=0.0, solubility_law=mat_r_law)
    vol_L = F.VolumeSubdomain(id=1, material=mat_L, locator=lambda x: x[0] <= 0.5)
    vol_R = F.VolumeSubdomain(id=2, material=mat_R, locator=lambda x: x[0] >= 0.5)
    left = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0.0))
    right = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 1.0))
    my_model.subdomains = [vol_L, vol_R, left, right]

    H = F.Species("H", subdomains=[vol_L, vol_R])
    my_model.species = [H]

    center_interface = F.Interface(
        id=5, subdomains=[vol_L, vol_R], penalty_term=penalty
    )
    my_model.interfaces = [center_interface]

    my_model.temperature = 100

    my_model.boundary_conditions = [
        F.FixedConcentrationBC(subdomain=left, value=1.0, species=H),
        F.FixedConcentrationBC(subdomain=right, value=0.0, species=H),
    ]

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

    my_model.exports = [
        F.VTXSpeciesExport(field=H, filename="mwe_L.bp", subdomain=vol_L),
        F.VTXSpeciesExport(field=H, filename="mwe_R.bp", subdomain=vol_R),
        F.VTXInterfaceResidualExport(
            field=H, filename="mwe_interface.bp", interface=center_interface
        ),
    ]

    my_model.initialise()
    my_model.run()

    max_residual = np.max(np.absolute(my_model.exports[2].function.x.array[:]))

    return max_residual


if __name__ == "__main__":
    penalty_terms = np.geomspace(1e-02, 1e05, num=10)
    mat_laws = ["sievert", "henry"]
    max_values_ss = []
    max_values_sh = []

    for penalty in penalty_terms:
        max_residual = sim(penalty=penalty, mat_r_law="sievert")
        print(
            f"Sievert-Sievert, Penalty term: {penalty:.2e}, Max interface residual: {max_residual:.2e}"
        )
        max_values_ss.append(max_residual)

    for penalty in penalty_terms:
        max_residual = sim(penalty=penalty, mat_r_law="henry")
        print(
            f"Sievert-Henry, Penalty term: {penalty:.2e}, Max interface residual: {max_residual:.2e}"
        )
        max_values_sh.append(max_residual)

    plt.figure()
    plt.plot(penalty_terms, max_values_ss, marker="o", label="Sievert-Sievert")
    plt.plot(penalty_terms, max_values_sh, marker="o", label="Sievert-Henry")
    plt.xscale("log")
    plt.yscale("log")
    plt.xlabel("Penalty term")
    plt.ylabel("Max interface residual")
    ax = plt.gca()
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    plt.legend()

    plt.show()

Related Issues

Motivation and Context

Type of Change

• 🐛 Bug fix (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 (no functional changes, no API changes)
• 📝 Documentation update
• ✅ Test update (adding missing tests or correcting existing tests)
• 🔧 Build/CI configuration change

Testing

• All existing tests pass locally (pytest)
• I have added new tests that prove my fix is effective or that my feature works

Code Quality Checklist

• My code follows the code style of this project (Ruff formatted: ruff format .)
• My code passes linting checks (ruff check .)
• I have performed a self-review of my own code
• I have commented my code, particularly in hard-to-understand areas

Documentation

• I have updated the documentation accordingly (if applicable)
• I have added docstrings to new functions/classes following the project conventions

Breaking Changes

Screenshots/Examples

Additional Notes

[codecov]

Codecov Report

❌ Patch coverage is 31.08108% with 51 lines in your changes missing coverage. Please review.
✅ Project coverage is 93.62%. Comparing base (5536c33) to head (a54ec6c).

Files with missing lines	Patch %	Lines
src/festim/exports/vtx.py	21.81%	43 Missing ⚠️
src/festim/hydrogen_transport_problem.py	12.50%	7 Missing ⚠️
src/festim/subdomain/interface.py	90.90%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1169      +/-   ##
==========================================
- Coverage   94.76%   93.62%   -1.15%     
==========================================
  Files          46       46              
  Lines        3534     3590      +56     
==========================================
+ Hits         3349     3361      +12     
- Misses        185      229      +44     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[jhdark]

@RemDelaporteMathurin So far, this only interpolates the values of the residual onto a function which lives on the interface, but maybe we could also add something additional to the logger, like obtaining the max absolute value from the function? Let me know what you think

[RemDelaporteMathurin]

i think using .x.array for this might be risky since it assumes they ordered wrt dofs in the exact same way, which might not be True.

the proper way to do it is to make a fem.Expression object with u_0/K_s - u_1/K_s and then interpolate it in self.function

[RemDelaporteMathurin]

• Logging: yes!
• We also need to be able to export not only the residual c/s - c/s but c/s itself