Fix penalty term

[RemDelaporteMathurin]

Proposed changes

Closes #1026

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)

Further comments

If this is a relatively large or complex change, kick off the discussion by explaining why you chose the solution you did and what alternatives you considered, etc...

[codecov]

Codecov Report

❌ Patch coverage is 88.23529% with 4 lines in your changes missing coverage. Please review.
✅ Project coverage is 95.06%. Comparing base (0783ddc) to head (a2958b9).

Files with missing lines	Patch %	Lines
src/festim/hydrogen_transport_problem.py	88.23%	4 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           fenicsx    #1030      +/-   ##
===========================================
+ Coverage    94.95%   95.06%   +0.10%     
===========================================
  Files           44       44              
  Lines         3112     3120       +8     
===========================================
+ Hits          2955     2966      +11     
+ Misses         157      154       -3     

☔ View full report in Codecov by Sentry.
📢 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.

[RemDelaporteMathurin]

@jorgensd

I have an issue testing this, sometimes one of the domains will be full of zeros with penalty method.

This reproduces the issue (on the fix-penalty branch).

import festim as F
import ufl
import numpy as np
import dolfinx
from mpi4py import MPI


def generate_mesh(n=20):
    def bottom_boundary(x):
        return np.isclose(x[1], 0.0)

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

    def half(x):
        return x[1] <= 0.5 + 1e-14

    mesh = dolfinx.mesh.create_unit_square(
        MPI.COMM_WORLD, n, n, dolfinx.mesh.CellType.triangle
    )

    # Split domain in half and set an interface tag of 5
    gdim = mesh.geometry.dim
    tdim = mesh.topology.dim
    fdim = tdim - 1
    top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary)
    bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom_boundary)
    num_facets_local = (
        mesh.topology.index_map(fdim).size_local
        + mesh.topology.index_map(fdim).num_ghosts
    )
    facets = np.arange(num_facets_local, dtype=np.int32)
    values = np.full_like(facets, 0, dtype=np.int32)
    values[top_facets] = 1
    values[bottom_facets] = 2

    bottom_cells = dolfinx.mesh.locate_entities(mesh, tdim, half)
    num_cells_local = (
        mesh.topology.index_map(tdim).size_local
        + mesh.topology.index_map(tdim).num_ghosts
    )
    cells = np.full(num_cells_local, 4, dtype=np.int32)
    cells[bottom_cells] = 3
    ct = dolfinx.mesh.meshtags(
        mesh, tdim, np.arange(num_cells_local, dtype=np.int32), cells
    )
    all_b_facets = dolfinx.mesh.compute_incident_entities(
        mesh.topology, ct.find(3), tdim, fdim
    )
    all_t_facets = dolfinx.mesh.compute_incident_entities(
        mesh.topology, ct.find(4), tdim, fdim
    )
    interface = np.intersect1d(all_b_facets, all_t_facets)
    values[interface] = 5

    mt = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, facets, values)
    return mesh, mt, ct


K_S_top = 3.0
K_S_bot = 6.0
D_top = 2.0
D_bot = 2.0

mesh, mt, ct = generate_mesh(20)

my_model = F.HydrogenTransportProblemDiscontinuous()
my_model.mesh = F.Mesh(mesh)
my_model.volume_meshtags = ct
my_model.facet_meshtags = mt

material_top = F.Material(D_0=D_top, E_D=0, K_S_0=K_S_top, E_K_S=0)
material_bottom = F.Material(D_0=D_bot, E_D=0, K_S_0=K_S_bot, E_K_S=0)

material_top.solubility_law = "sievert"
material_bottom.solubility_law = "sievert"

top_domain = F.VolumeSubdomain(4, material=material_top)
bottom_domain = F.VolumeSubdomain(3, material=material_bottom)

top_surface = F.SurfaceSubdomain(id=1)
bottom_surface = F.SurfaceSubdomain(id=2)
my_model.subdomains = [
    bottom_domain,
    top_domain,
    top_surface,
    bottom_surface,
]

my_model.interfaces = [
    F.Interface(5, (bottom_domain, top_domain), penalty_term=10),
]
my_model.surface_to_volume = {
    top_surface: top_domain,
    bottom_surface: bottom_domain,
}

H = F.Species("H", mobile=True)

my_model.species = [H]

for species in my_model.species:
    species.subdomains = [bottom_domain, top_domain]

my_model.boundary_conditions = [
    F.FixedConcentrationBC(top_surface, value=0, species=H),
    F.FixedConcentrationBC(bottom_surface, value=1, species=H),
]

x = ufl.SpatialCoordinate(mesh)

my_model.sources = []

my_model.temperature = 500.0  # lambda x: 300 + 10 * x[1] + 100 * x[0]

my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
my_model.exports = [
    F.VTXSpeciesExport(f"u_{subdomain.id}.bp", field=H, subdomain=subdomain)
    for subdomain in my_model.volume_subdomains
]

my_model.initialise()
my_model.run()

With my_model.method_interface = "nitsche" I get the correct result:

Any idea?

[RemDelaporteMathurin]

Update:

• @jorgensd found out that having a non-zero initial guess helps!
• I'm now struggling to test/use this on real cases with real magnitudes for solubilities/concentrations... turns out the previous penalty method works better in this specific case