{ "cells": [ { "cell_type": "markdown", "id": "e48786b3-d807-411a-b390-4e12d30d8ec0", "metadata": {}, "source": [ "# Deep dive: evaluate state using FEniCS" ] }, { "cell_type": "markdown", "id": "57b933ef-56e6-417d-a5c1-b7e51ca75b3c", "metadata": {}, "source": [ "In this notebook we solve the model via the finite element method using [FEniCSx](https://fenicsproject.org/)." ] }, { "cell_type": "code", "execution_count": 1, "id": "972befab-1105-4c27-9e6a-3ee7bb924f94", "metadata": {}, "outputs": [], "source": [ "import pathlib\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import ufl\n", "from basix.ufl import element, mixed_element\n", "from dolfinx import default_real_type, default_scalar_type, fem, mesh\n", "from dolfinx.fem.petsc import NonlinearProblem\n", "from dolfinx.io import XDMFFile\n", "from mpi4py import MPI\n", "from petsc4py.PETSc import ScalarType\n", "import plotly.graph_objects as go\n", "import plotly.io as pio\n", "from plotly.subplots import make_subplots" ] }, { "cell_type": "markdown", "id": "8603b6ba-4699-4350-b9c3-afc1dc92f332", "metadata": {}, "source": [ "The following line allow the plot at the end" ] }, { "cell_type": "code", "execution_count": 2, "id": "c01c68e0-371f-439b-a00b-3f6606de13c1", "metadata": {}, "outputs": [], "source": [ "pio.renderers.default = \"notebook_connected\"" ] }, { "cell_type": "markdown", "id": "dac73127-2504-465d-ac9b-64dbb8ab50a2", "metadata": {}, "source": [ "## Choice of geometry" ] }, { "cell_type": "markdown", "id": "d605b497-8581-4f60-beb2-3e41c567b04b", "metadata": {}, "source": [ "In this notebook we fix some of the parameters arbitrarily." ] }, { "cell_type": "markdown", "id": "0b7f612c-93c2-4012-a40b-e821438031e6", "metadata": {}, "source": [ "- The final time $T$ is 2 (representing seconds)." ] }, { "cell_type": "code", "execution_count": 3, "id": "90ebd715-94d2-4292-98e1-e6043123c158", "metadata": {}, "outputs": [], "source": [ "T = 2.0" ] }, { "cell_type": "markdown", "id": "82bca5eb-4aaf-4f66-bc36-e9682691e15c", "metadata": {}, "source": [ "- The space geometry (or domain) $\\Omega$ is the 1D unit interval, i.e. $(0,1)$." ] }, { "cell_type": "code", "execution_count": 4, "id": "f1c0cc8e-fece-495f-bc8a-8ae91fb5c1d7", "metadata": {}, "outputs": [], "source": [ "x0 = 0.0\n", "xL = 1.0" ] }, { "cell_type": "markdown", "id": "795f5985-3787-432e-84be-b259dbad7db1", "metadata": {}, "source": [ "- The diffusion coefficients $\\kappa_1$ and $\\kappa_2$ are constant at 1." ] }, { "cell_type": "code", "execution_count": 5, "id": "0985a736-b3a7-4235-9f3b-dac622206b0f", "metadata": {}, "outputs": [], "source": [ "def kappa_1(x):\n", " return np.ones_like(x[0])\n", "\n", "def kappa_2(x):\n", " return np.ones_like(x[0])" ] }, { "cell_type": "markdown", "id": "0c1c2e68-77b7-4758-938c-4d30535c5166", "metadata": {}, "source": [ "- The initial value $y_\\circ$ is constant at 5." ] }, { "cell_type": "code", "execution_count": 6, "id": "4be876ea-8a10-4dbe-b007-8f05b4b43cbb", "metadata": {}, "outputs": [], "source": [ "def y0(x):\n", " return 5 * np.ones_like(x[0])" ] }, { "cell_type": "markdown", "id": "c665facd-5f06-4f68-b451-111b61d5788e", "metadata": {}, "source": [ "- The control vector is $\\boldsymbol{\\mu} = [3, 3, 3, 3]$." ] }, { "cell_type": "code", "execution_count": 7, "id": "eeb4253d-9d98-4776-8016-f36bf183c86b", "metadata": {}, "outputs": [], "source": [ "mu = np.array([3.0, 3.0, 3.0, 3.0])" ] }, { "cell_type": "markdown", "id": "81203e1f-ce2c-4b79-b30f-26815c722e6a", "metadata": {}, "source": [ "- As the input function $u$ we take:\n", " $$u(t) := \\begin{cases} -3 & t \\le 4/3 \\\\ +3 & 4/3 < t \\le 2 \\end{cases}$$" ] }, { "cell_type": "code", "execution_count": 8, "id": "205d1d40-6b00-4cb1-863a-8a70e7974788", "metadata": {}, "outputs": [], "source": [ "def u_func(t):\n", " return -3 * (t <= 4 / 3) + 3 * (t > 4 / 3) * (t <= 2.0)" ] }, { "cell_type": "markdown", "id": "ce4fe36a-cad0-4d6f-8f8b-00b0591b1486", "metadata": {}, "source": [ "## Finite Element discretization" ] }, { "cell_type": "markdown", "id": "ba75ac0a-8431-4fd3-9d3e-2300b7938fda", "metadata": {}, "source": [ "To solve this problem numerically we use:" ] }, { "cell_type": "markdown", "id": "28545774-b3ad-4d58-93f9-de680c4d4f8b", "metadata": {}, "source": [ "- linear Lagrangian elements in space on 201 node points over the interval $\\Omega$." ] }, { "cell_type": "code", "execution_count": 9, "id": "07b97c63-559a-414a-8f91-76def8a536bf", "metadata": {}, "outputs": [], "source": [ "nx = 201\n", "degree = 1\n", "K = 201" ] }, { "cell_type": "markdown", "id": "7c48f68d-b0e5-4678-88ef-bdaab6fac917", "metadata": {}, "source": [ "- implicit Euler method in time in 201 time steps (counting initial and final time)" ] }, { "cell_type": "code", "execution_count": 10, "id": "82d8ccdf-7181-4369-a277-30703728cc9c", "metadata": {}, "outputs": [], "source": [ "dt = T / (K - 1)" ] }, { "cell_type": "markdown", "id": "c2de70c6-105b-4399-890d-fee3da895927", "metadata": {}, "source": [ "## Solution with FEniCS step by step" ] }, { "cell_type": "markdown", "id": "a9176fbf-d5cd-4072-938f-9d23671139c3", "metadata": {}, "source": [ "We first create the discretized space and the function space for the state $(y,q)$:" ] }, { "cell_type": "code", "execution_count": 11, "id": "185d1593-2b06-4132-8ed1-a683fbe2c445", "metadata": {}, "outputs": [], "source": [ "Omega = mesh.create_interval(\n", " comm=MPI.COMM_WORLD,\n", " nx=nx-1, # this function wants the number of intervals rather than the number of nodes\n", " points=[x0, xL],\n", ")\n", "\n", "P1 = element(\"Lagrange\", Omega.basix_cell(), degree)\n", "V = fem.functionspace(Omega, mixed_element([P1, P1]))" ] }, { "cell_type": "markdown", "id": "69609a11-16bf-4c95-9b50-ff35aa80084b", "metadata": {}, "source": [ "As $V$ is the solution space for the coupled state $(y,q)$, we also define the individual subspaces:" ] }, { "cell_type": "code", "execution_count": 12, "id": "affad53b-3555-4324-a507-fa1890d4eb07", "metadata": {}, "outputs": [], "source": [ "V_sub_0 = V.sub(0)\n", "Vy, Vy_to_V_sub_0 = V_sub_0.collapse()\n", "\n", "V_sub_1 = V.sub(1)\n", "Vq, Vq_to_V_sub_1 = V_sub_1.collapse()" ] }, { "cell_type": "markdown", "id": "cb37834b-6c81-4bfe-bc67-de139984d4ab", "metadata": {}, "source": [ "The solution at the current iteration $(y^n, q^n)$ is called `v` and the solution at the previous iteration $(y^{n-1}, q^{n-1})$ is called `v_old`. Here we initialize the two states in the joint function space." ] }, { "cell_type": "code", "execution_count": 13, "id": "0d5974be-4a38-49a3-a8ef-baa981144971", "metadata": {}, "outputs": [], "source": [ "v = fem.Function(V) # current solution\n", "v_old = fem.Function(V) # solution from previous converged step\n", "\n", "# Split mixed functions\n", "y, q = ufl.split(v)\n", "y_old, q_old = ufl.split(v_old)" ] }, { "cell_type": "markdown", "id": "c7f41197-172d-4469-820b-391afc9c04a2", "metadata": {}, "source": [ "The nonlinearity is defined as follow:" ] }, { "cell_type": "code", "execution_count": 14, "id": "7dae0356-6bdb-47a3-83a2-9e3370599451", "metadata": {}, "outputs": [], "source": [ "y_ = ufl.variable(y)\n", "q_ = ufl.variable(q)\n", "f = ufl.sqrt(y_) * ufl.sinh(q_)" ] }, { "cell_type": "markdown", "id": "3539db18-6af9-4297-a51f-fef0b63cb071", "metadata": {}, "source": [ "The initial condition should only be applied to the \"y-part\" of the state `v`. Here is how this is done." ] }, { "cell_type": "code", "execution_count": 15, "id": "efa10b10-fd77-4183-ba3e-e5c2510a123e", "metadata": {}, "outputs": [], "source": [ "v.x.array[:] = 0.0\n", "v.sub(0).interpolate(y0)\n", "v.x.scatter_forward()" ] }, { "cell_type": "markdown", "id": "ea017046-5733-4c68-92d0-9c097c7ee18b", "metadata": {}, "source": [ "As in the weak formulation the Dirichlet boundary conditions only determine the space where the solution is living, here they are collected in a list of conditions `bcs`:" ] }, { "cell_type": "code", "execution_count": 16, "id": "c6732f75-1482-4ac9-800d-0315bd829c49", "metadata": {}, "outputs": [], "source": [ "def on_left_boundary(x):\n", " return np.isclose(x[0], x0)\n", "\n", "\n", "boundary_dofs = fem.locate_dofs_geometrical(Vq, on_left_boundary)\n", "\n", "bc0 = fem.dirichletbc(default_scalar_type(0), Vq_to_V_sub_1[boundary_dofs], V_sub_1)\n", "bcs = [bc0]" ] }, { "cell_type": "markdown", "id": "8f65bab1-7071-48fe-b9b7-821a942d82a2", "metadata": {}, "source": [ "We interpolate the diffusion equations in the corresponding space:" ] }, { "cell_type": "code", "execution_count": 17, "id": "acff20bc-0449-4b39-84ca-dc659790c1e2", "metadata": {}, "outputs": [], "source": [ "k1 = fem.Function(Vy)\n", "k1.interpolate(kappa_1)\n", "\n", "k2 = fem.Function(Vq)\n", "k2.interpolate(kappa_2)" ] }, { "cell_type": "markdown", "id": "aa494be4-0e41-4d43-aa95-eb911154623f", "metadata": {}, "source": [ "We define the test functions and the nonlinear PDE operator `F` (where the two weak functions will be summed):" ] }, { "cell_type": "code", "execution_count": 18, "id": "da7a5029-dadb-48fa-bcf6-dbf88bd3c22b", "metadata": {}, "outputs": [], "source": [ "phi, psi = ufl.TestFunctions(V)\n", "b = fem.Constant(Omega, ScalarType(0.0))\n", "Fy = (\n", " ufl.inner(y, phi) * ufl.dx\n", " - ufl.inner(y_old, phi) * ufl.dx\n", " + dt * mu[0] * k1 * ufl.inner(ufl.grad(y), ufl.grad(phi)) * ufl.dx\n", " - dt * mu[1] * ufl.inner(f, phi) * ufl.dx\n", ")\n", "Fq = (\n", " mu[2] * k2 * ufl.inner(ufl.grad(q), ufl.grad(psi)) * ufl.dx\n", " + mu[3] * ufl.inner(f, psi) * ufl.dx\n", " + b * psi * ufl.ds\n", ")\n", "F = Fy + Fq" ] }, { "cell_type": "markdown", "id": "16d9c477-0d0b-4c38-b71a-493feedad9cd", "metadata": {}, "source": [ "Here we define the solver options and we initialize a `problem` object. This object is going to be solved in each iteration, while updating all the functions." ] }, { "cell_type": "code", "execution_count": 19, "id": "ec76ad9a-6ab2-4c04-8236-a891e4939924", "metadata": {}, "outputs": [], "source": [ "# Create nonlinear problem\n", "petsc_options = {\n", " \"snes_type\": \"newtonls\",\n", " \"snes_linesearch_type\": \"none\",\n", " \"snes_stol\": np.sqrt(np.finfo(default_real_type).eps) * 1e-2,\n", " \"snes_atol\": 0,\n", " \"snes_rtol\": 0,\n", " \"ksp_type\": \"preonly\",\n", " \"pc_type\": \"lu\",\n", " \"pc_factor_mat_solver_type\": \"mumps\",\n", "}\n", "problem = NonlinearProblem(\n", " F,\n", " v,\n", " bcs=bcs,\n", " petsc_options_prefix=\"fenics_FE_\",\n", " petsc_options=petsc_options,\n", ")" ] }, { "cell_type": "markdown", "id": "fa96eb8a-9653-452a-932a-fe650c614247", "metadata": {}, "source": [ "We collect data into the matrices `Y` and `Q`:\n", "\\begin{equation*}\n", " \\mathrm{Y} =\n", " \\left[\n", " \\begin{array}{c|c|c}\n", " y^1 & \\dots & y^K\n", " \\end{array}\n", " \\right]\n", " \\qquad\n", " \\mathrm{Q} = \n", " \\left[\n", " \\begin{array}{c|c|c}\n", " q^1 & \\dots & q^K\n", " \\end{array}\n", " \\right]\n", "\\end{equation*}\n", "These matrices will contain the columns of the coordinate arrays (in the FE basis)." ] }, { "cell_type": "code", "execution_count": 20, "id": "ab9dbcfd-5eb7-4cbe-b3e5-2276daec0702", "metadata": {}, "outputs": [], "source": [ "Y = np.zeros((nx, K-1))\n", "Q = np.zeros((nx, K-1))" ] }, { "cell_type": "markdown", "id": "a3eba9f0-94cc-490b-a341-5b1df7aa2225", "metadata": {}, "source": [ "The solver is advanced in time from to until a terminal time is reached" ] }, { "cell_type": "code", "execution_count": 21, "id": "073afcc2-35b0-4135-bca5-cf5353904abc", "metadata": { "editable": true, "scrolled": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "t_array = np.linspace(1, T, K) # array of time steps\n", "for k, t in enumerate(t_array[1:]): # we skip time=0\n", " v_old.x.array[:] = v.x.array\n", " b.value = u_func(t)\n", " problem.solve()\n", " Y[:, k] = v.x.array[Vy_to_V_sub_0]\n", " Q[:, k] = v.x.array[Vq_to_V_sub_1]" ] }, { "cell_type": "markdown", "id": "0f4d3fa4-5470-4bef-8a99-be173a7769bf", "metadata": {}, "source": [ "Next, let us plot the states in an interactive plot." ] }, { "cell_type": "code", "execution_count": 22, "id": "7a31d7fa-78d3-4460-8283-7b283221f814", "metadata": {}, "outputs": [ { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = go.Figure()\n", "\n", "# Array of x values\n", "x = Omega.geometry.x[:, 0]\n", "\n", "# Create subplots\n", "fig = make_subplots(\n", " rows=1,\n", " cols=2,\n", " column_widths=[0.5, 0.5],\n", ")\n", "\n", "# Add all $y^k$ and $q^k$ functions to trace\n", "for k, t in enumerate(t_array[1:]):\n", " fig.add_traces(\n", " [\n", " go.Scatter(\n", " x=x,\n", " y=Y[:,k],\n", " visible=False,\n", " name=f\"y(t={t:0.3f})\",\n", " ),\n", " go.Scatter(\n", " x=x,\n", " y=Q[:,k],\n", " visible=False,\n", " name=f\"q(t={t:0.3f})\",\n", " ),\n", " ],\n", " rows=[1,1],\n", " cols=[1,2],\n", " )\n", "\n", "# Make t=0 visible\n", "fig.data[0].visible = True\n", "fig.data[1].visible = True\n", "\n", "# Create and add slider\n", "steps = []\n", "for i in range(K-1):\n", " step = dict(\n", " method=\"update\",\n", " args=[{\"visible\": [False] * len(fig.data)},\n", " {\"title\": \"Slider switched to step: \" + str(i)}], # layout attribute\n", " )\n", " step[\"args\"][0][\"visible\"][2 * i] = True # Toggle i'th trace to \"visible\"\n", " step[\"args\"][0][\"visible\"][2 * i + 1] = True # Toggle i'th trace to \"visible\"\n", " steps.append(step)\n", "\n", "sliders = [dict(\n", " active=0,\n", " currentvalue={\"prefix\": \"Time: \"},\n", " pad={\"t\": 50},\n", " steps=steps\n", ")] \n", "\n", "# Add slider to figure\n", "fig.update_layout(\n", " sliders=sliders,\n", " yaxis={\"range\": [0, 7]},\n", " legend={\"x\": 0.5, \"y\":1.3, \"xanchor\": \"center\", \"orientation\":\"h\"},\n", ")\n", "\n", "# Show final interactive plot\n", "fig.show()" ] }, { "cell_type": "markdown", "id": "63a21f9e-559d-4eb5-901d-6eebfcc7fb57", "metadata": {}, "source": [ "We see that the state $y$ is \"inflated\" by the positive input $u$. When $u$ turns negative, $q$ changes to a more \"decrescent\" shape. This makes the state $y$ \"deflate\"." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.12.12" } }, "nbformat": 4, "nbformat_minor": 5 }