Alpha Calcualtions

blobid/plic/__init__.py

@@ -3,6 +3,8 @@
 :math:`\vec{n}\cdot\vec{x}=\alpha`
 """
 
-from .plic import get_normals, SUPPORTED_METHODS
+from .plic import get_normals, get_alpha, get_void_fraction, SUPPORTED_METHODS
 
-__all__ = ['get_normals', 'SUPPORTED_METHODS']
+from . import analytic_relations
+
+__all__ = ['get_normals', 'get_alpha', 'get_void_fraction', 'SUPPORTED_METHODS', 'analytic_relations']

blobid/plic/analytic_relations.py

@@ -0,0 +1,137 @@
+r"""
+Implements the analytical relationships described by Scardovelli and Zaleski[1]_.
+
+.. [1] R. Scardovelli & S. Zaleski, "Analytical Relations Connecting Linear Interfaces and Volume Fractions in
+   Rectangular Grids," Journal of Computational Physics, vol. 164, pp. 228-237, 2000.
+   [10.1006/jcph.2000.6567](https://doi.org/10.1006/jcph.2000.6567)
+"""
+import numpy as np
+
+
+def forward_problem(n1, n2, n3, alpha):
+    r"""
+    Solve the three-dimensional forward problem, :math:`f = f(\vec{n}, \alpha)`
+    """
+
+    # scale the slope and sort in ascending order
+    t = abs(n1) + abs(n2) + abs(n3)
+    m1, m2, m3 = sorted([abs(n1)/t, abs(n2)/t, abs(n3)/t])
+
+    # rotate so all normals are positive
+    if n1 < 0:
+        alpha = alpha - n1
+    if n2 < 0:
+        alpha = alpha - n2
+    if n3 < 0:
+        alpha = alpha - n3
+
+    # scale alpha
+    alpha = alpha / t
+
+    # Interface is outside the cell
+    if alpha >= 1:
+        return 1
+    if alpha <= 0:
+        return 0
+
+    # Interface is inside the cell
+    if alpha > 0.5:
+        return 1.0 - _forward_problem_scaled(m1, m2, m3, 1-alpha)
+    return _forward_problem_scaled(m1, m2, m3, alpha)
+
+
+def _forward_problem_scaled(m1, m2, m3, alpha):
+    assert 0 <= alpha <= 0.5
+    assert 0 <= m1 <= m2 <= m3
+
+    m12 = m1 + m2
+    p = 6*m1*m2*m3
+
+    V1 = (1 if m2 == 0 else m1/m2) * (m1/(6*m3))  # this form avoids the need for approximation
+
+    if alpha < m1:
+        return alpha**3 / p
+    if alpha < m2:
+        return alpha*(alpha-m1) / (2*m2*m3) + V1
+    if alpha < min(m12, m3):
+        return (alpha**2 * (3*m12 - alpha) + m1**2 * (m1 - 3*alpha) + m2**2 * (m2 - 3*alpha))/p
+    if m3 < m12:
+        return (alpha**2 * (3 - 2*alpha) + m1**2 * (m1 - 3*alpha) + m2**2 * (m2 - 3*alpha) + m3**2 * (m3 - 3*alpha))/p
+
+    return (2*alpha - m12)/(2*m3)
+
+
+def inverse_problem(n1, n2, n3, f):
+    r"""
+    Solve the three-dimensional inverse problem, :math:`\alpha = \alpha(\vec{n}, f)`
+
+    Only valid for :math:`f\in (0, 1)`
+
+    """
+    assert 0 < f < 1.0
+
+    # scale the slope and sort in ascending order
+    t = abs(n1) + abs(n2) + abs(n3)
+    m1, m2, m3 = sorted([abs(n1)/t, abs(n2)/t, abs(n3)/t])
+
+    # calculate the scaled alpha
+    if f > 0.5:
+        alpha = 1 - _inverse_problem_scaled(m1, m2, m3, 1 - f)
+    else:
+        alpha = _inverse_problem_scaled(m1, m2, m3, f)
+
+    # unscale alpha
+    alpha = alpha * t
+
+    # rotate to original orientation
+    if n1 < 0:
+        alpha = alpha + n1
+    if n2 < 0:
+        alpha = alpha + n2
+    if n3 < 0:
+        alpha = alpha + n3
+
+    return alpha
+
+
+def _inverse_problem_scaled(m1, m2, m3, V):
+    r"""
+    Solve the three-dimensional inverse problem
+    """
+
+    def cubic_root(c3, c2, c1, c0):
+        a0 = c0/c3
+        a1 = c1/c3
+        a2 = c2/c3
+
+        p = a1/3 - a2**2/9
+        q = (a1*a2 - 3*a0)/6 - a2**3/27
+
+        theta = np.arccos(q/np.sqrt(-p**3))/3
+
+        return np.sqrt(-p) * (np.sqrt(3)*np.sin(theta) - np.cos(theta)) - a2/3
+
+    assert 0 <= V <= 0.5
+    assert 0 <= m1 <= m2 <= m3
+
+    m12 = m1 + m2
+    p = 6*m1*m2*m3
+
+    V1 = (1 if m2 == 0 else m1/m2) * (m1/(6*m3))  # this form avoids the need for approximation
+    V2 = V1 + (m2 - m1)/(2 * m3)
+    V3 = (m3**2 * (3 * m12-m3) + m1**2 * (m1 - 3 * m3) + m2**2 * (m2 - 3 * m3)) / p if m3 < m12 else 0.5*m12/m3
+
+    if V < V1:
+        return (p*V)**(1/3)
+    if V < V2:
+        return 0.5 * (m1+np.sqrt(m1**2 + 8 * m2 * m3 * (V-V1)))
+    if V < V3:
+        return cubic_root(
+            -1, 3*m12, -3*(m1**2 + m2**2), m1**3 + m2**3 - p*V
+        )
+    if m3 < m12:
+        return cubic_root(
+            -2, 3, -3*(m1**2 + m2**2 + m3**2), m1**3 + m2**3 + m3**3 - p*V
+        )
+
+    return m3*V + m12/2

blobid/plic/plic.py

@@ -4,6 +4,8 @@
 import numpy as np
 import numpy.typing as npt
 
+from . import analytic_relations
+
 from .central_difference import get_normals as get_normals_CD
 from .nivira import get_normals as get_normals_WY
 
@@ -67,3 +69,71 @@ def get_normals(
             return get_normals_WY(void_fraction, normals_type)
         case _:
             raise ValueError(f"normals_method '{normals_method}' is not supported")
+
+
+def get_alpha(
+        void_fraction: np.ndarray,
+        normals: np.ndarray
+) -> np.ndarray:
+    r"""
+    Calculate the the interface intercept :math:`\alpha=\alpha(f, \vec{n})` using
+    [`blobid.plic.analytic_relations`](analytic_relations.html)
+
+    Parameters
+    ----------
+    void_fraction : ndarray[ni, nj, nk]
+        The void fractions :math:`f` for each grid cell.
+    normals : ndarray[3, ni, nj, nk]
+        The interface normals for each cell.
+
+    Returns
+    -------
+        alpha : ndarray[ni, nj, nk]
+            `alpha[i,j,k]` contains :math:`\alpha` for cell `[i,j,k]`.
+            If there is no interface or `normals` is undefined in the cell, :math:`\alpha=\text{NaN}`.
+    """
+
+    with np.nditer([void_fraction, normals[0, :, :, :], normals[1, :, :, :], normals[2, :, :, :], None],
+                   flags=['common_dtype'], op_dtypes=void_fraction.dtype) as itr:
+        for f, n1, n2, n3, a in itr:
+            if f <= 0 or f >= 1 or (n1 == 0 and n2 == 0 and n3 == 0):
+                # alpha is undefined
+                a[...] = np.nan
+            else:
+                a[...] = analytic_relations.inverse_problem(n1, n2, n3, f)
+
+        return itr.operands[-1]
+
+
+def get_void_fraction(
+        alpha: np.ndarray,
+        normals: np.ndarray
+) -> np.ndarray:
+    r"""
+    Calculate the the void fraction :math:`f=f(\alpha, \vec{n})` using
+    [`blobid.plic.analytic_relations`](analytic_relations.html)
+
+    Parameters
+    ----------
+    alpha : ndarray[ni, nj, nk]
+        The interface intercept :math:`\alpha` for each grid cell.
+    normals : ndarray[3, ni, nj, nk]
+        The interface normals for each cell.
+
+    Returns
+    -------
+        void_fraction : ndarray[ni, nj, nk]
+            `void_fraction[i,j,k]` contains :math:`f` for cell `[i,j,k]`.
+            If `normals` is undefined in the cell, :math:`f=\text{NaN}`.
+    """
+
+    with np.nditer([alpha, normals[0, :, :, :], normals[1, :, :, :], normals[2, :, :, :], None],
+                   flags=['common_dtype'], op_dtypes=alpha.dtype) as itr:
+        for a, n1, n2, n3, f in itr:
+            if np.isnan(a) or (n1 == 0 and n2 == 0 and n3 == 0):
+                # void fraction is undefined
+                f[...] = np.nan
+            else:
+                f[...] = analytic_relations.forward_problem(n1, n2, n3, a)
+
+        return itr.operands[-1]

tests/plic/test_analytic_relations.py

@@ -0,0 +1,64 @@
+import pytest
+
+import numpy as np
+import blobid.plic.analytic_relations as analytic_relations
+
+COUNT = 20000
+
+
+@pytest.fixture
+def slopes() -> np.ndarray:
+    """
+    List of slopes where 0 <= m1 <= m2 <= m2 < 1
+    """
+    m = np.random.rand(COUNT, 3)
+
+    m[:, 1] = m[:, 1] * m[:, 2]
+    m[:, 0] = m[:, 0] * m[:, 1]
+    m = m / np.sum(m, 1, keepdims=True)
+
+    return m
+
+
+@pytest.fixture
+def normals() -> np.ndarray:
+    """
+    List of normals
+    """
+    return (np.random.rand(COUNT, 3) * 3) - 1.5
+
+
+def test_round_trip(normals):
+    for n1, n2, n3, f in zip(normals[:, 0], normals[:, 1], normals[:, 2], np.random.rand(COUNT)*0.5):
+        alpha = analytic_relations.inverse_problem(n1, n2, n3, f)
+        f_test = analytic_relations.forward_problem(n1, n2, n3, alpha)
+
+        assert f_test == pytest.approx(f)
+
+
+def test_scaled_round_trip_forward(slopes):
+    for m1, m2, m3, V in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2], np.random.rand(COUNT)*0.5):
+        alpha = analytic_relations._inverse_problem_scaled(m1, m2, m3, V)
+        V_test = analytic_relations._forward_problem_scaled(m1, m2, m3, alpha)
+
+        assert V_test == pytest.approx(V)
+
+
+def test_scaled_round_trip_backward(slopes):
+    for m1, m2, m3, alpha in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2], np.random.rand(COUNT)*0.5):
+        V = analytic_relations._forward_problem_scaled(m1, m2, m3, alpha)
+        alpha_test = analytic_relations._inverse_problem_scaled(m1, m2, m3, V)
+
+        assert alpha_test == pytest.approx(alpha)
+
+
+def test_scaled_empty(slopes):
+    for m1, m2, m3 in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2]):
+        assert analytic_relations._forward_problem_scaled(m1, m2, m3, 0.0) == 0
+        assert analytic_relations._inverse_problem_scaled(m1, m2, m3, 0.0) == 0
+
+
+def test_scaled_half_full(slopes):
+    for m1, m2, m3 in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2]):
+        assert analytic_relations._forward_problem_scaled(m1, m2, m3, 0.5) == pytest.approx(0.5)
+        assert analytic_relations._inverse_problem_scaled(m1, m2, m3, 0.5) == pytest.approx(0.5)