From 8178701a8dc3400fdc6910bccde7001f7400464b Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Tue, 29 Jul 2025 16:39:02 -0400
Subject: [PATCH 01/17] First try at adding `get_alpha`

---
 blobid/plic/__init__.py |  4 +-
 blobid/plic/plic.py     | 81 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 83 insertions(+), 2 deletions(-)

diff --git a/blobid/plic/__init__.py b/blobid/plic/__init__.py
index 75c7461..5ed5b62 100644
--- a/blobid/plic/__init__.py
+++ b/blobid/plic/__init__.py
@@ -3,6 +3,6 @@
 :math:`\vec{n}\cdot\vec{x}=\alpha`
 """
 
-from .plic import get_normals, SUPPORTED_METHODS
+from .plic import get_normals, get_alpha, SUPPORTED_METHODS
 
-__all__ = ['get_normals', 'SUPPORTED_METHODS']
+__all__ = ['get_normals', 'get_alpha', 'SUPPORTED_METHODS']
diff --git a/blobid/plic/plic.py b/blobid/plic/plic.py
index f568e03..4314482 100644
--- a/blobid/plic/plic.py
+++ b/blobid/plic/plic.py
@@ -67,3 +67,84 @@ 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 normals :math:`\alpha=\alpha(f,\vec{n})`
+    """
+
+    def cubic_root(c3, c2, c1, c0):
+        """Solve cubic roots"""
+        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
+
+    def inverse_problem(m1, m2, m3, v):
+        """
+        This is a slightly modified version of the algorithm described by Scardovelli and Zaleski
+
+        Assumptions:
+        - m1 < m2 < m3
+        - m1 + m2 + m3 = 1
+        - 0 <= V < 0.5
+        """
+
+        m12 = m1 + m2
+        p = 6*m1*m2*m3
+
+        # using m1<=m2, we can satisfy the limit V1->0 as m1,m2->0 exactly
+        # This removes the need for the tolerance discussed by Scardovelli and Zaleski
+        V1 = (1 if m2 == 0 else m1/m2) * (m1/(6*m3))
+        V2 = V1 + (m2 - m1)/(2 * m3)
+        if m3 < m12:
+            V3 = (m3**2 * (3 * m12-m3) + m1**2 * (m1 - 3 * m3) + m2**2 * (m2 - 3 * m3)) / p
+        else:
+            V3 = 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
+
+    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, alpha in itr:
+            # cases where alpha is undefined
+            if f == 0 or f == 1 or (n1 == 0 and n2 == 0 and n3 == 0):
+                alpha[...] = np.nan
+                continue
+
+            # trivial case
+            if f == 0.5:
+                alpha[...] = 0.5
+                continue
+
+            # 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])
+
+            # solve inverse problem
+            alpha[...] = inverse_problem(m1, m2, m3, f if f < 0.5 else 1.0 - f)
+
+        return itr.operands[-1]

From 567ecbdc203b5caf9f1d3569b04af9d81422c3bd Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 12:09:43 -0400
Subject: [PATCH 02/17] Move analytic relations to own file

---
 blobid/plic/__init__.py           |  4 +-
 blobid/plic/analytic_relations.py | 85 +++++++++++++++++++++++++++++++
 blobid/plic/plic.py               | 57 ++-------------------
 3 files changed, 91 insertions(+), 55 deletions(-)
 create mode 100644 blobid/plic/analytic_relations.py

diff --git a/blobid/plic/__init__.py b/blobid/plic/__init__.py
index 5ed5b62..664fd05 100644
--- a/blobid/plic/__init__.py
+++ b/blobid/plic/__init__.py
@@ -5,4 +5,6 @@
 
 from .plic import get_normals, get_alpha, SUPPORTED_METHODS
 
