proteus.outgas.compaction¶

def darcy_velocity(
    porosity: float | np.ndarray,
    grain_size: float,
    delta_rho: float,
    gravity: float,
    melt_visc: float,
    *,
    buoyancy: str,
) -> float | np.ndarray:
    """Melt percolation speed relative to the matrix [m s-1].

    ``w_D = |delta_rho| g F(phi) / eta_melt``, the solver's separation velocity
    under ``buoyancy`` with the melt viscosity as drag. See the module docstring.

    Parameters
    ----------
    porosity : float or ndarray
        Volume fraction of melt [1].
    grain_size : float
        Crystal grain size [m], positive.
    delta_rho : float
        Solid minus melt density [kg m-3]; its magnitude is used.
    gravity : float
        Gravitational acceleration [m s-2].
    melt_visc : float
        Melt viscosity [Pa s], positive.
    buoyancy : str
        The solver's ``separation_buoyancy``, ``'mixture'`` or ``'solid'``.

    Raises
    ------
    ValueError
        If the melt viscosity is not positive, or ``buoyancy`` names no mode.
    RuntimeError
        If the speed is negative or not finite.
    """
    if melt_visc <= 0.0:
        raise ValueError(f'melt viscosity must be positive, got {melt_visc!r}')
    f = mobility_function(np.asarray(porosity, dtype=float), grain_size, buoyancy=buoyancy)
    w = abs(delta_rho) * gravity * f / melt_visc
    if not np.all(np.isfinite(w) & (w >= 0.0)):
        raise RuntimeError(
            f'percolation speed is negative or not finite at porosity {porosity!r}, gravity '
            f'{gravity!r} m/s2, density contrast {delta_rho!r} kg/m3 and melt viscosity '
            f'{melt_visc!r} Pa s'
        )
    return float(w) if np.ndim(porosity) == 0 else w

def drainage_integral(
    porosity_top: float,
    t_res: float,
    thickness: float,
    grain_size: float,
    delta_rho: float,
    gravity: float,
    melt_visc: float,
    mush_visc: float,
    *,
    buoyancy: str,
) -> tuple[float, float, float]:
    """Porosity a parcel retains after crossing the front [1].

    Integrates ``dphi/dt = -phi / max(L/w_D(phi), tau_s)`` from ``porosity_top``
    over the residence time, with ``w_D`` under ``buoyancy``. ``w_D`` falls as
    ``phi^2`` down to a porosity of about 3e-5. Below about 3e-6 the mobility tends to
    the Stokes residue of the tanh blend, a numerical floor, so the decay turns
    exponential, and after a residence far longer than a front's the porosity ends within
    the integrator's tolerance of zero. The result is clipped to ``[0, porosity_top]``;
    an end below 0 by more than the allowance 10 atol, or above ``porosity_top`` by more
    than 10 atol + rtol ``porosity_top``, raises ``RuntimeError``, as does a non-positive
    percolation time.

    Parameters
    ----------
    porosity_top : float
        Melt volume fraction where the parcel enters the front [1].
    t_res : float
        Residence time of the parcel in the front [s]. A non-positive or
        non-finite value returns ``porosity_top`` without integrating.
    thickness : float
        Front thickness [m].
    grain_size : float
        Crystal grain size [m].
    delta_rho : float
        Solid minus melt density across the front [kg m-3].
    gravity : float
        Gravitational acceleration [m s-2].
    melt_visc : float
        Melt viscosity [Pa s].
    mush_visc : float
        Shear viscosity of the mush [Pa s].
    buoyancy : str
        The solver's ``separation_buoyancy``, ``'mixture'`` or ``'solid'``.

    Returns
    -------
    tuple
        Retained volume fraction, the percolation time at the entry porosity
        [s], and the matrix deformation time [s].
    """
    from scipy.integrate import solve_ivp

    tau_s = matrix_time(mush_visc, delta_rho, gravity, thickness)
    w_top = darcy_velocity(
        porosity_top, grain_size, delta_rho, gravity, melt_visc, buoyancy=buoyancy
    )
    tau_d_top = thickness / w_top if w_top > 0.0 else float('inf')
    if not tau_d_top > 0.0:
        raise RuntimeError(
            f'percolation time {tau_d_top!r} s is not positive for a front {thickness!r} m thick'
        )

    if not np.isfinite(t_res) or t_res <= 0.0 or porosity_top <= 0.0:
        return float(max(porosity_top, 0.0)), tau_d_top, tau_s

    def rhs(_t, y):
        phi = max(float(y[0]), 0.0)
        w = darcy_velocity(phi, grain_size, delta_rho, gravity, melt_visc, buoyancy=buoyancy)
        tau_d = thickness / w if w > 0.0 else float('inf')
        # The slower process controls: melt leaves only if it can both
        # percolate through the pores and have the matrix deform to close
        # them, so the drainage time is max(tau_D, tau_s), not min.
        rate = max(tau_d, tau_s)
        if not np.isfinite(rate) or rate <= 0.0:
            return [0.0]
        return [-phi / rate]

    sol = solve_ivp(
        rhs,
        (0.0, float(t_res)),
        [float(porosity_top)],
        method='LSODA',
        rtol=_RTOL,
        atol=_ATOL,
    )
    if not sol.success:
        log.warning(
            'Drainage integral did not converge (%s); retaining the entry '
            'porosity as an upper bound.',
            sol.message,
        )
        return float(porosity_top), tau_d_top, tau_s
    retained = float(sol.y[0, -1])
    low, high = 10.0 * _ATOL, 10.0 * _ATOL + _RTOL * porosity_top
    if not -low <= retained <= porosity_top + high:
        raise RuntimeError(
            f'retained porosity {retained!r} lies outside [0, {porosity_top!r}] by more than '
            f'the allowance ({low!r} below, {high!r} above)'
        )
    return float(min(max(retained, 0.0), porosity_top)), tau_d_top, tau_s

