diff --git a/.github/workflows/check_tutorials.yml b/.github/workflows/check_tutorials.yml
index f8b3ea3a28..be0e7e6477 100644
--- a/.github/workflows/check_tutorials.yml
+++ b/.github/workflows/check_tutorials.yml
@@ -37,10 +37,10 @@ jobs:
     # Steps represent a sequence of tasks that will be executed as part of the job
     steps:
     # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     
     - name: Setup Python ${{ matrix.python-version}}
-      uses: actions/setup-python@v3
+      uses: actions/setup-python@v5
       with:
         # Version range or exact version of a Python version to use, using SemVer's version range syntax.
         python-version: ${{ matrix.python-version}}
@@ -50,14 +50,11 @@ jobs:
     # Cache pip dependencies
     # Source: https://github.com/actions/cache/blob/main/examples.md#python---pip
     - name: Cache dependencies
-      uses: actions/cache@v3
+      uses: actions/cache@v4
       id: cache
       with:
-        # Cache the full python environment, this is more efficient than just caching pip
-        # https://blog.allenai.org/python-caching-in-github-actions-e9452698e98d
         path: ${{ env.pythonLocation }}
-        # Hash both ordinary requirements and those specific for developmnet.
-        key: ${{ env.pythonLocation }}-${{ hashFiles('setup.py') }}-${{ hashFiles('**/requirements.txt', '**/requirements-dev.txt') }}
+        key: ${{ env.pythonLocation }}-${{ hashFiles('pyproject.toml') }}
 
     - name: Install external libraries
       run: |
@@ -65,22 +62,19 @@ jobs:
         sudo apt-get install libglu1-mesa libgeos-dev libffi-dev
         export LD_LIBRARY_PATH=/usr/local/lib64/:$LD_LIBRARY_PATH
 
-    - name: Install requirements
+    # Install the requirements, including the extras needed for testing, and PorePy
+    # itself.
+    - name: Install requirements and PorePy
       run: |
         pip install -U pip
-        # Install requirements, using eager updates to avoid stalled dependencies due to caching
-        # https://blog.allenai.org/python-caching-in-github-actions-e9452698e98d
-        pip install --upgrade --upgrade-strategy eager -r requirements.txt
+        pip install .[testing]
+        pip freeze
 
+    # Install jupyter - needed to run the tutorial notebooks.
     - name: Install jupyter
       run: |
         pip install jupyter 
 
-    - name: Install PorePy
-      run: |
-        pip install .
-        pip freeze       
-
     - name: Check tutorials
       if: ${{always()}}
       run: pytest tests -m tutorials
diff --git a/.github/workflows/run-pytest-all.yml b/.github/workflows/run-pytest-all.yml
index 6a97e605c7..6cec1f6f00 100644
--- a/.github/workflows/run-pytest-all.yml
+++ b/.github/workflows/run-pytest-all.yml
@@ -28,10 +28,10 @@ jobs:
     # Steps represent a sequence of tasks that will be executed as part of the job
     steps:
     # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     
     - name: Setup Python ${{ matrix.python-version}}
-      uses: actions/setup-python@v3
+      uses: actions/setup-python@v5
       with:
         # Version range or exact version of a Python version to use, using SemVer's version range syntax.
         python-version: ${{ matrix.python-version}}
@@ -41,14 +41,11 @@ jobs:
     # Cache pip dependencies
     # Source: https://github.com/actions/cache/blob/main/examples.md#python---pip
     - name: Cache dependencies
-      uses: actions/cache@v3
+      uses: actions/cache@v4
       id: cache
       with:
-        # Cache the full python environment, this is more efficient than just caching pip
-        # https://blog.allenai.org/python-caching-in-github-actions-e9452698e98d
         path: ${{ env.pythonLocation }}
-        # Hash requirements.
-        key: ${{ env.pythonLocation }}-${{ hashFiles('setup.py') }}-${{ hashFiles('**/requirements.txt') }}
+        key: ${{ env.pythonLocation }}-${{ hashFiles('pyproject.toml') }}
 
     - name: Install external libraries
       run: |
@@ -56,29 +53,18 @@ jobs:
         sudo apt-get install libglu1-mesa libgeos-dev libffi-dev
         export LD_LIBRARY_PATH=/usr/local/lib64/:$LD_LIBRARY_PATH
 
-    - name: Install requirements
+    # Install the requirements, including the extras needed for testing, and PorePy
+    # itself.
+    - name: Install requirements and PorePy
       run: |
         pip install -U pip
-        # Install requirements, using eager updates to avoid stalled dependencies due to caching
-        # https://blog.allenai.org/python-caching-in-github-actions-e9452698e98d
-        pip install --upgrade --upgrade-strategy eager -r requirements.txt
+        pip install .[testing]
+        pip freeze
 
     - name: Install pypardiso
       run: |
         pip install pypardiso
 
-    - name: Install jupyter
-      run: |
-        pip install jupyter 
-
-    - name: Install PorePy
-      run: |
-        pip install .
-        pip freeze       
-
     - name: all tests
       if: ${{always()}}
-      run: pytest --run-skipped -m "not tutorials"
-
-
-
+      run: pytest --run-skipped -m "not tutorials"
\ No newline at end of file
diff --git a/.github/workflows/run-pytest.yml b/.github/workflows/run-pytest.yml
index f74978d89d..eaacc504db 100644
--- a/.github/workflows/run-pytest.yml
+++ b/.github/workflows/run-pytest.yml
@@ -36,10 +36,10 @@ jobs:
     # Steps represent a sequence of tasks that will be executed as part of the job
     steps:
     # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     
     - name: Setup Python ${{ matrix.python-version}}
-      uses: actions/setup-python@v3
+      uses: actions/setup-python@v5
       with:
         # Version range or exact version of a Python version to use, using SemVer's version range syntax.
         python-version: ${{ matrix.python-version}}
@@ -49,14 +49,11 @@ jobs:
     # Cache pip dependencies
     # Source: https://github.com/actions/cache/blob/main/examples.md#python---pip
     - name: Cache dependencies
-      uses: actions/cache@v3
+      uses: actions/cache@v4
       id: cache
       with:
-        # Cache the full python environment, this is more efficient than just caching pip
-        # https://blog.allenai.org/python-caching-in-github-actions-e9452698e98d
         path: ${{ env.pythonLocation }}
-        # Hash requirements.
-        key: ${{ env.pythonLocation }}-${{ hashFiles('setup.py') }}-${{ hashFiles('**/requirements.txt') }}
+        key: ${{ env.pythonLocation }}-${{ hashFiles('pyproject.toml') }}
 
     - name: Install external libraries
       run: |
@@ -64,17 +61,13 @@ jobs:
         sudo apt-get install libglu1-mesa libgeos-dev libffi-dev
         export LD_LIBRARY_PATH=/usr/local/lib64/:$LD_LIBRARY_PATH
 
-    - name: Install requirements
+    # Install the requirements, including the extras needed for testing, and PorePy
+    # itself.
+    - name: Install requirements and PorePy
       run: |
         pip install -U pip
-        # Install requirements, using eager updates to avoid stalled dependencies due to caching
-        # https://blog.allenai.org/python-caching-in-github-actions-e9452698e98d
-        pip install --upgrade --upgrade-strategy eager -r requirements.txt
-
-    - name: Install PorePy
-      run: |
-        pip install .
-        pip freeze       
+        pip install .[testing]
+        pip freeze
 
     - name: Run tests
       if: ${{always()}}
diff --git a/.github/workflows/run-static-checks.yml b/.github/workflows/run-static-checks.yml
index a6a719bc48..6695b84904 100644
--- a/.github/workflows/run-static-checks.yml
+++ b/.github/workflows/run-static-checks.yml
@@ -1,4 +1,4 @@
-# Run static tests for code quality: Isort, black, flake8, mypy.
+# Run static tests for code quality: Isort, ruff (formatting and linting), mypy.
 
 name: Static tests
 
@@ -55,10 +55,10 @@ jobs:
     # Steps represent a sequence of tasks that will be executed as part of the job
     steps:
     # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it.
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     
     - name: Setup Python ${{ matrix.python-version}}
-      uses: actions/setup-python@v3
+      uses: actions/setup-python@v5
       with:
         # Version range or exact version of a Python version to use, using SemVer's version range syntax.
         python-version: ${{ matrix.python-version}}
@@ -68,43 +68,51 @@ jobs:
     # Cache pip dependencies
     # Source: https://github.com/actions/cache/blob/main/examples.md#python---pip
     - name: Cache dependencies
-      uses: actions/cache@v3
+      uses: actions/cache@v4
       id: cache
       with:
-        # Cache the full python environment, this is more efficient than just caching pip
-        # https://blog.allenai.org/python-caching-in-github-actions-e9452698e98d
         path: ${{ env.pythonLocation }}
-        key: ${{ env.pythonLocation }}-${{ hashFiles('setup.py') }}-${{ hashFiles('**/requirements.txt', '**/requirements-dev.txt') }}
+        key: ${{ env.pythonLocation }}-${{ hashFiles('pyproject.toml') }}
 
     - name: Cache mypy
-      uses: actions/cache@v3
+      uses: actions/cache@v4
       id: cache-mypy
       with:
         path: mypy-cache-${{ matrix.python-version}}
         key: ${{ runner.os }}-mypy
 
-    # Install the requirements. Mypy needs this to access stubs for numpy etc.
+    # Install the dependencies, which is needed to make mypy happy. Unravelling of the
+    # second command: The expression in the angular brackets is a python command that
+    # fetches the dependencies section in pyproject.toml; this is then made available to
+    # pip via -r. Option --no-bulid-isolation makes pip not use an isolated environment
+    # for the install (hence the installed packages are sure to be available also after
+    # the build), while --no-deps will only install the specified packages (in this
+    # case, the dependencies of PorePy). The third line similarly installs the
+    # development dependencies.
+    #
+    # Thanks, copilot!
     - name: Install requirements
       run: |
         pip install -U pip
-        # Install requirements, using eager updates to avoid stalled dependencies due to caching
-        # https://blog.allenai.org/python-caching-in-github-actions-e9452698e98d
-        pip install --upgrade --upgrade-strategy eager -r requirements.txt
+        pip install toml
+        pip install --no-build-isolation --no-deps -r <(python -c "import toml; print('\n'.join(toml.load('pyproject.toml')['project']['dependencies']))")
+        pip install --no-build-isolation --no-deps -r <(python -c "import toml; print('\n'.join(toml.load('pyproject.toml')['project']['optional-dependencies']['development']))")
         pip freeze
 
-    # Run black, isort, flake8, and mypy.
-    # Note that since these are static checks, we have not installed PorePy.
-    - name: black
+    # Run the various checks
+    - name: Ruff linting
       if: ${{ always() }}
-      run: black --check src
+      run: ruff check src
 
-    - name: isort
+    - name: Ruff formatting
       if: ${{always()}}
-      run: isort --check src
+      run: ruff format --check src
 
-    - name: flake8
+      # EK note to self: Ruff may introduce sorting of imports in the future, if so, we
+      # should use that instead of isort.
+    - name: isort 
       if: ${{ always() }}
-      run: flake8 src 
+      run: isort src 
 
     - name: mypy
       if: ${{ always() }}