-__all__ = ['get_normals', 'get_alpha', 'SUPPORTED_METHODS']
+from . import analytic_relations
+
+__all__ = ['get_normals', 'get_alpha', 'SUPPORTED_METHODS', 'analytic_relations']
diff --git a/blobid/plic/analytic_relations.py b/blobid/plic/analytic_relations.py
new file mode 100644
index 0000000..2a35bf5
--- /dev/null
+++ b/blobid/plic/analytic_relations.py
@@ -0,0 +1,85 @@
+r"""
+Implements the analytical relationships described by Scardovelli and Zaleski[1]_.
+
+These functions assume:
+.. math::
+    m_1 + m_2 + m_3=1
+.. math::
+    0\le m_1 \le m_2 \le m_3
+.. math::
+    \alpha \text{ and } V \in [0, 0.5]
+
+.. [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(m1, m2, m3, alpha):
+    r"""
+    Solve the three-dimensional forward problem
+    """
+    assert 0 <= alpha <= 0.5
+    assert 0 <= m1 <= m2 <= m3
+    assert np.isclose(m1 + m2 + m3, 1.0)
+
+    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(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
+    assert np.isclose(m1 + m2 + m3, 1.0)
+
+    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
diff --git a/blobid/plic/plic.py b/blobid/plic/plic.py
index 4314482..ced4ee8 100644
--- a/blobid/plic/plic.py
+++ b/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
 
@@ -73,59 +75,6 @@ def get_alpha(
         void_fraction: np.ndarray,
         normals: np.ndarray
 ) -> np.ndarray:
-    r"""
-    Calculate the the interface normals :math:`\alpha=\alpha(f,\vec{n})`
-    """
-
-    def cubic_root(c3, c2, c1, c0):
-        """Solve cubic roots"""
-        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
-
-    def inverse_problem(m1, m2, m3, v):
-        """
-        This is a slightly modified version of the algorithm described by Scardovelli and Zaleski
-
-        Assumptions:
-        - m1 < m2 < m3
-        - m1 + m2 + m3 = 1
-        - 0 <= V < 0.5
-        """
-
-        m12 = m1 + m2
-        p = 6*m1*m2*m3
-
-        # using m1<=m2, we can satisfy the limit V1->0 as m1,m2->0 exactly
-        # This removes the need for the tolerance discussed by Scardovelli and Zaleski
-        V1 = (1 if m2 == 0 else m1/m2) * (m1/(6*m3))
-        V2 = V1 + (m2 - m1)/(2 * m3)
-        if m3 < m12:
-            V3 = (m3**2 * (3 * m12-m3) + m1**2 * (m1 - 3 * m3) + m2**2 * (m2 - 3 * m3)) / p
-        else:
-            V3 = 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
 
     with np.nditer([void_fraction, normals[0, :, :, :], normals[1, :, :, :], normals[2, :, :, :], None],
                    flags=['common_dtype'], op_dtypes=void_fraction.dtype) as itr:
@@ -145,6 +94,6 @@ def inverse_problem(m1, m2, m3, v):
             m1, m2, m3 = sorted([abs(n1)/t, abs(n2)/t, abs(n3)/t])
 
             # solve inverse problem
-            alpha[...] = inverse_problem(m1, m2, m3, f if f < 0.5 else 1.0 - f)
+            alpha[...] = analytic_relations.inverse_problem(m1, m2, m3, f if f < 0.5 else 1.0 - f)
 
         return itr.operands[-1]

From 91ebba5845a842023da1f8b51e36f787f7d206b3 Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 12:38:37 -0400
Subject: [PATCH 03/17] Add tests of analytic relations

---
 tests/plic/test_analytic_relations.py | 47 +++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)
 create mode 100644 tests/plic/test_analytic_relations.py

diff --git a/tests/plic/test_analytic_relations.py b/tests/plic/test_analytic_relations.py
new file mode 100644
index 0000000..d328d1b
--- /dev/null
+++ b/tests/plic/test_analytic_relations.py
@@ -0,0 +1,47 @@
+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
+
+
+def test_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(m1, m2, m3, V)
+        V_test = analytic_relations.forward_problem(m1, m2, m3, alpha)
+
+        assert V_test == pytest.approx(V)
+
+
+def test_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(m1, m2, m3, alpha)
+        alpha_test = analytic_relations.inverse_problem(m1, m2, m3, V)
+
+        assert alpha_test == pytest.approx(alpha)
+
+
+def test_empty(slopes):
+    for m1, m2, m3 in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2]):
+        assert analytic_relations.forward_problem(m1, m2, m3, 0.0) == 0
+        assert analytic_relations.inverse_problem(m1, m2, m3, 0.0) == 0
+
+
+def test_half_full(slopes):
+    for m1, m2, m3 in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2]):
+        assert analytic_relations.forward_problem(m1, m2, m3, 0.5) == pytest.approx(0.5)
+        assert analytic_relations.inverse_problem(m1, m2, m3, 0.5) == pytest.approx(0.5)

From c867d4d7f000e3e26d04a609c676855e94252bdb Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 12:44:05 -0400
Subject: [PATCH 04/17] Support alpha and V in (0.5, 1.0]

---
 blobid/plic/analytic_relations.py     | 12 +++++++++---
 tests/plic/test_analytic_relations.py | 10 ++++++++--
 2 files changed, 17 insertions(+), 5 deletions(-)

diff --git a/blobid/plic/analytic_relations.py b/blobid/plic/analytic_relations.py
index 2a35bf5..425551f 100644
--- a/blobid/plic/analytic_relations.py
+++ b/blobid/plic/analytic_relations.py
@@ -7,7 +7,7 @@
 .. math::
     0\le m_1 \le m_2 \le m_3
 .. math::
-    \alpha \text{ and } V \in [0, 0.5]
+    \alpha \text{ and } V \in [0, 1.0]
 
 .. [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.
@@ -20,10 +20,13 @@ def forward_problem(m1, m2, m3, alpha):
     r"""
     Solve the three-dimensional forward problem
     """
-    assert 0 <= alpha <= 0.5
+    assert 0 <= alpha <= 1.0
     assert 0 <= m1 <= m2 <= m3
     assert np.isclose(m1 + m2 + m3, 1.0)
 
+    if alpha > 0.5:
+        return 1.0 - forward_problem(m1, m2, m3, 1-alpha)
+
     m12 = m1 + m2
     p = 6*m1*m2*m3
 
@@ -58,10 +61,13 @@ def cubic_root(c3, c2, c1, c0):
 
         return np.sqrt(-p) * (np.sqrt(3)*np.sin(theta) - np.cos(theta)) - a2/3
 
-    assert 0 <= V <= 0.5
+    assert 0 <= V <= 1.0
     assert 0 <= m1 <= m2 <= m3
     assert np.isclose(m1 + m2 + m3, 1.0)
 
+    if V > 0.5:
+        return 1.0 - inverse_problem(m1, m2, m3, 1-V)
+
     m12 = m1 + m2
     p = 6*m1*m2*m3
 
diff --git a/tests/plic/test_analytic_relations.py b/tests/plic/test_analytic_relations.py
index d328d1b..3408106 100644
--- a/tests/plic/test_analytic_relations.py
+++ b/tests/plic/test_analytic_relations.py
@@ -20,7 +20,7 @@ def slopes() -> np.ndarray:
 
 
 def test_round_trip_forward(slopes):
-    for m1, m2, m3, V in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2], np.random.rand(COUNT) * 0.5):
+    for m1, m2, m3, V in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2], np.random.rand(COUNT)):
         alpha = analytic_relations.inverse_problem(m1, m2, m3, V)
         V_test = analytic_relations.forward_problem(m1, m2, m3, alpha)
 
@@ -28,7 +28,7 @@ def test_round_trip_forward(slopes):
 
 
 def test_round_trip_backward(slopes):
-    for m1, m2, m3, alpha in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2], np.random.rand(COUNT) * 0.5):
+    for m1, m2, m3, alpha in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2], np.random.rand(COUNT)):
         V = analytic_relations.forward_problem(m1, m2, m3, alpha)
         alpha_test = analytic_relations.inverse_problem(m1, m2, m3, V)
 
@@ -45,3 +45,9 @@ def test_half_full(slopes):
     for m1, m2, m3 in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2]):
         assert analytic_relations.forward_problem(m1, m2, m3, 0.5) == pytest.approx(0.5)
         assert analytic_relations.inverse_problem(m1, m2, m3, 0.5) == pytest.approx(0.5)
+
+
+def test_full(slopes):
+    for m1, m2, m3 in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2]):
+        assert analytic_relations.forward_problem(m1, m2, m3, 1.0) == 1
+        assert analytic_relations.inverse_problem(m1, m2, m3, 1.0) == 1

From c67131806c204b91ceaf261e9e7ca8c73d5a8148 Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 13:22:40 -0400
Subject: [PATCH 05/17] Take care of scaling in analytic_relations

---
 blobid/plic/analytic_relations.py     | 71 ++++++++++++++++++++++-----
 blobid/plic/plic.py                   | 19 ++-----
 tests/plic/test_analytic_relations.py | 46 +++++++++++------
 3 files changed, 94 insertions(+), 42 deletions(-)

diff --git a/blobid/plic/analytic_relations.py b/blobid/plic/analytic_relations.py
index 425551f..08c2459 100644
--- a/blobid/plic/analytic_relations.py
+++ b/blobid/plic/analytic_relations.py
@@ -1,14 +1,6 @@
 r"""
 Implements the analytical relationships described by Scardovelli and Zaleski[1]_.
 