def locate_front(
    radius: np.ndarray,
    melt_fraction: np.ndarray,
    porosity: np.ndarray,
    rho_solid: np.ndarray,
    rho_melt: np.ndarray,
    rfront_loc: float,
    phi_min: float,
    n_front_min: int = 3,
    max_front_fraction: float = 1.0,
    r_floor: float | None = None,
) -> tuple[FrontGeometry | None, int]:
    """Locate the freezing front on the interior mesh.

    The front runs from the node where the solver melt fraction crosses
    ``rfront_loc`` at the top down to the node where the porosity falls to
    ``phi_min`` at the base. Using the solver's own melt fraction and its own
    rheological transition for the top means the front trapping uses is the
    front the solver uses for its rheology, rather than a second threshold that
    could disagree with it.

    A mantle crystallising from the bottom up is often still porous at its
    lowest node. The front then rests on the core-mantle boundary, which is
    impermeable, and that boundary is its base: the melt in it drains upward
    only, as it does from a front sitting on compacted cumulate.

    Parameters
    ----------
    r_floor : float or None
        Radius of the impermeable lower boundary of the mantle [m], the base
        of a front that is still porous at the lowest node. Defaults to the
        lowest node radius. Must not lie above it.

    Returns
    -------
    tuple
        The geometry and a branch code. The geometry is ``None`` only when no
        front exists. A geometry paired with :data:`BRANCH_GUARD` is returned
        when the front is present but split in two, too thin to resolve, too
        thick for the thin-front picture, reaching the surface node, or holding
        melt denser than the solid, with the cause in ``guard_reason``: the
        caller should then take the entry porosity as an upper bound rather than
        integrate.
    """
    r = np.asarray(radius, dtype=float)
    r_lo = float(r[0]) if r_floor is None else float(r_floor)
    if r_lo > float(r[0]):
        raise ValueError(
            f'r_floor = {r_lo:.6g} m lies above the lowest node at {float(r[0]):.6g} m; '
            'the lower boundary of the mantle cannot be inside the mesh'
        )
    phi_s = np.asarray(melt_fraction, dtype=float)
    por = np.asarray(porosity, dtype=float)
    mask = (por > phi_min) & (phi_s < rfront_loc)
    runs = _contiguous_runs(mask)
    if not runs:
        return None, BRANCH_NONE

    index = runs[0]
    top, base = int(index[-1]), int(index[0])
    drho_arr = np.asarray(rho_solid, dtype=float) - np.asarray(rho_melt, dtype=float)
    dense = bool(np.nanmin(drho_arr[index]) <= 0.0)

    geom = FrontGeometry(
        index=index,
        r_top=float(r[top]),
        r_base=float(r[base]),
        thickness=float(r[top] - r[base]),
        porosity_top=float(por[top]),
        delta_rho=float(np.nanmean(drho_arr[index])),
        rho_melt=float(np.nanmean(np.asarray(rho_melt, dtype=float)[index])),
        rho_solid_base=float(np.asarray(rho_solid, dtype=float)[base]),
        dense_melt=dense,
    )

    reasons = []
    if len(runs) > 1:
        reasons.append(f'the mush is split into {len(runs)} separate layers')
    if index.size < n_front_min:
        reasons.append(
            f'the front spans {index.size} nodes, fewer than the {n_front_min} '
            'needed to resolve it'
        )
    if top >= r.size - 1:
        reasons.append('the front reaches the surface node')
    if dense:
        reasons.append('the melt is denser than the solid inside the front')
    if reasons:
        geom.guard_reason = '; '.join(reasons)
        return geom, BRANCH_GUARD

    geom.r_top = _interp_crossing(r, phi_s, rfront_loc, top)
    if base == 0:
        geom.r_base = r_lo
    else:
        geom.r_base = _interp_crossing(r, por, phi_min, base - 1)
    geom.thickness = geom.r_top - geom.r_base
    span = float(r[-1] - r_lo)
    if geom.thickness <= 0.0:
        geom.guard_reason = 'the front has no positive thickness'
        return geom, BRANCH_GUARD
    if span > 0.0 and geom.thickness > max_front_fraction * span:
        geom.guard_reason = (
            f'the front spans {geom.thickness / span:.0%} of the mantle, above the '
            f'{max_front_fraction:.0%} a thin front allows'
        )
        return geom, BRANCH_GUARD
    return geom, BRANCH_DARCY

