From 4ea2b345d7dda1f2a599e0d693af5d4adc6dbbb5 Mon Sep 17 00:00:00 2001 From: Jakub Both Date: Tue, 27 Aug 2024 00:24:53 +0200 Subject: [PATCH 1/8] BUG: Fix lumped mass matrix for faces. --- src/darsia/utils/fv.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/darsia/utils/fv.py b/src/darsia/utils/fv.py index 383e233c..cc863d2b 100644 --- a/src/darsia/utils/fv.py +++ b/src/darsia/utils/fv.py @@ -63,8 +63,20 @@ 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. + volume_scaling = np.ones(grid.num_faces, dtype=float) + volume_scaling[grid.exterior_faces] = 0.5 + mass_matrix = sps.diags( + np.prod(grid.voxel_size) + * volume_scaling + * np.ones(grid.num_faces, dtype=float) ) else: raise NotImplementedError From dcd30006d1e6c899342a66a70af974992d433d7e Mon Sep 17 00:00:00 2001 From: Jakub Both Date: Tue, 27 Aug 2024 01:06:09 +0200 Subject: [PATCH 2/8] BUG: Previous idea was wrong - boundary faces are 0 anyhow. --- src/darsia/utils/fv.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/darsia/utils/fv.py b/src/darsia/utils/fv.py index cc863d2b..d0fbdd72 100644 --- a/src/darsia/utils/fv.py +++ b/src/darsia/utils/fv.py @@ -71,13 +71,9 @@ def __init__( # 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) - volume_scaling[grid.exterior_faces] = 0.5 - mass_matrix = sps.diags( - np.prod(grid.voxel_size) - * volume_scaling - * np.ones(grid.num_faces, dtype=float) - ) + mass_matrix = sps.diags(np.prod(grid.voxel_size) * volume_scaling) else: raise NotImplementedError From 4db0bb75081dbfefbe6eed70ce00f3a57ae1463b Mon Sep 17 00:00:00 2001 From: Jakub Both Date: Tue, 27 Aug 2024 01:19:34 +0200 Subject: [PATCH 3/8] TST: Adapt tests for fixed bugs. --- tests/unit/test_fv.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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(): From dc45e43985b21d8053fe3c8db6a339e00cd0ce11 Mon Sep 17 00:00:00 2001 From: Jakub Both Date: Tue, 27 Aug 2024 09:13:09 +0200 Subject: [PATCH 4/8] ENH: Extend cell to face capabilities for tensor quantities. --- src/darsia/utils/fv.py | 46 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 38 insertions(+), 8 deletions(-) diff --git a/src/darsia/utils/fv.py b/src/darsia/utils/fv.py index d0fbdd72..ce00e790 100644 --- a/src/darsia/utils/fv.py +++ b/src/darsia/utils/fv.py @@ -228,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: @@ -236,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 @@ -363,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: - # Iterate over normal directions + # 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)] + + else: + + raise NotImplementedError("Dimension not supported.") + + # Collect neighbouring cell quantities for orientation in range(grid.dim): + # Fetch faces faces = grid.faces[orientation] @@ -379,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": From 1f9981989b1eecd2d4910283caa386cd3236fbe2 Mon Sep 17 00:00:00 2001 From: Jakub Both Date: Tue, 27 Aug 2024 09:16:30 +0200 Subject: [PATCH 5/8] FIX: Use updated variable name. --- src/darsia/measure/wasserstein.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/darsia/measure/wasserstein.py b/src/darsia/measure/wasserstein.py index 07d8a753..1c95e25f 100644 --- a/src/darsia/measure/wasserstein.py +++ b/src/darsia/measure/wasserstein.py @@ -896,7 +896,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.") From cf28721861efe05b650b409435cacc4222015bdf Mon Sep 17 00:00:00 2001 From: Jakub Both Date: Wed, 28 Aug 2024 19:35:06 +0200 Subject: [PATCH 6/8] MAINT: Add petsc4py to build workflow --- .github/workflows/ci.yml | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 90a97d67..a48f0962 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -39,6 +39,13 @@ jobs: python-version: ${{ matrix.python-version}} architecture: x64 + - name: Install dependencies + run: | + sudo apt-get update + sudo apt-get install -y libhypre-dev + pip install numpy mpi4py + pip install petsc petsc4py + - name: Install DarSIA run: pip install .[dev] From 33a3453ffd5ca54fa2e509013db820826c08cc30 Mon Sep 17 00:00:00 2001 From: Jakub Both Date: Wed, 28 Aug 2024 19:56:17 +0200 Subject: [PATCH 7/8] MAINT: Extend build test to parallel versions of mumps and hypre. --- .github/workflows/ci.yml | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a48f0962..d7133e17 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -42,10 +42,33 @@ jobs: - name: Install dependencies run: | sudo apt-get update - sudo apt-get install -y libhypre-dev + 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] From d1b8955c6efe818e82b0969351c497cb3c2ccadc Mon Sep 17 00:00:00 2001 From: Jakub Both Date: Wed, 28 Aug 2024 20:37:05 +0200 Subject: [PATCH 8/8] MAINT: Use serial versions of mumps and hypre --- .github/workflows/ci.yml | 48 ++++++++++++++++++++++++++-------------- 1 file changed, 31 insertions(+), 17 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index d7133e17..22c82ab8 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -42,7 +42,7 @@ jobs: - name: Install dependencies run: | sudo apt-get update - sudo apt-get install -y mpich + sudo apt-get install -y libhypre-dev libmumps-seq-dev pip install numpy mpi4py pip install petsc petsc4py @@ -51,23 +51,37 @@ jobs: 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 + # Installation of petsc in parallel - - 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 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]