-These functions assume:
-.. math::
-    m_1 + m_2 + m_3=1
-.. math::
-    0\le m_1 \le m_2 \le m_3
-.. math::
-    \alpha \text{ and } V \in [0, 1.0]
-
 .. [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)
@@ -16,10 +8,36 @@
 import numpy as np
 
 
-def forward_problem(m1, m2, m3, alpha):
+def forward_problem(n1, n2, n3, alpha):
     r"""
-    Solve the three-dimensional forward problem
+    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
+
+    return _forward_problem_scaled(m1, m2, m3, alpha)
+
+
+def _forward_problem_scaled(m1, m2, m3, alpha):
     assert 0 <= alpha <= 1.0
     assert 0 <= m1 <= m2 <= m3
     assert np.isclose(m1 + m2 + m3, 1.0)
@@ -44,7 +62,36 @@ def forward_problem(m1, m2, m3, alpha):
     return (2*alpha - m12)/(2*m3)
 
 
-def inverse_problem(m1, m2, m3, V):
+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])
+
+    alpha = _inverse_problem_scaled(m1, m2, m3, f)
+
+    # scale 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
     """
@@ -66,7 +113,7 @@ def cubic_root(c3, c2, c1, c0):
     assert np.isclose(m1 + m2 + m3, 1.0)
 
     if V > 0.5:
