From dd542d49119c83059465efbe1885538a9f7531e6 Mon Sep 17 00:00:00 2001 From: ckhurana Date: Wed, 8 Oct 2025 10:47:52 -0400 Subject: [PATCH 01/17] boundary conditions chapter --- book/_toc.yml | 5 + .../boundary_conditions/concentration.ipynb | 302 ++++++++++++++++++ book/content/boundary_conditions/fluxes.ipynb | 42 +++ 3 files changed, 349 insertions(+) create mode 100644 book/content/boundary_conditions/concentration.ipynb create mode 100644 book/content/boundary_conditions/fluxes.ipynb diff --git a/book/_toc.yml b/book/_toc.yml index 6cac4fa6..40cd729b 100644 --- a/book/_toc.yml +++ b/book/_toc.yml @@ -21,6 +21,11 @@ parts: - file: content/material/material_htm - file: content/material/material_advanced + - caption: Boundary Conditions + chapters: + - file: content/boundary_conditions/concentration + - file: content/boundary_conditions/fluxes + - caption: Species & Reactions chapters: - file: content/species_reactions/species diff --git a/book/content/boundary_conditions/concentration.ipynb b/book/content/boundary_conditions/concentration.ipynb new file mode 100644 index 00000000..7109ce00 --- /dev/null +++ b/book/content/boundary_conditions/concentration.ipynb @@ -0,0 +1,302 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b24da17d", + "metadata": {}, + "source": [ + "# Concentration BCs #\n", + "\n", + "Boundary conditions (BCs) are essential to FESTIM simulations, as they describe the mathematical problem at the boundaries of the simulated domain. Read more about BCs _[here](https://festim.readthedocs.io/en/fenicsx/userguide/boundary_conditions.html)_.\n", + "\n", + "This tutorial goes over how to define concentration boundary conditions for hydrogen transport simulations.\n", + "\n", + "Objectives:\n", + "* Define a fixed concentration BC\n", + "* Choose which solubility law\n", + "* Plasma implantation approximation" + ] + }, + { + "cell_type": "markdown", + "id": "f018139f", + "metadata": {}, + "source": [ + "## Defining fixed concentration ##\n", + "\n", + "BCs need to be assigned to surfaces using FESTIM's `SurfaceSubdomain` class. To define the concentration of a specific species, we can use `FixedConcentrationBC`:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "c9b8fdc2", + "metadata": {}, + "outputs": [], + "source": [ + "from festim import FixedConcentrationBC, Species, SurfaceSubdomain\n", + "\n", + "boundary = SurfaceSubdomain(id=1)\n", + "H = Species(name=\"Hydrogen\")\n", + "\n", + "my_bc = FixedConcentrationBC(subdomain=boundary, value=10, species=H)" + ] + }, + { + "cell_type": "markdown", + "id": "0082375d", + "metadata": {}, + "source": [ + "The imposed concentration can be dependent on space, time and temperature, such as:\n", + "\n", + "$$ \n", + "\n", + "BC = 10 + x^2 + T(t+2)\n", + "\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7f47ff10", + "metadata": {}, + "outputs": [], + "source": [ + "my_custom_value = lambda x, t, T: 10 + x[0]**2 + T *(t + 2)\n", + "\n", + "my_bc = FixedConcentrationBC(subdomain=boundary, value=my_custom_value, species=H)" + ] + }, + { + "cell_type": "markdown", + "id": "b55d56c7", + "metadata": {}, + "source": [ + "Users can define the surface concentration using either Sieverts’ law, $c = S(T)\\sqrt P$, or Henry's law, $c=K_H P$, where $S(T)$ and $K_H$ denote the temperature-dependent Sieverts’ and Henry’s solubility coefficients, respectively, and $P$ is the partial pressure of the species on the surface. \n", + "\n", + "For Sieverts' law of solubility, we can use `festim.SievertsBC`:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "8c72893c", + "metadata": {}, + "outputs": [], + "source": [ + "from festim import SievertsBC, SurfaceSubdomain, Species\n", + "\n", + "boundary = SurfaceSubdomain(id=1)\n", + "H = Species(name=\"Hydrogen\")\n", + "\n", + "custom_pressure_value = lambda t: 2 + t\n", + "my_bc = SievertsBC(subdomain=3, S_0=2, E_S=0.1, species=H, pressure=custom_pressure_value)" + ] + }, + { + "cell_type": "markdown", + "id": "53ab8ee2", + "metadata": {}, + "source": [ + "Similarly, for Henry's law of solubility, we can use `festim.HenrysBC`:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "38b13107", + "metadata": {}, + "outputs": [], + "source": [ + "from festim import HenrysBC\n", + "\n", + "pressure_value = lambda t: 5 * t\n", + "my_bc = HenrysBC(subdomain=3, H_0=1.5, E_H=0.2, species=H, pressure=pressure_value)" + ] + }, + { + "cell_type": "markdown", + "id": "1f18a9f7", + "metadata": {}, + "source": [ + "## Plasma implantation ##" + ] + }, + { + "cell_type": "markdown", + "id": "5335a5cb", + "metadata": {}, + "source": [ + "We can also approximate plasma implantation using FESTIM's `ParticleSource` class, which is helpful in modeling thermo-desorption spectra (TDS) experiments or including the effect of plasma exposure on hydrogen transport. Learn more about how FESTIM approximates plasma implantation _[here](https://festim.readthedocs.io/en/fenicsx/theory.html)_.\n", + "\n", + "We can define a representative source term:\n", + "\n", + "$$ S_{ext} = \\varphi \\cdot f(x) $$\n", + "\n", + "where $\\varphi$ is the implantation flux and $f(x)$ is a Gaussian spatial distribution. For example, if we had an implantation flux $\\varphi = 1 \\mathrm{m}^{-2}\\mathrm{s}^{-1}$ and a distribution mean value of 1 $\\text{nm}$ and width of 3 $\\text{nm}$:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "6e467cbc", + "metadata": {}, + "outputs": [], + "source": [ + "import festim as F\n", + "import ufl\n", + "import numpy as np\n", + "\n", + "my_model = F.HydrogenTransportProblem()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a43e3993", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "ld: warning: duplicate -rpath '/Users/ckhurana/miniconda3/envs/festim-workshop/lib' ignored\n", + "ld: warning: duplicate -rpath '/Users/ckhurana/miniconda3/envs/festim-workshop/lib' ignored\n", + "ld: warning: duplicate -rpath '/Users/ckhurana/miniconda3/envs/festim-workshop/lib' ignored\n", + "ld: warning: duplicate -rpath '/Users/ckhurana/miniconda3/envs/festim-workshop/lib' ignored\n" + ] + } + ], + "source": [ + "\n", + "vertices = np.concatenate(\n", + " [\n", + " np.linspace(0, 30e-9, num=200),\n", + " np.linspace(30e-9, 3e-6, num=300),\n", + " np.linspace(3e-6, 20e-6, num=200),\n", + " ]\n", + ")\n", + "\n", + "my_model.mesh = F.Mesh1D(vertices)\n", + "\n", + "mat = F.Material(D_0=1e-7, E_D=0.2)\n", + "\n", + "volume_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 20e-6], material=mat)\n", + "left_boundary = F.SurfaceSubdomain1D(id=1, x=0)\n", + "right_boundary = F.SurfaceSubdomain1D(id=2, x=20e-6)\n", + "incident_flux = 1e12 # H/m2/s\n", + "\n", + "my_model.subdomains = [\n", + " volume_subdomain,\n", + " left_boundary,\n", + " right_boundary,\n", + "]\n", + "def gaussian_distribution(x, center, width):\n", + " return (\n", + " 1\n", + " / (width * (2 * ufl.pi) ** 0.5)\n", + " * ufl.exp(-0.5 * ((x[0] - center) / width) ** 2)\n", + " )\n", + "H = F.Species(\"H\")\n", + "my_model.species = [H]\n", + "\n", + "source_term = F.ParticleSource(\n", + " value=lambda x: incident_flux * gaussian_distribution(x, 10e-6, 20e-6),\n", + " volume=volume_subdomain,\n", + " species=H,\n", + ")\n", + "\n", + "my_model.sources = [source_term]\n", + "\n", + "my_model.boundary_conditions = [\n", + " F.FixedConcentrationBC(subdomain=left_boundary, value=0, species=H),\n", + " F.FixedConcentrationBC(subdomain=right_boundary, value=0, species=H),\n", + "]\n", + "\n", + "my_model.temperature = 400\n", + "my_model.settings = F.Settings(atol=1e10, rtol=1e-10, transient=False)\n", + "\n", + "profile_export = F.Profile1DExport(field=H,subdomain=volume_subdomain)\n", + "my_model.exports = [profile_export]\n", + "my_model.initialise()\n", + "my_model.run()\n", + "# " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "e966568c", + "metadata": {}, + "outputs": [], + "source": [ + "x = my_model.exports[0].x\n", + "c = my_model.exports[0].data " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "abbb180f", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGvCAYAAABW/q+QAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAI0NJREFUeJzt3X1UlGUC9/HfxMugJZOKgiQqtiaSbkcxEYu1TopiWe66qbmRdcyNbc1FTiff2ieyjqRb5rr4srqYdU6p2yLmORkrHRUt0dIFLSO3jJJVJsJVhrQFxfv5o8d5mhjQQUaZq+/nnPmDe677nvtquA7f7nnRZlmWJQAAAINcc7VPAAAAoLUROAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACME3y1T+BqOH/+vI4fP64OHTrIZrNd7dMBAACXwLIs1dbWKjo6Wtdc0/w1mp9k4Bw/flwxMTFX+zQAAEALVFRUqHv37s2O+UkGTocOHSR9/x8oPDz8Kp8NAAC4FC6XSzExMe6/4835SQbOhZelwsPDCRwAAALMpby9hDcZAwAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADDOFQmc5cuXKzY2VmFhYUpISNCuXbuaHV9UVKSEhASFhYWpd+/eWrlyZZNj169fL5vNpnHjxrXyWQMAgEDl98DZsGGDMjIyNG/ePJWUlCg5OVmpqak6evSo1/Hl5eUaM2aMkpOTVVJSorlz52rGjBnKy8trNParr77Sk08+qeTkZH9PAwAABBCbZVmWPx8gMTFRgwYN0ooVK9zb+vXrp3Hjxik7O7vR+FmzZmnz5s0qKytzb0tPT9eBAwdUXFzs3tbQ0KDhw4frkUce0a5du3Tq1Clt2rTpks7J5XLJ4XCopqZG4eHhLZ8cAAC4Ynz5++3XKzj19fXav3+/UlJSPLanpKRo9+7dXvcpLi5uNH7UqFHat2+fzp496942f/58denSRVOnTr3oedTV1cnlcnncAACAufwaONXV1WpoaFBkZKTH9sjISDmdTq/7OJ1Or+PPnTun6upqSdL777+v3NxcrV69+pLOIzs7Ww6Hw32LiYlpwWwAAECguCJvMrbZbB4/W5bVaNvFxl/YXltbqwcffFCrV69WRETEJT3+nDlzVFNT475VVFT4OAMAABBIgv158IiICAUFBTW6WlNVVdXoKs0FUVFRXscHBwerc+fOOnTokL788kuNHTvWff/58+clScHBwTp8+LBuvPFGj/3tdrvsdntrTAkAAAQAv17BCQ0NVUJCggoLCz22FxYWatiwYV73SUpKajR+69atGjx4sEJCQhQXF6ePPvpIpaWl7tu9996rO++8U6Wlpbz8BAAA/HsFR5IyMzOVlpamwYMHKykpSatWrdLRo0eVnp4u6fuXj44dO6bXXntN0vefmMrJyVFmZqamTZum4uJi5ebmat26dZKksLAw9e/f3+Mxrr/+eklqtB0AAPw0+T1wJk6cqBMnTmj+/PmqrKxU//79tWXLFvXs2VOSVFlZ6fGdOLGxsdqyZYtmzpypZcuWKTo6WkuXLtX48eP9faoAAMAQfv8enLaI78EBACDwtJnvwQEAALgaCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxrkigbN8+XLFxsYqLCxMCQkJ2rVrV7Pji4qKlJCQoLCwMPXu3VsrV670uH/16tVKTk5Wx44d1bFjR40YMUIffPCBP6cAAAACiN8DZ8OGDcrIyNC8efNUUlKi5ORkpaam6ujRo17Hl5eXa8yYMUpOTlZJSYnmzp2rGTNmKC8vzz1mx44deuCBB7R9+3YVFxerR48eSklJ0bFjx/w9HQAAEABslmVZ/nyAxMREDRo0SCtWrHBv69evn8aNG6fs7OxG42fNmqXNmzerrKzMvS09PV0HDhxQcXGx18doaGhQx44dlZOTo4ceeuii5+RyueRwOFRTU6Pw8PAWzAoAAFxpvvz99usVnPr6eu3fv18pKSke21NSUrR7926v+xQXFzcaP2rUKO3bt09nz571us+ZM2d09uxZderUyev9dXV1crlcHjcAAGAuvwZOdXW1GhoaFBkZ6bE9MjJSTqfT6z5Op9Pr+HPnzqm6utrrPrNnz9YNN9ygESNGeL0/OztbDofDfYuJiWnBbAAAQKC4Im8yttlsHj9bltVo28XGe9suSYsWLdK6deu0ceNGhYWFeT3enDlzVFNT475VVFT4OgUAABBAgv158IiICAUFBTW6WlNVVdXoKs0FUVFRXscHBwerc+fOHttffPFFLViwQO+++65+/vOfN3kedrtddru9hbMAAACBxq9XcEJDQ5WQkKDCwkKP7YWFhRo2bJjXfZKSkhqN37p1qwYPHqyQkBD3tj/96U967rnnVFBQoMGDB7f+yQMAgIDl95eoMjMz9be//U1r1qxRWVmZZs6cqaNHjyo9PV3S9y8f/fCTT+np6frqq6+UmZmpsrIyrVmzRrm5uXryySfdYxYtWqSnn35aa9asUa9eveR0OuV0OvXtt9/6ezoAACAA+PUlKkmaOHGiTpw4ofnz56uyslL9+/fXli1b1LNnT0lSZWWlx3fixMbGasuWLZo5c6aWLVum6OhoLV26VOPHj3ePWb58uerr6/XrX//a47GeeeYZZWVl+XtKAACgjfP79+C0RXwPDgAAgafNfA8OAADA1UDgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADDOFQmc5cuXKzY2VmFhYUpISNCuXbuaHV9UVKSEhASFhYWpd+/eWrlyZaMxeXl5io+Pl91uV3x8vPLz8/11+gAAIMD4PXA2bNigjIwMzZs3TyUlJUpOTlZqaqqOHj3qdXx5ebnGjBmj5ORklZSUaO7cuZoxY4by8vLcY4qLizVx4kSlpaXpwIEDSktL04QJE7R3715/TwcAAAQAm2VZlj8fIDExUYMGDdKKFSvc2/r166dx48YpOzu70fhZs2Zp8+bNKisrc29LT0/XgQMHVFxcLEmaOHGiXC6X3nnnHfeY0aNHq2PHjlq3bt1Fz8nlcsnhcKimpkbh4eGXMz0PlmXpu7MNrXY8AAACWbuQINlstlY7ni9/v4Nb7VG9qK+v1/79+zV79myP7SkpKdq9e7fXfYqLi5WSkuKxbdSoUcrNzdXZs2cVEhKi4uJizZw5s9GYJUuWeD1mXV2d6urq3D+7XK4WzObivjvboPj/80+/HBsAgEDzyfxRah/q19Rokl9foqqurlZDQ4MiIyM9tkdGRsrpdHrdx+l0eh1/7tw5VVdXNzumqWNmZ2fL4XC4bzExMS2dEgAACABXJKt+fHnKsqxmL1l5G//j7b4cc86cOcrMzHT/7HK5/BI57UKC9Mn8Ua1+XAAAAlG7kKCr9th+DZyIiAgFBQU1urJSVVXV6ArMBVFRUV7HBwcHq3Pnzs2OaeqYdrtddru9pdO4ZDab7apdigMAAP+fX1+iCg0NVUJCggoLCz22FxYWatiwYV73SUpKajR+69atGjx4sEJCQpod09QxAQDAT4vfLzdkZmYqLS1NgwcPVlJSklatWqWjR48qPT1d0vcvHx07dkyvvfaapO8/MZWTk6PMzExNmzZNxcXFys3N9fh01B/+8Af94he/0MKFC3Xffffprbfe0rvvvqv33nvP39MBAAABwO+BM3HiRJ04cULz589XZWWl+vfvry1btqhnz56SpMrKSo/vxImNjdWWLVs0c+ZMLVu2TNHR0Vq6dKnGjx/vHjNs2DCtX79eTz/9tP74xz/qxhtv1IYNG5SYmOjv6QAAgADg9+/BaYv89T04AADAf3z5+82/RQUAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACM49fAOXnypNLS0uRwOORwOJSWlqZTp041u49lWcrKylJ0dLTatWunO+64Q4cOHXLf/9///ldPPPGE+vbtq/bt26tHjx6aMWOGampq/DkVAAAQQPwaOJMnT1ZpaakKCgpUUFCg0tJSpaWlNbvPokWLtHjxYuXk5OjDDz9UVFSURo4cqdraWknS8ePHdfz4cb344ov66KOPtHbtWhUUFGjq1Kn+nAoAAAggNsuyLH8cuKysTPHx8dqzZ48SExMlSXv27FFSUpI+/fRT9e3bt9E+lmUpOjpaGRkZmjVrliSprq5OkZGRWrhwoR577DGvj/Xmm2/qwQcf1OnTpxUcHHzRc3O5XHI4HKqpqVF4ePhlzBIAAFwpvvz99tsVnOLiYjkcDnfcSNLQoUPlcDi0e/dur/uUl5fL6XQqJSXFvc1ut2v48OFN7iPJPdGm4qaurk4ul8vjBgAAzOW3wHE6neratWuj7V27dpXT6WxyH0mKjIz02B4ZGdnkPidOnNBzzz3X5NUdScrOzna/D8jhcCgmJuZSpwEAAAKQz4GTlZUlm83W7G3fvn2SJJvN1mh/y7K8bv+hH9/f1D4ul0t333234uPj9cwzzzR5vDlz5qimpsZ9q6iouJSpAgCAAHXxN6z8yPTp0zVp0qRmx/Tq1UsHDx7U119/3ei+b775ptEVmguioqIkfX8lp1u3bu7tVVVVjfapra3V6NGjdd111yk/P18hISFNno/dbpfdbm/2nAEAgDl8DpyIiAhFRERcdFxSUpJqamr0wQcfaMiQIZKkvXv3qqamRsOGDfO6T2xsrKKiolRYWKiBAwdKkurr61VUVKSFCxe6x7lcLo0aNUp2u12bN29WWFiYr9MAAAAG89t7cPr166fRo0dr2rRp2rNnj/bs2aNp06bpnnvu8fgEVVxcnPLz8yV9/9JURkaGFixYoPz8fH388cd6+OGH1b59e02ePFnS91duUlJSdPr0aeXm5srlcsnpdMrpdKqhocFf0wEAAAHE5ys4vnj99dc1Y8YM96ei7r33XuXk5HiMOXz4sMeX9D311FP67rvv9Pjjj+vkyZNKTEzU1q1b1aFDB0nS/v37tXfvXknSz372M49jlZeXq1evXn6cEQAACAR++x6ctozvwQEAIPC0ie/BAQAAuFoIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADG8WvgnDx5UmlpaXI4HHI4HEpLS9OpU6ea3ceyLGVlZSk6Olrt2rXTHXfcoUOHDjU5NjU1VTabTZs2bWr9CQAAgIDk18CZPHmySktLVVBQoIKCApWWliotLa3ZfRYtWqTFixcrJydHH374oaKiojRy5EjV1tY2GrtkyRLZbDZ/nT4AAAhQwf46cFlZmQoKCrRnzx4lJiZKklavXq2kpCQdPnxYffv2bbSPZVlasmSJ5s2bp1/96leSpFdffVWRkZF644039Nhjj7nHHjhwQIsXL9aHH36obt26+WsaAAAgAPntCk5xcbEcDoc7biRp6NChcjgc2r17t9d9ysvL5XQ6lZKS4t5mt9s1fPhwj33OnDmjBx54QDk5OYqKirroudTV1cnlcnncAACAufwWOE6nU127dm20vWvXrnI6nU3uI0mRkZEe2yMjIz32mTlzpoYNG6b77rvvks4lOzvb/T4gh8OhmJiYS50GAAAIQD4HTlZWlmw2W7O3ffv2SZLX98dYlnXR9838+P4f7rN582Zt27ZNS5YsueRznjNnjmpqaty3ioqKS94XAAAEHp/fgzN9+nRNmjSp2TG9evXSwYMH9fXXXze675tvvml0heaCCy83OZ1Oj/fVVFVVuffZtm2bjhw5ouuvv95j3/Hjxys5OVk7duxodFy73S673d7sOQMAAHP4HDgRERGKiIi46LikpCTV1NTogw8+0JAhQyRJe/fuVU1NjYYNG+Z1n9jYWEVFRamwsFADBw6UJNXX16uoqEgLFy6UJM2ePVuPPvqox34DBgzQyy+/rLFjx/o6HQAAYCC/fYqqX79+Gj16tKZNm6a//vWvkqTf/va3uueeezw+QRUXF6fs7Gz98pe/lM1mU0ZGhhYsWKA+ffqoT58+WrBggdq3b6/JkydL+v4qj7c3Fvfo0UOxsbH+mg4AAAggfgscSXr99dc1Y8YM96ei7r33XuXk5HiMOXz4sGpqatw/P/XUU/ruu+/0+OOP6+TJk0pMTNTWrVvVoUMHf54qAAAwiM2yLOtqn8SV5nK55HA4VFNTo/Dw8Kt9OgAA4BL48vebf4sKAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYJvtoncDVYliVJcrlcV/lMAADApbrwd/vC3/Hm/CQDp7a2VpIUExNzlc8EAAD4qra2Vg6Ho9kxNutSMsgw58+f1/Hjx9WhQwfZbLZWPbbL5VJMTIwqKioUHh7eqsduC0yfn2T+HJlf4DN9jswv8PlrjpZlqba2VtHR0brmmubfZfOTvIJzzTXXqHv37n59jPDwcGN/cSXz5yeZP0fmF/hMnyPzC3z+mOPFrtxcwJuMAQCAcQgcAABgHAKnldntdj3zzDOy2+1X+1T8wvT5SebPkfkFPtPnyPwCX1uY40/yTcYAAMBsXMEBAADGIXAAAIBxCBwAAGAcAgcAABiHwLmI5cuXKzY2VmFhYUpISNCuXbuaHV9UVKSEhASFhYWpd+/eWrlyZaMxeXl5io+Pl91uV3x8vPLz8/11+pfElzlu3LhRI0eOVJcuXRQeHq6kpCT985//9Bizdu1a2Wy2Rrf//e9//p6KV77Mb8eOHV7P/dNPP/UY15aeQ1/m9/DDD3ud38033+we05aev507d2rs2LGKjo6WzWbTpk2bLrpPoK1BX+cYaGvQ1/kF4hr0dY6BtA6zs7N16623qkOHDuratavGjRunw4cPX3S/trAOCZxmbNiwQRkZGZo3b55KSkqUnJys1NRUHT161Ov48vJyjRkzRsnJySopKdHcuXM1Y8YM5eXluccUFxdr4sSJSktL04EDB5SWlqYJEyZo7969V2paHnyd486dOzVy5Eht2bJF+/fv15133qmxY8eqpKTEY1x4eLgqKys9bmFhYVdiSh58nd8Fhw8f9jj3Pn36uO9rS8+hr/P785//7DGviooKderUSffff7/HuLby/J0+fVq33HKLcnJyLml8IK5BX+cYaGvQ1/ldEChrUPJ9joG0DouKivT73/9ee/bsUWFhoc6dO6eUlBSdPn26yX3azDq00KQhQ4ZY6enpHtvi4uKs2bNnex3/1FNPWXFxcR7bHnvsMWvo0KHunydMmGCNHj3aY8yoUaOsSZMmtdJZ+8bXOXoTHx9vPfvss+6fX3nlFcvhcLTWKV4WX+e3fft2S5J18uTJJo/Zlp7Dy33+8vPzLZvNZn355ZfubW3p+fshSVZ+fn6zYwJxDf7QpczRm7a8Bn/oUuYXaGvwx1ryHAbSOqyqqrIkWUVFRU2OaSvrkCs4Taivr9f+/fuVkpLisT0lJUW7d+/2uk9xcXGj8aNGjdK+fft09uzZZsc0dUx/askcf+z8+fOqra1Vp06dPLZ/++236tmzp7p376577rmn0f9dXgmXM7+BAweqW7duuuuuu7R9+3aP+9rKc9gaz19ubq5GjBihnj17emxvC89fSwTaGmwNbXkNXo5AWIOtJZDWYU1NjSQ1+n37obayDgmcJlRXV6uhoUGRkZEe2yMjI+V0Or3u43Q6vY4/d+6cqqurmx3T1DH9qSVz/LGXXnpJp0+f1oQJE9zb4uLitHbtWm3evFnr1q1TWFiYbrvtNn322Wetev4X05L5devWTatWrVJeXp42btyovn376q677tLOnTvdY9rKc3i5z19lZaXeeecdPfroox7b28rz1xKBtgZbQ1tegy0RSGuwNQTSOrQsS5mZmbr99tvVv3//Jse1lXX4k/zXxH1hs9k8frYsq9G2i43/8XZfj+lvLT2fdevWKSsrS2+99Za6du3q3j506FANHTrU/fNtt92mQYMG6S9/+YuWLl3aeid+iXyZX9++fdW3b1/3z0lJSaqoqNCLL76oX/ziFy06pr+19FzWrl2r66+/XuPGjfPY3taeP18F4hpsqUBZg74IxDV4OQJpHU6fPl0HDx7Ue++9d9GxbWEdcgWnCREREQoKCmpUk1VVVY2q84KoqCiv44ODg9W5c+dmxzR1TH9qyRwv2LBhg6ZOnaq///3vGjFiRLNjr7nmGt16661X/P88Lmd+PzR06FCPc28rz+HlzM+yLK1Zs0ZpaWkKDQ1tduzVev5aItDW4OUIhDXYWtrqGrxcgbQOn3jiCW3evFnbt29X9+7dmx3bVtYhgdOE0NBQJSQkqLCw0GN7YWGhhg0b5nWfpKSkRuO3bt2qwYMHKyQkpNkxTR3Tn1oyR+n7/2t8+OGH9cYbb+juu+++6ONYlqXS0lJ169btss/ZFy2d34+VlJR4nHtbeQ4vZ35FRUX6/PPPNXXq1Is+ztV6/loi0NZgSwXKGmwtbXUNXq5AWIeWZWn69OnauHGjtm3bptjY2Ivu02bWYau9XdlA69evt0JCQqzc3Fzrk08+sTIyMqxrr73W/U732bNnW2lpae7xX3zxhdW+fXtr5syZ1ieffGLl5uZaISEh1j/+8Q/3mPfff98KCgqyXnjhBausrMx64YUXrODgYGvPnj1XfH6W5fsc33jjDSs4ONhatmyZVVlZ6b6dOnXKPSYrK8sqKCiwjhw5YpWUlFiPPPKIFRwcbO3du7fNz+/ll1+28vPzrX//+9/Wxx9/bM2ePduSZOXl5bnHtKXn0Nf5XfDggw9aiYmJXo/Zlp6/2tpaq6SkxCopKbEkWYsXL7ZKSkqsr776yrIsM9agr3MMtDXo6/wCbQ1alu9zvCAQ1uHvfvc7y+FwWDt27PD4fTtz5ox7TFtdhwTORSxbtszq2bOnFRoaag0aNMjjo3FTpkyxhg8f7jF+x44d1sCBA63Q0FCrV69e1ooVKxod880337T69u1rhYSEWHFxcR4L92rwZY7Dhw+3JDW6TZkyxT0mIyPD6tGjhxUaGmp16dLFSklJsXbv3n0FZ+TJl/ktXLjQuvHGG62wsDCrY8eO1u233269/fbbjY7Zlp5DX39HT506ZbVr185atWqV1+O1pefvwkeGm/p9M2EN+jrHQFuDvs4vENdgS35PA2UdepuXJOuVV15xj2mr69D2/yYAAABgDN6DAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AAGjSzp07NXbsWEVHR8tms2nTpk1+fbysrCzZbDaPW1RUlM/HIXAAAECTTp8+rVtuuUU5OTlX7DFvvvlmVVZWum8fffSRz8cI9sN5AQAAQ6Smpio1NbXJ++vr6/X000/r9ddf16lTp9S/f38tXLhQd9xxR4sfMzg4uEVXbX6IKzgAAKDFHnnkEb3//vtav369Dh48qPvvv1+jR4/WZ5991uJjfvbZZ4qOjlZsbKwmTZqkL774wudj8G9RAQCAS2Kz2ZSfn69x48ZJko4cOaI+ffroP//5j6Kjo93jRowYoSFDhmjBggU+P8Y777yjM2fO6KabbtLXX3+t559/Xp9++qkOHTqkzp07X/JxuIIDAABa5F//+pcsy9JNN92k6667zn0rKirSkSNHJElffvllozcN//g2ffp09zFTU1M1fvx4DRgwQCNGjNDbb78tSXr11Vd9OjfegwMAAFrk/PnzCgoK0v79+xUUFORx33XXXSdJuuGGG1RWVtbscTp27Njkfddee60GDBjg80teBA4AAGiRgQMHqqGhQVVVVUpOTvY6JiQkRHFxcS1+jLq6OpWVlTV5/KYQOAAAoEnffvutPv/8c/fP5eXlKi0tVadOnXTTTTfpN7/5jR566CG99NJLGjhwoKqrq7Vt2zYNGDBAY8aM8fnxnnzySY0dO1Y9evRQVVWVnn/+eblcLk2ZMsWn4/AmYwAA0KQdO3bozjvvbLR9ypQpWrt2rc6ePavnn39er732mo4dO6bOnTsrKSlJzz77rAYMGODz402aNEk7d+5UdXW1unTpoqFDh+q5555TfHy8T8chcAAAgHH4FBUAADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4/xc5I//Q/+oy9QAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "c = np.array(c)\n", + "\n", + "plt.plot(x, c.squeeze())\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "13b71fb6", + "metadata": {}, + "source": [] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "7579b4a1", + "metadata": {}, + "outputs": [], + "source": [ + "c_flat = c[0, 0, :]\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "festim-workshop", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.18" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/book/content/boundary_conditions/fluxes.ipynb b/book/content/boundary_conditions/fluxes.ipynb new file mode 100644 index 00000000..da26f665 --- /dev/null +++ b/book/content/boundary_conditions/fluxes.ipynb @@ -0,0 +1,42 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "f633d832", + "metadata": {}, + "source": [ + "# Flux BCs #\n", + "\n", + "This tutorial goes over how to define flux boundary conditions in FESTIM.\n", + "\n", + "Objectives:\n", + "* Neuman: flux at boundaries\n", + "* Robin (based on concentrations with arbitrary expression): gradient is imposed as a function of the solutiin\n", + "* Surface reactions" + ] + }, + { + "cell_type": "markdown", + "id": "528d839b", + "metadata": {}, + "source": [ + "## Neumann BCs ##\n", + "\n", + "Users can impose Neumann boundary conditions to define the flux at boundaries. To impose a particle flux, we can use the " + ] + }, + { + "cell_type": "markdown", + "id": "7c26a601", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 3e3dc3ab84f47a8aa186c8b440c2938f6e18976a Mon Sep 17 00:00:00 2001 From: ckhurana Date: Tue, 14 Oct 2025 16:49:42 -0400 Subject: [PATCH 02/17] initial push for boundary condition chapters --- book/_toc.yml | 2 + .../boundary_conditions/concentration.ipynb | 302 ------------------ .../boundary_conditions/concentration.md | 169 ++++++++++ book/content/boundary_conditions/fluxes.ipynb | 42 --- book/content/boundary_conditions/fluxes.md | 162 ++++++++++ .../boundary_conditions/heat_transfer.md | 64 ++++ .../boundary_conditions/surface_reactions.md | 280 ++++++++++++++++ 7 files changed, 677 insertions(+), 344 deletions(-) delete mode 100644 book/content/boundary_conditions/concentration.ipynb create mode 100644 book/content/boundary_conditions/concentration.md delete mode 100644 book/content/boundary_conditions/fluxes.ipynb create mode 100644 book/content/boundary_conditions/fluxes.md create mode 100644 book/content/boundary_conditions/heat_transfer.md create mode 100644 book/content/boundary_conditions/surface_reactions.md diff --git a/book/_toc.yml b/book/_toc.yml index 40cd729b..30f4a6de 100644 --- a/book/_toc.yml +++ b/book/_toc.yml @@ -25,6 +25,8 @@ parts: chapters: - file: content/boundary_conditions/concentration - file: content/boundary_conditions/fluxes + - file: content/boundary_conditions/surface_reactions + - file: content/boundary_conditions/heat_transfer - caption: Species & Reactions chapters: diff --git a/book/content/boundary_conditions/concentration.ipynb b/book/content/boundary_conditions/concentration.ipynb deleted file mode 100644 index 7109ce00..00000000 --- a/book/content/boundary_conditions/concentration.ipynb +++ /dev/null @@ -1,302 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "b24da17d", - "metadata": {}, - "source": [ - "# Concentration BCs #\n", - "\n", - "Boundary conditions (BCs) are essential to FESTIM simulations, as they describe the mathematical problem at the boundaries of the simulated domain. Read more about BCs _[here](https://festim.readthedocs.io/en/fenicsx/userguide/boundary_conditions.html)_.\n", - "\n", - "This tutorial goes over how to define concentration boundary conditions for hydrogen transport simulations.\n", - "\n", - "Objectives:\n", - "* Define a fixed concentration BC\n", - "* Choose which solubility law\n", - "* Plasma implantation approximation" - ] - }, - { - "cell_type": "markdown", - "id": "f018139f", - "metadata": {}, - "source": [ - "## Defining fixed concentration ##\n", - "\n", - "BCs need to be assigned to surfaces using FESTIM's `SurfaceSubdomain` class. To define the concentration of a specific species, we can use `FixedConcentrationBC`:" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "id": "c9b8fdc2", - "metadata": {}, - "outputs": [], - "source": [ - "from festim import FixedConcentrationBC, Species, SurfaceSubdomain\n", - "\n", - "boundary = SurfaceSubdomain(id=1)\n", - "H = Species(name=\"Hydrogen\")\n", - "\n", - "my_bc = FixedConcentrationBC(subdomain=boundary, value=10, species=H)" - ] - }, - { - "cell_type": "markdown", - "id": "0082375d", - "metadata": {}, - "source": [ - "The imposed concentration can be dependent on space, time and temperature, such as:\n", - "\n", - "$$ \n", - "\n", - "BC = 10 + x^2 + T(t+2)\n", - "\n", - "$$" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "7f47ff10", - "metadata": {}, - "outputs": [], - "source": [ - "my_custom_value = lambda x, t, T: 10 + x[0]**2 + T *(t + 2)\n", - "\n", - "my_bc = FixedConcentrationBC(subdomain=boundary, value=my_custom_value, species=H)" - ] - }, - { - "cell_type": "markdown", - "id": "b55d56c7", - "metadata": {}, - "source": [ - "Users can define the surface concentration using either Sieverts’ law, $c = S(T)\\sqrt P$, or Henry's law, $c=K_H P$, where $S(T)$ and $K_H$ denote the temperature-dependent Sieverts’ and Henry’s solubility coefficients, respectively, and $P$ is the partial pressure of the species on the surface. \n", - "\n", - "For Sieverts' law of solubility, we can use `festim.SievertsBC`:" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "8c72893c", - "metadata": {}, - "outputs": [], - "source": [ - "from festim import SievertsBC, SurfaceSubdomain, Species\n", - "\n", - "boundary = SurfaceSubdomain(id=1)\n", - "H = Species(name=\"Hydrogen\")\n", - "\n", - "custom_pressure_value = lambda t: 2 + t\n", - "my_bc = SievertsBC(subdomain=3, S_0=2, E_S=0.1, species=H, pressure=custom_pressure_value)" - ] - }, - { - "cell_type": "markdown", - "id": "53ab8ee2", - "metadata": {}, - "source": [ - "Similarly, for Henry's law of solubility, we can use `festim.HenrysBC`:" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "id": "38b13107", - "metadata": {}, - "outputs": [], - "source": [ - "from festim import HenrysBC\n", - "\n", - "pressure_value = lambda t: 5 * t\n", - "my_bc = HenrysBC(subdomain=3, H_0=1.5, E_H=0.2, species=H, pressure=pressure_value)" - ] - }, - { - "cell_type": "markdown", - "id": "1f18a9f7", - "metadata": {}, - "source": [ - "## Plasma implantation ##" - ] - }, - { - "cell_type": "markdown", - "id": "5335a5cb", - "metadata": {}, - "source": [ - "We can also approximate plasma implantation using FESTIM's `ParticleSource` class, which is helpful in modeling thermo-desorption spectra (TDS) experiments or including the effect of plasma exposure on hydrogen transport. Learn more about how FESTIM approximates plasma implantation _[here](https://festim.readthedocs.io/en/fenicsx/theory.html)_.\n", - "\n", - "We can define a representative source term:\n", - "\n", - "$$ S_{ext} = \\varphi \\cdot f(x) $$\n", - "\n", - "where $\\varphi$ is the implantation flux and $f(x)$ is a Gaussian spatial distribution. For example, if we had an implantation flux $\\varphi = 1 \\mathrm{m}^{-2}\\mathrm{s}^{-1}$ and a distribution mean value of 1 $\\text{nm}$ and width of 3 $\\text{nm}$:" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "id": "6e467cbc", - "metadata": {}, - "outputs": [], - "source": [ - "import festim as F\n", - "import ufl\n", - "import numpy as np\n", - "\n", - "my_model = F.HydrogenTransportProblem()\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a43e3993", - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "ld: warning: duplicate -rpath '/Users/ckhurana/miniconda3/envs/festim-workshop/lib' ignored\n", - "ld: warning: duplicate -rpath '/Users/ckhurana/miniconda3/envs/festim-workshop/lib' ignored\n", - "ld: warning: duplicate -rpath '/Users/ckhurana/miniconda3/envs/festim-workshop/lib' ignored\n", - "ld: warning: duplicate -rpath '/Users/ckhurana/miniconda3/envs/festim-workshop/lib' ignored\n" - ] - } - ], - "source": [ - "\n", - "vertices = np.concatenate(\n", - " [\n", - " np.linspace(0, 30e-9, num=200),\n", - " np.linspace(30e-9, 3e-6, num=300),\n", - " np.linspace(3e-6, 20e-6, num=200),\n", - " ]\n", - ")\n", - "\n", - "my_model.mesh = F.Mesh1D(vertices)\n", - "\n", - "mat = F.Material(D_0=1e-7, E_D=0.2)\n", - "\n", - "volume_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 20e-6], material=mat)\n", - "left_boundary = F.SurfaceSubdomain1D(id=1, x=0)\n", - "right_boundary = F.SurfaceSubdomain1D(id=2, x=20e-6)\n", - "incident_flux = 1e12 # H/m2/s\n", - "\n", - "my_model.subdomains = [\n", - " volume_subdomain,\n", - " left_boundary,\n", - " right_boundary,\n", - "]\n", - "def gaussian_distribution(x, center, width):\n", - " return (\n", - " 1\n", - " / (width * (2 * ufl.pi) ** 0.5)\n", - " * ufl.exp(-0.5 * ((x[0] - center) / width) ** 2)\n", - " )\n", - "H = F.Species(\"H\")\n", - "my_model.species = [H]\n", - "\n", - "source_term = F.ParticleSource(\n", - " value=lambda x: incident_flux * gaussian_distribution(x, 10e-6, 20e-6),\n", - " volume=volume_subdomain,\n", - " species=H,\n", - ")\n", - "\n", - "my_model.sources = [source_term]\n", - "\n", - "my_model.boundary_conditions = [\n", - " F.FixedConcentrationBC(subdomain=left_boundary, value=0, species=H),\n", - " F.FixedConcentrationBC(subdomain=right_boundary, value=0, species=H),\n", - "]\n", - "\n", - "my_model.temperature = 400\n", - "my_model.settings = F.Settings(atol=1e10, rtol=1e-10, transient=False)\n", - "\n", - "profile_export = F.Profile1DExport(field=H,subdomain=volume_subdomain)\n", - "my_model.exports = [profile_export]\n", - "my_model.initialise()\n", - "my_model.run()\n", - "# " - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "e966568c", - "metadata": {}, - "outputs": [], - "source": [ - "x = my_model.exports[0].x\n", - "c = my_model.exports[0].data " - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "id": "abbb180f", - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGvCAYAAABW/q+QAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAI0NJREFUeJzt3X1UlGUC9/HfxMugJZOKgiQqtiaSbkcxEYu1TopiWe66qbmRdcyNbc1FTiff2ieyjqRb5rr4srqYdU6p2yLmORkrHRUt0dIFLSO3jJJVJsJVhrQFxfv5o8d5mhjQQUaZq+/nnPmDe677nvtquA7f7nnRZlmWJQAAAINcc7VPAAAAoLUROAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACME3y1T+BqOH/+vI4fP64OHTrIZrNd7dMBAACXwLIs1dbWKjo6Wtdc0/w1mp9k4Bw/flwxMTFX+zQAAEALVFRUqHv37s2O+UkGTocOHSR9/x8oPDz8Kp8NAAC4FC6XSzExMe6/4835SQbOhZelwsPDCRwAAALMpby9hDcZAwAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADDOFQmc5cuXKzY2VmFhYUpISNCuXbuaHV9UVKSEhASFhYWpd+/eWrlyZZNj169fL5vNpnHjxrXyWQMAgEDl98DZsGGDMjIyNG/ePJWUlCg5OVmpqak6evSo1/Hl5eUaM2aMkpOTVVJSorlz52rGjBnKy8trNParr77Sk08+qeTkZH9PAwAABBCbZVmWPx8gMTFRgwYN0ooVK9zb+vXrp3Hjxik7O7vR+FmzZmnz5s0qKytzb0tPT9eBAwdUXFzs3tbQ0KDhw4frkUce0a5du3Tq1Clt2rTpks7J5XLJ4XCopqZG4eHhLZ8cAAC4Ynz5++3XKzj19fXav3+/UlJSPLanpKRo9+7dXvcpLi5uNH7UqFHat2+fzp496942f/58denSRVOnTr3oedTV1cnlcnncAACAufwaONXV1WpoaFBkZKTH9sjISDmdTq/7OJ1Or+PPnTun6upqSdL777+v3NxcrV69+pLOIzs7Ww6Hw32LiYlpwWwAAECguCJvMrbZbB4/W5bVaNvFxl/YXltbqwcffFCrV69WRETEJT3+nDlzVFNT475VVFT4OAMAABBIgv158IiICAUFBTW6WlNVVdXoKs0FUVFRXscHBwerc+fOOnTokL788kuNHTvWff/58+clScHBwTp8+LBuvPFGj/3tdrvsdntrTAkAAAQAv17BCQ0NVUJCggoLCz22FxYWatiwYV73SUpKajR+69atGjx4sEJCQhQXF6ePPvpIpaWl7tu9996rO++8U6Wlpbz8BAAA/HsFR5IyMzOVlpamwYMHKykpSatWrdLRo0eVnp4u6fuXj44dO6bXXntN0vefmMrJyVFmZqamTZum4uJi5ebmat26dZKksLAw9e/f3+Mxrr/+eklqtB0AAPw0+T1wJk6cqBMnTmj+/PmqrKxU//79tWXLFvXs2VOSVFlZ6fGdOLGxsdqyZYtmzpypZcuWKTo6WkuXLtX48eP9faoAAMAQfv8enLaI78EBACDwtJnvwQEAALgaCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxrkigbN8+XLFxsYqLCxMCQkJ2rVrV7Pji4qKlJCQoLCwMPXu3VsrV670uH/16tVKTk5Wx44d1bFjR40YMUIffPCBP6cAAAACiN8DZ8OGDcrIyNC8efNUUlKi5ORkpaam6ujRo17Hl5eXa8yYMUpOTlZJSYnmzp2rGTNmKC8vzz1mx44deuCBB7R9+3YVFxerR48eSklJ0bFjx/w9HQAAEABslmVZ/nyAxMREDRo0SCtWrHBv69evn8aNG6fs7OxG42fNmqXNmzerrKzMvS09PV0HDhxQcXGx18doaGhQx44dlZOTo4ceeuii5+RyueRwOFRTU6Pw8PAWzAoAAFxpvvz99usVnPr6eu3fv18pKSke21NSUrR7926v+xQXFzcaP2rUKO3bt09nz571us+ZM2d09uxZderUyev9dXV1crlcHjcAAGAuvwZOdXW1GhoaFBkZ6bE9MjJSTqfT6z5Op9Pr+HPnzqm6utrrPrNnz9YNN9ygESNGeL0/OztbDofDfYuJiWnBbAAAQKC4Im8yttlsHj9bltVo28XGe9suSYsWLdK6deu0ceNGhYWFeT3enDlzVFNT475VVFT4OgUAABBAgv158IiICAUFBTW6WlNVVdXoKs0FUVFRXscHBwerc+fOHttffPFFLViwQO+++65+/vOfN3kedrtddru9hbMAAACBxq9XcEJDQ5WQkKDCwkKP7YWFhRo2bJjXfZKSkhqN37p1qwYPHqyQkBD3tj/96U967rnnVFBQoMGDB7f+yQMAgIDl95eoMjMz9be//U1r1qxRWVmZZs6cqaNHjyo9PV3S9y8f/fCTT+np6frqq6+UmZmpsrIyrVmzRrm5uXryySfdYxYtWqSnn35aa9asUa9eveR0OuV0OvXtt9/6ezoAACAA+PUlKkmaOHGiTpw4ofnz56uyslL9+/fXli1b1LNnT0lSZWWlx3fixMbGasuWLZo5c6aWLVum6OhoLV26VOPHj3ePWb58uerr6/XrX//a47GeeeYZZWVl+XtKAACgjfP79+C0RXwPDgAAgafNfA8OAADA1UDgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADDOFQmc5cuXKzY2VmFhYUpISNCuXbuaHV9UVKSEhASFhYWpd+/eWrlyZaMxeXl5io+Pl91uV3x8vPLz8/11+gAAIMD4PXA2bNigjIwMzZs3TyUlJUpOTlZqaqqOHj3qdXx5ebnGjBmj5ORklZSUaO7cuZoxY4by8vLcY4qLizVx4kSlpaXpwIEDSktL04QJE7R3715/TwcAAAQAm2VZlj8fIDExUYMGDdKKFSvc2/r166dx48YpOzu70fhZs2Zp8+bNKisrc29LT0/XgQMHVFxcLEmaOHGiXC6X3nnnHfeY0aNHq2PHjlq3bt1Fz8nlcsnhcKimpkbh4eGXMz0PlmXpu7MNrXY8AAACWbuQINlstlY7ni9/v4Nb7VG9qK+v1/79+zV79myP7SkpKdq9e7fXfYqLi5WSkuKxbdSoUcrNzdXZs2cVEhKi4uJizZw5s9GYJUuWeD1mXV2d6urq3D+7XK4WzObivjvboPj/80+/HBsAgEDzyfxRah/q19Rokl9foqqurlZDQ4MiIyM9tkdGRsrpdHrdx+l0eh1/7tw5VVdXNzumqWNmZ2fL4XC4bzExMS2dEgAACABXJKt+fHnKsqxmL1l5G//j7b4cc86cOcrMzHT/7HK5/BI57UKC9Mn8Ua1+XAAAAlG7kKCr9th+DZyIiAgFBQU1urJSVVXV6ArMBVFRUV7HBwcHq3Pnzs2OaeqYdrtddru9pdO4ZDab7apdigMAAP+fX1+iCg0NVUJCggoLCz22FxYWatiwYV73SUpKajR+69atGjx4sEJCQpod09QxAQDAT4vfLzdkZmYqLS1NgwcPVlJSklatWqWjR48qPT1d0vcvHx07dkyvvfaapO8/MZWTk6PMzExNmzZNxcXFys3N9fh01B/+8Af94he/0MKFC3Xffffprbfe0rvvvqv33nvP39MBAAABwO+BM3HiRJ04cULz589XZWWl+vfvry1btqhnz56SpMrKSo/vxImNjdWWLVs0c+ZMLVu2TNHR0Vq6dKnGjx/vHjNs2DCtX79eTz/9tP74xz/qxhtv1IYNG5SYmOjv6QAAgADg9+/BaYv89T04AADAf3z5+82/RQUAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACM49fAOXnypNLS0uRwOORwOJSWlqZTp041u49lWcrKylJ0dLTatWunO+64Q4cOHXLf/9///ldPPPGE+vbtq/bt26tHjx6aMWOGampq/DkVAAAQQPwaOJMnT1ZpaakKCgpUUFCg0tJSpaWlNbvPokWLtHjxYuXk5OjDDz9UVFSURo4cqdraWknS8ePHdfz4cb344ov66KOPtHbtWhUUFGjq1Kn+nAoAAAggNsuyLH8cuKysTPHx8dqzZ48SExMlSXv27FFSUpI+/fRT9e3bt9E+lmUpOjpaGRkZmjVrliSprq5OkZGRWrhwoR577DGvj/Xmm2/qwQcf1OnTpxUcHHzRc3O5XHI4HKqpqVF4ePhlzBIAAFwpvvz99tsVnOLiYjkcDnfcSNLQoUPlcDi0e/dur/uUl5fL6XQqJSXFvc1ut2v48OFN7iPJPdGm4qaurk4ul8vjBgAAzOW3wHE6neratWuj7V27dpXT6WxyH0mKjIz02B4ZGdnkPidOnNBzzz3X5NUdScrOzna/D8jhcCgmJuZSpwEAAAKQz4GTlZUlm83W7G3fvn2SJJvN1mh/y7K8bv+hH9/f1D4ul0t333234uPj9cwzzzR5vDlz5qimpsZ9q6iouJSpAgCAAHXxN6z8yPTp0zVp0qRmx/Tq1UsHDx7U119/3ei+b775ptEVmguioqIkfX8lp1u3bu7tVVVVjfapra3V6NGjdd111yk/P18hISFNno/dbpfdbm/2nAEAgDl8DpyIiAhFRERcdFxSUpJqamr0wQcfaMiQIZKkvXv3qqamRsOGDfO6T2xsrKKiolRYWKiBAwdKkurr61VUVKSFCxe6x7lcLo0aNUp2u12bN29WWFiYr9MAAAAG89t7cPr166fRo0dr2rRp2rNnj/bs2aNp06bpnnvu8fgEVVxcnPLz8yV9/9JURkaGFixYoPz8fH388cd6+OGH1b59e02ePFnS91duUlJSdPr0aeXm5srlcsnpdMrpdKqhocFf0wEAAAHE5ys4vnj99dc1Y8YM96ei7r33XuXk5HiMOXz4sMeX9D311FP67rvv9Pjjj+vkyZNKTEzU1q1b1aFDB0nS/v37tXfvXknSz372M49jlZeXq1evXn6cEQAACAR++x6ctozvwQEAIPC0ie/BAQAAuFoIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADG8WvgnDx5UmlpaXI4HHI4HEpLS9OpU6ea3ceyLGVlZSk6Olrt2rXTHXfcoUOHDjU5NjU1VTabTZs2bWr9CQAAgIDk18CZPHmySktLVVBQoIKCApWWliotLa3ZfRYtWqTFixcrJydHH374oaKiojRy5EjV1tY2GrtkyRLZbDZ/nT4AAAhQwf46cFlZmQoKCrRnzx4lJiZKklavXq2kpCQdPnxYffv2bbSPZVlasmSJ5s2bp1/96leSpFdffVWRkZF644039Nhjj7nHHjhwQIsXL9aHH36obt26+WsaAAAgAPntCk5xcbEcDoc7biRp6NChcjgc2r17t9d9ysvL5XQ6lZKS4t5mt9s1fPhwj33OnDmjBx54QDk5OYqKirroudTV1cnlcnncAACAufwWOE6nU127dm20vWvXrnI6nU3uI0mRkZEe2yMjIz32mTlzpoYNG6b77rvvks4lOzvb/T4gh8OhmJiYS50GAAAIQD4HTlZWlmw2W7O3ffv2SZLX98dYlnXR9838+P4f7rN582Zt27ZNS5YsueRznjNnjmpqaty3ioqKS94XAAAEHp/fgzN9+nRNmjSp2TG9evXSwYMH9fXXXze675tvvml0heaCCy83OZ1Oj/fVVFVVuffZtm2bjhw5ouuvv95j3/Hjxys5OVk7duxodFy73S673d7sOQMAAHP4HDgRERGKiIi46LikpCTV1NTogw8+0JAhQyRJe/fuVU1NjYYNG+Z1n9jYWEVFRamwsFADBw6UJNXX16uoqEgLFy6UJM2ePVuPPvqox34DBgzQyy+/rLFjx/o6HQAAYCC/fYqqX79+Gj16tKZNm6a//vWvkqTf/va3uueeezw+QRUXF6fs7Gz98pe/lM1mU0ZGhhYsWKA+ffqoT58+WrBggdq3b6/JkydL+v4qj7c3Fvfo0UOxsbH+mg4AAAggfgscSXr99dc1Y8YM96ei7r33XuXk5HiMOXz4sGpqatw/P/XUU/ruu+/0+OOP6+TJk0pMTNTWrVvVoUMHf54qAAAwiM2yLOtqn8SV5nK55HA4VFNTo/Dw8Kt9OgAA4BL48vebf4sKAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYhcAAAgHEIHAAAYBwCBwAAGIfAAQAAxiFwAACAcQgcAABgHAIHAAAYh8ABAADGIXAAAIBxCBwAAGAcAgcAABiHwAEAAMYJvtoncDVYliVJcrlcV/lMAADApbrwd/vC3/Hm/CQDp7a2VpIUExNzlc8EAAD4qra2Vg6Ho9kxNutSMsgw58+f1/Hjx9WhQwfZbLZWPbbL5VJMTIwqKioUHh7eqsduC0yfn2T+HJlf4DN9jswv8PlrjpZlqba2VtHR0brmmubfZfOTvIJzzTXXqHv37n59jPDwcGN/cSXz5yeZP0fmF/hMnyPzC3z+mOPFrtxcwJuMAQCAcQgcAABgHAKnldntdj3zzDOy2+1X+1T8wvT5SebPkfkFPtPnyPwCX1uY40/yTcYAAMBsXMEBAADGIXAAAIBxCBwAAGAcAgcAABiHwLmI5cuXKzY2VmFhYUpISNCuXbuaHV9UVKSEhASFhYWpd+/eWrlyZaMxeXl5io+Pl91uV3x8vPLz8/11+pfElzlu3LhRI0eOVJcuXRQeHq6kpCT985//9Bizdu1a2Wy2Rrf//e9//p6KV77Mb8eOHV7P/dNPP/UY15aeQ1/m9/DDD3ud38033+we05aev507d2rs2LGKjo6WzWbTpk2bLrpPoK1BX+cYaGvQ1/kF4hr0dY6BtA6zs7N16623qkOHDuratavGjRunw4cPX3S/trAOCZxmbNiwQRkZGZo3b55KSkqUnJys1NRUHT161Ov48vJyjRkzRsnJySopKdHcuXM1Y8YM5eXluccUFxdr4sSJSktL04EDB5SWlqYJEyZo7969V2paHnyd486dOzVy5Eht2bJF+/fv15133qmxY8eqpKTEY1x4eLgqKys9bmFhYVdiSh58nd8Fhw8f9jj3Pn36uO9rS8+hr/P785//7DGviooKderUSffff7/HuLby/J0+fVq33HKLcnJyLml8IK5BX+cYaGvQ1/ldEChrUPJ9joG0DouKivT73/9ee/bsUWFhoc6dO6eUlBSdPn26yX3azDq00KQhQ4ZY6enpHtvi4uKs2bNnex3/1FNPWXFxcR7bHnvsMWvo0KHunydMmGCNHj3aY8yoUaOsSZMmtdJZ+8bXOXoTHx9vPfvss+6fX3nlFcvhcLTWKV4WX+e3fft2S5J18uTJJo/Zlp7Dy33+8vPzLZvNZn355ZfubW3p+fshSVZ+fn6zYwJxDf7QpczRm7a8Bn/oUuYXaGvwx1ryHAbSOqyqqrIkWUVFRU2OaSvrkCs4Taivr9f+/fuVkpLisT0lJUW7d+/2uk9xcXGj8aNGjdK+fft09uzZZsc0dUx/askcf+z8+fOqra1Vp06dPLZ/++236tmzp7p376577rmn0f9dXgmXM7+BAweqW7duuuuuu7R9+3aP+9rKc9gaz19ubq5GjBihnj17emxvC89fSwTaGmwNbXkNXo5AWIOtJZDWYU1NjSQ1+n37obayDgmcJlRXV6uhoUGRkZEe2yMjI+V0Or3u43Q6vY4/d+6cqqurmx3T1DH9qSVz/LGXXnpJp0+f1oQJE9zb4uLitHbtWm3evFnr1q1TWFiYbrvtNn322Wetev4X05L5devWTatWrVJeXp42btyovn376q677tLOnTvdY9rKc3i5z19lZaXeeecdPfroox7b28rz1xKBtgZbQ1tegy0RSGuwNQTSOrQsS5mZmbr99tvVv3//Jse1lXX4k/zXxH1hs9k8frYsq9G2i43/8XZfj+lvLT2fdevWKSsrS2+99Za6du3q3j506FANHTrU/fNtt92mQYMG6S9/+YuWLl3aeid+iXyZX9++fdW3b1/3z0lJSaqoqNCLL76oX/ziFy06pr+19FzWrl2r66+/XuPGjfPY3taeP18F4hpsqUBZg74IxDV4OQJpHU6fPl0HDx7Ue++9d9GxbWEdcgWnCREREQoKCmpUk1VVVY2q84KoqCiv44ODg9W5c+dmxzR1TH9qyRwv2LBhg6ZOnaq///3vGjFiRLNjr7nmGt16661X/P88Lmd+PzR06FCPc28rz+HlzM+yLK1Zs0ZpaWkKDQ1tduzVev5aItDW4OUIhDXYWtrqGrxcgbQOn3jiCW3evFnbt29X9+7dmx3bVtYhgdOE0NBQJSQkqLCw0GN7YWGhhg0b5nWfpKSkRuO3bt2qwYMHKyQkpNkxTR3Tn1oyR+n7/2t8+OGH9cYbb+juu+++6ONYlqXS0lJ169btss/ZFy2d34+VlJR4nHtbeQ4vZ35FRUX6/PPPNXXq1Is+ztV6/loi0NZgSwXKGmwtbXUNXq5AWIeWZWn69OnauHGjtm3bptjY2Ivu02bWYau9XdlA69evt0JCQqzc3Fzrk08+sTIyMqxrr73W/U732bNnW2lpae7xX3zxhdW+fXtr5syZ1ieffGLl5uZaISEh1j/+8Q/3mPfff98KCgqyXnjhBausrMx64YUXrODgYGvPnj1XfH6W5fsc33jjDSs4ONhatmyZVVlZ6b6dOnXKPSYrK8sqKCiwjhw5YpWUlFiPPPKIFRwcbO3du7fNz+/ll1+28vPzrX//+9/Wxx9/bM2ePduSZOXl5bnHtKXn0Nf5XfDggw9aiYmJXo/Zlp6/2tpaq6SkxCopKbEkWYsXL7ZKSkqsr776yrIsM9agr3MMtDXo6/wCbQ1alu9zvCAQ1uHvfvc7y+FwWDt27PD4fTtz5ox7TFtdhwTORSxbtszq2bOnFRoaag0aNMjjo3FTpkyxhg8f7jF+x44d1sCBA63Q0FCrV69e1ooVKxod880337T69u1rhYSEWHFxcR4L92rwZY7Dhw+3JDW6TZkyxT0mIyPD6tGjhxUaGmp16dLFSklJsXbv3n0FZ+TJl/ktXLjQuvHGG62wsDCrY8eO1u233269/fbbjY7Zlp5DX39HT506ZbVr185atWqV1+O1pefvwkeGm/p9M2EN+jrHQFuDvs4vENdgS35PA2UdepuXJOuVV15xj2mr69D2/yYAAABgDN6DAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4BA4AAGjSzp07NXbsWEVHR8tms2nTpk1+fbysrCzZbDaPW1RUlM/HIXAAAECTTp8+rVtuuUU5OTlX7DFvvvlmVVZWum8fffSRz8cI9sN5AQAAQ6Smpio1NbXJ++vr6/X000/r9ddf16lTp9S/f38tXLhQd9xxR4sfMzg4uEVXbX6IKzgAAKDFHnnkEb3//vtav369Dh48qPvvv1+jR4/WZ5991uJjfvbZZ4qOjlZsbKwmTZqkL774wudj8G9RAQCAS2Kz2ZSfn69x48ZJko4cOaI+ffroP//5j6Kjo93jRowYoSFDhmjBggU+P8Y777yjM2fO6KabbtLXX3+t559/Xp9++qkOHTqkzp07X/JxuIIDAABa5F//+pcsy9JNN92k6667zn0rKirSkSNHJElffvllozcN//g2ffp09zFTU1M1fvx4DRgwQCNGjNDbb78tSXr11Vd9OjfegwMAAFrk/PnzCgoK0v79+xUUFORx33XXXSdJuuGGG1RWVtbscTp27Njkfddee60GDBjg80teBA4AAGiRgQMHqqGhQVVVVUpOTvY6JiQkRHFxcS1+jLq6OpWVlTV5/KYQOAAAoEnffvutPv/8c/fP5eXlKi0tVadOnXTTTTfpN7/5jR566CG99NJLGjhwoKqrq7Vt2zYNGDBAY8aM8fnxnnzySY0dO1Y9evRQVVWVnn/+eblcLk2ZMsWn4/AmYwAA0KQdO3bozjvvbLR9ypQpWrt2rc6ePavnn39er732mo4dO6bOnTsrKSlJzz77rAYMGODz402aNEk7d+5UdXW1unTpoqFDh+q5555TfHy8T8chcAAAgHH4FBUAADAOgQMAAIxD4AAAAOMQOAAAwDgEDgAAMA6BAwAAjEPgAAAA4xA4AADAOAQOAAAwDoEDAACMQ+AAAADjEDgAAMA4/xc5I//Q/+oy9QAAAABJRU5ErkJggg==", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "import matplotlib.pyplot as plt\n", - "\n", - "c = np.array(c)\n", - "\n", - "plt.plot(x, c.squeeze())\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "13b71fb6", - "metadata": {}, - "source": [] - }, - { - "cell_type": "code", - "execution_count": 9, - "id": "7579b4a1", - "metadata": {}, - "outputs": [], - "source": [ - "c_flat = c[0, 0, :]\n" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "festim-workshop", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.18" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/book/content/boundary_conditions/concentration.md b/book/content/boundary_conditions/concentration.md new file mode 100644 index 00000000..c1538d5d --- /dev/null +++ b/book/content/boundary_conditions/concentration.md @@ -0,0 +1,169 @@ +--- +jupytext: + formats: ipynb,md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.17.3 +kernelspec: + display_name: festim-workshop + language: python + name: python3 +--- + +# Concentration # + +Boundary conditions (BCs) are essential to FESTIM simulations, as they describe the mathematical problem at the boundaries of the simulated domain. Read more about BCs _[here](https://festim.readthedocs.io/en/fenicsx/userguide/boundary_conditions.html)_. + +This tutorial goes over how to define concentration boundary conditions for hydrogen transport simulations. + +Objectives: +* Define a fixed concentration BC +* Choose which solubility law (Sieverts' or Henry's) +* Solve a hydrogen transport problem with plasma implantation + ++++ + +## Defining fixed concentration ## + +BCs need to be assigned to surfaces using FESTIM's `SurfaceSubdomain` class. To define the concentration of a specific species, we can use `FixedConcentrationBC`: + +```{code-cell} ipython3 +from festim import FixedConcentrationBC, Species, SurfaceSubdomain + +boundary = SurfaceSubdomain(id=1) +H = Species(name="Hydrogen") + +my_bc = FixedConcentrationBC(subdomain=boundary, value=10, species=H) +``` + +The imposed concentration can be dependent on space, time and temperature, such as: + +$$ + +BC = 10 + x^2 + T(t+2) + +$$ + +```{code-cell} ipython3 +my_custom_value = lambda x, t, T: 10 + x[0]**2 + T *(t + 2) + +my_bc = FixedConcentrationBC(subdomain=boundary, value=my_custom_value, species=H) +``` + +## Choosing a solubility law ## + +Users can define the surface concentration using either Sieverts’ law, $c = S(T)\sqrt P$, or Henry's law, $c=K_H P$, where $S(T)$ and $K_H$ denote the temperature-dependent Sieverts’ and Henry’s solubility coefficients, respectively, and $P$ is the partial pressure of the species on the surface. + +For Sieverts' law of solubility, we can use `festim.SievertsBC`: + +```{code-cell} ipython3 +from festim import SievertsBC, SurfaceSubdomain, Species + +boundary = SurfaceSubdomain(id=1) +H = Species(name="Hydrogen") + +custom_pressure_value = lambda t: 2 + t +my_bc = SievertsBC(subdomain=3, S_0=2, E_S=0.1, species=H, pressure=custom_pressure_value) +``` + +Similarly, for Henry's law of solubility, we can use `festim.HenrysBC`: + +```{code-cell} ipython3 +from festim import HenrysBC + +pressure_value = lambda t: 5 * t +my_bc = HenrysBC(subdomain=3, H_0=1.5, E_H=0.2, species=H, pressure=pressure_value) +``` + +## Plasma implantation approximation ## + ++++ + +We can also approximate plasma implantation using FESTIM's `ParticleSource` class, which is helpful in modeling thermo-desorption spectra (TDS) experiments or including the effect of plasma exposure on hydrogen transport. Learn more about how FESTIM approximates plasma implantation _[here](https://festim.readthedocs.io/en/fenicsx/theory.html)_. + +Consider the following 1D plasma implantation problem, where we represent the plasma as a hydrogen source $S_{ext}$: + +$$ S_{ext} = \varphi \cdot f(x) $$ + +$$\varphi = 1\cdot 10^{13} \quad \mathrm{m}^{-2}\mathrm{s}^{-1}$$ + +where $\varphi$ is the implantation flux and $f(x)$ is a Gaussian spatial distribution (distribution mean value of 0.5 $\text{m}$ and width of 1 $\text{m}$). + +First, we setup a 1D mesh ranging from $ [0,1] $ and assign the subdomains and material: + +```{code-cell} ipython3 +import festim as F +import ufl +import numpy as np + +my_model = F.HydrogenTransportProblem() +vertices = np.linspace(0,1,2000) +my_model.mesh = F.Mesh1D(vertices) + +mat = F.Material(D_0=0.1, E_D=0.01) + +volume_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat) +left_boundary = F.SurfaceSubdomain1D(id=1, x=0) +right_boundary = F.SurfaceSubdomain1D(id=2, x=1) + +my_model.subdomains = [ + volume_subdomain, + left_boundary, + right_boundary, +] +``` + +Now, we define our `incident_flux` and `gaussian_distribution` function. We can use `ParticleSource` to represent the source term, and then add it to our model: + +```{code-cell} ipython3 +incident_flux = 1e13 + +def gaussian_distribution(x, center, width): + return ( + 1 + / (width * (2 * ufl.pi) ** 0.5) + * ufl.exp(-0.5 * ((x[0] - center) / width) ** 2) + ) +H = F.Species("H") +my_model.species = [H] + +source_term = F.ParticleSource( + value=lambda x: incident_flux * gaussian_distribution(x, .5, 1), + volume=volume_subdomain, + species=H, +) + +my_model.sources = [source_term] +``` + +Finally, we assign boundary conditions (zero concentration at $x=0$ and $x=1$) and solve our problem: + +```{code-cell} ipython3 +my_model.boundary_conditions = [ + F.FixedConcentrationBC(subdomain=left_boundary, value=0, species=H), + F.FixedConcentrationBC(subdomain=right_boundary, value=0, species=H), +] + +my_model.temperature = 400 +my_model.settings = F.Settings(atol=1e10, rtol=1e-10, transient=False) + +profile_export = F.Profile1DExport(field=H,subdomain=volume_subdomain) +my_model.exports = [profile_export] + +my_model.initialise() +my_model.run() +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +import matplotlib.pyplot as plt + +x = my_model.exports[0].x +c = my_model.exports[0].data[0][0] + +plt.plot(x, c) +plt.show() +``` diff --git a/book/content/boundary_conditions/fluxes.ipynb b/book/content/boundary_conditions/fluxes.ipynb deleted file mode 100644 index da26f665..00000000 --- a/book/content/boundary_conditions/fluxes.ipynb +++ /dev/null @@ -1,42 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "f633d832", - "metadata": {}, - "source": [ - "# Flux BCs #\n", - "\n", - "This tutorial goes over how to define flux boundary conditions in FESTIM.\n", - "\n", - "Objectives:\n", - "* Neuman: flux at boundaries\n", - "* Robin (based on concentrations with arbitrary expression): gradient is imposed as a function of the solutiin\n", - "* Surface reactions" - ] - }, - { - "cell_type": "markdown", - "id": "528d839b", - "metadata": {}, - "source": [ - "## Neumann BCs ##\n", - "\n", - "Users can impose Neumann boundary conditions to define the flux at boundaries. To impose a particle flux, we can use the " - ] - }, - { - "cell_type": "markdown", - "id": "7c26a601", - "metadata": {}, - "source": [] - } - ], - "metadata": { - "language_info": { - "name": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/book/content/boundary_conditions/fluxes.md b/book/content/boundary_conditions/fluxes.md new file mode 100644 index 00000000..4f1984b6 --- /dev/null +++ b/book/content/boundary_conditions/fluxes.md @@ -0,0 +1,162 @@ +--- +jupytext: + formats: ipynb,md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.17.3 +kernelspec: + display_name: festim-workshop + language: python + name: python3 +--- + +# Particle flux # + +This tutorial goes over how to define particle flux boundary conditions in FESTIM. + +Objectives: +* Learn how to define fluxes at boundaries +* Impose boundary fluxes as a function of the concentration +* Solve a problem with multi-species flux boundary conditions + ++++ + +## Defining flux at boundaries ## + +Users can impose the particle flux (Neumann BC) at boundaries using `ParticleFluxBC` class: + +```{code-cell} ipython3 +from festim import ParticleFluxBC, Species, SurfaceSubdomain + +boundary = SurfaceSubdomain(id=1) +H = Species(name="Hydrogen") + +my_flux_bc = ParticleFluxBC(subdomain=boundary, value=2, species=H) +``` + +## Defining concentration-dependant fluxes ## + +Similar to the concentration, the flux can be dependent on space, time and temperature. But for particle fluxes, the values can also be dependent on a species’ concentration. + +For example, let's define a hydrogen flux `J` that depends on the hydrogen concentration `c` and time `t`: + +$$ J(c,t) = 10t^2 + 2c $$ + +```{code-cell} ipython3 +from festim import ParticleFluxBC, Species, SurfaceSubdomain + +boundary = SurfaceSubdomain(id=1) +H = Species(name="Hydrogen") + +J = lambda t, c: 10*t**2 + 2*c + +my_flux_bc = ParticleFluxBC( + subdomain=boundary, + value=J, + species=H, + species_dependent_value={"c": H}, +) +``` + +## Multi-species flux boundary conditions ## + ++++ + +Users may also need to impose a flux boundary condition in multi-species problems where the flux depends on the concentrations of multiple species. + +Consider the following example with three species, A, B, and C, where the particle flux boundary condition depends on each species' concentration: + +$$ J(c_A, c_B, c_C) = 2c_A + 3c_b + 4c_C $$ + +We must first define each species using `Species` and then create the dictionary to be passed into `species_dependent_value`. The dictionary maps each argument in the custom flux function to the corresponding defined species: + +```{code-cell} ipython3 +import festim as F + +my_model = F.HydrogenTransportProblem() + +A = F.Species(name="A") +B = F.Species(name="B") +C = F.Species(name="C") +my_model.species = [A, B, C] + +my_custom_value = lambda c_A, c_B, c_C: 2*c_A + 3*c_B + 4*c_C +species_dependent_value = {"c_A": A, "c_B": B, "c_C": C} +``` + +Now, we create our 1D mesh and assign boundary conditions (flux BC on the left). The boundary condition `ParticleFluxBC` must be added for each species: + +```{code-cell} ipython3 +import numpy as np + +my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100)) + +D = 1 +mat = F.Material( + D_0={A: D, B: D, C: D}, + E_D={A: 0.01, B: 0.01, C: 0.01}, +) + +bulk = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat) +left = F.SurfaceSubdomain1D(id=1, x=0) +right = F.SurfaceSubdomain1D(id=2, x=1) + +my_model.subdomains = [bulk, left, right] +my_model.boundary_conditions = [ + F.ParticleFluxBC( + subdomain=left, + value=my_custom_value, + species=A, + species_dependent_value=species_dependent_value, + ), + F.ParticleFluxBC( + subdomain=left, + value=my_custom_value, + species=B, + species_dependent_value=species_dependent_value, + ), + F.ParticleFluxBC( + subdomain=left, + value=my_custom_value, + species=C, + species_dependent_value=species_dependent_value, + ), +] +``` + +```{note} +The diffusivity pre-factor `D_0` and activation energy `E_d` must be defined for each species in `Material`. Learn more about defining multi-species material properties __[here](https://festim-workshop.readthedocs.io/en/festim2/content/material/material_advanced.html#defining-species-dependent-material-properties)__. +``` + ++++ + +Finally, let's solve and plot the solution for each species: + +```{code-cell} ipython3 +my_model.temperature = 300 +my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) + +my_model.initialise() +my_model.run() +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +import matplotlib.pyplot as plt + +def plot_profile(species, **kwargs): + c = species.post_processing_solution.x.array[:] + x = species.post_processing_solution.function_space.mesh.geometry.x[:,0] + return plt.plot(x, c, **kwargs) + +for species in my_model.species: + plot_profile(species, label=species.name) + +plt.xlabel('Position') +plt.ylabel('Concentration') +plt.legend() +plt.show() +``` diff --git a/book/content/boundary_conditions/heat_transfer.md b/book/content/boundary_conditions/heat_transfer.md new file mode 100644 index 00000000..e2b89126 --- /dev/null +++ b/book/content/boundary_conditions/heat_transfer.md @@ -0,0 +1,64 @@ +--- +jupytext: + formats: ipynb,md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.17.3 +--- + +# Heat transfer # + +This tutorial goes over how to define boundary conditions for heat transfer simulations. + +Objectives: +* Define homogenous temperature boundary conditions +* Define heat flux boundary conditions + ++++ + +## Imposing homogenous temperature boundary conditions ## + +The temperature can be imposed on boundaries using `FixedTemperatureBC`: + +```{code-cell} +from festim import FixedTemperatureBC, SurfaceSubdomain + +boundary = SurfaceSubdomain(id=1) +my_bc = FixedTemperatureBC(subdomain=boundary, value=10) +``` + +To define the temperature as space or time dependent, a function can be passed to the `value` argument, such as: + +$$ \text{BC} = 10 + x^2 + t $$ + +```{code-cell} +from festim import FixedTemperatureBC, SurfaceSubdomain + +boundary = SurfaceSubdomain(id=1) +BC = lambda x, t: 10 + x[0]**2 + t + +my_bc = FixedTemperatureBC(subdomain=boundary, value=BC) +``` + +## Imposing heat flux boundary conditions ## + ++++ + +Heat fluxes can be imposed on boundaries using `HeatFluxBC`, which can depend on space, time, and temperature, such as: + +$$ \text{BC} = 2x + 10t + T $$ + +```{code-cell} +from festim import HeatFluxBC, SurfaceSubdomain + +boundary = SurfaceSubdomain(id=1) +BC = lambda x, t, T: 2 * x[0] + 10 * t + T + +my_flux_bc = HeatFluxBC(subdomain=boundary, value=BC) +``` + +```{note} +Read more about heat transfer settings __[here](https://festim-workshop.readthedocs.io/en/festim2/content/temperatures/temperatures_advanced)__. +``` diff --git a/book/content/boundary_conditions/surface_reactions.md b/book/content/boundary_conditions/surface_reactions.md new file mode 100644 index 00000000..15762910 --- /dev/null +++ b/book/content/boundary_conditions/surface_reactions.md @@ -0,0 +1,280 @@ +--- +jupytext: + formats: ipynb,md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.17.3 +kernelspec: + display_name: festim-workshop + language: python + name: python3 +--- + +# Surface reactions # + ++++ + +FESTIM V2.0 allows users to impose surface reactions on boundaries. Learn more about surface reactions __[here](https://festim-workshop.readthedocs.io/en/festim2/content/species_reactions/surface_reactions.html)__. + +Objectives: +* Impose a simple surface reaction boundary condition +* recombination and dissociation +* Multi-isotope + ++++ + +## Simple surface reaction BC ## + +A reaction is defined by specifying the reactants (`reactant`), gas pressure (`gas_pressure`), forward/backward rate constants (`k_r0` and `k_d0`), and rate activation energies (`E_kr` and `E_kd`). + +```{note} +Reaction rates are given as Arhennius laws. +``` + +To impose the folllowing surface reaction: + +$$ A + B \xrightarrow{k} C $$ + +```{code-cell} ipython3 +from festim import Species, SurfaceReactionBC, SurfaceSubdomain + +boundary = SurfaceSubdomain(id=1) +A = Species("A") +B = Species("B") +C = Species("C") + +my_bc = SurfaceReactionBC( + reactant=[A, B], + gas_pressure=1e5, + k_r0=1, + E_kr=0.1, + k_d0=1e-5, + E_kd=0.1, + subdomain=boundary, +) +``` + +## Recombination and Dissociation ## + +Hydrogen recombination/dissociation can also be modelled using `SurfaceReactionBC`, such as this reaction: + +$$ H + H \overset{K_r}{\underset{K_d}{\rightleftharpoons}} H_2 $$ + +```{code-cell} ipython3 +from festim import Species, SurfaceReactionBC, SurfaceSubdomain + +boundary = SurfaceSubdomain(id=1) +H = Species("H") + +my_bc = SurfaceReactionBC( + reactant=[H, H], + gas_pressure=1e5, + k_r0=1, + E_kr=0.1, + k_d0=1e-5, + E_kd=0.1, + subdomain=boundary, +) +``` + +## Multiple isotopes ## + +Surface reactions can involve multiple hydrogen isotopes, enabling the modelling of more complex interactions between species. For example, in a system with both mobile hydrogen and tritium, various molecular recombination pathways may occur at the surface, resulting in the formation of $H_2$, $T_2$, and $HT$: + +```{code-cell} ipython3 +from festim import Species, SurfaceReactionBC, SurfaceSubdomain + +boundary = SurfaceSubdomain(id=1) +H = Species("H") +T = Species("T") + +reac1 = SurfaceReactionBC( + reactant=[H, H], + gas_pressure=1e6, + k_r0=1.0, + E_kr=0.1, + k_d0=0.5, + E_kd=0.1, + subdomain=boundary, +) +reac2 = SurfaceReactionBC( + reactant=[T, T], + gas_pressure=1e6, + k_r0=1.0, + E_kr=0.1, + k_d0=0.5, + E_kd=0.1, + subdomain=boundary, +) +reac3 = SurfaceReactionBC( + reactant=[H, T], + gas_pressure=1e6, + k_r0=1.0, + E_kr=0.1, + k_d0=0.5, + E_kd=0.1, + subdomain=boundary, +) +``` + +## Isotopic exchange ## + +```{code-cell} ipython3 +import dolfinx.fem as fem +import numpy as np +import ufl + +import festim as F + +my_model = F.HydrogenTransportProblem() +my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 1000)) +my_mat = F.Material(name="mat", D_0=1, E_D=0) +vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat) +left = F.SurfaceSubdomain1D(id=1, x=0) +right = F.SurfaceSubdomain1D(id=2, x=1) + +my_model.subdomains = [vol, left, right] + +H = F.Species("H") +D = F.Species("D") +my_model.species = [H, D] + +my_model.temperature = 500 + +surface_reaction_hd = F.SurfaceReactionBC( + reactant=[H, D], + gas_pressure=0, + k_r0=0.01, + E_kr=0, + k_d0=0, + E_kd=0, + subdomain=right, +) + +surface_reaction_hh = F.SurfaceReactionBC( + reactant=[H, H], + gas_pressure=lambda t: ufl.conditional(ufl.gt(t, 1), 2, 0), + k_r0=0.02, + E_kr=0, + k_d0=0.03, + E_kd=0, + subdomain=right, +) + +surface_reaction_dd = F.SurfaceReactionBC( + reactant=[D, D], + gas_pressure=0, + k_r0=0.01, + E_kr=0, + k_d0=0, + E_kd=0, + subdomain=right, +) + +my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=2, species=H), + F.DirichletBC(subdomain=left, value=2, species=D), + surface_reaction_hd, + surface_reaction_hh, + surface_reaction_dd, +] + +H_flux_right = F.SurfaceFlux(H, right) +H_flux_left = F.SurfaceFlux(H, left) +D_flux_right = F.SurfaceFlux(D, right) +D_flux_left = F.SurfaceFlux(D, left) + +my_model.exports = [ + F.XDMFExport("test.xdmf", H), + H_flux_left, + H_flux_right, + D_flux_left, + D_flux_right, + HD_flux, + HH_flux, + DD_flux, +] + +my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=5, transient=True) + +my_model.settings.stepsize = 0.1 + +my_model.initialise() +my_model.run() + +import matplotlib.pyplot as plt + +plt.stackplot( + H_flux_left.t, + np.abs(H_flux_left.data), + np.abs(D_flux_left.data), + labels=["H_in", "D_in"], +) +plt.stackplot( + H_flux_right.t, + -np.abs(H_flux_right.data), + -np.abs(D_flux_right.data), + labels=["H_out", "D_out"], +) +plt.legend() +plt.xlabel("Time (s)") +plt.ylabel("Flux (atom/m^2/s)") +plt.figure() +plt.stackplot( + HD_flux.t, + np.abs(HH_flux.data), + np.abs(HD_flux.data), + np.abs(DD_flux.data), + labels=["HH", "HD", "DD"], +) +plt.legend(reverse=True) +plt.xlabel("Time (s)") +plt.ylabel("Flux (molecule/m^2/s)") + + +plt.figure() +plt.plot(H_flux_right.t, -np.array(H_flux_right.data), label="from gradient (H)") +plt.plot( + H_flux_right.t, + 2 * np.array(HH_flux.data) + np.array(HD_flux.data), + linestyle="--", + label="from reaction rates (H)", +) + +plt.plot(D_flux_right.t, -np.array(D_flux_right.data), label="from gradient (D)") +plt.plot( + D_flux_right.t, + 2 * np.array(DD_flux.data) + np.array(HD_flux.data), + linestyle="--", + label="from reaction rates (D)", +) +plt.xlabel("Time (s)") +plt.ylabel("Flux (atom/m^2/s)") +plt.legend() +plt.show() + +# check that H_flux_right == 2*HH_flux + HD_flux +H_flux_from_gradient = -np.array(H_flux_right.data) +H_flux_from_reac = 2 * np.array(HH_flux.data) + np.array(HD_flux.data) +assert np.allclose( + H_flux_from_gradient, + H_flux_from_reac, + rtol=0.5e-2, + atol=0.005, +) +# check that D_flux_right == 2*DD_flux + HD_flux +D_flux_from_gradient = -np.array(D_flux_right.data) +D_flux_from_reac = 2 * np.array(DD_flux.data) + np.array(HD_flux.data) +assert np.allclose( + D_flux_from_gradient, + D_flux_from_reac, + rtol=0.5e-2, + atol=0.005, +) +``` + +```{code-cell} ipython3 + +``` From 2e37c8b59adc87739c3bfd7e3c0b1b43a9629512 Mon Sep 17 00:00:00 2001 From: Chirag Date: Tue, 28 Oct 2025 12:33:17 -0400 Subject: [PATCH 03/17] added surface reactions and fluxes --- .../boundary_conditions/concentration.md | 2 +- book/content/boundary_conditions/fluxes.md | 71 +++- .../boundary_conditions/heat_transfer.md | 2 +- .../boundary_conditions/surface_reactions.md | 392 ++++++++++-------- 4 files changed, 276 insertions(+), 191 deletions(-) diff --git a/book/content/boundary_conditions/concentration.md b/book/content/boundary_conditions/concentration.md index c1538d5d..ec09ff61 100644 --- a/book/content/boundary_conditions/concentration.md +++ b/book/content/boundary_conditions/concentration.md @@ -5,7 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.17.3 + jupytext_version: 1.18.1 kernelspec: display_name: festim-workshop language: python diff --git a/book/content/boundary_conditions/fluxes.md b/book/content/boundary_conditions/fluxes.md index 4f1984b6..493deddc 100644 --- a/book/content/boundary_conditions/fluxes.md +++ b/book/content/boundary_conditions/fluxes.md @@ -5,7 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.17.3 + jupytext_version: 1.18.1 kernelspec: display_name: festim-workshop language: python @@ -64,13 +64,15 @@ my_flux_bc = ParticleFluxBC( +++ -Users may also need to impose a flux boundary condition in multi-species problems where the flux depends on the concentrations of multiple species. +Users may also need to impose a flux boundary condition which depends on the concentrations of multiple species. -Consider the following example with three species, A, B, and C, where the particle flux boundary condition depends on each species' concentration: +Consider the following 1D example with three species, $\mathrm{A}$, $\mathrm{B}$, and $\mathrm{C}$, with recombination fluxes $\phi_{AB}$ and $\phi_{ABC}$ that depend on the concentrations $\mathrm{c}$: -$$ J(c_A, c_B, c_C) = 2c_A + 3c_b + 4c_C $$ +$$ \phi_{AB} = -c_A c_B $$ -We must first define each species using `Species` and then create the dictionary to be passed into `species_dependent_value`. The dictionary maps each argument in the custom flux function to the corresponding defined species: +$$ \phi_{ABC}= -10 c_A c_B c_C$$ + +We must first define each species using `Species` and then create the dictionary to be passed into `species_dependent_value`. The dictionary maps each argument in the custom flux function to the corresponding defined species. We also define our custom functions for the fluxes: ```{code-cell} ipython3 import festim as F @@ -80,59 +82,88 @@ my_model = F.HydrogenTransportProblem() A = F.Species(name="A") B = F.Species(name="B") C = F.Species(name="C") -my_model.species = [A, B, C] -my_custom_value = lambda c_A, c_B, c_C: 2*c_A + 3*c_B + 4*c_C +my_model.species = [A, B, C] species_dependent_value = {"c_A": A, "c_B": B, "c_C": C} + +recombination_flux_AB = lambda c_A, c_B, c_C: -c_A*c_B +recombination_flux_ABC = lambda c_A, c_B, c_C: -10*c_A*c_B*c_C ``` -Now, we create our 1D mesh and assign boundary conditions (flux BC on the left). The boundary condition `ParticleFluxBC` must be added for each species: +Now, we create our 1D mesh and assign boundary conditions (recombination on the right, fixed concentration on the left). + +The boundary condition `ParticleFluxBC` must be added for each species: ```{code-cell} ipython3 import numpy as np my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100)) -D = 1 +D = 1e-2 mat = F.Material( - D_0={A: D, B: D, C: D}, + D_0={A: 8*D, B: 7*D, C: D}, E_D={A: 0.01, B: 0.01, C: 0.01}, ) bulk = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat) left = F.SurfaceSubdomain1D(id=1, x=0) right = F.SurfaceSubdomain1D(id=2, x=1) - my_model.subdomains = [bulk, left, right] + my_model.boundary_conditions = [ F.ParticleFluxBC( - subdomain=left, - value=my_custom_value, + subdomain=right, + value=recombination_flux_AB, species=A, species_dependent_value=species_dependent_value, ), F.ParticleFluxBC( - subdomain=left, - value=my_custom_value, + subdomain=right, + value=recombination_flux_AB, species=B, species_dependent_value=species_dependent_value, ), F.ParticleFluxBC( - subdomain=left, - value=my_custom_value, + subdomain=right, + value=recombination_flux_ABC, + species=A, + species_dependent_value=species_dependent_value, + ), + F.ParticleFluxBC( + subdomain=right, + value=recombination_flux_ABC, + species=B, + species_dependent_value=species_dependent_value, + ), + F.ParticleFluxBC( + subdomain=right, + value=recombination_flux_ABC, species=C, species_dependent_value=species_dependent_value, ), + F.FixedConcentrationBC(subdomain=left,value=1,species=A), + F.FixedConcentrationBC(subdomain=left,value=1,species=B), + F.FixedConcentrationBC(subdomain=left,value=1,species=C), +] + +right_flux_A = F.SurfaceFlux(field=A,surface=right) +right_flux_B = F.SurfaceFlux(field=B,surface=right) +right_flux_C = F.SurfaceFlux(field=C,surface=right) + +my_model.exports = [ + right_flux_A, + right_flux_B, + right_flux_C, ] ``` ```{note} -The diffusivity pre-factor `D_0` and activation energy `E_d` must be defined for each species in `Material`. Learn more about defining multi-species material properties __[here](https://festim-workshop.readthedocs.io/en/festim2/content/material/material_advanced.html#defining-species-dependent-material-properties)__. +The diffusivity pre-factor `D_0` and activation energy `E_D` must be defined for each species in `Material`. Learn more about defining multi-species material properties __[here](https://festim-workshop.readthedocs.io/en/festim2/content/material/material_advanced.html#defining-species-dependent-material-properties)__. ``` +++ -Finally, let's solve and plot the solution for each species: +Finally, let's solve and plot the profile for each species: ```{code-cell} ipython3 my_model.temperature = 300 @@ -160,3 +191,5 @@ plt.ylabel('Concentration') plt.legend() plt.show() ``` + +We see the higher recombination flux for $\mathrm{ABC}$ decreases the concentration of $\mathrm{C}$ throughout the material. diff --git a/book/content/boundary_conditions/heat_transfer.md b/book/content/boundary_conditions/heat_transfer.md index e2b89126..e78df302 100644 --- a/book/content/boundary_conditions/heat_transfer.md +++ b/book/content/boundary_conditions/heat_transfer.md @@ -5,7 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.17.3 + jupytext_version: 1.18.1 --- # Heat transfer # diff --git a/book/content/boundary_conditions/surface_reactions.md b/book/content/boundary_conditions/surface_reactions.md index 15762910..0a42c384 100644 --- a/book/content/boundary_conditions/surface_reactions.md +++ b/book/content/boundary_conditions/surface_reactions.md @@ -5,7 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.17.3 + jupytext_version: 1.18.1 kernelspec: display_name: festim-workshop language: python @@ -19,23 +19,48 @@ kernelspec: FESTIM V2.0 allows users to impose surface reactions on boundaries. Learn more about surface reactions __[here](https://festim-workshop.readthedocs.io/en/festim2/content/species_reactions/surface_reactions.html)__. Objectives: -* Impose a simple surface reaction boundary condition -* recombination and dissociation -* Multi-isotope +* Defining a simple surface reaction boundary condition +* Recombination and dissociation +* Isotopic exhange +* Complex isotopic exchange with multple hydrogenic species +++ -## Simple surface reaction BC ## +## Defining a simple surface reaction boundary condition ## -A reaction is defined by specifying the reactants (`reactant`), gas pressure (`gas_pressure`), forward/backward rate constants (`k_r0` and `k_d0`), and rate activation energies (`E_kr` and `E_kd`). +Surface reactions between adsorbed species and gas-phase products can be represented using the `SurfaceReactionBC` class. This boundary condition defines a reversible reaction of the form: -```{note} -Reaction rates are given as Arhennius laws. -``` +$$ +A + B \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad C +$$ + +where $\mathrm{A}$ and $\mathrm{B}$ are surface reactants and $\mathrm{C}$ is the product species in the gas phase. + +### Reaction Rate Formulation ### + +The forward and backward rate constants follow Arrhenius laws: + +$$ +K_r = k_{r0} e^{-E_{kr} / (k_B T)}, \qquad +K_d = k_{d0} e^{-E_{kd} / (k_B T)} +$$ + +The net surface reaction rate is given by: + +$$ +K = K_r c_A c_B - K_d P_C +$$ + +where: +- $\mathrm{k_{r0}}, \mathrm{k_{d0}}$ are rate pre-exponential constants +- $\mathrm{c_A}, \mathrm{c_B}$ are concentrations of reactant species $\mathrm{A}$ and $\mathrm{B}$ at the surface +- $\mathrm{P_C}$ is the partial pressure of the product species $\mathrm{C}$ +- $\mathrm{k_B}$ is the Boltzmann constant +- $\mathrm{T}$ is the surface temperature -To impose the folllowing surface reaction: +The flux of species $\mathrm{A}$ entering the surface is equal to $\mathrm{K}$ (if $\mathrm{A=B}$, the total particle flux entering the surface becomes $\mathrm{2K}$). -$$ A + B \xrightarrow{k} C $$ +We can use the `SurfaceReactionBC` class to impose the surface reaction above by specifying the reactants (`reactant`), gas pressure (`gas_pressure`), forward/backward rate constants (`k_r0` and `k_d0`), and rate activation energies (`E_kr` and `E_kd`): ```{code-cell} ipython3 from festim import Species, SurfaceReactionBC, SurfaceSubdomain @@ -50,17 +75,17 @@ my_bc = SurfaceReactionBC( gas_pressure=1e5, k_r0=1, E_kr=0.1, - k_d0=1e-5, - E_kd=0.1, + k_d0=0, + E_kd=0, subdomain=boundary, ) ``` -## Recombination and Dissociation ## +## Recombination and dissociation ## Hydrogen recombination/dissociation can also be modelled using `SurfaceReactionBC`, such as this reaction: -$$ H + H \overset{K_r}{\underset{K_d}{\rightleftharpoons}} H_2 $$ +$$ \mathrm{H + H} \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad \mathrm{H_2} $$ ```{code-cell} ipython3 from festim import Species, SurfaceReactionBC, SurfaceSubdomain @@ -79,202 +104,229 @@ my_bc = SurfaceReactionBC( ) ``` -## Multiple isotopes ## +## Isotopic exchange ## + +Isotopic exchange can occur where isotopes can swap positions with other isotopes, such as: + +$$\mathrm{T + H_2} \rightleftharpoons \mathrm{H} + \mathrm{HT}$$ + +These exchange processes occur through complex mechanisms and are important in understanding the behavior of hydrogenic species: + +If the $\mathrm{H_2}$ concentration is assumed much larger than $\mathrm{HT}$, the flux $\phi_T$ reduces to a first-order process in $\mathrm{T}$: -Surface reactions can involve multiple hydrogen isotopes, enabling the modelling of more complex interactions between species. For example, in a system with both mobile hydrogen and tritium, various molecular recombination pathways may occur at the surface, resulting in the formation of $H_2$, $T_2$, and $HT$: +$$ +\phi_T = +-\mathrm{K_{r0}} \exp\!\left(-\frac{\mathrm{E_{Kr}}}{\mathrm{k_B T}}\right) c_T^2 +-\mathrm{K_{r0}^\ast} \exp\!\left(-\frac{\mathrm{E_{Kr}^\ast}}{\mathrm{k_B T}}\right) c_{H_2} c_T +$$ + +Such fluxes can be implemented using `festim.ParticleFluxBC` with user-defined expressions. ```{code-cell} ipython3 -from festim import Species, SurfaceReactionBC, SurfaceSubdomain +import ufl +import festim as F -boundary = SurfaceSubdomain(id=1) -H = Species("H") -T = Species("T") +Kr_0 = 1.0 +E_Kr = 0.1 -reac1 = SurfaceReactionBC( - reactant=[H, H], - gas_pressure=1e6, - k_r0=1.0, - E_kr=0.1, - k_d0=0.5, - E_kd=0.1, - subdomain=boundary, -) -reac2 = SurfaceReactionBC( - reactant=[T, T], - gas_pressure=1e6, - k_r0=1.0, - E_kr=0.1, - k_d0=0.5, - E_kd=0.1, - subdomain=boundary, -) -reac3 = SurfaceReactionBC( - reactant=[H, T], - gas_pressure=1e6, - k_r0=1.0, - E_kr=0.1, - k_d0=0.5, - E_kd=0.1, - subdomain=boundary, -) +def my_custom_recombination_flux(c, T): + Kr_0_custom = 1.0 + E_Kr_custom = 0.5 # eV + h2_conc = 1e25 # assumed constant H2 concentration in + + recombination_flux = ( + -(Kr_0 * ufl.exp(-E_Kr / (F.k_B * T))) * c**2 + - (Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))) * h2_conc * c + ) + return recombination_flux ``` -## Isotopic exchange ## +```{note} +We can also use `F.SurfaceReactionBC` to define recombination fluxes, which we do in the next section for more complex isotopic exchanges. +``` ```{code-cell} ipython3 -import dolfinx.fem as fem import numpy as np -import ufl +my_model = F.HydrogenTransportProblem() +my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 100)) + +left_surf = F.SurfaceSubdomain1D(id=1, x=0) +right_surf = F.SurfaceSubdomain1D(id=2, x=1) + +material = F.Material(D_0=1e-2, E_D=0) + +vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=material) + +my_model.subdomains = [vol, left_surf, right_surf] +tritium = F.Species("T") +my_model.species = [tritium] + +my_custom_flux = F.ParticleFluxBC( + value=my_custom_recombination_flux, + subdomain=right_surf, + species_dependent_value={"c": tritium}, + species=tritium, +) + +my_model.boundary_conditions = [ + my_custom_flux, + F.FixedConcentrationBC(subdomain=left_surf, value=1, species=tritium), +] + +my_model.temperature = 300 +right_flux = F.SurfaceFlux(field=tritium, surface=right_surf) + +my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) + +my_model.initialise() +my_model.run() + +data_with_H2 = (tritium.solution.x.array) +``` + +We can compare the surface flux to the case where we have no isotopic exchange by removing the custom boundary condition: + +```{code-cell} ipython3 +my_model.boundary_conditions = [ + F.FixedConcentrationBC(subdomain=left_surf, value=1, species=tritium), +] + +my_model.initialise() +my_model.run() +data_no_H2 = (tritium.solution.x.array) +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +import matplotlib.pyplot as plt + +x = my_model.mesh.mesh.geometry.x[:, 0] +plt.plot(x, data_no_H2, label="No H_2") +plt.plot(x, data_with_H2, label="With H_2") + +plt.xlabel("x (m)") +plt.ylabel("Surface Flux (right)") +plt.legend() +plt.show() +``` + +We see that diffusion is present with isotopic exchange! + ++++ + +## Complex isotopic exchange with multple hydrogenic species ## + +Surface reactions can involve multiple hydrogen isotopes, allowing for the modeling of complex isotope-exchange mechanisms between species. For example, in a system with both mobile hydrogen and deuteriun, various molecular recombination pathways may occur at the surface, resulting in the formation of $H_2$, $D_2$, and $HD$: + +$$ \text{Reaction 1}: \mathrm{H+D} \rightarrow \mathrm{HD} \longrightarrow \phi_1 = K_{r1} c_H c_D - K_{d1}P_{HD} $$ +$$ \text{Reaction 2}: \mathrm{D+D} \rightarrow \mathrm{D_2} \longrightarrow \phi_2 = 2K_{r2} c_D^2 - K_{d2}P_{D2} $$ +$$ \text{Reaction 3}: \mathrm{D+H_2} \rightarrow \mathrm{HD + H} \longrightarrow \phi_3 = K_{r3} c_H c_D - K_{d3}P_{HD} $$ +$$ \text{Reaction 4}: \mathrm{H+H} \rightarrow \mathrm{H_2} \longrightarrow \phi_4 = 2K_{r4} c_H^2 - K_{d4}P_{H2} $$ + +Now consider the case where deuterium diffuses from left to right and reacts with background +$\mathrm{H_2}$, while $\mathrm{P_{HD}}$ and $\mathrm{P_{D_2}}$ are negligible at the surface. +Formation of $\mathrm{H}$ at the right boundary induces back-diffusion toward the left, +even though none existed initially. + +The boundary conditions for this scenario are: +$$ +c_D(x=0) = 1, \qquad c_H(x=0) = 0, \qquad P_{H2}(x=1) = \text{1000 Pa} +$$ + +First, let's define a 1D mesh ranging from $\mathrm{x=[0,1]}$: + +```{code-cell} ipython3 +import numpy as np import festim as F my_model = F.HydrogenTransportProblem() -my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 1000)) -my_mat = F.Material(name="mat", D_0=1, E_D=0) -vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat) -left = F.SurfaceSubdomain1D(id=1, x=0) -right = F.SurfaceSubdomain1D(id=2, x=1) +my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 100)) -my_model.subdomains = [vol, left, right] +left_surf = F.SurfaceSubdomain1D(id=1, x=0) +right_surf = F.SurfaceSubdomain1D(id=2, x=1) +material = F.Material(D_0=1e-2, E_D=0) +vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=material) + +my_model.subdomains = [vol, left_surf, right_surf] +``` + +Now, we define our species at recombination reactions using `SurfaceReactionBC`: + +```{code-cell} ipython3 H = F.Species("H") D = F.Species("D") my_model.species = [H, D] -my_model.temperature = 500 +H2 = F.SurfaceReactionBC( + reactant=[H, H], + gas_pressure=100000, + k_r0=1, + E_kr=0.1, + k_d0=1e-5, + E_kd=0.1, + subdomain=right_surf, +) -surface_reaction_hd = F.SurfaceReactionBC( +HD = F.SurfaceReactionBC( reactant=[H, D], gas_pressure=0, - k_r0=0.01, - E_kr=0, - k_d0=0, - E_kd=0, - subdomain=right, -) - -surface_reaction_hh = F.SurfaceReactionBC( - reactant=[H, H], - gas_pressure=lambda t: ufl.conditional(ufl.gt(t, 1), 2, 0), - k_r0=0.02, - E_kr=0, - k_d0=0.03, - E_kd=0, - subdomain=right, + k_r0=1, + E_kr=0.1, + k_d0=1e-5, + E_kd=0.1, + subdomain=right_surf, ) -surface_reaction_dd = F.SurfaceReactionBC( +D2 = F.SurfaceReactionBC( reactant=[D, D], gas_pressure=0, - k_r0=0.01, - E_kr=0, - k_d0=0, - E_kd=0, - subdomain=right, + k_r0=1, + E_kr=0.1, + k_d0=1e-5, + E_kd=0.1, + subdomain=right_surf, ) +``` -my_model.boundary_conditions = [ - F.DirichletBC(subdomain=left, value=2, species=H), - F.DirichletBC(subdomain=left, value=2, species=D), - surface_reaction_hd, - surface_reaction_hh, - surface_reaction_dd, -] +Finally, we add our boundary conditions and solve the steady-state problem: -H_flux_right = F.SurfaceFlux(H, right) -H_flux_left = F.SurfaceFlux(H, left) -D_flux_right = F.SurfaceFlux(D, right) -D_flux_left = F.SurfaceFlux(D, left) - -my_model.exports = [ - F.XDMFExport("test.xdmf", H), - H_flux_left, - H_flux_right, - D_flux_left, - D_flux_right, - HD_flux, - HH_flux, - DD_flux, +```{code-cell} ipython3 +my_model.boundary_conditions = [ + H2, + D2, + HD, + F.FixedConcentrationBC(subdomain=left_surf, value=1, species=D), + F.FixedConcentrationBC(subdomain=left_surf, value=0, species=H), ] -my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=5, transient=True) +my_model.temperature = 300 -my_model.settings.stepsize = 0.1 +my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) my_model.initialise() my_model.run() +``` + +```{code-cell} ipython3 +:tags: [hide-input] import matplotlib.pyplot as plt -plt.stackplot( - H_flux_left.t, - np.abs(H_flux_left.data), - np.abs(D_flux_left.data), - labels=["H_in", "D_in"], -) -plt.stackplot( - H_flux_right.t, - -np.abs(H_flux_right.data), - -np.abs(D_flux_right.data), - labels=["H_out", "D_out"], -) -plt.legend() -plt.xlabel("Time (s)") -plt.ylabel("Flux (atom/m^2/s)") -plt.figure() -plt.stackplot( - HD_flux.t, - np.abs(HH_flux.data), - np.abs(HD_flux.data), - np.abs(DD_flux.data), - labels=["HH", "HD", "DD"], -) -plt.legend(reverse=True) -plt.xlabel("Time (s)") -plt.ylabel("Flux (molecule/m^2/s)") - - -plt.figure() -plt.plot(H_flux_right.t, -np.array(H_flux_right.data), label="from gradient (H)") -plt.plot( - H_flux_right.t, - 2 * np.array(HH_flux.data) + np.array(HD_flux.data), - linestyle="--", - label="from reaction rates (H)", -) +def plot_profile(species, **kwargs): + c = species.post_processing_solution.x.array[:] + x = species.post_processing_solution.function_space.mesh.geometry.x[:,0] + return plt.plot(x, c, **kwargs) -plt.plot(D_flux_right.t, -np.array(D_flux_right.data), label="from gradient (D)") -plt.plot( - D_flux_right.t, - 2 * np.array(DD_flux.data) + np.array(HD_flux.data), - linestyle="--", - label="from reaction rates (D)", -) -plt.xlabel("Time (s)") -plt.ylabel("Flux (atom/m^2/s)") +for species in my_model.species: + plot_profile(species, label=species.name) + +plt.xlabel('Position') +plt.ylabel('Concentration') plt.legend() plt.show() - -# check that H_flux_right == 2*HH_flux + HD_flux -H_flux_from_gradient = -np.array(H_flux_right.data) -H_flux_from_reac = 2 * np.array(HH_flux.data) + np.array(HD_flux.data) -assert np.allclose( - H_flux_from_gradient, - H_flux_from_reac, - rtol=0.5e-2, - atol=0.005, -) -# check that D_flux_right == 2*DD_flux + HD_flux -D_flux_from_gradient = -np.array(D_flux_right.data) -D_flux_from_reac = 2 * np.array(DD_flux.data) + np.array(HD_flux.data) -assert np.allclose( - D_flux_from_gradient, - D_flux_from_reac, - rtol=0.5e-2, - atol=0.005, -) ``` -```{code-cell} ipython3 - -``` +We see that the background $\mathrm{H_2}$ reacts with the $\mathrm{D}$, removing the total amount of $\mathrm{D}$ from the surface. Conversely, the $\mathrm{H}$ diffuses from the surface towards the left. From 44b3887945284a4ad5da18406b48c16864d1cc56 Mon Sep 17 00:00:00 2001 From: Chirag Date: Tue, 28 Oct 2025 14:51:10 -0400 Subject: [PATCH 04/17] fixed some grammar errors --- book/content/boundary_conditions/surface_reactions.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/book/content/boundary_conditions/surface_reactions.md b/book/content/boundary_conditions/surface_reactions.md index 0a42c384..87645651 100644 --- a/book/content/boundary_conditions/surface_reactions.md +++ b/book/content/boundary_conditions/surface_reactions.md @@ -145,6 +145,10 @@ def my_custom_recombination_flux(c, T): We can also use `F.SurfaceReactionBC` to define recombination fluxes, which we do in the next section for more complex isotopic exchanges. ``` ++++ + +Let's work through a 1D example to illustrate the effect of isotopic exchange. We also define a fixed concentration on the left surface: + ```{code-cell} ipython3 import numpy as np @@ -231,6 +235,7 @@ Formation of $\mathrm{H}$ at the right boundary induces back-diffusion toward th even though none existed initially. The boundary conditions for this scenario are: + $$ c_D(x=0) = 1, \qquad c_H(x=0) = 0, \qquad P_{H2}(x=1) = \text{1000 Pa} $$ From a6a2e8951472d0d44a1ef25dc54334f13ec2f4c9 Mon Sep 17 00:00:00 2001 From: Chirag Date: Tue, 10 Feb 2026 08:57:29 -0500 Subject: [PATCH 05/17] re-organizing bcs chapter --- .../boundary_conditions/concentration.md | 157 +++++++- .../h_transport_advanced.ipynb | 0 .../h_transport_basic.ipynb | 341 ++++++++++++++++++ 3 files changed, 488 insertions(+), 10 deletions(-) create mode 100644 book/content/boundary_conditions/h_transport_advanced.ipynb create mode 100644 book/content/boundary_conditions/h_transport_basic.ipynb diff --git a/book/content/boundary_conditions/concentration.md b/book/content/boundary_conditions/concentration.md index ec09ff61..f02b76ae 100644 --- a/book/content/boundary_conditions/concentration.md +++ b/book/content/boundary_conditions/concentration.md @@ -19,37 +19,174 @@ Boundary conditions (BCs) are essential to FESTIM simulations, as they describe This tutorial goes over how to define concentration boundary conditions for hydrogen transport simulations. Objectives: +* Understand the mathematics behind a fixed concentration BC * Define a fixed concentration BC * Choose which solubility law (Sieverts' or Henry's) * Solve a hydrogen transport problem with plasma implantation +++ +## Understanding mathematics behind a fixed concentration BC ## + +A fixed concentration (Dirichlet) boundary condition prescribes the value of the mobile hydrogen isotope concentration at a boundary. This enforces the concentration to remain constant in time and space on the specified boundary, independent of the solution in the bulk. + +This boundary condition is typically used to represent surfaces in equilibrium with an infinite reservoir, imposed implantation conditions, or experimentally controlled concentrations. + +### Mathematical formulation + +On a boundary $\Gamma_D$, the mobile concentration satisfies + +$$ +c(\mathbf{x}, t) = c_0 \quad \text{for } \mathbf{x} \in \Gamma_D, +$$ + +where $c_0$ is the prescribed concentration value. + +In the weak formulation, this condition is enforced by directly constraining the degrees of freedom associated with the boundary. + ++++ + ## Defining fixed concentration ## -BCs need to be assigned to surfaces using FESTIM's `SurfaceSubdomain` class. To define the concentration of a specific species, we can use `FixedConcentrationBC`: +Users can prescribe a fixed concentration on a boundary using `FixedConcentrationBC`, which can depend on temperature, time, and space. + +Let us consider a 2D, steady-state, single-species example with the following *space-dependent* boundary conditions: + +$$ \text{Top:} \quad c = 2x + y $$ +$$ \text{Right:} \quad c = y^2 + 1 $$ ```{code-cell} ipython3 -from festim import FixedConcentrationBC, Species, SurfaceSubdomain +import festim as F +import numpy as np +from dolfinx.mesh import create_unit_square +from mpi4py import MPI -boundary = SurfaceSubdomain(id=1) -H = Species(name="Hydrogen") +top_subdomain = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1)) +right_subdomain = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1)) +left_subdomain = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0)) +bottom_subdomain = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[1], 0)) + +my_model = F.HydrogenTransportProblem() -my_bc = FixedConcentrationBC(subdomain=boundary, value=10, species=H) +my_model.mesh = F.Mesh(create_unit_square(MPI.COMM_WORLD, 10, 10)) + +H = F.Species("H") +my_model.species = [H] +mat = F.Material(D_0=1e-3, E_D=0) +my_model.material = mat +vol = F.VolumeSubdomain(id=5, material=mat) +my_model.subdomains = [top_subdomain, right_subdomain, left_subdomain, bottom_subdomain, vol] + +my_model.temperature = 400 +my_model.settings = F.Settings(atol=1e-8, rtol=1e-8, transient=False) ``` -The imposed concentration can be dependent on space, time and temperature, such as: +We can define functions for the top and right boundary conditions using `lambda` functions: -$$ +```{code-cell} ipython3 +top_bc_function = lambda x: 2.0 * x[0] + x[1] +right_bc_function = lambda x: x[1] ** 2 + 1.0 +``` -BC = 10 + x^2 + T(t+2) + -$$ ++++ + +To include these boundary conditions to our problem, we use `FixedConcentrationBC`. We must also specify which subdomain (boundary) each BC is applied to, as well as the corresponding species: + +```{code-cell} ipython3 +top_bc = F.FixedConcentrationBC( + subdomain=top_subdomain, + value=top_bc_function, + species=H +) + +right_bc = F.FixedConcentrationBC( + subdomain=right_subdomain, + value=right_bc_function, + species=H +) + +# left_bc = F.FixedConcentrationBC( +# subdomain=left_subdomain, +# value=.0, +# species=H +# ) + +# bottom_bc = F.FixedConcentrationBC( +# subdomain=bottom_subdomain, +# value=0.0, +# species=H +# ) +``` + +Finally, we add our BCs to `my_model.boundary_conditions` using a list, and then run the model: + +```{code-cell} ipython3 +my_model.boundary_conditions = [top_bc, right_bc] +my_model.initialise() +my_model.run() +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +import pyvista +from dolfinx import plot + +topology, cell_types, geometry = plot.vtk_mesh(H.post_processing_solution.function_space) +grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) +grid.point_data["c"] = H.post_processing_solution.x.array +grid.set_active_scalars("c") + +pyvista.start_xvfb() +pyvista.set_jupyter_backend("html") + +plotter = pyvista.Plotter() + +plotter.add_mesh(grid) +plotter.view_xy() +contours = grid.contour(isosurfaces=5, scalars="c") +plotter.add_mesh(contours, color="r", line_width=0.1) +if not pyvista.OFF_SCREEN: + plotter.show() +else: + figure = plotter.screenshot("concentration.png") +``` + + + ++++ + +```{note} + ++++ + +### Time and temperature dependent boundary conditions ### + ++++ + +Users can also impose concentration BCs that are dependent on space, time and temperature, such as: + +$$ c = 10 + x^2 + T(t+2) $$ + +To do so, we define a Python function: ```{code-cell} ipython3 my_custom_value = lambda x, t, T: 10 + x[0]**2 + T *(t + 2) -my_bc = FixedConcentrationBC(subdomain=boundary, value=my_custom_value, species=H) +boundary = F.SurfaceSubdomain(id=1) +my_bc = F.FixedConcentrationBC(subdomain=boundary, value=my_custom_value, species=H) ``` ## Choosing a solubility law ## diff --git a/book/content/boundary_conditions/h_transport_advanced.ipynb b/book/content/boundary_conditions/h_transport_advanced.ipynb new file mode 100644 index 00000000..e69de29b diff --git a/book/content/boundary_conditions/h_transport_basic.ipynb b/book/content/boundary_conditions/h_transport_basic.ipynb new file mode 100644 index 00000000..2f29a243 --- /dev/null +++ b/book/content/boundary_conditions/h_transport_basic.ipynb @@ -0,0 +1,341 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4740a699", + "metadata": {}, + "source": [ + "# Hydrogen transport: basic boundary conditions\n", + "\n", + "This section discusses how to implement basic boundary conditions (BCs) for hydrogen transport problems. Boundary conditions are essential to FESTIM simulations, as they describe the mathematical problem at the boundaries of the simulated domain. Read more about BCs _[here](https://festim.readthedocs.io/en/fenicsx/userguide/boundary_conditions.html)_.\n", + "\n", + "Objectives:\n", + "* Learn mathematical formulation for concentration and flux boundary conditions\n", + "* Implement fixed concentration BCs\n", + "* Implement particle flux BCs" + ] + }, + { + "cell_type": "markdown", + "id": "8dbd044f", + "metadata": {}, + "source": [ + "## Understanding math behind concentration and flux boundary conditions\n", + "\n", + "### Fixed concentration\n", + "\n", + "A fixed concentration (Dirichlet) boundary condition prescribes the value of the mobile hydrogen isotope concentration at a boundary. This enforces the concentration to remain constant in time and space on the specified boundary, independent of the solution in the bulk.\n", + "\n", + "This boundary condition is typically used to represent surfaces in equilibrium with an infinite reservoir, imposed implantation conditions, or experimentally controlled concentrations.\n", + "\n", + "#### Mathematical formulation\n", + "\n", + "On a boundary $\\Gamma_D$, the mobile concentration satisfies\n", + "\n", + "$$\n", + "c(\\mathbf{x}, t) = c_0 \\quad \\text{for } \\mathbf{x} \\in \\Gamma_D,\n", + "$$\n", + "\n", + "where $c_0$ is the prescribed concentration value.\n", + "\n", + "In the weak formulation, this condition is enforced by directly constraining the degrees of freedom associated with the boundary.\n", + "\n", + "### Particle flux boundary conditions\n", + "\n", + "A particle flux (Neumann) boundary condition prescribes the normal flux of mobile hydrogen isotopes across a boundary. Unlike fixed concentration conditions, flux boundary conditions do not directly constrain the concentration value at the surface. Instead, they control the rate at which particles enter or leave the domain.\n", + "\n", + "Flux boundary conditions are commonly used to represent implantation from a plasma, outgassing to vacuum, permeation through a surface, or symmetry boundaries where no net transport occurs.\n", + "\n", + "#### Mathematical formulation\n", + "\n", + "The particle flux $\\mathbf{J}$ is typically given by Fick’s law,\n", + "\n", + "$$\n", + "\\mathbf{J} = -D \\nabla c,\n", + "$$\n", + "\n", + "where $D$ is the diffusion coefficient and $c$ is the mobile concentration.\n", + "\n", + "On a boundary $\\Gamma_N$, a prescribed normal flux is imposed as\n", + "\n", + "$$\n", + "\\mathbf{J} \\cdot \\mathbf{n} = g(\\mathbf{x}, t) \\quad \\text{for } \\mathbf{x} \\in \\Gamma_N,\n", + "$$\n", + "\n", + "where:\n", + "- $\\mathbf{n}$ is the outward unit normal vector,\n", + "- $g(\\mathbf{x}, t)$ is the imposed particle flux (positive for outward flux, negative for inward flux).\n", + "\n", + "A special case is the **zero-flux (symmetry) boundary condition**, given by\n", + "\n", + "$$\n", + "\\mathbf{J} \\cdot \\mathbf{n} = 0 \\quad \\text{on } \\Gamma,\n", + "$$\n", + "\n", + "which implies no net particle transport across the boundary.\n", + "\n", + "#### Weak formulation\n", + "\n", + "In the weak form, flux boundary conditions appear naturally as surface integrals after integration by parts. For a test function $v$, the boundary contribution is\n", + "\n", + "$$\n", + "\\int_{\\Gamma_N} g(\\mathbf{x}, t)\\, v \\, \\mathrm{d}\\Gamma.\n", + "$$\n", + "\n", + "Because of this, flux boundary conditions are often referred to as *natural boundary conditions* and do not require explicit modification of the solution space, unlike Dirichlet conditions.\n" + ] + }, + { + "cell_type": "markdown", + "id": "885b1855", + "metadata": {}, + "source": [ + "## Imposing a fixed concentration BC \n", + "\n", + "Users can prescribe a fixed concentration on a boundary using `FixedConcentrationBC`, which can depend on temperature, time, and space.\n", + "\n", + "Let us consider a 2D, steady-state, single-species example with the following *space-dependent* boundary conditions:\n", + "\n", + "$$ \\text{Top:} \\quad c = 2x + y $$\n", + "$$ \\text{Right:} \\quad c = y^2 + 1 $$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea3994c1", + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "import festim as F\n", + "import numpy as np\n", + "from dolfinx.mesh import create_unit_square\n", + "from mpi4py import MPI\n", + "\n", + "top_subdomain = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1))\n", + "right_subdomain = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1))\n", + "left_subdomain = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0))\n", + "bottom_subdomain = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[1], 0))\n", + "\n", + "my_model = F.HydrogenTransportProblem()\n", + "\n", + "my_model.mesh = F.Mesh(create_unit_square(MPI.COMM_WORLD, 10, 10))\n", + "\n", + "H = F.Species(\"H\")\n", + "my_model.species = [H]\n", + "mat = F.Material(D_0=1e-3, E_D=0)\n", + "my_model.material = mat\n", + "vol = F.VolumeSubdomain(id=5, material=mat)\n", + "my_model.subdomains = [top_subdomain, right_subdomain, left_subdomain, bottom_subdomain, vol]\n", + "\n", + "my_model.temperature = 400\n", + "my_model.settings = F.Settings(atol=1e-8, rtol=1e-8, transient=False)" + ] + }, + { + "cell_type": "markdown", + "id": "2bc407e0", + "metadata": {}, + "source": [ + "We can define functions for the top and right boundary conditions using `lambda` functions:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "683f609b", + "metadata": {}, + "outputs": [], + "source": [ + "top_bc_function = lambda x: 2.0 * x[0] + x[1]\n", + "right_bc_function = lambda x: x[1] ** 2 + 1.0" + ] + }, + { + "cell_type": "markdown", + "id": "2552b567", + "metadata": {}, + "source": [ + "```{note}\n", + "`x[0]`, `x[1]`, `x[2]` corresponse to *(x, y, z)* in cartesian space.\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "fee11d25", + "metadata": {}, + "source": [ + "To include these boundary conditions to our problem, we use `FixedConcentrationBC`. We must also specify which subdomain (boundary) each BC is applied to, as well as the corresponding species:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "372c20cb", + "metadata": {}, + "outputs": [], + "source": [ + "top_bc = F.FixedConcentrationBC(\n", + " subdomain=top_subdomain,\n", + " value=top_bc_function,\n", + " species=H\n", + ")\n", + "\n", + "right_bc = F.FixedConcentrationBC(\n", + " subdomain=right_subdomain,\n", + " value=right_bc_function,\n", + " species=H\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "71012707", + "metadata": {}, + "outputs": [], + "source": [ + "my_model.boundary_conditions = [top_bc, right_bc]\n", + "my_model.initialise()\n", + "my_model.run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af055320", + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "import pyvista\n", + "from dolfinx import plot\n", + "\n", + "topology, cell_types, geometry = plot.vtk_mesh(H.post_processing_solution.function_space)\n", + "grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)\n", + "grid.point_data[\"c\"] = H.post_processing_solution.x.array\n", + "grid.set_active_scalars(\"c\")\n", + "\n", + "pyvista.start_xvfb()\n", + "pyvista.set_jupyter_backend(\"html\")\n", + "\n", + "plotter = pyvista.Plotter()\n", + "\n", + "plotter.add_mesh(grid)\n", + "plotter.view_xy()\n", + "contours = grid.contour(isosurfaces=5, scalars=\"c\")\n", + "plotter.add_mesh(contours, color=\"r\", line_width=0.1)\n", + "if not pyvista.OFF_SCREEN:\n", + " plotter.show()\n", + "else:\n", + " figure = plotter.screenshot(\"concentration.png\")" + ] + }, + { + "cell_type": "markdown", + "id": "5e3b21c6", + "metadata": {}, + "source": [ + "```{note}\n", + "\n", + "In this example, we did not define any boundary conditions for the left and bottom surfaces. If no BC is set on a boundary, it is assumed that the flux is null (symmetry boundary condition):\n", + "\n", + "$$\n", + "\\mathbf{J} \\cdot \\mathbf{n} = 0 \\quad \\text{on } \\Gamma\n", + "$$\n", + "\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "c2a7b30e", + "metadata": {}, + "source": [ + "### Time and temperature dependent boundary conditions ###" + ] + }, + { + "cell_type": "markdown", + "id": "eebf6d4e", + "metadata": {}, + "source": [ + "Users can also impose concentration BCs that are dependent on space, time and temperature, such as:\n", + "\n", + "$$ c = 10 + x^2 + T(t+2) $$\n", + "\n", + "To do so, we define a Python function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7e12721b", + "metadata": {}, + "outputs": [], + "source": [ + "my_custom_value = lambda x, t, T: 10 + x[0]**2 + T *(t + 2)\n", + "\n", + "boundary = F.SurfaceSubdomain(id=1)\n", + "my_bc = F.FixedConcentrationBC(subdomain=boundary, value=my_custom_value, species=H)" + ] + }, + { + "cell_type": "markdown", + "id": "ffb1284e", + "metadata": {}, + "source": [ + "```{tip}\n", + "This is equivalent to:\n", + "`\n", + "def my_custom_value(x, t, T):\n", + " return 10 + x[0]**2 + T *(t + 2)\n", + "`\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "2619c3a3", + "metadata": {}, + "source": [ + "### Multi-species fixed concentration BC ##" + ] + }, + { + "cell_type": "markdown", + "id": "520174bb", + "metadata": {}, + "source": [ + "## Imposing a particle flux BC" + ] + }, + { + "cell_type": "markdown", + "id": "71b75467", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "festim-workshop", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.10.18" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From dd4ba2a1b8dc5d62d862c725147b86fec91a2bf1 Mon Sep 17 00:00:00 2001 From: ck768 Date: Tue, 10 Feb 2026 13:32:15 -0500 Subject: [PATCH 06/17] moving fluxes, concentration bcs to basic hydrogen transport chapter --- book/_toc.yml | 4 +- .../h_transport_advanced.ipynb | 0 .../h_transport_basic.ipynb | 341 ------------- .../boundary_conditions/h_transport_basic.md | 454 ++++++++++++++++++ 4 files changed, 455 insertions(+), 344 deletions(-) delete mode 100644 book/content/boundary_conditions/h_transport_advanced.ipynb delete mode 100644 book/content/boundary_conditions/h_transport_basic.ipynb create mode 100644 book/content/boundary_conditions/h_transport_basic.md diff --git a/book/_toc.yml b/book/_toc.yml index ee441bbf..4e99a860 100644 --- a/book/_toc.yml +++ b/book/_toc.yml @@ -23,9 +23,7 @@ parts: - caption: Boundary Conditions chapters: - - file: content/boundary_conditions/concentration - - file: content/boundary_conditions/fluxes - - file: content/boundary_conditions/surface_reactions + - file: content/boundary_conditions/h_transport_basic - file: content/boundary_conditions/heat_transfer - caption: Species & Reactions diff --git a/book/content/boundary_conditions/h_transport_advanced.ipynb b/book/content/boundary_conditions/h_transport_advanced.ipynb deleted file mode 100644 index e69de29b..00000000 diff --git a/book/content/boundary_conditions/h_transport_basic.ipynb b/book/content/boundary_conditions/h_transport_basic.ipynb deleted file mode 100644 index 2f29a243..00000000 --- a/book/content/boundary_conditions/h_transport_basic.ipynb +++ /dev/null @@ -1,341 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "4740a699", - "metadata": {}, - "source": [ - "# Hydrogen transport: basic boundary conditions\n", - "\n", - "This section discusses how to implement basic boundary conditions (BCs) for hydrogen transport problems. Boundary conditions are essential to FESTIM simulations, as they describe the mathematical problem at the boundaries of the simulated domain. Read more about BCs _[here](https://festim.readthedocs.io/en/fenicsx/userguide/boundary_conditions.html)_.\n", - "\n", - "Objectives:\n", - "* Learn mathematical formulation for concentration and flux boundary conditions\n", - "* Implement fixed concentration BCs\n", - "* Implement particle flux BCs" - ] - }, - { - "cell_type": "markdown", - "id": "8dbd044f", - "metadata": {}, - "source": [ - "## Understanding math behind concentration and flux boundary conditions\n", - "\n", - "### Fixed concentration\n", - "\n", - "A fixed concentration (Dirichlet) boundary condition prescribes the value of the mobile hydrogen isotope concentration at a boundary. This enforces the concentration to remain constant in time and space on the specified boundary, independent of the solution in the bulk.\n", - "\n", - "This boundary condition is typically used to represent surfaces in equilibrium with an infinite reservoir, imposed implantation conditions, or experimentally controlled concentrations.\n", - "\n", - "#### Mathematical formulation\n", - "\n", - "On a boundary $\\Gamma_D$, the mobile concentration satisfies\n", - "\n", - "$$\n", - "c(\\mathbf{x}, t) = c_0 \\quad \\text{for } \\mathbf{x} \\in \\Gamma_D,\n", - "$$\n", - "\n", - "where $c_0$ is the prescribed concentration value.\n", - "\n", - "In the weak formulation, this condition is enforced by directly constraining the degrees of freedom associated with the boundary.\n", - "\n", - "### Particle flux boundary conditions\n", - "\n", - "A particle flux (Neumann) boundary condition prescribes the normal flux of mobile hydrogen isotopes across a boundary. Unlike fixed concentration conditions, flux boundary conditions do not directly constrain the concentration value at the surface. Instead, they control the rate at which particles enter or leave the domain.\n", - "\n", - "Flux boundary conditions are commonly used to represent implantation from a plasma, outgassing to vacuum, permeation through a surface, or symmetry boundaries where no net transport occurs.\n", - "\n", - "#### Mathematical formulation\n", - "\n", - "The particle flux $\\mathbf{J}$ is typically given by Fick’s law,\n", - "\n", - "$$\n", - "\\mathbf{J} = -D \\nabla c,\n", - "$$\n", - "\n", - "where $D$ is the diffusion coefficient and $c$ is the mobile concentration.\n", - "\n", - "On a boundary $\\Gamma_N$, a prescribed normal flux is imposed as\n", - "\n", - "$$\n", - "\\mathbf{J} \\cdot \\mathbf{n} = g(\\mathbf{x}, t) \\quad \\text{for } \\mathbf{x} \\in \\Gamma_N,\n", - "$$\n", - "\n", - "where:\n", - "- $\\mathbf{n}$ is the outward unit normal vector,\n", - "- $g(\\mathbf{x}, t)$ is the imposed particle flux (positive for outward flux, negative for inward flux).\n", - "\n", - "A special case is the **zero-flux (symmetry) boundary condition**, given by\n", - "\n", - "$$\n", - "\\mathbf{J} \\cdot \\mathbf{n} = 0 \\quad \\text{on } \\Gamma,\n", - "$$\n", - "\n", - "which implies no net particle transport across the boundary.\n", - "\n", - "#### Weak formulation\n", - "\n", - "In the weak form, flux boundary conditions appear naturally as surface integrals after integration by parts. For a test function $v$, the boundary contribution is\n", - "\n", - "$$\n", - "\\int_{\\Gamma_N} g(\\mathbf{x}, t)\\, v \\, \\mathrm{d}\\Gamma.\n", - "$$\n", - "\n", - "Because of this, flux boundary conditions are often referred to as *natural boundary conditions* and do not require explicit modification of the solution space, unlike Dirichlet conditions.\n" - ] - }, - { - "cell_type": "markdown", - "id": "885b1855", - "metadata": {}, - "source": [ - "## Imposing a fixed concentration BC \n", - "\n", - "Users can prescribe a fixed concentration on a boundary using `FixedConcentrationBC`, which can depend on temperature, time, and space.\n", - "\n", - "Let us consider a 2D, steady-state, single-species example with the following *space-dependent* boundary conditions:\n", - "\n", - "$$ \\text{Top:} \\quad c = 2x + y $$\n", - "$$ \\text{Right:} \\quad c = y^2 + 1 $$" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ea3994c1", - "metadata": { - "tags": [ - "hide-input" - ] - }, - "outputs": [], - "source": [ - "import festim as F\n", - "import numpy as np\n", - "from dolfinx.mesh import create_unit_square\n", - "from mpi4py import MPI\n", - "\n", - "top_subdomain = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1))\n", - "right_subdomain = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1))\n", - "left_subdomain = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0))\n", - "bottom_subdomain = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[1], 0))\n", - "\n", - "my_model = F.HydrogenTransportProblem()\n", - "\n", - "my_model.mesh = F.Mesh(create_unit_square(MPI.COMM_WORLD, 10, 10))\n", - "\n", - "H = F.Species(\"H\")\n", - "my_model.species = [H]\n", - "mat = F.Material(D_0=1e-3, E_D=0)\n", - "my_model.material = mat\n", - "vol = F.VolumeSubdomain(id=5, material=mat)\n", - "my_model.subdomains = [top_subdomain, right_subdomain, left_subdomain, bottom_subdomain, vol]\n", - "\n", - "my_model.temperature = 400\n", - "my_model.settings = F.Settings(atol=1e-8, rtol=1e-8, transient=False)" - ] - }, - { - "cell_type": "markdown", - "id": "2bc407e0", - "metadata": {}, - "source": [ - "We can define functions for the top and right boundary conditions using `lambda` functions:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "683f609b", - "metadata": {}, - "outputs": [], - "source": [ - "top_bc_function = lambda x: 2.0 * x[0] + x[1]\n", - "right_bc_function = lambda x: x[1] ** 2 + 1.0" - ] - }, - { - "cell_type": "markdown", - "id": "2552b567", - "metadata": {}, - "source": [ - "```{note}\n", - "`x[0]`, `x[1]`, `x[2]` corresponse to *(x, y, z)* in cartesian space.\n", - "```" - ] - }, - { - "cell_type": "markdown", - "id": "fee11d25", - "metadata": {}, - "source": [ - "To include these boundary conditions to our problem, we use `FixedConcentrationBC`. We must also specify which subdomain (boundary) each BC is applied to, as well as the corresponding species:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "372c20cb", - "metadata": {}, - "outputs": [], - "source": [ - "top_bc = F.FixedConcentrationBC(\n", - " subdomain=top_subdomain,\n", - " value=top_bc_function,\n", - " species=H\n", - ")\n", - "\n", - "right_bc = F.FixedConcentrationBC(\n", - " subdomain=right_subdomain,\n", - " value=right_bc_function,\n", - " species=H\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "71012707", - "metadata": {}, - "outputs": [], - "source": [ - "my_model.boundary_conditions = [top_bc, right_bc]\n", - "my_model.initialise()\n", - "my_model.run()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "af055320", - "metadata": { - "tags": [ - "hide-input" - ] - }, - "outputs": [], - "source": [ - "import pyvista\n", - "from dolfinx import plot\n", - "\n", - "topology, cell_types, geometry = plot.vtk_mesh(H.post_processing_solution.function_space)\n", - "grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)\n", - "grid.point_data[\"c\"] = H.post_processing_solution.x.array\n", - "grid.set_active_scalars(\"c\")\n", - "\n", - "pyvista.start_xvfb()\n", - "pyvista.set_jupyter_backend(\"html\")\n", - "\n", - "plotter = pyvista.Plotter()\n", - "\n", - "plotter.add_mesh(grid)\n", - "plotter.view_xy()\n", - "contours = grid.contour(isosurfaces=5, scalars=\"c\")\n", - "plotter.add_mesh(contours, color=\"r\", line_width=0.1)\n", - "if not pyvista.OFF_SCREEN:\n", - " plotter.show()\n", - "else:\n", - " figure = plotter.screenshot(\"concentration.png\")" - ] - }, - { - "cell_type": "markdown", - "id": "5e3b21c6", - "metadata": {}, - "source": [ - "```{note}\n", - "\n", - "In this example, we did not define any boundary conditions for the left and bottom surfaces. If no BC is set on a boundary, it is assumed that the flux is null (symmetry boundary condition):\n", - "\n", - "$$\n", - "\\mathbf{J} \\cdot \\mathbf{n} = 0 \\quad \\text{on } \\Gamma\n", - "$$\n", - "\n", - "```" - ] - }, - { - "cell_type": "markdown", - "id": "c2a7b30e", - "metadata": {}, - "source": [ - "### Time and temperature dependent boundary conditions ###" - ] - }, - { - "cell_type": "markdown", - "id": "eebf6d4e", - "metadata": {}, - "source": [ - "Users can also impose concentration BCs that are dependent on space, time and temperature, such as:\n", - "\n", - "$$ c = 10 + x^2 + T(t+2) $$\n", - "\n", - "To do so, we define a Python function:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7e12721b", - "metadata": {}, - "outputs": [], - "source": [ - "my_custom_value = lambda x, t, T: 10 + x[0]**2 + T *(t + 2)\n", - "\n", - "boundary = F.SurfaceSubdomain(id=1)\n", - "my_bc = F.FixedConcentrationBC(subdomain=boundary, value=my_custom_value, species=H)" - ] - }, - { - "cell_type": "markdown", - "id": "ffb1284e", - "metadata": {}, - "source": [ - "```{tip}\n", - "This is equivalent to:\n", - "`\n", - "def my_custom_value(x, t, T):\n", - " return 10 + x[0]**2 + T *(t + 2)\n", - "`\n", - "```" - ] - }, - { - "cell_type": "markdown", - "id": "2619c3a3", - "metadata": {}, - "source": [ - "### Multi-species fixed concentration BC ##" - ] - }, - { - "cell_type": "markdown", - "id": "520174bb", - "metadata": {}, - "source": [ - "## Imposing a particle flux BC" - ] - }, - { - "cell_type": "markdown", - "id": "71b75467", - "metadata": {}, - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "festim-workshop", - "language": "python", - "name": "python3" - }, - "language_info": { - "name": "python", - "version": "3.10.18" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/book/content/boundary_conditions/h_transport_basic.md b/book/content/boundary_conditions/h_transport_basic.md new file mode 100644 index 00000000..e85dfc8e --- /dev/null +++ b/book/content/boundary_conditions/h_transport_basic.md @@ -0,0 +1,454 @@ +--- +jupytext: + formats: ipynb,md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.17.3 +kernelspec: + display_name: festim-workshop + language: python + name: python3 +--- + +# Hydrogen transport: basic boundary conditions + +This section discusses how to implement basic boundary conditions (BCs) for hydrogen transport problems. Boundary conditions are essential to FESTIM simulations, as they describe the mathematical problem at the boundaries of the simulated domain. Read more about BCs _[here](https://festim.readthedocs.io/en/fenicsx/userguide/boundary_conditions.html)_. + +Objectives: +* Learn mathematical formulation for concentration and flux boundary conditions +* Implement fixed concentration boundary conditions +* Implement particle flux boundary conditions + ++++ + +## Understanding math behind concentration and flux boundary conditions + +### Fixed concentration + +A fixed concentration (Dirichlet) boundary condition prescribes the value of the mobile hydrogen isotope concentration at a boundary. This enforces the concentration to remain constant in time and space on the specified boundary, independent of the solution in the bulk. + +This boundary condition is typically used to represent surfaces in equilibrium with an infinite reservoir, imposed implantation conditions, or experimentally controlled concentrations. + +#### Mathematical formulation + +On a boundary $\Gamma_D$, the mobile concentration satisfies + +$$ +c(\mathbf{x}, t) = c_0 \quad \text{for } \mathbf{x} \in \Gamma_D, +$$ + +where $c_0$ is the prescribed concentration value. + +In the weak formulation, this condition is enforced by directly constraining the degrees of freedom associated with the boundary. + +### Particle flux boundary conditions + +A particle flux (Neumann) boundary condition prescribes the normal flux of mobile hydrogen isotopes across a boundary. Unlike fixed concentration conditions, flux boundary conditions do not directly constrain the concentration value at the surface. Instead, they control the rate at which particles enter or leave the domain. + +Flux boundary conditions are commonly used to represent implantation from a plasma, outgassing to vacuum, permeation through a surface, or symmetry boundaries where no net transport occurs. + +#### Mathematical formulation + +The particle flux $\mathbf{J}$ is typically given by Fick’s law, + +$$ +\mathbf{J} = -D \nabla c, +$$ + +where $D$ is the diffusion coefficient and $c$ is the mobile concentration. + +On a boundary $\Gamma_N$, a prescribed normal flux is imposed as + +$$ +\mathbf{J} \cdot \mathbf{n} = g(\mathbf{x}, t) \quad \text{for } \mathbf{x} \in \Gamma_N, +$$ + +where: +- $\mathbf{n}$ is the outward unit normal vector, +- $g(\mathbf{x}, t)$ is the imposed particle flux (positive for outward flux, negative for inward flux). + +A special case is the **zero-flux (symmetry) boundary condition**, given by + +$$ +\mathbf{J} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma, +$$ + +which implies no net particle transport across the boundary. + +#### Weak formulation + +In the weak form, flux boundary conditions appear naturally as surface integrals after integration by parts. For a test function $v$, the boundary contribution is + +$$ +\int_{\Gamma_N} g(\mathbf{x}, t)\, v \, \mathrm{d}\Gamma. +$$ + +Because of this, flux boundary conditions are often referred to as *natural boundary conditions* and do not require explicit modification of the solution space, unlike Dirichlet conditions. + ++++ + +## Imposing fixed concentration boundary conditions + +Users can prescribe a fixed concentration on a boundary using `FixedConcentrationBC`, which can depend on temperature, time, and space. + +Let us consider a 2D, steady-state, single-species example with the following *space-dependent* boundary conditions: + +$$ \text{Top:} \quad c = 2x + y $$ +$$ \text{Right:} \quad c = y^2 + 1 $$ + +```{code-cell} ipython3 +:tags: [hide-input] + +import festim as F +import numpy as np +from dolfinx.mesh import create_unit_square +from mpi4py import MPI + +my_model = F.HydrogenTransportProblem() +my_model.mesh = F.Mesh(create_unit_square(MPI.COMM_WORLD, 10, 10)) +H = F.Species("H") +my_model.species = [H] +mat = F.Material(D_0=1e-3, E_D=0) +my_model.material = mat +my_model.temperature = 400 +my_model.settings = F.Settings(atol=1e-8, rtol=1e-8, transient=False) +``` + +We can define the top and right boundary conditions using `lambda` functions: + +```{code-cell} ipython3 +top_bc_function = lambda x: 2.0 * x[0] + x[1] +right_bc_function = lambda x: x[1] ** 2 + 1.0 +``` + +```{note} +`x[0]`, `x[1]`, `x[2]` corresponse to *(x, y, z)* in cartesian space. +``` + ++++ + +Now, let's create the boundaries to assign the BCs: + +```{code-cell} ipython3 +top_subdomain = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1)) +right_subdomain = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1)) +left_subdomain = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0)) +bottom_subdomain = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[1], 0)) +vol = F.VolumeSubdomain(id=5, material=mat) + +my_model.subdomains = [top_subdomain, right_subdomain, left_subdomain, bottom_subdomain, vol] +``` + +To include these boundary conditions to our problem, we use `FixedConcentrationBC`. We must also specify which subdomain (boundary) each BC is applied to, as well as the corresponding species: + +```{code-cell} ipython3 +top_bc = F.FixedConcentrationBC( + subdomain=top_subdomain, + value=top_bc_function, + species=H +) + +right_bc = F.FixedConcentrationBC( + subdomain=right_subdomain, + value=right_bc_function, + species=H +) +``` + +```{code-cell} ipython3 +my_model.boundary_conditions = [top_bc, right_bc] +my_model.initialise() +my_model.run() +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +import pyvista +from dolfinx import plot + +topology, cell_types, geometry = plot.vtk_mesh(H.post_processing_solution.function_space) +grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) +grid.point_data["c"] = H.post_processing_solution.x.array +grid.set_active_scalars("c") + +# pyvista.start_xvfb() +# pyvista.set_jupyter_backend("html") + +plotter = pyvista.Plotter() + +plotter.add_mesh(grid) +plotter.view_xy() +contours = grid.contour(isosurfaces=5, scalars="c") +plotter.add_mesh(contours, color="r", line_width=0.1) +if not pyvista.OFF_SCREEN: + plotter.show() +else: + figure = plotter.screenshot("concentration.png") +``` + +```{note} + +In this example, we did not define any boundary conditions for the left and bottom surfaces. If no BC is set on a boundary, it is assumed that the flux is null (symmetry boundary condition): + +$$ +\mathbf{J} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma +$$ + +``` + ++++ + +### Adding time and temperature dependent boundary conditions ### + ++++ + +Users can also impose concentration BCs that are dependent on space $x$, time $t$, and temperature $T$, such as: + +$$ c = 10 + x^2 + T(t+2) $$ + +To do so, we define a Python function: + +```{code-cell} ipython3 +def my_custom_value(x, t, T): + return 10 + x[0]**2 + T *(t + 2) + +boundary = F.SurfaceSubdomain(id=1) +my_bc = F.FixedConcentrationBC(subdomain=boundary, value=my_custom_value, species=H) +``` + +```{tip} +This is equivalent to: +` +my_custom_value = lambda x, t, T: 10 + x[0]**2 + T *(t + 2) + +` +``` + ++++ + +### Multi-species fixed concentration BC ## + +Users may also want to add BCs for multiple species in their model. For example, if our model included both deuterium and hydrogen and we wanted to impose separate fixed concentrations on a boundary, we need to define each species and then one BC for each: + +```{code-cell} ipython3 +import festim as F + +H = F.Species("H") +D = F.Species("D") +boundary = F.SurfaceSubdomain(id=1) + +top_H = F.FixedConcentrationBC( + subdomain=boundary, + value=5.0, + species=H +) + +top_D = F.FixedConcentrationBC( + subdomain=boundary, + value=10.0, + species=D +) +``` + +## Imposing particle flux boundary conditions + ++++ + +Users can impose the particle flux at boundaries using `ParticleFluxBC` class. At minimum, this BC requires a boundary, value, and species to be added: + +```{code-cell} ipython3 +import festim as F + +boundary = F.SurfaceSubdomain(id=1) +H = F.Species(name="H") + +my_flux_bc = F.ParticleFluxBC(subdomain=boundary, value=2, species=H) +``` + +To use this BC to your problem, simply add it to `boundary_conditions` attribute of your problem: + +```{code-cell} ipython3 +my_model = F.HydrogenTransportProblem() +my_model.boundary_conditions = [my_flux_bc] +``` + +### Species-dependent fluxes + +Similar to the concentration, the flux can be dependent on space, time and temperature. But for particle fluxes, the values can also be dependent on a species’ concentration. + +For example, let's define a hydrogen flux `J` that depends on the hydrogen concentration `c` and time `t`: + +$$ J(c,t) = 10t^2 + 2c $$ + +```{code-cell} ipython3 +import festim as F + +my_model = F.HydrogenTransportProblem() +boundary = F.SurfaceSubdomain(id=1) +H = F.Species(name="H") + +J = lambda t, c: 10*t**2 + 2*c +``` + +To add this BC, we need to create a dictionary that maps the concentration variable `c` in our custom function to our species `H`: + +```{code-cell} ipython3 +mapping_dict = {"c": H} +``` + +Now, we add our `ParticleFluxBC`, also utilizing the `species_dependent_value` argument: + +```{code-cell} ipython3 +my_flux_bc = F.ParticleFluxBC( + subdomain=boundary, + value=J, + species=H, + species_dependent_value=mapping_dict, +) + +my_model.boundary_conditions = [my_flux_bc] +``` + +```{note} +For multiple species, the dictionary + ++++ + +### Multi-species flux boundary conditions + +Users may also need to impose a flux boundary condition which depends on the concentrations of multiple species. + +Consider the following 1D example with three species, $\mathrm{A}$, $\mathrm{B}$, and $\mathrm{C}$. The boundary conditions include recombination fluxes $\phi_{AB}$ and $\phi_{ABC}$ that depend on the concentrations $\mathrm{c}$ (on the right) and fixed concentrations (on the left): + +$$ \mathrm{Right}: \phi_{AB} = -c_A c_B $$ + +$$ \mathrm{Right}: \phi_{ABC}= -10 c_A c_B c_C$$ + +$$ \mathrm{Left}: c_A = c_B = c_C = 1 $$ + +We must first define each species using `Species` and then create the dictionary to be passed into `species_dependent_value`. The dictionary maps each argument in the custom flux function to the corresponding defined species. We also define our custom functions for the fluxes: + +```{code-cell} ipython3 +import festim as F + +my_model = F.HydrogenTransportProblem() + +A = F.Species(name="A") +B = F.Species(name="B") +C = F.Species(name="C") + +my_model.species = [A, B, C] +species_dependent_value = {"c_A": A, "c_B": B, "c_C": C} + +recombination_flux_AB = lambda c_A, c_B, c_C: -c_A*c_B +recombination_flux_ABC = lambda c_A, c_B, c_C: -10*c_A*c_B*c_C +``` + +Now, we create our 1D mesh and assign boundary conditions (recombination on the right, fixed concentration on the left). + +The boundary condition `ParticleFluxBC` must be added for each species: + +```{code-cell} ipython3 +import numpy as np + +my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100)) + +D = 1e-2 +mat = F.Material( + D_0={A: 8*D, B: 7*D, C: D}, + E_D={A: 0.01, B: 0.01, C: 0.01}, +) + +bulk = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat) +left = F.SurfaceSubdomain1D(id=1, x=0) +right = F.SurfaceSubdomain1D(id=2, x=1) +my_model.subdomains = [bulk, left, right] + +my_model.boundary_conditions = [ + F.ParticleFluxBC( + subdomain=right, + value=recombination_flux_AB, + species=A, + species_dependent_value=species_dependent_value, + ), + F.ParticleFluxBC( + subdomain=right, + value=recombination_flux_AB, + species=B, + species_dependent_value=species_dependent_value, + ), + F.ParticleFluxBC( + subdomain=right, + value=recombination_flux_ABC, + species=A, + species_dependent_value=species_dependent_value, + ), + F.ParticleFluxBC( + subdomain=right, + value=recombination_flux_ABC, + species=B, + species_dependent_value=species_dependent_value, + ), + F.ParticleFluxBC( + subdomain=right, + value=recombination_flux_ABC, + species=C, + species_dependent_value=species_dependent_value, + ), + F.FixedConcentrationBC(subdomain=left,value=1,species=A), + F.FixedConcentrationBC(subdomain=left,value=1,species=B), + F.FixedConcentrationBC(subdomain=left,value=1,species=C), +] + +right_flux_A = F.SurfaceFlux(field=A,surface=right) +right_flux_B = F.SurfaceFlux(field=B,surface=right) +right_flux_C = F.SurfaceFlux(field=C,surface=right) + +my_model.exports = [ + right_flux_A, + right_flux_B, + right_flux_C, +] +``` + +```{note} +The diffusivity pre-factor `D_0` and activation energy `E_D` must be defined for each species in `Material`. Learn more about defining multi-species material properties [here](../material/material_advanced.md). +``` + ++++ + +Finally, let's solve and plot the profile for each species: + +```{code-cell} ipython3 +:tags: [hide-input] + +my_model.temperature = 300 +my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) + +my_model.initialise() +my_model.run() +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +import matplotlib.pyplot as plt + +def plot_profile(species, **kwargs): + c = species.post_processing_solution.x.array[:] + x = species.post_processing_solution.function_space.mesh.geometry.x[:,0] + return plt.plot(x, c, **kwargs) + +for species in my_model.species: + plot_profile(species, label=species.name) + +plt.xlabel('Position') +plt.ylabel('Concentration') +plt.legend() +plt.show() +``` + +We see the higher recombination flux for $\mathrm{ABC}$ decreases the concentration of $\mathrm{C}$ throughout the material. From 8eda15e95c6d581a0eb96129d13014b405f68b47 Mon Sep 17 00:00:00 2001 From: ck768 Date: Tue, 10 Feb 2026 13:36:23 -0500 Subject: [PATCH 07/17] fixes --- book/content/boundary_conditions/concentration.md | 4 ++-- book/content/boundary_conditions/fluxes.md | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/book/content/boundary_conditions/concentration.md b/book/content/boundary_conditions/concentration.md index f02b76ae..e6154277 100644 --- a/book/content/boundary_conditions/concentration.md +++ b/book/content/boundary_conditions/concentration.md @@ -1,11 +1,11 @@ --- jupytext: - formats: ipynb,md:myst + formats: md:myst,ipynb text_representation: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.18.1 + jupytext_version: 1.17.3 kernelspec: display_name: festim-workshop language: python diff --git a/book/content/boundary_conditions/fluxes.md b/book/content/boundary_conditions/fluxes.md index 493deddc..b5a6d9c9 100644 --- a/book/content/boundary_conditions/fluxes.md +++ b/book/content/boundary_conditions/fluxes.md @@ -1,11 +1,11 @@ --- jupytext: - formats: ipynb,md:myst + formats: md:myst,ipynb text_representation: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.18.1 + jupytext_version: 1.17.3 kernelspec: display_name: festim-workshop language: python From a7af26712e4d2ef2e40ec5dbab14b9696dbfd800 Mon Sep 17 00:00:00 2001 From: ck768 Date: Wed, 11 Feb 2026 08:47:37 -0500 Subject: [PATCH 08/17] moved to basic and advanced sections --- .../boundary_conditions/concentration.md | 20 ++---- .../h_transport_advanced.md | 11 +++ .../boundary_conditions/h_transport_basic.md | 72 +++++++++++-------- environment.yml | 27 +++---- 4 files changed, 74 insertions(+), 56 deletions(-) create mode 100644 book/content/boundary_conditions/h_transport_advanced.md diff --git a/book/content/boundary_conditions/concentration.md b/book/content/boundary_conditions/concentration.md index e6154277..fa54ccd3 100644 --- a/book/content/boundary_conditions/concentration.md +++ b/book/content/boundary_conditions/concentration.md @@ -1,6 +1,6 @@ --- jupytext: - formats: md:myst,ipynb + formats: ipynb,md:myst text_representation: extension: .md format_name: myst @@ -156,26 +156,20 @@ else: figure = plotter.screenshot("concentration.png") ``` - +In this example, we did not define any boundary conditions for the left and bottom surfaces. If no BC is set on a boundary, it is assumed that the flux is null (symmetry boundary condition): -+++ +$$ +\mathbf{J} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma +$$ -```{note} +``` +++ ### Time and temperature dependent boundary conditions ### -+++ - Users can also impose concentration BCs that are dependent on space, time and temperature, such as: $$ c = 10 + x^2 + T(t+2) $$ diff --git a/book/content/boundary_conditions/h_transport_advanced.md b/book/content/boundary_conditions/h_transport_advanced.md new file mode 100644 index 00000000..2f98edf0 --- /dev/null +++ b/book/content/boundary_conditions/h_transport_advanced.md @@ -0,0 +1,11 @@ +--- +jupytext: + formats: ipynb,md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.19.1 +--- + +# Hydrogen transport: advanced diff --git a/book/content/boundary_conditions/h_transport_basic.md b/book/content/boundary_conditions/h_transport_basic.md index e85dfc8e..c7a8703a 100644 --- a/book/content/boundary_conditions/h_transport_basic.md +++ b/book/content/boundary_conditions/h_transport_basic.md @@ -5,7 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.17.3 + jupytext_version: 1.19.1 kernelspec: display_name: festim-workshop language: python @@ -124,7 +124,7 @@ right_bc_function = lambda x: x[1] ** 2 + 1.0 ``` ```{note} -`x[0]`, `x[1]`, `x[2]` corresponse to *(x, y, z)* in cartesian space. +`x[0]`, `x[1]`, `x[2]` corresponse to *(x, y, z)* in cartesian space. Additionally, boundary conditions can also be given as Python functions (see *{ref}`label-time-temp-concentration`* to learn more.) ``` +++ @@ -157,6 +157,8 @@ right_bc = F.FixedConcentrationBC( ) ``` +We can add these to our problem, we pass them into `boundary_conditions` as a list and then run the simulation: + ```{code-cell} ipython3 my_model.boundary_conditions = [top_bc, right_bc] my_model.initialise() @@ -174,15 +176,14 @@ grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) grid.point_data["c"] = H.post_processing_solution.x.array grid.set_active_scalars("c") -# pyvista.start_xvfb() -# pyvista.set_jupyter_backend("html") +pyvista.start_xvfb() +pyvista.set_jupyter_backend("html") plotter = pyvista.Plotter() plotter.add_mesh(grid) plotter.view_xy() contours = grid.contour(isosurfaces=5, scalars="c") -plotter.add_mesh(contours, color="r", line_width=0.1) if not pyvista.OFF_SCREEN: plotter.show() else: @@ -201,6 +202,7 @@ $$ +++ +(label-time-temp-concentration)= ### Adding time and temperature dependent boundary conditions ### +++ @@ -220,18 +222,19 @@ my_bc = F.FixedConcentrationBC(subdomain=boundary, value=my_custom_value, specie ``` ```{tip} -This is equivalent to: +This custom function is equivalent to: + ` my_custom_value = lambda x, t, T: 10 + x[0]**2 + T *(t + 2) - ` -``` +++ ### Multi-species fixed concentration BC ## -Users may also want to add BCs for multiple species in their model. For example, if our model included both deuterium and hydrogen and we wanted to impose separate fixed concentrations on a boundary, we need to define each species and then one BC for each: +Users may also want to add BCs for multiple species in their model. + +For example, if we wanted to impose separate fixed concentrations for deuterium and hydrogen on a boundary, we need to define each species and then one BC for each: ```{code-cell} ipython3 import festim as F @@ -312,24 +315,20 @@ my_flux_bc = F.ParticleFluxBC( my_model.boundary_conditions = [my_flux_bc] ``` -```{note} -For multiple species, the dictionary - -+++ - ### Multi-species flux boundary conditions Users may also need to impose a flux boundary condition which depends on the concentrations of multiple species. Consider the following 1D example with three species, $\mathrm{A}$, $\mathrm{B}$, and $\mathrm{C}$. The boundary conditions include recombination fluxes $\phi_{AB}$ and $\phi_{ABC}$ that depend on the concentrations $\mathrm{c}$ (on the right) and fixed concentrations (on the left): -$$ \mathrm{Right}: \phi_{AB} = -c_A c_B $$ +$$ \text{Right (recombination flux):} \quad \phi_{AB} = -c_A c_B $$ + +$$ \text{Right (recombination flux):} \quad \phi_{ABC} = -10 c_A c_B c_C$$ -$$ \mathrm{Right}: \phi_{ABC}= -10 c_A c_B c_C$$ +$$ \text{Left (fixed concentration):} \quad c_A = c_B = c_C = 1 $$ -$$ \mathrm{Left}: c_A = c_B = c_C = 1 $$ -We must first define each species using `Species` and then create the dictionary to be passed into `species_dependent_value`. The dictionary maps each argument in the custom flux function to the corresponding defined species. We also define our custom functions for the fluxes: +First, let us define each species using `Species`: ```{code-cell} ipython3 import festim as F @@ -341,15 +340,22 @@ B = F.Species(name="B") C = F.Species(name="C") my_model.species = [A, B, C] +``` + +Now, we create the dictionary to be passed into `species_dependent_value`; this maps each argument in the custom flux function to the corresponding defined species: + +```{code-cell} ipython3 species_dependent_value = {"c_A": A, "c_B": B, "c_C": C} +``` +We also define our custom functions for the fluxes: + +```{code-cell} ipython3 recombination_flux_AB = lambda c_A, c_B, c_C: -c_A*c_B recombination_flux_ABC = lambda c_A, c_B, c_C: -10*c_A*c_B*c_C ``` -Now, we create our 1D mesh and assign boundary conditions (recombination on the right, fixed concentration on the left). - -The boundary condition `ParticleFluxBC` must be added for each species: +Now, we create our 1D mesh and corresponding boundaries: ```{code-cell} ipython3 import numpy as np @@ -366,7 +372,17 @@ bulk = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat) left = F.SurfaceSubdomain1D(id=1, x=0) right = F.SurfaceSubdomain1D(id=2, x=1) my_model.subdomains = [bulk, left, right] +``` +```{note} +The diffusivity pre-factor `D_0` and activation energy `E_D` must be defined for each species in `Material`. Learn more about defining multi-species material properties [here](../material/material_advanced.md). Here we choose arbitrarily different diffusivity coefficients to see transport between species. +``` + ++++ + +Let's assign boundary conditions (recombination on the right, fixed concentration on the left). The boundary condition `ParticleFluxBC` must be added for each species: + +```{code-cell} ipython3 my_model.boundary_conditions = [ F.ParticleFluxBC( subdomain=right, @@ -402,7 +418,11 @@ my_model.boundary_conditions = [ F.FixedConcentrationBC(subdomain=left,value=1,species=B), F.FixedConcentrationBC(subdomain=left,value=1,species=C), ] +``` +We can export the flux for each species on the right using `SurfaceFlux` (see [post-processing derived values](../post_process/derived.md) to learn more about exporting fluxes): + +```{code-cell} ipython3 right_flux_A = F.SurfaceFlux(field=A,surface=right) right_flux_B = F.SurfaceFlux(field=B,surface=right) right_flux_C = F.SurfaceFlux(field=C,surface=right) @@ -414,12 +434,6 @@ my_model.exports = [ ] ``` -```{note} -The diffusivity pre-factor `D_0` and activation energy `E_D` must be defined for each species in `Material`. Learn more about defining multi-species material properties [here](../material/material_advanced.md). -``` - -+++ - Finally, let's solve and plot the profile for each species: ```{code-cell} ipython3 @@ -430,10 +444,6 @@ my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) my_model.initialise() my_model.run() -``` - -```{code-cell} ipython3 -:tags: [hide-input] import matplotlib.pyplot as plt diff --git a/environment.yml b/environment.yml index 469ed3be..b23d27f2 100644 --- a/environment.yml +++ b/environment.yml @@ -2,22 +2,25 @@ name: festim-workshop channels: - conda-forge - defaults + dependencies: + - python<3.14 + - pip=23.3 + - setuptools<82 + - wheel<0.45 - jupyter-book<2.0 - - python<3.14 # som pyvista incompatibility issues with python 3.14 - - jupytext # necessary to open MyST files using BinderHub https://jupyterbook.org/en/stable/interactive/launchbuttons.html#launchbuttons-binder - - sphinx-tags - - matplotlib - - meshio[all] - - pip>=20.1 + - jupytext - ipykernel - nbconvert + - sphinx-tags + - matplotlib + - scipy + - shapely + - meshio - fenics-dolfinx - - festim=2.0b1 - python-gmsh - - shapely - - scipy + - festim=2.0b1 - pip: - - pyparsing - - h-transport-materials==0.16 - - pyvista[all,trame] # trame cannot be installed on conda \ No newline at end of file + - h-transport-materials==0.16 + - pyparsing + - pyvista[all,trame]<0.47.0 From 38dde951fe687625e7ab8b7ca33035b14ca29394 Mon Sep 17 00:00:00 2001 From: ck768 Date: Thu, 12 Feb 2026 15:31:52 -0500 Subject: [PATCH 09/17] separating examples to new section --- book/_toc.yml | 2 + book/content/boundary_conditions/examples.md | 13 + .../h_transport_advanced.md | 297 ++++++++++++++++++ 3 files changed, 312 insertions(+) create mode 100644 book/content/boundary_conditions/examples.md diff --git a/book/_toc.yml b/book/_toc.yml index 4e99a860..cfa3ebc7 100644 --- a/book/_toc.yml +++ b/book/_toc.yml @@ -24,6 +24,8 @@ parts: - caption: Boundary Conditions chapters: - file: content/boundary_conditions/h_transport_basic + - file: content/boundary_conditions/h_transport_advanced + - file: content/boundary_conditions/examples - file: content/boundary_conditions/heat_transfer - caption: Species & Reactions diff --git a/book/content/boundary_conditions/examples.md b/book/content/boundary_conditions/examples.md new file mode 100644 index 00000000..57e9f2e7 --- /dev/null +++ b/book/content/boundary_conditions/examples.md @@ -0,0 +1,13 @@ +--- +jupytext: + formats: ipynb,md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.19.1 +--- + +# Boundary condition: examples + ++++ diff --git a/book/content/boundary_conditions/h_transport_advanced.md b/book/content/boundary_conditions/h_transport_advanced.md index 2f98edf0..c445a045 100644 --- a/book/content/boundary_conditions/h_transport_advanced.md +++ b/book/content/boundary_conditions/h_transport_advanced.md @@ -6,6 +6,303 @@ jupytext: format_name: myst format_version: 0.13 jupytext_version: 1.19.1 +kernelspec: + display_name: festim-workshop + language: python + name: python3 --- # Hydrogen transport: advanced + +This section discusses how to implement advanced boundary conditions (BCs) for hydrogen transport problems. + +Objectives: +* Learn mathematical formulation for solubility laws and surface reactions +* Choose solubility laws (Henry's or Siervert's) +* Implement surface reaction boundary conditions +* Model isotopic exchange as a recombination flux + ++++ + +## Learn mathematical formulation for solubility laws and surface reactions + ++++ + +### Solubility laws + +A solubility law prescribes the equilibrium relationship between the hydrogen concentration in a solid and the surrounding environment. Two common solubility laws for hydrogen are **Henry’s law** and **Sieverts’ law**: + +#### Henry's law + +Henry’s law assumes a linear relationship between the hydrogen concentration in the solid and the partial pressure of hydrogen in the gas: + +$$ +c_s = k_H \, P_H +$$ + +where: +- $c_s$ is the surface concentration, +- $P_H$ is the hydrogen partial pressure in the surrounding environment, +- $k_H$ is the Henry’s law constant (material and temperature-dependent). + +```{tip} +Users interested in molten salt interactions with hydrogen will usually need to use *Henry's Law*. +``` + +#### Sieverts' law + +Sieverts’ law applies when hydrogen dissolves in metals as atomic hydrogen. The solubility is proportional to the square root of the hydrogen partial pressure (due to diatomic hydrogen dissociating into atoms in the solid): + +$$ +c_s = k_S \, \sqrt{P_H} +$$ + +where: +- $k_S$ is the Sieverts’ constant (which depends on temperature and the metal) + +```{tip} +Users interested in liquid metal interactions with hydrogen will usually need to use **Sievert's Law**. +``` + +#### Weak formulation + +In the weak form, solubility is imposed as a **Robin-type boundary condition**: + +$$ +\int_{\Gamma_s} k_s \left(c - c_s \right) v \, \mathrm{d}\Gamma +$$ + +where $k_s$ is a surface reaction rate constant controlling the exchange between bulk and surface concentration. This allows modeling of adsorption/desorption and surface-limited transport at material boundaries. + +--- + ++++ + +### Simple surface reaction + +Surface reactions between adsorbed species and gas-phase products can be represented using the `SurfaceReactionBC` class. This boundary condition defines a reversible reaction of the form: + +$$ +A + B \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad C +$$ + +where: +- $\mathrm{A}, \mathrm{B}$ are surface reactants +- $\mathrm{C}$ is the product species in the gas phase +- $\mathrm{K}_r$ and $\mathrm{K}_d$ are forward and backward reaction rates, respectively + +#### Mathematical formulation + +Let: +- $c_A, c_B$ be the surface concentrations of $\mathrm{A}$ and $\mathrm{B}$ +- $P_C$ be the partial pressure of $\mathrm{C}$ +- $T$ the surface temperature +- $k_B$ the Boltzmann constant + +The forward and backward reaction rates follow Arrhenius laws: + +$$ +K_r = k_{r0} \exp\!\left(-\frac{E_{kr}}{k_B T}\right), \qquad +K_d = k_{d0} \exp\!\left(-\frac{E_{kd}}{k_B T}\right) +$$ + +The net surface reaction rate is: + +$$ +K = K_r \, c_A c_B - K_d \, P_C +$$ + +The resulting flux of species $\mathrm{A}$ into the surface is: + +$$ +\mathbf{J}_A \cdot \mathbf{n} = K +$$ + +If $\mathrm{A=B}$, the total flux becomes: + +$$ +\mathbf{J}_A \cdot \mathbf{n} = 2 K +$$ + +where $\mathbf{n}$ is the outward normal vector at the surface. + +#### Weak form contribution + +In the weak formulation, the surface reaction contributes as a Robin-type term: + +$$ +\int_{\Gamma_s} \mathbf{J}_A \cdot \mathbf{n} \, v \, d\Gamma = +\int_{\Gamma_s} K \, v \, d\Gamma +$$ + +where $v$ is a test function. + +--- + ++++ + +### Recombination and dissociation + +Hydrogen recombination or dissociation can be modeled with `SurfaceReactionBC`, e.g.: + +$$ +\mathrm{H + H} \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad \mathrm{H_2} +$$ + +#### Mathematical formulation + +Let: +- $c_H$ be the surface concentration of atomic hydrogen +- $P_{H_2}$ the partial pressure of molecular hydrogen + +The net reaction rate is: + +$$ +K = K_r \, c_H^2 - K_d \, P_{H_2} +$$ + +The corresponding flux of atomic hydrogen is: + +$$ +\mathbf{J}_H \cdot \mathbf{n} = - 2 K +$$ + +#### Weak form contribution + +The weak form contribution of recombination flux is: + +$$ +\int_{\Gamma_s} \mathbf{J}_H \cdot \mathbf{n} \, v \, d\Gamma = +- \int_{\Gamma_s} 2 K \, v \, d\Gamma +$$ + +--- + ++++ + +### Isotopic exchange + +Isotopic exchange occurs when hydrogenic isotopes swap positions, for example: + +$$ +\mathrm{T + H_2} \rightleftharpoons \mathrm{H} + \mathrm{HT} +$$ + +#### Mathematical formulation + +Let: +- $c_T, c_{H_2}, c_{HT}$ be the surface concentrations of tritium, molecular hydrogen, and hydrogen-tritium +- $K_{r0}, K_{r0}^\ast$ the pre-exponential factors +- $E_{Kr}, E_{Kr}^\ast$ the activation energies + +If $c_{H_2} \gg c_{HT}$, the flux of tritium is approximated as: + +$$ +\phi_T = +- K_{r0} \exp\!\left(-\frac{E_{Kr}}{k_B T}\right) c_T^2 +- K_{r0}^\ast \exp\!\left(-\frac{E_{Kr}^\ast}{k_B T}\right) c_{H_2} c_T +$$ + +#### Weak form contribution + +The weak form contribution for tritium flux is: + +$$ +\int_{\Gamma_s} \phi_T \, v \, d\Gamma = +- \int_{\Gamma_s} \Bigg[ +K_{r0} \exp\!\left(-\frac{E_{Kr}}{k_B T}\right) c_T^2 ++ K_{r0}^\ast \exp\!\left(-\frac{E_{Kr}^\ast}{k_B T}\right) c_{H_2} c_T +\Bigg] v \, d\Gamma +$$ + +These fluxes can be implemented in FESTIM using `ParticleFluxBC` with user-defined expressions for each reaction term (see [here](examples.ipynb) to learn more). + ++++ + +## Implementing simple surface reaction boundary conditions + +Users can impose a surface reaction at boundary using `SurfaceReactionBC`. + +First, create a boundary to impose the BC and the reactant species: + +```{code-cell} ipython3 +import festim as F + +boundary = F.SurfaceSubdomain(id=1) +A = F.Species("A") +B = F.Species("B") +``` + +Now, we add these as arguments to the `SurfaceReactionBC` class. We'll also need to assign a `gas_pressure` (corresponding to the partial pressure of the product species), and the corresponding rate (`k_r0`, `k_d0`) and energy (`E_kr`, `E_dr`) coefficients (see above to learn more about these parameters): + +```{code-cell} ipython3 +my_bc = F.SurfaceReactionBC( + reactant=[A, B], + gas_pressure=1e5, + k_r0=1, + E_kr=0.1, + k_d0=0, + E_kd=0, + subdomain=boundary, +) +``` + +Finally, add the BC to your problem's `boundary_conditions` attribute using a list: + +```{code-cell} ipython3 +my_model = F.HydrogenTransportProblem() +my_model.boundary_conditions = [my_bc] +``` + +### Modeling recombination and dissociation + +Recombination and dissociation can also be modeled using `SurfaceReactionBC`, where the forward and backward rates of this reaction correspond to recombination and dissociation, respectively. + +To model the reaction: + +$$ \mathrm{H} + \mathrm{H} \rightleftharpoons \mathrm{H_2}$$ + +where $ \text{Species A} = \text{Species B} = \text{H} $, assign your `reactants` list accordingly: + +```{code-cell} ipython3 +H = F.Species("H") + +my_recombination_bc = F.SurfaceReactionBC( + reactant=[A, A], + gas_pressure=1e5, + k_r0=1, + E_kr=0.1, + k_d0=1e-5, + E_kd=0.1, + subdomain=boundary, +) +``` + +## Modeling isotopic exchange as a recombination flux + +Isotopic exchange reactions can be modeled as user-defined expressions using `ufl`. + +Assuuming a large, constant concentration of molecular hydrogen $H_2$ at a boundary, we can define a recombination flux using rate and energy coefficients: + +```{code-cell} ipython3 +import ufl +import festim as F + +Kr_0 = 1.0 +E_Kr = 0.1 + +def my_custom_recombination_flux(c, T): + Kr_0_custom = 1.0 + E_Kr_custom = 0.5 # eV + h2_conc = 1e25 # assumed constant H2 concentration in + + recombination_flux = ( + -(Kr_0 * ufl.exp(-E_Kr / (F.k_B * T))) * c**2 + - (Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))) * h2_conc * c + ) + return recombination_flux +``` + +```{note} +For more complex isotopic exchanges, we can also use `SurfaceReactionBC` to define recombination fluxes. See [NOTE] to learn more. +``` From cb63c6eba2c2576a8a6e50c2ed53d23b60c68811 Mon Sep 17 00:00:00 2001 From: Chirag Date: Fri, 13 Feb 2026 16:47:26 -0500 Subject: [PATCH 10/17] first pass at bcs with more integrated examples --- book/content/boundary_conditions/examples.md | 290 ++++++++++++++++- .../h_transport_advanced.md | 306 ++++++++++++++---- .../boundary_conditions/h_transport_basic.md | 3 +- 3 files changed, 537 insertions(+), 62 deletions(-) diff --git a/book/content/boundary_conditions/examples.md b/book/content/boundary_conditions/examples.md index 57e9f2e7..bde78add 100644 --- a/book/content/boundary_conditions/examples.md +++ b/book/content/boundary_conditions/examples.md @@ -6,8 +6,296 @@ jupytext: format_name: myst format_version: 0.13 jupytext_version: 1.19.1 +kernelspec: + display_name: festim-workshop + language: python + name: python3 --- -# Boundary condition: examples +# Advanced examples + +This section discusses a few advanced examples using various boundary conditions to help users understand how BCs can help model their systems. + +Objective: +* Model plasma implantation using volumetric sources and approximations +* Model complex multi-species isotopic exchange surface reactions + ++++ + +## Modeling plasma implantation + +We can model plasma implantation using FESTIM's `ParticleSource` class, which is a class used to define volumetric source terms. This is helpful in modeling thermo-desorption spectra (TDS) experiments or including the effect of plasma exposure on hydrogen transport. + +Consider the following 1D plasma implantation problem, where we represent the plasma as a hydrogen source $S_{ext}$ implanted on a tungsten material that is 20 microns thick: + +$$ S_{ext} = \varphi_{imp} \cdot f(x) $$ +$$ L = 20 \mu \mathrm{m}$$ +$$\varphi_{imp} = 10^{13} \quad \mathrm{m}^{-2}\mathrm{s}^{-1}$$ + +where $\varphi_{imp}$ is the implantation flux and $f(x)$ is a Gaussian spatial distribution (distribution mean value represents implantation depth). + +First, we setup a 1D mesh ranging from $ [0,L] $ and assign the subdomains and material properties for tungsten: + +```{code-cell} ipython3 +import festim as F +import ufl +import numpy as np + +L = 20e-6 +my_model = F.HydrogenTransportProblem() +vertices = np.linspace(0, L, 1000) +my_model.mesh = F.Mesh1D(vertices) + +tungsten = F.Material( + D_0=4.1e-07, # m2/s + E_D=0.39, # eV +) +volume_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=tungsten) +left_boundary = F.SurfaceSubdomain1D(id=1, x=0) +right_boundary = F.SurfaceSubdomain1D(id=2, x=L) + +my_model.subdomains = [ + volume_subdomain, + left_boundary, + right_boundary, +] +``` + +Now, we define our `incident_flux` and `gaussian_distribution` function. Here, we define the mean implantation depth $R_p$ and distribution width to be a few nanometers long: + +$$ R_p = 4 \times 10^{-9} \mathrm{m} $$ +$$ \mathrm{width} = 2.5 \times 10^{-9} \mathrm{m} $$ + +```{code-cell} ipython3 +incident_flux = 1e13 +Rp = 4e-9 +width = 2.5e-9 + +def gaussian_distribution(x, center, width): + return ( + 1 + / (width * (2 * ufl.pi) ** 0.5) + * ufl.exp(-0.5 * ((x[0] - center) / width) ** 2) + ) +``` + +We can define our species and use `ParticleSource` to represent the source term, and then add it to our model: + +```{code-cell} ipython3 +H = F.Species("H") +my_model.species = [H] + +source_term = F.ParticleSource( + value=lambda x: incident_flux * gaussian_distribution(x, Rp, width), + volume=volume_subdomain, + species=H, +) + +my_model.sources = [source_term] +``` + +Finally, we assign boundary conditions (zero concentration at $x=0$ and $x=L$) and solve our problem: + +```{code-cell} ipython3 +my_model.boundary_conditions = [ + F.FixedConcentrationBC(subdomain=left_boundary, value=0, species=H), + F.FixedConcentrationBC(subdomain=right_boundary, value=0, species=H), +] + +my_model.temperature = 400 +my_model.settings = F.Settings(atol=1e10, rtol=1e-10, transient=False) + +profile_export = F.Profile1DExport(field=H,subdomain=volume_subdomain) +my_model.exports = [profile_export] + +my_model.initialise() +my_model.run() +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +import matplotlib.pyplot as plt + +x = my_model.exports[0].x +c = my_model.exports[0].data[0] + +plt.plot(x, c) +plt.xlabel("x") +plt.ylabel("c") +plt.show() +``` + +We see that a huge spike in concentration in the first few nanometers of tungesten, where the implantation range is focused. +++ + +### Approximating plasma implantation using fixed concentration boundary conditions + +If recombination is fast enough, the spike shown above can be approximated as a fixed concentration boundary condition that mainly drives diffusion across the material. Learn more about the plasma implantation approximation approach _[here](https://festim.readthedocs.io/en/fenicsx/theory.html#plasma-implantation-approximation)_. + +To see how we might approximate this, let's define a maximum concentration to set on the left boundary, representing the spike from the implantation: + +$$ c_m = \frac{R_p \cdot \varphi_{imp}}{D} $$ + +where $\varphi_{imp}$ is the implantation flux and $R_p$ is the implantation depth, both which we defined earlier, and $D$ is the material diffusivity: + +```{code-cell} ipython3 +D = tungsten.D_0 * np.exp(-tungsten.E_D / (F.k_B * my_model.temperature)) +c_m = incident_flux * Rp / D +``` + +Now, we'll change our boundary conditions to represent the implantation as a fixed concentration on the left boundary, and remove the source term from our problem: + +```{code-cell} ipython3 +my_model.boundary_conditions = [ + F.FixedConcentrationBC(subdomain=left_boundary, value=c_m, species=H), + F.FixedConcentrationBC(subdomain=right_boundary, value=0, species=H), +] + +my_model.sources = [] +profile_export = F.Profile1DExport(field=H,subdomain=volume_subdomain) +my_model.exports = [profile_export] + +my_model.initialise() +my_model.run() +``` + +Let's compare the profiles from the approximation to the volumetric source results: + +```{code-cell} ipython3 +:tags: [hide-input] + +import matplotlib.pyplot as plt + +x2 = my_model.exports[0].x +c2 = my_model.exports[0].data[0] + +plt.plot(x, c, label="With volumetric source") +plt.plot(x2, c2, label="Approximation") + +plt.xlabel("x") +plt.ylabel("c") +plt.legend() +plt.show() +``` + +Using the approximation is computationally less expensive, and still provides similar diffusion profiles. + ++++ + +## Complex isotopic exchange with multple hydrogenic species ## + +Surface reactions can involve multiple hydrogen isotopes, allowing for the modeling of complex isotope-exchange mechanisms between species. For example, in a system with both mobile hydrogen and deuteriun, various molecular recombination pathways may occur at the surface, resulting in the formation of $H_2$, $D_2$, and $HD$: + +$$ \text{Reaction 1}: \mathrm{H+D} \rightarrow \mathrm{HD} \longrightarrow \phi_1 = K_{r1} c_H c_D - K_{d1}P_{HD} $$ +$$ \text{Reaction 2}: \mathrm{D+D} \rightarrow \mathrm{D_2} \longrightarrow \phi_2 = 2K_{r2} c_D^2 - K_{d2}P_{D2} $$ +$$ \text{Reaction 3}: \mathrm{D+H_2} \rightarrow \mathrm{HD + H} \longrightarrow \phi_3 = K_{r3} c_H c_D - K_{d3}P_{HD} $$ +$$ \text{Reaction 4}: \mathrm{H+H} \rightarrow \mathrm{H_2} \longrightarrow \phi_4 = 2K_{r4} c_H^2 - K_{d4}P_{H2} $$ + +Now consider the case where deuterium diffuses from left to right and reacts with background +$\mathrm{H_2}$, while $\mathrm{P_{HD}}$ and $\mathrm{P_{D_2}}$ are negligible at the surface. +Formation of $\mathrm{H}$ at the right boundary induces back-diffusion toward the left, +even though none existed initially. + +The boundary conditions for this scenario are: + +$$ +c_D(x=0) = 1, \qquad c_H(x=0) = 0, \qquad P_{H2}(x=1) = \text{1000 Pa} +$$ + +First, let's define a 1D mesh ranging from $\mathrm{x=[0,1]}$: + +```{code-cell} ipython3 +import numpy as np +import festim as F + +my_model = F.HydrogenTransportProblem() +my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 100)) + +left_surf = F.SurfaceSubdomain1D(id=1, x=0) +right_surf = F.SurfaceSubdomain1D(id=2, x=1) + +material = F.Material(D_0=1e-2, E_D=0) +vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=material) + +my_model.subdomains = [vol, left_surf, right_surf] +``` + +Now, we define our species at recombination reactions using `SurfaceReactionBC`: + +```{code-cell} ipython3 +H = F.Species("H") +D = F.Species("D") +my_model.species = [H, D] + +H2 = F.SurfaceReactionBC( + reactant=[H, H], + gas_pressure=100000, + k_r0=1, + E_kr=0.1, + k_d0=1e-5, + E_kd=0.1, + subdomain=right_surf, +) + +HD = F.SurfaceReactionBC( + reactant=[H, D], + gas_pressure=0, + k_r0=1, + E_kr=0.1, + k_d0=1e-5, + E_kd=0.1, + subdomain=right_surf, +) + +D2 = F.SurfaceReactionBC( + reactant=[D, D], + gas_pressure=0, + k_r0=1, + E_kr=0.1, + k_d0=1e-5, + E_kd=0.1, + subdomain=right_surf, +) +``` + +Finally, we add our boundary conditions and solve the steady-state problem: + +```{code-cell} ipython3 +my_model.boundary_conditions = [ + H2, + D2, + HD, + F.FixedConcentrationBC(subdomain=left_surf, value=1, species=D), + F.FixedConcentrationBC(subdomain=left_surf, value=0, species=H), +] + +my_model.temperature = 300 + +my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) + +my_model.initialise() +my_model.run() +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +import matplotlib.pyplot as plt + +def plot_profile(species, **kwargs): + c = species.post_processing_solution.x.array[:] + x = species.post_processing_solution.function_space.mesh.geometry.x[:,0] + return plt.plot(x, c, **kwargs) + +for species in my_model.species: + plot_profile(species, label=species.name) + +plt.xlabel('Position') +plt.ylabel('Concentration') +plt.legend() +plt.show() +``` + +We see that the background $\mathrm{H_2}$ reacts with the $\mathrm{D}$, removing the total amount of $\mathrm{D}$ from the surface. Conversely, the $\mathrm{H}$ diffuses from the surface towards the left. diff --git a/book/content/boundary_conditions/h_transport_advanced.md b/book/content/boundary_conditions/h_transport_advanced.md index c445a045..84f8ffd7 100644 --- a/book/content/boundary_conditions/h_transport_advanced.md +++ b/book/content/boundary_conditions/h_transport_advanced.md @@ -17,14 +17,13 @@ kernelspec: This section discusses how to implement advanced boundary conditions (BCs) for hydrogen transport problems. Objectives: -* Learn mathematical formulation for solubility laws and surface reactions * Choose solubility laws (Henry's or Siervert's) * Implement surface reaction boundary conditions * Model isotopic exchange as a recombination flux +++ -## Learn mathematical formulation for solubility laws and surface reactions +## Understanding background for solubility laws and surface reactions +++ @@ -74,11 +73,48 @@ $$ where $k_s$ is a surface reaction rate constant controlling the exchange between bulk and surface concentration. This allows modeling of adsorption/desorption and surface-limited transport at material boundaries. ++++ + +### Imposing a solubility law + +Users can define the surface concentration using either Sieverts’ law. For Sieverts' law of solubility, we can use `SievertsBC`: + +```{code-cell} ipython3 +import festim as F + +boundary = F.SurfaceSubdomain(id=1) +H = F.Species(name="Hydrogen") + +custom_pressure_value = lambda t: 2 + t +my_sieverts_bc = F.SievertsBC(subdomain=3, S_0=2, E_S=0.1, species=H, pressure=custom_pressure_value) +``` + +Similarly, for Henry's law of solubility, we can use `HenrysBC`: + +```{code-cell} ipython3 +pressure_value = lambda t: 5 * t +my_henrys_bc = F.HenrysBC(subdomain=3, H_0=1.5, E_H=0.2, species=H, pressure=pressure_value) +``` + +To use these BCs in your simulation, add them to `boundary_conditions`: + +```{code-cell} ipython3 +my_model = F.HydrogenTransportProblem() +my_model.boundary_conditions = [my_sieverts_bc, my_henrys_bc] +``` + +```{note} +Using `HenrysBC` or `SievertsBC` is just a convenient way for users to define a `FixedConcentrationBC`, using the solubility law to calculate the concentration. +``` + ++++ + + --- +++ -### Simple surface reaction +## Simple surface reaction Surface reactions between adsorbed species and gas-phase products can be represented using the `SurfaceReactionBC` class. This boundary condition defines a reversible reaction of the form: @@ -91,7 +127,7 @@ where: - $\mathrm{C}$ is the product species in the gas phase - $\mathrm{K}_r$ and $\mathrm{K}_d$ are forward and backward reaction rates, respectively -#### Mathematical formulation +### Mathematical formulation Let: - $c_A, c_B$ be the surface concentrations of $\mathrm{A}$ and $\mathrm{B}$ @@ -137,11 +173,48 @@ $$ where $v$ is a test function. ++++ + +### Implementing surface reaction boundary conditions + +Users can impose a surface reaction on a boundary using `SurfaceReactionBC`. + +First, create a boundary to impose the BC and the reactant species: + +```{code-cell} ipython3 +import festim as F + +boundary = F.SurfaceSubdomain(id=1) +A = F.Species("A") +B = F.Species("B") +``` + +Now, we add these as arguments to the `SurfaceReactionBC` class. We'll also need to assign a `gas_pressure` (corresponding to the partial pressure of the product species), and the forward/backward rate (`k_r0`, `k_d0`) and energy (`E_kr`, `E_dr`) coefficients (see above to learn more about these parameters): + +```{code-cell} ipython3 +my_bc = F.SurfaceReactionBC( + reactant=[A, B], + gas_pressure=1e5, + k_r0=1, + E_kr=0.1, + k_d0=0, + E_kd=0, + subdomain=boundary, +) +``` + +Finally, add the BC to your problem's `boundary_conditions` attribute using a list: + +```{code-cell} ipython3 +my_model = F.HydrogenTransportProblem() +my_model.boundary_conditions = [my_bc] +``` + --- +++ -### Recombination and dissociation +## Recombination and dissociation Hydrogen recombination or dissociation can be modeled with `SurfaceReactionBC`, e.g.: @@ -149,7 +222,7 @@ $$ \mathrm{H + H} \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad \mathrm{H_2} $$ -#### Mathematical formulation +### Mathematical formulation Let: - $c_H$ be the surface concentration of atomic hydrogen @@ -176,11 +249,37 @@ $$ - \int_{\Gamma_s} 2 K \, v \, d\Gamma $$ ++++ + --- +++ -### Isotopic exchange +### Modeling recombination and dissociation + +Recombination and dissociation can also be modeled using `SurfaceReactionBC`, where the forward and backward rates of this reaction correspond to recombination and dissociation, respectively. + +To model the reaction: + +$$ \mathrm{H} + \mathrm{H} \rightleftharpoons \mathrm{H_2}$$ + +where $ \text{Species A} = \text{Species B} = \text{H} $, assign your `reactants` list accordingly: + +```{code-cell} ipython3 +H = F.Species("H") + +my_recombination_bc = F.SurfaceReactionBC( + reactant=[A, A], + gas_pressure=1e5, + k_r0=1, + E_kr=0.1, + k_d0=1e-5, + E_kd=0.1, + subdomain=boundary, +) +``` + +## Isotopic exchange Isotopic exchange occurs when hydrogenic isotopes swap positions, for example: @@ -188,7 +287,7 @@ $$ \mathrm{T + H_2} \rightleftharpoons \mathrm{H} + \mathrm{HT} $$ -#### Mathematical formulation +### Mathematical formulation Let: - $c_T, c_{H_2}, c_{HT}$ be the surface concentrations of tritium, molecular hydrogen, and hydrogen-tritium @@ -219,90 +318,179 @@ These fluxes can be implemented in FESTIM using `ParticleFluxBC` with user-defin +++ -## Implementing simple surface reaction boundary conditions +### Modeling isotopic exchange as a recombination flux -Users can impose a surface reaction at boundary using `SurfaceReactionBC`. +Isotopic exchange reactions can be modeled as user-defined expressions using `ufl`. -First, create a boundary to impose the BC and the reactant species: +Assuming a large, constant concentration of molecular hydrogen $H_2$ at a boundary, we can define a recombination flux using rate and energy coefficients: ```{code-cell} ipython3 +import ufl import festim as F -boundary = F.SurfaceSubdomain(id=1) -A = F.Species("A") -B = F.Species("B") +Kr_0 = 1.0 +E_Kr = 0.1 + +def my_custom_recombination_flux(c, T): + Kr_0_custom = 1.0 + E_Kr_custom = 0.5 # eV + h2_conc = 1e25 # assumed constant H2 concentration in + + recombination_flux = ( + -(Kr_0 * ufl.exp(-E_Kr / (F.k_B * T))) * c**2 + - (Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))) * h2_conc * c + ) + return recombination_flux +``` + +```{note} +For more complex isotopic exchanges, we can also use `SurfaceReactionBC` to define recombination fluxes. See [modeling isotopic exchange for multiple species](examples.md) to learn more. ``` -Now, we add these as arguments to the `SurfaceReactionBC` class. We'll also need to assign a `gas_pressure` (corresponding to the partial pressure of the product species), and the corresponding rate (`k_r0`, `k_d0`) and energy (`E_kr`, `E_dr`) coefficients (see above to learn more about these parameters): ++++ + +Let's work through a 1D example to illustrate the effect of isotopic exchange. We'll consider tritium diffusion from left to right, with a large partial pressure of $H_2$ at the right boundary (where isotopic exchange can occur). + +First, we'll set up our mesh and materials: ```{code-cell} ipython3 -my_bc = F.SurfaceReactionBC( - reactant=[A, B], - gas_pressure=1e5, - k_r0=1, - E_kr=0.1, - k_d0=0, - E_kd=0, - subdomain=boundary, +import numpy as np +from dolfinx.mesh import create_unit_square +from mpi4py import MPI + +my_model = F.HydrogenTransportProblem() +my_model.mesh = F.Mesh(create_unit_square(MPI.COMM_WORLD, 300, 300)) + +right_subdomain = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1)) +left_subdomain = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0)) + +material = F.Material(D_0=1e-2, E_D=0) + +vol = F.VolumeSubdomain(id=5, material=material) +my_model.subdomains = [vol, left_subdomain, right_subdomain] +``` + +Now we'll add our species `tritium` to our model and define a custom flux to represent isotopic exchange. We use `ParticleFluxBC` and the custom flux function we defined above to define this BC on the right boundary: + +```{code-cell} ipython3 +tritium = F.Species("T") +my_model.species = [tritium] + +my_custom_flux = F.ParticleFluxBC( + value=my_custom_recombination_flux, + subdomain=right_subdomain, + species_dependent_value={"c": tritium}, + species=tritium, ) ``` -Finally, add the BC to your problem's `boundary_conditions` attribute using a list: +For diffusion across the mesh, we also define a fixed concentration on the left surface: ```{code-cell} ipython3 -my_model = F.HydrogenTransportProblem() -my_model.boundary_conditions = [my_bc] +my_model.boundary_conditions = [ + my_custom_flux, + F.FixedConcentrationBC(subdomain=left_subdomain, value=1e20, species=tritium), +] ``` -### Modeling recombination and dissociation +We'll export the flux using `SurfaceFlux` (see [exporting derived quantities to learn more](../post_process/derived.md)) and save this data to a variable `data_with_H2`: -Recombination and dissociation can also be modeled using `SurfaceReactionBC`, where the forward and backward rates of this reaction correspond to recombination and dissociation, respectively. +```{code-cell} ipython3 +surface_flux = F.SurfaceFlux(field=tritium, surface=right_subdomain) +``` -To model the reaction: +```{code-cell} ipython3 +my_model.temperature = 300 +my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, stepsize=50, final_time=500) +my_model.exports = [surface_flux] +my_model.initialise() +my_model.run() -$$ \mathrm{H} + \mathrm{H} \rightleftharpoons \mathrm{H_2}$$ +data_with_H2 = (surface_flux.data) +``` -where $ \text{Species A} = \text{Species B} = \text{H} $, assign your `reactants` list accordingly: +```{code-cell} ipython3 +:tags: [hide-input, hide-output] + +import pyvista +from dolfinx import plot + +topology, cell_types, geometry = plot.vtk_mesh(tritium.post_processing_solution.function_space) +grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) +grid.point_data["c"] = tritium.post_processing_solution.x.array +grid.set_active_scalars("c") +``` + +We can compare the surface flux to the case where we have no isotopic exchange by removing the custom boundary condition and saving the results to `data_no_H2`: ```{code-cell} ipython3 -H = F.Species("H") +my_model.boundary_conditions = [ + F.FixedConcentrationBC(subdomain=left_subdomain, value=1e20, species=tritium), +] + +my_model.exports = [surface_flux] +my_model.initialise() +my_model.run() +data_no_H2 = (surface_flux.data) +``` -my_recombination_bc = F.SurfaceReactionBC( - reactant=[A, A], - gas_pressure=1e5, - k_r0=1, - E_kr=0.1, - k_d0=1e-5, - E_kd=0.1, - subdomain=boundary, -) +```{code-cell} ipython3 +:tags: [hide-input, hide-output] + +topology, cell_types, geometry = plot.vtk_mesh(tritium.post_processing_solution.function_space) +grid2 = pyvista.UnstructuredGrid(topology, cell_types, geometry) +grid2.point_data["c"] = tritium.post_processing_solution.x.array +grid2.set_active_scalars("c") ``` -## Modeling isotopic exchange as a recombination flux +Let's plot both `data_with_H2` and `data_no_H2` to see the concentraton profile for each case: -Isotopic exchange reactions can be modeled as user-defined expressions using `ufl`. +```{code-cell} ipython3 +:tags: [hide-input] + +import matplotlib.pyplot as plt -Assuuming a large, constant concentration of molecular hydrogen $H_2$ at a boundary, we can define a recombination flux using rate and energy coefficients: +x = my_model.mesh.mesh.geometry.x[:, 0] +t = surface_flux.t +plt.plot(t, data_no_H2, label="No H_2") +plt.plot(t, data_with_H2, label="With H_2") + +plt.xlabel("Time (s)") +plt.ylabel("Surface Flux (right)") +plt.legend() +plt.show() +``` + +We see that the flux is higher with isotopic exchange, as we'd expect. Let's also take a look at the concentration profiles (top with isotopic exchange and bottom without): ```{code-cell} ipython3 -import ufl -import festim as F +:tags: [hide-input] -Kr_0 = 1.0 -E_Kr = 0.1 +pyvista.set_jupyter_backend("html") -def my_custom_recombination_flux(c, T): - Kr_0_custom = 1.0 - E_Kr_custom = 0.5 # eV - h2_conc = 1e25 # assumed constant H2 concentration in +plotter = pyvista.Plotter() - recombination_flux = ( - -(Kr_0 * ufl.exp(-E_Kr / (F.k_B * T))) * c**2 - - (Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))) * h2_conc * c - ) - return recombination_flux +plotter.add_mesh(grid) +plotter.view_xy() +if not pyvista.OFF_SCREEN: + plotter.show() +else: + figure = plotter.screenshot("concentration.png") ``` -```{note} -For more complex isotopic exchanges, we can also use `SurfaceReactionBC` to define recombination fluxes. See [NOTE] to learn more. +```{code-cell} ipython3 +:tags: [hide-input] + +pyvista.set_jupyter_backend("html") + +plotter = pyvista.Plotter() + +plotter.add_mesh(grid2) +plotter.view_xy() +if not pyvista.OFF_SCREEN: + plotter.show() +else: + figure = plotter.screenshot("concentration.png") ``` + +The results without isotopic exchange show virtually no diffusion for a given inlet concentration, indicating that isotopic exchange helps enchance diffusion! diff --git a/book/content/boundary_conditions/h_transport_basic.md b/book/content/boundary_conditions/h_transport_basic.md index c7a8703a..cf4cbc7c 100644 --- a/book/content/boundary_conditions/h_transport_basic.md +++ b/book/content/boundary_conditions/h_transport_basic.md @@ -12,7 +12,7 @@ kernelspec: name: python3 --- -# Hydrogen transport: basic boundary conditions +# Hydrogen transport: basic This section discusses how to implement basic boundary conditions (BCs) for hydrogen transport problems. Boundary conditions are essential to FESTIM simulations, as they describe the mathematical problem at the boundaries of the simulated domain. Read more about BCs _[here](https://festim.readthedocs.io/en/fenicsx/userguide/boundary_conditions.html)_. @@ -176,7 +176,6 @@ grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) grid.point_data["c"] = H.post_processing_solution.x.array grid.set_active_scalars("c") -pyvista.start_xvfb() pyvista.set_jupyter_backend("html") plotter = pyvista.Plotter() From 494228f9348ea3fe808499a0e46abe4141b69671 Mon Sep 17 00:00:00 2001 From: Chirag Date: Fri, 13 Feb 2026 17:04:47 -0500 Subject: [PATCH 11/17] fixed code cells and removed old files --- .../boundary_conditions/concentration.md | 300 ---------------- book/content/boundary_conditions/fluxes.md | 195 ---------- .../h_transport_advanced.md | 2 +- .../boundary_conditions/h_transport_basic.md | 3 +- .../boundary_conditions/surface_reactions.md | 337 ------------------ 5 files changed, 2 insertions(+), 835 deletions(-) delete mode 100644 book/content/boundary_conditions/concentration.md delete mode 100644 book/content/boundary_conditions/fluxes.md delete mode 100644 book/content/boundary_conditions/surface_reactions.md diff --git a/book/content/boundary_conditions/concentration.md b/book/content/boundary_conditions/concentration.md deleted file mode 100644 index fa54ccd3..00000000 --- a/book/content/boundary_conditions/concentration.md +++ /dev/null @@ -1,300 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.17.3 -kernelspec: - display_name: festim-workshop - language: python - name: python3 ---- - -# Concentration # - -Boundary conditions (BCs) are essential to FESTIM simulations, as they describe the mathematical problem at the boundaries of the simulated domain. Read more about BCs _[here](https://festim.readthedocs.io/en/fenicsx/userguide/boundary_conditions.html)_. - -This tutorial goes over how to define concentration boundary conditions for hydrogen transport simulations. - -Objectives: -* Understand the mathematics behind a fixed concentration BC -* Define a fixed concentration BC -* Choose which solubility law (Sieverts' or Henry's) -* Solve a hydrogen transport problem with plasma implantation - -+++ - -## Understanding mathematics behind a fixed concentration BC ## - -A fixed concentration (Dirichlet) boundary condition prescribes the value of the mobile hydrogen isotope concentration at a boundary. This enforces the concentration to remain constant in time and space on the specified boundary, independent of the solution in the bulk. - -This boundary condition is typically used to represent surfaces in equilibrium with an infinite reservoir, imposed implantation conditions, or experimentally controlled concentrations. - -### Mathematical formulation - -On a boundary $\Gamma_D$, the mobile concentration satisfies - -$$ -c(\mathbf{x}, t) = c_0 \quad \text{for } \mathbf{x} \in \Gamma_D, -$$ - -where $c_0$ is the prescribed concentration value. - -In the weak formulation, this condition is enforced by directly constraining the degrees of freedom associated with the boundary. - -+++ - -## Defining fixed concentration ## - -Users can prescribe a fixed concentration on a boundary using `FixedConcentrationBC`, which can depend on temperature, time, and space. - -Let us consider a 2D, steady-state, single-species example with the following *space-dependent* boundary conditions: - -$$ \text{Top:} \quad c = 2x + y $$ -$$ \text{Right:} \quad c = y^2 + 1 $$ - -```{code-cell} ipython3 -import festim as F -import numpy as np -from dolfinx.mesh import create_unit_square -from mpi4py import MPI - -top_subdomain = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1)) -right_subdomain = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1)) -left_subdomain = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0)) -bottom_subdomain = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[1], 0)) - -my_model = F.HydrogenTransportProblem() - -my_model.mesh = F.Mesh(create_unit_square(MPI.COMM_WORLD, 10, 10)) - -H = F.Species("H") -my_model.species = [H] -mat = F.Material(D_0=1e-3, E_D=0) -my_model.material = mat -vol = F.VolumeSubdomain(id=5, material=mat) -my_model.subdomains = [top_subdomain, right_subdomain, left_subdomain, bottom_subdomain, vol] - -my_model.temperature = 400 -my_model.settings = F.Settings(atol=1e-8, rtol=1e-8, transient=False) -``` - -We can define functions for the top and right boundary conditions using `lambda` functions: - -```{code-cell} ipython3 -top_bc_function = lambda x: 2.0 * x[0] + x[1] -right_bc_function = lambda x: x[1] ** 2 + 1.0 -``` - - - -+++ - -To include these boundary conditions to our problem, we use `FixedConcentrationBC`. We must also specify which subdomain (boundary) each BC is applied to, as well as the corresponding species: - -```{code-cell} ipython3 -top_bc = F.FixedConcentrationBC( - subdomain=top_subdomain, - value=top_bc_function, - species=H -) - -right_bc = F.FixedConcentrationBC( - subdomain=right_subdomain, - value=right_bc_function, - species=H -) - -# left_bc = F.FixedConcentrationBC( -# subdomain=left_subdomain, -# value=.0, -# species=H -# ) - -# bottom_bc = F.FixedConcentrationBC( -# subdomain=bottom_subdomain, -# value=0.0, -# species=H -# ) -``` - -Finally, we add our BCs to `my_model.boundary_conditions` using a list, and then run the model: - -```{code-cell} ipython3 -my_model.boundary_conditions = [top_bc, right_bc] -my_model.initialise() -my_model.run() -``` - -```{code-cell} ipython3 -:tags: [hide-input] - -import pyvista -from dolfinx import plot - -topology, cell_types, geometry = plot.vtk_mesh(H.post_processing_solution.function_space) -grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) -grid.point_data["c"] = H.post_processing_solution.x.array -grid.set_active_scalars("c") - -pyvista.start_xvfb() -pyvista.set_jupyter_backend("html") - -plotter = pyvista.Plotter() - -plotter.add_mesh(grid) -plotter.view_xy() -contours = grid.contour(isosurfaces=5, scalars="c") -plotter.add_mesh(contours, color="r", line_width=0.1) -if not pyvista.OFF_SCREEN: - plotter.show() -else: - figure = plotter.screenshot("concentration.png") -``` - -```{note} - -In this example, we did not define any boundary conditions for the left and bottom surfaces. If no BC is set on a boundary, it is assumed that the flux is null (symmetry boundary condition): - -$$ -\mathbf{J} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma -$$ - -``` - -+++ - -### Time and temperature dependent boundary conditions ### - -Users can also impose concentration BCs that are dependent on space, time and temperature, such as: - -$$ c = 10 + x^2 + T(t+2) $$ - -To do so, we define a Python function: - -```{code-cell} ipython3 -my_custom_value = lambda x, t, T: 10 + x[0]**2 + T *(t + 2) - -boundary = F.SurfaceSubdomain(id=1) -my_bc = F.FixedConcentrationBC(subdomain=boundary, value=my_custom_value, species=H) -``` - -## Choosing a solubility law ## - -Users can define the surface concentration using either Sieverts’ law, $c = S(T)\sqrt P$, or Henry's law, $c=K_H P$, where $S(T)$ and $K_H$ denote the temperature-dependent Sieverts’ and Henry’s solubility coefficients, respectively, and $P$ is the partial pressure of the species on the surface. - -For Sieverts' law of solubility, we can use `festim.SievertsBC`: - -```{code-cell} ipython3 -from festim import SievertsBC, SurfaceSubdomain, Species - -boundary = SurfaceSubdomain(id=1) -H = Species(name="Hydrogen") - -custom_pressure_value = lambda t: 2 + t -my_bc = SievertsBC(subdomain=3, S_0=2, E_S=0.1, species=H, pressure=custom_pressure_value) -``` - -Similarly, for Henry's law of solubility, we can use `festim.HenrysBC`: - -```{code-cell} ipython3 -from festim import HenrysBC - -pressure_value = lambda t: 5 * t -my_bc = HenrysBC(subdomain=3, H_0=1.5, E_H=0.2, species=H, pressure=pressure_value) -``` - -## Plasma implantation approximation ## - -+++ - -We can also approximate plasma implantation using FESTIM's `ParticleSource` class, which is helpful in modeling thermo-desorption spectra (TDS) experiments or including the effect of plasma exposure on hydrogen transport. Learn more about how FESTIM approximates plasma implantation _[here](https://festim.readthedocs.io/en/fenicsx/theory.html)_. - -Consider the following 1D plasma implantation problem, where we represent the plasma as a hydrogen source $S_{ext}$: - -$$ S_{ext} = \varphi \cdot f(x) $$ - -$$\varphi = 1\cdot 10^{13} \quad \mathrm{m}^{-2}\mathrm{s}^{-1}$$ - -where $\varphi$ is the implantation flux and $f(x)$ is a Gaussian spatial distribution (distribution mean value of 0.5 $\text{m}$ and width of 1 $\text{m}$). - -First, we setup a 1D mesh ranging from $ [0,1] $ and assign the subdomains and material: - -```{code-cell} ipython3 -import festim as F -import ufl -import numpy as np - -my_model = F.HydrogenTransportProblem() -vertices = np.linspace(0,1,2000) -my_model.mesh = F.Mesh1D(vertices) - -mat = F.Material(D_0=0.1, E_D=0.01) - -volume_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat) -left_boundary = F.SurfaceSubdomain1D(id=1, x=0) -right_boundary = F.SurfaceSubdomain1D(id=2, x=1) - -my_model.subdomains = [ - volume_subdomain, - left_boundary, - right_boundary, -] -``` - -Now, we define our `incident_flux` and `gaussian_distribution` function. We can use `ParticleSource` to represent the source term, and then add it to our model: - -```{code-cell} ipython3 -incident_flux = 1e13 - -def gaussian_distribution(x, center, width): - return ( - 1 - / (width * (2 * ufl.pi) ** 0.5) - * ufl.exp(-0.5 * ((x[0] - center) / width) ** 2) - ) -H = F.Species("H") -my_model.species = [H] - -source_term = F.ParticleSource( - value=lambda x: incident_flux * gaussian_distribution(x, .5, 1), - volume=volume_subdomain, - species=H, -) - -my_model.sources = [source_term] -``` - -Finally, we assign boundary conditions (zero concentration at $x=0$ and $x=1$) and solve our problem: - -```{code-cell} ipython3 -my_model.boundary_conditions = [ - F.FixedConcentrationBC(subdomain=left_boundary, value=0, species=H), - F.FixedConcentrationBC(subdomain=right_boundary, value=0, species=H), -] - -my_model.temperature = 400 -my_model.settings = F.Settings(atol=1e10, rtol=1e-10, transient=False) - -profile_export = F.Profile1DExport(field=H,subdomain=volume_subdomain) -my_model.exports = [profile_export] - -my_model.initialise() -my_model.run() -``` - -```{code-cell} ipython3 -:tags: [hide-input] - -import matplotlib.pyplot as plt - -x = my_model.exports[0].x -c = my_model.exports[0].data[0][0] - -plt.plot(x, c) -plt.show() -``` diff --git a/book/content/boundary_conditions/fluxes.md b/book/content/boundary_conditions/fluxes.md deleted file mode 100644 index b5a6d9c9..00000000 --- a/book/content/boundary_conditions/fluxes.md +++ /dev/null @@ -1,195 +0,0 @@ ---- -jupytext: - formats: md:myst,ipynb - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.17.3 -kernelspec: - display_name: festim-workshop - language: python - name: python3 ---- - -# Particle flux # - -This tutorial goes over how to define particle flux boundary conditions in FESTIM. - -Objectives: -* Learn how to define fluxes at boundaries -* Impose boundary fluxes as a function of the concentration -* Solve a problem with multi-species flux boundary conditions - -+++ - -## Defining flux at boundaries ## - -Users can impose the particle flux (Neumann BC) at boundaries using `ParticleFluxBC` class: - -```{code-cell} ipython3 -from festim import ParticleFluxBC, Species, SurfaceSubdomain - -boundary = SurfaceSubdomain(id=1) -H = Species(name="Hydrogen") - -my_flux_bc = ParticleFluxBC(subdomain=boundary, value=2, species=H) -``` - -## Defining concentration-dependant fluxes ## - -Similar to the concentration, the flux can be dependent on space, time and temperature. But for particle fluxes, the values can also be dependent on a species’ concentration. - -For example, let's define a hydrogen flux `J` that depends on the hydrogen concentration `c` and time `t`: - -$$ J(c,t) = 10t^2 + 2c $$ - -```{code-cell} ipython3 -from festim import ParticleFluxBC, Species, SurfaceSubdomain - -boundary = SurfaceSubdomain(id=1) -H = Species(name="Hydrogen") - -J = lambda t, c: 10*t**2 + 2*c - -my_flux_bc = ParticleFluxBC( - subdomain=boundary, - value=J, - species=H, - species_dependent_value={"c": H}, -) -``` - -## Multi-species flux boundary conditions ## - -+++ - -Users may also need to impose a flux boundary condition which depends on the concentrations of multiple species. - -Consider the following 1D example with three species, $\mathrm{A}$, $\mathrm{B}$, and $\mathrm{C}$, with recombination fluxes $\phi_{AB}$ and $\phi_{ABC}$ that depend on the concentrations $\mathrm{c}$: - -$$ \phi_{AB} = -c_A c_B $$ - -$$ \phi_{ABC}= -10 c_A c_B c_C$$ - -We must first define each species using `Species` and then create the dictionary to be passed into `species_dependent_value`. The dictionary maps each argument in the custom flux function to the corresponding defined species. We also define our custom functions for the fluxes: - -```{code-cell} ipython3 -import festim as F - -my_model = F.HydrogenTransportProblem() - -A = F.Species(name="A") -B = F.Species(name="B") -C = F.Species(name="C") - -my_model.species = [A, B, C] -species_dependent_value = {"c_A": A, "c_B": B, "c_C": C} - -recombination_flux_AB = lambda c_A, c_B, c_C: -c_A*c_B -recombination_flux_ABC = lambda c_A, c_B, c_C: -10*c_A*c_B*c_C -``` - -Now, we create our 1D mesh and assign boundary conditions (recombination on the right, fixed concentration on the left). - -The boundary condition `ParticleFluxBC` must be added for each species: - -```{code-cell} ipython3 -import numpy as np - -my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100)) - -D = 1e-2 -mat = F.Material( - D_0={A: 8*D, B: 7*D, C: D}, - E_D={A: 0.01, B: 0.01, C: 0.01}, -) - -bulk = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat) -left = F.SurfaceSubdomain1D(id=1, x=0) -right = F.SurfaceSubdomain1D(id=2, x=1) -my_model.subdomains = [bulk, left, right] - -my_model.boundary_conditions = [ - F.ParticleFluxBC( - subdomain=right, - value=recombination_flux_AB, - species=A, - species_dependent_value=species_dependent_value, - ), - F.ParticleFluxBC( - subdomain=right, - value=recombination_flux_AB, - species=B, - species_dependent_value=species_dependent_value, - ), - F.ParticleFluxBC( - subdomain=right, - value=recombination_flux_ABC, - species=A, - species_dependent_value=species_dependent_value, - ), - F.ParticleFluxBC( - subdomain=right, - value=recombination_flux_ABC, - species=B, - species_dependent_value=species_dependent_value, - ), - F.ParticleFluxBC( - subdomain=right, - value=recombination_flux_ABC, - species=C, - species_dependent_value=species_dependent_value, - ), - F.FixedConcentrationBC(subdomain=left,value=1,species=A), - F.FixedConcentrationBC(subdomain=left,value=1,species=B), - F.FixedConcentrationBC(subdomain=left,value=1,species=C), -] - -right_flux_A = F.SurfaceFlux(field=A,surface=right) -right_flux_B = F.SurfaceFlux(field=B,surface=right) -right_flux_C = F.SurfaceFlux(field=C,surface=right) - -my_model.exports = [ - right_flux_A, - right_flux_B, - right_flux_C, -] -``` - -```{note} -The diffusivity pre-factor `D_0` and activation energy `E_D` must be defined for each species in `Material`. Learn more about defining multi-species material properties __[here](https://festim-workshop.readthedocs.io/en/festim2/content/material/material_advanced.html#defining-species-dependent-material-properties)__. -``` - -+++ - -Finally, let's solve and plot the profile for each species: - -```{code-cell} ipython3 -my_model.temperature = 300 -my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) - -my_model.initialise() -my_model.run() -``` - -```{code-cell} ipython3 -:tags: [hide-input] - -import matplotlib.pyplot as plt - -def plot_profile(species, **kwargs): - c = species.post_processing_solution.x.array[:] - x = species.post_processing_solution.function_space.mesh.geometry.x[:,0] - return plt.plot(x, c, **kwargs) - -for species in my_model.species: - plot_profile(species, label=species.name) - -plt.xlabel('Position') -plt.ylabel('Concentration') -plt.legend() -plt.show() -``` - -We see the higher recombination flux for $\mathrm{ABC}$ decreases the concentration of $\mathrm{C}$ throughout the material. diff --git a/book/content/boundary_conditions/h_transport_advanced.md b/book/content/boundary_conditions/h_transport_advanced.md index 84f8ffd7..621954a3 100644 --- a/book/content/boundary_conditions/h_transport_advanced.md +++ b/book/content/boundary_conditions/h_transport_advanced.md @@ -314,7 +314,7 @@ K_{r0} \exp\!\left(-\frac{E_{Kr}}{k_B T}\right) c_T^2 \Bigg] v \, d\Gamma $$ -These fluxes can be implemented in FESTIM using `ParticleFluxBC` with user-defined expressions for each reaction term (see [here](examples.ipynb) to learn more). +These fluxes can be implemented in FESTIM using `ParticleFluxBC` with user-defined expressions for each reaction term, as shown below. +++ diff --git a/book/content/boundary_conditions/h_transport_basic.md b/book/content/boundary_conditions/h_transport_basic.md index cf4cbc7c..f477a75d 100644 --- a/book/content/boundary_conditions/h_transport_basic.md +++ b/book/content/boundary_conditions/h_transport_basic.md @@ -226,8 +226,7 @@ This custom function is equivalent to: ` my_custom_value = lambda x, t, T: 10 + x[0]**2 + T *(t + 2) ` - -+++ +``` ### Multi-species fixed concentration BC ## diff --git a/book/content/boundary_conditions/surface_reactions.md b/book/content/boundary_conditions/surface_reactions.md deleted file mode 100644 index 87645651..00000000 --- a/book/content/boundary_conditions/surface_reactions.md +++ /dev/null @@ -1,337 +0,0 @@ ---- -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 ---- - -# Surface reactions # - -+++ - -FESTIM V2.0 allows users to impose surface reactions on boundaries. Learn more about surface reactions __[here](https://festim-workshop.readthedocs.io/en/festim2/content/species_reactions/surface_reactions.html)__. - -Objectives: -* Defining a simple surface reaction boundary condition -* Recombination and dissociation -* Isotopic exhange -* Complex isotopic exchange with multple hydrogenic species - -+++ - -## Defining a simple surface reaction boundary condition ## - -Surface reactions between adsorbed species and gas-phase products can be represented using the `SurfaceReactionBC` class. This boundary condition defines a reversible reaction of the form: - -$$ -A + B \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad C -$$ - -where $\mathrm{A}$ and $\mathrm{B}$ are surface reactants and $\mathrm{C}$ is the product species in the gas phase. - -### Reaction Rate Formulation ### - -The forward and backward rate constants follow Arrhenius laws: - -$$ -K_r = k_{r0} e^{-E_{kr} / (k_B T)}, \qquad -K_d = k_{d0} e^{-E_{kd} / (k_B T)} -$$ - -The net surface reaction rate is given by: - -$$ -K = K_r c_A c_B - K_d P_C -$$ - -where: -- $\mathrm{k_{r0}}, \mathrm{k_{d0}}$ are rate pre-exponential constants -- $\mathrm{c_A}, \mathrm{c_B}$ are concentrations of reactant species $\mathrm{A}$ and $\mathrm{B}$ at the surface -- $\mathrm{P_C}$ is the partial pressure of the product species $\mathrm{C}$ -- $\mathrm{k_B}$ is the Boltzmann constant -- $\mathrm{T}$ is the surface temperature - -The flux of species $\mathrm{A}$ entering the surface is equal to $\mathrm{K}$ (if $\mathrm{A=B}$, the total particle flux entering the surface becomes $\mathrm{2K}$). - -We can use the `SurfaceReactionBC` class to impose the surface reaction above by specifying the reactants (`reactant`), gas pressure (`gas_pressure`), forward/backward rate constants (`k_r0` and `k_d0`), and rate activation energies (`E_kr` and `E_kd`): - -```{code-cell} ipython3 -from festim import Species, SurfaceReactionBC, SurfaceSubdomain - -boundary = SurfaceSubdomain(id=1) -A = Species("A") -B = Species("B") -C = Species("C") - -my_bc = SurfaceReactionBC( - reactant=[A, B], - gas_pressure=1e5, - k_r0=1, - E_kr=0.1, - k_d0=0, - E_kd=0, - subdomain=boundary, -) -``` - -## Recombination and dissociation ## - -Hydrogen recombination/dissociation can also be modelled using `SurfaceReactionBC`, such as this reaction: - -$$ \mathrm{H + H} \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad \mathrm{H_2} $$ - -```{code-cell} ipython3 -from festim import Species, SurfaceReactionBC, SurfaceSubdomain - -boundary = SurfaceSubdomain(id=1) -H = Species("H") - -my_bc = SurfaceReactionBC( - reactant=[H, H], - gas_pressure=1e5, - k_r0=1, - E_kr=0.1, - k_d0=1e-5, - E_kd=0.1, - subdomain=boundary, -) -``` - -## Isotopic exchange ## - -Isotopic exchange can occur where isotopes can swap positions with other isotopes, such as: - -$$\mathrm{T + H_2} \rightleftharpoons \mathrm{H} + \mathrm{HT}$$ - -These exchange processes occur through complex mechanisms and are important in understanding the behavior of hydrogenic species: - -If the $\mathrm{H_2}$ concentration is assumed much larger than $\mathrm{HT}$, the flux $\phi_T$ reduces to a first-order process in $\mathrm{T}$: - -$$ -\phi_T = --\mathrm{K_{r0}} \exp\!\left(-\frac{\mathrm{E_{Kr}}}{\mathrm{k_B T}}\right) c_T^2 --\mathrm{K_{r0}^\ast} \exp\!\left(-\frac{\mathrm{E_{Kr}^\ast}}{\mathrm{k_B T}}\right) c_{H_2} c_T -$$ - -Such fluxes can be implemented using `festim.ParticleFluxBC` with user-defined expressions. - -```{code-cell} ipython3 -import ufl -import festim as F - -Kr_0 = 1.0 -E_Kr = 0.1 - -def my_custom_recombination_flux(c, T): - Kr_0_custom = 1.0 - E_Kr_custom = 0.5 # eV - h2_conc = 1e25 # assumed constant H2 concentration in - - recombination_flux = ( - -(Kr_0 * ufl.exp(-E_Kr / (F.k_B * T))) * c**2 - - (Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))) * h2_conc * c - ) - return recombination_flux -``` - -```{note} -We can also use `F.SurfaceReactionBC` to define recombination fluxes, which we do in the next section for more complex isotopic exchanges. -``` - -+++ - -Let's work through a 1D example to illustrate the effect of isotopic exchange. We also define a fixed concentration on the left surface: - -```{code-cell} ipython3 -import numpy as np - -my_model = F.HydrogenTransportProblem() -my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 100)) - -left_surf = F.SurfaceSubdomain1D(id=1, x=0) -right_surf = F.SurfaceSubdomain1D(id=2, x=1) - -material = F.Material(D_0=1e-2, E_D=0) - -vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=material) - -my_model.subdomains = [vol, left_surf, right_surf] -tritium = F.Species("T") -my_model.species = [tritium] - -my_custom_flux = F.ParticleFluxBC( - value=my_custom_recombination_flux, - subdomain=right_surf, - species_dependent_value={"c": tritium}, - species=tritium, -) - -my_model.boundary_conditions = [ - my_custom_flux, - F.FixedConcentrationBC(subdomain=left_surf, value=1, species=tritium), -] - -my_model.temperature = 300 -right_flux = F.SurfaceFlux(field=tritium, surface=right_surf) - -my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) - -my_model.initialise() -my_model.run() - -data_with_H2 = (tritium.solution.x.array) -``` - -We can compare the surface flux to the case where we have no isotopic exchange by removing the custom boundary condition: - -```{code-cell} ipython3 -my_model.boundary_conditions = [ - F.FixedConcentrationBC(subdomain=left_surf, value=1, species=tritium), -] - -my_model.initialise() -my_model.run() -data_no_H2 = (tritium.solution.x.array) -``` - -```{code-cell} ipython3 -:tags: [hide-input] - -import matplotlib.pyplot as plt - -x = my_model.mesh.mesh.geometry.x[:, 0] -plt.plot(x, data_no_H2, label="No H_2") -plt.plot(x, data_with_H2, label="With H_2") - -plt.xlabel("x (m)") -plt.ylabel("Surface Flux (right)") -plt.legend() -plt.show() -``` - -We see that diffusion is present with isotopic exchange! - -+++ - -## Complex isotopic exchange with multple hydrogenic species ## - -Surface reactions can involve multiple hydrogen isotopes, allowing for the modeling of complex isotope-exchange mechanisms between species. For example, in a system with both mobile hydrogen and deuteriun, various molecular recombination pathways may occur at the surface, resulting in the formation of $H_2$, $D_2$, and $HD$: - -$$ \text{Reaction 1}: \mathrm{H+D} \rightarrow \mathrm{HD} \longrightarrow \phi_1 = K_{r1} c_H c_D - K_{d1}P_{HD} $$ -$$ \text{Reaction 2}: \mathrm{D+D} \rightarrow \mathrm{D_2} \longrightarrow \phi_2 = 2K_{r2} c_D^2 - K_{d2}P_{D2} $$ -$$ \text{Reaction 3}: \mathrm{D+H_2} \rightarrow \mathrm{HD + H} \longrightarrow \phi_3 = K_{r3} c_H c_D - K_{d3}P_{HD} $$ -$$ \text{Reaction 4}: \mathrm{H+H} \rightarrow \mathrm{H_2} \longrightarrow \phi_4 = 2K_{r4} c_H^2 - K_{d4}P_{H2} $$ - -Now consider the case where deuterium diffuses from left to right and reacts with background -$\mathrm{H_2}$, while $\mathrm{P_{HD}}$ and $\mathrm{P_{D_2}}$ are negligible at the surface. -Formation of $\mathrm{H}$ at the right boundary induces back-diffusion toward the left, -even though none existed initially. - -The boundary conditions for this scenario are: - -$$ -c_D(x=0) = 1, \qquad c_H(x=0) = 0, \qquad P_{H2}(x=1) = \text{1000 Pa} -$$ - -First, let's define a 1D mesh ranging from $\mathrm{x=[0,1]}$: - -```{code-cell} ipython3 -import numpy as np -import festim as F - -my_model = F.HydrogenTransportProblem() -my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 100)) - -left_surf = F.SurfaceSubdomain1D(id=1, x=0) -right_surf = F.SurfaceSubdomain1D(id=2, x=1) - -material = F.Material(D_0=1e-2, E_D=0) -vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=material) - -my_model.subdomains = [vol, left_surf, right_surf] -``` - -Now, we define our species at recombination reactions using `SurfaceReactionBC`: - -```{code-cell} ipython3 -H = F.Species("H") -D = F.Species("D") -my_model.species = [H, D] - -H2 = F.SurfaceReactionBC( - reactant=[H, H], - gas_pressure=100000, - k_r0=1, - E_kr=0.1, - k_d0=1e-5, - E_kd=0.1, - subdomain=right_surf, -) - -HD = F.SurfaceReactionBC( - reactant=[H, D], - gas_pressure=0, - k_r0=1, - E_kr=0.1, - k_d0=1e-5, - E_kd=0.1, - subdomain=right_surf, -) - -D2 = F.SurfaceReactionBC( - reactant=[D, D], - gas_pressure=0, - k_r0=1, - E_kr=0.1, - k_d0=1e-5, - E_kd=0.1, - subdomain=right_surf, -) -``` - -Finally, we add our boundary conditions and solve the steady-state problem: - -```{code-cell} ipython3 -my_model.boundary_conditions = [ - H2, - D2, - HD, - F.FixedConcentrationBC(subdomain=left_surf, value=1, species=D), - F.FixedConcentrationBC(subdomain=left_surf, value=0, species=H), -] - -my_model.temperature = 300 - -my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) - -my_model.initialise() -my_model.run() -``` - -```{code-cell} ipython3 -:tags: [hide-input] - -import matplotlib.pyplot as plt - -def plot_profile(species, **kwargs): - c = species.post_processing_solution.x.array[:] - x = species.post_processing_solution.function_space.mesh.geometry.x[:,0] - return plt.plot(x, c, **kwargs) - -for species in my_model.species: - plot_profile(species, label=species.name) - -plt.xlabel('Position') -plt.ylabel('Concentration') -plt.legend() -plt.show() -``` - -We see that the background $\mathrm{H_2}$ reacts with the $\mathrm{D}$, removing the total amount of $\mathrm{D}$ from the surface. Conversely, the $\mathrm{H}$ diffuses from the surface towards the left. From 62521126a52b10ac530ff1e9d8eca2b07924240e Mon Sep 17 00:00:00 2001 From: Chirag Date: Fri, 20 Feb 2026 10:47:07 -0500 Subject: [PATCH 12/17] refined mesh near implantation for example --- book/content/boundary_conditions/examples.md | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/book/content/boundary_conditions/examples.md b/book/content/boundary_conditions/examples.md index bde78add..09663363 100644 --- a/book/content/boundary_conditions/examples.md +++ b/book/content/boundary_conditions/examples.md @@ -43,8 +43,22 @@ import numpy as np L = 20e-6 my_model = F.HydrogenTransportProblem() -vertices = np.linspace(0, L, 1000) + +vertices = np.concatenate( + [ + np.linspace(0, 30e-9, num=200), + np.linspace(3e-6, L, num=500), + ] +) + my_model.mesh = F.Mesh1D(vertices) +``` + +```{note} +We break up the mesh region into two regions so we can refine the region closer to the implantation depth (defined below) +``` + +```{code-cell} ipython3 tungsten = F.Material( D_0=4.1e-07, # m2/s From a5e171960b67b4b611d58b848829b52df8e94cb8 Mon Sep 17 00:00:00 2001 From: ck768 Date: Mon, 23 Feb 2026 13:45:21 -0500 Subject: [PATCH 13/17] addressed changes from review. moved basic BC theory to intro section, removed weak formulation discussions) --- book/_toc.yml | 1 + book/content/boundary_conditions/examples.md | 31 ++++- .../h_transport_advanced.md | 120 ++++-------------- .../boundary_conditions/h_transport_basic.md | 99 +++------------ 4 files changed, 74 insertions(+), 177 deletions(-) diff --git a/book/_toc.yml b/book/_toc.yml index f77c850b..29efa335 100644 --- a/book/_toc.yml +++ b/book/_toc.yml @@ -23,6 +23,7 @@ parts: - caption: Boundary Conditions chapters: + - file: content/boundary_conditions/bc_intro - file: content/boundary_conditions/h_transport_basic - file: content/boundary_conditions/h_transport_advanced - file: content/boundary_conditions/examples diff --git a/book/content/boundary_conditions/examples.md b/book/content/boundary_conditions/examples.md index 09663363..2d782e40 100644 --- a/book/content/boundary_conditions/examples.md +++ b/book/content/boundary_conditions/examples.md @@ -59,7 +59,6 @@ We break up the mesh region into two regions so we can refine the region closer ``` ```{code-cell} ipython3 - tungsten = F.Material( D_0=4.1e-07, # m2/s E_D=0.39, # eV @@ -245,7 +244,7 @@ my_model.species = [H, D] H2 = F.SurfaceReactionBC( reactant=[H, H], - gas_pressure=100000, + gas_pressure=1000, k_r0=1, E_kr=0.1, k_d0=1e-5, @@ -274,6 +273,33 @@ D2 = F.SurfaceReactionBC( ) ``` +Now, let's add our isotopic exchange reaction using `ParticleFluxBC` (see [here](h_transport_advanced.md) to learn more about defining isotopic exchange fluxes): + +```{code-cell} ipython3 +import ufl + +Kr_0 = 1.0 +E_Kr = 0.1 + +def my_custom_recombination_flux(c, T): + Kr_0_custom = 1.0 + E_Kr_custom = 0.5 # eV + h2_conc = 1e25 # assumed constant H2 concentration in + + recombination_flux = ( + -(Kr_0 * ufl.exp(-E_Kr / (F.k_B * T))) * c**2 + - (Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))) * h2_conc * c + ) + return recombination_flux + +HDH = F.ParticleFluxBC( + value=my_custom_recombination_flux, + subdomain=right_surf, + species_dependent_value={"c": D}, + species=D, +) +``` + Finally, we add our boundary conditions and solve the steady-state problem: ```{code-cell} ipython3 @@ -281,6 +307,7 @@ my_model.boundary_conditions = [ H2, D2, HD, + HDH, F.FixedConcentrationBC(subdomain=left_surf, value=1, species=D), F.FixedConcentrationBC(subdomain=left_surf, value=0, species=H), ] diff --git a/book/content/boundary_conditions/h_transport_advanced.md b/book/content/boundary_conditions/h_transport_advanced.md index 621954a3..8c9b8dee 100644 --- a/book/content/boundary_conditions/h_transport_advanced.md +++ b/book/content/boundary_conditions/h_transport_advanced.md @@ -63,16 +63,6 @@ where: Users interested in liquid metal interactions with hydrogen will usually need to use **Sievert's Law**. ``` -#### Weak formulation - -In the weak form, solubility is imposed as a **Robin-type boundary condition**: - -$$ -\int_{\Gamma_s} k_s \left(c - c_s \right) v \, \mathrm{d}\Gamma -$$ - -where $k_s$ is a surface reaction rate constant controlling the exchange between bulk and surface concentration. This allows modeling of adsorption/desorption and surface-limited transport at material boundaries. - +++ ### Imposing a solubility law @@ -109,7 +99,6 @@ Using `HenrysBC` or `SievertsBC` is just a convenient way for users to define a +++ - --- +++ @@ -162,16 +151,25 @@ $$ where $\mathbf{n}$ is the outward normal vector at the surface. -#### Weak form contribution +### Recombination and dissociation -In the weak formulation, the surface reaction contributes as a Robin-type term: +If $\text{Species A = Species B}$, recombination or dissociation can also be modeled with `SurfaceReactionBC`, e.g.: $$ -\int_{\Gamma_s} \mathbf{J}_A \cdot \mathbf{n} \, v \, d\Gamma = -\int_{\Gamma_s} K \, v \, d\Gamma +\mathrm{H + H} \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad \mathrm{H_2} $$ -where $v$ is a test function. +where + +$$ +K = K_r \, c_H^2 - K_d \, P_{H_2} +$$ + +and corresponding flux of atomic hydrogen is: + +$$ +\mathbf{J}_H \cdot \mathbf{n} = - 2 K +$$ +++ @@ -189,7 +187,7 @@ A = F.Species("A") B = F.Species("B") ``` -Now, we add these as arguments to the `SurfaceReactionBC` class. We'll also need to assign a `gas_pressure` (corresponding to the partial pressure of the product species), and the forward/backward rate (`k_r0`, `k_d0`) and energy (`E_kr`, `E_dr`) coefficients (see above to learn more about these parameters): +Now, we add these as arguments to the `SurfaceReactionBC` class. We'll also need to assign a `gas_pressure` (corresponding to the partial pressure of the product species), and the forward/backward rate (`k_r0`, `k_d0`) and energy (`E_kr`, `E_kd`) coefficients (see above to learn more about these parameters): ```{code-cell} ipython3 my_bc = F.SurfaceReactionBC( @@ -210,44 +208,9 @@ my_model = F.HydrogenTransportProblem() my_model.boundary_conditions = [my_bc] ``` ---- - -+++ - -## Recombination and dissociation - -Hydrogen recombination or dissociation can be modeled with `SurfaceReactionBC`, e.g.: - -$$ -\mathrm{H + H} \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad \mathrm{H_2} -$$ - -### Mathematical formulation - -Let: -- $c_H$ be the surface concentration of atomic hydrogen -- $P_{H_2}$ the partial pressure of molecular hydrogen - -The net reaction rate is: - -$$ -K = K_r \, c_H^2 - K_d \, P_{H_2} -$$ - -The corresponding flux of atomic hydrogen is: - -$$ -\mathbf{J}_H \cdot \mathbf{n} = - 2 K -$$ - -#### Weak form contribution - -The weak form contribution of recombination flux is: - -$$ -\int_{\Gamma_s} \mathbf{J}_H \cdot \mathbf{n} \, v \, d\Gamma = -- \int_{\Gamma_s} 2 K \, v \, d\Gamma -$$ +```{note} +Using `SurfaceReactionBC` is just a convenient way for users to define surface fluxes, which can be done manually (as shown below in the isotopic exchange example). +``` +++ @@ -255,35 +218,13 @@ $$ +++ -### Modeling recombination and dissociation - -Recombination and dissociation can also be modeled using `SurfaceReactionBC`, where the forward and backward rates of this reaction correspond to recombination and dissociation, respectively. - -To model the reaction: - -$$ \mathrm{H} + \mathrm{H} \rightleftharpoons \mathrm{H_2}$$ - -where $ \text{Species A} = \text{Species B} = \text{H} $, assign your `reactants` list accordingly: - -```{code-cell} ipython3 -H = F.Species("H") - -my_recombination_bc = F.SurfaceReactionBC( - reactant=[A, A], - gas_pressure=1e5, - k_r0=1, - E_kr=0.1, - k_d0=1e-5, - E_kd=0.1, - subdomain=boundary, -) -``` - ## Isotopic exchange Isotopic exchange occurs when hydrogenic isotopes swap positions, for example: $$ +\mathrm{T + T} \rightleftharpoons \mathrm{T_2} +\\ \mathrm{T + H_2} \rightleftharpoons \mathrm{H} + \mathrm{HT} $$ @@ -302,18 +243,6 @@ $$ - K_{r0}^\ast \exp\!\left(-\frac{E_{Kr}^\ast}{k_B T}\right) c_{H_2} c_T $$ -#### Weak form contribution - -The weak form contribution for tritium flux is: - -$$ -\int_{\Gamma_s} \phi_T \, v \, d\Gamma = -- \int_{\Gamma_s} \Bigg[ -K_{r0} \exp\!\left(-\frac{E_{Kr}}{k_B T}\right) c_T^2 -+ K_{r0}^\ast \exp\!\left(-\frac{E_{Kr}^\ast}{k_B T}\right) c_{H_2} c_T -\Bigg] v \, d\Gamma -$$ - These fluxes can be implemented in FESTIM using `ParticleFluxBC` with user-defined expressions for each reaction term, as shown below. +++ @@ -344,7 +273,7 @@ def my_custom_recombination_flux(c, T): ``` ```{note} -For more complex isotopic exchanges, we can also use `SurfaceReactionBC` to define recombination fluxes. See [modeling isotopic exchange for multiple species](examples.md) to learn more. +For more complex isotopic exchanges, we can also use `SurfaceReactionBC` add other reactions. See [modeling isotopic exchange for multiple species](examples.md) to learn more. ``` +++ @@ -443,7 +372,7 @@ grid2.point_data["c"] = tritium.post_processing_solution.x.array grid2.set_active_scalars("c") ``` -Let's plot both `data_with_H2` and `data_no_H2` to see the concentraton profile for each case: +Let's plot both `data_with_H2` and `data_no_H2` to see the concentration profile for each case: ```{code-cell} ipython3 :tags: [hide-input] @@ -457,6 +386,7 @@ plt.plot(t, data_with_H2, label="With H_2") plt.xlabel("Time (s)") plt.ylabel("Surface Flux (right)") +plt.yscale("log") plt.legend() plt.show() ``` @@ -484,8 +414,8 @@ else: pyvista.set_jupyter_backend("html") plotter = pyvista.Plotter() - -plotter.add_mesh(grid2) +vmin, vmax = grid["c"].min(), grid2["c"].max() +plotter.add_mesh(grid2, clim=[vmin, vmax]) plotter.view_xy() if not pyvista.OFF_SCREEN: plotter.show() diff --git a/book/content/boundary_conditions/h_transport_basic.md b/book/content/boundary_conditions/h_transport_basic.md index f477a75d..4ca7d6ab 100644 --- a/book/content/boundary_conditions/h_transport_basic.md +++ b/book/content/boundary_conditions/h_transport_basic.md @@ -17,78 +17,11 @@ kernelspec: This section discusses how to implement basic boundary conditions (BCs) for hydrogen transport problems. Boundary conditions are essential to FESTIM simulations, as they describe the mathematical problem at the boundaries of the simulated domain. Read more about BCs _[here](https://festim.readthedocs.io/en/fenicsx/userguide/boundary_conditions.html)_. Objectives: -* Learn mathematical formulation for concentration and flux boundary conditions * Implement fixed concentration boundary conditions * Implement particle flux boundary conditions +++ -## Understanding math behind concentration and flux boundary conditions - -### Fixed concentration - -A fixed concentration (Dirichlet) boundary condition prescribes the value of the mobile hydrogen isotope concentration at a boundary. This enforces the concentration to remain constant in time and space on the specified boundary, independent of the solution in the bulk. - -This boundary condition is typically used to represent surfaces in equilibrium with an infinite reservoir, imposed implantation conditions, or experimentally controlled concentrations. - -#### Mathematical formulation - -On a boundary $\Gamma_D$, the mobile concentration satisfies - -$$ -c(\mathbf{x}, t) = c_0 \quad \text{for } \mathbf{x} \in \Gamma_D, -$$ - -where $c_0$ is the prescribed concentration value. - -In the weak formulation, this condition is enforced by directly constraining the degrees of freedom associated with the boundary. - -### Particle flux boundary conditions - -A particle flux (Neumann) boundary condition prescribes the normal flux of mobile hydrogen isotopes across a boundary. Unlike fixed concentration conditions, flux boundary conditions do not directly constrain the concentration value at the surface. Instead, they control the rate at which particles enter or leave the domain. - -Flux boundary conditions are commonly used to represent implantation from a plasma, outgassing to vacuum, permeation through a surface, or symmetry boundaries where no net transport occurs. - -#### Mathematical formulation - -The particle flux $\mathbf{J}$ is typically given by Fick’s law, - -$$ -\mathbf{J} = -D \nabla c, -$$ - -where $D$ is the diffusion coefficient and $c$ is the mobile concentration. - -On a boundary $\Gamma_N$, a prescribed normal flux is imposed as - -$$ -\mathbf{J} \cdot \mathbf{n} = g(\mathbf{x}, t) \quad \text{for } \mathbf{x} \in \Gamma_N, -$$ - -where: -- $\mathbf{n}$ is the outward unit normal vector, -- $g(\mathbf{x}, t)$ is the imposed particle flux (positive for outward flux, negative for inward flux). - -A special case is the **zero-flux (symmetry) boundary condition**, given by - -$$ -\mathbf{J} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma, -$$ - -which implies no net particle transport across the boundary. - -#### Weak formulation - -In the weak form, flux boundary conditions appear naturally as surface integrals after integration by parts. For a test function $v$, the boundary contribution is - -$$ -\int_{\Gamma_N} g(\mathbf{x}, t)\, v \, \mathrm{d}\Gamma. -$$ - -Because of this, flux boundary conditions are often referred to as *natural boundary conditions* and do not require explicit modification of the solution space, unlike Dirichlet conditions. - -+++ - ## Imposing fixed concentration boundary conditions Users can prescribe a fixed concentration on a boundary using `FixedConcentrationBC`, which can depend on temperature, time, and space. @@ -116,19 +49,6 @@ my_model.temperature = 400 my_model.settings = F.Settings(atol=1e-8, rtol=1e-8, transient=False) ``` -We can define the top and right boundary conditions using `lambda` functions: - -```{code-cell} ipython3 -top_bc_function = lambda x: 2.0 * x[0] + x[1] -right_bc_function = lambda x: x[1] ** 2 + 1.0 -``` - -```{note} -`x[0]`, `x[1]`, `x[2]` corresponse to *(x, y, z)* in cartesian space. Additionally, boundary conditions can also be given as Python functions (see *{ref}`label-time-temp-concentration`* to learn more.) -``` - -+++ - Now, let's create the boundaries to assign the BCs: ```{code-cell} ipython3 @@ -141,6 +61,19 @@ vol = F.VolumeSubdomain(id=5, material=mat) my_model.subdomains = [top_subdomain, right_subdomain, left_subdomain, bottom_subdomain, vol] ``` +We can define the top and right boundary conditions using `lambda` functions: + +```{code-cell} ipython3 +top_bc_function = lambda x: 2.0 * x[0] + x[1] +right_bc_function = lambda x: x[1] ** 2 + 1.0 +``` + +```{note} +`x[0]`, `x[1]`, `x[2]` correspond to the first, second, and third space coordinate ( *(x, y, z)* in cartesian space). See *{ref}`label-time-temp-concentration`* to learn more. +``` + ++++ + To include these boundary conditions to our problem, we use `FixedConcentrationBC`. We must also specify which subdomain (boundary) each BC is applied to, as well as the corresponding species: ```{code-cell} ipython3 @@ -294,6 +227,12 @@ H = F.Species(name="H") J = lambda t, c: 10*t**2 + 2*c ``` +```{note} +A flux boundary condition dependent on the actual solution is called a Robin boundary condition. +``` + ++++ + To add this BC, we need to create a dictionary that maps the concentration variable `c` in our custom function to our species `H`: ```{code-cell} ipython3 From fe6fd92e0fba3bc74d1ae9178e78d654060cdc34 Mon Sep 17 00:00:00 2001 From: ck768 Date: Mon, 23 Feb 2026 13:55:23 -0500 Subject: [PATCH 14/17] minor fixes and forgot intro file --- book/content/boundary_conditions/bc_intro.md | 79 +++++++++++++++++++ book/content/boundary_conditions/examples.md | 2 +- .../h_transport_advanced.md | 4 +- 3 files changed, 82 insertions(+), 3 deletions(-) create mode 100644 book/content/boundary_conditions/bc_intro.md diff --git a/book/content/boundary_conditions/bc_intro.md b/book/content/boundary_conditions/bc_intro.md new file mode 100644 index 00000000..b21ac271 --- /dev/null +++ b/book/content/boundary_conditions/bc_intro.md @@ -0,0 +1,79 @@ +--- +jupytext: + formats: ipynb,md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.19.1 +--- + +# Fixed value and flux boundary condition mathematics + +This section discusses the math behind fixed temperature/concentration and flux boundary conditions (BCs). + ++++ + +## Fixed concentration/temperature boundary conditions + +A fixed concentration (Dirichlet) boundary condition prescribes the value of the mobile hydrogen isotope concentration at a boundary. This enforces the concentration to remain constant in time and space on the specified boundary, independent of the solution in the bulk. This also applies to a fixed temperature boundary condition. + +For hydrogen transport, this boundary condition is typically used to represent surfaces in equilibrium with a gas, imposed implantation conditions, or experimentally controlled concentrations. Fixed temperatures are commonly used to model infinite heat reservoirs or sinks. + +### Mathematical formulation + +On a boundary $\Gamma_D$, the mobile concentration satisfies + +$$ +c(\mathbf{x}, t) = c_0 \quad \text{for } \mathbf{x} \in \Gamma_D, +$$ + +where $c_0$ is the prescribed concentration value. + +This condition is enforced by directly constraining the degrees of freedom associated with the boundary. + +## Flux boundary conditions + +A particle flux (Neumann) boundary condition prescribes the normal flux of mobile hydrogen isotopes across a boundary. Unlike fixed concentration conditions, flux boundary conditions do not directly constrain the concentration value at the surface. Instead, they control the rate at which particles enter or leave the domain. + +Flux boundary conditions are commonly used to represent implantation from a plasma, outgassing to vacuum, permeation through a surface, or symmetry boundaries where no net transport occurs. + +This also applies for heat flux boundary conditions. + +### Mathematical formulation + +The particle flux $\mathbf{J}$ is typically given by Fick’s law, + +$$ +\mathbf{J} = -D \nabla c, +$$ + +where $D$ is the diffusion coefficient and $c$ is the mobile concentration. + +On a boundary $\Gamma_N$, a prescribed normal flux is imposed as + +$$ +\mathbf{J} \cdot \mathbf{n} = g(\mathbf{x}, t) \quad \text{for } \mathbf{x} \in \Gamma_N, +$$ + +where: +- $\mathbf{n}$ is the outward unit normal vector, +- $g(\mathbf{x}, t)$ is the imposed particle flux (positive for outward flux, negative for inward flux). + +A special case is the **zero-flux (symmetry) boundary condition**, given by + +$$ +\mathbf{J} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma, +$$ + +which implies no net particle transport across the boundary. + +#### Weak formulation + +In the weak form, flux boundary conditions appear naturally as surface integrals after integration by parts. For a test function $v$, the boundary contribution is + +$$ +\int_{\Gamma_N} g(\mathbf{x}, t)\, v \, \mathrm{d}\Gamma. +$$ + +Because of this, flux boundary conditions are often referred to as *natural boundary conditions* and do not require explicit modification of the solution space, unlike Dirichlet conditions. diff --git a/book/content/boundary_conditions/examples.md b/book/content/boundary_conditions/examples.md index 2d782e40..63d3b72e 100644 --- a/book/content/boundary_conditions/examples.md +++ b/book/content/boundary_conditions/examples.md @@ -12,7 +12,7 @@ kernelspec: name: python3 --- -# Advanced examples +# Hydrogen transport: examples This section discusses a few advanced examples using various boundary conditions to help users understand how BCs can help model their systems. diff --git a/book/content/boundary_conditions/h_transport_advanced.md b/book/content/boundary_conditions/h_transport_advanced.md index 8c9b8dee..f76d8db3 100644 --- a/book/content/boundary_conditions/h_transport_advanced.md +++ b/book/content/boundary_conditions/h_transport_advanced.md @@ -330,7 +330,7 @@ surface_flux = F.SurfaceFlux(field=tritium, surface=right_subdomain) ```{code-cell} ipython3 my_model.temperature = 300 -my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, stepsize=50, final_time=500) +my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, stepsize=50, final_time=200) my_model.exports = [surface_flux] my_model.initialise() my_model.run() @@ -372,7 +372,7 @@ grid2.point_data["c"] = tritium.post_processing_solution.x.array grid2.set_active_scalars("c") ``` -Let's plot both `data_with_H2` and `data_no_H2` to see the concentration profile for each case: +Let's plot both `data_with_H2` and `data_no_H2` to see the flux versus time for each case: ```{code-cell} ipython3 :tags: [hide-input] From 25c6f105e3f23d7ba0e055cc106987c106a11ef3 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 23 Feb 2026 16:43:30 -0500 Subject: [PATCH 15/17] improved isotopic exchange case --- .../h_transport_advanced.md | 250 +++++++++++------- 1 file changed, 151 insertions(+), 99 deletions(-) diff --git a/book/content/boundary_conditions/h_transport_advanced.md b/book/content/boundary_conditions/h_transport_advanced.md index f76d8db3..2de84c30 100644 --- a/book/content/boundary_conditions/h_transport_advanced.md +++ b/book/content/boundary_conditions/h_transport_advanced.md @@ -1,11 +1,11 @@ --- jupytext: - formats: ipynb,md:myst + formats: md:myst,ipynb text_representation: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.19.1 + jupytext_version: 1.20.0 kernelspec: display_name: festim-workshop language: python @@ -222,157 +222,242 @@ Using `SurfaceReactionBC` is just a convenient way for users to define surface f Isotopic exchange occurs when hydrogenic isotopes swap positions, for example: -$$ +```{math} +:label: t2_recomb \mathrm{T + T} \rightleftharpoons \mathrm{T_2} -\\ +``` + +```{math} +:label: isotopic_exchange \mathrm{T + H_2} \rightleftharpoons \mathrm{H} + \mathrm{HT} -$$ +``` ### Mathematical formulation Let: -- $c_T, c_{H_2}, c_{HT}$ be the surface concentrations of tritium, molecular hydrogen, and hydrogen-tritium -- $K_{r0}, K_{r0}^\ast$ the pre-exponential factors -- $E_{Kr}, E_{Kr}^\ast$ the activation energies +- $c_T$ be the surface concentrations of tritium +- $P_\mathrm{H_2}, P_\mathrm{HT}$ the partial pressures of $\mathrm{H_2}$ and $\mathrm{HT}$ +- $K_{r}$ the rate coefficient of reaction {eq}`t2_recomb` +- $K_{r}^\ast$ the rate coefficient of reaction {eq}`:label: isotopic_exchange` + +$$ +K_{r} = K_{r0} \exp\!\left(-\frac{E_{Kr}}{k_B T}\right) +$$ + +$$ +K_{r}^\ast = K_{r0}^\ast \exp\!\left(-\frac{E_{Kr}^\ast}{k_B T}\right) +$$ If $c_{H_2} \gg c_{HT}$, the flux of tritium is approximated as: $$ -\phi_T = -- K_{r0} \exp\!\left(-\frac{E_{Kr}}{k_B T}\right) c_T^2 -- K_{r0}^\ast \exp\!\left(-\frac{E_{Kr}^\ast}{k_B T}\right) c_{H_2} c_T +\phi_T = \phi_\mathrm{recombination} + \phi_\mathrm{exchange} +$$ +with +$$ +\phi_\mathrm{recombination} = - K_{r} c_T^2 \\ +\phi_\mathrm{exchange} = - K_{r}^\ast P_{H_2} c_T $$ These fluxes can be implemented in FESTIM using `ParticleFluxBC` with user-defined expressions for each reaction term, as shown below. +++ -### Modeling isotopic exchange as a recombination flux - -Isotopic exchange reactions can be modeled as user-defined expressions using `ufl`. - -Assuming a large, constant concentration of molecular hydrogen $H_2$ at a boundary, we can define a recombination flux using rate and energy coefficients: - -```{code-cell} ipython3 -import ufl -import festim as F - -Kr_0 = 1.0 -E_Kr = 0.1 - -def my_custom_recombination_flux(c, T): - Kr_0_custom = 1.0 - E_Kr_custom = 0.5 # eV - h2_conc = 1e25 # assumed constant H2 concentration in +```{note} +$K_{r0}^\ast$ and $K_{r0}$ have different dimensions since $P_\mathrm{H_2}$ is a pressure. +``` - recombination_flux = ( - -(Kr_0 * ufl.exp(-E_Kr / (F.k_B * T))) * c**2 - - (Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))) * h2_conc * c - ) - return recombination_flux +```{note} +This example is not tracking the dissolution and transport of protium in the metal. And assumes that the partial presure of $\mathrm{H_2}$ is constant. ``` ++++ + ```{note} For more complex isotopic exchanges, we can also use `SurfaceReactionBC` add other reactions. See [modeling isotopic exchange for multiple species](examples.md) to learn more. ``` +++ +### Implementation + ++++ + Let's work through a 1D example to illustrate the effect of isotopic exchange. We'll consider tritium diffusion from left to right, with a large partial pressure of $H_2$ at the right boundary (where isotopic exchange can occur). First, we'll set up our mesh and materials: ```{code-cell} ipython3 import numpy as np -from dolfinx.mesh import create_unit_square -from mpi4py import MPI my_model = F.HydrogenTransportProblem() -my_model.mesh = F.Mesh(create_unit_square(MPI.COMM_WORLD, 300, 300)) +my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 100)) -right_subdomain = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1)) -left_subdomain = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0)) +right_subdomain = F.SurfaceSubdomain1D(id=2, x=1) +left_subdomain = F.SurfaceSubdomain1D(id=3, x=0) material = F.Material(D_0=1e-2, E_D=0) -vol = F.VolumeSubdomain(id=5, material=material) +vol = F.VolumeSubdomain1D(id=5, material=material, borders=[0, 1]) my_model.subdomains = [vol, left_subdomain, right_subdomain] ``` +Isotopic exchange reactions can be modeled as user-defined expressions using `ufl`. + +```{code-cell} ipython3 +import ufl +import festim as F + +Kr_0_custom = 1 +E_Kr_custom = 0.0 # eV +P_H2 = 10 # assumed constant H2 concentration in + +def isotopic_exchange_flux_func(c, T): + return - (Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))) * P_H2 * c +``` + Now we'll add our species `tritium` to our model and define a custom flux to represent isotopic exchange. We use `ParticleFluxBC` and the custom flux function we defined above to define this BC on the right boundary: ```{code-cell} ipython3 tritium = F.Species("T") my_model.species = [tritium] -my_custom_flux = F.ParticleFluxBC( - value=my_custom_recombination_flux, +isotopic_exchange_flux = F.ParticleFluxBC( + value=isotopic_exchange_flux_func, subdomain=right_subdomain, species_dependent_value={"c": tritium}, species=tritium, ) ``` +```{code-cell} ipython3 +t2_recomb = F.SurfaceReactionBC( + reactant=[tritium, tritium], + gas_pressure=0, + k_r0=1e-22, + E_kr=0, + k_d0=0, + E_kd=0, + subdomain=right_subdomain, +) +``` + For diffusion across the mesh, we also define a fixed concentration on the left surface: ```{code-cell} ipython3 my_model.boundary_conditions = [ - my_custom_flux, + isotopic_exchange_flux, + t2_recomb, F.FixedConcentrationBC(subdomain=left_subdomain, value=1e20, species=tritium), ] ``` -We'll export the flux using `SurfaceFlux` (see [exporting derived quantities to learn more](../post_process/derived.md)) and save this data to a variable `data_with_H2`: +We'll export the flux using `SurfaceFlux` (see [exporting derived quantities to learn more](../post_process/derived.md)) and save this data to a variable `data_with_H2`. We also create a `Profile1DExport` to be able to plot the concentration profile later: ```{code-cell} ipython3 surface_flux = F.SurfaceFlux(field=tritium, surface=right_subdomain) + +profile = F.Profile1DExport(field=tritium, subdomain=vol) ``` ```{code-cell} ipython3 my_model.temperature = 300 -my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, stepsize=50, final_time=200) -my_model.exports = [surface_flux] +my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, stepsize=1, final_time=120) + +my_model.exports = [surface_flux, profile] my_model.initialise() my_model.run() -data_with_H2 = (surface_flux.data) +data_with_H2 = surface_flux.data.copy() ``` +Let's plot the evolution of the tritium concentration profile: + ```{code-cell} ipython3 -:tags: [hide-input, hide-output] +:tags: [hide-input] -import pyvista -from dolfinx import plot +import matplotlib.pyplot as plt -topology, cell_types, geometry = plot.vtk_mesh(tritium.post_processing_solution.function_space) -grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) -grid.point_data["c"] = tritium.post_processing_solution.x.array -grid.set_active_scalars("c") +import matplotlib.animation as animation +from IPython.display import HTML + +# Turn off interactive plotting to prevent static display +plt.ioff() + +# Create the figure and axis +fig, ax = plt.subplots(figsize=(10, 6)) + +# Set up the plot limits and labels +ax.set_ylabel("Tritium concentration") +ax.set_xlabel("Position") + +# Initialize empty line objects for the animation +line_current, = ax.plot([], [], color="black", linewidth=2, zorder=1000) +lines_previous = [] + +def animate(frame): + # Clear previous iteration lines + for line in lines_previous: + line.remove() + lines_previous.clear() + + # Plot all previous iterations in light grey + for i in range(frame): + line_prev, = ax.plot(profile.x, profile.data[i], color="lightgrey", linewidth=0.8, alpha=0.6) + lines_previous.append(line_prev) + + # Plot current iteration in black + line_current.set_data(profile.x, profile.data[frame]) + + ax.set_title(f"Time: {profile.t[frame]:.2f}") + + return [line_current] + lines_previous + +# Create animation +anim = animation.FuncAnimation( + fig, animate, frames=len(profile.data), interval=100, blit=False, repeat=True +) + +# Close the figure to prevent static display +plt.close(fig) + +# Display only the animation +HTML(anim.to_jshtml()) ``` -We can compare the surface flux to the case where we have no isotopic exchange by removing the custom boundary condition and saving the results to `data_no_H2`: +We can compare the surface flux to the case where we have no isotopic exchange by removing the isotopic exchange boundary condition and saving the results to `data_no_H2`: ```{code-cell} ipython3 my_model.boundary_conditions = [ F.FixedConcentrationBC(subdomain=left_subdomain, value=1e20, species=tritium), + t2_recomb, ] -my_model.exports = [surface_flux] +profile = F.Profile1DExport(field=tritium, subdomain=vol) +my_model.exports = [surface_flux, profile] my_model.initialise() my_model.run() -data_no_H2 = (surface_flux.data) +data_no_H2 = surface_flux.data.copy() ``` +Without isotopic exchange the gradient is much lower, indicating that isotopic exchange helps enhance surface outgassing of tritium! + ```{code-cell} ipython3 -:tags: [hide-input, hide-output] +:tags: [hide-input] + +# Create animation +anim = animation.FuncAnimation( + fig, animate, frames=len(profile.data), interval=100, blit=False, repeat=True +) -topology, cell_types, geometry = plot.vtk_mesh(tritium.post_processing_solution.function_space) -grid2 = pyvista.UnstructuredGrid(topology, cell_types, geometry) -grid2.point_data["c"] = tritium.post_processing_solution.x.array -grid2.set_active_scalars("c") +# Close the figure to prevent static display +plt.close(fig) + +# Display only the animation +HTML(anim.to_jshtml()) ``` -Let's plot both `data_with_H2` and `data_no_H2` to see the flux versus time for each case: +Let's plot both `data_with_H2` and `data_no_H2` to see the temporal evolution of the surface flux for each case: ```{code-cell} ipython3 :tags: [hide-input] @@ -382,45 +467,12 @@ import matplotlib.pyplot as plt x = my_model.mesh.mesh.geometry.x[:, 0] t = surface_flux.t plt.plot(t, data_no_H2, label="No H_2") -plt.plot(t, data_with_H2, label="With H_2") +plt.plot(t, data_with_H2,label="With H_2") -plt.xlabel("Time (s)") +plt.xlabel("Time") plt.ylabel("Surface Flux (right)") -plt.yscale("log") -plt.legend() +plt.legend(reverse=True) plt.show() ``` -We see that the flux is higher with isotopic exchange, as we'd expect. Let's also take a look at the concentration profiles (top with isotopic exchange and bottom without): - -```{code-cell} ipython3 -:tags: [hide-input] - -pyvista.set_jupyter_backend("html") - -plotter = pyvista.Plotter() - -plotter.add_mesh(grid) -plotter.view_xy() -if not pyvista.OFF_SCREEN: - plotter.show() -else: - figure = plotter.screenshot("concentration.png") -``` - -```{code-cell} ipython3 -:tags: [hide-input] - -pyvista.set_jupyter_backend("html") - -plotter = pyvista.Plotter() -vmin, vmax = grid["c"].min(), grid2["c"].max() -plotter.add_mesh(grid2, clim=[vmin, vmax]) -plotter.view_xy() -if not pyvista.OFF_SCREEN: - plotter.show() -else: - figure = plotter.screenshot("concentration.png") -``` - -The results without isotopic exchange show virtually no diffusion for a given inlet concentration, indicating that isotopic exchange helps enchance diffusion! +We see that the flux is higher with isotopic exchange, as we'd expect. From d903923f5f8fbd151bb6f2d7c1801c2fcff86fa9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 23 Feb 2026 17:14:17 -0500 Subject: [PATCH 16/17] completed isotopic exchange example --- book/content/boundary_conditions/examples.md | 87 ++++++++++++-------- 1 file changed, 53 insertions(+), 34 deletions(-) diff --git a/book/content/boundary_conditions/examples.md b/book/content/boundary_conditions/examples.md index 63d3b72e..979631f6 100644 --- a/book/content/boundary_conditions/examples.md +++ b/book/content/boundary_conditions/examples.md @@ -1,11 +1,11 @@ --- jupytext: - formats: ipynb,md:myst + formats: md:myst,ipynb text_representation: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.19.1 + jupytext_version: 1.20.0 kernelspec: display_name: festim-workshop language: python @@ -197,14 +197,26 @@ Using the approximation is computationally less expensive, and still provides si +++ -## Complex isotopic exchange with multple hydrogenic species ## +## Complex isotopic exchange with multiple hydrogenic species ## Surface reactions can involve multiple hydrogen isotopes, allowing for the modeling of complex isotope-exchange mechanisms between species. For example, in a system with both mobile hydrogen and deuteriun, various molecular recombination pathways may occur at the surface, resulting in the formation of $H_2$, $D_2$, and $HD$: -$$ \text{Reaction 1}: \mathrm{H+D} \rightarrow \mathrm{HD} \longrightarrow \phi_1 = K_{r1} c_H c_D - K_{d1}P_{HD} $$ -$$ \text{Reaction 2}: \mathrm{D+D} \rightarrow \mathrm{D_2} \longrightarrow \phi_2 = 2K_{r2} c_D^2 - K_{d2}P_{D2} $$ -$$ \text{Reaction 3}: \mathrm{D+H_2} \rightarrow \mathrm{HD + H} \longrightarrow \phi_3 = K_{r3} c_H c_D - K_{d3}P_{HD} $$ -$$ \text{Reaction 4}: \mathrm{H+H} \rightarrow \mathrm{H_2} \longrightarrow \phi_4 = 2K_{r4} c_H^2 - K_{d4}P_{H2} $$ +$$ \text{Reaction 1}: \mathrm{H+D} \rightleftharpoons \mathrm{HD} \longrightarrow R_1 = K_{r1} c_H c_D - K_{d1}P_{HD} $$ +$$ \text{Reaction 2}: \mathrm{D+D} \rightleftharpoons \mathrm{D_2} \longrightarrow R_2 = K_{r2} c_D^2 - K_{d2}P_{D2} $$ +$$ \text{Reaction 3}: \mathrm{H+H} \rightleftharpoons \mathrm{H_2} \longrightarrow R_3 = K_{r3} c_H^2 - K_{d3}P_{H2} $$ + +$$ \text{Reaction 4}: \mathrm{D+H_2} \rightleftharpoons \mathrm{HD + H} \longrightarrow R_4 = K_{r4} P_{H2} c_D - K_{d4}P_{HD} c_H $$ + + +From this reaction scheme, the surface fluxes of H and D are: + +$$ +\varphi_\mathrm{H} = -R_1 - 2 R_3 + R_4 +$$ + +$$ +\varphi_\mathrm{D} = -R_1 - 2 R_2 - R_4 +$$ Now consider the case where deuterium diffuses from left to right and reacts with background $\mathrm{H_2}$, while $\mathrm{P_{HD}}$ and $\mathrm{P_{D_2}}$ are negligible at the surface. @@ -242,9 +254,11 @@ H = F.Species("H") D = F.Species("D") my_model.species = [H, D] -H2 = F.SurfaceReactionBC( - reactant=[H, H], - gas_pressure=1000, +P_h2 = 1000 + +reaction_1_bc = F.SurfaceReactionBC( + reactant=[H, D], + gas_pressure=0, k_r0=1, E_kr=0.1, k_d0=1e-5, @@ -252,8 +266,8 @@ H2 = F.SurfaceReactionBC( subdomain=right_surf, ) -HD = F.SurfaceReactionBC( - reactant=[H, D], +reaction_2_bc = F.SurfaceReactionBC( + reactant=[D, D], gas_pressure=0, k_r0=1, E_kr=0.1, @@ -262,9 +276,9 @@ HD = F.SurfaceReactionBC( subdomain=right_surf, ) -D2 = F.SurfaceReactionBC( - reactant=[D, D], - gas_pressure=0, +reaction_3_bc = F.SurfaceReactionBC( + reactant=[H, H], + gas_pressure=P_h2, k_r0=1, E_kr=0.1, k_d0=1e-5, @@ -276,38 +290,43 @@ D2 = F.SurfaceReactionBC( Now, let's add our isotopic exchange reaction using `ParticleFluxBC` (see [here](h_transport_advanced.md) to learn more about defining isotopic exchange fluxes): ```{code-cell} ipython3 -import ufl +Kr_0_custom = 10000.0 +E_Kr_custom = 0.5 # eV -Kr_0 = 1.0 -E_Kr = 0.1 +def K_exchange(T): + return Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T)) -def my_custom_recombination_flux(c, T): - Kr_0_custom = 1.0 - E_Kr_custom = 0.5 # eV - h2_conc = 1e25 # assumed constant H2 concentration in +def isotopic_exchange_D(c_D, T): + return - K_exchange(T) * P_h2 * c_D + +def isotopic_exchange_H(c_D, T): + return + K_exchange(T) * P_h2 * c_D - recombination_flux = ( - -(Kr_0 * ufl.exp(-E_Kr / (F.k_B * T))) * c**2 - - (Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))) * h2_conc * c - ) - return recombination_flux -HDH = F.ParticleFluxBC( - value=my_custom_recombination_flux, +reaction_4_D = F.ParticleFluxBC( + value=isotopic_exchange_D, subdomain=right_surf, - species_dependent_value={"c": D}, + species_dependent_value={"c_D": D}, species=D, ) + +reaction_4_H = F.ParticleFluxBC( + value=isotopic_exchange_H, + subdomain=right_surf, + species_dependent_value={"c_D": D}, + species=H, +) ``` Finally, we add our boundary conditions and solve the steady-state problem: ```{code-cell} ipython3 my_model.boundary_conditions = [ - H2, - D2, - HD, - HDH, + reaction_1_bc, + reaction_2_bc, + reaction_3_bc, + reaction_4_D, + reaction_4_H, F.FixedConcentrationBC(subdomain=left_surf, value=1, species=D), F.FixedConcentrationBC(subdomain=left_surf, value=0, species=H), ] From 6415423d31178a9fd8aeda3869e7365f66dd0e5d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 23 Feb 2026 21:38:06 -0500 Subject: [PATCH 17/17] fixed eq ref --- book/content/boundary_conditions/h_transport_advanced.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/book/content/boundary_conditions/h_transport_advanced.md b/book/content/boundary_conditions/h_transport_advanced.md index 2de84c30..78563756 100644 --- a/book/content/boundary_conditions/h_transport_advanced.md +++ b/book/content/boundary_conditions/h_transport_advanced.md @@ -238,7 +238,7 @@ Let: - $c_T$ be the surface concentrations of tritium - $P_\mathrm{H_2}, P_\mathrm{HT}$ the partial pressures of $\mathrm{H_2}$ and $\mathrm{HT}$ - $K_{r}$ the rate coefficient of reaction {eq}`t2_recomb` -- $K_{r}^\ast$ the rate coefficient of reaction {eq}`:label: isotopic_exchange` +- $K_{r}^\ast$ the rate coefficient of reaction {eq}`isotopic_exchange` $$ K_{r} = K_{r0} \exp\!\left(-\frac{E_{Kr}}{k_B T}\right)