Merge pull request #878 from RemDelaporteMathurin/discontinuity-generic

Interface concentration discontinuities (fenicsx)

.github/dependabot.yml

@@ -0,0 +1,15 @@
+# To get started with Dependabot version updates, you'll need to specify which
+# package ecosystems to update and where the package manifests are located.
+# Please see the documentation for all configuration options:
+# https://docs.github.com/code-security/dependabot/dependabot-version-updates/configuration-options-for-the-dependabot.yml-file
+
+version: 2
+updates:
+  - package-ecosystem: "github-actions" # See documentation for possible values
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"
+  - package-ecosystem: "pip" # See documentation for possible values
+    directory: "./"
+    schedule:
+      interval: "weekly"

.github/workflows/benchmark.yml

@@ -4,16 +4,16 @@ on: [pull_request, push]
 jobs:
   benchmark:
     runs-on: ubuntu-latest
-    container: dolfinx/dolfinx:v0.7.3
+    container: dolfinx/dolfinx:v0.9.0
 
     steps:
-    - name: Checkout code
-      uses: actions/checkout@v2
+      - name: Checkout code
+        uses: actions/checkout@v2
 
-    - name: Install local package and dependencies
-      run: |
-        pip install .[test]
+      - name: Install local package and dependencies
+        run: |
+          pip install .[test]
 
-    - name: Run benchmark
-      run: |
-        python3 test/benchmark.py
+      - name: Run benchmark
+        run: |
+          python3 test/benchmark.py

.github/workflows/code_formatting.yml

@@ -6,5 +6,25 @@ jobs:
   lint:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
-      - uses: psf/black@stable
+      - uses: actions/checkout@v4
+
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+      - name: Install linting tools
+        run: pip install mypy ruff
+
+      - name: ruff format
+        run: |
+          ruff format --check .
+
+      - name: ruff check
+        continue-on-error: true
+        run: |
+          ruff check .
+
+      - name: mypy
+        continue-on-error: true
+        run: |
+          python -m mypy .

.gitignore

@@ -105,3 +105,9 @@ venv.bak/
 
 # version file
 **_version.py
+
+# Input/Output files
+*.bp
+*.xdmf
+*.h5
+*.txt

docs/source/api/index.rst

@@ -0,0 +1,26 @@
+FESTIM API
+==========
+
+
+.. automodule:: festim.boundary_conditions
+    :members:
+    :show-inheritance:
+    :exclude-members: __weakref__
+    :private-members:
+    :inherited-members:
+
+
+.. .. automodule:: festim.hydrogen_transport_problem
+..     :members:
+..     :show-inheritance:
+..     :exclude-members: __weakref__
+..     :private-members:
+..     :inherited-members:
+
+
+
+.. .. automodule:: festim.heat_transfer_problem
+..     :members:
+..     :show-inheritance:
+..     :exclude-members: __weakref__
+..     :inherited-members:

docs/source/index.rst

@@ -34,6 +34,7 @@ Contents
    userguide/index
    devguide/index
    publications
+   api/index
 
 
 Indices and tables

convergence_rates.py → examples/convergence_rates.py

@@ -1,9 +1,11 @@
-import festim as F
-import numpy as np
-from dolfinx import fem
-import ufl
 import mpi4py.MPI as MPI
+
 import matplotlib.pyplot as plt
+import numpy as np
+import ufl
+from dolfinx import fem
+
+import festim as F
 
 
 def source_from_exact_solution(
@@ -58,19 +60,27 @@ def error_L2(u_computed, u_exact, degree_raise=3):
 
 
 def run(N):
-    exact_solution = lambda x, t: 2 * x[0] ** 2 + 20 * t
+    def exact_solution(x, t):
+        return 2 * x[0] ** 2 + 20 * t
+
+    def density(T):
+        return 0.2 * T + 2
 
-    density = lambda T: 0.2 * T + 2
-    heat_capacity = lambda T: 0.2 * T + 3
-    thermal_conductivity = lambda T: 0.1 * T + 4
+    def heat_capacity(T):
+        return 0.2 * T + 3
+
+    def thermal_conductivity(T):
+        return 0.1 * T + 4
 
     mms_source_from_sp = source_from_exact_solution(
         exact_solution,
         density=lambda x, t: density(exact_solution(x, t)),
         heat_capacity=lambda x, t: heat_capacity(exact_solution(x, t)),
         thermal_conductivity=lambda x, t: thermal_conductivity(exact_solution(x, t)),
     )
-    mms_source = lambda x, t: mms_source_from_sp((x[0], None, None), t)
+
+    def mms_source(x, t):
+        return mms_source_from_sp((x[0], None, None), t)
 
     my_problem = F.HeatTransferProblem()
 
@@ -110,18 +120,19 @@ def run(N):
     my_problem.settings.stepsize = F.Stepsize(0.05)
 
     my_problem.exports = [
-        F.VTXExportForTemperature(filename="test_transient_heat_transfer.bp")
+        F.VTXSpeciesExport(filename="test_transient_heat_transfer.bp")
     ]
 
     my_problem.initialise()
     my_problem.run()
 
     computed_solution = my_problem.u
-    final_time_sim = (
-        my_problem.t.value
-    )  # we use the exact final time of the simulation which may differ from the one specified in the settings
+    # we use the exact final time of the simulation which may differ from the one specified in the settings
+    final_time_sim = my_problem.t.value
+
+    def exact_solution_end(x):
+        return exact_solution(x, final_time_sim)
 
-    exact_solution_end = lambda x: exact_solution(x, final_time_sim)
     L2_error = error_L2(computed_solution, exact_solution_end)
     return L2_error
 

examples/multi_isotope_trapping_example.py

@@ -1,6 +1,6 @@
-import festim as F
 import numpy as np
 
+import festim as F
 
 my_model = F.HydrogenTransportProblem()
 

examples/multi_material_1d.py

@@ -0,0 +1,102 @@
+import numpy as np
+
+import festim as F
+
+my_model = F.HTransportProblemDiscontinuous()
+
+interface_1 = 0.5
+interface_2 = 0.7
+
+# for some reason if the mesh isn't fine enough then I have a random SEGV error
+N = 2000
+vertices = np.concatenate(
+    [
+        np.linspace(0, interface_1, num=N),
+        np.linspace(interface_1, interface_2, num=N),
+        np.linspace(interface_2, 1, num=N),
+    ]
+)
+
+my_model.mesh = F.Mesh1D(vertices)
+
+material_left = F.Material(D_0=2.0, E_D=0)
+material_mid = F.Material(D_0=2.0, E_D=0)
+material_right = F.Material(D_0=2.0, E_D=0)
+
+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=[vertices[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, vertices[-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 = [
+    F.Interface(6, (left_domain, middle_domain)),
+    F.Interface(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 my_model.species:
+    species.subdomains = my_model.volume_subdomains
+
+
+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=1e-5, transient=False)
+
+my_model.exports = [
+    F.VTXSpeciesExport(filename=f"u_{subdomain.id}.bp", field=H, subdomain=subdomain)
+    for subdomain in my_model.volume_subdomains
+]
+
+my_model.initialise()
+my_model.run()