-        return 1.0 - inverse_problem(m1, m2, m3, 1-V)
+        return 1.0 - _inverse_problem_scaled(m1, m2, m3, 1-V)
 
     m12 = m1 + m2
     p = 6*m1*m2*m3
diff --git a/blobid/plic/plic.py b/blobid/plic/plic.py
index ced4ee8..01648c5 100644
--- a/blobid/plic/plic.py
+++ b/blobid/plic/plic.py
@@ -79,21 +79,10 @@ def get_alpha(
     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, alpha in itr:
-            # cases where alpha is undefined
-            if f == 0 or f == 1 or (n1 == 0 and n2 == 0 and n3 == 0):
+            if f <= 0 or f >= 1 or (n1 == 0 and n2 == 0 and n3 == 0):
+                # alpha is undefined
                 alpha[...] = np.nan
-                continue
-
-            # trivial case
-            if f == 0.5:
-                alpha[...] = 0.5
-                continue
-
-            # 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])
-
-            # solve inverse problem
-            alpha[...] = analytic_relations.inverse_problem(m1, m2, m3, f if f < 0.5 else 1.0 - f)
+            else:
+                alpha[...] = analytic_relations.inverse_problem(n1, n2, n3, f)     
 
         return itr.operands[-1]
diff --git a/tests/plic/test_analytic_relations.py b/tests/plic/test_analytic_relations.py
index 3408106..86dfc23 100644
--- a/tests/plic/test_analytic_relations.py
+++ b/tests/plic/test_analytic_relations.py
@@ -19,35 +19,51 @@ def slopes() -> np.ndarray:
     return m
 
 
-def test_round_trip_forward(slopes):
+@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)):
+        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)):
-        alpha = analytic_relations.inverse_problem(m1, m2, m3, V)
-        V_test = analytic_relations.forward_problem(m1, m2, m3, alpha)
+        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_round_trip_backward(slopes):
+def test_scaled_round_trip_backward(slopes):
     for m1, m2, m3, alpha in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2], np.random.rand(COUNT)):
-        V = analytic_relations.forward_problem(m1, m2, m3, alpha)
-        alpha_test = analytic_relations.inverse_problem(m1, m2, m3, V)
+        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_empty(slopes):
+def test_scaled_empty(slopes):
     for m1, m2, m3 in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2]):
-        assert analytic_relations.forward_problem(m1, m2, m3, 0.0) == 0
-        assert analytic_relations.inverse_problem(m1, m2, m3, 0.0) == 0
+        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_half_full(slopes):
+def test_scaled_half_full(slopes):
     for m1, m2, m3 in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2]):
-        assert analytic_relations.forward_problem(m1, m2, m3, 0.5) == pytest.approx(0.5)
-        assert analytic_relations.inverse_problem(m1, m2, m3, 0.5) == pytest.approx(0.5)
+        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)
 
 
-def test_full(slopes):
+def test_scaled_full(slopes):
     for m1, m2, m3 in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2]):
-        assert analytic_relations.forward_problem(m1, m2, m3, 1.0) == 1
-        assert analytic_relations.inverse_problem(m1, m2, m3, 1.0) == 1
+        assert analytic_relations._forward_problem_scaled(m1, m2, m3, 1.0) == 1
+        assert analytic_relations._inverse_problem_scaled(m1, m2, m3, 1.0) == 1

From 2d69c3ccb991af60657e82a8b42b78b9352a2b0a Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 13:28:16 -0400
Subject: [PATCH 06/17] Avoid recursion

---
 blobid/plic/analytic_relations.py     | 21 +++++++++++----------
 tests/plic/test_analytic_relations.py | 12 +++---------
 2 files changed, 14 insertions(+), 19 deletions(-)

diff --git a/blobid/plic/analytic_relations.py b/blobid/plic/analytic_relations.py
index 08c2459..44b01cf 100644
--- a/blobid/plic/analytic_relations.py
+++ b/blobid/plic/analytic_relations.py
@@ -34,17 +34,17 @@ def forward_problem(n1, n2, n3, alpha):
     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 <= 1.0
+    assert 0 <= alpha <= 0.5
     assert 0 <= m1 <= m2 <= m3
     assert np.isclose(m1 + m2 + m3, 1.0)
 
-    if alpha > 0.5:
-        return 1.0 - forward_problem(m1, m2, m3, 1-alpha)
-
     m12 = m1 + m2
     p = 6*m1*m2*m3
 
@@ -75,9 +75,13 @@ def inverse_problem(n1, n2, n3, f):
     t = abs(n1) + abs(n2) + abs(n3)
     m1, m2, m3 = sorted([abs(n1)/t, abs(n2)/t, abs(n3)/t])
 
