Gas enclosure tutorials

CLAUDE.md

@@ -0,0 +1,79 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## What this is
+
+FESTIM-workshop is a **Jupyter Book tutorial** (published at festim-workshop.readthedocs.io) for
+[FESTIM](https://festim.readthedocs.io), a finite-element hydrogen-transport / heat-transfer code
+built on FEniCSx (dolfinx). This repo is teaching material — collections of executable notebooks —
+not a Python package. There is nothing to `pip install` from here and no unit test suite; "correct"
+means the book builds and every notebook runs top to bottom without error or warning.
+
+## Environment
+
+Everything runs in the `festim-workshop` conda env (dolfinx, FEniCSx, and pip-only deps cannot be
+installed with plain pip):
+
+```bash
+conda env create -f environment.yml   # first time
+conda activate festim-workshop
+```
+
+Key pinned deps: `festim=2.1`, `fenics-dolfinx`, `h-transport-materials` (material property DB),
+`foam2dolfinx` (OpenFOAM→dolfinx CFD coupling), `openmc2dolfinx` (neutronics), `autoemulate`+`torch`
+(ML surrogates), `python-gmsh`, `pyvista`/`trame` (3D viz). Python is pinned `<3.13`.
+
+## Build the book
+
+```bash
+jupyter-book build book          # outputs to book/_build/html/index.html
+```
+
+Read the Docs reproduces the build via a `pre_build` step then Sphinx — see `.readthedocs.yml`.
+Two build behaviors matter when editing content:
+
+- **Notebooks are force-executed on every build** (`execute_notebooks: force`, `timeout: -1` in
+  `book/_config.yml`). A change that makes any cell raise breaks the build. Expect long build times.
+- **`fail_on_warning: true`** (in `.readthedocs.yml`): a Sphinx warning — a broken cross-reference,
+  a missing image, a bad `codeautolink` — fails the RTD build even if notebooks execute fine.
+
+Run a single notebook without a full book build to check it executes:
+
+```bash
+jupyter nbconvert --to notebook --execute book/content/applications/task02.ipynb --stdout > /dev/null
+# MyST-markdown pages: convert first with jupytext, or open in Jupyter (they are paired notebooks)
+jupytext --to notebook book/content/material/material_basics.md -o /tmp/out.ipynb
+```
+
+## Content structure
+
+- `book/intro.md` — landing page; `book/_config.yml` — Jupyter Book config; `book/_toc.yml` — table
+  of contents. **Any new page must be registered in `book/_toc.yml` or it won't appear.**
+- `book/content/` is organized by concept, mirroring FESTIM's building blocks:
+  `meshes/`, `material/`, `boundary_conditions/`, `species_reactions/`, `initial_conditions/`,
+  `temperatures/`, `post_process/`, `misc/`, and `applications/` (end-to-end `taskNN` cases plus
+  `multiphysics/` CFD+neutronics and `ml/` surrogate/active-learning).
+- `book/content/misc/festim_from_scratch.md` + `dolfinx_*.md` show the same physics written in raw
+  dolfinx rather than the FESTIM API — useful when explaining what FESTIM does under the hood.
+
+### Page formats — two interchangeable notebook types
+
+Every tutorial page is an **executable notebook** in one of two formats; treat both as code:
+
+1. `.ipynb` — standard Jupyter notebooks (most `applications/taskNN`).
+2. `.md` with a `jupytext` YAML header (`formats: ipynb,md:myst`) — MyST-markdown notebooks. These
+   are notebooks stored as text and are executed by the build exactly like `.ipynb`. Editing the
+   Python inside a fenced ```` ```{code-cell} ```` block is editing runnable code.
+
+There are no plain prose `.md` files — every `.md` under `book/content/` is a MyST notebook.
+
+## Conventions
+
+- API objects in code cells auto-link to upstream docs via `sphinx-codeautolink` + `intersphinx`
+  (festim, dolfinx, numpy, matplotlib…). Reference real, importable names so links resolve; a bad
+  reference is a build warning.
+- `.bp`, `.h5`/`.xdmf`, `.msh` files committed next to notebooks are simulation outputs / meshes that
+  cells read or write — don't assume they're stale artifacts before checking the notebook that uses them.
+- Note: `.github/workflows/test_notebooks.yml` still targets an old flat `tasks/` directory that no
+  longer exists (content moved under `book/content/applications/`); that workflow is stale.

book/_static/new-badge.js

@@ -7,6 +7,12 @@
 // Badges the page's entry in the left-hand navigation sidebar.
 const NEW_PAGES = [
   "content/applications/ml",
+  "content/applications/task11",
+  "content/gas_enclosures/index",
+  "content/gas_enclosures/enclosure_intro",
+  "content/gas_enclosures/coupling",
+  "content/gas_enclosures/openings",
+  "content/gas_enclosures/geometries",
 ];
 
 // Badges a single section at its heading and in the right-hand "Contents" panel,

book/_toc.yml

@@ -49,7 +49,14 @@ parts:
       sections:
       - file: content/temperatures/temperatures_basic
       - file: content/temperatures/temperatures_advanced
-
+
+    - file: content/gas_enclosures/index
+      sections:
+      - file: content/gas_enclosures/enclosure_intro
+      - file: content/gas_enclosures/coupling
+      - file: content/gas_enclosures/openings
+      - file: content/gas_enclosures/geometries
+
     - file: content/post_process/intro
       sections:
       - file: content/post_process/derived
@@ -80,6 +87,7 @@ parts:
     - file: content/applications/task08
     - file: content/applications/task06
     - file: content/applications/task10
+    - file: content/applications/task11
     - file: content/applications/microstructure
     - file: content/applications/multiphysics/index
       sections:

book/content/boundary_conditions/h_transport_advanced.md

@@ -560,3 +560,73 @@ plt.show()
 ```{note}
 Increasing the penalty improves how well the target value is enforced, but very large values can degrade the conditioning of the linear system. In practice, the penalty should be chosen large enough to satisfy the constraint to the desired accuracy without harming convergence.
 ```
+
++++
+
+## Concentration-dependent inputs
+
+Since FESTIM 2.2, the `value` of a `festim.ParticleSource` or a `festim.ParticleFluxBC`
+can depend on the concentration of *another* species, in addition to space `x`, time `t` and temperature
+`T`. You declare the dependency with `species_dependent_value`, a dict that maps each species-argument name
+in your callable to a `festim.Species`. FESTIM binds those arguments to the live species solution,
+so the term is evaluated, and differentiated in the Newton solve, against the coupled concentration field.
+
+As an example, we drive a species `A` in from the left boundary and produce a second species `B` locally at
+a rate proportional to the concentration of `A`, i.e. $S_B = 0.5\,c_A$:
+
+```{code-cell} ipython3
+import festim as F
+import numpy as np
+from dolfinx.mesh import create_unit_interval
+from mpi4py import MPI
+
+model = F.HydrogenTransportProblem()
+model.mesh = F.Mesh(create_unit_interval(MPI.COMM_WORLD, 40))
+material = F.Material(D_0=1.0, E_D=0.0)
+vol = F.VolumeSubdomain(id=1, material=material)
+left = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[0], 0.0))
+model.subdomains = [vol, left]
+
+A = F.Species("A")
+B = F.Species("B")
+model.species = [A, B]
+model.temperature = 500.0
+
+# the argument name "c_A" in the callable is bound to species A's concentration
+model.sources = [
+    F.ParticleSource(
+        value=lambda c_A: 0.5 * c_A,
+        volume=vol,
+        species=B,
+        species_dependent_value={"c_A": A},
+    )
+]
+model.boundary_conditions = [F.FixedConcentrationBC(subdomain=left, value=1.0, species=A)]
+model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=True,
+                            final_time=0.3, stepsize=0.02)
+
+model.initialise()
+model.run()
+```
+
+`B` appears wherever `A` has diffused, without any explicit coupling term written by hand:
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import matplotlib.pyplot as plt
+
+x = A.post_processing_solution.function_space.tabulate_dof_coordinates()[:, 0]
+order = np.argsort(x)
+plt.plot(x[order], A.post_processing_solution.x.array[order], label="A (driven in)")
+plt.plot(x[order], B.post_processing_solution.x.array[order], label="B (produced from A)")
+plt.xlabel("x")
+plt.ylabel("concentration")
+plt.legend()
+plt.grid(alpha=0.3)
+plt.show()
+```
+
+The same `species_dependent_value` argument works on `festim.ParticleFluxBC`, letting a surface
+flux depend on a species' concentration, for example `value=lambda c: k * c` with
+`species_dependent_value={"c": my_species}`.

