Merge branch 'develop' into maint-omit-redundant-residual-computations

.github/workflows/check_tutorials.yml

@@ -1,5 +1,5 @@
-# Run unit, integration, and functional tests.
-name: Pytest
+# Run checks on the tutorials
+name: Check tutorials
 
 # Controls when the action will run. Triggers the workflow on push or pull request
 # events for the main and develop branches
@@ -29,8 +29,8 @@ jobs:
     # https://help.github.com/en/actions/language-and-framework-guides/using-python-with-github-actions#specifying-a-python-version
     strategy:
       matrix:
-        # Only check the latest version of 
-        python-version:  ["3.12"]
+        # Run the test on the latest supported version of Python
+        python-version:  ["3.13"]
       # Complete all versions in matrix even if one fails.
       fail-fast: false
 
@@ -62,9 +62,7 @@ jobs:
     - name: Install external libraries
       run: |
         # Various packages that must be installed 
-        # Install libffi v7. This was (seemed to be) necessary in Nov 2022. The reason could be that
-        # ubuntu-latest was updated, libffi was removed, and scipy complained.
-        sudo apt-get install libglu1-mesa libgeos-dev libffi7
+        sudo apt-get install libglu1-mesa libgeos-dev libffi-dev
         export LD_LIBRARY_PATH=/usr/local/lib64/:$LD_LIBRARY_PATH
 
     - name: Install requirements

.github/workflows/run-pytest-all.yml

@@ -21,7 +21,7 @@ jobs:
     # https://help.github.com/en/actions/language-and-framework-guides/using-python-with-github-actions#specifying-a-python-version
     strategy:
       matrix:
-        python-version:  ["3.12"]  # run only the last python version
+        python-version:  ["3.13"]  # run only the last python version
       # Complete all versions in matrix even if one fails.
       fail-fast: false
 
@@ -53,9 +53,7 @@ jobs:
     - name: Install external libraries
       run: |
         # Various packages that must be installed 
-        # Install libffi v7. This was (seemed to be) necessary in Nov 2022. The reason could be that
-        # ubuntu-latest was updated, libffi was removed, and scipy complained.
-        sudo apt-get install libglu1-mesa libgeos-dev libffi7
+        sudo apt-get install libglu1-mesa libgeos-dev libffi-dev
         export LD_LIBRARY_PATH=/usr/local/lib64/:$LD_LIBRARY_PATH
 
     - name: Install requirements

.github/workflows/run-pytest.yml

@@ -29,7 +29,7 @@ jobs:
     # https://help.github.com/en/actions/language-and-framework-guides/using-python-with-github-actions#specifying-a-python-version
     strategy:
       matrix:
-        python-version:  ["3.11", "3.12"]
+        python-version:  ["3.11", "3.12", "3.13"]
       # Complete all versions in matrix even if one fails.
       fail-fast: false
 
@@ -61,9 +61,7 @@ jobs:
     - name: Install external libraries
       run: |
         # Various packages that must be installed 
-        # Install libffi v7. This was (seemed to be) necessary in Nov 2022. The reason could be that
-        # ubuntu-latest was updated, libffi was removed, and scipy complained.
-        sudo apt-get install libglu1-mesa libgeos-dev libffi7
+        sudo apt-get install libglu1-mesa libgeos-dev libffi-dev
         export LD_LIBRARY_PATH=/usr/local/lib64/:$LD_LIBRARY_PATH
 
     - name: Install requirements

.github/workflows/run-static-checks.yml

@@ -48,7 +48,7 @@ jobs:
     strategy:
       matrix:
         # Run the tests only on the latest supported Python version.
-        python-version:  ["3.12"]
+        python-version:  ["3.13"]
       # Complete all versions in matrix even if one fails.
       fail-fast: false
 

dockerfiles/base/Dockerfile

@@ -6,7 +6,7 @@
 
 # Base the Docker image on the official Python image.
 # For more information, please refer to https://aka.ms/vscode-docker-python