-    alpha = _inverse_problem_scaled(m1, m2, m3, f)
+    # 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)
 
-    # scale alpha
+    # unscale alpha
     alpha = alpha * t
 
     # rotate to original orientation
@@ -108,13 +112,10 @@ def cubic_root(c3, c2, c1, c0):
 
         return np.sqrt(-p) * (np.sqrt(3)*np.sin(theta) - np.cos(theta)) - a2/3
 
-    assert 0 <= V <= 1.0
+    assert 0 <= V <= 0.5
     assert 0 <= m1 <= m2 <= m3
     assert np.isclose(m1 + m2 + m3, 1.0)
 
-    if V > 0.5:
-        return 1.0 - _inverse_problem_scaled(m1, m2, m3, 1-V)
-
     m12 = m1 + m2
     p = 6*m1*m2*m3
 
diff --git a/tests/plic/test_analytic_relations.py b/tests/plic/test_analytic_relations.py
index 86dfc23..9364a52 100644
--- a/tests/plic/test_analytic_relations.py
+++ b/tests/plic/test_analytic_relations.py
@@ -28,7 +28,7 @@ def normals() -> np.ndarray:
 
 
 def test_round_trip(normals):
-    for n1, n2, n3, f in zip(normals[:, 0], normals[:, 1], normals[:, 2], np.random.rand(COUNT)):
+    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)
 
@@ -36,7 +36,7 @@ def test_round_trip(normals):
 
 
 def test_scaled_round_trip_forward(slopes):
-    for m1, m2, m3, V in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2], np.random.rand(COUNT)):
+    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)
 
@@ -44,7 +44,7 @@ def test_scaled_round_trip_forward(slopes):
 
 
 def test_scaled_round_trip_backward(slopes):
-    for m1, m2, m3, alpha in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2], np.random.rand(COUNT)):
+    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)
 
@@ -61,9 +61,3 @@ 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)
-
-
-def test_scaled_full(slopes):
-    for m1, m2, m3 in zip(slopes[:, 0], slopes[:, 1], slopes[:, 2]):
-        assert analytic_relations._forward_problem_scaled(m1, m2, m3, 1.0) == 1
-        assert analytic_relations._inverse_problem_scaled(m1, m2, m3, 1.0) == 1

From c4490a2fdf671b3cf4ae9b8e1030de183ba60ca6 Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 14:07:41 -0400
Subject: [PATCH 07/17] Added get_void_fraction function and round-trip test

---
 blobid/plic/__init__.py                       |  4 ++--
 blobid/plic/analytic_relations.py             |  2 +-
 blobid/plic/plic.py                           | 23 ++++++++++++++++---
 .../{test_get_normals.py => test_plic.py}     | 19 +++++++++++++++
 4 files changed, 42 insertions(+), 6 deletions(-)
 rename tests/plic/{test_get_normals.py => test_plic.py} (64%)

diff --git a/blobid/plic/__init__.py b/blobid/plic/__init__.py
index 664fd05..5dc742f 100644
--- a/blobid/plic/__init__.py
+++ b/blobid/plic/__init__.py
@@ -3,8 +3,8 @@
 :math:`\vec{n}\cdot\vec{x}=\alpha`
 """
 
-from .plic import get_normals, get_alpha, SUPPORTED_METHODS
+from .plic import get_normals, get_alpha, get_void_fraction, SUPPORTED_METHODS
 
 from . import analytic_relations
 
-__all__ = ['get_normals', 'get_alpha', 'SUPPORTED_METHODS', 'analytic_relations']
+__all__ = ['get_normals', 'get_alpha', 'get_void_fraction', 'SUPPORTED_METHODS', 'analytic_relations']
diff --git a/blobid/plic/analytic_relations.py b/blobid/plic/analytic_relations.py
index 44b01cf..0e536de 100644
--- a/blobid/plic/analytic_relations.py
+++ b/blobid/plic/analytic_relations.py
@@ -69,7 +69,7 @@ def inverse_problem(n1, n2, n3, f):
     Only valid for :math:`f\in (0, 1)`
 
     """
-    assert 0 <= f <= 1.0
+    assert 0 < f < 1.0
 
     # scale the slope and sort in ascending order
     t = abs(n1) + abs(n2) + abs(n3)
diff --git a/blobid/plic/plic.py b/blobid/plic/plic.py
index 01648c5..e7be7a7 100644
--- a/blobid/plic/plic.py
+++ b/blobid/plic/plic.py
@@ -78,11 +78,28 @@ def get_alpha(
 
     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, alpha in 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
-                alpha[...] = np.nan
+                a[...] = np.nan
             else:
-                alpha[...] = analytic_relations.inverse_problem(n1, n2, n3, f)     
+                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:
+
+    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]
diff --git a/tests/plic/test_get_normals.py b/tests/plic/test_plic.py
similarity index 64%
rename from tests/plic/test_get_normals.py
rename to tests/plic/test_plic.py
index b5829f1..b7a508a 100644
--- a/tests/plic/test_get_normals.py
+++ b/tests/plic/test_plic.py
@@ -1,5 +1,7 @@
 import hashlib
 
