From dd542d49119c83059465efbe1885538a9f7531e6 Mon Sep 17 00:00:00 2001 From: ckhurana Date: Wed, 8 Oct 2025 10:47:52 -0400 Subject: [PATCH 1/4] 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 6cac4fa..40cd729 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 0000000..7109ce0 --- /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 0000000..da26f66 --- /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 2/4] 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 40cd729..30f4a6d 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 7109ce0..0000000 --- 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 0000000..c1538d5 --- /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 da26f66..0000000 --- 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 0000000..4f1984b --- /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 0000000..e2b8912 --- /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 0000000..1576291 --- /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 3/4] 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 c1538d5..ec09ff6 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 4f1984b..493dedd 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 e2b8912..e78df30 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 1576291..0a42c38 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 4/4] 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 0a42c38..8764565 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} $$