diff --git a/.github/actions/install-dependencies/action.yml b/.github/actions/install-dependencies/action.yml index c49f3d98..c7db29e5 100644 --- a/.github/actions/install-dependencies/action.yml +++ b/.github/actions/install-dependencies/action.yml @@ -6,4 +6,4 @@ runs: - name: Install apt dependencies and upgrade pip shell: bash -el {0} run: | - apt-get update && apt-get install -y libxrender1 + apt-get update && apt-get install -y libxrender1 libgl1-mesa-dev mesa-utils gzip \ No newline at end of file diff --git a/.github/workflows/book_stable.yml b/.github/workflows/book_stable.yml index 576d30a1..b4e8ae5b 100644 --- a/.github/workflows/book_stable.yml +++ b/.github/workflows/book_stable.yml @@ -20,7 +20,7 @@ jobs: env: PYVISTA_OFF_SCREEN: false - PYVISTA_JUPYTER_BACKEND: "html" + PYVISTA_JUPYTER_BACKEND: html steps: - uses: actions/checkout@v5 @@ -31,10 +31,10 @@ jobs: - name: Install book deps run: | python3 -m pip install --break-system-packages -U pip setuptools pkgconfig poetry-core - python3 -m pip install --no-binary=h5py .[netgen] + python3 -m pip install --no-binary=h5py --no-build-isolation .[netgen] - name: Build the book - run: jupyter-book build . -W + run: jupyter-book build . - uses: actions/upload-artifact@v4 if: always() diff --git a/.github/workflows/test_stable.yml b/.github/workflows/test_stable.yml index b0be68bb..c9bff84f 100644 --- a/.github/workflows/test_stable.yml +++ b/.github/workflows/test_stable.yml @@ -5,19 +5,18 @@ on: workflow_call: pull_request: branches: ["release"] -env: - HDF5_MPI: "ON" - HDF5_DIR: "/usr/local/" - DISPLAY: ":99.0" - DEB_PYTHON_INSTALL_LAYOUT: deb_system - LIBGL_ALWAYS_SOFTWARE: 1 jobs: test: runs-on: ubuntu-latest container: ghcr.io/fenics/dolfinx/lab:stable env: + HDF5_MPI: "ON" + HDF5_DIR: "/usr/local/" + DEB_PYTHON_INSTALL_LAYOUT: deb_system PYVISTA_OFF_SCREEN: true + PYVISTA_JUPYTER_BACKEND: html + LIBGL_ALWAYS_SOFTWARE: 1 # Steps represent a sequence of tasks that will be executed as part of the job steps: @@ -30,7 +29,7 @@ jobs: - name: Install additional deps run: | python3 -m pip install -U pip setuptools pkgconfig poetry-core - python3 -m pip install --no-binary=h5py .[netgen] + python3 -m pip install --no-binary=h5py --no-build-isolation .[netgen] - name: Test complex notebooks in parallel working-directory: chapter1 diff --git a/chapter1/fundamentals_code.ipynb b/chapter1/fundamentals_code.ipynb index f06dc269..a8e7a586 100644 --- a/chapter1/fundamentals_code.ipynb +++ b/chapter1/fundamentals_code.ipynb @@ -107,8 +107,8 @@ "id": "3", "metadata": {}, "source": [ - "Note that in addition to give how many elements we would like to have in each direction.\n", - "We also have to supply the _MPI-communicator_.\n", + "Note that in addition to give how many elements we would like to have in each direction,\n", + "we also have to supply the _MPI-communicator_.\n", "This is to specify how we would like the program to behave in parallel.\n", "If we supply {py:data}`MPI.COMM_WORLD` we create a single mesh,\n", "whose data is distributed over the number of processors we would like to use.\n", @@ -124,11 +124,11 @@ " Once the mesh has been created, we can create the finite element function space $V$.\n", "The finite element function space does not need to be the same as the one used to describe the mesh.\n", "DOLFINx supports a wide range of arbitrary order finite element function spaces, see:\n", - "[Supported elements in DOLFINx](https://defelement.org/lists/implementations/basix.ufl.html)\n", + "[Supported elements in DOLFINx](https://defelement.org/lists/implementations/basix_ufl.html)\n", "for an extensive list.\n", "To create a function space, we need to specify what mesh the space is defined on,\n", - "what element famil the space is based on, and the degree of the element.\n", - "These can for instance be defned through a tuple `(\"family\", degree)`, as shown below" + "what element family the space is based on, and the degree of the element.\n", + "These can for instance be defined through a tuple `(\"family\", degree)`, as shown below" ] }, { @@ -248,7 +248,7 @@ "In FEniCSx, we do not specify boundary conditions as part of the function space,\n", "so it is sufficient to use a common space for the trial and test function.\n", "\n", - "We use the {py:mod}`Unified Form Language` (UFL) to specify the varitional formulations.\n", + "We use the {py:mod}`Unified Form Language` (UFL) to specify the variational formulations.\n", "See {cite}`fundamentals-ufl2014` for more details." ] }, @@ -296,7 +296,7 @@ "However, if we would like to change this parameter later in the simulation,\n", "we would have to redefine our variational formulation.\n", "The {py:attr}`dolfinx.fem.Constant.value` allows us to update the value in $f$ by using `f.value=5`.\n", - "Additionally, by indicating that $f$ is a constant, we speed of compilation of the variational\n", + "Additionally, by indicating that $f$ is a constant, we speed up compilation of the variational\n", "formulations required for the created linear system.\n", "```\n", "\n", diff --git a/chapter1/fundamentals_code.py b/chapter1/fundamentals_code.py index dffc9b32..4c8d0c20 100644 --- a/chapter1/fundamentals_code.py +++ b/chapter1/fundamentals_code.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.18.1 +# jupytext_version: 1.19.1 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -97,8 +97,8 @@ domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral) # - -# Note that in addition to give how many elements we would like to have in each direction. -# We also have to supply the _MPI-communicator_. +# Note that in addition to give how many elements we would like to have in each direction, +# we also have to supply the _MPI-communicator_. # This is to specify how we would like the program to behave in parallel. # If we supply {py:data}`MPI.COMM_WORLD` we create a single mesh, # whose data is distributed over the number of processors we would like to use. @@ -175,7 +175,7 @@ # In FEniCSx, we do not specify boundary conditions as part of the function space, # so it is sufficient to use a common space for the trial and test function. # -# We use the {py:mod}`Unified Form Language` (UFL) to specify the varitional formulations. +# We use the {py:mod}`Unified Form Language` (UFL) to specify the variational formulations. # See {cite}`fundamentals-ufl2014` for more details. @@ -200,7 +200,7 @@ # However, if we would like to change this parameter later in the simulation, # we would have to redefine our variational formulation. # The {py:attr}`dolfinx.fem.Constant.value` allows us to update the value in $f$ by using `f.value=5`. -# Additionally, by indicating that $f$ is a constant, we speed of compilation of the variational +# Additionally, by indicating that $f$ is a constant, we speed up compilation of the variational # formulations required for the created linear system. # ``` # diff --git a/chapter2/amr.py b/chapter2/amr.py index 836e6a51..2f2fd2eb 100644 --- a/chapter2/amr.py +++ b/chapter2/amr.py @@ -95,6 +95,7 @@ # Again, we visualize the curved mesh with pyvista. # + tags=["hide-input"] +plotter = pyvista.Plotter() curved_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(curved_mesh)) curved_grid.cell_data["ct"] = ct.values plotter.add_mesh( diff --git a/chapter2/hyperelasticity.ipynb b/chapter2/hyperelasticity.ipynb index 1271a4bb..818d3549 100644 --- a/chapter2/hyperelasticity.ipynb +++ b/chapter2/hyperelasticity.ipynb @@ -437,7 +437,7 @@ "for n in range(1, 10):\n", " T.value[2] = n * tval0\n", " problem.solve()\n", - " converged = problem.solver.getConvergedReason() > 0\n", + " converged = problem.solver.getConvergedReason()\n", " num_its = problem.solver.getIterationNumber()\n", " assert converged > 0, f\"Solver did not converge with reason {converged}.\"\n", "\n", diff --git a/chapter2/hyperelasticity.py b/chapter2/hyperelasticity.py index a088f592..b3fb1069 100644 --- a/chapter2/hyperelasticity.py +++ b/chapter2/hyperelasticity.py @@ -209,7 +209,7 @@ def right(x): for n in range(1, 10): T.value[2] = n * tval0 problem.solve() - converged = problem.solver.getConvergedReason() > 0 + converged = problem.solver.getConvergedReason() num_its = problem.solver.getIterationNumber() assert converged > 0, f"Solver did not converge with reason {converged}." diff --git a/chapter2/navierstokes.md b/chapter2/navierstokes.md index 336562cc..938115c3 100644 --- a/chapter2/navierstokes.md +++ b/chapter2/navierstokes.md @@ -68,11 +68,11 @@ As mentioned in the note in [Linear elasticity implementation](./linearelasticit $w_i=\sum_{j}\left(u_j\frac{\partial}{\partial x_j}\right)u_i = \sum_j u_j\frac{\partial u_i}{\partial x_j}$. This term can be implemented in FEniCSx as either `grad(u)*u`, since this expression becomes $\sum_j\frac{\partial u_i}{\partial x_j}u_j$, or as `dot(u, nabla_grad(u))` since this -expression becomes $\sum_i u_i\frac{\partial u_j}{x_i}$. We will use the notation `dot(u, nabla_grad(u))` below since it corresponds more closely to the standard notation $u\cdot \nabla u$. +expression becomes $\sum_i u_i\frac{\partial u_j}{\partial x_i}$. We will use the notation `dot(u, nabla_grad(u))` below since it corresponds more closely to the standard notation $u\cdot \nabla u$. ``` We now move on to the second step in our splitting scheme for the incompressible Navier-Stokes equations. In the first step, we computed the *tentative velocity* $u^*$ based on the pressure from the previous time step. -We may now use the computed tentative velocity to compute the new pressure $p^n$: +We may now use the computed tentative velocity to compute the new pressure $p^{n+1}$: ```{math} :label: ipcs-two \langle \nabla p^{n+1}, \nabla q \rangle = \langle \nabla p^n, \nabla q\rangle - \frac{\rho}{\Delta t}\langle \nabla \cdot u^*, q\rangle. diff --git a/chapter3/em.ipynb b/chapter3/em.ipynb index 4a983877..23a1998c 100644 --- a/chapter3/em.ipynb +++ b/chapter3/em.ipynb @@ -44,7 +44,7 @@ "\n", "which holds for an isotropic linear magnetic medium.\n", "Here, $\\mu$ is the magnetic permability of the material.\n", - "Now, since $B$ is solenodial (divergence free) accoording to Maxwell's equations, we known that $B$ must be the curl of some vector field $A$. This field is called the magnetic vector potential. Since the problem is static and thus $\\frac{\\partial D}{\\partial t}=0$, it follows that\n", + "Now, since $B$ is solenodial (divergence free) according to Maxwell's equations, we known that $B$ must be the curl of some vector field $A$. This field is called the magnetic vector potential. Since the problem is static and thus $\\frac{\\partial D}{\\partial t}=0$, it follows that\n", "\n", "$$\n", "J = \\nabla \\times H = \\nabla \\times(\\mu^{-1} B)=\\nabla \\times (\\mu^{-1}\\nabla \\times A ) = -\\nabla \\cdot (\\mu^{-1}\\nabla A).\n", diff --git a/chapter3/em.py b/chapter3/em.py index e8dd582f..6efe178a 100644 --- a/chapter3/em.py +++ b/chapter3/em.py @@ -53,7 +53,7 @@ # # which holds for an isotropic linear magnetic medium. # Here, $\mu$ is the magnetic permability of the material. -# Now, since $B$ is solenodial (divergence free) accoording to Maxwell's equations, we known that $B$ must be the curl of some vector field $A$. This field is called the magnetic vector potential. Since the problem is static and thus $\frac{\partial D}{\partial t}=0$, it follows that +# Now, since $B$ is solenodial (divergence free) according to Maxwell's equations, we known that $B$ must be the curl of some vector field $A$. This field is called the magnetic vector potential. Since the problem is static and thus $\frac{\partial D}{\partial t}=0$, it follows that # # $$ # J = \nabla \times H = \nabla \times(\mu^{-1} B)=\nabla \times (\mu^{-1}\nabla \times A ) = -\nabla \cdot (\mu^{-1}\nabla A). diff --git a/chapter3/subdomains.ipynb b/chapter3/subdomains.ipynb index c43de222..250e3490 100644 --- a/chapter3/subdomains.ipynb +++ b/chapter3/subdomains.ipynb @@ -5,11 +5,13 @@ "metadata": {}, "source": [ "# Defining subdomains for different materials\n", + "\n", "Author: Jørgen S. Dokken\n", "\n", "Solving PDEs in domains made up of different materials is a frequently encountered task. In FEniCSx, we handle these problems by defining a Discontinous cell-wise constant function.\n", "Such a function can be created over any mesh in the following way\n", - "## Subdomains on built-in meshes" + "\n", + "## Subdomains on built-in meshes\n" ] }, { @@ -56,7 +58,7 @@ }, "source": [ "We will use a simple example with two materials in two dimensions to demonstrate the idea. The whole domain will be $\\Omega=[0,1]\\times[0,1]$, which consists of two subdomains\n", - "$\\Omega_0=[0,1]\\times [0,1/2]$ and $\\Omega_1=[0,1]\\times[1/2, 1]$. We start by creating two python functions, where each returns `True` if the input coordinate is inside its domain." + "$\\Omega_0=[0,1]\\times [0,1/2]$ and $\\Omega_1=[0,1]\\times[1/2, 1]$. We start by creating two python functions, where each returns `True` if the input coordinate is inside its domain.\n" ] }, { @@ -84,14 +86,16 @@ "$$\n", "-\\nabla \\cdot [\\kappa (x,y)\\nabla u(x, y)]= 1 \\qquad \\text{in } \\Omega,\n", "$$\n", + "\n", "$$\n", "u=u_D=1 \\qquad \\text{on } \\partial\\Omega_D=[0,y], y\\in[0,1]\n", "$$\n", + "\n", "$$\n", "-\\frac{\\partial u}{\\partial n}=0 \\qquad \\text{on } \\partial\\Omega\\setminus \\partial\\Omega_D\n", "$$\n", "\n", - "Our next step is to define $\\kappa$" + "Our next step is to define $\\kappa$\n" ] }, { @@ -112,7 +116,7 @@ "In the previous code block, we found which cells (triangular elements) satisfy the condition for being in $\\Omega_0, \\Omega_1$. As the $DG-0$ function contains only one degree of freedom per cell, there is a one to one mapping between the cell indicies and the degrees of freedom. We let $\\kappa=\\begin{cases}\n", "1 &\\text{if } x\\in\\Omega_0\\\\\n", "0.1& \\text{if } x\\in\\Omega_1\\\\\n", - "\\end{cases}$" + "\\end{cases}$\n" ] }, { @@ -129,7 +133,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We are now ready to define our variational formulation and Dirichlet boundary condition after using integration by parts" + "We are now ready to define our variational formulation and Dirichlet boundary condition after using integration by parts\n" ] }, { @@ -151,7 +155,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We can now solve and visualize the solution of the problem" + "We can now solve and visualize the solution of the problem\n" ] }, { @@ -217,7 +221,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We clearly observe different behavior in the two regions, which both have the same Dirichlet boundary condition on the left side, where $x=0$." + "We clearly observe different behavior in the two regions, which both have the same Dirichlet boundary condition on the left side, where $x=0$.\n" ] }, { @@ -227,7 +231,8 @@ }, "source": [ "## Interpolation with Python-function\n", - "As we saw in the first approach, in many cases, we can use the geometrical coordinates to determine which coefficient we should use. Using the unstructured mesh from the previous example, we illustrate an alternative approach using interpolation:" + "\n", + "As we saw in the first approach, in many cases, we can use the geometrical coordinates to determine which coefficient we should use. Using the unstructured mesh from the previous example, we illustrate an alternative approach using interpolation:\n" ] }, { @@ -261,7 +266,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We verify this by assembling the error between this new function and the old one" + "We verify this by assembling the error between this new function and the old one\n" ] }, { @@ -280,7 +285,8 @@ "metadata": {}, "source": [ "## Subdomains defined from external mesh data\n", - "Let us now consider the same problem, but using GMSH to generate the mesh and subdomains. We will then in turn show how to use this data to generate discontinuous functions in DOLFINx." + "\n", + "Let us now consider the same problem, but using GMSH to generate the mesh and subdomains. We will then in turn show how to use this data to generate discontinuous functions in DOLFINx.\n" ] }, { @@ -329,7 +335,8 @@ "metadata": {}, "source": [ "## Read in MSH files with DOLFINx\n", - "You can read in MSH files with DOLFINx, which will read them in on a single process, and then distribute them over the available ranks in the MPI communicator." + "\n", + "You can read in MSH files with DOLFINx, which will read them in on a single process, and then distribute them over the available ranks in the MPI communicator.\n" ] }, { @@ -353,14 +360,17 @@ }, "source": [ "## Convert msh-files to XDMF using meshio\n", + "\n", "We will use `meshio` to read in the `msh` file, and convert it to a more suitable IO format. Meshio requires `h5py`, and can be installed on linux with the following commands:\n", + "\n", "```{code}\n", "export HDF5_MPI=\"ON\"\n", "export CC=mpicc\n", "export HDF5_DIR=\"/usr/lib/x86_64-linux-gnu/hdf5/mpich/\"\n", "pip3 install --no-cache-dir --no-binary=h5py h5py meshio\n", "```\n", - "We start by creating a convenience function for extracting data for a single cell type, and creating a new `meshio.Mesh`." + "\n", + "We start by creating a convenience function for extracting data for a single cell type, and creating a new `meshio.Mesh`.\n" ] }, { @@ -385,7 +395,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "This function returns a meshio mesh, including physical markers for the given type. The `prune_z` argument is for cases where we want to use two dimensional meshes. The last coordinate in the mesh (as it is generated in a 3D space) has to be removed for DOLFINx to consider this as a two dimensional geometry." + "This function returns a meshio mesh, including physical markers for the given type. The `prune_z` argument is for cases where we want to use two dimensional meshes. The last coordinate in the mesh (as it is generated in a 3D space) has to be removed for DOLFINx to consider this as a two dimensional geometry.\n" ] }, { @@ -401,8 +411,8 @@ " # Create and save one file for the mesh, and one file for the facets\n", " triangle_mesh = create_mesh(msh, \"triangle\", prune_z=True)\n", " line_mesh = create_mesh(msh, \"line\", prune_z=True)\n", - " meshio.write(\"mesh.xdmf\", triangle_mesh)\n", - " meshio.write(\"mt.xdmf\", line_mesh)\n", + " meshio.write(\"mesh.xdmf\", triangle_mesh, compression=None)\n", + " meshio.write(\"mt.xdmf\", line_mesh, compression=None)\n", "MPI.COMM_WORLD.barrier()" ] }, @@ -412,7 +422,7 @@ "source": [ "We have now written the mesh and the cell markers to one file, and the facet markers in a separate file. We can now read this data in DOLFINx using `XDMFFile.read_mesh` and `XDMFFile.read_meshtags`. The `dolfinx.MeshTags` stores the index of the entity, along with the value of the marker in two one dimensional arrays.\n", "\n", - "Note that we have generated and written the mesh on only one processor. However, the `xdmf`-format supports parallel IO, and we can thus read the mesh in parallel." + "Note that we have generated and written the mesh on only one processor. However, the `xdmf`-format supports parallel IO, and we can thus read the mesh in parallel.\n" ] }, { @@ -433,7 +443,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We have now read in the mesh and corresponding cell and facet data. We can now create our discontinuous function `kappa` as follows" + "We have now read in the mesh and corresponding cell and facet data. We can now create our discontinuous function `kappa` as follows\n" ] }, { @@ -454,7 +464,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We can also efficiently use the facet data `ft` to create the Dirichlet boundary condition" + "We can also efficiently use the facet data `ft` to create the Dirichlet boundary condition\n" ] }, { @@ -475,7 +485,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We can now solve the problem in a similar fashion as above" + "We can now solve the problem in a similar fashion as above\n" ] }, { diff --git a/chapter3/subdomains.py b/chapter3/subdomains.py index 33baff2e..3f368778 100644 --- a/chapter3/subdomains.py +++ b/chapter3/subdomains.py @@ -14,11 +14,14 @@ # --- # # Defining subdomains for different materials +# # Author: Jørgen S. Dokken # # Solving PDEs in domains made up of different materials is a frequently encountered task. In FEniCSx, we handle these problems by defining a Discontinous cell-wise constant function. # Such a function can be created over any mesh in the following way +# # ## Subdomains on built-in meshes +# # + from dolfinx import default_scalar_type @@ -54,6 +57,7 @@ # We will use a simple example with two materials in two dimensions to demonstrate the idea. The whole domain will be $\Omega=[0,1]\times[0,1]$, which consists of two subdomains # $\Omega_0=[0,1]\times [0,1/2]$ and $\Omega_1=[0,1]\times[1/2, 1]$. We start by creating two python functions, where each returns `True` if the input coordinate is inside its domain. +# # + @@ -74,14 +78,17 @@ def Omega_1(x): # $$ # -\nabla \cdot [\kappa (x,y)\nabla u(x, y)]= 1 \qquad \text{in } \Omega, # $$ +# # $$ # u=u_D=1 \qquad \text{on } \partial\Omega_D=[0,y], y\in[0,1] # $$ +# # $$ # -\frac{\partial u}{\partial n}=0 \qquad \text{on } \partial\Omega\setminus \partial\Omega_D # $$ # # Our next step is to define $\kappa$ +# kappa = Function(Q) cells_0 = locate_entities(mesh, mesh.topology.dim, Omega_0) @@ -91,11 +98,13 @@ def Omega_1(x): # 1 &\text{if } x\in\Omega_0\\ # 0.1& \text{if } x\in\Omega_1\\ # \end{cases}$ +# kappa.x.array[cells_0] = np.full_like(cells_0, 1, dtype=default_scalar_type) kappa.x.array[cells_1] = np.full_like(cells_1, 0.1, dtype=default_scalar_type) -# We are now ready to define our variational formulation and Dirichlet boundary condition after using integration by parts +# We are now ready to define our variational formulation and Dirichlet boundary condition after using integration by parts +# V = functionspace(mesh, ("Lagrange", 1)) u, v = TrialFunction(V), TestFunction(V) @@ -106,6 +115,7 @@ def Omega_1(x): bcs = [dirichletbc(default_scalar_type(1), dofs, V)] # We can now solve and visualize the solution of the problem +# # + problem = LinearProblem( @@ -152,9 +162,12 @@ def Omega_1(x): # We clearly observe different behavior in the two regions, which both have the same Dirichlet boundary condition on the left side, where $x=0$. +# # ## Interpolation with Python-function +# # As we saw in the first approach, in many cases, we can use the geometrical coordinates to determine which coefficient we should use. Using the unstructured mesh from the previous example, we illustrate an alternative approach using interpolation: +# def eval_kappa(x): @@ -172,13 +185,16 @@ def eval_kappa(x): kappa2.interpolate(eval_kappa) # We verify this by assembling the error between this new function and the old one +# # Difference in kappa's error = mesh.comm.allreduce(assemble_scalar(form((kappa - kappa2) ** 2 * dx))) print(error) # ## Subdomains defined from external mesh data +# # Let us now consider the same problem, but using GMSH to generate the mesh and subdomains. We will then in turn show how to use this data to generate discontinuous functions in DOLFINx. +# gmsh.initialize() proc = MPI.COMM_WORLD.rank @@ -215,7 +231,9 @@ def eval_kappa(x): gmsh.finalize() # ## Read in MSH files with DOLFINx +# # You can read in MSH files with DOLFINx, which will read them in on a single process, and then distribute them over the available ranks in the MPI communicator. +# mesh_data = gmshio.read_from_msh("mesh.msh", MPI.COMM_WORLD, gdim=2) mesh = mesh_data.mesh @@ -225,14 +243,18 @@ def eval_kappa(x): facet_markers = mesh_data.facet_tags # ## Convert msh-files to XDMF using meshio +# # We will use `meshio` to read in the `msh` file, and convert it to a more suitable IO format. Meshio requires `h5py`, and can be installed on linux with the following commands: +# # ```{code} # export HDF5_MPI="ON" # export CC=mpicc # export HDF5_DIR="/usr/lib/x86_64-linux-gnu/hdf5/mpich/" # pip3 install --no-cache-dir --no-binary=h5py h5py meshio # ``` +# # We start by creating a convenience function for extracting data for a single cell type, and creating a new `meshio.Mesh`. +# def create_mesh(mesh, cell_type, prune_z=False): @@ -248,6 +270,7 @@ def create_mesh(mesh, cell_type, prune_z=False): # This function returns a meshio mesh, including physical markers for the given type. The `prune_z` argument is for cases where we want to use two dimensional meshes. The last coordinate in the mesh (as it is generated in a 3D space) has to be removed for DOLFINx to consider this as a two dimensional geometry. +# if proc == 0: # Read in mesh @@ -256,13 +279,14 @@ def create_mesh(mesh, cell_type, prune_z=False): # Create and save one file for the mesh, and one file for the facets triangle_mesh = create_mesh(msh, "triangle", prune_z=True) line_mesh = create_mesh(msh, "line", prune_z=True) - meshio.write("mesh.xdmf", triangle_mesh) - meshio.write("mt.xdmf", line_mesh) + meshio.write("mesh.xdmf", triangle_mesh, compression=None) + meshio.write("mt.xdmf", line_mesh, compression=None) MPI.COMM_WORLD.barrier() # We have now written the mesh and the cell markers to one file, and the facet markers in a separate file. We can now read this data in DOLFINx using `XDMFFile.read_mesh` and `XDMFFile.read_meshtags`. The `dolfinx.MeshTags` stores the index of the entity, along with the value of the marker in two one dimensional arrays. # # Note that we have generated and written the mesh on only one processor. However, the `xdmf`-format supports parallel IO, and we can thus read the mesh in parallel. +# with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf: mesh = xdmf.read_mesh(name="Grid") @@ -272,6 +296,7 @@ def create_mesh(mesh, cell_type, prune_z=False): ft = xdmf.read_meshtags(mesh, name="Grid") # We have now read in the mesh and corresponding cell and facet data. We can now create our discontinuous function `kappa` as follows +# Q = functionspace(mesh, ("DG", 0)) kappa = Function(Q) @@ -281,6 +306,7 @@ def create_mesh(mesh, cell_type, prune_z=False): kappa.x.array[top_cells] = np.full_like(top_cells, 0.1, dtype=default_scalar_type) # We can also efficiently use the facet data `ft` to create the Dirichlet boundary condition +# V = functionspace(mesh, ("Lagrange", 1)) u_bc = Function(V) @@ -290,6 +316,7 @@ def create_mesh(mesh, cell_type, prune_z=False): bcs = [dirichletbc(default_scalar_type(1), left_dofs, V)] # We can now solve the problem in a similar fashion as above +# # + u, v = TrialFunction(V), TestFunction(V) diff --git a/chapter4/compiler_parameters.ipynb b/chapter4/compiler_parameters.ipynb index 39feeb91..760f3437 100644 --- a/chapter4/compiler_parameters.ipynb +++ b/chapter4/compiler_parameters.ipynb @@ -82,7 +82,8 @@ "id": "4", "metadata": {}, "source": [ - "We start by considering the different levels of optimization that the C compiler can use on the optimized code. A list of optimization options and explanations can be found [here](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html)" + "We start by considering the different levels of optimization that the C compiler can use on the optimized code.\n", + "A list of optimization options and explanations can be found [here](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html)" ] }, { @@ -100,7 +101,9 @@ "id": "6", "metadata": {}, "source": [ - "The next option we can choose is if we want to compile the code with `-march=native` or not. This option enables instructions for the local machine, and can give different results on different systems. More information can be found [here](https://gcc.gnu.org/onlinedocs/gcc/AArch64-Options.html#g_t-march-and--mcpu-Feature-Modifiers)" + "The next option we can choose is if we want to compile the code with `-march=native` or not.\n", + "This option enables instructions for the local machine, and can give different results on different systems.\n", + "More information can be found [here](https://gcc.gnu.org/onlinedocs/gcc/AArch64-Options.html#g_t-march-and--mcpu-Feature-Modifiers)" ] }, { diff --git a/chapter4/compiler_parameters.py b/chapter4/compiler_parameters.py index 2ba44ff5..97e5c33c 100644 --- a/chapter4/compiler_parameters.py +++ b/chapter4/compiler_parameters.py @@ -64,11 +64,14 @@ def compile_form(space: str, degree: int, jit_options: Dict): return end - start -# We start by considering the different levels of optimization that the C compiler can use on the optimized code. A list of optimization options and explanations can be found [here](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html) +# We start by considering the different levels of optimization that the C compiler can use on the optimized code. +# A list of optimization options and explanations can be found [here](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html) optimization_options = ["-O1", "-O2", "-O3", "-Ofast"] -# The next option we can choose is if we want to compile the code with `-march=native` or not. This option enables instructions for the local machine, and can give different results on different systems. More information can be found [here](https://gcc.gnu.org/onlinedocs/gcc/AArch64-Options.html#g_t-march-and--mcpu-Feature-Modifiers) +# The next option we can choose is if we want to compile the code with `-march=native` or not. +# This option enables instructions for the local machine, and can give different results on different systems. +# More information can be found [here](https://gcc.gnu.org/onlinedocs/gcc/AArch64-Options.html#g_t-march-and--mcpu-Feature-Modifiers) march_native = [True, False] diff --git a/chapter4/newton-solver.ipynb b/chapter4/newton-solver.ipynb index f979e173..750832f7 100644 --- a/chapter4/newton-solver.ipynb +++ b/chapter4/newton-solver.ipynb @@ -184,11 +184,24 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "0d5ce361", "metadata": {}, "outputs": [], "source": [ "solver = PETSc.KSP().create(mesh.comm)\n", + "solver.setType(\"preonly\")\n", + "solver.getPC().setType(\"lu\")\n", + "solver.getPC().setFactorSolverType(\"mumps\")\n", + "solver.setErrorIfNotConverged(True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ "solver.setOperators(A)\n", "du = dolfinx.fem.Function(V)" ] @@ -257,7 +270,7 @@ "\n", " # Compute norm of update\n", " correction_norm = du.x.petsc_vec.norm(0)\n", - " print(f\"Iteration {i}: Correction norm {correction_norm}\")\n", + " PETSc.Sys.Print(f\"Iteration {i}: Correction norm {correction_norm}\")\n", " if correction_norm < 1e-10:\n", " break\n", " solutions[i, :] = uh.x.array[sort_order]" @@ -279,7 +292,7 @@ "outputs": [], "source": [ "dolfinx.fem.petsc.assemble_vector(L, residual)\n", - "print(f\"Final residual {L.norm(0)}\")\n", + "PETSc.Sys.Print(f\"Final residual {L.norm(0)}\")\n", "A.destroy()\n", "L.destroy()\n", "solver.destroy()" @@ -314,7 +327,7 @@ " u_ex = root(x)\n", " L2_error = dolfinx.fem.form(ufl.inner(uh - u_ex, uh - u_ex) * ufl.dx)\n", " global_L2 = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_error), op=MPI.SUM)\n", - " print(f\"L2-error (root {j}) {np.sqrt(global_L2)}\")\n", + " PETSc.Sys.Print(f\"L2-error (root {j}) {np.sqrt(global_L2)}\")\n", "\n", " kwargs = {} if j == 0 else {\"label\": \"u_exact\"}\n", " plt.plot(x_spacing, root(x_spacing.reshape(1, -1)), *args, **kwargs)\n", @@ -410,7 +423,11 @@ "A = dolfinx.fem.petsc.create_matrix(jacobian)\n", "L = dolfinx.fem.petsc.create_vector(dolfinx.fem.extract_function_spaces(residual))\n", "solver = PETSc.KSP().create(mesh.comm)\n", - "solver.setOperators(A)" + "solver.setOperators(A)\n", + "solver.setType(\"preonly\")\n", + "solver.getPC().setType(\"lu\")\n", + "solver.getPC().setFactorSolverType(\"mumps\")\n", + "solver.setErrorIfNotConverged(True)" ] }, { @@ -461,15 +478,20 @@ " dolfinx.fem.petsc.assemble_matrix(A, jacobian, bcs=[bc])\n", " A.assemble()\n", " dolfinx.fem.petsc.assemble_vector(L, residual)\n", + "\n", + " # Compute b - alpha * J(u_D-u_(i-1))\n", + " dolfinx.fem.petsc.apply_lifting(\n", + " L, [jacobian], [[bc]], x0=[uh.x.petsc_vec], alpha=-1.0\n", + " )\n", " L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)\n", - " L.scale(-1)\n", "\n", - " # Compute b - J(u_D-u_(i-1))\n", - " dolfinx.fem.petsc.apply_lifting(L, [jacobian], [[bc]], x0=[uh.x.petsc_vec], alpha=1)\n", - " # Set du|_bc = u_{i-1}-u_D\n", - " dolfinx.fem.petsc.set_bc(L, [bc], uh.x.petsc_vec, 1.0)\n", + " # Set du|_bc = - (u_{i-1}-u_D)\n", + " dolfinx.fem.petsc.set_bc(L, [bc], uh.x.petsc_vec, -1.0)\n", " L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)\n", "\n", + " # Compute negative residual\n", + " L.scale(-1)\n", + "\n", " # Solve linear problem\n", " solver.solve(L, du.x.petsc_vec)\n", " du.x.scatter_forward()\n", @@ -486,8 +508,10 @@ " np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(error), op=MPI.SUM))\n", " )\n", " du_norm.append(correction_norm)\n", - "\n", - " print(f\"Iteration {i}: Correction norm {correction_norm}, L2 error: {L2_error[-1]}\")\n", + " PETSc.Sys.Print(\n", + " f\"Iteration {i}: Correction norm {correction_norm}, L2 error: {L2_error[-1]}\",\n", + " flush=True,\n", + " )\n", " if correction_norm < 1e-10:\n", " break" ] @@ -542,8 +566,7 @@ "outputs": [], "source": [ "error_max = domain.comm.allreduce(np.max(np.abs(uh.x.array - u_D.x.array)), op=MPI.MAX)\n", - "if domain.comm.rank == 0:\n", - " print(f\"Error_max: {error_max:.2e}\")" + "PETSc.Sys.Print(f\"Error_max: {error_max:.2e}\")" ] }, { diff --git a/chapter4/newton-solver.py b/chapter4/newton-solver.py index 882ddae9..80b460af 100644 --- a/chapter4/newton-solver.py +++ b/chapter4/newton-solver.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.18.1 +# jupytext_version: 1.19.1 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -99,6 +99,11 @@ def root_1(x): # Next, we create the linear solver and the vector to hold `du`. solver = PETSc.KSP().create(mesh.comm) +solver.setType("preonly") +solver.getPC().setType("lu") +solver.getPC().setFactorSolverType("mumps") +solver.setErrorIfNotConverged(True) + solver.setOperators(A) du = dolfinx.fem.Function(V) @@ -139,7 +144,7 @@ def root_1(x): # Compute norm of update correction_norm = du.x.petsc_vec.norm(0) - print(f"Iteration {i}: Correction norm {correction_norm}") + PETSc.Sys.Print(f"Iteration {i}: Correction norm {correction_norm}") if correction_norm < 1e-10: break solutions[i, :] = uh.x.array[sort_order] @@ -147,7 +152,7 @@ def root_1(x): # We now compute the magnitude of the residual. dolfinx.fem.petsc.assemble_vector(L, residual) -print(f"Final residual {L.norm(0)}") +PETSc.Sys.Print(f"Final residual {L.norm(0)}") A.destroy() L.destroy() solver.destroy() @@ -167,7 +172,7 @@ def root_1(x): u_ex = root(x) L2_error = dolfinx.fem.form(ufl.inner(uh - u_ex, uh - u_ex) * ufl.dx) global_L2 = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_error), op=MPI.SUM) - print(f"L2-error (root {j}) {np.sqrt(global_L2)}") + PETSc.Sys.Print(f"L2-error (root {j}) {np.sqrt(global_L2)}") kwargs = {} if j == 0 else {"label": "u_exact"} plt.plot(x_spacing, root(x_spacing.reshape(1, -1)), *args, **kwargs) @@ -227,6 +232,10 @@ def u_exact(x): L = dolfinx.fem.petsc.create_vector(dolfinx.fem.extract_function_spaces(residual)) solver = PETSc.KSP().create(mesh.comm) solver.setOperators(A) +solver.setType("preonly") +solver.getPC().setType("lu") +solver.getPC().setFactorSolverType("mumps") +solver.setErrorIfNotConverged(True) # Since this problem has strong Dirichlet conditions, we need to apply lifting to the right hand side of our Newton problem. # We previously had that we wanted to solve the system: @@ -262,15 +271,20 @@ def u_exact(x): dolfinx.fem.petsc.assemble_matrix(A, jacobian, bcs=[bc]) A.assemble() dolfinx.fem.petsc.assemble_vector(L, residual) + + # Compute b - alpha * J(u_D-u_(i-1)) + dolfinx.fem.petsc.apply_lifting( + L, [jacobian], [[bc]], x0=[uh.x.petsc_vec], alpha=-1.0 + ) L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) - L.scale(-1) - # Compute b - J(u_D-u_(i-1)) - dolfinx.fem.petsc.apply_lifting(L, [jacobian], [[bc]], x0=[uh.x.petsc_vec], alpha=1) - # Set du|_bc = u_{i-1}-u_D - dolfinx.fem.petsc.set_bc(L, [bc], uh.x.petsc_vec, 1.0) + # Set du|_bc = - (u_{i-1}-u_D) + dolfinx.fem.petsc.set_bc(L, [bc], uh.x.petsc_vec, -1.0) L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD) + # Compute negative residual + L.scale(-1) + # Solve linear problem solver.solve(L, du.x.petsc_vec) du.x.scatter_forward() @@ -287,8 +301,10 @@ def u_exact(x): np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(error), op=MPI.SUM)) ) du_norm.append(correction_norm) - - print(f"Iteration {i}: Correction norm {correction_norm}, L2 error: {L2_error[-1]}") + PETSc.Sys.Print( + f"Iteration {i}: Correction norm {correction_norm}, L2 error: {L2_error[-1]}", + flush=True, + ) if correction_norm < 1e-10: break @@ -315,8 +331,7 @@ def u_exact(x): # We compute the max error and plot the solution error_max = domain.comm.allreduce(np.max(np.abs(uh.x.array - u_D.x.array)), op=MPI.MAX) -if domain.comm.rank == 0: - print(f"Error_max: {error_max:.2e}") +PETSc.Sys.Print(f"Error_max: {error_max:.2e}") u_topology, u_cell_types, u_geometry = dolfinx.plot.vtk_mesh(V) u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry) diff --git a/docker/Dockerfile b/docker/Dockerfile index ff9f0c2d..a6b6b545 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -39,9 +39,13 @@ RUN python3 -m pip install vtk ADD pyproject.toml /tmp/pyproject.toml WORKDIR /tmp # As we install netgen from source we don't need the netgen deps here -RUN python3 -m pip install --no-cache-dir --no-binary=h5py -v . +RUN python3 -m pip install --no-cache-dir --no-binary=h5py --no-build-isolation -v . RUN python3 -m pip cache purge -RUN python3 -m pip install ngsPETSc --no-deps +RUN python3 -m pip install ngsPETSc --no-build-isolation --no-deps + +ENV PYVISTA_OFF_SCREEN="false" +ENV PYVISTA_JUPYTER_BACKEND="trame" +ENV LIBGL_ALWAYS_SOFTWARE=1 ENV PYVISTA_OFF_SCREEN="false" ENV PYVISTA_JUPYTER_BACKEND="trame" diff --git a/pyproject.toml b/pyproject.toml index 399e85f4..ee1104ec 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,7 +6,7 @@ build-backend = "setuptools.build_meta" name = "DOLFINx_Tutorial" version = "0.10.0" dependencies = [ - "jupyter-book", + "jupyter-book<2.0", "meshio", "h5py", "seaborn", diff --git a/references.bib b/references.bib index 476f5a61..eb8f2718 100644 --- a/references.bib +++ b/references.bib @@ -1,7 +1,7 @@ @book{Langtangen_Mardal_FEM_2019, author = {Langtangen, Hans Petter and Mardal, Kent-Andre}, - title = {Introduction to Numerical Methods for Variational Problems}, + title = {{Introduction to Numerical Methods for Variational Problems}}, year = {2019}, publisher = {Springer International Publishing}, address = {Cham}, @@ -11,7 +11,7 @@ @book{Langtangen_Mardal_FEM_2019 @article{ufl2014, author = {Aln\ae{}s, Martin S. and Logg, Anders and \O{}lgaard, Kristian B. and Rognes, Marie E. and Wells, Garth N.}, - title = {Unified Form Language: A Domain-Specific Language for Weak Formulations of Partial Differential Equations}, + title = {{Unified Form Language: A Domain-Specific Language for Weak Formulations of Partial Differential Equations}}, year = {2014}, issue_date = {February 2014}, publisher = {Association for Computing Machinery}, @@ -95,7 +95,7 @@ @article{Guermond1999 @book{Langtangen2016, author = {Langtangen, Hans Petter}, - title = {A Primer on Scientific Programming with Python}, + title = {{A Primer on Scientific Programming with Python}}, year = {2016}, publisher = {Springer Berlin Heidelberg}, address = {Berlin, Heidelberg}, @@ -106,7 +106,7 @@ @book{Langtangen2016 @book{FenicsTutorial, author = {Langtangen, Hans Petter and Logg, Anders}, - title = {Solving PDEs in Python: The FEniCS Tutorial I}, + title = {{Solving PDEs in Python: The FEniCS Tutorial I}}, year = {2016}, publisher = {Springer International Publishing}, address = {Cham}, @@ -144,11 +144,11 @@ @article{Nitsche1971 @misc{rehor2025pctools, - title = {{FEniCSx-pctools: Tools for PETSc block linear algebra preconditioning in FEniCSx}}, - author = {Martin Řehoř and Jack S. Hale}, - year = {2025}, - eprint = {2402.02523}, - archiveprefix = {arXiv}, - primaryclass = {cs.MS}, - doi = {10.48550/arXiv.2402.02523} -} \ No newline at end of file + author = {Řehoř, Martin and Hale, Jack S.}, + doi = {10.5334/jors.494}, + journal = {Journal of Open Research Software}, + keyword = {en}, + month = {Sep}, + title = {FEniCSx-pctools: Tools for PETSc Block Linear Algebra Preconditioning in FEniCSx}, + year = {2025} +}