+import pytest
+
 import numpy as np
 import blobid.plic as plic
 
@@ -41,3 +43,20 @@ def test_end_to_end(fs_vof):
 
         hash_object = hashlib.md5(normals.tobytes(order='C'))
         assert hash_object.hexdigest() == signature
+
+
+def test_round_trip(fs_vof_medium):
+    normals = plic.get_normals(fs_vof_medium.astype(np.float32), 'CD')
+    undefined_normal = np.all(normals == 0, axis=0)
+
+    void_fraction = fs_vof_medium[1:-1, 1:-1, 1:-1].astype(np.float32)
+
+    alpha = plic.get_alpha(void_fraction, normals)
+    void_fraction_test = plic.get_void_fraction(alpha, normals)
+
+    with np.nditer([void_fraction, void_fraction_test, undefined_normal]) as itr:
+        for f, f_test, undef in itr:
+            if np.isnan(f_test):
+                assert undef or f == 0 or f == 1
+            else:
+                assert f_test == pytest.approx(f, abs=1e-6)

From 1d9ff5302ff57126ef16ad2aa4ba67b75dc0407c Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 14:19:46 -0400
Subject: [PATCH 08/17] Add documentation

---
 blobid/plic/plic.py | 34 ++++++++++++++++++++++++++++++++++
 1 file changed, 34 insertions(+)

diff --git a/blobid/plic/plic.py b/blobid/plic/plic.py
index e7be7a7..bbc7c96 100644
--- a/blobid/plic/plic.py
+++ b/blobid/plic/plic.py
@@ -75,6 +75,23 @@ 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:
@@ -92,6 +109,23 @@ 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:

From 003085f5ed55c0954a16559bebf7d4d4debb1293 Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 14:25:57 -0400
Subject: [PATCH 09/17] Add alpha to end_to_end hash test

---
 tests/plic/test_plic.py | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/tests/plic/test_plic.py b/tests/plic/test_plic.py
index b7a508a..8b4ce11 100644
--- a/tests/plic/test_plic.py
+++ b/tests/plic/test_plic.py
@@ -34,15 +34,20 @@ def assert_undefined(norm):
 
 def test_end_to_end(fs_vof):
     expected_result = [
-        ['CD', '599820a35d442746d95e8cd2f06fbc61'],
-        ['WY', '2aeb743728db4ecefda7ce7527d836b4']
+        ['CD', '599820a35d442746d95e8cd2f06fbc61', '969c062702d97ebc363ba20dbde987c0'],
+        ['WY', '2aeb743728db4ecefda7ce7527d836b4', '97a7522dc051332d0bc0d1020b9a6f85']
     ]
 
-    for method, signature in expected_result:
+    for method, normals_signature, alpha_signature in expected_result:
         normals = plic.get_normals(fs_vof.astype(np.float32), method)
 
         hash_object = hashlib.md5(normals.tobytes(order='C'))
-        assert hash_object.hexdigest() == signature
+        assert hash_object.hexdigest() == normals_signature
+
+        alpha = plic.get_alpha(fs_vof[1:-1, 1:-1, 1:-1].astype(np.float32), normals)
+
+        hash_object = hashlib.md5(alpha.tobytes(order='C'))
+        assert hash_object.hexdigest() == alpha_signature
 
 
 def test_round_trip(fs_vof_medium):

From a2a8e38ad92ce9f7efb7824ec46aae47b1a7e75f Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 14:34:11 -0400
Subject: [PATCH 10/17] Remove expensive assert

---
 blobid/plic/analytic_relations.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/blobid/plic/analytic_relations.py b/blobid/plic/analytic_relations.py
index 0e536de..fa924f8 100644
--- a/blobid/plic/analytic_relations.py
+++ b/blobid/plic/analytic_relations.py
@@ -43,7 +43,6 @@ def forward_problem(n1, n2, n3, alpha):
 def _forward_problem_scaled(m1, m2, m3, alpha):
     assert 0 <= alpha <= 0.5
     assert 0 <= m1 <= m2 <= m3
-    assert np.isclose(m1 + m2 + m3, 1.0)
 
     m12 = m1 + m2
     p = 6*m1*m2*m3
@@ -114,7 +113,6 @@ def cubic_root(c3, c2, c1, c0):
 
     assert 0 <= V <= 0.5
     assert 0 <= m1 <= m2 <= m3
-    assert np.isclose(m1 + m2 + m3, 1.0)
 
     m12 = m1 + m2
     p = 6*m1*m2*m3

From c1a2b1d6a7caa26c3dbb550b38ef7c98ae770306 Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 14:39:16 -0400
Subject: [PATCH 11/17] Add benchmarks for plic

