diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..493a5f9 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,34 @@ +name: CI + +on: + push: + branches: ["**"] + pull_request: + branches: ["**"] + +jobs: + build-and-test: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v4 + + - name: Install MPI + run: sudo apt-get update && sudo apt-get install -y openmpi-bin libopenmpi-dev + + - name: Build + run: make -C src/ + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.x" + + - name: Install Python dependencies + run: pip install -r requirements-test.txt + + - name: Run unit tests + run: pytest tests/ -v + + - name: Run integration tests + run: pytest tests/ -v --run-integration diff --git a/.gitignore b/.gitignore index bc97f63..d9ab01f 100644 --- a/.gitignore +++ b/.gitignore @@ -1,34 +1,53 @@ -# POREMAPS +# ============================================================ +# POREMAPS — project-specific simulation files +# ============================================================ + +# Binary input geometry files (keep only benchmark examples) *.raw !simulations/99_benchmarks/01_poiseuille/ptube_54_54_50_vs_2e-05.raw !simulations/99_benchmarks/02_channel_flow/rect_channel_24_14_50_vs_5e-05.raw !simulations/99_benchmarks/03_spherepackings/SPHERE3D_bcc_NX100_NY100_NZ100_PHI0.3_DIAMETER0.0008745237084764591_VSIZE1e-05.raw !simulations/99_benchmarks/03_spherepackings/SPHERE3D_fcc_NX100_NY100_NZ100_PHI0.3_DIAMETER0.000693979234383925_VSIZE1e-05.raw + +# Solver input files (keep only benchmark examples and the template) *.inp !simulations/99_benchmarks/01_poiseuille/input_ptube_54_54_50_vs_2e-05.inp !simulations/99_benchmarks/02_channel_flow/input_rect_channel_24_14_50_vs_5e-05.raw !simulations/99_benchmarks/03_spherepackings/input_SPHERE3D_bcc_NX100_NY100_NZ100_PHI0.3_DIAMETER0.0008745237084764591_VSIZE1e-05.inp !simulations/99_benchmarks/03_spherepackings/input_SPHERE3D_fcc_NX100_NY100_NZ100_PHI0.3_DIAMETER0.000693979234383925_VSIZE1e-05.inp !input_template.inp + +# Solver log and binary output files *.log + +# VTK/ParaView output files (generated by fields2vtu.py) +*.vtu +*.pvtu +*.vts +*.pvts + +# Compiled solver binary (legacy name kept for safety) STOKES_SOLVER* -POREMAPS* +bin/POREMAPS +# ============================================================ +# C/C++ build artifacts +# ============================================================ -# Prerequisites +# Dependency files *.d -# Compiled Object files +# Compiled object files *.slo *.lo *.o *.obj -# Precompiled Headers +# Precompiled headers *.gch *.pch -# Compiled Dynamic libraries +# Shared libraries *.so *.dylib *.dll @@ -37,7 +56,7 @@ POREMAPS* *.mod *.smod -# Compiled Static libraries +# Static libraries *.lai *.la *.a @@ -47,3 +66,31 @@ POREMAPS* *.exe *.out *.app + +# ============================================================ +# Python +# ============================================================ +__pycache__/ +*.pyc +*.pyo +*.pyd +.pytest_cache/ +.coverage +htmlcov/ + +# ============================================================ +# Editors and IDEs +# ============================================================ +.vscode/ +.idea/ +*.swp +*.swo +*~ +.clangd/ +compile_commands.json + +# ============================================================ +# Operating system +# ============================================================ +.DS_Store +Thumbs.db diff --git a/LICENSE.md b/LICENSE.md index f06d9bf..63e260f 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -1,4 +1,4 @@ -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in diff --git a/README.md b/README.md index 7caa4d7..4482c04 100644 --- a/README.md +++ b/README.md @@ -1,35 +1,69 @@ [![POREMAPS](doc/poremaps_logo_grey.png)](https://www.mib.uni-stuttgart.de/cont/) -**POREMAPS** is a Finite Difference based PORous Media Anisotropic Permeability Solver for Stokes flow. + +**POREMAPS** is a Finite Difference based PORous Media Anisotropic Permeability Solver for Stokes flow. [![Identifier](https://img.shields.io/badge/doi-10.18419%2Fdarus--3676-d45815)](https://doi.org/10.18419/darus-3676) [![Identifier](https://img.shields.io/badge/Publication-blue)](https://doi.org/10.69631/ipj.v2i1nr39) -## How to build: -We suggest to build yourself a local version of [OpenMPI](https://www.open-mpi.org/). -The code was tested with the following versions: -- Linux: - 1. `openmpi-3.1.5` - 2. `openmpi-4.1.5` - 3. `openmpi-5.0.1` -- Windows - 1. `MS-MPI v10.1.1` - -Open `makefile` and edit `CLINKER` path to your `mpiCC`. The build command creates an executable `POREMAPS` in the `bin` directory. +Given a binarized 3-D pore geometry (`.raw` file), POREMAPS solves the Stokes equations on a Cartesian finite-difference grid using MPI domain decomposition and computes the anisotropic intrinsic permeability tensor of the sample. + +--- + +## Contents + +- [How to build](#how-to-build) +- [How to run](#how-to-run) +- [Testing](#testing) +- [Boundary conditions](#boundary-conditions) +- [Input data](#input-data) +- [Output files](#output-files) +- [Aspects of parallelization](#aspects-of-parallelization) +- [Compute the permeability tensor](#compute-the-permeability-tensor) +- [License](#license) +- [How to cite](#how-to-cite) + +--- + +## How to build + +POREMAPS requires an MPI C++ compiler (`mpiCC`). We recommend building a local version of [OpenMPI](https://www.open-mpi.org/). Tested versions: + +| Platform | MPI implementation | +|---|---| +| Linux | `openmpi-3.1.5`, `openmpi-4.1.5`, `openmpi-5.0.1` | +| Windows 10 | `MS-MPI v10.1.1` | + +**Linux** — if `mpiCC` is on your `PATH` (standard OpenMPI installation), just run: ```shell -linux@fastmachine:~$ cd PATH/TO/POREMAPS/src/ -linux@fastmachine:~$ vim makefile # edit CLINKER -linux@fastmachine:~$ make +cd PATH/TO/POREMAPS/src/ +make ``` -To build it on MS Windows (MS Windows 10), we recommend to create a Visual Studio (v16.8.3) project from the existing code and link to the corresponding MS-MPI (v10.1.1) via the project properties. Tested setup in brackets. + +If your MPI compiler is installed at a non-standard location, pass it on the command line: +```shell +make MPICXX=/path/to/your/mpiCC +``` + +A debug build (no optimisation, debug symbols) is also available: +```shell +make debug +``` + +The executable `POREMAPS` is placed in the `bin/` directory. + +**Windows** — create a Visual Studio (v16.8.3) project from the source files and link MS-MPI (v10.1.1) via the project properties. + +--- ## How to run -Copy the executable `POREMAPS` to the folder with input file and geometry and run by + ```shell -linux@fastmachine:~$ mylocal_mpirun -np 8 POREMAPS my_inputfile +mpirun -np 8 ./bin/POREMAPS my_inputfile.inp ``` -Parameters in input file `input_template.inp`: -```cpp +The input file must follow the exact format shown in `input_template.inp`. All parameters are read in order; the field names are ignored by the parser. + +```text dom_decomposition 0 0 0 boundary_method 0 geometry_file_name ptube_54_54_50_vs_2e-05.raw @@ -45,127 +79,219 @@ porosity 1.0 dom_interest 0 0 0 0 0 0 write_output 1 1 0 0 ``` -Parameters explained: -1. **dom_decomposition**: Number of Ranks in $\mathbf{e}_1-$, $\mathbf{e}_2-$ and $\mathbf{e}_3-$ direction. If MPI should do this by itself set it to `0 0 0`. -2. **boundary_method**: Information see below. -3. **geometry_file_name**: File name of the input geometry (`.raw`-file). -4. **size_x_y_z**: Size of the input geometry in voxel. -5. **voxel_size**: Voxel size in meter. -6. **max_iter**: Maximum number of iterations. -7. **it_eval**: Frequency to evaluate convergence criterion. -8. **it_write**: Frequency to write log file. -9. **log_file_name**: Name of log file. -10. **solving_algorithm**: Algorithm see `solver.cc` for more information. -11. **eps**: Convergence criterion. -12. **porosity**: Porosity of the sample or domain of interest. If porosity equals -1.0 the solver computes it for the whole domain (not the effective porosity). Specify the porosity if, e.g., solid frames are used. -13. **dom_interest**: Start and End of domain of interest for all three directions in voxel. Solver gives the permeability of this subdomain and for whole domain if set. If all entries equal 0 only permeability of the whole domain is computed. -14. **write_output** Defines which `.raw`-file output to be written. Four entries for velocity (all three components), pressure, neighborhood and domain decomposition in that order. If 1, field is written. - - -## Boundary Conditions - -Choose the **boundary_method** in the input file according to the requirements: - -| Value | Description | -|:-----:|-------------| -| 0 | periodic in all three directions | -| 1 | periodic in direction of main pressure gradient, slip condition in the two other directions | -| 2 | periodic in direction of main pressure gradient, no-slip condition in the two other directions | -| 3 | non-periodic in direction of main pressure gradient, slip condition in the two other directions | -| 4 | non-periodic in direction of main pressure gradient, no-slip condition in the two other directions | + +### Parameter reference + +| Parameter | Value(s) | Description | +|---|---|---| +| `dom_decomposition` | `int int int` | Number of MPI ranks in $\mathbf{e}_1$, $\mathbf{e}_2$, $\mathbf{e}_3$. Set to `0 0 0` to let MPI choose automatically. | +| `boundary_method` | `0`–`4` | Physical boundary setup; see [Boundary conditions](#boundary-conditions). | +| `geometry_file_name` | `string` | Path to the input geometry file (`.raw`). | +| `size_x_y_z` | `int int int` | Global domain size in voxels $(n_x,\ n_y,\ n_z)$. | +| `voxel_size` | `float` | Voxel edge length in metres. | +| `max_iter` | `int` | Maximum number of pseudo-time iterations. | +| `it_eval` | `int` | Evaluate the convergence criterion every `it_eval` iterations. | +| `it_write` | `int` | Write a log entry every `it_write` iterations (must be a multiple of `it_eval`). | +| `log_file_name` | `string` | Path for the output log file. | +| `solving_algorithm` | `1`, `2`, or `3` | `1` = Jacobi, `2` = Gauss–Seidel, `3` = SOR. | +| `eps` | `float` | Convergence threshold; iteration stops when the relative permeability change falls below this value. | +| `porosity` | `float` | Porosity of the sample. Set to `-1.0` to have the solver compute it from the geometry (not the effective porosity). Specify explicitly when using solid frames or sub-domain evaluation. | +| `dom_interest` | `x_min x_max y_min y_max z_min z_max` | Voxel index bounds of a sub-domain of interest. The solver reports additional permeability values for this region. Set all entries to `0` to skip sub-domain evaluation. | +| `write_output` | `int int int int` | Controls which field files are written: velocity (`velx/y/z`), pressure, voxel neighbourhood, domain decomposition — in that order. `1` = write, `0` = skip. | + +--- + +## Testing + +A Python test suite covering the core numerical kernels lives in the `tests/` directory. The tests verify the finite-difference stencils, domain decomposition arithmetic, boundary condition assignment, and permeability and convergence formulae via reference implementations that mirror the C++ source without requiring MPI. + +**Dependencies** — `numpy` and `pytest`: +```shell +pip install numpy pytest +``` + +**Run all unit tests** from the repository root: +```shell +pytest tests/ +``` + +**Run with verbose output:** +```shell +pytest tests/ -v +``` + +**Integration tests** invoke the compiled `POREMAPS` binary on small generated geometries. They require `mpirun` on `PATH` and the binary to be built. They are opt-in: +```shell +pytest tests/ --run-integration +``` + +| Module | Functions tested | Source file | +|---|---|---| +| `test_solver.py` | `split_number`, `get_2nd_derivation` | `solver.cc` | +| `test_geometry.py` | `eval_geometry`, `get_dom_limits`, `get_proc_porosity` | `geometry.cc` | +| `test_evaluation.py` | `compute_permeability`, `compute_convergence` | `evaluation.cc` | +| `test_initialize.py` | `determine_comm_pattern`, `set_halos_initial` | `initialize.cc` | +| `test_integration.py` | Full binary run via `mpirun` | — (opt-in) | + +--- + +## Boundary conditions + +The main pressure gradient and flow direction is always $\mathbf{e}_3$ (z). Choose `boundary_method` according to the physical setup: + +| Value | z-direction | x- and y-directions | +|:-----:|---|---| +| `0` | periodic | periodic | +| `1` | periodic | slip walls | +| `2` | periodic | no-slip walls | +| `3` | non-periodic (pressure-driven) | slip walls | +| `4` | non-periodic (pressure-driven) | no-slip walls | + +--- ## Input data -8-bit RAW datasets are unprocessed, binary image files where each pixel's intensity is stored as an 8-bit value (ranging from 0 to 255). In POREMAPS we use binarized data resulting in a dataset consisting of zeros (fluid voxel) and ones (solid voxels). These datasets contain only pixel data without any additional headers or metadata, meaning they lack information about dimensions or image format. As a result, the user must manually specify these parameters when opening or processing the files. The solver expects to read data stored Fortran-like index order. -In `python` we suggest to use libraries such as `numpy` and `matplotlib` or `PIL` to open the file and display slices of the 3D-image. See the `python` script `fields2vtu.py` for an example. The mentioned script also helps to convert the `.raw` data into `.vtu` format to open them with [ParaView](https://www.paraview.org/). -For the input dataset, permeability, meaning at least one available flow path through the material, must be ensured (this can be tested using a labeling algorithm `scipy.ndimage.label`), where only face-to-face connections are considered permeable. Beyond that, POREMAPS does not undertake any changes to the domains. If, for example, mirrored domains are to be used, this must be done in preprocessing by the user. We also recommend using `python`, `numpy` for this. + +POREMAPS reads 8-bit binary (`.raw`) geometry files. Each byte encodes one voxel: `0` = fluid, `1` = solid. There is no file header — dimensions must be provided via `size_x_y_z` in the input file. + +**Byte order:** data must be stored with x varying fastest (Fortran index order). When preparing input files with NumPy, use a `(nz, ny, nx)` shaped array written in C order, which gives x-fastest layout: + +```python +import numpy as np +geom = np.zeros((nz, ny, nx), dtype=np.uint8) # 0 = fluid, 1 = solid +geom.tofile("geometry.raw") # default C order → x varies fastest +``` + +**Connectivity:** the geometry must contain at least one connected flow path from the inlet to the outlet face (face-to-face connectivity). This can be verified with: + +```python +from scipy.ndimage import label +labeled, n = label(geom == 0) # label connected fluid regions +``` + +Preprocessing (mirroring, padding, solid frame addition) must be done before passing the file to POREMAPS. We recommend NumPy for this. + +To visualise `.raw` output fields, use the included `fields2vtu.py` script to convert them to `.vtu` format for [ParaView](https://www.paraview.org/). + +--- ## Output files -The solver creates three files for the final velocity field (`velx_*.raw`, `vely_*.raw`, `velz_*.raw`) ( $\mathbf{e}_1-$, $\mathbf{e}_2-$, $\mathbf{e}_3-$ direction ) and one for the final pressure field (`press_*.raw`). Additionaly, the determined voxel neighborhood cases (`voxel_neighborhood_*.raw`) and domain decomposition (`domain_decomp_*.raw`) are created. All files are in raw image format (`double` for all velocity and pressure files, `int` for neighborhood and domain decomposition). Using the file `fields2vtu.py` allows a direct conversion into a `*.vtu` file which can be visualized, e.g., with [ParaView](https://www.paraview.org/). All intermediate and final computational results for the permeability are included in the file `permeability_*.log`. The different entries of the `permeability_*.log` file are briefly described below: -1. **iteration**: Number of the iteration. -2. **conv**: Convergence criterion. -3. **TPS**: Computed time steps per second. -4. **wmax_velz**: Maximum fluid voxel velocity (`z`-component) $[\mathrm{m} / \mathrm{s}]$. -5. **wmean_velz**: Mean fluid voxel velocity (`z`-component) $[\mathrm{m} / \mathrm{s}]$. -6. **k13**, **k23**, **k33**: Intrinsic permeabilities $[\mathrm{m}^2]$ in the given spatial directions. Pressure gradient and flux are computed based on the given domain of interest (see input file). If no domain of interest is given, $k_{13}$ = $k_{23}$ = $k_{33}$ = $0.0$. -7. **wk13**, **wk23**, **wk33**: Intrinsic permeabilities $[\mathrm{m}^2]$ in the given spatial directions. Pressure gradient and flux are computed based on the whole domain. +### Field files (`.raw`) + +Written when the corresponding `write_output` flag is set to `1`: + +| File pattern | Data type | Content | +|---|---|---| +| `velx_*.raw`, `vely_*.raw`, `velz_*.raw` | `double` | Final velocity field ($\mathbf{e}_1$, $\mathbf{e}_2$, $\mathbf{e}_3$ components) | +| `press_*.raw` | `double` | Final pressure field | +| `voxel_neighborhood_*.raw` | `int` | Voxel neighbourhood case codes | +| `domain_decomp_*.raw` | `int` | MPI domain decomposition map | + +Use `fields2vtu.py` to convert any of these to `.vtu` for visualisation in [ParaView](https://www.paraview.org/). + +### Log file (`.log`) + +One row per write event; columns in order: + +| Column | Unit | Description | +|---|---|---| +| `iteration` | — | Iteration number | +| `conv` | — | Relative permeability change (convergence criterion) | +| `TPS` | 1/s | Time steps computed per wall-clock second | +| `wmax_velz` | m/s | Maximum fluid-voxel z-velocity (whole domain) | +| `wmean_velz` | m/s | Mean fluid-voxel z-velocity (whole domain) | +| `k13`, `k23`, `k33` | m² | Permeability components for the `dom_interest` sub-domain (`0.0` if not set) | +| `wk13`, `wk23`, `wk33` | m² | Permeability components for the whole domain | +--- ## Aspects of parallelization -There is no simple answer to the question of how many cores can be recommended for specific domain sizes. We recommend at least one subdomain of $50^3$ voxels (for large porosities) and maximum $100^3$ voxels (for small porosities) per core. This should be an appropriate consideration, but is basically more of a rule of thumb than a requirement. In principle, significantly larger subdomains can be processed per core, and if possible, for example, 8 different simulations can be computed on one desktop computer. If there are concerns about parallelization, one can run a short simulation of a problem ($1000$ time steps) and take a look at the TPS values in the log file. This might also useful/interesting when experimenting with fixed domain decompositions. In our paper, see below, you can also find core numbers we used for larger simulations ($>1500^3$ voxels). +A useful rule of thumb is to target **50³–100³ voxels per MPI rank** (smaller sub-domains for high-porosity samples, larger for low-porosity ones). Larger sub-domains are technically possible and may be preferable when memory bandwidth is the bottleneck. + +To assess parallel efficiency for a given problem, run a short simulation (a few hundred iterations) and inspect the TPS value in the log. Experimenting with manual `dom_decomposition` settings versus the automatic (`0 0 0`) decomposition can also be worthwhile for non-cubic domains. Core counts used for larger simulations (> 1500³ voxels) are reported in the paper below. + +--- ## Compute the permeability tensor -The $z-$ or $\mathbf{e}_3-$ direction is the direction of the main pressure gradient. Rotate the domain (e.g. with `numpy.transpose()`) and rotate all results respectively to get the entries in a column. -Example: -original domain: **k13** -> $k_{13}$, **k23** -> $k_{23}$, **k33** -> $k_{33}$ +The pressure gradient is always applied along $\mathbf{e}_3$ (z). To obtain off-diagonal tensor entries, rotate the domain and run the solver once per orientation. Each run fills one column of the permeability tensor: -`np.transpose(domain, (2, 0, 1))`: **k13** -> $k_{32}$, **k23** -> $k_{12}$, **k33** -> $k_{22}$ +| `numpy.transpose` call | Output column | +|---|---| +| *(no transpose — original)* | $k_{13},\ k_{23},\ k_{33}$ | +| `np.transpose(domain, (2, 0, 1))` | $k_{32},\ k_{12},\ k_{22}$ | +| `np.transpose(domain, (1, 2, 0))` | $k_{21},\ k_{31},\ k_{11}$ | -`np.transpose(domain, (1, 2, 0))`: **k13** -> $k_{21}$, **k23** -> $k_{31}$, **k33** -> $k_{11}$ +Remember to rotate the result vectors by the same permutation to map log-file entries back to physical tensor indices. + +--- ## License -The solver is licensed under the terms and conditions of the MIT License version 3 or - at your option - any later -version. The License can be [found online](https://opensource.org/license/mit/) or in the LICENSE.md file -provided in the topmost directory of source code tree. +POREMAPS is released under the [MIT License](https://opensource.org/license/mit/). See `LICENSE.md` in the repository root for the full text. + +--- ## How to cite -The solver is research software and developed at a research institute. Please cite **specific releases** according to [**DaRUS**](https://doi.org/10.18419/darus-3676) version. +Please cite **specific releases** via [DaRUS](https://doi.org/10.18419/darus-3676). -If you are using poremaps in scientific publications and in the academic context, please cite our publications: +If you use POREMAPS in a scientific publication, please also cite: ```bib @article{Krach2025a, - author={Krach, David and Ruf, Matthias and Steeb, Holger}, - title={{POREMAPS}: A finite difference based Porous Media Anisotropic Permeability Solver for Stokes flow}, - DOI={10.69631/ipj.v2i1nr39}, - journal={InterPore Journal}, - pages={IPJ260225–7}, - volume={2}, - number={1}, - month={Feb.}, - year={2025}, - place={De Bilt, The Netherlands} + author = {Krach, David and Ruf, Matthias and Steeb, Holger}, + title = {{POREMAPS}: A finite difference based Porous Media Anisotropic Permeability Solver for Stokes flow}, + DOI = {10.69631/ipj.v2i1nr39}, + journal = {InterPore Journal}, + pages = {IPJ260225--7}, + volume = {2}, + number = {1}, + month = {Feb.}, + year = {2025}, + place = {De Bilt, The Netherlands} } ``` ```bib @data{Krach2024a, - author = {Krach, David and Ruf, Matthias and Steeb, Holger}, + author = {Krach, David and Ruf, Matthias and Steeb, Holger}, publisher = {DaRUS}, - title = {{POREMAPS 1.0.0: Code, Benchmarks, Applications}}, - year = {2024}, - version = {V1}, - doi = {10.18419/darus-3676}, - url = {https://doi.org/10.18419/darus-3676} + title = {{POREMAPS 1.0.0: Code, Benchmarks, Applications}}, + year = {2024}, + version = {V1}, + doi = {10.18419/darus-3676}, + url = {https://doi.org/10.18419/darus-3676} } ``` -## Solver is used in following publications +--- + +## Solver is used in the following publications [![Identifier](https://img.shields.io/badge/Publication_ADWR_Krach_et.al._(2025)-blue)](https://doi.org/10.1016/j.advwatres.2024.104860) ```bib @article{Krach2025b, - title = {A novel geometry-informed drag term formulation for pseudo-3D Stokes simulations with varying apertures}, + title = {A novel geometry-informed drag term formulation for pseudo-3D Stokes simulations with varying apertures}, journal = {Advances in Water Resources}, - volume = {195}, - year = {2025}, - doi = {https://doi.org/10.1016/j.advwatres.2024.104860}, - author = {David Krach and Felix Weinhardt and Mingfeng Wang and Martin Schneider and Holger Class and Holger Steeb}, + volume = {195}, + year = {2025}, + doi = {https://doi.org/10.1016/j.advwatres.2024.104860}, + author = {David Krach and Felix Weinhardt and Mingfeng Wang and Martin Schneider and Holger Class and Holger Steeb}, keywords = {Porous media, Stokes flow, Biomineralization, Microfluidics, Image-based simulations, Computational efficiency versus accuracy} } - ``` +--- + ## Developer -- [David Krach](https://www.mib.uni-stuttgart.de/institute/team/Krach/) E-mail: [david.krach@mib.uni-stuttgart.de](mailto:david.krach@mib.uni-stuttgart.de) -- [Matthias Ruf](https://www.mib.uni-stuttgart.de/institute/team/Ruf-00001/) E-mail: [matthias.ruf@mib.uni-stuttgart.de](mailto:matthias.ruf@mib.uni-stuttgart.de) +- [David Krach](https://www.mib.uni-stuttgart.de/institute/team/Krach/) — [david.krach@mib.uni-stuttgart.de](mailto:david.krach@mib.uni-stuttgart.de) +- [Matthias Ruf](https://www.mib.uni-stuttgart.de/institute/team/Ruf-00001/) — [matthias.ruf@mib.uni-stuttgart.de](mailto:matthias.ruf@mib.uni-stuttgart.de) ## Contact -- [Software Support Institute of Applied Mechanics](mailto:software@mib.uni-stuttgart.de) -- [Data Support Institute of Applied Mechanics](mailto:data@mib.uni-stuttgart.de) + +- [Software Support — Institute of Applied Mechanics](mailto:software@mib.uni-stuttgart.de) +- [Data Support — Institute of Applied Mechanics](mailto:data@mib.uni-stuttgart.de) diff --git a/fields2vtu.py b/fields2vtu.py index 764b13f..cdd87d8 100644 --- a/fields2vtu.py +++ b/fields2vtu.py @@ -1,6 +1,6 @@ #!/usr/bin/python3 -# Copyright 2024 David Krach, Matthias Ruf +# Copyright 2024-2026 David Krach, Matthias Ruf # Permission is hereby granted, free of charge, to any person obtaining a copy of # this software and associated documentation files (the “Software”), to deal in diff --git a/requirements-test.txt b/requirements-test.txt new file mode 100644 index 0000000..1cc18fc --- /dev/null +++ b/requirements-test.txt @@ -0,0 +1,2 @@ +numpy +pytest diff --git a/simulations/README.md b/simulations/README.md index 58f463d..e571351 100644 --- a/simulations/README.md +++ b/simulations/README.md @@ -1,14 +1,65 @@ # Stokes Solver Benchmarks -### Poiseuille Flow in a Tube -Compare the velocity in main flow direction $v_3$ to the analytical solution: -$$v_3(r) = \frac{\Delta p}{4L\mu} (R^2 -r^2)$$ -See [Batchelor, 1967](https://doi.org/10.1017/CBO9780511800955). - -### Rectangular Channel Flow -Compare Flow $Q_i$ to analytical solution, cf. [White and Majdalani, 2006](https://books.google.de/books?id=fl6wPwAACAAJ). - -### Sphere-packings -Compare normalized permeability to Kozeny-Carman equation: -$$k_1^{KC} = \frac{D^2}{c_{KC}} \frac{\phi^3}{(1-\phi)^2}$$ -See [Carman, 1997](https://doi.org/10.1016/S0263-8762(97)80003-2) and [Kozeny, 1927](https://books.google.de/books?id=yERGGwAACAAJ). +This directory contains the three benchmark cases used to validate POREMAPS against known analytical solutions. + +## Folder structure + +``` +simulations/ +└── 99_benchmarks/ + ├── 01_poiseuille/ # Hagen-Poiseuille flow in a periodic tube + ├── 02_channel_flow/ # Pressure-driven flow in a rectangular channel + └── 03_spherepackings/ # Isotropic sphere-packing arrays (BCC, FCC) +``` + +Each subfolder contains the geometry (`.raw`), input file (`.inp`), and — for the Poiseuille case — pre-computed output fields and a log file. + +## How to run + +Build the binary first (see [../README.md](../README.md)), then run a benchmark from its subfolder, e.g.: + +```shell +cd simulations/99_benchmarks/01_poiseuille +mpirun -np 4 ../../../bin/POREMAPS input_ptube_54_54_50_vs_2e-05.inp +``` + +--- + +## Benchmark 1 — Hagen-Poiseuille Flow in a Tube (`01_poiseuille/`) + +**Geometry:** circular tube embedded in a 54 × 54 × 50 voxel periodic domain (voxel size $\Delta x = 2 \times 10^{-5}$ m). + +**Analytical reference:** for a circular tube of radius $R$ under an axial pressure gradient $\Delta p / L$ and dynamic viscosity $\mu$, the Hagen-Poiseuille velocity profile is: +$$v_3(r) = \frac{\Delta p}{4L\mu} (R^2 - r^2)$$ +See Batchelor, 1967. + +**What to check:** compare the computed $v_3$ field (from `velz_*.raw`) to the parabolic profile above, evaluated at the radial distance $r$ of each fluid voxel from the tube centre. + +**Included output:** pre-computed velocity, pressure, voxel-neighbourhood, and domain-decomposition fields are included alongside the log file. + +--- + +## Benchmark 2 — Rectangular Channel Flow (`02_channel_flow/`) + +**Geometry:** full rectangular cross-section channel, 24 × 14 × 50 voxels (voxel size $\Delta x = 5 \times 10^{-5}$ m; physical cross-section 1.2 mm × 0.7 mm). + +**Analytical reference:** the volumetric flow rate $Q$ through a rectangular duct has a closed-form Fourier-series solution that depends on the aspect ratio of the cross-section and the applied pressure gradient. See White and Majdalani, 2006. + +**What to check:** integrate the computed $v_3$ field (`velz_*.raw`) over the cross-section and compare the resulting flow rate $Q$ to the analytical value. + +--- + +## Benchmark 3 — Sphere Packings (`03_spherepackings/`) + +**Geometries:** two ordered sphere-packing arrays at porosity $\phi = 0.3$ (solid volume fraction 0.7), each 100 × 100 × 100 voxels (voxel size $\Delta x = 1 \times 10^{-5}$ m): + +| Input file | Lattice | Sphere diameter $D$ | +|---|---|---| +| `input_SPHERE3D_bcc_…` | Body-centred cubic (BCC) | 8.745 × 10⁻⁴ m | +| `input_SPHERE3D_fcc_…` | Face-centred cubic (FCC) | 6.940 × 10⁻⁴ m | + +**Analytical reference:** the Kozeny-Carman equation relates the intrinsic permeability $k$ to sphere diameter and porosity: +$$k^{KC} = \frac{D^2}{c_{KC}} \frac{\phi^3}{(1-\phi)^2}$$ +where $c_{KC}$ is the Kozeny-Carman constant. See Carman, 1997 and Kozeny, 1927. + +**What to check:** compare the whole-domain permeability $k_{33}$ (column `wk33` in the log file) to $k^{KC}$ computed from the geometry parameters above. diff --git a/src/boundary_conditions.cc b/src/boundary_conditions.cc index 14f11a7..897da2b 100644 --- a/src/boundary_conditions.cc +++ b/src/boundary_conditions.cc @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in @@ -51,8 +51,8 @@ void reapply_pressure_gradient(bool*** proc_geom, int* cart_coords, int* dims) { - int i, j, k; - int zmax_procs = dims[2] - 1; // max number of process in z dir in cart communicator + int j, k; + int zmax_procs = dims[2] - 1; // index of last rank in the z direction int lim = 2; if (cart_coords[2] == 0){ // reset inlet pressure for (j = lim; j < new_proc_size[1] - lim; j++){ diff --git a/src/boundary_conditions.h b/src/boundary_conditions.h index b7571f5..adfc4b5 100644 --- a/src/boundary_conditions.h +++ b/src/boundary_conditions.h @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in diff --git a/src/calc_flow.cc b/src/calc_flow.cc index d7ef449..53fed7f 100644 --- a/src/calc_flow.cc +++ b/src/calc_flow.cc @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in @@ -78,19 +78,15 @@ void calc_flow( bool*** proc_geom, char* log_file) { - int i, j, k; // - int iteration = 0; // Set current iteration to zero - double perm_conv = 1.0 + eps; // Convergence Criterion - double start_it0, end_it0, elapsed_seconds, secs_it_sum; // Timer related vars - double phi_global; // Porosity of the sample or domain of interest - double permeability = 0.0; // Voxel based K - double prev_permeability = 0.0; // k of previous timestep - double perm_diff = eps + 1.0; // k difference t+1 and t - double wmax_velz = 0.0, wmean_velz = 0.0; // storing velocity information - double k13 = 0.0, k23 = 0.0, k33 = 0.0; // storing kii information domain of interest - double wk13 = 0.0, wk23 = 0.0, wk33 = 0.0; // storing kii information whole domain - // main pressure gradient direction - bool firstline = true; // flag if first line in logfile should be written + int iteration = 0; // current iteration counter + double perm_conv = 1.0 + eps; // convergence criterion (initialised above threshold) + double start_it0, end_it0, elapsed_seconds; // per-iteration timing (rank 0 only) + double phi_global; // domain porosity + double permeability = 0.0; // mean z-velocity based permeability (voxel units) + double wmax_velz = 0.0, wmean_velz = 0.0; // whole-domain velocity diagnostics + double k13 = 0.0, k23 = 0.0, k33 = 0.0; // permeability tensor, domain of interest [m^2] + double wk13 = 0.0, wk23 = 0.0, wk33 = 0.0; // permeability tensor, whole domain [m^2] + bool firstline = true; // controls whether to write the header line in the log file // if porosity is given as -1.0 --> compute it by get_porosity if (porosity == -1.0){ @@ -119,8 +115,6 @@ void calc_flow( bool*** proc_geom, communicate_halos(vel_z, new_proc_size, neighbors, MPI_DOUBLE, comm_cart, bc, 1); // Set Dirichlet Boundarys for the pressure again reapply_pressure_gradient(proc_geom, size, new_proc_size, press, comm_cart, cart_coords, dims); - - MPI_Barrier(comm_cart); // evaluate convergence if (iteration%it_eval == 0 && iteration > 1000){ @@ -148,9 +142,7 @@ void calc_flow( bool*** proc_geom, phi_global); if ( rank == 0 ){ - // compute relevent properties - end_it0 = MPI_Wtime(); elapsed_seconds = end_it0 - start_it0; secs_it_sum += elapsed_seconds; - // Write logfile + end_it0 = MPI_Wtime(); elapsed_seconds = end_it0 - start_it0; write_logfile( iteration, k13, k23, k33, wk13, wk23, wk33, @@ -163,10 +155,9 @@ void calc_flow( bool*** proc_geom, } } } - MPI_Barrier(comm_cart); if (rank == 0 && iteration%100 == 0){ - end_it0 = MPI_Wtime(); elapsed_seconds = end_it0 - start_it0; secs_it_sum += elapsed_seconds; + end_it0 = MPI_Wtime(); elapsed_seconds = end_it0 - start_it0; printf("Iteration: %i\t Convergence: %e\t", iteration, perm_conv); printf("Time: %f\t TPS: %f\n", elapsed_seconds, (1./elapsed_seconds)); } diff --git a/src/calc_flow.h b/src/calc_flow.h index 1258c7c..f6e434f 100644 --- a/src/calc_flow.h +++ b/src/calc_flow.h @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in diff --git a/src/constants.h b/src/constants.h index fff58f9..f14ccea 100644 --- a/src/constants.h +++ b/src/constants.h @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in @@ -32,16 +32,14 @@ OTHER DEALINGS IN THE SOFTWARE. * Contents: * * Changelog: Version: 1.0.0 - -file_name_char : Number of bytes reserved in the memory for saving a char series -pass : Empty function to pass not yet specified functions -omega_sor : Relaxation parameter for the SOR solver -Re : Reynolds number -ndims : Number of dimensions of Cartesian grid -c2 : Dimensionless speed of sound -ndims : Number of spatial dimensions -dt : Timestep size play with this number to decrease - convergence time + * + * file_name_char : max bytes for a filename string + * pass : no-op placeholder macro + * omega_sor : relaxation factor for the SOR solver (1.0 = Gauss-Seidel) + * Re : dimensionless viscosity (Reynolds-like parameter) + * ndims : number of spatial dimensions of the Cartesian grid + * c2 : dimensionless speed-of-sound squared (artificial compressibility) + * dt : pseudo-time step size; reduce to improve stability ***********************************************************************/ #ifndef file_name_char diff --git a/src/evaluation.cc b/src/evaluation.cc index 1771b58..287f77b 100644 --- a/src/evaluation.cc +++ b/src/evaluation.cc @@ -1,28 +1,28 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf -Permission is hereby granted, free of charge, to any person obtaining a copy of -this software and associated documentation files (the “Software”), to deal in -the Software without restriction, including without limitation the rights to use, -copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the -Software, and to permit persons to whom the Software is furnished to do so, +Permission is hereby granted, free of charge, to any person obtaining a copy of +this software and associated documentation files (the "Software"), to deal in +the Software without restriction, including without limitation the rights to use, +copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the +Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: -The above copyright notice and this permission notice shall be included in all +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. -THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, -EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES -OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND -NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT -HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, -WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - * Authors: David Krach + * Authors: David Krach * Date: 2021/12 * Contact: david.krach@mib.uni-stuttgart.de * @@ -31,9 +31,15 @@ OTHER DEALINGS IN THE SOFTWARE. * * double compute_convergence: * compute convergence criterion - * Contents: - * + * Contents: + * * Changelog: Version: 1.0.0 + * + * krach 2025: removed redundant MPI_Barrier calls (MPI_Allreduce + * is itself a collective and does not need surrounding barriers). + * Batched multiple scalar MPI_Allreduce calls into single calls + * over arrays, reducing collective operations in + * compute_permeability from 13 to 5. ***********************************************************************/ @@ -45,23 +51,23 @@ OTHER DEALINGS IN THE SOFTWARE. #include "constants.h" // ********************************************************************** -void compute_permeability( bool*** proc_geom, - int* size, - int* new_proc_size, - double*** vel_x, - double*** vel_y, +void compute_permeability( bool*** proc_geom, + int* size, + int* new_proc_size, + double*** vel_x, + double*** vel_y, double*** vel_z, - double*** press, + double*** press, double* k13, double* k23, double* k33, double* wk13, double* wk23, double* wk33, - double* wmax_velz, - double* wmean_velz, + double* wmax_velz, + double* wmean_velz, MPI_Comm comm_cart, - int* starts, + int* starts, int* dom_interest, double voxelsize, double porosity) @@ -69,10 +75,10 @@ void compute_permeability( bool*** proc_geom, int i, j, k; int gi, gj, gk; // local coordinates int lim = 2; // to consider halos - double velx_sum = 0.0, vely_sum = 0.0, velz_sum = 0.0; + double velx_sum = 0.0, vely_sum = 0.0, velz_sum = 0.0; double wvelx_sum = 0.0, wvely_sum = 0.0, wvelz_sum = 0.0; double wvelz_max = 0.0; // store velocities - double p0sum = 0.0, pendsum = 0.0; // stores pressuers + double p0sum = 0.0, pendsum = 0.0; // stores pressures double pgrad_dofi = 1.0; // pressure gradient of subdomain int p0flcount = 0, pendflcount = 0; // voxel counter int flcount = 0, wflcount = 0; // voxel counter @@ -94,15 +100,15 @@ void compute_permeability( bool*** proc_geom, } } } - // if true --> compute properties of subdomain - if (dom_interest[0] == 0 && dom_interest[1] == 0 && - dom_interest[2] == 0 && dom_interest[3] == 0 && + // if true --> compute properties of subdomain + if (dom_interest[0] == 0 && dom_interest[1] == 0 && + dom_interest[2] == 0 && dom_interest[3] == 0 && dom_interest[4] == 0 && dom_interest[5] == 0){ flcount = 0; eval_dofi = false; } else { - // Check wether the subdomain of the rank is within the + // Check wether the subdomain of the rank is within the // the domain domain of interest for ( i = lim; i < new_proc_size[2] - lim; i++ ){ for ( j = lim; j < new_proc_size[1] - lim; j++ ){ @@ -113,7 +119,7 @@ void compute_permeability( bool*** proc_geom, gk = starts[0] + k - lim; if ( gk < dom_interest[0] || gk > dom_interest[1] || gj < dom_interest[2] || gj > dom_interest[3] || - gi < dom_interest[4] || gk > dom_interest[5] ){ + gi < dom_interest[4] || gi > dom_interest[5] ){ continue; } else { @@ -143,7 +149,7 @@ void compute_permeability( bool*** proc_geom, } // Make sure that there is no random value in storage and set all - // relevant vels to zero if there is no fluid voxel present on the + // relevant vels to zero if there is no fluid voxel present on the // subdomain of the rank if ( wflcount == 0 ){ wvelx_sum = 0.0; @@ -156,38 +162,31 @@ void compute_permeability( bool*** proc_geom, vely_sum = 0.0; velz_sum = 0.0; } + if ( p0flcount == 0 ){ p0sum = 0.0; } + if ( pendflcount == 0 ){ pendsum = 0.0; } - if ( p0flcount == 0 ){ - p0sum = 0.0; - } - if ( pendflcount == 0 ){ - pendsum = 0.0; - } - - // Communicate domain of interest properties if necessary - MPI_Barrier(comm_cart); - - // int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm) + // Communicate domain of interest properties if necessary. + // MPI_Allreduce is a collective that synchronises all ranks implicitly; + // no surrounding MPI_Barrier calls are needed. if ( eval_dofi == true ){ - MPI_Allreduce(MPI_IN_PLACE, &velx_sum, 1, MPI_DOUBLE, MPI_SUM, comm_cart); - MPI_Allreduce(MPI_IN_PLACE, &vely_sum, 1, MPI_DOUBLE, MPI_SUM, comm_cart); - MPI_Allreduce(MPI_IN_PLACE, &velz_sum, 1, MPI_DOUBLE, MPI_SUM, comm_cart); - MPI_Allreduce(MPI_IN_PLACE, &p0sum, 1, MPI_DOUBLE, MPI_SUM, comm_cart); - MPI_Allreduce(MPI_IN_PLACE, &pendsum, 1, MPI_DOUBLE, MPI_SUM, comm_cart); - MPI_Allreduce(MPI_IN_PLACE, &flcount, 1, MPI_INT, MPI_SUM, comm_cart); - MPI_Allreduce(MPI_IN_PLACE, &p0flcount, 1, MPI_INT, MPI_SUM, comm_cart); - MPI_Allreduce(MPI_IN_PLACE, &pendflcount, 1, MPI_INT, MPI_SUM, comm_cart); + // Batch all SUM-doubles into one call and SUM-ints into one call. + double dsum_buf[5] = {velx_sum, vely_sum, velz_sum, p0sum, pendsum}; + int isum_buf[3] = {flcount, p0flcount, pendflcount}; + MPI_Allreduce(MPI_IN_PLACE, dsum_buf, 5, MPI_DOUBLE, MPI_SUM, comm_cart); + MPI_Allreduce(MPI_IN_PLACE, isum_buf, 3, MPI_INT, MPI_SUM, comm_cart); + velx_sum = dsum_buf[0]; vely_sum = dsum_buf[1]; + velz_sum = dsum_buf[2]; p0sum = dsum_buf[3]; + pendsum = dsum_buf[4]; + flcount = isum_buf[0]; p0flcount = isum_buf[1]; + pendflcount = isum_buf[2]; } - MPI_Barrier(comm_cart); - - MPI_Allreduce(MPI_IN_PLACE, &wvelx_sum, 1, MPI_DOUBLE, MPI_SUM, comm_cart); - MPI_Allreduce(MPI_IN_PLACE, &wvely_sum, 1, MPI_DOUBLE, MPI_SUM, comm_cart); - MPI_Allreduce(MPI_IN_PLACE, &wvelz_sum, 1, MPI_DOUBLE, MPI_SUM, comm_cart); + // Whole-domain reductions: batch the three SUM-doubles, keep MAX separate. + double wdsum_buf[3] = {wvelx_sum, wvely_sum, wvelz_sum}; + MPI_Allreduce(MPI_IN_PLACE, wdsum_buf, 3, MPI_DOUBLE, MPI_SUM, comm_cart); MPI_Allreduce(MPI_IN_PLACE, &wflcount, 1, MPI_INT, MPI_SUM, comm_cart); MPI_Allreduce(MPI_IN_PLACE, &wvelz_max, 1, MPI_DOUBLE, MPI_MAX, comm_cart); - - MPI_Barrier(comm_cart); + wvelx_sum = wdsum_buf[0]; wvely_sum = wdsum_buf[1]; wvelz_sum = wdsum_buf[2]; // Compute permeabilities of domain of interest in real dimensions if ( eval_dofi == true ){ @@ -206,21 +205,18 @@ void compute_permeability( bool*** proc_geom, *wk13 = (wvelx_sum/wflcount) * (1.0/Re) * porosity * pow(voxelsize,2); *wk23 = (wvely_sum/wflcount) * (1.0/Re) * porosity * pow(voxelsize,2); *wk33 = (wvelz_sum/wflcount) * (1.0/Re) * porosity * pow(voxelsize,2); - - *wmean_velz = wvelz_sum/wflcount * voxelsize; - *wmax_velz = wvelz_max * voxelsize; - - MPI_Barrier(comm_cart); + *wmean_velz = wvelz_sum/wflcount * voxelsize; + *wmax_velz = wvelz_max * voxelsize; } -double compute_convergence( bool*** proc_geom, - int* size, - int* new_proc_size, - double*** vel_z, +double compute_convergence( bool*** proc_geom, + int* size, + int* new_proc_size, + double*** vel_z, double* permeability, MPI_Comm comm_cart) { @@ -244,30 +240,23 @@ double compute_convergence( bool*** proc_geom, } } - - // Make sure that there is no random value in storage and set all - // relevant vels to zero if there is no fluid voxel present on the + // relevant vels to zero if there is no fluid voxel present on the // subdomain of the rank - if (flcount == 0){ - velz_sum = 0.0; - } - - MPI_Barrier(comm_cart); - - // int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm) - MPI_Allreduce(MPI_IN_PLACE, &velz_sum, 1, MPI_DOUBLE, MPI_SUM, comm_cart); - MPI_Allreduce(MPI_IN_PLACE, &flcount, 1, MPI_INT, MPI_SUM, comm_cart); + if (flcount == 0){ velz_sum = 0.0; } - MPI_Barrier(comm_cart); + // Batch the two SUM reductions into one call each (double + int). + // MPI_Allreduce is itself a synchronising collective; no barrier needed. + double cbuf[2] = {velz_sum, (double)flcount}; + MPI_Allreduce(MPI_IN_PLACE, cbuf, 2, MPI_DOUBLE, MPI_SUM, comm_cart); + velz_sum = cbuf[0]; + flcount = (int)cbuf[1]; mean_velz = velz_sum/flcount; current_permeability = mean_velz/Re; // in vx^2; // Compute convergence criterion conv = fabs(current_permeability - *permeability)/current_permeability; *permeability = current_permeability; - MPI_Barrier(comm_cart); return conv; - -} \ No newline at end of file +} diff --git a/src/evaluation.h b/src/evaluation.h index 242ae2c..6208f5d 100644 --- a/src/evaluation.h +++ b/src/evaluation.h @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in @@ -26,7 +26,7 @@ OTHER DEALINGS IN THE SOFTWARE. * Date: 2021/12 * Contact: david.krach@mib.uni-stuttgart.de * - * Purpose: Header for evaluate.cc + * Purpose: Header for evaluation.cc * * Contents: * diff --git a/src/geometry.cc b/src/geometry.cc index 26ca90f..41e24ad 100644 --- a/src/geometry.cc +++ b/src/geometry.cc @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in @@ -176,7 +176,7 @@ void eval_geometry( bool*** proc_geom, int*** voxel_neighborhood) { // See Bentz et al 2007. - int i, j, k, l, direction; + int i, j, k, direction; int i_, j_, k_; int shift_factor = 0; @@ -354,10 +354,10 @@ double get_porosity( bool*** proc_geom, } } - MPI_Allreduce(MPI_IN_PLACE, &fluid_count, 1, MPI_DOUBLE, MPI_SUM, comm_cart); - MPI_Allreduce(MPI_IN_PLACE, &vox_count , 1, MPI_DOUBLE, MPI_SUM, comm_cart); - - MPI_Barrier(comm_cart); + // Batch both SUM reductions into one call. + double gbuf[2] = {fluid_count, vox_count}; + MPI_Allreduce(MPI_IN_PLACE, gbuf, 2, MPI_DOUBLE, MPI_SUM, comm_cart); + fluid_count = gbuf[0]; vox_count = gbuf[1]; g_porosity = fluid_count/vox_count; return g_porosity; diff --git a/src/geometry.h b/src/geometry.h index b0b49d0..3b95a74 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in diff --git a/src/initialize.cc b/src/initialize.cc index 8a96146..e2bf8a0 100644 --- a/src/initialize.cc +++ b/src/initialize.cc @@ -1,56 +1,57 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf -Permission is hereby granted, free of charge, to any person obtaining a copy of -this software and associated documentation files (the “Software”), to deal in -the Software without restriction, including without limitation the rights to use, -copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the -Software, and to permit persons to whom the Software is furnished to do so, +Permission is hereby granted, free of charge, to any person obtaining a copy of +this software and associated documentation files (the "Software"), to deal in +the Software without restriction, including without limitation the rights to use, +copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the +Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: -The above copyright notice and this permission notice shall be included in all +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. -THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, -EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES -OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND -NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT -HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, -WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - * Authors: David Krach + * Authors: David Krach * Date: 2021/12 * Contact: david.krach@mib.uni-stuttgart.de * * Purpose: void allocate_initial_field: - * Initialize fields and commuincate halos - * + * Allocate and initialise velocity and pressure fields, + * then communicate pressure halos. + * * void determine_comm_pattern: - * Get ranks communication pattern depending on - * domain decomposition and BCs + * Set per-rank bc[] array that controls which directions + * communicate (MPI) and which apply a wall condition. * * void set_halos_initial: - * Set vel values in halos - * Contents: - * + * Zero velocity halo layers for no-slip boundary ranks. + * Contents: + * * Changelog: Version: 1.0.0 - -Each rank stores a array with 3 ints with holds information on how to -communicate in each spatial direction [x, y, z] = {0, 1, ..., 6} - -options: 0: communicate in both directions - 1: communicate in +1, slip condition in -1 - 2: communicate in +1, no-slip condition in -1 - 3: communicate in -1, slip condition in +1 - 4: communicate in -1, no-slip condition in +1 - 5: slip condition in -1 and +1 - 6: no-slip condition in -1 and +1 - -************************************************************************/ + * + * Each rank stores bc[3], one entry per spatial direction [x, y, z]. + * Values encode which wall (if any) the rank borders and the wall type: + * + * 0: communicate in both directions (interior rank) + * 1: communicate in +1, slip condition in -1 + * 2: communicate in +1, no-slip condition in -1 + * 3: communicate in -1, slip condition in +1 + * 4: communicate in -1, no-slip condition in +1 + * 5: slip condition in both -1 and +1 + * 6: no-slip condition in both -1 and +1 + * + ************************************************************************/ // HEADER *************************************************************** @@ -65,22 +66,21 @@ options: 0: communicate in both directions // ********************************************************************** -void allocate_initial_field(bool*** proc_geom, - int* size, - int* new_proc_size, - double*** press, +void allocate_initial_field(bool*** proc_geom, + int* size, + int* new_proc_size, + double*** press, double*** vel_x, - double*** vel_y, - double*** vel_z, - MPI_Comm comm_cart, - int* neighbors, - int* starts, + double*** vel_y, + double*** vel_z, + MPI_Comm comm_cart, + int* neighbors, + int* starts, unsigned int* bc) { - int i, j , k; - double gradp; + int i, j, k; - // allocate velocities and set initial value + // Allocate velocity and pressure arrays; set all fields to zero/initial values. for (i = 0; i < new_proc_size[2]; i++){ vel_x[i] = (double**)malloc(new_proc_size[1]*sizeof(double*)); vel_y[i] = (double**)malloc(new_proc_size[1]*sizeof(double*)); @@ -92,28 +92,22 @@ void allocate_initial_field(bool*** proc_geom, vel_z[i][j] = (double*)malloc(new_proc_size[0]*sizeof(double)); press[i][j] = (double*)malloc(new_proc_size[0]*sizeof(double)); for (k = 0; k < new_proc_size[0]; k++){ - if (proc_geom[i][j][k] == 0){ // initiate velocity for fluid voxels - vel_z[i][j][k] = 0.0; - } - else { - vel_z[i][j][k] = 0.0; - } vel_x[i][j][k] = 0.0; vel_y[i][j][k] = 0.0; + vel_z[i][j][k] = 0.0; press[i][j][k] = -1.0; - } } } - // initialize pressure gradient + + // Initialise pressure with a linear gradient along z (flow direction). + // Fluid voxels get a value proportional to their z position; solid voxels get 0. int lim = 2; for (i = lim; i < new_proc_size[2] - lim; i++){ for (j = lim; j < new_proc_size[1] - lim; j++){ for (k = lim; k < new_proc_size[0] - lim; k++){ - // distinguish between solid and fluid voxels if (proc_geom[i][j][k] == 0){ - // press[i][j][k] = (starts[2]+i-lim)+2; - press[i][j][k] = (size[2]-starts[2]-i)+lim+2; + press[i][j][k] = (size[2] - starts[2] - i) + lim + 2; } else { press[i][j][k] = 0.0; @@ -121,154 +115,142 @@ void allocate_initial_field(bool*** proc_geom, } } } - // MPI_Datatype etype = MPI_DOUBLE; - // Communicate pressure values in halos - communicate_halos(press, new_proc_size, neighbors, MPI_DOUBLE, comm_cart, bc, 0); + // Communicate pressure values into halo layers. + communicate_halos(press, new_proc_size, neighbors, MPI_DOUBLE, comm_cart, bc, 0); } -void determine_comm_pattern(int* dims, - int rank, - int* cart_coords, - unsigned int* bc, - int bc_method, +void determine_comm_pattern(int* dims, + int rank, + int* cart_coords, + unsigned int* bc, + int bc_method, int* neighbors) { /* -Possible input options in bc_method - -options: 0: requires: domain periodic in all spatial directions - communication: all ranks in all directions - set: nothing - - 1: requires: domain periodic in z direction - communication: all ranks in z directions , non-boundary interfaces - set: slip condition on x,y-boundary - - 2: requires: domain periodic in z direction - communication: all ranks in z directions , non-boundary interfaces - set: no-slip condition on x,y-boundary - - 3: requires: nothing - communication: non-boundary interfaces - set: slip condition on x,y-boundary , continouus flow in z direction - - 4: requires: nothing - communication: non-boundary interfaces - set: no-slip condition on x,y-boundary , continouus flow in z direction - -These are the 7 options to determine dependend on case and position of the rank in cart_coord + * bc_method selects the physical setup: + * + * 0: fully periodic in all directions + * → all ranks communicate in all directions + * + * 1: periodic in z; slip walls in x and y + * → all ranks communicate in z; boundary ranks in x/y apply slip + * + * 2: periodic in z; no-slip walls in x and y + * → all ranks communicate in z; boundary ranks in x/y apply no-slip + * + * 3: non-periodic; slip walls in x and y; pressure-driven flow in z + * → boundary ranks in all directions apply slip + * + * 4: non-periodic; no-slip walls in x and y; pressure-driven flow in z + * → boundary ranks in x/y apply no-slip; boundary ranks in z apply slip + */ -options: 0: communicate in both directions - 1: communicate in +1, slip condition in -1 - 2: communicate in +1, no-slip condition in -1 - 3: communicate in -1, slip condition in +1 - 4: communicate in -1, no-slip condition in +1 - 5: slip condition in -1 and +1 - 6: no-slip condition in -1 and +1 -*/ + int i; - int i, j, k; + // Default: all ranks communicate in both directions (interior behaviour). + // Boundary ranks below override their relevant entry. + for (i = 0; i < ndims; i++){ bc[i] = 0; } if (bc_method == 0){ - for (i = 0; i < ndims; i++){bc[i] = 0;} + // Fully periodic — bc stays 0 for all ranks and all directions. } - + else if (bc_method == 1){ - bc[2] = 0; - if (dims[0] == 1){ bc[0] = 5;} + // Periodic z; slip x and y walls. + bc[2] = 0; // z is periodic: all ranks communicate + + if (dims[0] == 1){ bc[0] = 5; } // single rank in x: slip on both sides else if (dims[0] > 1){ - if (cart_coords[0] == 0){ bc[0] = 1;} - if (cart_coords[0] + 1 == dims[0]){ bc[0] = 3;} + if (cart_coords[0] == 0) { bc[0] = 1; } // first rank in x: slip at -x + if (cart_coords[0] + 1 == dims[0]) { bc[0] = 3; } // last rank in x: slip at +x } - if (dims[1] == 1){ bc[1] = 5;} + if (dims[1] == 1){ bc[1] = 5; } else if (dims[1] > 1){ - if (cart_coords[1] == 0){ bc[1] = 1;} - if (cart_coords[1] + 1 == dims[1]){ bc[1] = 3;} - + if (cart_coords[1] == 0) { bc[1] = 1; } + if (cart_coords[1] + 1 == dims[1]) { bc[1] = 3; } } } else if (bc_method == 2){ + // Periodic z; no-slip x and y walls. bc[2] = 0; - if (dims[0] == 1){ bc[0] = 6;} + + if (dims[0] == 1){ bc[0] = 6; } else if (dims[0] > 1){ - if (cart_coords[0] == 0){ bc[0] = 2;} - if (cart_coords[0] + 1 == dims[0]){ bc[0] = 4;} + if (cart_coords[0] == 0) { bc[0] = 2; } + if (cart_coords[0] + 1 == dims[0]) { bc[0] = 4; } } - if (dims[1] == 1){ bc[1] = 6;} + if (dims[1] == 1){ bc[1] = 6; } else if (dims[1] > 1){ - if (cart_coords[1] == 0){ bc[1] = 2;} - if (cart_coords[1] + 1 == dims[1]){ bc[1] = 4;} - + if (cart_coords[1] == 0) { bc[1] = 2; } + if (cart_coords[1] + 1 == dims[1]) { bc[1] = 4; } } } else if (bc_method == 3){ - if (dims[0] == 1){ bc[0] = 5;} + // Non-periodic; slip walls in all directions. + if (dims[0] == 1){ bc[0] = 5; } else if (dims[0] > 1){ - if (cart_coords[0] == 0){ bc[0] = 1;} - if (cart_coords[0] + 1 == dims[0]){ bc[0] = 3;} + if (cart_coords[0] == 0) { bc[0] = 1; } + if (cart_coords[0] + 1 == dims[0]) { bc[0] = 3; } } - if (dims[1] == 1){ bc[1] = 5;} + if (dims[1] == 1){ bc[1] = 5; } else if (dims[1] > 1){ - if (cart_coords[1] == 0){ bc[1] = 1;} - if (cart_coords[1] + 1 == dims[1]){ bc[1] = 3;} - + if (cart_coords[1] == 0) { bc[1] = 1; } + if (cart_coords[1] + 1 == dims[1]) { bc[1] = 3; } } - if (dims[2] == 1){ bc[2] = 5;} + if (dims[2] == 1){ bc[2] = 5; } else if (dims[2] > 1){ - if (cart_coords[2] == 0){ bc[2] = 1;} - if (cart_coords[2] + 1 == dims[2]){ bc[2] = 3;} - + if (cart_coords[2] == 0) { bc[2] = 1; } + if (cart_coords[2] + 1 == dims[2]) { bc[2] = 3; } } } else if (bc_method == 4){ - if (dims[0] == 1){ bc[0] = 6;} + // Non-periodic; no-slip x and y walls; slip z walls (pressure-driven). + if (dims[0] == 1){ bc[0] = 6; } else if (dims[0] > 1){ - if (cart_coords[0] == 0){ bc[0] = 2;} - if (cart_coords[0] + 1 == dims[0]){ bc[0] = 4;} + if (cart_coords[0] == 0) { bc[0] = 2; } + if (cart_coords[0] + 1 == dims[0]) { bc[0] = 4; } } - if (dims[1] == 1){ bc[1] = 6;} + if (dims[1] == 1){ bc[1] = 6; } else if (dims[1] > 1){ - if (cart_coords[1] == 0){ bc[1] = 2;} - if (cart_coords[1] + 1 == dims[1]){ bc[1] = 4;} - + if (cart_coords[1] == 0) { bc[1] = 2; } + if (cart_coords[1] + 1 == dims[1]) { bc[1] = 4; } } - if (dims[2] == 1){ bc[2] = 5;} // 5 here is right ! + // z direction uses slip (not no-slip) even for bc_method 4 — the pressure + // gradient is imposed via reapply_pressure_gradient, not a wall condition. + if (dims[2] == 1){ bc[2] = 5; } else if (dims[2] > 1){ - if (cart_coords[2] == 0){ bc[2] = 2;} - if (cart_coords[2] + 1 == dims[2]){ bc[2] = 4;} - + if (cart_coords[2] == 0) { bc[2] = 2; } + if (cart_coords[2] + 1 == dims[2]) { bc[2] = 4; } } } +} -} - - -void set_halos_initial( double*** var, - int* new_proc_size, +void set_halos_initial( double*** var, + int* new_proc_size, unsigned int* bc) { - /* - initially set halos to zero - */ - int i, j, k; + // Zero the halo layers that correspond to no-slip wall boundaries. + // Only called for velocity fields; pressure halos are handled separately. + int i, j; + // X halos: indices 0,1 (left wall) and nx-2,nx-1 (right wall). for (i = 0; i < new_proc_size[2]; i++){ for (j = 0; j < new_proc_size[1]; j++){ if (bc[0] == 2 || bc[0] == 6){ var[i][j][0] = 0.0; - var[i][j][0] = 0.0; + var[i][j][1] = 0.0; // fixed: was erroneously set to [0] twice } if (bc[0] == 4 || bc[0] == 6){ var[i][j][new_proc_size[0]-1] = 0.0; @@ -277,6 +259,7 @@ void set_halos_initial( double*** var, } } + // Y halos: indices 0,1 (bottom wall) and ny-2,ny-1 (top wall). for (i = 0; i < new_proc_size[2]; i++){ for (j = 0; j < new_proc_size[0]; j++){ if (bc[1] == 2 || bc[1] == 6){ @@ -290,18 +273,17 @@ void set_halos_initial( double*** var, } } + // Z halos: indices 0,1 (front wall) and nz-2,nz-1 (back wall). for (i = 0; i < new_proc_size[1]; i++){ for (j = 0; j < new_proc_size[0]; j++){ - if (bc[3] == 2 || bc[3] == 6){ + if (bc[2] == 2 || bc[2] == 6){ // fixed: was bc[3] (out of bounds) var[0][i][j] = 0.0; var[1][i][j] = 0.0; } - if (bc[3] == 4 || bc[3] == 6){ + if (bc[2] == 4 || bc[2] == 6){ // fixed: was bc[3] (out of bounds) var[new_proc_size[2]-1][i][j] = 0.0; var[new_proc_size[2]-2][i][j] = 0.0; } } } } - - diff --git a/src/initialize.h b/src/initialize.h index 1929fbe..3b5ccc3 100644 --- a/src/initialize.h +++ b/src/initialize.h @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in @@ -26,9 +26,7 @@ OTHER DEALINGS IN THE SOFTWARE. * Date: 2021/02 * Contact: david.krach@mib.uni-stuttgart.de * - * Purpose: - * - * Contents: Header for initialize.cc + * Purpose: Header for initialize.cc * * Changelog: Version: 1.0.0 ***********************************************************************/ diff --git a/src/main_stokes_solver.cc b/src/main_stokes_solver.cc index 2f6cae4..b024907 100644 --- a/src/main_stokes_solver.cc +++ b/src/main_stokes_solver.cc @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in @@ -144,7 +144,7 @@ int main (int argc, char** argv){ // Define input file // If no command line argument is given the solver assumes the input file is called "input_param.inp" if (argc == 1){input_fn = "input_param.inp";} - else if (argc = 2){input_fn = (argv[1]);} + else if (argc == 2){input_fn = (argv[1]);} if (rank_glob == 0){ if (argc == 1){printf("No input file specified. Using input_param.inp\n");} diff --git a/src/makefile b/src/makefile index 863a9c4..2d984c0 100644 --- a/src/makefile +++ b/src/makefile @@ -1,33 +1,79 @@ -OBJECTS= read_problem_params.o geometry.o output.o parallelization.o initialize.o calc_flow.o solver.o boundary_conditions.o evaluation.o -CLINKER=path/to/my/mpiCC +# ----------------------------------------------------------------------- +# POREMAPS — Finite Difference Stokes Solver +# ----------------------------------------------------------------------- +# Override the compiler and flags on the command line if needed, e.g.: +# make MPICXX=mpic++ CXXFLAGS="-O2 -std=c++17" +# ----------------------------------------------------------------------- -FD_stokes_solver: $(OBJECTS) main_stokes_solver.cc - -$(CLINKER) -o ../bin/POREMAPS $(OBJECTS) main_stokes_solver.cc -lm - rm -f $(OBJECTS) +MPICXX ?= mpiCC +CXXFLAGS = -O3 -std=c++14 -Wall +LDFLAGS = -lm +TARGET = ../bin/POREMAPS -read_problem_params.o: read_problem_params.h read_problem_params.cc - -$(CLINKER) -o read_problem_params.o -c read_problem_params.cc +SRCS = read_problem_params.cc \ + geometry.cc \ + output.cc \ + parallelization.cc \ + initialize.cc \ + calc_flow.cc \ + solver.cc \ + boundary_conditions.cc \ + evaluation.cc \ + main_stokes_solver.cc -geometry.o: constants.h geometry.h geometry.cc - -$(CLINKER) -o geometry.o -c geometry.cc -lm +OBJS = $(SRCS:.cc=.o) -output.o: output.h output.cc - -$(CLINKER) -o output.o -c output.cc -lm +# ----------------------------------------------------------------------- +.PHONY: all clean debug -parallelization.o: parallelization.h parallelization.cc - -$(CLINKER) -o parallelization.o -c parallelization.cc +all: $(TARGET) -initialize.o: initialize.h initialize.cc - -$(CLINKER) -o initialize.o -c initialize.cc -lm +debug: CXXFLAGS = -O0 -g -std=c++14 -Wall +debug: $(TARGET) -calc_flow.o: constants.h solver.h boundary_conditions.h output.h calc_flow.h calc_flow.cc - -$(CLINKER) -o calc_flow.o -c calc_flow.cc -lm +# Link +$(TARGET): $(OBJS) + $(MPICXX) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) -solver.o: constants.h solver.h solver.cc - -$(CLINKER) -o solver.o -c solver.cc +# Compile (pattern rule covers all translation units) +%.o: %.cc + $(MPICXX) $(CXXFLAGS) -c -o $@ $< -boundary_conditions.o: constants.h boundary_conditions.h boundary_conditions.cc - -$(CLINKER) -o boundary_conditions.o -c boundary_conditions.cc +# ----------------------------------------------------------------------- +# Header dependencies (keep in sync with #includes in each .cc file) +# ----------------------------------------------------------------------- +read_problem_params.o: read_problem_params.cc read_problem_params.h -evaluation.o: constants.h evaluation.h evaluation.cc - -$(CLINKER) -o evaluation.o -c evaluation.cc +geometry.o: geometry.cc geometry.h \ + constants.h parallelization.h + +output.o: output.cc output.h \ + constants.h + +parallelization.o: parallelization.cc parallelization.h \ + constants.h + +initialize.o: initialize.cc initialize.h \ + constants.h parallelization.h + +calc_flow.o: calc_flow.cc calc_flow.h \ + constants.h parallelization.h solver.h \ + boundary_conditions.h geometry.h \ + evaluation.h output.h + +solver.o: solver.cc solver.h \ + constants.h + +boundary_conditions.o: boundary_conditions.cc boundary_conditions.h \ + constants.h + +evaluation.o: evaluation.cc evaluation.h \ + constants.h + +main_stokes_solver.o: main_stokes_solver.cc \ + constants.h read_problem_params.h geometry.h \ + output.h calc_flow.h parallelization.h initialize.h + +# ----------------------------------------------------------------------- +clean: + rm -f $(OBJS) $(TARGET) diff --git a/src/output.cc b/src/output.cc index 1ef9d61..ab012cb 100644 --- a/src/output.cc +++ b/src/output.cc @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in @@ -91,18 +91,19 @@ void write_output_raw( double*** var, } MPI_File_set_view(fh, write_offset, etype, filetype, "native", MPI_INFO_NULL); - MPI_File_write_at_all( fh, - write_offset, - temp_var, - proc_size[2]*proc_size[1]*proc_size[0], - etype, + MPI_File_write_at_all( fh, + write_offset, + temp_var, + proc_size[2]*proc_size[1]*proc_size[0], + etype, &write_status); delete [] temp_var; - + MPI_File_close(&fh); + MPI_Type_free(&filetype); } -void write_output_raw( int*** var, +void write_output_raw( int*** var, char* file_name, int* size, int* proc_size, @@ -139,19 +140,19 @@ void write_output_raw( int*** var, } } MPI_File_set_view(fh, write_offset, etype, filetype, "native", MPI_INFO_NULL); - MPI_File_write_at_all( fh, - write_offset, - temp_var, - proc_size[2]*proc_size[1]*proc_size[0], - etype, + MPI_File_write_at_all( fh, + write_offset, + temp_var, + proc_size[2]*proc_size[1]*proc_size[0], + etype, &write_status); delete [] temp_var; - + MPI_File_close(&fh); + MPI_Type_free(&filetype); } - void write_domain_decomposition(int*** domain_decomposition, char* file_name, int* size, diff --git a/src/output.h b/src/output.h index 21acbf4..ecae214 100644 --- a/src/output.h +++ b/src/output.h @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in diff --git a/src/parallelization.cc b/src/parallelization.cc index b178040..4ed64c9 100644 --- a/src/parallelization.cc +++ b/src/parallelization.cc @@ -1,28 +1,28 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf -Permission is hereby granted, free of charge, to any person obtaining a copy of -this software and associated documentation files (the “Software”), to deal in -the Software without restriction, including without limitation the rights to use, -copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the -Software, and to permit persons to whom the Software is furnished to do so, +Permission is hereby granted, free of charge, to any person obtaining a copy of +this software and associated documentation files (the "Software"), to deal in +the Software without restriction, including without limitation the rights to use, +copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the +Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: -The above copyright notice and this permission notice shall be included in all +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. -THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, -EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES -OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND -NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT -HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, -WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - * Authors: David Krach + * Authors: David Krach * Date: 2021/12 * Contact: david.krach@mib.uni-stuttgart.de * @@ -31,20 +31,27 @@ OTHER DEALINGS IN THE SOFTWARE. * * void communicate_halos: * pressure and velocity - * all 3 directions, if no communication is + * all 3 directions, if no communication is * needed due to BCs or small number of ranks - * MPI_PROC_NULL means that there is no neighbour + * MPI_PROC_NULL means that there is no neighbour * (defined as -2 in MPI4) - * Contents: - * + * Contents: + * * Changelog: Version: 1.0.0 krach 2021/01/14 added different cases for bc - added pindicator since pressure-bc have - to be treated different than velocity bc + added pindicator since pressure-bc have + to be treated different than velocity bc in case of no slip conditions if pindicator == 0 --> var is pressure if pindicator == 1 --> var is velocity + + krach 2026: replaced element-by-element MPI_Sendrecv loops + with pack/send/unpack pattern. + Each direction now requires 2 MPI_Sendrecv calls + (6 total) instead of O(ny*nz) calls. + Also fixed MPI derived datatype leak in + communicate_geom_halos. ***********************************************************************/ @@ -57,143 +64,255 @@ OTHER DEALINGS IN THE SOFTWARE. // ********************************************************************** -void communicate_geom_halos(bool*** var, - int* new_proc_size, - int* neighbors, - MPI_Datatype etype, +void communicate_geom_halos(bool*** var, + int* new_proc_size, + int* neighbors, + MPI_Datatype etype, MPI_Comm comm_cart) { - int i, j, k; - MPI_Status srstatus[12]; // MPI_Status object for all communication dirs - MPI_Datatype zsend, ysend, xsend; // derived datatypes for communication - - int zcount = new_proc_size[0]; - int ycount = new_proc_size[0]; - int xcount = 1; - - MPI_Type_contiguous(zcount, etype, &zsend); - MPI_Type_commit(&zsend); - MPI_Type_contiguous(ycount, etype, &ysend); - MPI_Type_commit(&ysend); - MPI_Type_contiguous(xcount, etype, &xsend); - MPI_Type_commit(&xsend); - - // comunicate in x direction <--> size[0], neighbor[0]: -1 in xdir, neighbor[1]: +1 in xdir - // communication bool by bool is improveable - for (i = 0; i < new_proc_size[2]; i++){ - for (j = 0; j < new_proc_size[1]; j++){ - MPI_Sendrecv(&var[i][j][2], 1, xsend, neighbors[0], 0, - &var[i][j][new_proc_size[0]-2], 1, xsend, neighbors[1], 0, comm_cart, &srstatus[0]); - MPI_Sendrecv(&var[i][j][3], 1, xsend, neighbors[0], 1, - &var[i][j][new_proc_size[0]-1], 1, xsend, neighbors[1], 1, comm_cart, &srstatus[1]); - MPI_Sendrecv(&var[i][j][new_proc_size[0]-4], 1, xsend, neighbors[1], 2, - &var[i][j][0], 1, xsend, neighbors[0], 2, comm_cart, &srstatus[2]); - MPI_Sendrecv(&var[i][j][new_proc_size[0]-3], 1, xsend, neighbors[1], 3, - &var[i][j][1], 1, xsend, neighbors[0], 3, comm_cart, &srstatus[3]); - } - } + int i, j, k; + MPI_Status status; + + int nx = new_proc_size[0]; + int ny = new_proc_size[1]; + int nz = new_proc_size[2]; + + // --- X direction --- + // Send x=2,3 to left neighbor (neighbors[0]); recv from right (neighbors[1]) into x=nx-2,nx-1 + // Send x=nx-4,nx-3 to right neighbor (neighbors[1]); recv from left (neighbors[0]) into x=0,1 + { + int sz = nz * ny * 2; + bool *sl = new bool[sz], *sr = new bool[sz]; + bool *rl = new bool[sz], *rr = new bool[sz]; - // communicate in y direction <--> size[1], neighbor[2]: -1 in ydir, neighbor[3]: +1 in ydir - for (i = 0; i < new_proc_size[2]; i++){ - // synchronous blocking sendrecvs - MPI_Sendrecv(var[i][2], 1, ysend, neighbors[2], 4, - var[i][new_proc_size[1]-2], 1, ysend, neighbors[3], 4, comm_cart, &srstatus[4]); - MPI_Sendrecv(var[i][3], 1, ysend, neighbors[2], 5, - var[i][new_proc_size[1]-1], 1, ysend, neighbors[3], 5, comm_cart, &srstatus[5]); - MPI_Sendrecv(var[i][new_proc_size[1]-4], 1, ysend, neighbors[3], 6, - var[i][0], 1, ysend, neighbors[2], 6, comm_cart, &srstatus[6]); - MPI_Sendrecv(var[i][new_proc_size[1]-3], 1, ysend, neighbors[3], 7, - var[i][1], 1, ysend, neighbors[2], 7, comm_cart, &srstatus[7]); + for (i = 0; i < nz; i++) + for (j = 0; j < ny; j++) { + int idx = 2*(i*ny + j); + sl[idx] = var[i][j][2]; sl[idx+1] = var[i][j][3]; + sr[idx] = var[i][j][nx-4]; sr[idx+1] = var[i][j][nx-3]; + } + + MPI_Sendrecv(sl, sz, etype, neighbors[0], 0, + rr, sz, etype, neighbors[1], 0, comm_cart, &status); + MPI_Sendrecv(sr, sz, etype, neighbors[1], 1, + rl, sz, etype, neighbors[0], 1, comm_cart, &status); + + if (neighbors[1] != MPI_PROC_NULL) + for (i = 0; i < nz; i++) + for (j = 0; j < ny; j++) { + int idx = 2*(i*ny + j); + var[i][j][nx-2] = rr[idx]; var[i][j][nx-1] = rr[idx+1]; + } + + if (neighbors[0] != MPI_PROC_NULL) + for (i = 0; i < nz; i++) + for (j = 0; j < ny; j++) { + int idx = 2*(i*ny + j); + var[i][j][0] = rl[idx]; var[i][j][1] = rl[idx+1]; + } + + delete[] sl; delete[] sr; delete[] rl; delete[] rr; } - // communicate in z direction <--> size[2], neighbor[4]: -1 in zdir, neighbor[5]: +1 in zdir - for (i = 0; i < new_proc_size[1]; i++){ - // synchronous blocking sendrecvs - MPI_Sendrecv(var[2][i], 1, zsend, neighbors[4], 8, - var[new_proc_size[2]-2][i], 1, zsend, neighbors[5], 8, comm_cart, &srstatus[8]); - MPI_Sendrecv(var[3][i], 1, zsend, neighbors[4], 9, - var[new_proc_size[2]-1][i], 1, zsend, neighbors[5], 9, comm_cart, &srstatus[9]); - MPI_Sendrecv(var[new_proc_size[2]-4][i], 1, zsend, neighbors[5], 10, - var[0][i], 1, zsend, neighbors[4], 10, comm_cart, &srstatus[10]); - MPI_Sendrecv(var[new_proc_size[2]-3][i], 1, zsend, neighbors[5], 11, - var[1][i], 1, zsend, neighbors[4], 11, comm_cart, &srstatus[11]); + // --- Y direction --- + // Send y=2,3 to bottom neighbor (neighbors[2]); recv from top (neighbors[3]) into y=ny-2,ny-1 + // Send y=ny-4,ny-3 to top neighbor (neighbors[3]); recv from bottom (neighbors[2]) into y=0,1 + { + int sz = nz * nx * 2; + bool *sl = new bool[sz], *sr = new bool[sz]; + bool *rl = new bool[sz], *rr = new bool[sz]; + + for (i = 0; i < nz; i++) + for (k = 0; k < nx; k++) { + int idx = 2*(i*nx + k); + sl[idx] = var[i][2][k]; sl[idx+1] = var[i][3][k]; + sr[idx] = var[i][ny-4][k]; sr[idx+1] = var[i][ny-3][k]; + } + + MPI_Sendrecv(sl, sz, etype, neighbors[2], 2, + rr, sz, etype, neighbors[3], 2, comm_cart, &status); + MPI_Sendrecv(sr, sz, etype, neighbors[3], 3, + rl, sz, etype, neighbors[2], 3, comm_cart, &status); + + if (neighbors[3] != MPI_PROC_NULL) + for (i = 0; i < nz; i++) + for (k = 0; k < nx; k++) { + int idx = 2*(i*nx + k); + var[i][ny-2][k] = rr[idx]; var[i][ny-1][k] = rr[idx+1]; + } + + if (neighbors[2] != MPI_PROC_NULL) + for (i = 0; i < nz; i++) + for (k = 0; k < nx; k++) { + int idx = 2*(i*nx + k); + var[i][0][k] = rl[idx]; var[i][1][k] = rl[idx+1]; + } + + delete[] sl; delete[] sr; delete[] rl; delete[] rr; } + // --- Z direction --- + // Send z=2,3 to front neighbor (neighbors[4]); recv from back (neighbors[5]) into z=nz-2,nz-1 + // Send z=nz-4,nz-3 to back neighbor (neighbors[5]); recv from front (neighbors[4]) into z=0,1 + { + int sz = ny * nx * 2; + bool *sl = new bool[sz], *sr = new bool[sz]; + bool *rl = new bool[sz], *rr = new bool[sz]; + + for (j = 0; j < ny; j++) + for (k = 0; k < nx; k++) { + int idx = 2*(j*nx + k); + sl[idx] = var[2][j][k]; sl[idx+1] = var[3][j][k]; + sr[idx] = var[nz-4][j][k]; sr[idx+1] = var[nz-3][j][k]; + } + + MPI_Sendrecv(sl, sz, etype, neighbors[4], 4, + rr, sz, etype, neighbors[5], 4, comm_cart, &status); + MPI_Sendrecv(sr, sz, etype, neighbors[5], 5, + rl, sz, etype, neighbors[4], 5, comm_cart, &status); + + if (neighbors[5] != MPI_PROC_NULL) + for (j = 0; j < ny; j++) + for (k = 0; k < nx; k++) { + int idx = 2*(j*nx + k); + var[nz-2][j][k] = rr[idx]; var[nz-1][j][k] = rr[idx+1]; + } + + if (neighbors[4] != MPI_PROC_NULL) + for (j = 0; j < ny; j++) + for (k = 0; k < nx; k++) { + int idx = 2*(j*nx + k); + var[0][j][k] = rl[idx]; var[1][j][k] = rl[idx+1]; + } + + delete[] sl; delete[] sr; delete[] rl; delete[] rr; + } } -void communicate_halos( double*** var, - int* new_proc_size, - int* neighbors, - MPI_Datatype etype, - MPI_Comm comm_cart, - unsigned int* bc, +void communicate_halos( double*** var, + int* new_proc_size, + int* neighbors, + MPI_Datatype etype, + MPI_Comm comm_cart, + unsigned int* bc, int pindicator) { - int i, j, k; - MPI_Status srstatus[12]; // MPI_Status object for all communication dirs - MPI_Datatype zsend, ysend, xsend; // derived datatypes for communication - - int zcount = new_proc_size[0]; - int ycount = new_proc_size[0]; - int xcount = 1; - // int zmax_procs = dims[2] - 1; // max number of process in z dir in cart communicator - - // Definition of Datatypes to send/recv - MPI_Type_contiguous(zcount, etype, &zsend); - MPI_Type_commit(&zsend); - MPI_Type_contiguous(ycount, etype, &ysend); - MPI_Type_commit(&ysend); - MPI_Type_contiguous(xcount, etype, &xsend); - MPI_Type_commit(&xsend); - - // comunicate in x direction <--> size[0], neighbor[0]: -1 in xdir, neighbor[1]: +1 in xdir - // communication bool by bool is improveable - for (i = 0; i < new_proc_size[2]; i++){ - for (j = 0; j < new_proc_size[1]; j++){ - MPI_Sendrecv(&var[i][j][2], 1, xsend, neighbors[0], 0, - &var[i][j][new_proc_size[0]-2], 1, xsend, neighbors[1], 0, comm_cart, &srstatus[0]); - MPI_Sendrecv(&var[i][j][3], 1, xsend, neighbors[0], 1, - &var[i][j][new_proc_size[0]-1], 1, xsend, neighbors[1], 1, comm_cart, &srstatus[1]); - MPI_Sendrecv(&var[i][j][new_proc_size[0]-4], 1, xsend, neighbors[1], 2, - &var[i][j][0], 1, xsend, neighbors[0], 2, comm_cart, &srstatus[2]); - MPI_Sendrecv(&var[i][j][new_proc_size[0]-3], 1, xsend, neighbors[1], 3, - &var[i][j][1], 1, xsend, neighbors[0], 3, comm_cart, &srstatus[3]); - } - } + int i, j, k; + MPI_Status status; + + int nx = new_proc_size[0]; + int ny = new_proc_size[1]; + int nz = new_proc_size[2]; + + // --- X direction --- + // Send x=2,3 to left neighbor (neighbors[0]); recv from right (neighbors[1]) into x=nx-2,nx-1 + // Send x=nx-4,nx-3 to right neighbor (neighbors[1]); recv from left (neighbors[0]) into x=0,1 + { + int sz = nz * ny * 2; + double *sl = new double[sz], *sr = new double[sz]; + double *rl = new double[sz], *rr = new double[sz]; + + for (i = 0; i < nz; i++) + for (j = 0; j < ny; j++) { + int idx = 2*(i*ny + j); + sl[idx] = var[i][j][2]; sl[idx+1] = var[i][j][3]; + sr[idx] = var[i][j][nx-4]; sr[idx+1] = var[i][j][nx-3]; + } - // communicate in y direction <--> size[1], neighbor[2]: -1 in ydir, neighbor[3]: +1 in ydir - for (i = 0; i < new_proc_size[2]; i++){ - // synchronous blocking sendrecvs - MPI_Sendrecv(var[i][2], 1, ysend, neighbors[2], 4, - var[i][new_proc_size[1]-2], 1, ysend, neighbors[3], 4, comm_cart, &srstatus[4]); - MPI_Sendrecv(var[i][3], 1, ysend, neighbors[2], 5, - var[i][new_proc_size[1]-1], 1, ysend, neighbors[3], 5, comm_cart, &srstatus[5]); - MPI_Sendrecv(var[i][new_proc_size[1]-4], 1, ysend, neighbors[3], 6, - var[i][0], 1, ysend, neighbors[2], 6, comm_cart, &srstatus[6]); - MPI_Sendrecv(var[i][new_proc_size[1]-3], 1, ysend, neighbors[3], 7, - var[i][1], 1, ysend, neighbors[2], 7, comm_cart, &srstatus[7]); + MPI_Sendrecv(sl, sz, etype, neighbors[0], 0, + rr, sz, etype, neighbors[1], 0, comm_cart, &status); + MPI_Sendrecv(sr, sz, etype, neighbors[1], 1, + rl, sz, etype, neighbors[0], 1, comm_cart, &status); + + if (neighbors[1] != MPI_PROC_NULL) + for (i = 0; i < nz; i++) + for (j = 0; j < ny; j++) { + int idx = 2*(i*ny + j); + var[i][j][nx-2] = rr[idx]; var[i][j][nx-1] = rr[idx+1]; + } + + if (neighbors[0] != MPI_PROC_NULL) + for (i = 0; i < nz; i++) + for (j = 0; j < ny; j++) { + int idx = 2*(i*ny + j); + var[i][j][0] = rl[idx]; var[i][j][1] = rl[idx+1]; + } + + delete[] sl; delete[] sr; delete[] rl; delete[] rr; } - // communicate in z direction <--> size[2], neighbor[4]: -1 in zdir, neighbor[5]: +1 in zdir - for (i = 0; i < new_proc_size[1]; i++){ - - // synchronous blocking sendrecvs - MPI_Sendrecv(var[2][i], 1, zsend, neighbors[4], 8, - var[new_proc_size[2]-2][i], 1, zsend, neighbors[5], 8, comm_cart, &srstatus[8]); - MPI_Sendrecv(var[3][i], 1, zsend, neighbors[4], 9, - var[new_proc_size[2]-1][i], 1, zsend, neighbors[5], 9, comm_cart, &srstatus[9]); - MPI_Sendrecv(var[new_proc_size[2]-4][i], 1, zsend, neighbors[5], 10, - var[0][i], 1, zsend, neighbors[4], 10, comm_cart, &srstatus[10]); - MPI_Sendrecv(var[new_proc_size[2]-3][i], 1, zsend, neighbors[5], 11, - var[1][i], 1, zsend, neighbors[4], 11, comm_cart, &srstatus[11]); - + // --- Y direction --- + // Send y=2,3 to bottom neighbor (neighbors[2]); recv from top (neighbors[3]) into y=ny-2,ny-1 + // Send y=ny-4,ny-3 to top neighbor (neighbors[3]); recv from bottom (neighbors[2]) into y=0,1 + { + int sz = nz * nx * 2; + double *sl = new double[sz], *sr = new double[sz]; + double *rl = new double[sz], *rr = new double[sz]; + + for (i = 0; i < nz; i++) + for (k = 0; k < nx; k++) { + int idx = 2*(i*nx + k); + sl[idx] = var[i][2][k]; sl[idx+1] = var[i][3][k]; + sr[idx] = var[i][ny-4][k]; sr[idx+1] = var[i][ny-3][k]; + } + + MPI_Sendrecv(sl, sz, etype, neighbors[2], 2, + rr, sz, etype, neighbors[3], 2, comm_cart, &status); + MPI_Sendrecv(sr, sz, etype, neighbors[3], 3, + rl, sz, etype, neighbors[2], 3, comm_cart, &status); + + if (neighbors[3] != MPI_PROC_NULL) + for (i = 0; i < nz; i++) + for (k = 0; k < nx; k++) { + int idx = 2*(i*nx + k); + var[i][ny-2][k] = rr[idx]; var[i][ny-1][k] = rr[idx+1]; + } + + if (neighbors[2] != MPI_PROC_NULL) + for (i = 0; i < nz; i++) + for (k = 0; k < nx; k++) { + int idx = 2*(i*nx + k); + var[i][0][k] = rl[idx]; var[i][1][k] = rl[idx+1]; + } + + delete[] sl; delete[] sr; delete[] rl; delete[] rr; } - // Free the datatypes - MPI_Type_free(&zsend); - MPI_Type_free(&ysend); - MPI_Type_free(&xsend); + // --- Z direction --- + // Send z=2,3 to front neighbor (neighbors[4]); recv from back (neighbors[5]) into z=nz-2,nz-1 + // Send z=nz-4,nz-3 to back neighbor (neighbors[5]); recv from front (neighbors[4]) into z=0,1 + { + int sz = ny * nx * 2; + double *sl = new double[sz], *sr = new double[sz]; + double *rl = new double[sz], *rr = new double[sz]; -} + for (j = 0; j < ny; j++) + for (k = 0; k < nx; k++) { + int idx = 2*(j*nx + k); + sl[idx] = var[2][j][k]; sl[idx+1] = var[3][j][k]; + sr[idx] = var[nz-4][j][k]; sr[idx+1] = var[nz-3][j][k]; + } + + MPI_Sendrecv(sl, sz, etype, neighbors[4], 4, + rr, sz, etype, neighbors[5], 4, comm_cart, &status); + MPI_Sendrecv(sr, sz, etype, neighbors[5], 5, + rl, sz, etype, neighbors[4], 5, comm_cart, &status); + if (neighbors[5] != MPI_PROC_NULL) + for (j = 0; j < ny; j++) + for (k = 0; k < nx; k++) { + int idx = 2*(j*nx + k); + var[nz-2][j][k] = rr[idx]; var[nz-1][j][k] = rr[idx+1]; + } + + if (neighbors[4] != MPI_PROC_NULL) + for (j = 0; j < ny; j++) + for (k = 0; k < nx; k++) { + int idx = 2*(j*nx + k); + var[0][j][k] = rl[idx]; var[1][j][k] = rl[idx+1]; + } + + delete[] sl; delete[] sr; delete[] rl; delete[] rr; + } +} diff --git a/src/parallelization.h b/src/parallelization.h index 1162b37..fc9c689 100644 --- a/src/parallelization.h +++ b/src/parallelization.h @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in diff --git a/src/read_problem_params.cc b/src/read_problem_params.cc index 036273a..e5cfdba 100644 --- a/src/read_problem_params.cc +++ b/src/read_problem_params.cc @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in @@ -27,12 +27,12 @@ OTHER DEALINGS IN THE SOFTWARE. * Contact: david.krach@mib.uni-stuttgart.de * * Purpose: void read_problem_params: - * Read problem parameters. File is alwasy input_param.inp + * Read problem parameters. File is always input_param.inp * if not specified otherwise. * * void read_parallel_params: - * Read parallel parameters, same functiontype, reads only - * parameters needed for domain decompostion + * Read parallel parameters, same function type, reads only + * parameters needed for domain decomposition * Contents: * * Changelog: Version: 1.0.0 @@ -97,7 +97,7 @@ void read_problem_params( char* file_name, fprintf(stderr, "Setting EPS to a value smaller equal than zero."); } if (*it_write % *it_eval){ - fprintf(stderr, "Set writing and evaluating freq to proper values (it_write mod it_eval =! 0)."); + fprintf(stderr, "Set writing and evaluating freq to proper values (it_write mod it_eval != 0)."); } if (*porosity > 1.0 || (*porosity < -0.0 && *porosity != -1.0)){ fprintf(stderr, "Porosity has to be 0.0 < porosity <= 1.0 or -1.0 to be evaluated."); diff --git a/src/read_problem_params.h b/src/read_problem_params.h index e761ba7..a67d72d 100644 --- a/src/read_problem_params.h +++ b/src/read_problem_params.h @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in @@ -27,7 +27,7 @@ OTHER DEALINGS IN THE SOFTWARE. * Date: 2021/02 * Contact: david.krach@mib.uni-stuttgart.de * - * Purpose: Header file to read_problem_params.cc + * Purpose: Header for read_problem_params.cc * * Contents: * diff --git a/src/solver.cc b/src/solver.cc index 34ff014..3cf3710 100644 --- a/src/solver.cc +++ b/src/solver.cc @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in @@ -63,8 +63,6 @@ void get_2nd_derivation(int neighborhood_case, double &d2vdx2) { - double v1, v2, v3, v4; // velocity definition see sketch - if (neighborhood_case == 1) { // d2vdx2 = -8.0 * v2; d2vdx2 = -8.0 * cvel; @@ -177,7 +175,7 @@ void central_diff5( bool*** proc_geom, get_2nd_derivation(case_y, vel_x[i][j-2][k], vel_x[i][j-1][k], vel_x[i][j][k], vel_x[i][j+1][k], vel_x[i][j+2][k], d2vxdy2); get_2nd_derivation(case_y, vel_z[i][j-2][k], vel_z[i][j-1][k], vel_z[i][j][k], vel_z[i][j+1][k], vel_z[i][j+2][k], d2vzdy2); - // perpendicular direction y + // perpendicular direction z get_2nd_derivation(case_z, vel_x[i-2][j][k], vel_x[i-1][j][k], vel_x[i][j][k], vel_x[i+1][j][k], vel_x[i+2][j][k], d2vxdz2); get_2nd_derivation(case_z, vel_y[i-2][j][k], vel_y[i-1][j][k], vel_y[i][j][k], vel_y[i+1][j][k], vel_y[i+2][j][k], d2vydz2); @@ -214,16 +212,11 @@ void central_diff5( bool*** proc_geom, vel_y[i][j][k] = (1.0 - omega_sor) * vel_y[i][j][k] + omega_sor * (-dp_dy + vely_sum / Re) * dt; vel_z[i][j][k] = (1.0 - omega_sor) * vel_z[i][j][k] + omega_sor * (-dp_dz + velz_sum / Re) * dt; } - // if solver is 2, 3 the pressure is computed with updated vel_*[][][] + // Solvers 2 and 3 use already-updated velocities for the pressure step. if (solver == 2 || solver == 3){ v_sum = (vel_z[i+1][j][k] - vel_z[i][j][k] + vel_y[i][j+1][k] - vel_y[i][j][k] + vel_x[i][j][k+1] - vel_x[i][j][k]); } - if (solver == 1 || solver == 2 || solver == 3){ - press[i][j][k] = press[i][j][k] - c2 * v_sum * dt; - } - // else if (solver == 3){ - // press[i][j][k] = (1.0 - omega_sor) * press[i][j][k] - c2 * omega_sor * v_sum * dt; - // } + press[i][j][k] = press[i][j][k] - c2 * v_sum * dt; } } } diff --git a/src/solver.h b/src/solver.h index 2f6e6a2..3386c4a 100644 --- a/src/solver.h +++ b/src/solver.h @@ -1,7 +1,7 @@ /************************************************************************ Parallel Finite Difference Solver for Stokes Equations in Porous Media -Copyright 2024 David Krach, Matthias Ruf +Copyright 2024-2026 David Krach, Matthias Ruf Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..a51409f --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,79 @@ +# SPDX-License-Identifier: MIT +# +# Copyright 2024-2026 David Krach, Matthias Ruf +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +"""Shared pytest fixtures for the POREMAPS test suite.""" + +import numpy as np +import pytest + +from helpers import make_proc_geom + + +def pytest_addoption(parser): + parser.addoption( + "--run-integration", + action="store_true", + default=False, + help="Also run integration tests that invoke the POREMAPS binary.", + ) + + +@pytest.fixture +def run_integration(request): + return request.config.getoption("--run-integration") + + +# --------------------------------------------------------------------------- +# Small geometry fixtures +# --------------------------------------------------------------------------- + +@pytest.fixture +def full_fluid_geom(): + """8×8×8 all-fluid interior wrapped in solid halos; shape (12, 12, 12).""" + interior = np.zeros((8, 8, 8), dtype=bool) + return make_proc_geom(interior) + + +@pytest.fixture +def channel_z_geom(): + """10×6×6 channel open in z; solid walls at x=0/x_max and y=0/y_max. + + Interior shape (nz=10, ny=6, nx=6); fluid occupies y=1..4, x=1..4. + Full array shape (14, 10, 10) after halos. + """ + interior = np.zeros((10, 6, 6), dtype=bool) + interior[:, 0, :] = True # solid y-wall (bottom) + interior[:, -1, :] = True # solid y-wall (top) + interior[:, :, 0] = True # solid x-wall (left) + interior[:, :, -1] = True # solid x-wall (right) + return make_proc_geom(interior) + + +@pytest.fixture +def single_solid_voxel_geom(): + """6×6×6 mostly-fluid interior with one solid voxel at position (3,3,3). + + Full array shape (10, 10, 10) after halos. + """ + interior = np.zeros((6, 6, 6), dtype=bool) + interior[3, 3, 3] = True + return make_proc_geom(interior) diff --git a/tests/helpers.py b/tests/helpers.py new file mode 100644 index 0000000..996c3ed --- /dev/null +++ b/tests/helpers.py @@ -0,0 +1,427 @@ +# SPDX-License-Identifier: MIT +# +# Copyright 2024-2026 David Krach, Matthias Ruf +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +"""Pure-Python reference implementations of POREMAPS C++ kernels. + +Used by the pytest suite to verify numerical logic independently of the +MPI-parallel C++ code. + +Array conventions +----------------- +All 3-D arrays use numpy shape (nz+4, ny+4, nx+4), matching the C++ +convention where the first index is z, second is y, third is x. The +physical (inner) domain occupies ``[HALO:-HALO]`` in every dimension. +Geometry encoding: 0 = fluid, non-zero = solid. +""" + +import numpy as np + +# Physical constants — must match src/constants.h +Re = 0.01 +c2 = 1.5e6 +dt = 1.0e-4 +omega_sor = 1.2 +HALO = 2 # ghost-layer thickness on each face + + +# --------------------------------------------------------------------------- +# solver.cc +# --------------------------------------------------------------------------- + +def split_number(number): + """Decode a 3-digit base-10 neighbourhood code into (case_x, case_y, case_z). + + Mirrors solver.cc::split_number(). + """ + case_z = number % 10 + number //= 10 + case_y = number % 10 + number //= 10 + case_x = number % 10 + return case_x, case_y, case_z + + +def get_2nd_derivation(case, llvel, lvel, cvel, rvel, rrvel): + """Return the FD approximation to d²v/dx² for the given neighbourhood case. + + Mirrors solver.cc::get_2nd_derivation(). + + Parameters + ---------- + case : int + Neighbourhood case 1–6. + llvel, lvel, cvel, rvel, rrvel : float + Stencil values: second-left, left, centre, right, second-right. + """ + if case == 1: + return -8.0 * cvel + if case == 2: + return 8.0 / 3.0 * rvel - 16.0 / 3.0 * cvel + if case == 3: + return 8.0 / 3.0 * lvel - 16.0 / 3.0 * cvel + if case == 4: + return 2.0 * rvel - 5.0 * cvel - 1.0 / 5.0 * rrvel + if case == 5: + return 2.0 * lvel - 1.0 / 5.0 * llvel - 5.0 * cvel + if case == 6: + return lvel - 2.0 * cvel + rvel + raise ValueError(f"neighbourhood case must be 1–6, got {case}") + + +# --------------------------------------------------------------------------- +# geometry.cc +# --------------------------------------------------------------------------- + +def get_dom_limits(size, dims, cart_coords): + """Compute per-rank domain limits. + + Mirrors geometry.cc::get_dom_limits(). + + Parameters + ---------- + size : sequence[3] Global domain size [nx, ny, nz]. + dims : sequence[3] Cartesian decomposition [dx, dy, dz]. + cart_coords : sequence[3] This rank's coordinates in the Cartesian grid. + + Returns + ------- + proc_size : list[3] Owned voxels per direction (after extra-slice assignment). + new_proc_size : list[3] proc_size + 2*HALO. + starts : list[3] Global start indices (based on base proc_size). + ends : list[3] starts + base proc_size. + """ + base = [size[i] // dims[i] for i in range(3)] + addslice = [size[i] - dims[i] * base[i] for i in range(3)] + + starts = [0, 0, 0] + ends = [0, 0, 0] + for i in range(3): + if cart_coords[i] == 0: + cum = 0 + elif 0 < cart_coords[i] <= addslice[i]: + cum = cart_coords[i] + else: + cum = addslice[i] + starts[i] = cart_coords[i] * base[i] + cum + ends[i] = starts[i] + base[i] + + # Extra slice for early ranks (matches C++ order: ends computed before update) + proc_size = [base[i] + (1 if cart_coords[i] < addslice[i] else 0) + for i in range(3)] + new_proc_size = [p + 2 * HALO for p in proc_size] + return proc_size, new_proc_size, starts, ends + + +def eval_geometry(proc_geom): + """Encode voxel neighbourhood into a 3-digit base-10 code for each fluid voxel. + + Mirrors geometry.cc::eval_geometry() (Bentz et al. 2007). + + Parameters + ---------- + proc_geom : ndarray, bool, shape (nz, ny, nx) + 0 = fluid, non-zero = solid (including halos). + + Returns + ------- + nhood : ndarray, int32, same shape + Neighbourhood code: 0 for unconstrained fluid, positive 3-digit code for + near-wall fluid, -1 for solid voxels. Only interior voxels are processed. + """ + nhood = np.where(proc_geom == 0, 0, -1).astype(np.int32) + nz, ny, nx = proc_geom.shape + lim = HALO + + # (shift vector in [z,y,x], shift_factor) + directions = [ + ((0, 0, 1), 100), # x direction + ((0, 1, 0), 10), # y direction + ((1, 0, 0), 1), # z direction + ] + + for (di, dj, dk), sf in directions: + for i in range(lim, nz - lim): + for j in range(lim, ny - lim): + for k in range(lim, nx - lim): + if proc_geom[i, j, k] != 0: + continue # skip solid voxels + + lv = 0 if proc_geom[i - di, j - dj, k - dk] == 0 else 1 + rv = 0 if proc_geom[i + di, j + dj, k + dk] == 0 else 1 + llv = -1 + rrv = -1 + if lv == 0: + llv = 0 if proc_geom[i - 2*di, j - 2*dj, k - 2*dk] == 0 else 1 + if rv == 0: + rrv = 0 if proc_geom[i + 2*di, j + 2*dj, k + 2*dk] == 0 else 1 + + if lv == 1 and rv == 1: + case = 1 + elif lv == 1 and rv == 0 and rrv == 1: + case = 2 + elif llv == 1 and lv == 0 and rv == 1: + case = 3 + elif lv == 1 and rv == 0 and rrv == 0: + case = 4 + elif llv == 0 and lv == 0 and rv == 1: + case = 5 + elif lv == 0 and rv == 0: + case = 6 + else: + case = 0 # unclassified; should not occur in valid geometry + + nhood[i, j, k] += case * sf + + return nhood + + +def get_proc_porosity(proc_geom): + """Return the porosity of the interior domain (halos excluded). + + Mirrors geometry.cc::get_proc_porosity(). + """ + interior = proc_geom[HALO:-HALO, HALO:-HALO, HALO:-HALO] + fluid_count = int(np.count_nonzero(interior == 0)) + vox_count = interior.size + return fluid_count / vox_count + + +# --------------------------------------------------------------------------- +# initialize.cc +# --------------------------------------------------------------------------- + +def determine_comm_pattern(dims, cart_coords, bc_method): + """Return bc[3] for the given rank in the Cartesian topology. + + Mirrors initialize.cc::determine_comm_pattern(). + + bc values per direction: + 0: communicate both ways (interior rank) + 1: communicate +1, slip at -1 + 2: communicate +1, no-slip at -1 + 3: communicate -1, slip at +1 + 4: communicate -1, no-slip at +1 + 5: slip at both -1 and +1 + 6: no-slip at both -1 and +1 + + Parameters + ---------- + dims : sequence[3] Cartesian topology dimensions [dx, dy, dz]. + cart_coords : sequence[3] This rank's Cartesian coordinates. + bc_method : int 0–4 physical setup selector. + + Returns + ------- + bc : list[3] + """ + bc = [0, 0, 0] + + if bc_method == 0: + pass # fully periodic — bc stays 0 + + elif bc_method in (1, 2): + # Periodic z; slip (method 1) or no-slip (method 2) walls in x and y. + lo = 1 if bc_method == 1 else 2 # first-rank code + hi = 3 if bc_method == 1 else 4 # last-rank code + both = 5 if bc_method == 1 else 6 # single-rank code + bc[2] = 0 # z is always periodic + for ax in (0, 1): + if dims[ax] == 1: + bc[ax] = both + elif dims[ax] > 1: + if cart_coords[ax] == 0: + bc[ax] = lo + if cart_coords[ax] + 1 == dims[ax]: + bc[ax] = hi + + elif bc_method == 3: + # Non-periodic; slip walls in all three directions. + for ax in range(3): + if dims[ax] == 1: + bc[ax] = 5 + elif dims[ax] > 1: + if cart_coords[ax] == 0: + bc[ax] = 1 + if cart_coords[ax] + 1 == dims[ax]: + bc[ax] = 3 + + elif bc_method == 4: + # Non-periodic; no-slip x/y walls; no-slip z boundaries. + for ax in (0, 1): + if dims[ax] == 1: + bc[ax] = 6 + elif dims[ax] > 1: + if cart_coords[ax] == 0: + bc[ax] = 2 + if cart_coords[ax] + 1 == dims[ax]: + bc[ax] = 4 + if dims[2] == 1: + bc[2] = 5 + elif dims[2] > 1: + if cart_coords[2] == 0: + bc[2] = 2 + if cart_coords[2] + 1 == dims[2]: + bc[2] = 4 + + return bc + + +def set_halos_initial(var, bc): + """Zero velocity halo layers that correspond to no-slip wall boundaries (in-place). + + Mirrors initialize.cc::set_halos_initial(). + + Parameters + ---------- + var : ndarray, float, shape (nz, ny, nx) + bc : sequence[3] Boundary condition codes [x, y, z]. + """ + # X halos (last axis) + if bc[0] in (2, 6): + var[:, :, 0] = 0.0 + var[:, :, 1] = 0.0 + if bc[0] in (4, 6): + var[:, :, -1] = 0.0 + var[:, :, -2] = 0.0 + + # Y halos (middle axis) + if bc[1] in (2, 6): + var[:, 0, :] = 0.0 + var[:, 1, :] = 0.0 + if bc[1] in (4, 6): + var[:, -1, :] = 0.0 + var[:, -2, :] = 0.0 + + # Z halos (first axis) + if bc[2] in (2, 6): + var[0, :, :] = 0.0 + var[1, :, :] = 0.0 + if bc[2] in (4, 6): + var[-1, :, :] = 0.0 + var[-2, :, :] = 0.0 + + +# --------------------------------------------------------------------------- +# evaluation.cc +# --------------------------------------------------------------------------- + +def compute_permeability_whole_domain(vel_x, vel_y, vel_z, proc_geom, + voxelsize, porosity): + """Compute whole-domain permeability components (single rank, no MPI). + + Mirrors the whole-domain section of evaluation.cc::compute_permeability(). + + Returns + ------- + wk13, wk23, wk33 : float Permeability tensor components [m²]. + wmean_velz : float Mean z-velocity [m/s]. + wmax_velz : float Maximum z-velocity [m/s]. + """ + sl = (slice(HALO, -HALO),) * 3 + fluid = proc_geom[sl] == 0 + + if not fluid.any(): + return 0.0, 0.0, 0.0, 0.0, 0.0 + + wflcount = int(fluid.sum()) + wvelx_sum = float(vel_x[sl][fluid].sum()) + wvely_sum = float(vel_y[sl][fluid].sum()) + wvelz_sum = float(vel_z[sl][fluid].sum()) + wvelz_max = float(vel_z[sl][fluid].max()) + + wk13 = (wvelx_sum / wflcount) * (1.0 / Re) * porosity * voxelsize ** 2 + wk23 = (wvely_sum / wflcount) * (1.0 / Re) * porosity * voxelsize ** 2 + wk33 = (wvelz_sum / wflcount) * (1.0 / Re) * porosity * voxelsize ** 2 + + wmean_velz = wvelz_sum / wflcount * voxelsize + wmax_velz = wvelz_max * voxelsize + + return wk13, wk23, wk33, wmean_velz, wmax_velz + + +def compute_convergence(vel_z, proc_geom, permeability): + """Compute the convergence criterion (single rank, no MPI). + + Mirrors evaluation.cc::compute_convergence(). + + Parameters + ---------- + vel_z : ndarray, float, shape (nz, ny, nx) + proc_geom : ndarray, bool, shape (nz, ny, nx) + permeability : float Permeability from the previous iteration. + + Returns + ------- + conv : float Relative change |k_new - k_old| / k_new. + new_perm : float Updated permeability value. + """ + sl = (slice(HALO, -HALO),) * 3 + fluid = proc_geom[sl] == 0 + + if not fluid.any(): + return 0.0, permeability + + velz_sum = float(vel_z[sl][fluid].sum()) + flcount = int(fluid.sum()) + + mean_velz = velz_sum / flcount + current_permeability = mean_velz / Re + if current_permeability != 0.0: + conv = abs(current_permeability - permeability) / current_permeability + else: + conv = 0.0 + + return conv, current_permeability + + +# --------------------------------------------------------------------------- +# Geometry construction helpers +# --------------------------------------------------------------------------- + +def make_proc_geom(geom_interior, halo_value=1): + """Wrap an interior geometry array with HALO solid layers on every face. + + Parameters + ---------- + geom_interior : ndarray, shape (nz, ny, nx) 0=fluid, non-zero=solid. + halo_value : scalar Value to fill halo cells with (default 1 = solid). + + Returns + ------- + ndarray, same dtype, shape (nz+4, ny+4, nx+4). + """ + nz, ny, nx = geom_interior.shape + out = np.full((nz + 2*HALO, ny + 2*HALO, nx + 2*HALO), + halo_value, dtype=geom_interior.dtype) + out[HALO:-HALO, HALO:-HALO, HALO:-HALO] = geom_interior + return out + + +def write_geometry(geom_interior, path): + """Write a geometry to a binary .raw file (uint8, x-fastest / Fortran order). + + Parameters + ---------- + geom_interior : ndarray, shape (nz, ny, nx) 0=fluid, 1=solid. + path : str or Path-like. + """ + geom_interior.astype(np.uint8).ravel(order='C').tofile(path) diff --git a/tests/test_evaluation.py b/tests/test_evaluation.py new file mode 100644 index 0000000..b1ae516 --- /dev/null +++ b/tests/test_evaluation.py @@ -0,0 +1,203 @@ +# SPDX-License-Identifier: MIT +# +# Copyright 2024-2026 David Krach, Matthias Ruf +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +"""Tests for evaluation.cc functions: compute_permeability, compute_convergence.""" + +import numpy as np +import pytest + +from helpers import ( + compute_permeability_whole_domain, + compute_convergence, + make_proc_geom, + HALO, + Re, +) + + +def _fluid_geom(nz=6, ny=6, nx=6): + return make_proc_geom(np.zeros((nz, ny, nx), dtype=bool)) + + +def _vel(value, shape): + return np.full(shape, value, dtype=float) + + +class TestComputePermeabilityWholeDomain: + def test_uniform_velz_gives_correct_k33(self): + geom = _fluid_geom() + s = geom.shape + dx = 1e-6 + phi = 1.0 + wk13, wk23, wk33, _, _ = compute_permeability_whole_domain( + np.zeros(s), np.zeros(s), _vel(1.0, s), geom, dx, phi + ) + assert wk33 == pytest.approx((1.0 / Re) * phi * dx**2) + assert wk13 == pytest.approx(0.0) + assert wk23 == pytest.approx(0.0) + + def test_uniform_velx_gives_correct_k13(self): + geom = _fluid_geom() + s = geom.shape + dx = 1e-6 + wk13, wk23, wk33, _, _ = compute_permeability_whole_domain( + _vel(1.0, s), np.zeros(s), np.zeros(s), geom, dx, 1.0 + ) + assert wk13 == pytest.approx((1.0 / Re) * dx**2) + assert wk23 == pytest.approx(0.0) + assert wk33 == pytest.approx(0.0) + + def test_zero_velocity_gives_zero_permeability(self): + geom = _fluid_geom() + s = geom.shape + z = np.zeros(s) + wk13, wk23, wk33, wmean, wmax = compute_permeability_whole_domain( + z, z, z, geom, 1e-6, 1.0 + ) + assert wk13 == pytest.approx(0.0) + assert wk23 == pytest.approx(0.0) + assert wk33 == pytest.approx(0.0) + assert wmean == pytest.approx(0.0) + assert wmax == pytest.approx(0.0) + + def test_porosity_scales_permeability_linearly(self): + geom = _fluid_geom() + s = geom.shape + vz = _vel(1.0, s) + dx = 1e-6 + _, _, k33_p05, _, _ = compute_permeability_whole_domain( + np.zeros(s), np.zeros(s), vz, geom, dx, 0.5 + ) + _, _, k33_p10, _, _ = compute_permeability_whole_domain( + np.zeros(s), np.zeros(s), vz, geom, dx, 1.0 + ) + assert k33_p05 == pytest.approx(k33_p10 * 0.5) + + def test_voxelsize_scales_permeability_as_squared(self): + geom = _fluid_geom() + s = geom.shape + vz = _vel(1.0, s) + _, _, k33_1, _, _ = compute_permeability_whole_domain( + np.zeros(s), np.zeros(s), vz, geom, 1.0, 1.0 + ) + _, _, k33_2, _, _ = compute_permeability_whole_domain( + np.zeros(s), np.zeros(s), vz, geom, 2.0, 1.0 + ) + assert k33_2 == pytest.approx(k33_1 * 4.0) + + def test_mean_and_max_velz_scaled_by_voxelsize(self): + geom = _fluid_geom(4, 4, 4) + s = geom.shape + vz = np.zeros(s) + vz[HALO:-HALO, HALO:-HALO, HALO:-HALO] = 3.0 + dx = 2.0 + _, _, _, wmean, wmax = compute_permeability_whole_domain( + np.zeros(s), np.zeros(s), vz, geom, dx, 1.0 + ) + assert wmean == pytest.approx(3.0 * dx) + assert wmax == pytest.approx(3.0 * dx) + + def test_max_velz_is_maximum_not_mean(self): + # Half of interior has vel=1, other half has vel=2 → max=2, mean=1.5 + geom = _fluid_geom(4, 4, 4) + s = geom.shape + vz = np.zeros(s) + nz, ny, nx = 4, 4, 4 + vz[HALO:HALO+nz//2, HALO:-HALO, HALO:-HALO] = 1.0 + vz[HALO+nz//2:-HALO, HALO:-HALO, HALO:-HALO] = 2.0 + _, _, _, wmean, wmax = compute_permeability_whole_domain( + np.zeros(s), np.zeros(s), vz, geom, 1.0, 1.0 + ) + assert wmax == pytest.approx(2.0) + assert wmean == pytest.approx(1.5) + + def test_all_solid_returns_zeros(self): + geom = make_proc_geom(np.ones((6, 6, 6), dtype=bool)) + s = geom.shape + ones = np.ones(s) + wk13, wk23, wk33, wmean, wmax = compute_permeability_whole_domain( + ones, ones, ones, geom, 1e-6, 0.0 + ) + assert wk13 == 0.0 and wk23 == 0.0 and wk33 == 0.0 + assert wmean == 0.0 and wmax == 0.0 + + +class TestComputeConvergence: + def test_zero_change_gives_zero_convergence(self): + geom = _fluid_geom() + s = geom.shape + vel_z = _vel(1.0, s) + # Set old permeability to exactly what the current field will give. + old_perm = 1.0 / Re + conv, new_perm = compute_convergence(vel_z, geom, old_perm) + assert conv == pytest.approx(0.0) + assert new_perm == pytest.approx(old_perm) + + def test_permeability_updated_on_each_call(self): + geom = _fluid_geom() + s = geom.shape + _, p1 = compute_convergence(_vel(1.0, s), geom, 0.0) + _, p2 = compute_convergence(_vel(2.0, s), geom, p1) + assert p1 == pytest.approx(1.0 / Re) + assert p2 == pytest.approx(2.0 / Re) + + def test_convergence_is_relative_change(self): + geom = _fluid_geom() + s = geom.shape + old = 1.0 / Re + new = 2.0 / Re + conv, _ = compute_convergence(_vel(2.0, s), geom, old) + assert conv == pytest.approx(abs(new - old) / new) + + def test_from_zero_gives_convergence_one(self): + # |new - 0| / new = 1 + geom = _fluid_geom() + s = geom.shape + conv, _ = compute_convergence(_vel(2.0, s), geom, 0.0) + assert conv == pytest.approx(1.0) + + def test_velocity_doubled_doubles_permeability(self): + geom = _fluid_geom() + s = geom.shape + _, p1 = compute_convergence(_vel(1.0, s), geom, 0.0) + _, p2 = compute_convergence(_vel(2.0, s), geom, 0.0) + assert p2 == pytest.approx(2.0 * p1) + + def test_all_solid_returns_unchanged_permeability(self): + geom = make_proc_geom(np.ones((6, 6, 6), dtype=bool)) + s = geom.shape + old = 0.42 + conv, new_perm = compute_convergence(_vel(1.0, s), geom, old) + assert conv == pytest.approx(0.0) + assert new_perm == pytest.approx(old) + + def test_convergence_approaches_zero_over_iterations(self): + # Simulate successive iterations with slowly changing velocity. + geom = _fluid_geom() + s = geom.shape + perm = 0.0 + vel = 1.0 + for _ in range(10): + conv, perm = compute_convergence(_vel(vel, s), geom, perm) + vel *= 1.001 # tiny change each "iteration" + # After many iterations with tiny changes convergence should be small. + assert conv < 0.01 diff --git a/tests/test_geometry.py b/tests/test_geometry.py new file mode 100644 index 0000000..7ab4d5c --- /dev/null +++ b/tests/test_geometry.py @@ -0,0 +1,191 @@ +# SPDX-License-Identifier: MIT +# +# Copyright 2024-2026 David Krach, Matthias Ruf +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +"""Tests for geometry.cc functions: eval_geometry, get_dom_limits, get_proc_porosity.""" + +import numpy as np +import pytest + +from helpers import ( + eval_geometry, + get_dom_limits, + get_proc_porosity, + split_number, + make_proc_geom, + HALO, +) + + +class TestGetDomLimits: + def test_single_rank_exact_division(self): + size, dims, coords = [10, 8, 6], [1, 1, 1], [0, 0, 0] + ps, nps, starts, ends = get_dom_limits(size, dims, coords) + assert ps == size + assert nps == [s + 4 for s in size] + assert starts == [0, 0, 0] + assert ends == size + + def test_two_ranks_even_split(self): + # 10 voxels into 2 ranks of 5 each; no remainder + size, dims = [10, 10, 10], [2, 1, 1] + ps0, nps0, st0, _ = get_dom_limits(size, dims, [0, 0, 0]) + ps1, nps1, st1, _ = get_dom_limits(size, dims, [1, 0, 0]) + assert ps0[0] == 5 and ps1[0] == 5 + assert st0[0] == 0 and st1[0] == 5 + assert nps0[0] == 9 and nps1[0] == 9 + + def test_two_ranks_odd_split(self): + # 11 voxels → rank 0 gets 6, rank 1 gets 5 + size, dims = [11, 8, 8], [2, 1, 1] + ps0, _, st0, _ = get_dom_limits(size, dims, [0, 0, 0]) + ps1, _, st1, _ = get_dom_limits(size, dims, [1, 0, 0]) + assert ps0[0] + ps1[0] == 11 + assert ps0[0] == 6 and ps1[0] == 5 + # starts[1] accounts for rank 0's extra slice: 1*base + 1 = 6 + assert st1[0] == 6 + + def test_starts_cover_full_domain_three_equal_ranks(self): + # 12 voxels into 3 ranks of 4 each + size, dims = [12, 12, 12], [3, 1, 1] + starts = [get_dom_limits(size, dims, [c, 0, 0])[2][0] for c in range(3)] + assert starts == [0, 4, 8] + + def test_new_proc_size_has_halos(self): + size, dims = [10, 10, 10], [1, 1, 1] + _, nps, _, _ = get_dom_limits(size, dims, [0, 0, 0]) + assert nps == [s + 4 for s in size] + + def test_three_uneven_ranks_sizes_sum_to_global(self): + # 10 voxels into 3: base=3, remainder=1 → rank 0 gets 4, others get 3 + size, dims = [10, 8, 8], [3, 1, 1] + sizes = [get_dom_limits(size, dims, [c, 0, 0])[0][0] for c in range(3)] + assert sum(sizes) == size[0] + assert sizes[0] == 4 and sizes[1] == 3 and sizes[2] == 3 + + def test_y_and_z_also_split(self): + # Decompose all three dimensions + size, dims = [6, 9, 12], [2, 3, 4] + total_x = sum(get_dom_limits(size, dims, [c, 0, 0])[0][0] for c in range(2)) + total_y = sum(get_dom_limits(size, dims, [0, c, 0])[0][1] for c in range(3)) + total_z = sum(get_dom_limits(size, dims, [0, 0, c])[0][2] for c in range(4)) + assert total_x == size[0] + assert total_y == size[1] + assert total_z == size[2] + + +class TestGetProcPorosity: + def test_all_fluid(self, full_fluid_geom): + assert get_proc_porosity(full_fluid_geom) == pytest.approx(1.0) + + def test_all_solid_interior(self): + interior = np.ones((6, 6, 6), dtype=bool) + geom = make_proc_geom(interior) + assert get_proc_porosity(geom) == pytest.approx(0.0) + + def test_half_solid_half_fluid(self): + interior = np.zeros((6, 6, 6), dtype=bool) + interior[:, :, :3] = True # solid on one x-half + geom = make_proc_geom(interior) + assert get_proc_porosity(geom) == pytest.approx(0.5) + + def test_halos_are_excluded(self): + # Halos are solid by default; a fully-fluid interior still gives φ=1. + interior = np.zeros((4, 4, 4), dtype=bool) + geom = make_proc_geom(interior, halo_value=1) + assert get_proc_porosity(geom) == pytest.approx(1.0) + + def test_channel_porosity(self, channel_z_geom): + # Interior 10×6×6; walls at y=0, y=5, x=0, x=5 (4 fluid rows in each) + # fluid fraction = (4*4) / (6*6) = 16/36 + assert get_proc_porosity(channel_z_geom) == pytest.approx(16 / 36) + + def test_single_solid_voxel(self, single_solid_voxel_geom): + total = 6 ** 3 + fluid = total - 1 + assert get_proc_porosity(single_solid_voxel_geom) == pytest.approx(fluid / total) + + +class TestEvalGeometry: + def test_solid_voxels_stay_minus_one(self): + interior = np.zeros((6, 6, 6), dtype=bool) + interior[3, 3, 3] = True # one solid voxel + geom = make_proc_geom(interior) + nhood = eval_geometry(geom) + solid_mask = geom != 0 + assert (nhood[solid_mask] == -1).all() + + def test_fluid_voxels_have_nonnegative_code(self, full_fluid_geom): + nhood = eval_geometry(full_fluid_geom) + interior = (slice(HALO, -HALO),) * 3 + fluid = full_fluid_geom[interior] == 0 + assert (nhood[interior][fluid] >= 0).all() + + def test_isolated_fluid_voxel_all_directions_case1(self): + # Single fluid voxel with solid on all sides → case 1 in every direction + interior = np.ones((7, 7, 7), dtype=bool) + interior[3, 3, 3] = False # one fluid voxel + geom = make_proc_geom(interior) + nhood = eval_geometry(geom) + code = nhood[3 + HALO, 3 + HALO, 3 + HALO] + cx, cy, cz = split_number(code) + assert cx == 1 and cy == 1 and cz == 1 + + def test_fully_fluid_interior_voxel_all_directions_case6(self, full_fluid_geom): + # A voxel well away from the halo boundary should have code 666 + nhood = eval_geometry(full_fluid_geom) + nz, ny, nx = full_fluid_geom.shape + mid = (nz // 2, ny // 2, nx // 2) + assert nhood[mid] == 6 * 100 + 6 * 10 + 6 + + def test_first_fluid_x_in_channel_is_near_wall(self, channel_z_geom): + # First fluid voxel in x (interior x=1, array x=HALO+1=3) has: + # left = solid wall (case_x has wall-left component) + # right = fluid, next-right = fluid → case 4 + nhood = eval_geometry(channel_z_geom) + nz, ny, _ = channel_z_geom.shape + i = nz // 2 + j = ny // 2 + k = HALO + 1 # first fluid x-index in full array + cx, _, _ = split_number(nhood[i, j, k]) + assert cx == 4 + + def test_all_solid_geometry_all_minus_one(self): + interior = np.ones((6, 6, 6), dtype=bool) + geom = make_proc_geom(interior) + nhood = eval_geometry(geom) + assert (nhood[HALO:-HALO, HALO:-HALO, HALO:-HALO] == -1).all() + + def test_two_fluid_voxels_adjacent_case2_and_case3(self): + # Two adjacent fluid voxels in x surrounded by solid: + # left voxel: wall on left (lvox=1), fluid on right (rvox=0), wall beyond (rrvox=1) → case_x=2 + # right voxel: fluid on left (lvox=0), wall on right (rvox=1), wall beyond (llvox=1) → case_x=3 + interior = np.ones((5, 5, 5), dtype=bool) + interior[2, 2, 1] = False # left fluid + interior[2, 2, 2] = False # right fluid + geom = make_proc_geom(interior) + nhood = eval_geometry(geom) + # left voxel: array index (2+HALO, 2+HALO, 1+HALO) = (4, 4, 3) + cx_left, _, _ = split_number(nhood[2 + HALO, 2 + HALO, 1 + HALO]) + # right voxel: array index (4, 4, 4) + cx_right, _, _ = split_number(nhood[2 + HALO, 2 + HALO, 2 + HALO]) + assert cx_left == 2 + assert cx_right == 3 diff --git a/tests/test_initialize.py b/tests/test_initialize.py new file mode 100644 index 0000000..00e62ed --- /dev/null +++ b/tests/test_initialize.py @@ -0,0 +1,254 @@ +# SPDX-License-Identifier: MIT +# +# Copyright 2024-2026 David Krach, Matthias Ruf +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +"""Tests for initialize.cc functions: determine_comm_pattern, set_halos_initial.""" + +import numpy as np +import pytest + +from helpers import determine_comm_pattern, set_halos_initial, HALO + + +class TestDetermineCommPattern: + """Test bc[] array produced by determine_comm_pattern for each bc_method.""" + + # ------------------------------------------------------------------ + # bc_method 0: fully periodic — all ranks, all directions → bc = 0 + # ------------------------------------------------------------------ + def test_method0_single_rank(self): + assert determine_comm_pattern([1, 1, 1], [0, 0, 0], 0) == [0, 0, 0] + + def test_method0_interior_rank(self): + assert determine_comm_pattern([4, 4, 4], [2, 2, 2], 0) == [0, 0, 0] + + def test_method0_boundary_rank_still_zero(self): + assert determine_comm_pattern([4, 4, 4], [0, 0, 0], 0) == [0, 0, 0] + + # ------------------------------------------------------------------ + # bc_method 1: periodic z; slip walls in x and y + # ------------------------------------------------------------------ + def test_method1_single_rank_in_all_dims(self): + bc = determine_comm_pattern([1, 1, 1], [0, 0, 0], 1) + assert bc[0] == 5 and bc[1] == 5 and bc[2] == 0 + + def test_method1_z_is_always_zero(self): + for coords in ([0, 0, 0], [1, 1, 1], [2, 2, 2]): + bc = determine_comm_pattern([3, 3, 3], coords, 1) + assert bc[2] == 0 + + def test_method1_first_rank_in_x(self): + bc = determine_comm_pattern([3, 1, 1], [0, 0, 0], 1) + assert bc[0] == 1 + + def test_method1_last_rank_in_x(self): + bc = determine_comm_pattern([3, 1, 1], [2, 0, 0], 1) + assert bc[0] == 3 + + def test_method1_interior_rank_x_zero(self): + bc = determine_comm_pattern([3, 1, 1], [1, 0, 0], 1) + assert bc[0] == 0 + + def test_method1_first_rank_in_y(self): + bc = determine_comm_pattern([1, 4, 1], [0, 0, 0], 1) + assert bc[1] == 1 + + def test_method1_last_rank_in_y(self): + bc = determine_comm_pattern([1, 4, 1], [0, 3, 0], 1) + assert bc[1] == 3 + + # ------------------------------------------------------------------ + # bc_method 2: periodic z; no-slip walls in x and y + # ------------------------------------------------------------------ + def test_method2_single_rank(self): + bc = determine_comm_pattern([1, 1, 1], [0, 0, 0], 2) + assert bc[0] == 6 and bc[1] == 6 and bc[2] == 0 + + def test_method2_first_rank_in_x(self): + bc = determine_comm_pattern([2, 1, 1], [0, 0, 0], 2) + assert bc[0] == 2 + + def test_method2_last_rank_in_x(self): + bc = determine_comm_pattern([2, 1, 1], [1, 0, 0], 2) + assert bc[0] == 4 + + def test_method2_interior_rank_zero(self): + bc = determine_comm_pattern([3, 1, 1], [1, 0, 0], 2) + assert bc[0] == 0 + + def test_method2_z_is_always_zero(self): + bc = determine_comm_pattern([1, 1, 3], [0, 0, 1], 2) + assert bc[2] == 0 + + # ------------------------------------------------------------------ + # bc_method 3: non-periodic; slip walls in all directions + # ------------------------------------------------------------------ + def test_method3_single_rank_all(self): + bc = determine_comm_pattern([1, 1, 1], [0, 0, 0], 3) + assert bc == [5, 5, 5] + + def test_method3_first_rank_z(self): + bc = determine_comm_pattern([1, 1, 3], [0, 0, 0], 3) + assert bc[2] == 1 + + def test_method3_last_rank_z(self): + bc = determine_comm_pattern([1, 1, 3], [0, 0, 2], 3) + assert bc[2] == 3 + + def test_method3_interior_all_zero(self): + bc = determine_comm_pattern([3, 3, 3], [1, 1, 1], 3) + assert bc == [0, 0, 0] + + # ------------------------------------------------------------------ + # bc_method 4: non-periodic; no-slip x/y; no-slip z boundaries + # ------------------------------------------------------------------ + def test_method4_single_rank_x_y_noslip(self): + bc = determine_comm_pattern([1, 1, 1], [0, 0, 0], 4) + assert bc[0] == 6 and bc[1] == 6 + + def test_method4_single_rank_z_code(self): + # Single rank in z gets code 5 (as per C++ code) + bc = determine_comm_pattern([1, 1, 1], [0, 0, 0], 4) + assert bc[2] == 5 + + def test_method4_first_rank_z(self): + bc = determine_comm_pattern([1, 1, 3], [0, 0, 0], 4) + assert bc[2] == 2 + + def test_method4_last_rank_z(self): + bc = determine_comm_pattern([1, 1, 3], [0, 0, 2], 4) + assert bc[2] == 4 + + def test_method4_interior_all_zero(self): + bc = determine_comm_pattern([3, 3, 3], [1, 1, 1], 4) + assert bc == [0, 0, 0] + + # ------------------------------------------------------------------ + # All interior ranks are unaffected by bc_method (remain 0) + # ------------------------------------------------------------------ + def test_interior_rank_unaffected_by_bc_method(self): + for method in range(5): + bc = determine_comm_pattern([5, 5, 5], [2, 2, 2], method) + assert bc == [0, 0, 0], f"bc_method={method} changed interior rank bc" + + # ------------------------------------------------------------------ + # Default initialisation: bc starts at 0 for all ranks + # ------------------------------------------------------------------ + def test_default_is_zero_before_overrides(self): + # bc_method 0 leaves everything at 0 + bc = determine_comm_pattern([4, 4, 4], [1, 2, 2], 0) + assert bc == [0, 0, 0] + + +class TestSetHalosInitial: + def _ones(self, shape=(12, 12, 12)): + return np.ones(shape, dtype=float) + + def test_zero_bc_leaves_array_unchanged(self): + var = self._ones() + ref = var.copy() + set_halos_initial(var, [0, 0, 0]) + np.testing.assert_array_equal(var, ref) + + def test_slip_bcs_leave_array_unchanged(self): + # bc values 1, 3, 5 are slip — they must NOT zero any halo layers + for slip_bc in ( + [1, 0, 0], [3, 0, 0], [5, 0, 0], + [0, 1, 0], [0, 3, 0], [0, 5, 0], + [0, 0, 1], [0, 0, 3], [0, 0, 5], + [5, 5, 5], + ): + var = self._ones() + ref = var.copy() + set_halos_initial(var, slip_bc) + assert (var == ref).all(), f"bc={slip_bc} unexpectedly modified array" + + # X halos (axis 2, last dimension) + def test_no_slip_x_left_zeros_columns_0_and_1(self): + var = self._ones() + set_halos_initial(var, [2, 0, 0]) + assert (var[:, :, 0] == 0.0).all() + assert (var[:, :, 1] == 0.0).all() + assert (var[:, :, -1] == 1.0).all() # right side untouched + + def test_no_slip_x_right_zeros_last_two_columns(self): + var = self._ones() + set_halos_initial(var, [4, 0, 0]) + assert (var[:, :, -1] == 0.0).all() + assert (var[:, :, -2] == 0.0).all() + assert (var[:, :, 0] == 1.0).all() # left side untouched + + def test_no_slip_x_both_zeros_all_x_halos(self): + var = self._ones() + set_halos_initial(var, [6, 0, 0]) + assert (var[:, :, 0] == 0.0).all() + assert (var[:, :, 1] == 0.0).all() + assert (var[:, :, -1] == 0.0).all() + assert (var[:, :, -2] == 0.0).all() + + # Y halos (axis 1, middle dimension) + def test_no_slip_y_left_zeros_rows_0_and_1(self): + var = self._ones() + set_halos_initial(var, [0, 2, 0]) + assert (var[:, 0, :] == 0.0).all() + assert (var[:, 1, :] == 0.0).all() + assert (var[:, -1, :] == 1.0).all() + + def test_no_slip_y_both_zeros_all_y_halos(self): + var = self._ones() + set_halos_initial(var, [0, 6, 0]) + assert (var[:, 0, :] == 0.0).all() + assert (var[:, 1, :] == 0.0).all() + assert (var[:, -1, :] == 0.0).all() + assert (var[:, -2, :] == 0.0).all() + + # Z halos (axis 0, first dimension) + def test_no_slip_z_left_zeros_planes_0_and_1(self): + var = self._ones() + set_halos_initial(var, [0, 0, 2]) + assert (var[0, :, :] == 0.0).all() + assert (var[1, :, :] == 0.0).all() + assert (var[-1, :, :] == 1.0).all() + + def test_no_slip_z_both_zeros_all_z_halos(self): + var = self._ones() + set_halos_initial(var, [0, 0, 6]) + assert (var[0, :, :] == 0.0).all() + assert (var[1, :, :] == 0.0).all() + assert (var[-1, :, :] == 0.0).all() + assert (var[-2, :, :] == 0.0).all() + + def test_interior_never_zeroed_even_with_full_no_slip(self): + var = self._ones() + set_halos_initial(var, [6, 6, 6]) + interior = var[HALO:-HALO, HALO:-HALO, HALO:-HALO] + assert (interior == 1.0).all() + + def test_all_no_slip_zeros_exactly_halo_layers(self): + # With bc=[6,6,6], exactly the HALO=2 outer layers on each face should be 0. + var = self._ones(shape=(10, 10, 10)) + set_halos_initial(var, [6, 6, 6]) + # Inner 6×6×6 block should be 1.0 + assert (var[2:-2, 2:-2, 2:-2] == 1.0).all() + # Corner/edge/face halos should be 0.0 (spot checks) + assert var[0, 0, 0] == 0.0 + assert var[1, 5, 5] == 0.0 + assert var[9, 9, 9] == 0.0 diff --git a/tests/test_integration.py b/tests/test_integration.py new file mode 100644 index 0000000..4e7b3d2 --- /dev/null +++ b/tests/test_integration.py @@ -0,0 +1,187 @@ +# SPDX-License-Identifier: MIT +# +# Copyright 2024-2026 David Krach, Matthias Ruf +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +"""Integration tests: invoke the POREMAPS binary on small generated geometries. + +These tests are skipped unless ``--run-integration`` is passed to pytest. +They also require: + * The binary compiled at ../bin/POREMAPS (run ``make`` in src/). + * ``mpirun`` available on PATH. + +Example +------- + pytest tests/ --run-integration +""" + +import os +import shutil +import subprocess +import tempfile +import textwrap + +import numpy as np +import pytest + +BINARY = os.path.abspath( + os.path.join(os.path.dirname(__file__), "..", "bin", "POREMAPS") +) + + +# --------------------------------------------------------------------------- +# Module-level skip fixture applied to every test in this file +# --------------------------------------------------------------------------- + +@pytest.fixture(autouse=True) +def _check_prerequisites(run_integration): + if not run_integration: + pytest.skip("pass --run-integration to enable binary tests") + if not os.path.isfile(BINARY): + pytest.skip(f"binary not found at {BINARY} — run 'make' in src/ first") + if shutil.which("mpirun") is None: + pytest.skip("mpirun not found on PATH") + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _write_geom(path, nx, ny, nz): + """Write an all-fluid uint8 raw geometry file (x-fastest order).""" + geom = np.zeros((nz, ny, nx), dtype=np.uint8) + geom.ravel(order='C').tofile(path) + + +def _write_inp(path, geom_path, log_path, nx, ny, nz, + bc_method=0, voxelsize=1e-6, max_iter=100, eps=1e-12): + """Write a POREMAPS input parameter file.""" + text = textwrap.dedent(f"""\ + dom_decomposition 0 0 0 + boundary_method {bc_method} + geometry_file_name {geom_path} + size_x_y_z {nx} {ny} {nz} + voxel_size {voxelsize} + max_iter {max_iter} + it_eval 50 + it_write 50 + log_file_name {log_path} + solving_algorithm 2 + eps {eps} + porosity -1.0 + dom_interest 0 0 0 0 0 0 + write_output 0 0 0 0 + """) + with open(path, "w") as f: + f.write(text) + + +def _run(tmpdir, nx=8, ny=8, nz=8, nranks=1, **inp_kwargs): + """Set up geometry + input files, run the binary, return (result, log_path).""" + geom_path = os.path.join(tmpdir, "geom.raw") + log_path = os.path.join(tmpdir, "run.log") + inp_path = os.path.join(tmpdir, "input.inp") + + _write_geom(geom_path, nx, ny, nz) + _write_inp(inp_path, geom_path, log_path, nx, ny, nz, **inp_kwargs) + + result = subprocess.run( + ["mpirun", "-n", str(nranks), BINARY, inp_path], + capture_output=True, + text=True, + cwd=tmpdir, + timeout=120, + ) + return result, log_path + + +# --------------------------------------------------------------------------- +# Tests +# --------------------------------------------------------------------------- + +class TestIntegration: + def test_binary_exits_successfully(self): + with tempfile.TemporaryDirectory() as tmpdir: + result, _ = _run(tmpdir, nx=8, ny=8, nz=8, max_iter=100) + assert result.returncode == 0, ( + f"POREMAPS exited with code {result.returncode}\n" + f"stdout:\n{result.stdout}\nstderr:\n{result.stderr}" + ) + + def test_log_file_is_created(self): + # The solver only writes the log after iteration 1000; use max_iter > 1000. + with tempfile.TemporaryDirectory() as tmpdir: + result, log_path = _run(tmpdir, nx=8, ny=8, nz=8, max_iter=1100) + assert result.returncode == 0 + assert os.path.isfile(log_path), "Log file was not created by the solver." + + def test_log_file_contains_data(self): + with tempfile.TemporaryDirectory() as tmpdir: + result, log_path = _run(tmpdir, nx=8, ny=8, nz=8, max_iter=1100) + assert result.returncode == 0 + with open(log_path) as f: + content = f.read() + assert len(content.strip()) > 0, "Log file is empty." + + def test_log_contains_numeric_values(self): + with tempfile.TemporaryDirectory() as tmpdir: + result, log_path = _run(tmpdir, nx=8, ny=8, nz=8, max_iter=1100) + assert result.returncode == 0 + with open(log_path) as f: + lines = [ln.strip() for ln in f if ln.strip()] + assert lines, "Log file has no non-empty lines." + # At least the last data line should be parseable as floats. + # Log format is comma-separated: "iter, val, val, ..." + data_lines = [ln for ln in lines if not ln.startswith('#')] + assert data_lines, "Log file has no data lines (only a header)." + try: + values = [float(v) for v in data_lines[-1].split(',')] + except ValueError: + pytest.fail(f"Last log line is not numeric: {data_lines[-1]!r}") + assert len(values) > 0 + + def test_two_mpi_ranks(self): + """Solver must complete without error when split across 2 MPI ranks.""" + with tempfile.TemporaryDirectory() as tmpdir: + result, _ = _run(tmpdir, nx=8, ny=8, nz=8, nranks=2, max_iter=50) + assert result.returncode == 0, ( + f"2-rank run failed\nstdout:\n{result.stdout}\nstderr:\n{result.stderr}" + ) + + def test_permeability_positive_all_fluid(self): + """The whole-domain z-permeability (wk33) in the log must be positive.""" + with tempfile.TemporaryDirectory() as tmpdir: + result, log_path = _run( + tmpdir, nx=8, ny=8, nz=10, + max_iter=1100, eps=1e-8, voxelsize=1e-5, + ) + assert result.returncode == 0 + with open(log_path) as f: + lines = [ln.strip() for ln in f if ln.strip()] + assert lines + data_lines = [ln for ln in lines if not ln.startswith('#')] + assert data_lines, "Log file has no data lines (only a header)." + try: + values = [float(v) for v in data_lines[-1].split(',')] + except ValueError: + pytest.skip("Could not parse log file format — check log output") + assert any(v > 0 for v in values), ( + f"Expected at least one positive value in log line: {data_lines[-1]}" + ) diff --git a/tests/test_solver.py b/tests/test_solver.py new file mode 100644 index 0000000..37ba205 --- /dev/null +++ b/tests/test_solver.py @@ -0,0 +1,164 @@ +# SPDX-License-Identifier: MIT +# +# Copyright 2024-2026 David Krach, Matthias Ruf +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +"""Tests for split_number and get_2nd_derivation (src/solver.cc).""" + +import pytest + +from helpers import split_number, get_2nd_derivation + + +class TestSplitNumber: + def test_zero(self): + assert split_number(0) == (0, 0, 0) + + def test_single_digit_z(self): + assert split_number(6) == (0, 0, 6) + + def test_single_digit_y(self): + assert split_number(30) == (0, 3, 0) + + def test_single_digit_x(self): + assert split_number(200) == (2, 0, 0) + + def test_all_ones(self): + assert split_number(111) == (1, 1, 1) + + def test_all_sixes(self): + assert split_number(666) == (6, 6, 6) + + def test_mixed(self): + assert split_number(246) == (2, 4, 6) + assert split_number(135) == (1, 3, 5) + + def test_round_trip(self): + """Encoding cases as x*100+y*10+z and decoding must recover originals.""" + for x in range(0, 7): + for y in range(0, 7): + for z in range(0, 7): + code = x * 100 + y * 10 + z + assert split_number(code) == (x, y, z) + + def test_only_case_z(self): + cx, cy, cz = split_number(1) + assert (cx, cy, cz) == (0, 0, 1) + + def test_only_case_x(self): + cx, cy, cz = split_number(600) + assert (cx, cy, cz) == (6, 0, 0) + + +class TestGet2ndDerivation: + """Verify each neighbourhood case against the finite-difference formulae.""" + + # ------------------------------------------------------------------ + # Case 1: both neighbours are walls → -8 * cvel + # ------------------------------------------------------------------ + def test_case1_unit_centre(self): + assert get_2nd_derivation(1, 0.0, 0.0, 1.0, 0.0, 0.0) == pytest.approx(-8.0) + + def test_case1_zero_centre(self): + assert get_2nd_derivation(1, 5.0, 5.0, 0.0, 5.0, 5.0) == pytest.approx(0.0) + + def test_case1_ignores_neighbours(self): + # llvel/lvel/rvel/rrvel are irrelevant for case 1 + r1 = get_2nd_derivation(1, 99.0, 99.0, 2.0, 99.0, 99.0) + assert r1 == pytest.approx(-16.0) + + # ------------------------------------------------------------------ + # Case 2: wall on left, one fluid cell to the right + # ------------------------------------------------------------------ + def test_case2_formula(self): + result = get_2nd_derivation(2, 0.0, 0.0, 1.0, 3.0, 0.0) + expected = 8.0 / 3.0 * 3.0 - 16.0 / 3.0 * 1.0 + assert result == pytest.approx(expected) + + def test_case2_zero_result(self): + # 8/3 * rvel == 16/3 * cvel → rvel = 2 * cvel + assert get_2nd_derivation(2, 0.0, 0.0, 1.0, 2.0, 0.0) == pytest.approx(0.0) + + # ------------------------------------------------------------------ + # Case 3: wall on right, one fluid cell to the left + # ------------------------------------------------------------------ + def test_case3_formula(self): + result = get_2nd_derivation(3, 0.0, 3.0, 1.0, 0.0, 0.0) + expected = 8.0 / 3.0 * 3.0 - 16.0 / 3.0 * 1.0 + assert result == pytest.approx(expected) + + def test_case2_case3_mirror_symmetry(self): + """Swapping left and right with equal values gives the same result.""" + r2 = get_2nd_derivation(2, 0.0, 0.0, 2.0, 4.0, 0.0) + r3 = get_2nd_derivation(3, 0.0, 4.0, 2.0, 0.0, 0.0) + assert r2 == pytest.approx(r3) + + # ------------------------------------------------------------------ + # Case 4: wall on left, two fluid cells to the right + # ------------------------------------------------------------------ + def test_case4_formula(self): + result = get_2nd_derivation(4, 0.0, 0.0, 1.0, 2.0, 3.0) + expected = 2.0 * 2.0 - 5.0 * 1.0 - 1.0 / 5.0 * 3.0 + assert result == pytest.approx(expected) + + # ------------------------------------------------------------------ + # Case 5: wall on right, two fluid cells to the left + # ------------------------------------------------------------------ + def test_case5_formula(self): + result = get_2nd_derivation(5, 3.0, 2.0, 1.0, 0.0, 0.0) + expected = 2.0 * 2.0 - 1.0 / 5.0 * 3.0 - 5.0 * 1.0 + assert result == pytest.approx(expected) + + def test_case4_case5_mirror_symmetry(self): + """Mirroring stencil values between cases 4 and 5 gives equal results.""" + r4 = get_2nd_derivation(4, 0.0, 0.0, 2.0, 4.0, 6.0) + r5 = get_2nd_derivation(5, 6.0, 4.0, 2.0, 0.0, 0.0) + assert r4 == pytest.approx(r5) + + # ------------------------------------------------------------------ + # Case 6: fluid on both sides — standard 3-point central difference + # ------------------------------------------------------------------ + def test_case6_quadratic_field(self): + # v(x) = x²: v(-1)=1, v(0)=0, v(1)=1 → d²v/dx² = 2 + assert get_2nd_derivation(6, 0.0, 1.0, 0.0, 1.0, 0.0) == pytest.approx(2.0) + + def test_case6_linear_field_zero_laplacian(self): + # v(x) = x: v(-1)=-1, v(0)=0, v(1)=1 → d²v/dx² = 0 + assert get_2nd_derivation(6, 0.0, -1.0, 0.0, 1.0, 0.0) == pytest.approx(0.0) + + def test_case6_constant_field_zero_laplacian(self): + assert get_2nd_derivation(6, 5.0, 5.0, 5.0, 5.0, 5.0) == pytest.approx(0.0) + + def test_case6_ignores_second_neighbours(self): + # llvel and rrvel are unused in case 6 + r1 = get_2nd_derivation(6, 99.0, 1.0, 2.0, 3.0, 99.0) + r2 = get_2nd_derivation(6, -1.0, 1.0, 2.0, 3.0, -1.0) + assert r1 == pytest.approx(r2) + + # ------------------------------------------------------------------ + # Error handling + # ------------------------------------------------------------------ + def test_invalid_case_zero_raises(self): + with pytest.raises(ValueError): + get_2nd_derivation(0, 1.0, 1.0, 1.0, 1.0, 1.0) + + def test_invalid_case_seven_raises(self): + with pytest.raises(ValueError): + get_2nd_derivation(7, 1.0, 1.0, 1.0, 1.0, 1.0)