Merge branch 'main' into petsc4py_ksp

.github/workflows/ci.yml

@@ -39,6 +39,50 @@ jobs:
         python-version: ${{ matrix.python-version}}
         architecture: x64
 
+    - name: Install dependencies
+      run: |
+        sudo apt-get update
+        sudo apt-get install -y libhypre-dev libmumps-seq-dev
+        pip install numpy mpi4py
+        pip install petsc petsc4py
+
+    - name: Set environment variables
+      run: |
+        echo "PETSC_DIR=/usr/local/petsc" >> $GITHUB_ENV
+        echo "PETSC_ARCH=arch-linux-c-opt" >> $GITHUB_ENV
+
+    # Installation of petsc in parallel
+
+    # - name: Install dependencies
+    #   run: |
+    #     sudo apt-get update
+    #     sudo apt-get install -y mpich
+    #     pip install numpy mpi4py
+    #     pip install petsc petsc4py
+
+    # - name: Set environment variables
+    #   run: |
+    #     echo "PETSC_DIR=/usr/local/petsc" >> $GITHUB_ENV
+    #     echo "PETSC_ARCH=arch-linux-c-opt" >> $GITHUB_ENV
+
+    # - name: Install Hypre
+    #   run: |
+    #     git clone https://github.com/hypre-space/hypre.git
+    #     cd hypre/src
+    #     ./configure --with-MPI
+    #     make
+    #     sudo make install
+
+    # - name: Install MUMPS
+    #   run: |
+    #     git clone https://github.com/scivision/mumps.git
+    #     cd mumps
+    #     mkdir build
+    #     cd build
+    #     cmake -G "Ninja" -DMUMPS_parallel=on ..
+    #     cmake --build .
+    #     sudo make install
+
     - name: Install DarSIA
       run: pip install .[dev]
 

src/darsia/measure/wasserstein.py

@@ -987,7 +987,7 @@ def _compute_face_weight(self, flat_flux: np.ndarray) -> np.ndarray:
 
             # Combine weights**2 / |weight * flux| on faces
             face_weights = harm_avg_face_weights**2 / norm_weighted_face_flux
-            face_weights_inv = 1.0 / face_weights
+            face_weights_inv = norm_weighted_face_flux / harm_avg_face_weights**2
         else:
             raise ValueError(f"Mobility mode {self.mobility_mode} not supported.")
 

src/darsia/utils/fv.py

@@ -63,9 +63,17 @@ def __init__(
             )
         elif mode == "faces":
             if lumping:
-                mass_matrix = 0.5 * sps.diags(
-                    np.prod(grid.voxel_size) * np.ones(grid.num_faces, dtype=float)
-                )
+                # For a Cartesian grid, the lumped mass matrix is diagonal and
+                # contains the area half of the cell to left and right of the face on the
+                # diagonal. This origins from integration of RT0 basis functions over the
+                # cells connected to the face. The mass matrix is lumped, i.e. the
+                # diagonal entries are the sum of the off-diagonal entries.
+                # After all, this is identical to the arithmetic avg of cell volumes
+                # connected to the face. This requires careful handling of the boundary faces.
+                # NOTE: This assumes equidistant grid spacing.
+                # NOTE: Fluxes on boundary are 0 and excluded.
+                volume_scaling = np.ones(grid.num_faces, dtype=float)
+                mass_matrix = sps.diags(np.prod(grid.voxel_size) * volume_scaling)
             else:
                 raise NotImplementedError
 
@@ -220,15 +228,17 @@ def __init__(self, grid: darsia.Grid) -> None:
         self.num_tangential_directions = grid.dim - 1
         self.grid = grid
 
-    def __call__(self, normal_flux: np.ndarray, concatenate: bool = True) -> np.ndarray:
+    def __call__(
+        self, normal_flux: np.ndarray, concatenate: bool = True
+    ) -> list[np.ndarray] | np.ndarray:
         """Apply the operator to the normal fluxes.
 
         Args:
             normal_flux (np.ndarray): normal fluxes
             concatenate (bool, optional): whether to concatenate the tangential fluxes
 
         Returns:
-            np.ndarray: tangential fluxes
+            np.ndarray or list of arrays: tangential fluxes
 
         """
         # Apply the operator to the normal fluxes
@@ -355,24 +365,52 @@ def cell_to_face_average(
         np.ndarray: face-based quantity
 
     """
-    # Collect cell quantities for each face
+
+    # Prepare arrays
     neighbouring_cell_values = np.zeros((grid.num_faces, 2), dtype=float)
     face_qty = np.zeros(grid.num_faces, dtype=float)
 
-    # Flatten cell quantities
-    flat_cell_qty = cell_qty.ravel("F")
+    # Apply normal projection of cell values onto face orientations (axis-aligned)
+    if len(cell_qty.shape) == grid.dim or (
+        len(cell_qty.shape) == grid.dim + 1 and cell_qty.shape[-1] == 1
+    ):
+
+        # Flatten cell quantities
+        single_flat_cell_qty = cell_qty.ravel("F")
+        flat_cell_qty = [single_flat_cell_qty for _ in range(grid.dim)]
+
+    elif len(cell_qty.shape) == grid.dim + 1 and cell_qty.shape[-1] == grid.dim:
+
+        # Flatten each component of the cell quantities
+        flat_cell_qty = [cell_qty[..., i].ravel("F") for i in range(grid.dim)]
+
+    elif len(cell_qty.shape) == grid.dim + 2 and cell_qty.shape[-2:] == (
+        grid.dim,
+        grid.dim,
+    ):
+        # Flatten each diagonal component of the cell quantities (from normal projection)
+        flat_cell_qty = [cell_qty[..., i, i].ravel("F") for i in range(grid.dim)]
 
-    # Iterate over normal directions
+    else:
+
+        raise NotImplementedError("Dimension not supported.")
+
+    # Collect neighbouring cell quantities
     for orientation in range(grid.dim):
+
         # Fetch faces
         faces = grid.faces[orientation]
 
         # Fetch neighbouring cells
         neighbouring_cells = grid.connectivity[faces]
 
         # Fetch cell quantities
-        neighbouring_cell_values[faces, 0] = flat_cell_qty[neighbouring_cells[:, 0]]
-        neighbouring_cell_values[faces, 1] = flat_cell_qty[neighbouring_cells[:, 1]]
+        neighbouring_cell_values[faces, 0] = flat_cell_qty[orientation][
+            neighbouring_cells[:, 0]
+        ]
+        neighbouring_cell_values[faces, 1] = flat_cell_qty[orientation][
+            neighbouring_cells[:, 1]
+        ]
 
     # Perform averaging
     if mode == "arithmetic":

tests/unit/test_fv.py

@@ -101,7 +101,7 @@ def test_mass_face_2d():
 
     # Check diagonal values
     assert len(np.unique(np.diag(mass))) == 1
-    assert np.isclose(mass[0, 0], 0.5 * 0.125)
+    assert all([np.isclose(mass[i, i], 0.5 * 0.25) for i in range(20)])
 
 
 def test_mass_face_3d():
@@ -116,7 +116,7 @@ def test_mass_face_3d():
 
     # Check diagonal values
     assert len(np.unique(np.diag(mass))) == 1
-    assert np.isclose(mass[0, 0], 0.5 * 0.25)
+    assert all([np.isclose(mass[i, i], 0.5 * 0.25 * 2) for i in range(60)])
 
 
 def test_tangential_reconstruction_2d_1():