{ "cells": [ { "cell_type": "markdown", "id": "46b14116", "metadata": { "slideshow": { "slide_type": "skip" }, "toc": true }, "source": [ "

Table of Contents

\n", "
" ] }, { "cell_type": "markdown", "id": "ff47678c", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# A tutorial on solving singularly perturbed problems in Firedrake\n", "\n", "**21st Workshop on Numerical Methods for Problems with Layer Phenomena**\n", " \n", "[Niall Madden](https://www.niallmadden.ie/) and Sean Tobin, April 2025\n", "\n", "You can access this notebook at:\n", "* Jupyter notebook: https://www.niallmadden.ie/SPPs_Firedrake.ipynb\n", "* HTML version: https://www.niallmadden.ie/SPPs_Firedrake.html" ] }, { "cell_type": "markdown", "id": "04c10b1c", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "In this session, we introduce how to solve one-dimensional and two-dimensional singularly perturbed differential equations (SPDEs) on $(0,1)^d$, $d=1,2$ in Firedrake.\n", "\n", "For a very basic introduction to **Firedrake**, see https://niallmadden.ie/IntroToFiredrake. \n", "For more advanced topics, see https://www.firedrakeproject.org/documentation.html\n", "\n", "While this session will include the very fundamentals of using Firedrake for solving PDEs, there will be an emphasis on topics particular to SPDEs:\n", "* constructing layer-adapted meshes\n", "* error estimation in various norms\n", "* convergence tables\n", "* extension to higher dimensions" ] }, { "cell_type": "markdown", "id": "4b955393", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## What is Firedrake?\n", "\n", "According to the documentation: _Firedrake is an automated system for the solution of partial differential equations using the finite element method (FEM). Firedrake uses sophisticated code generation to provide mathematicians, scientists, and engineers with a very high productivity way to create sophisticated high performance simulations._\n", "\n", "Crucial:\n", "* Available as a Python module;\n", "* You can use one of the many built-in mesh generators, or define a mesh (by hand, or import from, say, gmsh);\n", "* Define a elements, and a function space on that mesh;\n", "* **Express the weak form of a PDE**\n", "* Call a solver\n", "* Compute errors and/or visualise solutions.\n" ] }, { "cell_type": "markdown", "id": "91331a42", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Background\n", "\n", "Firedrake began in 2014 is as a \"reimagining\" of `FEniCS`. \n", "* It uses the FEniCS _Unified Form Language_ (UFL) to express PDEs. \n", "\n", "* Boasts sophisticated, programmable solvers through \"_seamless_\" coupling with PETSc.\n", "\n", "* Triangular, quadrilateral, and tetrahedral unstructured meshes.\n", "\n", "* Vast array of finite element spaces.\n", "\n", "Etc. \n", "\n", "(To the user, Firedrake may seem a little closer the FEniCS than FEniCSx)." ] }, { "cell_type": "markdown", "id": "a474dbff", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Getting started on Google's colab\n", "\n", "For serious work, you should install Firedrake on your own system. That is not too hard to do: as of last month, there is a `pip install` approach.\n", "\n", "But to get started, we'll try Google colab: https://colab.google/\n", "\n", "Note the incantation needed to install Firedrake. Can take several minutes to run:" ] }, { "cell_type": "code", "execution_count": null, "id": "9cf89eda", "metadata": {}, "outputs": [], "source": [ "try:\n", " from firedrake import *\n", "except ImportError:\n", " !wget \"https://fem-on-colab.github.io/releases/firedrake-install-release-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n", " from firedrake import *" ] }, { "cell_type": "markdown", "id": "769aef83", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "We'll also use the `numpy`, `pandas`, and `matplotlib` modules:" ] }, { "cell_type": "code", "execution_count": null, "id": "777ef1ab", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "# %matplotlib notebook # may not work on colab\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "73903791", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## A one-dimensional problem" ] }, { "cell_type": "markdown", "id": "e325caa9", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "The first problem we'll solve is this ODE: _find $u(x)$ such that_\n", "$$-\\varepsilon^2 u'' + r(x)u = f(x) \\text{ on } (0,1), \\qquad\\text{ with } u(0)=u(1)=0.$$\n", "\n", "The weak form is $$\\underbrace{\\varepsilon^2 (u',v') + (ru,v)}_{a(u,v)} = \\underbrace{(f,v)}_{l(v)},$$\n", "where here $(u,v):=\\int_0^1 u(x)v(x) dx$, is the usual $L^2$ inner product.\n", "\n", "The specific example we'll consider will have $r(x) \\equiv 1$ and $f(x)=e^x$.\n", "\n", "We'll solve this problem using a standard $P_1$ (continuous) Galerkin method.\n" ] }, { "cell_type": "markdown", "id": "bb9ac524", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Problem data\n", "I like to use `lambda` functions for defining functions" ] }, { "cell_type": "code", "execution_count": null, "id": "b04162f8", "metadata": {}, "outputs": [], "source": [ "epsilon = Constant(0.1)\n", "r = Constant(\"1.0\")\n", "f = lambda x : exp(x)\n", "u_e = lambda x : exp(x)/(1-epsilon**2) + (exp((x-1)/epsilon)*(exp(1.0)-exp(-1.0/epsilon)) + \\\n", " exp(-x/epsilon)*(1 - exp((epsilon-1)/epsilon)))/(1-epsilon**2)*(exp(-2/epsilon) -1) " ] }, { "cell_type": "markdown", "id": "55f760f6", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Define the mesh\n", "\n", "We'll use a uniform mesh with $N=16$ intervals." ] }, { "cell_type": "code", "execution_count": null, "id": "fc83b797", "metadata": {}, "outputs": [], "source": [ "N = 16 \n", "mesh1D = UnitIntervalMesh(N)" ] }, { "cell_type": "markdown", "id": "acc5b172", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "It can be a good idea to check the mesh properties (especially in higher dimensions, or with hierarchically constructed meshes):" ] }, { "cell_type": "code", "execution_count": null, "id": "5365c6c4", "metadata": {}, "outputs": [], "source": [ "print(f\"mesh1D has {mesh1D.num_vertices()} vertices and {mesh1D.num_cells()} edges.\")" ] }, { "cell_type": "markdown", "id": "bf7a859e", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### The FEM space\n", "\n", "We define a function space on our mesh. The basic syntax is \\\n", "`V = FunctionSpace( , , )`\n", "\n", "In the following example we'll use Continuous Galerkin (\"CG\" = \"Lagrange\") elements of degree 1.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2003e980", "metadata": {}, "outputs": [], "source": [ "V = FunctionSpace(mesh1D, \"CG\", 1)" ] }, { "cell_type": "markdown", "id": "958422ff", "metadata": {}, "source": [ "While such spaces are used, primarily, for defining trial and test functions, and functions for storing solution, they can be used to represent any function you like (via interpolation), and preform operations on them (arithmetic, differential, integral, etc)." ] }, { "cell_type": "markdown", "id": "23ed05b5", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Here we'll define functions on the space `V` and use that to represent the true solution. To do that we'll\n", "* extract the spatial coordinates, `x`, which has only one member: `x[0]`\n", "* express the solution in terms of that." ] }, { "cell_type": "code", "execution_count": null, "id": "e0878b76", "metadata": {}, "outputs": [], "source": [ "u_true = Function(V)\n", "x = SpatialCoordinate(V)\n", "u_true.interpolate( u_e(x[0]) )\n", "plot(u_true)" ] }, { "cell_type": "markdown", "id": "5b9b3cc8", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### The Weak Form\n", "\n", "Define trial and test functions:" ] }, { "cell_type": "code", "execution_count": null, "id": "1eabe3b1", "metadata": {}, "outputs": [], "source": [ "u = TrialFunction(V)\n", "v = TestFunction(V)" ] }, { "cell_type": "markdown", "id": "8e3e94dc", "metadata": {}, "source": [ "We now define the bilinear and linear forms\n", "$$ a(u,v) := \\varepsilon^2(u',v') + (ru,v) \n", "\\qquad and \\quad \n", "l(v) := (f,v).\n", "$$\n", "* one can represent $u'(x)$ and $v'(x)$ as `u.dx(0)` and `v.dx(0)`, where here the index says we are differentiating with respect to the first (and only) dimension.\n", "* `dx` is also used to represents the (symbolic) integral operator." ] }, { "cell_type": "code", "execution_count": null, "id": "59d304ce", "metadata": {}, "outputs": [], "source": [ "a = (epsilon**2)*u.dx(0)*v.dx(0)*dx + r*u*v*dx\n", "l = f(x[0])*v*dx" ] }, { "cell_type": "markdown", "id": "20ffa5c9", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Boundary conditions\n", "\n", "It is easy to impose any reasonable boundary conditions you can think of. Here we'll use set homogeneous Dirichlet: " ] }, { "cell_type": "code", "execution_count": null, "id": "ff748ab2", "metadata": {}, "outputs": [], "source": [ "BCs = DirichletBC(V, 0.0, \"on_boundary\")" ] }, { "cell_type": "markdown", "id": "e3a70ca7", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Solve\n", "\n", "And we are ready to solve! Define a function to store the solution, and call the solver:" ] }, { "cell_type": "code", "execution_count": null, "id": "c4a05a6e", "metadata": {}, "outputs": [], "source": [ "uN = Function(V)\n", "solve(a == l, uN, bcs=BCs)\n", "plot(uN, label=\"$u^N$\")\n", "plt.legend()" ] }, { "cell_type": "markdown", "id": "4f581bea", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Compute errors\n", "\n", "We'll define a function $d(x)=u(x)-u^N(x)$ and compute the $L_2$, energy, and $\\max$ norms of that.\n", "* For $L_2$-type norms, apply the `dx` operator to form an integral, and then `assemble` that to evaluate.\n", "* For maximum pointwise error, choose a set of points (e.g., the mesh points) and evaluate `d` at those, using the `at` method. (There are other approaches)." ] }, { "cell_type": "code", "execution_count": null, "id": "224af739", "metadata": {}, "outputs": [], "source": [ "d = Function(V).interpolate(uN-u_true)\n", "\n", "L2_error = sqrt(assemble( d**2*dx)) # same as L2=\n", "H1_error = sqrt(assemble( d.dx(0)**2*dx))\n", "a_error = sqrt( float(epsilon)*H1_error**2.0 + L2_error**2.0)\n", "\n", "xp = np.linspace(0,1,N+1) # points at which we compute errors\n", "max_pw_error = max(np.abs(d.at(xp)))\n", "\n", "print(f\"Errors: L2={L2_error : 8.3e}, H1={H1_error : 8.3e}, Energy={a_error : 8.3e}, Max={max_pw_error: 8.3e}\")" ] }, { "cell_type": "markdown", "id": "1d5759cc", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Plot solutions\n", "\n", "While we can plot solutions using the `plot()` function, which is an interface to `matplotlib`. But if you are familiar with `matplotlib` it might be better to use it directly:" ] }, { "cell_type": "code", "execution_count": null, "id": "62f2d2b1", "metadata": {}, "outputs": [], "source": [ "x_vals = mesh1D.coordinates.dat.data[:]\n", "\n", "plot(u_true, label='$u(x)$') # true\n", "plt.plot(x_vals, uN.at(x_vals),'r--o', label='$u^N$') # use for small N\n", "plt.xlabel(\"x\")\n", "plt.legend()" ] }, { "cell_type": "markdown", "id": "6b4cab4d", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## A fitted mesh in 1D\n", "\n", "We now what to compute solutions (and corresponding errors) on a fitted mesh. At this stage, it makes sense to wrap our solver in a function, so we can call it for different $\\varepsilon$ and meshes." ] }, { "cell_type": "markdown", "id": "027c3861", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "### Solver as a function\n", "\n", "This function also returns both the FEM and true solution for the provided `epsilon`, on the same function space." ] }, { "cell_type": "code", "execution_count": null, "id": "373950a6", "metadata": {}, "outputs": [], "source": [ "def Solve_1DRD(epsilon, my_mesh, my_degree):\n", " V = FunctionSpace(my_mesh, \"CG\", my_degree)\n", " x = SpatialCoordinate(V)\n", " u_true = Function(V)\n", " u_true.interpolate( u_e(x[0]) )\n", " u = TrialFunction(V)\n", " v = TestFunction(V)\n", " a = (epsilon**2)*u.dx(0)*v.dx(0)*dx + r*u*v*dx\n", " l = f(x[0])*v*dx\n", " uN = Function(V)\n", " BCs = DirichletBC(V, 0.0, \"on_boundary\")\n", " solve(a == l, uN, bcs=BCs)\n", " return(uN, u_true)" ] }, { "cell_type": "markdown", "id": "d3b91831", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Let's check if that works. We'll solve for $\\varepsilon^2 = 10^{-4}$ and $N=16$; expect oscillations." ] }, { "cell_type": "code", "execution_count": null, "id": "330af1b4", "metadata": {}, "outputs": [], "source": [ "N = np.int32(16)\n", "epsilon = 1.0e-2\n", "mesh1D = UnitIntervalMesh(N)\n", "uN,uE = Solve_1DRD(epsilon, mesh1D, 1)\n", "x_vals = mesh1D.coordinates.dat.data[:]\n", "plot(uE, label='$u(x)$') # true\n", "plt.plot(x_vals, uN.at(x_vals),'r--o', label='$u^N$') # use for small N\n", "plt.legend()" ] }, { "cell_type": "markdown", "id": "ec4263c4", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### A 1D fitted mesh\n", "\n", "In this example we'll construct a simple mesh piecewise uniform (Shishkin) mesh, and solve on that. There are various ways to do this. In this example, we construct a uniform mesh, and then change the mesh values.\n", "\n", "First define the mesh points as a `numpy` array." ] }, { "cell_type": "code", "execution_count": null, "id": "de570ea6", "metadata": {}, "outputs": [], "source": [ "N=16\n", "p=2 # degree of the element\n", "tau = min( (p+1)*epsilon*np.log(float(N)), 0.25) # mesh transition point\n", "mesh1 = np.linspace(0.0, tau, int(N/4)+1) # left mesh; assuming N divisible by 4\n", "mesh2 = np.linspace(tau, 1-tau, int(N/2)+1) # interior\n", "mesh3 = np.linspace(1-tau, 1.0, int(N/4)+1) # right mesh\n", "\n", "Smesh_points = np.concatenate( (mesh1, mesh2[1:-1], mesh3) ) " ] }, { "cell_type": "markdown", "id": "db588600", "metadata": {}, "source": [ "Now make a (Firedrake) mesh, and overwrite the mesh points" ] }, { "cell_type": "code", "execution_count": null, "id": "4ed6c01e", "metadata": {}, "outputs": [], "source": [ "Smesh1D = UnitIntervalMesh(N)\n", "Smesh1D.coordinates.dat.data[:]=Smesh_points" ] }, { "cell_type": "markdown", "id": "850fab63", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Let's test:" ] }, { "cell_type": "code", "execution_count": null, "id": "9bdd0944", "metadata": {}, "outputs": [], "source": [ "uN,uE = Solve_1DRD(epsilon, Smesh1D, p)" ] }, { "cell_type": "code", "execution_count": null, "id": "5fef8a62", "metadata": {}, "outputs": [], "source": [ "plot(uE, label='$u(x)$') # true\n", "plt.plot(Smesh_points, uN.at(Smesh_points),'r--o', label='$u^N$') # use for small N\n", "plt.legend()" ] }, { "cell_type": "markdown", "id": "43219afa", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Of course, we'd like to check the errors. We'll define a function to do that:" ] }, { "cell_type": "code", "execution_count": null, "id": "f57b2e4e", "metadata": {}, "outputs": [], "source": [ "def ErrorNorms(u1, u2):\n", " d = Function(V).interpolate(u1-u2)\n", " L2 = sqrt(assemble( d**2*dx)) # same as L2=\n", " H1 = sqrt(assemble( d.dx(0)**2*dx))\n", " Energy = sqrt( float(epsilon)*H1**2.0 + L2**2.0)\n", " xp = np.linspace(0,1,N+1) # points at which we compute errors\n", " Max = max(np.abs(d.at(xp)))\n", " return (L2, H1, Energy, Max)" ] }, { "cell_type": "code", "execution_count": null, "id": "5553fab2", "metadata": {}, "outputs": [], "source": [ "L2, H1, Energy, Max = ErrorNorms(uN, uE)\n", "print(f\"epsilon={epsilon:0.1e}, N={N:4d}, Errors: L2={L2:.3e}, H1={H1:.3e}, Energy={Energy:.3e}, max={Max:.3e}\")" ] }, { "cell_type": "markdown", "id": "034a3108", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Convergence tables\n", "\n", "Many studies of (robust) methods for singularly perturbed problems feature convergence tables showing how the errors depend on both $N$ and $\\varepsilon$.\n", "\n", "Here is how we can do that (though it is not Firedrake-specific). \n", "\n", "First: a function for defining a mesh." ] }, { "cell_type": "code", "execution_count": null, "id": "960f4e21", "metadata": {}, "outputs": [], "source": [ "def ShishkinMesh(epsilon, N, p):\n", " tau = min( (p+1)*epsilon*np.log(float(N)), 0.25) # mesh transition point\n", " mesh1 = np.linspace(0.0, tau, int(N/4)+1) # left mesh; assuming N divisible by 4\n", " mesh2 = np.linspace(tau, 1-tau, int(N/2)+1) # interior\n", " mesh3 = np.linspace(1-tau, 1.0, int(N/4)+1) # right mesh\n", " Smesh_points = np.concatenate( (mesh1, mesh2[1:-1], mesh3) ) \n", " mesh1d = UnitIntervalMesh(N)\n", " mesh1d.coordinates.dat.data[:]=Smesh_points\n", " return(mesh1d)" ] }, { "cell_type": "markdown", "id": "02923374", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Now we'll solve the problem, and compute errors for \n", "* $\\varepsilon = 10^{-1}, 10^{-2}, \\dots, 10^{-6}$\n", "* $N = 32, 64, \\dots, 1024$" ] }, { "cell_type": "code", "execution_count": null, "id": "2551c62f", "metadata": {}, "outputs": [], "source": [ "p = 1\n", "Epsilons = 10.0**np.arange(-1,-7,-1)\n", "Ns = 2**np.arange(5,11)\n", "\n", "L2 = np.zeros( (len(Epsilons), len(Ns)) )\n", "Energy = np.zeros_like(L2)\n", "Max = np.zeros_like(L2)\n", "\n", "e=-1\n", "for epsilon in [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]:\n", " e+=1\n", " n=-1\n", " for N in Ns:\n", " n+=1\n", " Smesh1 = ShishkinMesh(epsilon, N, p)\n", " uN, uE = Solve_1DRD(epsilon, Smesh1, p)\n", " L2[e,n], H1, Energy[e,n], Max[e,n] = ErrorNorms(uN, uE)" ] }, { "cell_type": "markdown", "id": "7054f685", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "We'll display the table using a `pandas` data frame. You could also easily make a log-log plot." ] }, { "cell_type": "code", "execution_count": null, "id": "330ca46d", "metadata": {}, "outputs": [], "source": [ "pd.set_option('display.float_format', lambda x: '%10.2e' % x)\n", "table = pd.DataFrame(Max, index=Epsilons, columns=Ns)\n", "print('Maximum pointwise errors')\n", "print(table)" ] }, { "cell_type": "markdown", "id": "6e657c9f", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "The mandatory log-log plot:" ] }, { "cell_type": "code", "execution_count": null, "id": "e461ab94", "metadata": {}, "outputs": [], "source": [ "table.T.plot(loglog=True, style='--o')" ] }, { "cell_type": "markdown", "id": "aaa7b881", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## A two-dimensional problem\n", "\n", "We'll finish with the two-dimensional problem: _find $u=u(x,y)$ such that_\n", "$$-\\varepsilon^2 \\Delta u + r(x,y)u = f(x,y) \\text{ on } (0,1)^2,$$\n", "with homogeneous Dirichlet boundary conditions.\n", "\n", "The specific example we'll consider will have $r(x,y) = \\sqrt{1+x+y}$ \n", "and $f(x,y)=e^{x-y}$.\n", "\n" ] }, { "cell_type": "markdown", "id": "4c564d56", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### The mesh\n", "\n", "We'll solve with quadrilateral elements, on an $N \\times N$ tensor product (Shishkin) mesh. Handily, the `TensorRectangleMesh` method can take two numpy arrays as inputs, and generate the corresponding tensor product mesh." ] }, { "cell_type": "code", "execution_count": null, "id": "a7b03ddf", "metadata": {}, "outputs": [], "source": [ "N = 32\n", "p = 1 # degree of the element\n", "\n", "tau = min( (p+1)*epsilon*np.log(float(N)), 0.25) # mesh transition point\n", "mesh1 = np.linspace(0.0, tau, int(N/4)+1) # left mesh; assuming N divisible by 4\n", "mesh2 = np.linspace(tau, 1-tau, int(N/2)+1) # interior\n", "mesh3 = np.linspace(1-tau, 1.0, int(N/4)+1) # right mesh\n", "xp = np.concatenate( (mesh1, mesh2[1:-1], mesh3) )\n", "yp = xp\n", "\n", "mesh2D = TensorRectangleMesh(xp, yp, quadrilateral=True)\n", "print(f\"mesh2D has {mesh2D.num_vertices()} vertices and {mesh2D.num_cells()} edges.\")" ] }, { "cell_type": "markdown", "id": "a4715132", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### The function space and weak form.\n", "\n", "Since we are now in 2D, we can use the `grad` operator and `inner()` method. (We could have used these in 1D too, if we wished)." ] }, { "cell_type": "code", "execution_count": null, "id": "bfdf66f5", "metadata": {}, "outputs": [], "source": [ "V = FunctionSpace(mesh2D, \"CG\", p)\n", "x,y = SpatialCoordinate(V)\n", "u = TrialFunction(V)\n", "v = TestFunction(V)\n", "\n", "a = (epsilon**2)*inner(grad(u), grad(v))*dx + inner(sqrt(1+x+y)*u,v)*dx\n", "l = exp(x+y)*v*dx\n", "uN = Function(V)\n", "BCs = DirichletBC(V, 0.0, \"on_boundary\")\n", "solve(a == l, uN, bcs=BCs)" ] }, { "cell_type": "markdown", "id": "9b83b33e", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Plotting\n", "\n", "Finally we plot the solution. Plotting is a little limited (I think) and it is better to export for use in some dedicated tool, such as Paraview. But we'll make do:" ] }, { "cell_type": "code", "execution_count": null, "id": "c08caa82", "metadata": {}, "outputs": [], "source": [ "fig = plt.figure()\n", "axes = fig.add_subplot(projection='3d')\n", "trisurf(uN, axes=axes);" ] }, { "cell_type": "markdown", "id": "4b29cd34", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "To make an interactive plot, first make a 2D np array of solution values:" ] }, { "cell_type": "code", "execution_count": null, "id": "9de93c0e", "metadata": {}, "outputs": [], "source": [ "[X,Y] = np.meshgrid(xp,yp)\n", "uN_vals = np.zeros_like(X)\n", "for i in range(len(xp)):\n", " for j in range(len(yp)):\n", " uN_vals[i,j] = uN.at(xp[i],yp[j])" ] }, { "cell_type": "markdown", "id": "5e697ef6", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "The use `plot_surface` from matplotlib: (note: might not work on colab)" ] }, { "cell_type": "code", "execution_count": null, "id": "e11d6dcb", "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(8, 6))\n", "ax = fig.add_subplot(111, projection='3d')\n", "surf = ax.plot_surface(X,Y, uN_vals, cmap='bwr')\n", "fig.colorbar(surf)\n", "ax.set_xlabel('x', labelpad=20)\n", "ax.set_ylabel('y', labelpad=20)" ] }, { "cell_type": "markdown", "id": "05c764fb", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Other things Firedrake can and cannot do\n", "\n", "1. Many other elements are implemented; see https://www.firedrakeproject.org/variational-problems.html#supported-finite-elements. These include elemenst on triangles, tets, quads, and hex elements. \n", "2. dG elements are implemented, but only on conforming meshes.\n", "3. The 'irksome' module adds support for RK-methods in time (and solver support).\n", "4. Geometric multigrid \n", "5. Visualisation is limited (just an interface to `matplotlib`: use of Paraview recommended).\n", "6. Nonlinear problems are represented in a natural way, but care is needed with solver options.\n" ] }, { "cell_type": "markdown", "id": "5aba9e92", "metadata": {}, "source": [ "**The End**" ] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "FD25i", "language": "python", "name": "fd25i" }, "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.11.6" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": true, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }