Longitudinal_velocity_gradient_moments.ipynb

Kernel: Python 3 (CoCalc)

import sympy as sp
from sympy import Function, Integer
from sympy.functions.combinatorial.factorials import factorial, factorial2
from collections import Counter

class TraceS(Function):
    nargs = 1
    @classmethod
    def eval(cls, k):
        # keep symbolic (no automatic evaluation)
        return None


class AngleAverage(Function):
    nargs = 1
    def _sympystr(self, printer):
        return f"<{printer.doprint(self.args[0])}>"


def integer_partitions(n: int, max_part: int | None = None):
    """Yield integer partitions of n as nonincreasing lists."""
    if max_part is None or max_part > n:
        max_part = n
    if n == 0:
        yield []
        return
    for first in range(min(max_part, n), 0, -1):
        for rest in integer_partitions(n - first, first):
            yield [first] + rest


# =========================================================
# MomentPolynomial[n] (Bell-polynomial-like combinatorial sum)
# =========================================================

def moment_polynomial(n: int):
    if n < 0:
        raise ValueError("n must be nonnegative")
    if n == 0:
        return sp.Integer(1)

    total = sp.Integer(0)

    for part in integer_partitions(n):
        tally = Counter(part)  # k -> multiplicity m
        klist = list(tally.keys())
        mklist = [tally[k] for k in klist]

        # coeff = n! / ( Π(k^m) Π(m!) ) * 2^(n - Σ m)
        denom = sp.Integer(1)
        for k, m in zip(klist, mklist):
            denom *= (sp.Integer(k) ** m) * factorial(m)

        coeff = factorial(n) / denom
        coeff *= sp.Integer(2) ** (n - sum(mklist))

        term = coeff
        for k, m in zip(klist, mklist):
            term *= TraceS(Integer(k)) ** m

        total += term

    return sp.simplify(total)


# =========================================================
# IsotropicMoment[n] with options:
#   incompressible: set TraceS(1)=0
#   reduce_ch: Cayley–Hamilton reduction for 3x3 to express Tr(S^k), k>3
# =========================================================

def isotropic_moment(n: int, *, incompressible: bool = False, reduce_ch: bool = False):
    if n <= 0:
        raise ValueError("n must be positive")

    Cn = sp.Rational(1, 1) / factorial2(2 * n + 1)  # 1/(2n+1)!!
    poly = moment_polynomial(n)
    expr = poly

    if reduce_ch:
        # elementary symmetric polynomials in eigenvalues:
        e1 = TraceS(1)
        e2 = -(e1**2 - TraceS(2)) / 2
        e3 = (e1**3 - 3 * TraceS(2) * e1 + 2 * TraceS(3)) / 6

        # power sums p[k] = Tr(S^k), CH recurrence:
        p_cache = {1: TraceS(1), 2: TraceS(2), 3: TraceS(3)}

        def p(k: int):
            k = int(k)
            if k in p_cache:
                return p_cache[k]
            val = sp.expand(e1 * p(k - 1) + e2 * p(k - 2) + e3 * p(k - 3))
            p_cache[k] = val
            return val

        # replace TraceS(k) for k>3
        ks = [int(a.args[0]) for a in poly.atoms(TraceS)]
        maxK = max(ks) if ks else 0
        subs = {TraceS(k): p(k) for k in range(4, maxK + 1)}
        expr = sp.expand(poly.subs(subs))

    if incompressible:
        expr = sp.expand(expr.subs({TraceS(1): sp.Integer(0)}))

    return sp.simplify(Cn * expr)


# =========================================================
# Formatting: wrap each monomial in <...>
# =========================================================

def format_moment(expr):
    expanded = sp.expand(expr)
    terms = expanded.as_ordered_terms() if isinstance(expanded, sp.Add) else [expanded]

    out = sp.Integer(0)
    for term in terms:
        # "coeff" = set all TraceS(*) -> 1
        subs = {a: 1 for a in term.atoms(TraceS)}
        coeff = sp.simplify(term.subs(subs))
        rest = sp.simplify(term / coeff) if coeff != 0 else term
        out += coeff * AngleAverage(rest)

    return sp.factor_terms(out)


# =========================================================
# Public wrapper: A11nmoment[n]
# default: compressible, ReduceCH ON
# =========================================================

def A11nmoment(n: int, *, incompressible: bool = False, reduce_ch: bool = True):
    raw = isotropic_moment(n, incompressible=incompressible, reduce_ch=reduce_ch)
    return format_moment(raw)


# ----------------- Examples -----------------

if __name__ == "__main__":
    # n=3, compressible
    print("n=3 compressible:")
    print(A11nmoment(3, incompressible=False))
    # -> ( <TraceS(1)**3> + 6<TraceS(1)TraceS(2)> + 8<TraceS(3)> ) / 105

    # n=4, compressible (CH reduction ON by default)
    print("\nn=4 compressible (ReduceCH=True):")
    print(A11nmoment(4))

    # n=4, incompressible
    print("\nn=4 incompressible:")
    print(A11nmoment(4, incompressible=True))

    # n=4, incompressible
    print("\nn=4 incompressible:")
    print(A11nmoment(4, incompressible=True))

n=3 compressible:
(6*<TraceS(1)*TraceS(2)> + <TraceS(1)**3> + 8*<TraceS(3)>)/105

n=4 compressible (ReduceCH=True):
(32*<TraceS(1)*TraceS(3)> - 12*<TraceS(1)**2*TraceS(2)> + 3*<TraceS(1)**4> + 12*<TraceS(2)**2>)/315

n=4 incompressible:
4*<TraceS(2)**2>/105