---
 tests/plic/test_plic.py | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/tests/plic/test_plic.py b/tests/plic/test_plic.py
index 8b4ce11..54c1e98 100644
--- a/tests/plic/test_plic.py
+++ b/tests/plic/test_plic.py
@@ -65,3 +65,22 @@ def test_round_trip(fs_vof_medium):
                 assert undef or f == 0 or f == 1
             else:
                 assert f_test == pytest.approx(f, abs=1e-6)
+
+
+def test_bench_get_alpha(benchmark, fs_vof_small):
+    normals = plic.get_normals(fs_vof_small.astype(np.float32), 'WY')
+    void_fraction = fs_vof_small[1:-1, 1:-1, 1:-1].astype(np.float32)
+
+    benchmark(plic.get_alpha,
+              void_fraction=void_fraction,
+              normals=normals)
+
+
+def test_bench_get_void_fraction(benchmark, fs_vof_small):
+    normals = plic.get_normals(fs_vof_small.astype(np.float32), 'WY')
+    void_fraction = fs_vof_small[1:-1, 1:-1, 1:-1].astype(np.float32)
+    alpha = plic.get_alpha(void_fraction, normals)
+
+    benchmark(plic.get_void_fraction,
+              alpha=alpha,
+              normals=normals)

From e52c428ae52cb492d11283999da4bc7fe35d984a Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 14:58:11 -0400
Subject: [PATCH 12/17] Use njit

---
 blobid/plic/analytic_relations.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/blobid/plic/analytic_relations.py b/blobid/plic/analytic_relations.py