diff --git a/Readme.md b/Readme.md
index 288711c4da..656401eee1 100644
--- a/Readme.md
+++ b/Readme.md
@@ -1,6 +1,7 @@
 ![Pytest](https://github.com/pmgbergen/porepy/actions/workflows/run-pytest.yml/badge.svg)
 ![Pytest including slow](https://github.com/pmgbergen/porepy/actions/workflows/run-pytest-all.yml/badge.svg)
-![Mypy, black, isort, flake8](https://github.com/pmgbergen/porepy/actions/workflows/run-static-checks.yml/badge.svg)
+![Mypy, ruff, isort](https://github.com/pmgbergen/porepy/actions/workflows/run-static-checks.yml/badge.svg)
+![Tutorials](https://github.com/pmgbergen/porepy/actions/workflows/check_tutorials.yml/badge.svg)
 [![DOI](https://zenodo.org/badge/89228838.svg)](https://zenodo.org/badge/latestdoi/89228838)
 [![License: GPL v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
 
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000000..495817cc43
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,81 @@
+[build-system]
+requires = ["setuptools>=42", "wheel"]
+build-backend = "setuptools.build_meta"
+
+# Project information.
+[project]
+name = "porepy"
+version = "1.10.0"
+description = "Simulation tool for fractured and deformable porous media"
+readme = "Readme.md"
+license = { text = "GPL" }
+keywords = ["porous media simulation fractures deformable"]
+maintainers = [
+    { name = "Eirik Keilegavlen", email = "Eirik.Keilegavlen@uib.no" },
+    { name = "Ivar Stefansson", email = "Ivar.Stefansson@uib.no" },
+    { name = "Jakub Both", email = "Jakub.Both@uib.no" }
+]
+urls = { "Homepage" = "https://github.com/pmgbergen/porepy" }
+classifiers = [
+    "Programming Language :: Python",
+    "Programming Language :: Python :: 3",
+    "License :: OSI Approved :: GNU General Public License (GPL)",
+    "Operating System :: OS Independent"
+]
+dependencies = [
+    "deepdiff",
+    "future",
+    "gmsh",
+    "matplotlib < 3.11",
+    "meshio",
+    "networkx",
+    "numba",
+    "numpy < 2.2",
+    "scipy",
+    "seaborn",
+    "shapely",
+    "six",
+    "sympy",
+    "typing_extensions",
+]
+
+[project.optional-dependencies]
+development = [
+    "isort",
+    "ruff",
+    "mypy",
+    "mypy-extensions"
+]
+testing = [
+    "pytest >= 4.6",
+    "pytest-cov",
+    "pytest-runner"
+]
+
+# Information for the build system.
+[tool.setuptools.packages.find]
+where = ["src"]
+
+[tool.setuptools.package-data]
+"porepy" = [
+    "py.typed",
+    "applications/md_grids/gmsh_file_library/**/*.csv",
+    "applications/md_grids/gmsh_file_library/**/*.geo"
+]
+
+# Testing.
+[tool.pytest.ini_options]
+python_files = "tests/*.py"
+addopts = "--cov=src/porepy --cov-report term-missing -p no:warnings"
+
+# Static analysis.
+[tool.ruff.lint]
+ignore = [
+    "E731",  # do not assign a lambda expression, use a def
+    "E266",  # too many leading '#' for block comment
+]
+exclude = [
+    "src/porepy/__init__.py",
+    "src/porepy/numerics/ad/__init__.py",
+    "src/porepy/compositional/__init__.py"
+]
diff --git a/requirements.txt b/requirements.txt
deleted file mode 100644
index 9b5a77f404..0000000000
--- a/requirements.txt
+++ /dev/null
@@ -1,22 +0,0 @@
-future
-gmsh
-matplotlib < 3.11
-meshio >= 5.0
-networkx
-numba >= 0.57
-numpy
-requests
-scipy
-seaborn
-shapely
-six
-sympy
-typing_extensions
-black == 24.*
-deepdiff
-flake8
-isort
-mypy
-pytest >= 4.6
-pytest-cov
-pytest-runner
diff --git a/setup.cfg b/setup.cfg
deleted file mode 100644
index ab65e82880..0000000000
--- a/setup.cfg
+++ /dev/null
@@ -1,40 +0,0 @@
-[metadata]
-description-file = Readme.md
-
-[aliases]
-test=pytest
-
-[tool:pytest]
-python_files = tests/*.py
-addopts = --cov=src/porepy --cov-report term-missing -p no:warnings
-
-[flake8]
-max-line-length = 95
-ignore = 
-  # F541: f-string is missing placeholders
-  F541,
-  # W503: line break before binary operator
-  W503,
-  # E731 do not assign a lambda expression, use a def
-  E731,
-  # E203 whitespace before ':'
-  E203,
-  # E266 too many leading '#' for block comment
-  E266,
-  # allow to use \dot
-  W605
-  # Multiple statements on each line - this conflicts with black's preferences
-  E701
-  E704
-exclude = 
-  src/porepy/__init__.py,
-  src/porepy/numerics/ad/__init__.py,
-per-file-ignores =
-  # F401: imported but unused on files __init__.py
-  __init__.py:F401
-  # F403: from .module import *, ignore because of definition of __all__
-  src/porepy/compositional/__init__.py:F403,F401
-  
-[mypy]
-warn_unused_configs = True
-plugins = numpy.typing.mypy_plugin 
diff --git a/setup.py b/setup.py
index 8346773d17..b3e9a79cc3 100644
--- a/setup.py
+++ b/setup.py
@@ -1,30 +1,5 @@
 """Set-up file for PorePy for installations usins ``pip install .``"""
-from setuptools import find_packages, setup
 
+from setuptools import setup
 
-with open("requirements.txt") as f:
-    required = f.read().splitlines()
-
-
-setup(
-    name="porepy",
-    url="https://github.com/pmgbergen/porepy",
-    version="1.10.0",
-    license="GPL",
-    keywords=["porous media simulation fractures deformable"],
-    install_requires=required,
-    description="Simulation tool for fractured and deformable porous media",
-    maintainer="Eirik Keilegavlen",
-    maintainer_email="Eirik.Keilegavlen@uib.no",
-    platforms=["Linux", "Windows", "Mac OS-X"],
-    package_data={
-        "porepy": [
-            "py.typed",
-            "applications/md_grids/gmsh_file_library/**/*.csv",
-            "applications/md_grids/gmsh_file_library/**/*.geo",
-        ],
-    },
-    packages=find_packages("src"),
-    package_dir={"": "src"},
-    zip_safe=False,
-)
+setup()
diff --git a/src/porepy/__init__.py b/src/porepy/__init__.py
index 85e1336ca5..bd5b69a974 100644
--- a/src/porepy/__init__.py
+++ b/src/porepy/__init__.py
@@ -1,4 +1,4 @@
-"""   PorePy.
+"""PorePy.
 
 Root directory for the PorePy package. Contains the following sub-packages:
 
diff --git a/src/porepy/applications/boundary_conditions/__init__.py b/src/porepy/applications/boundary_conditions/__init__.py
index 43ea7a9b9e..cf9b1ccbd6 100644
--- a/src/porepy/applications/boundary_conditions/__init__.py
+++ b/src/porepy/applications/boundary_conditions/__init__.py
@@ -1 +1,3 @@
 from . import model_boundary_conditions
+
+__all__ = ["model_boundary_conditions"]
diff --git a/src/porepy/applications/test_utils/__init__.py b/src/porepy/applications/test_utils/__init__.py
index 80c72cf922..03df7e557b 100644
--- a/src/porepy/applications/test_utils/__init__.py
+++ b/src/porepy/applications/test_utils/__init__.py
@@ -15,3 +15,14 @@
     reference_sparse_arrays,
     vtk,
 )
+
+__all__ = [
+    "arrays",
+    "common_xpfa_tests",
+    "grids",
+    "models",
+    "partial_discretization",
+    "reference_dense_arrays",
+    "reference_sparse_arrays",
+    "vtk",
+]
diff --git a/src/porepy/applications/test_utils/arrays.py b/src/porepy/applications/test_utils/arrays.py
index e94dcc3851..b2d9ded6d8 100644
--- a/src/porepy/applications/test_utils/arrays.py
+++ b/src/porepy/applications/test_utils/arrays.py
@@ -1,5 +1,4 @@
-"""Test helpers for arrays and matrices.
-"""
+"""Test helpers for arrays and matrices."""
 
 import numpy as np
 import scipy.sparse as sps
diff --git a/src/porepy/applications/test_utils/common_xpfa_tests.py b/src/porepy/applications/test_utils/common_xpfa_tests.py
index e4f8295313..c6073e607c 100644
--- a/src/porepy/applications/test_utils/common_xpfa_tests.py
+++ b/src/porepy/applications/test_utils/common_xpfa_tests.py
@@ -56,7 +56,7 @@ def perturb_grid(g, rate, dx):
 
 
 def create_grid_mpfa_mpsa_reproduce_known_values(
-    grid_type: Literal["cart", "simplex"]
+    grid_type: Literal["cart", "simplex"],
 ) -> tuple[pp.Grid, pp.Grid]:
     """Create grids for the tests that mpfa and mpsa reproduce known values.
 
@@ -822,7 +822,7 @@ def test_linear_pressure_dirichlet_conditions_perturbed_grid(self):
 
 
 def test_split_discretization_into_subproblems(
-    discr_class: Union[pp.Mpfa, pp.Mpsa, pp.Biot]
+    discr_class: Union[pp.Mpfa, pp.Mpsa, pp.Biot],
 ):
     """Test that the discretization matrices produced by Mpxa are the same if they
     are split into subproblems or not.
diff --git a/src/porepy/applications/test_utils/grids.py b/src/porepy/applications/test_utils/grids.py
index 5c42e3cfc3..d565910954 100644
--- a/src/porepy/applications/test_utils/grids.py
+++ b/src/porepy/applications/test_utils/grids.py
@@ -1,5 +1,4 @@
-"""Helper methods to compare grids.
-"""
+"""Helper methods to compare grids."""
 
 import numpy as np
 
diff --git a/src/porepy/applications/test_utils/vtk.py b/src/porepy/applications/test_utils/vtk.py
index a25829be43..d75e7fa3d7 100644
--- a/src/porepy/applications/test_utils/vtk.py
+++ b/src/porepy/applications/test_utils/vtk.py
@@ -1,5 +1,4 @@
-"""Test helpers for the vtk files.
-"""
+"""Test helpers for the vtk files."""
 
 import shutil
 import xml.etree.ElementTree as ET
diff --git a/src/porepy/applications/test_utils/well_models.py b/src/porepy/applications/test_utils/well_models.py
index cc9884dbd3..7750cc9ad5 100644
--- a/src/porepy/applications/test_utils/well_models.py
+++ b/src/porepy/applications/test_utils/well_models.py
@@ -1,5 +1,4 @@
-"""Contains code for setting up a simple but non-trivial model with a well.
-"""
+"""Contains code for setting up a simple but non-trivial model with a well."""
 
 from typing import Literal
 
@@ -9,7 +8,6 @@
 
 
 class OneVerticalWell(pp.PorePyModel):
-
     def set_well_network(self) -> None:
         """Assign well network class."""
         points = np.array([[0.5, 0.5], [0.5, 0.5], [0.2, 1]])
@@ -154,7 +152,6 @@ def bc_type_fourier_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
 
 
 class WellPermeability(pp.constitutive_laws.CubicLawPermeability):
-
     def permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         """Permeability [m^2].
 
diff --git a/src/porepy/compositional/base.py b/src/porepy/compositional/base.py
index ead671c812..43a27b5991 100644
--- a/src/porepy/compositional/base.py
+++ b/src/porepy/compositional/base.py
@@ -98,7 +98,6 @@ class Component:
     """
 
     def __init__(self, *args, **kwargs) -> None:
-
         self.name: str = str(kwargs.get("name", "unnamed_component"))
         """Name of the component. Can be named by providing a keyword argument 'name'
         when instantiating."""
@@ -164,7 +163,6 @@ class Compound(Component, Generic[ComponentLike]):
     """
 
     def __init__(self, *args, **kwargs) -> None:
-
         self.molar_mass: pp.number
         if "molar_mass" not in kwargs:
             raise ValueError(
@@ -260,7 +258,7 @@ def compound_molar_mass(self, domains: pp.SubdomainsOrBoundaries) -> pp.ad.Opera
             # fluid component.
             if not hasattr(pc, "molar_mass"):
                 raise TypeError(
-                    f"Cannot assemble compound molar mass: Active tracer of type"
+                    "Cannot assemble compound molar mass: Active tracer of type"
                     + f" {type(pc)} has no attribute `molar_mass`."
                 )
             else:
@@ -478,7 +476,6 @@ def __init__(
         state: PhysicalState,
         name: str,
     ) -> None:
-
         self._ref_component_index: int = 0
         """See :meth:`reference_component_index`."""
 
@@ -734,7 +731,6 @@ def __init__(
         components: list[ComponentLike],
         phases: list[PhaseLike],
     ) -> None:
-
         self._ref_phase_index: int = 0
         """See :meth:`reference_phase_index`."""
         self._ref_component_index: int = 0
@@ -961,7 +957,6 @@ def density(self, domains: pp.SubdomainsOrBoundaries) -> pp.ad.Operator:
         """
 
         if self.num_phases > 1:
-
             op = pp.ad.sum_operator_list(
                 [
                     phase.saturation(domains) * phase.density(domains)
@@ -1012,7 +1007,6 @@ def specific_enthalpy(self, domains: pp.SubdomainsOrBoundaries) -> pp.ad.Operato
         """
 
         if self.num_phases > 1:
-
             op = pp.ad.sum_operator_list(
                 [
                     phase.fraction(domains) * phase.specific_enthalpy(domains)
@@ -1051,7 +1045,6 @@ def thermal_conductivity(
 
         """
         if self.num_phases > 1:
-
             op = pp.ad.sum_operator_list(
                 [
                     phase.saturation(domains) * phase.thermal_conductivity(domains)
@@ -1061,7 +1054,6 @@ def thermal_conductivity(
             )
 
         else:
-
             op = self.reference_phase.thermal_conductivity(domains)
             op.set_name("fluid_thermal_conductivity")
 
diff --git a/src/porepy/compositional/compositional_mixins.py b/src/porepy/compositional/compositional_mixins.py
index c45a1ab39f..6cf0f89713 100644
--- a/src/porepy/compositional/compositional_mixins.py
+++ b/src/porepy/compositional/compositional_mixins.py
@@ -783,9 +783,9 @@ def tracer_fraction(
             A callable which returns the tracer fraction for a given set of domains.
 
         """
-        assert (
-            tracer in compound.active_tracers
-        ), f"Solute {tracer.name} not in compound {compound.name}"
+        assert tracer in compound.active_tracers, (
+            f"Solute {tracer.name} not in compound {compound.name}"
+        )
 
         if self.has_independent_tracer_fraction(tracer, compound):
             fraction = self._fraction_factory(
@@ -909,9 +909,9 @@ def extended_fraction(
             A callable which returns the extended fraction for a given set of domains.
 
         """
-        assert (
-            component in phase
-        ), f"Component {component.name} not in phase {phase.name}"
+        assert component in phase, (
+            f"Component {component.name} not in phase {phase.name}"
+        )
 
         fraction: DomainFunctionType
 
@@ -971,9 +971,9 @@ def partial_fraction(
             A callable which returns the extended fraction for a given set of domains.
 
         """
-        assert (
-            component in phase
-        ), f"Component {component.name} not in phase {phase.name}"
+        assert component in phase, (
+            f"Component {component.name} not in phase {phase.name}"
+        )
 
         fraction: DomainFunctionType
 
diff --git a/src/porepy/examples/__init__.py b/src/porepy/examples/__init__.py
index 3ec1f4cce7..bdee8d8c92 100644
--- a/src/porepy/examples/__init__.py
+++ b/src/porepy/examples/__init__.py
@@ -16,3 +16,11 @@
 from .flow_benchmark_3d_case_3 import FlowBenchmark3dCase3Model
 from .mandel_biot import MandelExactSolution, MandelSolutionStrategy
 from .terzaghi_biot import TerzaghiExactSolution, TerzaghiSetup
+
+__all__ = [
+    "FlowBenchmark3dCase3Model",
+    "MandelExactSolution",
+    "MandelSolutionStrategy",
+    "TerzaghiExactSolution",
+    "TerzaghiSetup",
+]
diff --git a/src/porepy/examples/flow_benchmark_3d_case_3.py b/src/porepy/examples/flow_benchmark_3d_case_3.py
index 88de040a56..00d20cf2de 100644
--- a/src/porepy/examples/flow_benchmark_3d_case_3.py
+++ b/src/porepy/examples/flow_benchmark_3d_case_3.py
@@ -68,7 +68,6 @@ def set_geometry(self) -> None:
 
 
 class IntersectionPermeability(Permeability):
-
     def intersection_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         """Constant intersection permeability.
 
diff --git a/src/porepy/examples/mandel_biot.py b/src/porepy/examples/mandel_biot.py
index 5b7893251e..d20f984463 100644
--- a/src/porepy/examples/mandel_biot.py
+++ b/src/porepy/examples/mandel_biot.py
@@ -1251,7 +1251,6 @@ def grid_type(self) -> Literal["simplex", "cartesian", "tensor_grid"]:
 
 # -----> Boundary conditions
 class MandelBoundaryConditionsMechanicsTimeDependent(pp.PorePyModel):
-
     exact_sol: MandelExactSolution
     """Exact solution object."""
 
@@ -1333,7 +1332,6 @@ def bc_values_displacement(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
 
 
 class MandelBoundaryConditionsSinglePhaseFlow(pp.PorePyModel):
-
     def bc_type_darcy_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         """Define boundary condition types for the Darcy flux.
 
diff --git a/src/porepy/fracs/fracture_network.py b/src/porepy/fracs/fracture_network.py
index 57eb56c8f3..57c879844a 100644
--- a/src/porepy/fracs/fracture_network.py
+++ b/src/porepy/fracs/fracture_network.py
@@ -6,14 +6,14 @@
 from __future__ import annotations
 
 import warnings
-from typing import List, Optional, Union
+from typing import Optional, Union, cast
 
 import porepy as pp
 from porepy.fracs.fracture_network_2d import FractureNetwork2d
 from porepy.fracs.fracture_network_3d import FractureNetwork3d
 
 # Custom typings
-FractureList = Union[List[pp.LineFracture], List[pp.PlaneFracture]]
+FractureList = Union[list[pp.LineFracture], list[pp.PlaneFracture]]
 FractureNetwork = Union[FractureNetwork2d, FractureNetwork3d]
 
 
@@ -127,12 +127,14 @@ def create_fracture_network(
     # check above, it is safe to ignore the argument type
     if dim == 2:
         fracture_network_2d = FractureNetwork2d(
-            fractures=fracs, domain=domain, tol=tol  # type: ignore[arg-type]
+            fractures=cast(list[pp.LineFracture], fracs),
+            domain=domain,
+            tol=tol,  # type: ignore[arg-type]
         )
         return fracture_network_2d
     else:
         fracture_network_3d = FractureNetwork3d(
-            fractures=fracs,  # type: ignore[arg-type]
+            fractures=cast(list[pp.PlaneFracture], fracs),
             domain=domain,
             tol=tol,
             run_checks=run_checks,
diff --git a/src/porepy/fracs/fracture_network_2d.py b/src/porepy/fracs/fracture_network_2d.py
index dce24e782a..fb4cd16ccd 100644
--- a/src/porepy/fracs/fracture_network_2d.py
+++ b/src/porepy/fracs/fracture_network_2d.py
@@ -506,9 +506,9 @@ def _find_and_split_intersections(self, constraints: np.ndarray) -> None:
 
         tags = np.zeros((2, edges.shape[1]), dtype=int)
 
-        tags[0][
-            np.logical_not(self.tags["boundary"])
-        ] = GmshInterfaceTags.FRACTURE.value
+        tags[0][np.logical_not(self.tags["boundary"])] = (
+            GmshInterfaceTags.FRACTURE.value
+        )
         tags[0][self.tags["boundary"]] = GmshInterfaceTags.DOMAIN_BOUNDARY_LINE.value
         tags[0][constraints] = GmshInterfaceTags.AUXILIARY_LINE.value
 
diff --git a/src/porepy/fracs/fracture_network_3d.py b/src/porepy/fracs/fracture_network_3d.py
index e275ce57f3..3c406858e3 100644
--- a/src/porepy/fracs/fracture_network_3d.py
+++ b/src/porepy/fracs/fracture_network_3d.py
@@ -1182,7 +1182,9 @@ def _remove_edge_intersections(
             # somewhat to obtain a more robust algorithm. Not sure about how to do
             # this consistent.
             p_new, edges_new, tags = pp.intersections.split_intersecting_segments_2d(
-                p_2d, edges_2d, tol=self.tol  # type: ignore[misc]
+                p_2d,
+                edges_2d,
+                tol=self.tol,  # type: ignore[misc]
             )
             # Then, patch things up by converting new points to 3D,
 
diff --git a/src/porepy/fracs/meshing.py b/src/porepy/fracs/meshing.py
index 1208ff0c95..3bb67eb1b6 100644
--- a/src/porepy/fracs/meshing.py
+++ b/src/porepy/fracs/meshing.py
@@ -240,7 +240,9 @@ def _tag_faces(grids, check_highest_dim=True):
 
     # Assume only one grid of highest dimension
     if check_highest_dim:
-        assert len(grids[0]) == 1, "Must be exactly" "1 grid of dim: " + str(len(grids))
+        assert len(grids[0]) == 1, "There must be exactly 1 grid of dim: " + str(
+            len(grids)
+        )
 
     for g_h in np.atleast_1d(grids[0]):
         bnd_faces = g_h.get_all_boundary_faces()
diff --git a/src/porepy/fracs/simplex.py b/src/porepy/fracs/simplex.py
index 1600a542fd..226b989221 100644
--- a/src/porepy/fracs/simplex.py
+++ b/src/porepy/fracs/simplex.py
@@ -1,5 +1,5 @@
 """Module with functionality to create simplex grids with fractures from a
-tesselation generated by gmsh. """
+tesselation generated by gmsh."""
 
 from __future__ import annotations
 
diff --git a/src/porepy/fracs/split_grid.py b/src/porepy/fracs/split_grid.py
index 8761a67fed..0df261a9a3 100644
--- a/src/porepy/fracs/split_grid.py
+++ b/src/porepy/fracs/split_grid.py
@@ -558,7 +558,7 @@ def update_cell_connectivity(
     # domain. Check that each face added two cells:
     if sum(left_cell) * 2 != left_cell.size:
         raise ValueError(
-            "Fractures must either be" "on boundary or completely inside domain"
+            "Fractures must either be on boundary or completely inside domain"
         )
 
     # We create a cell_faces mapping for the new faces. This will be added on the end
diff --git a/src/porepy/fracs/structured.py b/src/porepy/fracs/structured.py
index 083fda968a..d06c58257c 100644
--- a/src/porepy/fracs/structured.py
+++ b/src/porepy/fracs/structured.py
@@ -213,9 +213,9 @@ def _create_lower_dim_grids_3d(
         is_xy_frac = np.allclose(f[2, 0], f[2])
         is_xz_frac = np.allclose(f[1, 0], f[1])
         is_yz_frac = np.allclose(f[0, 0], f[0])
-        assert (
-            is_xy_frac + is_xz_frac + is_yz_frac == 1
-        ), "Fracture must align to x-, y- or z-axis"
+        assert is_xy_frac + is_xz_frac + is_yz_frac == 1, (
+            "Fracture must align to x-, y- or z-axis."
+        )
 
         # snap to grid
         if physdims is None:
@@ -356,10 +356,10 @@ def _create_lower_dim_grids_3d(
         e_pt = pts[:, edges[1, e]]
         nodes = _find_nodes_on_line(g_3d, nx, s_pt, e_pt)
         loc_coord = g_3d.nodes[:, nodes]
-        assert (
-            loc_coord.shape[1] > 1
-        ), "1d grid in intersection should span\
+        assert loc_coord.shape[1] > 1, (
+            "1d grid in intersection should span\
             more than one node"
+        )
         g = msh_2_grid.create_embedded_line_grid(loc_coord, nodes)
         g_1d.append(g)
 
diff --git a/src/porepy/fracs/utils.py b/src/porepy/fracs/utils.py
index cf6f91376e..5f8c77b343 100644
--- a/src/porepy/fracs/utils.py
+++ b/src/porepy/fracs/utils.py
@@ -1,5 +1,5 @@
 """This module contains (frontend) utility functions related to fractures and their
-meshing. """
+meshing."""
 
 from __future__ import annotations
 
diff --git a/src/porepy/geometry/intersections.py b/src/porepy/geometry/intersections.py
index 7d388d4a78..c3e4073ba3 100644
--- a/src/porepy/geometry/intersections.py
+++ b/src/porepy/geometry/intersections.py
@@ -2062,7 +2062,9 @@ def _min_max_coord(coord):
                 is_ccw = np.array(
                     [
                         pp.geometry_property_checks.is_ccw_polyline(
-                            start[:, i], middle[:, i], end[:, i]  # type:ignore
+                            start[:, i],  # type: ignore
+                            middle[:, i],  # type: ignore
+                            end[:, i],  # type:ignore
                         )
                         for i in range(poly.shape[1])  # type:ignore
                     ]
diff --git a/src/porepy/geometry/map_geometry.py b/src/porepy/geometry/map_geometry.py
index 06e984abad..452d42733f 100644
--- a/src/porepy/geometry/map_geometry.py
+++ b/src/porepy/geometry/map_geometry.py
@@ -10,7 +10,7 @@
 
 
 def force_point_collinearity(
-    pts: np.ndarray[Any, np.dtype[np.float64]]
+    pts: np.ndarray[Any, np.dtype[np.float64]],
 ) -> np.ndarray[Any, np.dtype[np.float64]]:
     """Given a set of points, return them aligned on a line.
 
@@ -555,7 +555,7 @@ def compute_normal(
 
 
 def compute_normals_1d(
-    pts: np.ndarray[Any, np.dtype[np.float64]]
+    pts: np.ndarray[Any, np.dtype[np.float64]],
 ) -> np.ndarray[Any, np.dtype[np.float64]]:
     """Compute the normals of a set of points aligned along a 1d line.
 
diff --git a/src/porepy/grids/coarsening.py b/src/porepy/grids/coarsening.py
index 69a91f9cb9..196190c8f2 100644
--- a/src/porepy/grids/coarsening.py
+++ b/src/porepy/grids/coarsening.py
@@ -64,9 +64,7 @@ def coarsen(
         if method_kwargs.get("if_seeds", False):
             seeds = generate_seeds(g)
         matrix = _tpfa_matrix(g)
-        partition = create_partition(
-            matrix, g, seeds=seeds, **method_kwargs
-        )  # type: ignore
+        partition = create_partition(matrix, g, seeds=seeds, **method_kwargs)  # type: ignore
 
     else:
         raise ValueError(f"Undefined method `{method}` for coarsening algorithm.")
diff --git a/src/porepy/grids/md_grid.py b/src/porepy/grids/md_grid.py
index 9842d4f667..a6919793d1 100644
--- a/src/porepy/grids/md_grid.py
+++ b/src/porepy/grids/md_grid.py
@@ -976,7 +976,7 @@ def __repr__(self) -> str:
                 for sd in self.subdomains(dim=dim):
                     num_sd += 1
                     nc += sd.num_cells
-                s += f"{num_sd} grids of dimension {dim}" f" with in total {nc} cells\n"
+                s += f"{num_sd} grids of dimension {dim} with in total {nc} cells\n"
         if self.num_interfaces() > 0:
             for dim in range(self.dim_max(), self.dim_min(), -1):
                 num_intf = 0
diff --git a/src/porepy/grids/mdg_generation.py b/src/porepy/grids/mdg_generation.py
index 67fce4a0b1..6558057652 100644
--- a/src/porepy/grids/mdg_generation.py
+++ b/src/porepy/grids/mdg_generation.py
@@ -57,7 +57,7 @@ def _validate_args_types(
 
 
 def _validate_grid_type_value(
-    grid_type: Literal["simplex", "cartesian", "tensor_grid"]
+    grid_type: Literal["simplex", "cartesian", "tensor_grid"],
 ):
     """Validates grid_type value.
 
diff --git a/src/porepy/grids/mortar_grid.py b/src/porepy/grids/mortar_grid.py
index 7e9aff066e..afc71dec97 100644
--- a/src/porepy/grids/mortar_grid.py
+++ b/src/porepy/grids/mortar_grid.py
@@ -160,7 +160,7 @@ def __init__(
         """Nodes in the side grids with ``shape=(3, num_nodes)``."""
 
         # Set projections
-        if not (primary_secondary is None):
+        if primary_secondary is not None:
             self._init_projections(primary_secondary, face_duplicate_ind)
             self._set_projections()
 
diff --git a/src/porepy/grids/refinement.py b/src/porepy/grids/refinement.py
index 6f94c941f9..f8fde8a7ca 100644
--- a/src/porepy/grids/refinement.py
+++ b/src/porepy/grids/refinement.py
@@ -577,9 +577,9 @@ def structured_refinement(
     data = np.ones(indices.size)
     coarse_fine = sps.csc_matrix((data, indices, indptr))
 
-    assert (
-        indices.size == g_ref.num_cells
-    ), "Every fine cell should be inside exactly one coarse cell"
+    assert indices.size == g_ref.num_cells, (
+        "Every fine cell should be inside exactly one coarse cell"
+    )
     return coarse_fine
 
 
diff --git a/src/porepy/models/compositional_flow.py b/src/porepy/models/compositional_flow.py
index 1bc2287596..6d7b83089b 100644
--- a/src/porepy/models/compositional_flow.py
+++ b/src/porepy/models/compositional_flow.py
@@ -258,9 +258,9 @@ def set_equations(self) -> None:
         using the fractional flow framework."""
         super().set_equations()
 
-        assert self.params[
-            "rediscretize_darcy_flux"
-        ], "Model params['rediscretize_darcy_flux'] must be flagged as True."
+        assert self.params["rediscretize_darcy_flux"], (
+            "Model params['rediscretize_darcy_flux'] must be flagged as True."
+        )
 
     def fluid_flux(self, domains: pp.SubdomainsOrBoundaries) -> pp.ad.Operator:
         """The fluid flux is given solely by the :attr:`darcy_flux`, assuming the total
@@ -1560,9 +1560,9 @@ def before_nonlinear_iteration(self) -> None:
         self.update_thermodynamic_properties_of_phases()
         # NOTE: Mypy complaints about trivial body of protocol.
         # But this is a mixin. We assert it is indeed a solution strategy and proceed.
-        assert isinstance(
-            self, pp.SolutionStrategy
-        ), "This is a mixin. Require SolutionStrategy as base."
+        assert isinstance(self, pp.SolutionStrategy), (
+            "This is a mixin. Require SolutionStrategy as base."
+        )
         super().before_nonlinear_iteration()  # type:ignore[safe-super]
 
     def after_nonlinear_convergence(self) -> None:
@@ -1574,9 +1574,9 @@ def after_nonlinear_convergence(self) -> None:
         The progression is performed after the super-call.
 
         """
-        assert isinstance(
-            self, pp.SolutionStrategy
-        ), "This is a mixin. Require SolutionStrategy as base."
+        assert isinstance(self, pp.SolutionStrategy), (
+            "This is a mixin. Require SolutionStrategy as base."
+        )
         super().after_nonlinear_convergence()  # type:ignore[safe-super]
 
         subdomains = self.mdg.subdomains()
@@ -1620,9 +1620,9 @@ class SolutionStrategyNonlinearMPFA(pp.PorePyModel):
     """See :class:`~porepy.models.constitutive_laws.DarcysLaw`."""
 
     def __init__(self, params: Optional[dict] = None) -> None:
-        assert isinstance(
-            self, pp.SolutionStrategy
-        ), "This is a mixin. Require SolutionStrategy as base."
+        assert isinstance(self, pp.SolutionStrategy), (
+            "This is a mixin. Require SolutionStrategy as base."
+        )
         super().__init__(params)  # type:ignore[safe-super]
 
         self._nonlinear_flux_discretizations: list[pp.ad._ad_utils.MergedOperator] = []
@@ -1646,9 +1646,9 @@ def set_nonlinear_discretizations(self) -> None:
         """After the super-call, this method adds the
         :meth:`fourier_flux_discretization` and the :meth:`darcy_flux_discretization`
         to the update framework using :meth:`add_nonlinear_flux_discretization`."""
-        assert isinstance(
-            self, pp.SolutionStrategy
-        ), "This is a mixin. Require SolutionStrategy as base."
+        assert isinstance(self, pp.SolutionStrategy), (
+            "This is a mixin. Require SolutionStrategy as base."
+        )
         super().set_nonlinear_discretizations()  # type:ignore[safe-super]
 
         subdomains = self.mdg.subdomains()
@@ -1689,9 +1689,9 @@ def before_nonlinear_iteration(self) -> None:
 
         """
         self.rediscretize_fluxes()
-        assert isinstance(
-            self, pp.SolutionStrategy
-        ), "This is a mixin. Require SolutionStrategy as base."
+        assert isinstance(self, pp.SolutionStrategy), (
+            "This is a mixin. Require SolutionStrategy as base."
+        )
         super().before_nonlinear_iteration()  # type:ignore[safe-super]
 
 
@@ -1891,9 +1891,9 @@ def solve_linear_system(self) -> np.ndarray:
 
         # NOTE mypy complaints about trivial body of protocol.
         # But this is a mixin. We assert it is indeed a solution strategy and proceed
-        assert isinstance(
-            self, pp.SolutionStrategy
-        ), "This is a mixin. Require SolutionStrategy as base."
+        assert isinstance(self, pp.SolutionStrategy), (
+            "This is a mixin. Require SolutionStrategy as base."
+        )
         sol = super().solve_linear_system()  # type:ignore[safe-super]
 
         if self.params.get("reduce_linear_system", False):
diff --git a/src/porepy/models/mass_and_energy_balance.py b/src/porepy/models/mass_and_energy_balance.py
index 1a4bacf563..294d5330fd 100644
--- a/src/porepy/models/mass_and_energy_balance.py
+++ b/src/porepy/models/mass_and_energy_balance.py
@@ -1,4 +1,4 @@
-"""Combine single-physics models into coupled mass and energy balance equations. """
+"""Combine single-physics models into coupled mass and energy balance equations."""
 
 from __future__ import annotations
 
diff --git a/src/porepy/models/poromechanics.py b/src/porepy/models/poromechanics.py
index d441c592ba..d5e0cd2454 100644
--- a/src/porepy/models/poromechanics.py
+++ b/src/porepy/models/poromechanics.py
@@ -145,9 +145,9 @@ def set_discretization_parameters(self) -> None:
                 "scalar_vector_mappings", {}
             )
             scalar_vector_mappings[self.darcy_keyword] = self.biot_tensor([sd])
-            data[pp.PARAMETERS][self.stress_keyword][
-                "scalar_vector_mappings"
-            ] = scalar_vector_mappings
+            data[pp.PARAMETERS][self.stress_keyword]["scalar_vector_mappings"] = (
+                scalar_vector_mappings
+            )
 
     def _is_nonlinear_problem(self) -> bool:
         """The coupled problem is nonlinear."""
diff --git a/src/porepy/models/run_models.py b/src/porepy/models/run_models.py
index fecfaaed9b..d19c4424fc 100644
--- a/src/porepy/models/run_models.py
+++ b/src/porepy/models/run_models.py
@@ -1,4 +1,4 @@
-""" This module contains functions to run stationary and time-dependent models."""
+"""This module contains functions to run stationary and time-dependent models."""
 
 from __future__ import annotations
 
diff --git a/src/porepy/models/solution_strategy.py b/src/porepy/models/solution_strategy.py
index 15682453b2..af9bd5e754 100644
--- a/src/porepy/models/solution_strategy.py
+++ b/src/porepy/models/solution_strategy.py
@@ -289,9 +289,9 @@ def set_materials(self):
             dict[str, pp.Constants], self.params.get("material_constants", {})
         )
         # If the user provided material constants, assert they are in dictionary form
-        assert isinstance(
-            constants, dict
-        ), "model.params['material_constants'] must be a dictionary."
+        assert isinstance(constants, dict), (
+            "model.params['material_constants'] must be a dictionary."
+        )
 
         # Use standard models for fluid, solid and numerical constants if not provided.
         # Otherwise get the given constants.
diff --git a/src/porepy/models/thermoporomechanics.py b/src/porepy/models/thermoporomechanics.py
index 9482239944..db10c0ec17 100644
--- a/src/porepy/models/thermoporomechanics.py
+++ b/src/porepy/models/thermoporomechanics.py
@@ -175,9 +175,9 @@ def set_discretization_parameters(self) -> None:
                 self.solid_thermal_expansion_tensor([sd])
             )
             scalar_vector_mappings[self.darcy_keyword] = self.biot_tensor([sd])
-            data[pp.PARAMETERS][self.stress_keyword][
-                "scalar_vector_mappings"
-            ] = scalar_vector_mappings
+            data[pp.PARAMETERS][self.stress_keyword]["scalar_vector_mappings"] = (
+                scalar_vector_mappings
+            )
 
     def set_nonlinear_discretizations(self) -> None:
         """Collect discretizations for nonlinear terms."""
diff --git a/src/porepy/numerics/ad/__init__.py b/src/porepy/numerics/ad/__init__.py
index d930125e59..7f1677d9c3 100644
--- a/src/porepy/numerics/ad/__init__.py
+++ b/src/porepy/numerics/ad/__init__.py
@@ -1,4 +1,4 @@
-""" Init file for all AD functionality.
+"""Init file for all AD functionality.
 
 They should all be accessible through a calling
    >>> import porepy as pp
diff --git a/src/porepy/numerics/ad/_ad_parser.py b/src/porepy/numerics/ad/_ad_parser.py
index 147094353f..58859ec839 100644
--- a/src/porepy/numerics/ad/_ad_parser.py
+++ b/src/porepy/numerics/ad/_ad_parser.py
@@ -327,9 +327,9 @@ def _evaluate_single(
 
             case Operations.evaluate:
                 # Operator functions should have at least 1 child (themselves).
-                assert (
-                    len(child_values) >= 1
-                ), "Operator functions must have at least 1 child."
+                assert len(child_values) >= 1, (
+                    "Operator functions must have at least 1 child."
+                )
                 assert hasattr(op, "func"), (
                     f"Operators with operation {operation} must have a functional"
                     + " representation `func` implemented as a callable member."
diff --git a/src/porepy/numerics/ad/forward_mode.py b/src/porepy/numerics/ad/forward_mode.py
index 2a65475d6a..b01058c544 100644
--- a/src/porepy/numerics/ad/forward_mode.py
+++ b/src/porepy/numerics/ad/forward_mode.py
@@ -1,5 +1,4 @@
-"""The module contains the data class for forward mode automatic differentiation.
-"""
+"""The module contains the data class for forward mode automatic differentiation."""
 
 from __future__ import annotations
 
diff --git a/src/porepy/numerics/ad/grid_operators.py b/src/porepy/numerics/ad/grid_operators.py
index ce123476e5..db653f57e1 100644
--- a/src/porepy/numerics/ad/grid_operators.py
+++ b/src/porepy/numerics/ad/grid_operators.py
@@ -1,4 +1,4 @@
-""" Ad representation of grid-related quantities needed to write equations. The classes
+"""Ad representation of grid-related quantities needed to write equations. The classes
 defined here are mainly wrappers that constructs Ad matrices based on grid information.
 
 """
diff --git a/src/porepy/numerics/ad/operator_functions.py b/src/porepy/numerics/ad/operator_functions.py
index a592383ec6..b79504edcf 100644
--- a/src/porepy/numerics/ad/operator_functions.py
+++ b/src/porepy/numerics/ad/operator_functions.py
@@ -95,9 +95,9 @@ def __call__(self, *args: pp.ad.Operator) -> pp.ad.Operator:
             to make the numerical function available during parsing (see :meth:`parse`).
 
         """
-        assert (
-            len(args) > 0
-        ), "Operator functions must be called with at least 1 argument."
+        assert len(args) > 0, (
+            "Operator functions must be called with at least 1 argument."
+        )
 
         op = Operator(
             name=f"{self.name}{[a.name for a in args]}",
@@ -319,9 +319,9 @@ def get_values(self, *args: float | np.ndarray | AdArray) -> float | np.ndarray:
         return result.val if isinstance(result, AdArray) else result
 
     def get_jacobian(self, *args: float | np.ndarray | AdArray) -> sps.spmatrix:
-        assert any(
-            isinstance(a, AdArray) for a in args
-        ), "No Ad arrays passed as arguments."
+        assert any(isinstance(a, AdArray) for a in args), (
+            "No Ad arrays passed as arguments."
+        )
         result = self._func(*args)
         assert isinstance(result, AdArray)
         return result.jac
diff --git a/src/porepy/numerics/ad/operators.py b/src/porepy/numerics/ad/operators.py
index dfc5578f2d..56b896796a 100644
--- a/src/porepy/numerics/ad/operators.py
+++ b/src/porepy/numerics/ad/operators.py
@@ -1736,22 +1736,22 @@ def __init__(self, variables: list[Variable]) -> None:
 
         # check assumptions
         if len(variables) > 0:
-            assert (
-                len(set(time_indices)) == 1
-            ), "Cannot create md-variable from variables at different time steps."
+            assert len(set(time_indices)) == 1, (
+                "Cannot create md-variable from variables at different time steps."
+            )
             # NOTE both must be unique for all sub-variables, to avoid md-variables
             # having sub-variables at different iterate states.
             # Both current value, and most recent previous iterate have iterate index 0,
             # hence the need to check the size of the current_iter set.
-            assert (
-                len(set(iter_indices)) == 1 and len(set(current_iter)) == 1
-            ), "Cannot create md-variable from variables at different iterates."
-            assert (
-                len(set(names)) == 1
-            ), "Cannot create md-variable from variables with different names."
-            assert len(set(domains)) == len(
-                domains
-            ), "Cannot create md-variable from variables with overlapping domains."
+            assert len(set(iter_indices)) == 1 and len(set(current_iter)) == 1, (
+                "Cannot create md-variable from variables at different iterates."
+            )
+            assert len(set(names)) == 1, (
+                "Cannot create md-variable from variables with different names."
+            )
+            assert len(set(domains)) == len(domains), (
+                "Cannot create md-variable from variables with overlapping domains."
+            )
         # Default values for empty md variable
         else:
             time_indices = [-1]
diff --git a/src/porepy/numerics/discretization.py b/src/porepy/numerics/discretization.py
index f41287e47b..9e8ae89ec9 100644
--- a/src/porepy/numerics/discretization.py
+++ b/src/porepy/numerics/discretization.py
@@ -1,4 +1,4 @@
-""" Module contains the abstract superclass for all discretizations."""
+"""Module contains the abstract superclass for all discretizations."""
 
 import abc
 from typing import Dict, Union
diff --git a/src/porepy/numerics/fracture_deformation/propagate_fracture.py b/src/porepy/numerics/fracture_deformation/propagate_fracture.py
index 6278e5135f..baae01bac2 100644
--- a/src/porepy/numerics/fracture_deformation/propagate_fracture.py
+++ b/src/porepy/numerics/fracture_deformation/propagate_fracture.py
@@ -624,9 +624,9 @@ def _update_connectivity_fracture_grid(
             sd_primary.tags["domain_boundary_nodes"][fi], axis=0
         )
         sd_secondary.tags["tip_faces"][new_face_indices_l] = ~domain_boundary_faces
-        sd_secondary.tags["domain_boundary_faces"][
-            new_face_indices_l
-        ] = domain_boundary_faces
+        sd_secondary.tags["domain_boundary_faces"][new_face_indices_l] = (
+            domain_boundary_faces
+        )
 
         # Expand array of face-nodes in sd_secondary
         all_faces_l = np.append(all_faces_l, faces_l[:, ~exist], axis=1)
diff --git a/src/porepy/numerics/fv/biot.py b/src/porepy/numerics/fv/biot.py
index ff1e12bdf9..3797647a5e 100644
--- a/src/porepy/numerics/fv/biot.py
+++ b/src/porepy/numerics/fv/biot.py
@@ -1,4 +1,4 @@
-""" Modules contains discretization of poro-elasticity by the multi-point stress
+"""Modules contains discretization of poro-elasticity by the multi-point stress
 approximation.
 
 The discretization scheme is described in
diff --git a/src/porepy/numerics/fv/fvutils.py b/src/porepy/numerics/fv/fvutils.py
index 3ece21e5d0..e1c771fcd0 100644
--- a/src/porepy/numerics/fv/fvutils.py
+++ b/src/porepy/numerics/fv/fvutils.py
@@ -353,7 +353,7 @@ def find_active_indices(
 
 
 def parse_partition_arguments(
-    partition_arguments: Optional[dict[str, int]] = None
+    partition_arguments: Optional[dict[str, int]] = None,
 ) -> tuple[int | None, int | None]:
     """Parse arguments related to the splitting of discretization into subproblems.
 
@@ -477,9 +477,13 @@ def subproblems(
         raise ValueError("Either max_memory or num_subproblems must be given")
 
     if num_part == 1:
-        yield sd, np.arange(sd.num_faces), np.arange(sd.num_cells), np.arange(
-            sd.num_cells
-        ), np.arange(sd.num_faces)
+        yield (
+            sd,
+            np.arange(sd.num_faces),
+            np.arange(sd.num_cells),
+            np.arange(sd.num_cells),
+            np.arange(sd.num_faces),
+        )
 
     else:
         # Since MPxA discretizations are based on interaction regions (cells in the dual
diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py
index b3c977ddb3..7c19d68752 100644
--- a/src/porepy/numerics/fv/mpfa.py
+++ b/src/porepy/numerics/fv/mpfa.py
@@ -1,6 +1,4 @@
-"""Implementation of the multi-point flux approximation O-method.
-
-"""
+"""Implementation of the multi-point flux approximation O-method."""
 
 from __future__ import annotations
 
@@ -1553,7 +1551,7 @@ def _create_bound_rhs(
         elif num_bound == 0:  # all of them are empty
             neu_rob_dir_ind = neu_rob_ind
         else:
-            raise ValueError("Boundary values should be either Dirichlet or " "Neumann")
+            raise ValueError("Boundary values should be either Dirichlet or Neumann")
 
         num_subfno = subcell_topology.num_subfno_unique
 
diff --git a/src/porepy/numerics/fv/tpsa.py b/src/porepy/numerics/fv/tpsa.py
index c6be13ccee..1c5a7a4e01 100644
--- a/src/porepy/numerics/fv/tpsa.py
+++ b/src/porepy/numerics/fv/tpsa.py
@@ -259,7 +259,6 @@ class Tpsa:
     """
 
     def __init__(self, keyword: str) -> None:
-
         self.keyword: str = keyword
         """Keyword used to identify the parameter dictionary."""
 
@@ -677,9 +676,7 @@ def discretize(self, sd: Grid, data: dict) -> None:
                 )
                 + t_shear_rob  # Contribution from Robin boundary conditions.
             )
-        ).reshape(
-            (nd, nf), order="F"
-        )  #
+        ).reshape((nd, nf), order="F")
 
         # Discretize the stress-displacement relation.
         stress, bound_stress = self._vector_laplace_matrices(
diff --git a/src/porepy/numerics/linalg/__init__.py b/src/porepy/numerics/linalg/__init__.py
index a9a2076fa1..3bc99d885b 100644
--- a/src/porepy/numerics/linalg/__init__.py
+++ b/src/porepy/numerics/linalg/__init__.py
@@ -1,2 +1 @@
-"""Directory is intended to contain various linear algebra modules.
-"""
+"""Directory is intended to contain various linear algebra modules."""
diff --git a/src/porepy/numerics/linalg/matrix_operations.py b/src/porepy/numerics/linalg/matrix_operations.py
index 97d979ffaa..f972da8839 100644
--- a/src/porepy/numerics/linalg/matrix_operations.py
+++ b/src/porepy/numerics/linalg/matrix_operations.py
@@ -659,7 +659,6 @@ def invert_diagonal_blocks_numba(a: sps.csr_matrix, size: np.ndarray) -> np.ndar
             parallel=True,
         )
         def inv_compiled_function(is_csr_q, data, indices, indptr, sz):
-
             # Construction of simple data structures (low complexity)
             # Indices for block positions, flattened inverse block positions and nonzeros
             # Expanded block positions
diff --git a/src/porepy/numerics/vem/__init__.py b/src/porepy/numerics/vem/__init__.py
index 3c4716afa1..cc97d1ba15 100644
--- a/src/porepy/numerics/vem/__init__.py
+++ b/src/porepy/numerics/vem/__init__.py
@@ -1,2 +1 @@
-""" Implementation of methods related to the virtual element method.
-"""
+"""Implementation of methods related to the virtual element method."""
diff --git a/src/porepy/numerics/vem/hybrid.py b/src/porepy/numerics/vem/hybrid.py
index a057599d39..c727ce8583 100644
--- a/src/porepy/numerics/vem/hybrid.py
+++ b/src/porepy/numerics/vem/hybrid.py
@@ -3,6 +3,7 @@
 
 @author: fumagalli, alessio
 """
+
 import numpy as np
 import scipy.sparse as sps
 from numpy.linalg import solve
diff --git a/src/porepy/params/data.py b/src/porepy/params/data.py
index de839c3e13..ab7a986155 100644
--- a/src/porepy/params/data.py
+++ b/src/porepy/params/data.py
@@ -1,4 +1,4 @@
-r""" Contains class for storing data / parameters associated with a grid.
+r"""Contains class for storing data / parameters associated with a grid.
 
 At present, the Parameters class is a simple wrapper around a dictionary.
 
diff --git a/src/porepy/params/parameter_dictionaries.py b/src/porepy/params/parameter_dictionaries.py
index 700f8e6c95..36114239d2 100644
--- a/src/porepy/params/parameter_dictionaries.py
+++ b/src/porepy/params/parameter_dictionaries.py
@@ -1,4 +1,4 @@
-""" Parameter dictionaries.
+"""Parameter dictionaries.
 
 Here, we store various parameter dictionaries with "sensible" (typically unitary or
 zero) default values for the parameters required by the discretization objects.
diff --git a/src/porepy/params/tensor.py b/src/porepy/params/tensor.py
index f7895ec9a8..6ced3143f9 100644
--- a/src/porepy/params/tensor.py
+++ b/src/porepy/params/tensor.py
@@ -46,8 +46,7 @@ def __init__(
 
         if np.any(kxx < 0):
             raise ValueError(
-                "Tensor is not positive definite because of "
-                "components in x-direction"
+                "Tensor is not positive definite because of components in x-direction"
             )
 
         perm[0, 0, ::] = kxx
@@ -64,8 +63,7 @@ def __init__(
         # Onsager's principle - tensor should be positive definite
         if np.any((kxx * kyy - kxy * kxy) < 0):
             raise ValueError(
-                "Tensor is not positive definite because of "
-                "components in y-direction"
+                "Tensor is not positive definite because of components in y-direction"
             )
 
         perm[1, 0, ::] = kxy
@@ -92,8 +90,7 @@ def __init__(
             < 0
         ):
             raise ValueError(
-                "Tensor is not positive definite because of "
-                "components in z-direction"
+                "Tensor is not positive definite because of components in z-direction"
             )
 
         perm[2, 0, ::] = kxz
@@ -268,9 +265,9 @@ def __str__(self) -> str:
         # below acknowledges the full fourth-order nature of this tensor (as opposed to
         # the storage format in self.values, which joins two and two dimensions).
         if self.values.shape[:2] == (4, 4):
-            s += f"Each cell has a tensor of shape (2, 2, 2, 2)"
+            s += "Each cell has a tensor of shape (2, 2, 2, 2)"
         elif self.values.shape[:2] == (9, 9):
-            s += f"Each cell has a tensor of shape (3, 3, 3, 3)"
+            s += "Each cell has a tensor of shape (3, 3, 3, 3)"
         else:
             # In EK's understanding, this should never happen. Give a fair description
             # of the situation, and hope the user knows what is going on.
diff --git a/src/porepy/utils/graph.py b/src/porepy/utils/graph.py
index d046a8876b..bb5ebcd61a 100644
--- a/src/porepy/utils/graph.py
+++ b/src/porepy/utils/graph.py
@@ -47,7 +47,7 @@ def color_nodes(self):
             else:
                 return
         raise RuntimeWarning(
-            "number of regions can not be greater than " "number of nodes"
+            "number of regions can not be greater than number of nodes"
         )
 
     def bfs(self, start, color):
diff --git a/src/porepy/utils/grid_utils.py b/src/porepy/utils/grid_utils.py
index c8aeb47a10..5d33bde320 100644
--- a/src/porepy/utils/grid_utils.py
+++ b/src/porepy/utils/grid_utils.py
@@ -1,5 +1,4 @@
-"""Module contains various utility functions for working with grids.
-"""
+"""Module contains various utility functions for working with grids."""
 
 import numpy as np
 import scipy.sparse as sps
diff --git a/src/porepy/utils/interpolation_tables.py b/src/porepy/utils/interpolation_tables.py
index 1a1839a101..670ee3f60a 100644
--- a/src/porepy/utils/interpolation_tables.py
+++ b/src/porepy/utils/interpolation_tables.py
@@ -1,4 +1,4 @@
-""" The module contains interpolation tables, intended for use in function
+"""The module contains interpolation tables, intended for use in function
 evaluations. Specifically, the motivation is to facilitate the parametrization
 framework described in
 
diff --git a/src/porepy/utils/mcolon.py b/src/porepy/utils/mcolon.py
index 6f894261b4..7ba40f7ac8 100644
--- a/src/porepy/utils/mcolon.py
+++ b/src/porepy/utils/mcolon.py
@@ -1,4 +1,4 @@
-""" Efficient numpy.arange for arrays of start and end indices.
+"""Efficient numpy.arange for arrays of start and end indices.
 
 Acknowledgements:
     The functions are a python translation of the corresponding matlab
@@ -43,7 +43,7 @@ def mcolon(lo, hi):
         hi = hi * np.ones(lo.size, dtype=np.int64)
     if lo.size != hi.size:
         raise ValueError(
-            "Low and high should have same number of elements, " "or a single item "
+            "Low and high should have same number of elements, or a single item "
         )
 
     i = hi >= lo + 1
diff --git a/src/porepy/utils/permutations.py b/src/porepy/utils/permutations.py
index a6b40f39d5..98a3847545 100644
--- a/src/porepy/utils/permutations.py
+++ b/src/porepy/utils/permutations.py
@@ -1,5 +1,4 @@
-""" Utility function for permutation of numbers.
-"""
+"""Utility function for permutation of numbers."""
 
 
 def multinary_permutations(base, length):
diff --git a/src/porepy/utils/setmembership.py b/src/porepy/utils/setmembership.py
index 4abe903102..d263d944b8 100644
--- a/src/porepy/utils/setmembership.py
+++ b/src/porepy/utils/setmembership.py
@@ -10,7 +10,9 @@
 from scipy.spatial import KDTree
 
 
-def unique_rows(data: np.ndarray[Any, np.dtype[np.float64]]) -> Tuple[
+def unique_rows(
+    data: np.ndarray[Any, np.dtype[np.float64]],
+) -> Tuple[
     np.ndarray[Any, np.dtype[np.float64]],
     np.ndarray[Any, np.dtype[np.int64]],
     np.ndarray[Any, np.dtype[np.int64]],
diff --git a/src/porepy/utils/sort_points.py b/src/porepy/utils/sort_points.py
index 504d0a10ac..273538dfce 100644
--- a/src/porepy/utils/sort_points.py
+++ b/src/porepy/utils/sort_points.py
@@ -1,5 +1,4 @@
-""" Functions to sort points and edges belonging to geometric objects.
-"""
+"""Functions to sort points and edges belonging to geometric objects."""
 
 from __future__ import annotations
 
diff --git a/src/porepy/viz/diagnostics_mixin.py b/src/porepy/viz/diagnostics_mixin.py
index bf46470dad..e296441323 100644
--- a/src/porepy/viz/diagnostics_mixin.py
+++ b/src/porepy/viz/diagnostics_mixin.py
@@ -1,6 +1,4 @@
-"""Module for diagnostics of PorePy's models built on seaborn.
-
-"""
+"""Module for diagnostics of PorePy's models built on seaborn."""
 
 from __future__ import annotations