-FROM python:3.12-slim
+FROM python:3.13-slim
 
 
 # Keeps Python from generating .pyc files in the container

src/porepy/__init__.py

@@ -226,12 +226,6 @@
 from porepy.fracs import meshing, fracture_importer
 from porepy.grids import coarsening, partition, refinement
 from porepy.numerics import displacement_correlation
-from porepy.utils.default_domains import (
-    CubeDomain,
-    SquareDomain,
-    UnitSquareDomain,
-    UnitCubeDomain,
-)
 
 # Applications
 from porepy.applications.md_grids import (

src/porepy/applications/profiling/run_profiling.py

@@ -198,7 +198,7 @@ def run_model_with_tracer(args, model) -> None:
             refinement.
             - min_duration (int): Minimum duration in microseconds for a function to be
             recorded by VizTracer.
-            - keep_output (bool): Whether to keep the output file after viewing it.
+            - clear_output (bool): Whether to clear the output file after viewing it.
         model: The model to be run and profiled.
 
     Raises:
@@ -249,7 +249,7 @@ def run_model_with_tracer(args, model) -> None:
     results_path = pathlib.Path(__file__).parent / save_file
     tracer.save(str(results_path))
     subprocess.run(["vizviewer", "--port", "9002", results_path])
-    if not args.keep_output:
+    if args.clear_output:
         results_path.unlink()
 
 
@@ -288,10 +288,10 @@ def run_model_with_tracer(args, model) -> None:
         + " named after the chosen physics, geometry, and grid refinement.",
     )
     parser.add_argument(
-        "--keep_output",
+        "--clear_output",
         action="store_true",
-        default=True,
-        help="Keep viztracer output after running.",
+        help="Clear viztracer output after running. By default, the log files are"
+        + " kept.",
     )
     parser.add_argument(
         "--min_duration",

src/porepy/applications/test_utils/benchmarks.py

@@ -108,7 +108,7 @@ def effective_normal_permeability(
         projection = pp.ad.MortarProjections(self.mdg, subdomains, interfaces, dim=1)
 
         normal_gradient = pp.ad.Scalar(2) * (
-            projection.secondary_to_mortar_avg
+            projection.secondary_to_mortar_avg()
             @ self.aperture(subdomains) ** pp.ad.Scalar(-1)
         )
 

src/porepy/examples/mandel_biot.py

@@ -23,6 +23,7 @@
 
 from __future__ import annotations
 
+import warnings
 from dataclasses import dataclass
 from typing import Callable, Literal, Union, cast
 
@@ -1022,6 +1023,18 @@ def _plot_horizontal_flux(self, color_map: mcolors.ListedColormap) -> None:
             color_map: listed color map object.
 
         """
+        if self.grid_type() == "cartesian":
+            # A division by zero was discovered while running the Mandel-problem on a
+            # Cartesian grid. The issue was related to plotting, and the current
+            # solution is to skip the troublesome plotting if the grid type is
+            # Cartesian. The GitHub-issue for this is found at:
+            # https://github.com/pmgbergen/porepy/issues/1137
+            warnings.warn(
+                """Division by x-component of the normal vector for internal faces of
+                cells adjacent to the southern boundary causes division by zero for
+                cartesian grids."""
+            )
+            return
         sd = self.mdg.subdomains()[0]
         xf = sd.face_centers[0]
         nx = sd.face_normals[0]

src/porepy/fracs/fracture_network_3d.py

@@ -1460,8 +1460,8 @@ def _points_2_plane(
         # local numbering of points
         edges_2d = np.empty_like(edges_loc)
         for ei in range(edges_loc.shape[1]):
-            edges_2d[0, ei] = np.argwhere(p_ind_loc == edges_loc[0, ei])
-            edges_2d[1, ei] = np.argwhere(p_ind_loc == edges_loc[1, ei])
+            edges_2d[0, ei] = np.argmax(p_ind_loc == edges_loc[0, ei])
+            edges_2d[1, ei] = np.argmax(p_ind_loc == edges_loc[1, ei])
 
         assert edges_2d[:2].max() < p_loc.shape[1]
 
@@ -1662,7 +1662,8 @@ def impose_external_boundary(
                 # points should not be sorted. Splitting of non-convex fractures into
                 # convex subparts is handled below.
                 new_frac = pp.PlaneFracture(constrained_polys[sub_i], sort_points=False)
-                self.add(new_frac)
+                new_frac.index = len(self.fractures)
+                self.fractures.append(new_frac)
                 ind_map = np.hstack((ind_map, fi))
 
         # Now, we may modify the fractures for two reasons:
@@ -1989,7 +1990,9 @@ def new_pairs_and_angles(ei):
                 # thus a linear ordering should be fine also for a subpolygon.
                 verts = np.unique(triangles[tris])
                 # To be sure, check the convexity of the polygon.
-                self.add(pp.PlaneFracture(f.pts[:, verts], check_convexity=False))
+                new_frac = pp.PlaneFracture(f.pts[:, verts], check_convexity=False)
+                new_frac.index = len(self.fractures)
+                self.fractures.append(new_frac)
                 ind_map = np.hstack((ind_map, fi))
 
             # Finally, increase pointer to ind_map array
@@ -2011,7 +2014,9 @@ def new_pairs_and_angles(ei):
         boundary_tags = self.tags.get("boundary", [False] * len(self.fractures))
         if keep_box:
             for pnt in self.domain.polytope:
-                self.add(pp.PlaneFracture(pnt))
+                new_frac = pp.PlaneFracture(pnt)
+                new_frac.index = len(self.fractures)
+                self.fractures.append(new_frac)
                 boundary_tags.append(True)
         self.tags["boundary"] = boundary_tags
 

src/porepy/geometry/constrain_geometry.py

@@ -48,12 +48,7 @@ def lines_by_polygon(
 
     """
     import shapely.geometry as shapely_geometry
-    import shapely.speedups as shapely_speedups
 
-    try:
-        shapely_speedups.enable()
-    except AttributeError:
-        pass
     # it stores the points after the intersection
     int_pts = np.empty((2, 0))
     # define the polygon

src/porepy/geometry/intersections.py

@@ -1729,12 +1729,6 @@ def triangulations(
     # to other cell shapes. This would require more general data structures, but should
     # not be too much of an effort.
     import shapely.geometry as shapely_geometry
-    import shapely.speedups as shapely_speedups
-
-    try:
-        shapely_speedups.enable()
-    except AttributeError:
-        pass
 
     n_1 = t_1.shape[1]
     n_2 = t_2.shape[1]
@@ -1893,13 +1887,6 @@ def surface_tessellations(
 
     # local imports
     import shapely.geometry as shapely_geometry
-    import shapely.speedups as shapely_speedups
-
-    try:
-        shapely_speedups.enable()
-    except AttributeError:
-        # Nothing to do here, but this may be slow.
-        pass
 
     def _min_max_coord(coord):
         # Convenience function to get max and minimum coordinates for a set of polygons

src/porepy/grids/simplex.py

@@ -12,7 +12,7 @@
 
 import numpy as np
 import scipy.sparse as sps
-import scipy.spatial.qhull
+import scipy.spatial
 
 from porepy.grids.grid import Grid
 from porepy.utils import accumarray, setmembership
@@ -55,7 +55,7 @@ def __init__(
         self.dim = 2
 
         if tri is None:
-            triangulation = scipy.spatial.qhull.Delaunay(p.transpose())
+            triangulation = scipy.spatial.Delaunay(p.transpose())
             tri = triangulation.simplices
             tri = tri.transpose()
 
@@ -240,7 +240,7 @@ def __init__(
         # Transform points to column vector if necessary (scipy.Delaunay requires this
         # format)
         if tet is None:
-            tesselation = scipy.spatial.qhull.Delaunay(p.transpose())
+            tesselation = scipy.spatial.Delaunay(p.transpose())
             tet = tesselation.simplices.transpose()
 
         if name is None: