Introduce FrictionFactorModel protocol with built‑in implementations and user extensibility

[dad616610]

Summary

This PR replaces the hard‑coded friction factor calculations with a FrictionFactorModel protocol, making it easy to use built‑in models (Nikuradse, Swamee‑Jain, Colebrook) or supply custom ones. A RegimeAwareFrictionFactorModel can switch between models based on the Reynolds number. The old string‑based interface is still supported.

Note for maintainers: Please squash this PR when merging.

Key changes

• New FrictionFactorModel protocol – any callable implementing compute_lambda_and_dlambda_dm() works with the pipe flow solver.

• Unified computation – calc_lambda and calc_der_lambda are merged into a single method that returns both $\lambda$ and $\frac{\mathrm{d}\lambda}{\mathrm{d}\dot{m}}$, reusing intermediate results and keeping the interface simple.

• Built‑in models – Nikuradse, SwameeJain, Colebrook, and RegimeAwareFrictionFactorModel.

• Unified Nikuradse formulations – the previously separate compressible and incompressible Nikuradse functions were mathematically equivalent¹; now only one implementation remains. Though I had to relax the tolerance in the tests (see commit 9960d71).

• Consistent pre‑filtering – zero‑flow and zero‑length pipes are now filtered once before any model is called, instead of being handled inside each model.

• Removed Numba support from friction factor models (rationale²: easier custom numba‑models by users, and no clean way to integrate use_numba automatically). This can be revisited later.

• Documentation – new pipeflow/friction_factor_models page; the pipe component page now refers to it.

• Backward compatibility – string names like "nikuradse" are converted to the appropriate model instance internally.

Known limitations / follow‑ups

• The Nikuradse derivative currently uses $|m|$ instead of the physically correct odd form $-\frac{64}{Re \cdot m}$. Fixing it slightly degrades convergence; this will be addressed separately.

• The Colebrook error message still mentions max_iter_colebrook; a follow‑up can adjust it for the new API.

• Numba support for friction models can be re‑added if performance demands it; custom numba‑based models can be used now.

Impact

• Enables users to implement and use any friction factor model without modifying pandapipes internals.

• Potentially closes add friction model "Hofer" #683 and Implemented new friction factor calculation spanning entire flow regime #767 as these friction factor models can now be implemented from user-side.

---

Footnotes

Proof of mathematical equivalence

$$ \begin{align} &\frac{1}{\left( -2\log_{10}\left( \frac{k}{3.71 D} \right) \right)^2} \\ =& \frac{1}{\left( -2\left( \log_{10}\left( \frac{k}{D} \right) - \log_{10}\left( 3.71 \right) \right) \right)^2} \\ =& \frac{1}{\left( -2 \log_{10}\left( \frac{k}{D} \right) + 2 \log_{10}\left( 3.71 \right) \right)^2} \\ =& \frac{1}{\left( 2 \log_{10}\left( \frac{D}{k} \right) + 2 \log_{10}\left( 3.71 \right) \right)^2} \\ =& \frac{1}{\left( 2 \log_{10}\left( \frac{D}{k} \right) + 1.13875 \right)^2} \\ \approx& \frac{1}{\left( 2 \log_{10}\left( \frac{D}{k} \right) + 1.14 \right)^2} \end{align} $$

On numba support removal

Numba starts beating pure-numpy functions with >10k elements -- below they're almost equal. I ran the benchmark tests for different sizes of elements. Each function was run 100 times per round, with 7 rounds. The best (minimum) time across the 7 rounds is shown in seconds.

Nikuradse results

n = 100
        nikuradse: 0.0004
        nikuradse_nb: 0.0001
        nikuradse_nb_par: 0.0011
n = 1_000
        nikuradse: 0.0014
        nikuradse_nb: 0.0010
        nikuradse_nb_par: 0.0010
n = 10_000
        nikuradse: 0.0078
        nikuradse_nb: 0.0061
        nikuradse_nb_par: 0.0013
n = 100_000
        nikuradse: 0.0843
        nikuradse_nb: 0.0606
        nikuradse_nb_par: 0.0037
n = 1_000_000
        nikuradse: 1.2289
        nikuradse_nb: 0.6955
        nikuradse_nb_par: 0.1360

Swamee-Jain results

n = 100
        swameejain: 0.0008
        swameejain_nb: 0.0004
        swameejain_nb_par: 0.0003
n = 1_000
        swameejain: 0.0020
        swameejain_nb: 0.0024
        swameejain_nb_par: 0.0024
n = 10_000
        swameejain: 0.0154
        swameejain_nb: 0.0231
        swameejain_nb_par: 0.0231
n = 100_000
        swameejain: 0.2044
        swameejain_nb: 0.2314
        swameejain_nb_par: 0.2307
n = 1_000_000
        swameejain: 2.3740
        swameejain_nb: 2.3845
        swameejain_nb_par: 2.3856

Colebrook results

n = 100
        colebrook: 0.0068
        colebrook_nb: 0.0014
        colebrook_nb_par: 0.0012
n = 1_000
        colebrook: 0.0181
        colebrook_nb: 0.0136
        colebrook_nb_par: 0.0016
n = 10_000
        colebrook: 0.1417
        colebrook_nb: 0.1378
        colebrook_nb_par: 0.0059
n = 100_000
        colebrook: 1.3551
        colebrook_nb: 1.3814
        colebrook_nb_par: 0.0496
n = 1_000_000
        colebrook: 27.5238
        colebrook_nb: 13.8841
        colebrook_nb_par: 0.584

Benchmarking code

def nikuradse(k_over_D, re):
    return 64 / re + 1 / (-2 * np.log10(k_over_D / 3.71)) ** 2


@nb.njit
def nikuradse_nb(k_over_D, re):
    n = re.size
    res = np.empty(n, dtype=np.float64)
    for i in range(n):
        res[i] = 64 / re[i] + 1 / (-2 * np.log10(k_over_D[i] / 3.71)) ** 2
    return res


@nb.njit(parallel=True)
def nikuradse_nb_par(k_over_D, re):
    n = re.size
    res = np.empty(n, dtype=np.float64)
    for i in nb.prange(n):
        res[i] = 64 / re[i] + 1 / (-2 * np.log10(k_over_D[i] / 3.71)) ** 2
    return res


def swameejain(k_over_D, re):
    inv_re_09 = 1 / re**0.9
    inner_log_term = k_over_D / 3.7 + 5.74 * inv_re_09
    log_term = np.log(inner_log_term)
    log_squared = log_term * log_term

    # a = 0.25 * ln(10)**2
    a = 1.325474527619599502640416597148504422899
    return a / log_squared


@nb.njit
def swameejain_nb(k_over_D, re):
    n = re.size
    res = np.empty(n, dtype=np.float64)
    a = 1.325474527619599502640416597148504422899

    for i in range(n):
        inv_re_09 = 1 / re[i] ** 0.9
        inner_log_term = k_over_D[i] / 3.7 + 5.74 * inv_re_09
        log_term = np.log(inner_log_term)
        log_squared = log_term * log_term

        res[i] = a / log_squared
    return res


@nb.njit
def swameejain_nb_par(k_over_D, re):
    n = re.size
    res = np.empty(n, dtype=np.float64)
    a = 1.325474527619599502640416597148504422899

    for i in nb.prange(n):
        inv_re_09 = 1 / re[i] ** 0.9
        inner_log_term = k_over_D[i] / 3.7 + 5.74 * inv_re_09
        log_term = np.log(inner_log_term)
        log_squared = log_term * log_term

        res[i] = a / log_squared
    return res


max_iter = 10000
tolerance = 1e-4


def colebrook(k_over_D, re):
    lambda_prev = 1 / (-2 * np.log10(k_over_D / 3.71)) ** 2
    lambda_curr = lambda_prev

    a = k_over_D / 3.71
    b = 2.51 / re
    # 1 / ln(10)
    inv_ln10 = 0.4342944819032518276511289189166050822944
    for _ in range(max_iter):
        inv_lambda_sqrt = 1 / np.sqrt(lambda_curr)
        inner_log_term = a + b * inv_lambda_sqrt
        cubed_inv_lambda_sqrt = inv_lambda_sqrt**3

        f = inv_lambda_sqrt + 2 * np.log10(inner_log_term)
        df = (
            -0.5 * cubed_inv_lambda_sqrt
            - b * cubed_inv_lambda_sqrt * inv_ln10 / inner_log_term
        )

        lambda_curr = lambda_prev - f / df
        if np.all(np.abs(lambda_curr - lambda_prev) < tolerance):
            break
        lambda_prev = lambda_curr
    return lambda_curr


@nb.njit
def colebrook_nb(k_over_D, re):
    n = re.size
    res = np.empty(n, dtype=np.float64)

    # 1 / ln(10)
    inv_ln10 = 0.4342944819032518276511289189166050822944
    for i in range(n):
        lambda_prev = 1 / (-2 * np.log10(k_over_D[i] / 3.71)) ** 2
        lambda_curr = lambda_prev

        a = k_over_D[i] / 3.71
        b = 2.51 / re[i]

        for _ in range(max_iter):
            inv_lambda_sqrt = 1 / np.sqrt(lambda_curr)
            inner_log_term = a + b * inv_lambda_sqrt
            cubed_inv_lambda_sqrt = inv_lambda_sqrt**3

            f = inv_lambda_sqrt + 2 * np.log10(inner_log_term)
            df = (
                -0.5 * cubed_inv_lambda_sqrt
                - b * cubed_inv_lambda_sqrt * inv_ln10 / inner_log_term
            )

            lambda_curr = lambda_prev - f / df
            if np.abs(lambda_curr - lambda_prev) < tolerance:
                res[i] = lambda_curr
                break
            lambda_prev = lambda_curr

    return res


@nb.njit(parallel=True)
def colebrook_nb_par(k_over_D, re):
    n = re.size
    res = np.empty(n, dtype=np.float64)

    # 1 / ln(10)
    inv_ln10 = 0.4342944819032518276511289189166050822944
    for i in nb.prange(n):
        lambda_prev = 1 / (-2 * np.log10(k_over_D[i] / 3.71)) ** 2
        lambda_curr = lambda_prev

        a = k_over_D[i] / 3.71
        b = 2.51 / re[i]

        for _ in range(max_iter):
            inv_lambda_sqrt = 1 / np.sqrt(lambda_curr)
            inner_log_term = a + b * inv_lambda_sqrt
            cubed_inv_lambda_sqrt = inv_lambda_sqrt**3

            f = inv_lambda_sqrt + 2 * np.log10(inner_log_term)
            df = (
                -0.5 * cubed_inv_lambda_sqrt
                - b * cubed_inv_lambda_sqrt * inv_ln10 / inner_log_term
            )

            lambda_curr = lambda_prev - f / df
            if np.abs(lambda_curr - lambda_prev) < tolerance:
                res[i] = lambda_curr
                break
            lambda_prev = lambda_curr

    return res


def bench_fs(sizes, fs):
    rng = np.random.default_rng(seed=42)
    for n in sizes:
        k = rng.uniform(1e-5, 1e-2, size=n)
        d = rng.uniform(20, 800, size=n)
        k_over_D = k / d
        re = rng.uniform(2500, 8000, size=n)

        print(f"{n = :_}")
        for f in fs:
            res = r(
                partial(f, k_over_D, re),
                repeat=7,
                number=100,
                globals={**globals(), **locals()},
            )
            print(f"\t{f.__name__}: {min(res):.4f}")


def bench(sizes):
    bench_fs(sizes, (nikuradse, nikuradse_nb, nikuradse_nb_par))
    print()
    bench_fs(sizes, (swameejain, swameejain_nb, swameejain_nb_par))
    print()
    bench_fs(sizes, (colebrook, colebrook_nb, colebrook_nb_par))


if __name__ == "__main__":
    sizes = [100, 1000, 10000, 100_000, 1_000_000]
    bench(sizes)

[dad616610]

Ah, I forgot to mention #754 as the main reason behind this PR. Now a composite friction factor model can be created without modifying the pandapipes internals and now we have a built-in RegimeAwareFrictionFactorModel that distinguishes between regimes: laminar, transient and turbulent. If that's too much, I guess either users can create their own FrictionFactorModel just for laminar and turbulent, or we might alter the current RegimeAwareFrictionFactorModel to respect only 2 regimes, instead of 3.

I've also improved the performance of Colebrook (up to ~1.76x in these benchmarks, see: 859f7f4).

Benchmarking results

colebrook uses a**3, while colebrook_fast uses a * a * a. Each function was run 100 times per round, with 7 rounds. The best (minimum) time across the 7 rounds is shown in seconds.

n = 10
        colebrook: 0.0070
        colebrook_fast: 0.0068, speedup: 1.0223
n = 100
        colebrook: 0.0081
        colebrook_fast: 0.0073, speedup: 1.1147
n = 1_000
        colebrook: 0.0183
        colebrook_fast: 0.0120, speedup: 1.5202
n = 10_000
        colebrook: 0.1578
        colebrook_fast: 0.0894, speedup: 1.7643
n = 100_000
        colebrook: 1.4826
        colebrook_fast: 0.8518, speedup: 1.7406
n = 300_000
        colebrook: 8.2176
        colebrook_fast: 6.3043, speedup: 1.3035

Code to reproduce

from functools import partial
from timeit import repeat as r

import numpy as np

max_iter = 10_000
tolerance = 1e-4


def colebrook(k_over_D, re):
    lambda_prev = 1 / (-2 * np.log10(k_over_D / 3.71)) ** 2
    lambda_curr = None

    a = k_over_D / 3.71
    b = 2.51 / re
    # 1 / ln(10)
    inv_ln10 = 0.4342944819032518276511289189166050822944
    for _ in range(max_iter):
        inv_lambda_sqrt = 1 / np.sqrt(lambda_prev)
        inner_log_term = a + b * inv_lambda_sqrt
        cubed_inv_lambda_sqrt = inv_lambda_sqrt**3

        f = inv_lambda_sqrt + 2 * np.log10(inner_log_term)
        df = (
            -0.5 * cubed_inv_lambda_sqrt
            - b * cubed_inv_lambda_sqrt * inv_ln10 / inner_log_term
        )
        lambda_curr = lambda_prev - f / df
        if np.all(np.abs(lambda_curr - lambda_prev) < tolerance):
            break
        lambda_prev = lambda_curr
    return lambda_curr


def colebrook_fast(k_over_D, re):
    lambda_prev = 1 / (-2 * np.log10(k_over_D / 3.71)) ** 2
    lambda_curr = None

    a = k_over_D / 3.71
    b = 2.51 / re
    # 1 / ln(10)
    inv_ln10 = 0.4342944819032518276511289189166050822944
    for _ in range(max_iter):
        inv_lambda_sqrt = 1 / np.sqrt(lambda_prev)
        inner_log_term = a + b * inv_lambda_sqrt
        cubed_inv_lambda_sqrt = inv_lambda_sqrt * inv_lambda_sqrt * inv_lambda_sqrt

        f = inv_lambda_sqrt + 2 * np.log10(inner_log_term)
        df = (
            -0.5 * cubed_inv_lambda_sqrt
            - b * cubed_inv_lambda_sqrt * inv_ln10 / inner_log_term
        )
        lambda_curr = lambda_prev - f / df
        if np.all(np.abs(lambda_curr - lambda_prev) < tolerance):
            break
        lambda_prev = lambda_curr
    return lambda_curr


def bench(sizes):
    fs = (colebrook, colebrook_fast)
    rng = np.random.default_rng(seed=42)
    for n in sizes:
        k = rng.uniform(1e-5, 1e-2, size=n)
        d = rng.uniform(20, 800, size=n)
        k_over_D = k / d
        re = rng.uniform(2500, 8000, size=n)

        print(f"{n = :_}")
        for i, f in enumerate(fs):
            res = r(
                partial(f, k_over_D, re),
                repeat=7,
                number=100,
                globals={**globals(), **locals()},
            )
            if i == 0:
                baseline = min(res)
                print(f"\t{f.__name__}: {min(res):.4f}")
            else:
                print(
                    f"\t{f.__name__}: {min(res):.4f}, speedup: {baseline / min(res):.4f}"
                )


bench(sizes=np.array([10, 100, 1_000, 10_000, 100_000, 300_000]))