diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 90a97d67..22c82ab8 100644 --- a/.github/workflows/ci.yml +++ b/.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] diff --git a/src/darsia/measure/wasserstein.py b/src/darsia/measure/wasserstein.py index dd1f2f20..63bf6753 100644 --- a/src/darsia/measure/wasserstein.py +++ b/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.") diff --git a/src/darsia/utils/fv.py b/src/darsia/utils/fv.py index 383e233c..ce00e790 100644 --- a/src/darsia/utils/fv.py +++ b/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,7 +228,9 @@ 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: @@ -228,7 +238,7 @@ def __call__(self, normal_flux: np.ndarray, concatenate: bool = True) -> np.ndar 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,15 +365,39 @@ 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] @@ -371,8 +405,12 @@ def cell_to_face_average( 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": diff --git a/tests/unit/test_fv.py b/tests/unit/test_fv.py index 3f5f8e29..c1d0d04d 100644 --- a/tests/unit/test_fv.py +++ b/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():