book/content/gas_enclosures/coupling.md

@@ -0,0 +1,191 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.18.1
+kernelspec:
+  display_name: festim-workshop
+  language: python
+  name: python3
+---
+
+# Coupling modes
+
+An enclosure can be coupled to the solid in two ways:
+
+* **Flux coupling**: a `festim.SurfaceReactionBC` whose `gas_pressure` is the `GasSpecies`.
+  The gas exchanges with the solid through an explicit surface reaction
+  ($2\,\mathrm{H} \rightleftharpoons \mathrm{H_2}$), as in [](enclosure_intro.md). Use this when the
+  recombination/dissociation kinetics matter.
+* **Dirichlet (concentration) coupling**: a `festim.SievertsBC` or `festim.HenrysBC`
+  whose `pressure` is the `GasSpecies`. The surface concentration is tied *instantaneously* to the gas
+  pressure through the solubility law ($c = S\sqrt{P}$ for Sieverts, $c = K_H P$ for Henry). Use this
+  when the surface is assumed to be at equilibrium with the gas.
+
+This page covers the Dirichlet mode and verifies it against analytical solutions.
+
+```{admonition} The pressure lives in a real function space
+:class: note
+Because the pressure is a single global unknown, it cannot be interpolated into a nodal field, so these
+boundary conditions can only be enforced **weakly**, with Nitsche's method. Pass `enforce_weakly=True`
+and a `penalty` coefficient.
+```
+
++++
+
+## A helper to build the slab + enclosure
+
+We reuse the same setup for both verifications: a 1D slab whose left surface is coupled to a gas
+enclosure by Henry's law, optionally with a perfect sink ($c=0$) on the right surface.
+
+```{code-cell} ipython3
+import numpy as np
+import festim as F
+
+
+def build(length, area, volume, D_0=1e-1, T=500.0, P0=1e5, K_H=1e15,
+          penalty=100, sink=False, final_time=4000.0, dt=20.0, n=60):
+    model = F.HydrogenTransportProblemDiscontinuous()
+    model.mesh = F.Mesh1D(np.linspace(0.0, length, num=n + 1))
+
+    material = F.Material(name="mat", D_0=D_0, E_D=0.0)
+    vol = F.VolumeSubdomain1D(id=1, borders=[0.0, length], material=material)
+    left = F.SurfaceSubdomain1D(id=1, x=0.0)
+    right = F.SurfaceSubdomain1D(id=2, x=length)
+    model.subdomains = [vol, left, right]
+
+    H = F.Species("H", subdomains=[vol])
+    model.species = [H]
+    model.temperature = T
+
+    H2 = F.GasSpecies(name="H2", initial_pressure=P0)
+    model.enclosures = [
+        F.Enclosure(volume=volume, species=[H2], temperature=T, surfaces={left: area})
+    ]
+
+    # Henry's law coupling: the surface concentration follows c = K_H * P(t)
+    model.boundary_conditions = [
+        F.HenrysBC(subdomain=left, H_0=K_H, E_H=0.0, pressure=H2, species=H,
+                   enforce_weakly=True, penalty=penalty)
+    ]
+    if sink:
+        model.boundary_conditions.append(
+            F.FixedConcentrationBC(subdomain=right, value=0.0, species=H)
+        )
+
+    model.initial_conditions = [F.InitialConcentration(value=0.0, species=H, volume=vol)]
+
+    pressure = F.GasPressure(field=H2)
+    model.exports = [pressure]
+
+    # atol must match the magnitude of the residual: with c ~ K_H*P ~ 1e20 the residual
+    # bottoms out well above the 1e-8 that would suit an order-1 problem.
+    model.settings = F.Settings(atol=1e8, rtol=1e-8, transient=True,
+                                final_time=final_time, stepsize=F.Stepsize(dt))
+    model.show_progress_bar = False
+    return model, H2, pressure
+```
+
+## Verification 1: a closed enclosure reaches equilibrium
+
+With no sink, the enclosure and the slab share a fixed number of hydrogen atoms. At equilibrium the
+concentration is uniform and in equilibrium with the gas ($c = K_H P_\infty$ everywhere), so conservation
+of atoms gives
+
+$$
+A\,\ell\,K_H\,P_\infty + \frac{P_\infty V}{k_B T} = \frac{P_0 V}{k_B T}
+\qquad\Longrightarrow\qquad
+P_\infty = \frac{P_0}{1 + A\,\ell\,K_H\,k_B T / V}.
+$$
+
+```{code-cell} ipython3
+length, area, V_enc, T, P0, K_H = 2.0, 0.25, 1e-6, 500.0, 1e5, 1e15
+
+model, H2, pressure = build(length, area, V_enc, T=T, P0=P0, K_H=K_H)
+model.initialise()
+model.run()
+
+kT = F.k_B_SI * T
+expected = P0 / (1 + area * length * K_H * kT / V_enc)
+
+print(f"final pressure : {H2.value:.6e} Pa")
+print(f"analytical     : {expected:.6e} Pa")
+assert abs(H2.value - expected) / expected < 1e-4
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import matplotlib.pyplot as plt
+
+plt.plot(pressure.t, pressure.data, label="FESTIM")
+plt.axhline(expected, color="k", ls="--", label="analytical equilibrium")
+plt.xlabel("time (s)")
+plt.ylabel("H2 pressure (Pa)")
+plt.legend()
+plt.grid(alpha=0.3)
+plt.show()
+```
+
+## Verification 2: the TMAP permeation case (issue #996)
+
+Now add a perfect sink on the far side. Hydrogen coupled in through Henry's law on the left permeates
+across the slab and is pumped away on the right, so the enclosure pressure decays. Separating variables
+gives the eigenvalue equation
+
+$$
+\lambda \tan\lambda = \beta, \qquad \beta = \frac{k_B T\,K_H\,A\,\ell}{V},
+$$
+
+and at late times the pressure decays as $\exp(-D\,\lambda_1^2\,t/\ell^2)$, where $\lambda_1$ is the
+first root.
+
+```{code-cell} ipython3
+import scipy.optimize
+
+D_0 = 1e-1
+beta = F.k_B_SI * T * K_H * area * length / V_enc
+lambda_1 = scipy.optimize.brentq(
+    lambda lam: lam * np.tan(lam) - beta, 1e-12, np.pi / 2 - 1e-9
+)
+continuous_rate = D_0 * lambda_1**2 / length**2
+print(f"beta = {beta:.3f},  lambda_1 = {lambda_1:.4f},  decay rate = {continuous_rate:.3e} 1/s")
+```
+
+```{code-cell} ipython3
+model, H2, pressure = build(length, area, V_enc, D_0=D_0, T=T, P0=P0, K_H=K_H,
+                            sink=True, final_time=300.0, dt=2.0)
+model.initialise()
+model.run()
+
+t = np.array(pressure.t)
+P = np.array(pressure.data)
+
+# fit the rate where the first eigenmode dominates (higher modes gone, before the solve
+# stalls on tolerances)
+window = (P < 1e-2 * P[0]) & (P > 1e-6 * P[0])
+fitted_rate = -np.polyfit(t[window], np.log(P[window]), 1)[0]
+
+# backward Euler integrates dP/dt = -r P as P_{n+1} = P_n/(1 + r*dt), a log-slope of
+# -ln(1 + r*dt)/dt; compare against that discrete rate rather than the continuous one.
+dt = 2.0
+discrete_rate = np.log(1 + continuous_rate * dt) / dt
+print(f"fitted rate   : {fitted_rate:.4e} 1/s")
+print(f"discrete rate : {discrete_rate:.4e} 1/s")
+assert abs(fitted_rate - discrete_rate) / discrete_rate < 5e-3
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+plt.semilogy(t, P, label="FESTIM")
+plt.semilogy(t, P[0] * np.exp(-fitted_rate * t), "k--", label="fitted first mode")
+plt.xlabel("time (s)")
+plt.ylabel("H2 pressure (Pa)")
+plt.legend()
+plt.grid(alpha=0.3, which="both")
+plt.show()
+```