def matrix_time(mush_visc: float, delta_rho: float, gravity: float, thickness: float) -> float:
    """Matrix deformation time [s], McKenzie (1984) compaction viscosity.

    ``tau_s = mu_s / (delta_rho g L)``. The compaction pressure ``delta_rho g L``
    drives a volumetric strain rate against the bulk viscosity ``zeta = mu_s/phi``,
    and closing a porosity ``phi`` takes ``phi zeta / (delta_rho g L)``, which is
    independent of porosity.

    Parameters
    ----------
    mush_visc : float
        Shear viscosity of the mush [Pa s].
    delta_rho : float
        Solid minus melt density [kg m-3]; its magnitude is used.
    gravity : float
        Gravitational acceleration [m s-2].
    thickness : float
        Front thickness [m].

    Returns
    -------
    float
        Deformation time [s]; infinite when ``|delta_rho| g L`` is not positive
        or not finite.
    """
    denom = abs(delta_rho) * gravity * thickness
    if denom <= 0.0 or not np.isfinite(denom):
        return float('inf')
    return float(mush_visc) / denom

def mobility_function(
    porosity: float | np.ndarray, grain_size: float, *, buoyancy: str
) -> float | np.ndarray:
    """Permeability over porosity, ``F = K / phi`` [m2], of the solver's separation law.

    ``aragog.eos.mobility_function``: the three regimes of Bower et al. (2018)
    section 2.1, Blake-Kozeny-Carman ``a^2 phi^2 / ((1-phi)^2 * 1000)``,
    Rumpf-Gupte ``a^2 phi^4.5 * 5/7`` and Stokes ``a^2 * 2/9``, blended with tanh
    weights; under ``buoyancy = 'mixture'`` the first two carry ``1 - phi``.

    Parameters
    ----------
    porosity : float or ndarray
        Volume fraction of melt [1].
    grain_size : float
        Crystal grain size [m].
    buoyancy : str
        The solver's ``separation_buoyancy``, ``'mixture'`` or ``'solid'``.

    Raises
    ------
    ValueError
        If ``grain_size`` is not positive, which would give no permeability, or
        ``buoyancy`` names no mode.
    """
    # Imported here: the solver package takes seconds to load.
    from aragog.eos import mobility_function as solver_mobility

    if grain_size <= 0.0:
        raise ValueError(f'grain_size must be positive, got {grain_size!r}')
    f = solver_mobility(np.asarray(porosity, dtype=float), float(grain_size), buoyancy)
    return float(f) if np.ndim(porosity) == 0 else np.asarray(f, dtype=float)

def porosity_from_densities(
    rho: np.ndarray,
    rho_solid: np.ndarray,
    rho_melt: np.ndarray,
    melt_fraction: np.ndarray,
) -> np.ndarray:
    """Volume fraction of melt from the mixture and end-member densities [1].

    Inside the mush, ``phi = (rho_s - rho) / (rho_s - rho_l)``, the lever rule
    on density. The interior solver builds the mush density from the same two
    phase-boundary densities, so the rule returns the melt volume fraction
    exactly whichever phase is denser: where the tabulated melt is denser than
    the solid, numerator and denominator are both negative and their ratio is
    still the positive melt volume. The solver's own separation flux floors the
    denominator there instead, to stop settling where the melt reads as the
    denser phase; that floor would read such a mush as solid rock and cut the
    front short, so it is not applied here.

    A node that is fully solid or fully molten takes its melt fraction, 0 or 1.
    The solver evaluates its density on the single-phase table at the node's own
    entropy rather than as a two-phase mixture, so the lever rule does not hold
    there: a cold solid denser than the solid at the phase boundary would
    otherwise read as porous wherever the melt is the denser phase. Where the
    two phase-boundary densities are equal the rule is 0/0, and mass and volume
    fractions coincide, so the melt fraction is used there too. The result is
    clipped to [0, 1].

    Parameters
    ----------
    rho : ndarray
        Mixture density [kg m-3].
    rho_solid, rho_melt : ndarray
        Solid and melt densities at the phase boundary [kg m-3].
    melt_fraction : ndarray
        Melt mass fraction of the interior solver [1].
    """
    rho_s, rho_l, rho_mix, phi_s = np.broadcast_arrays(
        np.asarray(rho_solid, dtype=float),
        np.asarray(rho_melt, dtype=float),
        np.asarray(rho, dtype=float),
        np.asarray(melt_fraction, dtype=float),
    )
    drho = rho_s - rho_l
    two_phase = (phi_s > 0.0) & (phi_s < 1.0) & (drho != 0.0)
    with np.errstate(divide='ignore', invalid='ignore'):
        lever = (rho_s - rho_mix) / drho
    return np.clip(np.where(two_phase, lever, phi_s), 0.0, 1.0)

def volume_to_mass_fraction(f_vol: float, rho_melt: float, rho_solid: float) -> float:
    """Convert a trapped melt volume fraction to the mass fraction the budget uses.

    ``F_tl`` is a volume fraction in the compaction physics, where the
    disaggregation limit is a geometric packing threshold, and a mass fraction in
    the volatile budget. The two differ by the melt-to-mixture density ratio,
    about 0.93 at a ten percent density contrast and a porosity of 0.3. The
    conversion is applied once, here, to every fraction the front scheme hands
    to the budget.

    Parameters
    ----------
    f_vol : float
        Melt volume fraction [1]; clipped to [0, 1].
    rho_melt, rho_solid : float
        Melt and solid densities [kg m-3].

    Returns
    -------
    float
        Melt mass fraction [1]. Where the mixture density is not positive or not
        finite, the clipped volume fraction is returned unchanged.
    """
    f = float(np.clip(f_vol, 0.0, 1.0))
    denom = f * rho_melt + (1.0 - f) * rho_solid
    if denom <= 0.0 or not np.isfinite(denom):
        return f
    return float(f * rho_melt / denom)