index fa924f8..aa89f31 100644
--- a/blobid/plic/analytic_relations.py
+++ b/blobid/plic/analytic_relations.py
@@ -6,6 +6,7 @@
    [10.1006/jcph.2000.6567](https://doi.org/10.1006/jcph.2000.6567)
 """
 import numpy as np
+from .._numba_support import njit
 
 
 def forward_problem(n1, n2, n3, alpha):
@@ -40,6 +41,7 @@ def forward_problem(n1, n2, n3, alpha):
     return _forward_problem_scaled(m1, m2, m3, alpha)
 
 
+@njit
 def _forward_problem_scaled(m1, m2, m3, alpha):
     assert 0 <= alpha <= 0.5
     assert 0 <= m1 <= m2 <= m3
@@ -94,6 +96,7 @@ def inverse_problem(n1, n2, n3, f):
     return alpha
 
 
+@njit
 def _inverse_problem_scaled(m1, m2, m3, V):
     r"""
     Solve the three-dimensional inverse problem

From 3290a74140f2e37c5618b350fc82648ce8d64275 Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 15:47:02 -0400
Subject: [PATCH 13/17] Format cleanup

---
 tests/plic/test_analytic_relations.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/plic/test_analytic_relations.py b/tests/plic/test_analytic_relations.py
index 9364a52..49b2e32 100644
--- a/tests/plic/test_analytic_relations.py
+++ b/tests/plic/test_analytic_relations.py
@@ -5,6 +5,7 @@
 
 COUNT = 20000
 
+
 @pytest.fixture
 def slopes() -> np.ndarray:
     """

From a530d40fd9b9c39c3fce57f434c703e7eb6864bf Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 15:56:17 -0400
Subject: [PATCH 14/17] Use double precision for test_round_trip

---
 tests/plic/test_plic.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/plic/test_plic.py b/tests/plic/test_plic.py
index 54c1e98..592cb41 100644
--- a/tests/plic/test_plic.py
+++ b/tests/plic/test_plic.py
@@ -51,10 +51,10 @@ def test_end_to_end(fs_vof):
 
 
 def test_round_trip(fs_vof_medium):
-    normals = plic.get_normals(fs_vof_medium.astype(np.float32), 'CD')
+    normals = plic.get_normals(fs_vof_medium.astype(np.float64), 'CD')
     undefined_normal = np.all(normals == 0, axis=0)
 
-    void_fraction = fs_vof_medium[1:-1, 1:-1, 1:-1].astype(np.float32)
+    void_fraction = fs_vof_medium[1:-1, 1:-1, 1:-1].astype(np.float64)
 
     alpha = plic.get_alpha(void_fraction, normals)
     void_fraction_test = plic.get_void_fraction(alpha, normals)
@@ -64,7 +64,7 @@ def test_round_trip(fs_vof_medium):
             if np.isnan(f_test):
                 assert undef or f == 0 or f == 1
             else:
-                assert f_test == pytest.approx(f, abs=1e-6)
+                assert f_test == pytest.approx(f)
 
 
 def test_bench_get_alpha(benchmark, fs_vof_small):

From e4de765564f4d38059719922dd210b8079db60a4 Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 15:58:15 -0400
Subject: [PATCH 15/17] Revert "Use njit"

This reverts commit e52c428ae52cb492d11283999da4bc7fe35d984a.
---
 blobid/plic/analytic_relations.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/blobid/plic/analytic_relations.py b/blobid/plic/analytic_relations.py
index aa89f31..fa924f8 100644
--- a/blobid/plic/analytic_relations.py
+++ b/blobid/plic/analytic_relations.py
@@ -6,7 +6,6 @@
    [10.1006/jcph.2000.6567](https://doi.org/10.1006/jcph.2000.6567)
 """
 import numpy as np
-from .._numba_support import njit
 
 
 def forward_problem(n1, n2, n3, alpha):
@@ -41,7 +40,6 @@ def forward_problem(n1, n2, n3, alpha):
     return _forward_problem_scaled(m1, m2, m3, alpha)
 
 
-@njit
 def _forward_problem_scaled(m1, m2, m3, alpha):
     assert 0 <= alpha <= 0.5
     assert 0 <= m1 <= m2 <= m3
@@ -96,7 +94,6 @@ def inverse_problem(n1, n2, n3, f):
     return alpha
 
 
-@njit
 def _inverse_problem_scaled(m1, m2, m3, V):
     r"""
     Solve the three-dimensional inverse problem

From 63638e63baf86e274110b234cb1d3757c0b74bf7 Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 16:02:40 -0400
Subject: [PATCH 16/17] Allow some differences in alpha

---
 tests/plic/test_plic.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/tests/plic/test_plic.py b/tests/plic/test_plic.py
index 592cb41..010100e 100644
--- a/tests/plic/test_plic.py
+++ b/tests/plic/test_plic.py
@@ -34,8 +34,8 @@ def assert_undefined(norm):
 
 def test_end_to_end(fs_vof):
     expected_result = [
-        ['CD', '599820a35d442746d95e8cd2f06fbc61', '969c062702d97ebc363ba20dbde987c0'],
-        ['WY', '2aeb743728db4ecefda7ce7527d836b4', '97a7522dc051332d0bc0d1020b9a6f85']
+        ['CD', '599820a35d442746d95e8cd2f06fbc61', '04e255859bc9db0720dbacca5556b300'],
+        ['WY', '2aeb743728db4ecefda7ce7527d836b4', '0bbcb27a94c4758e0f05f88d526ef4fc']
     ]
 
     for method, normals_signature, alpha_signature in expected_result:
@@ -46,6 +46,9 @@ def test_end_to_end(fs_vof):
 
         alpha = plic.get_alpha(fs_vof[1:-1, 1:-1, 1:-1].astype(np.float32), normals)
 
+        # cast alpha to lower precision to allow some floating point differences
+        alpha = alpha.astype(np.float16)
+
         hash_object = hashlib.md5(alpha.tobytes(order='C'))
         assert hash_object.hexdigest() == alpha_signature
 

From 465ed9fe59afe44c3529a506e4a2bf492d6181fd Mon Sep 17 00:00:00 2001
From: Declan Gaylo <dgaylo@mit.edu>
Date: Mon, 11 Aug 2025 16:08:42 -0400
Subject: [PATCH 17/17] Remove hash test of alpha

---
 tests/plic/test_plic.py | 14 +++-----------
 1 file changed, 3 insertions(+), 11 deletions(-)

diff --git a/tests/plic/test_plic.py b/tests/plic/test_plic.py
index 010100e..3c1a6ee 100644
--- a/tests/plic/test_plic.py
+++ b/tests/plic/test_plic.py
@@ -34,24 +34,16 @@ def assert_undefined(norm):
 
 def test_end_to_end(fs_vof):
     expected_result = [
-        ['CD', '599820a35d442746d95e8cd2f06fbc61', '04e255859bc9db0720dbacca5556b300'],
-        ['WY', '2aeb743728db4ecefda7ce7527d836b4', '0bbcb27a94c4758e0f05f88d526ef4fc']
+        ['CD', '599820a35d442746d95e8cd2f06fbc61'],
+        ['WY', '2aeb743728db4ecefda7ce7527d836b4']
     ]
 
-    for method, normals_signature, alpha_signature in expected_result:
+    for method, normals_signature in expected_result:
         normals = plic.get_normals(fs_vof.astype(np.float32), method)
 
         hash_object = hashlib.md5(normals.tobytes(order='C'))
         assert hash_object.hexdigest() == normals_signature
 
-        alpha = plic.get_alpha(fs_vof[1:-1, 1:-1, 1:-1].astype(np.float32), normals)
-
-        # cast alpha to lower precision to allow some floating point differences
-        alpha = alpha.astype(np.float16)
-
-        hash_object = hashlib.md5(alpha.tobytes(order='C'))
-        assert hash_object.hexdigest() == alpha_signature
-
 
 def test_round_trip(fs_vof_medium):
     normals = plic.get_normals(fs_vof_medium.astype(np.float64), 'CD')