[JoaoBentoNeto]

Summary

This Pull Request introduces improvements to the lbpm_greyscale_simulator routine, including bug fixes, corrections to the physical formulation, and new features.

---

Bug Fixes

• fix: correct absolute permeability calculation: Corrected the calculation of Darcy velocity used to determine absolute permeability, ensuring proper averaging in greyscale media with partially porous voxels.

The absolute permeability is calculated by:

$$\Large \begin{aligned} K = \frac{\langle u \rangle \nu}{g - \frac{\Delta p}{L}} \end{aligned}$$

where $\langle u\rangle$ is the domain-averaged Darcy velocity, $\nu$ is the kinematic viscosity, $g$ is the body force, $\Delta p$ is the pressure difference between the inlet and outlet, and $L$ is the domain length in the streamwise direction.

By definition, the Darcy velocity is obtained by averaging the local velocity over the total volume of the porous medium:

$$\Large \begin{aligned} \langle u \rangle &= \frac{1}{V_{\mathrm{total}}} \int_{V_{\mathrm{total}}} u ~\mathrm{d}V \\ &= \frac{1}{V_{\text{total}}} \sum_i^{N_{\text{total}}} \left(\int_{V_{\text{i}}} u ~\mathrm{d}V \right) \\ &= \frac{1}{V_{\mathrm{total}}} \sum_{i=1}^{N_{\mathrm{total}}} \left(\langle u\rangle_i V_{\mathrm{voxel}}\right)\\ &= \frac{1}{N_{\mathrm{total}}} \sum_{i=1}^{N_{\mathrm{total}}}\langle u\rangle_i \end{aligned}$$

where $\langle u\rangle_i$ is the voxel-averaged Darcy velocity, $V_{\mathrm{voxel}}$ is the voxel volume, and $N_{\mathrm{total}}$ is the total number of voxels.

The previous implementation attempted to calculate the Darcy velocity by averaging the voxel velocities over the pore voxels and subsequently multiplying the result by the total porosity:

$$\Large \begin{aligned} \langle u \rangle_{\mathrm{previous}} &= \left( \frac{1}{N_{\mathrm{pores}}} \sum_{i=1}^{N_{\mathrm{total}}}\langle u\rangle_i \right)\phi_{\mathrm{total}}\\ &= \left(\frac{1}{N_{\text{pores}}} \sum_i^{N_{\text{total}}} \langle u \rangle_i \right) \frac{\sum_i^{N_{\text{total}}} \phi_i}{N_{\text{total}}} \\ &= \left( \frac{1}{N_{\mathrm{total}}} \sum_{i=1}^{N_{\mathrm{total}}}\langle u\rangle_i \right) \left( \frac{\sum_{i=1}^{N_{\mathrm{total}}}\phi_i} {N_{\mathrm{pores}}} \right) \end{aligned}$$

This expression is equivalent to the volume-averaged definition only when all non-solid voxels have unit porosity ($\phi_i=1$). Therefore, the previous approach is generally incorrect for greyscale media, where voxels may have intermediate porosity values ($0<\phi_i\leq1$).

Therefore, the calculation has been replaced with the direct volume average of the voxel Darcy velocities, eliminating the incorrect porosity correction.

• fix: corrected the 'GreyPorosity' calculation: Excluded communication halo cells from the porosity calculation by moving the relevant operations inside the conditional statement that checks whether a voxel doesn't belongs to the halo. Also removed unnecessary variables.

• fix: corrected the permeability unit for the raw input file: Added a unit conversion for permeability values read from raw input files, ensuring consistency with the label-based property assignment. Both input methods now uses permeability in squared micrometers.

• fix: resolve memory leak

• fix: replace raw arrays with std::vector for file assignments: Replaced raw arrays with std::vector for file-based property assignments to be consistent with the modifiction used to solve the memory leak..

Features

• feature: add boolean toggle for Forchheimer non-linear drag term: Introduced a Boolean Forchheimer flag to control the inclusion of the non-linear drag force in the flow calculations.

• feature: add VTK (.vti/.pvti) output support

• improvement: add exception for undeclared voxel labels: Added an exception to handle voxels containing labels that have not been declared in the input configuration.

• improvement: enhanced log output: Enhanced the simulation log to display the current time step, convergence error, and absolute permeability in both square micrometers (µm²) and millidarcies (mD).

• feature: add ViscosityRatio parameter: Clarified the nomenclature of viscosity-related variables, distinguishing the parameters associated with the Darcy drag term (viscosity, rlx, and tau) from those associated with the Brinkman effective viscosity (viscosity_eff, rlx_eff, and tau_eff). Introduced the floating-point parameter ViscosityRatio to control the ratio between effective and dynamic viscosities.

Refactoring and Style

• refactor: enhance error handling and validation checks: Consolidated the length checks for ComponentLabels, PorosityList, and PermeabilityList into a single validation block at the beginning of the routine. Moved the physical property validation checks ($0<\phi\leq1$ and $k>0$) from the main domain loop to the input-list processing loop.

• refactor: unify porosity and permeability assignment loops: Combined the porosity and permeability assignment loops into a single loop to improve readability and reduce redundant iterations.

• style: apply clang-format to greyscale files: Applied the project's C/C++ formatting configuration using clang-format to ensure consistency with the rest of the codebase.