Kinetic energy fraction (#329)

* Coded a function to compute the efficiency factor in Hydrodynamics

* Coded it for the template model as well

* Made a test for the efficiency factor

* Replaced simps by simpson

* Try something to make the tests work

* added a test for the efficiency factor in the 2-step

* Added the formula in the docstring

---------

Co-authored-by: Jorinde van de Vis <jorindevandevis@outlook.com>

Author: benoitlaurent96

src/WallGo/hydrodynamics.py

@@ -7,7 +7,7 @@
 import numpy as np
 import numpy.typing as npt
 from scipy.optimize import root_scalar, root, minimize_scalar
-from scipy.integrate import solve_ivp
+from scipy.integrate import solve_ivp, simpson
 from .exceptions import WallGoError
 from .thermodynamics import Thermodynamics
 from .hydrodynamicsTemplateModel import HydrodynamicsTemplateModel
@@ -484,7 +484,7 @@ def matching(
         return vp, vm, Tp, Tm
 
     def shockDE(
-        self, v: float, xiAndT: np.ndarray
+        self, v: float, xiAndT: np.ndarray, shockWave: bool=True
     ) -> Tuple[npt.ArrayLike, npt.ArrayLike]:
         r"""
         Hydrodynamic equations for the self-similar coordinate :math:`\xi = r/t` and
@@ -498,6 +498,9 @@ def shockDE(
         xiAndT : array
             Values of the self-similar coordinate :math:`\xi = r/t` and
             the temperature :math:`T`
+        shockWave : bool, optional
+            If True, the integration happens in the shock wave. If False, it happens in
+            the rarefaction wave. Default is True.
 
         Returns
         -------
@@ -515,20 +518,19 @@ def shockDE(
                 {"v": v, "xi": xi, "T": T},
             )
 
+        if shockWave:
+            csq = self.thermodynamics.csqHighT(T)
+        else:
+            csq = self.thermodynamics.csqLowT(T)
         eq1 = (
             gammaSq(v)
             * (1.0 - v * xi)
-            * (boostVelocity(xi, v) ** 2 / self.thermodynamics.csqHighT(T) - 1.0)
+            * (boostVelocity(xi, v) ** 2 / csq - 1.0)
             * xi
             / 2.0
             / v
         )
-        eq2 = (
-            self.thermodynamics.wHighT(T)
-            / self.thermodynamics.dpHighT(T)
-            * gammaSq(v)
-            * boostVelocity(xi, v)
-        )
+        eq2 = T * gammaSq(v) * boostVelocity(xi, v)
         return [eq1, eq2]
 
     def solveHydroShock(self, vw: float, vp: float, Tp: float) -> float:
@@ -558,7 +560,7 @@ def shock(v: float, xiAndT: np.ndarray | list) -> float:
             xi, T = xiAndT
             return float(boostVelocity(xi, v) * xi - self.thermodynamics.csqHighT(T))
 
-        shock.terminal = True  # What's happening here?
+        shock.terminal = True
         xi0T0 = [vw, Tp]
         vpcent = boostVelocity(vw, vp)
         if shock(vpcent, xi0T0) > 0:
@@ -903,6 +905,85 @@ def shock(
         )
         return float(sol.root)
 
+    def efficiencyFactor(self, vw: float) -> float:
+        r"""
+        Computes the efficiency factor
+        :math:`\kappa=\frac{4}{v_w^3 \alpha_n w_n}\int d\xi \xi^2 v^2\gamma^2 w`.
+
+        Parameters
+        ----------
+        vw : float
+            Wall velocity.
+
+        Returns
+        -------
+        float
+            Value of the efficiency factor :math:`\kappa`.
+
+        """
+        # Separates the efficiency factor into a contribution from the shock wave and
+        # the rarefaction wave.
+        kappaSW = 0.0
+        kappaRW = 0.0
+
+        vp, vm, Tp, Tm = self.findMatching(vw)
+
+        # If deflagration or hybrid, computes the shock wave contribution
+        if vw < self.vJ:
+            def shock(v: float, xiAndT: np.ndarray | list) -> float:
+                xi, T = xiAndT
+                return float(boostVelocity(xi, v)*xi - self.thermodynamics.csqHighT(T))
+
+            shock.terminal = True
+            xi0T0 = [vw, Tp]
+            vpcent = boostVelocity(vw, vp)
+            if shock(vpcent, xi0T0) < 0 and vw != vp:
+                # Integrate the shock wave
+                solShock = solve_ivp(
+                    self.shockDE,
+                    [vpcent, 1e-10],
+                    xi0T0,
+                    events=shock,
+                    rtol=self.rtol,
+                    atol=0,
+                )  # solve differential equation all the way from v = v+ to v = 0
+                vPlasma = solShock.t
+                xi = solShock.y[0]
+                T = solShock.y[1]
+                enthalpy = np.array([self.thermodynamics.wHighT(t) for t in T])
+
+                # Integrate the solution to get kappa
+                kappaSW = 4 * simpson(
+                    y=xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy,
+                    x=xi
+                ) / (vw**3*self.thermodynamics.wHighT(self.Tnucl)*self.template.alN)
+
+        # If hybrid or detonation, computes the rarefaction wave contribution
+        if vw**2 > self.thermodynamics.csqLowT(Tm):
+            xi0T0 = [vw, Tm]
+            vmcent = boostVelocity(vw, vm)
+            # Integrate the rarefaction wave
+            solRarefaction = solve_ivp(
+                self.shockDE,
+                [vmcent, 1e-10],
+                xi0T0,
+                rtol=self.rtol,
+                atol=0,
+                args=(False,)
+            )  # solve differential equation all the way from v = v- to v = 0
+            vPlasma = solRarefaction.t
+            xi = solRarefaction.y[0]
+            T = solRarefaction.y[1]
+            enthalpy = np.array([self.thermodynamics.wLowT(t) for t in T])
+
+            # Integrate the solution to get kappa
+            kappaRW = -4 * simpson(
+                y=xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy,
+                x=xi
+            ) / (vw**3*self.thermodynamics.wHighT(self.Tnucl)*self.template.alN)
+
+        return kappaSW + kappaRW
+
     def _mappingT(self, TpTm: list[float]) -> list[float]:
         """
         Maps the variables Tp and Tm, which are constrained to

src/WallGo/hydrodynamicsTemplateModel.py

@@ -6,10 +6,10 @@
 import warnings
 import logging
 import numpy as np
-from scipy.integrate import solve_ivp
+from scipy.integrate import solve_ivp, simpson
 from scipy.optimize import root_scalar, minimize_scalar, OptimizeResult
 from .exceptions import WallGoError
-from .helpers import gammaSq
+from .helpers import gammaSq, boostVelocity
 from .thermodynamics import Thermodynamics
 
 
@@ -328,23 +328,31 @@ def solveAlpha(self, vw: float, constraint: bool = True) -> float:
         except ValueError as exc:
             raise WallGoError("alpha can not be found", data={"vw": vw}) from exc
 
-    def _dxiAndWdv(self, v: float, xiAndW: np.ndarray) -> np.ndarray:
+    def _dxiAndWdv(
+            self, v: float, xiAndW: np.ndarray, shockWave: bool=True
+        ) -> np.ndarray:
         """
         Fluid equations in the shock wave as a function of v.
         """
         xi, w = xiAndW
+        if shockWave:
+            csq = self.cs2
+        else:
+            csq = self.cb2
         muXiV = (xi - v) / (1 - xi * v)
-        dwdv = w * (1 + 1 / self.cs2) * muXiV / (1 - v**2)
+        dwdv = w * (1 + 1 / csq) * muXiV / (1 - v**2)
         if v != 0:
             dxidv = (
-                xi * (1 - v * xi) * (muXiV**2 / self.cs2 - 1) / (2 * v * (1 - v**2))
+                xi * (1 - v * xi) * (muXiV**2 / csq - 1) / (2 * v * (1 - v**2))
             )
         else:
             # If v = 0, dxidv is set to a very high value
             dxidv = 1e50
         return np.array([dxidv, dwdv])
 
-    def integratePlasma(self, v0: float, vw: float, wp: float) -> OptimizeResult:
+    def integratePlasma(
+            self, v0: float, vw: float, wp: float, shockWave: bool=True
+        ) -> OptimizeResult:
         """
         Integrates the fluid equations in the shock wave until it reaches the shock
         front.
@@ -358,6 +366,9 @@ def integratePlasma(self, v0: float, vw: float, wp: float) -> OptimizeResult:
             Wall velocity.
         wp : float
             Enthalpy just in front of the wall.
+        shockWave : bool, optional
+            If True, the integration happens in the shock wave. If False, it happens in
+            the rarefaction wave. Default is True.
 
         Returns
         -------
@@ -366,15 +377,16 @@ def integratePlasma(self, v0: float, vw: float, wp: float) -> OptimizeResult:
 
         """
 
-        def event(v: float, xiAndW: np.ndarray) -> float:
+        def event(v: float, xiAndW: np.ndarray, shockWave: bool=True) -> float:
             # Function that is 0 at the shock wave front. Is used by solve_ivp to stop
             # the integration
             xi = xiAndW[0]
             return float((xi * (xi - v) / (1 - xi * v) - self.cs2) * v)
 
-        event.terminal = True
+        event.terminal = shockWave
         sol = solve_ivp(
-            self._dxiAndWdv, (v0, 1e-10), [vw, wp], events=event, rtol=self.rtol, atol=0
+            self._dxiAndWdv, (v0, 1e-10), [vw, wp],
+            events=event, rtol=self.rtol/10, atol=0, args=(shockWave,)
         )
         return sol
 
@@ -673,3 +685,58 @@ def detonationVAndT(self, vw: float) -> tuple[float, ...]:
         vm = (part + np.sqrt(part**2 - 4 * self.cb2 * vp**2)) / (2 * vp)
         Tm = self._findTm(vm, vp, self.Tnucl)
         return vp, vm, self.Tnucl, Tm
+
+    def efficiencyFactor(self, vw: float) -> float:
+        r"""
+        Computes the efficiency factor
+        :math:`\kappa=\frac{4}{v_w^3 \alpha_n w_n}\int d\xi \xi^2 v^2\gamma^2 w`..
+
+        Parameters
+        ----------
+        vw : float
+            Wall velocity.
+
+        Returns
+        -------
+        float
+            Value of the efficiency factor :math:`\kappa`.
+
+        """
+        # Separates the efficiency factor into a contribution from the shock wave and
+        # the rarefaction wave.
+        kappaSW = 0.0
+        kappaRW = 0.0
+
+        vp, vm, Tp, Tm = self.findMatching(vw)
+
+        # Computes the enthalpies (in units where w_s(Tn)=1)
+        wp = (Tp/self.Tnucl)**self.mu
+        wm = gammaSq(vp)*vp*wp/(gammaSq(vm)*vm)
+
+        # If deflagration or hybrid, computes the shock wave contribution
+        if vw < self.vJ:
+            solShock = self.integratePlasma(boostVelocity(vw, vp), vw, wp)
+            vPlasma = solShock.t
+            xi = solShock.y[0]
+            enthalpy = solShock.y[1]
+
+            # Integrate the solution to get kappa
+            kappaSW = 4 * simpson(
+                y=xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy,
+                x=xi
+            ) / (vw**3 * self.alN)
+
+        # If hybrid or detonation, computes the rarefaction wave contribution
+        if vw > self.cb:
+            solRarefaction = self.integratePlasma(boostVelocity(vw, vm), vw, wm, False)
+            vPlasma = solRarefaction.t
+            xi = solRarefaction.y[0]
+            enthalpy = solRarefaction.y[1]
+
+            # Integrate the solution to get kappa
+            kappaRW = -4 * simpson(
+                y=xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy,
+                x=xi
+            ) / (vw**3 * self.alN)
+
+        return kappaSW + kappaRW

tests/test_HydroTemplateModel.py

@@ -118,6 +118,22 @@ def test_findvwLTE():
         res2[i] = hydroTemplate.findvwLTE()
     np.testing.assert_allclose(res1,res2,rtol = 10**-4,atol = 0)
 
+def test_efficiencyFactor():
+    res1,res2 = np.zeros(N),np.zeros(N)
+    psiN = 1-0.5*rng.random(N)
+    alN = (1-psiN)/3+10**(-3*rng.random(N)-0.5)
+    cs2 = 1/4+(1/3-1/4)*rng.random(N)
+    cb2 = cs2-(1/3-1/4)*rng.random(N)
+    for i in range(N):
+        model = TestModelTemplate(alN[i],psiN[i],cb2[i],cs2[i],1,1)
+        hydrodynamics = WallGo.Hydrodynamics(model,tmax,tmin,1e-8,1e-8)
+        hydroTemplate = WallGo.HydrodynamicsTemplateModel(model)
+        vMin = max(hydroTemplate.vMin, hydrodynamics.vMin)
+        vw = vMin + (1-vMin)*rng.random()
+        res1[i] = hydrodynamics.efficiencyFactor(vw)
+        res2[i] = hydroTemplate.efficiencyFactor(vw)
+    np.testing.assert_allclose(res1,res2,rtol = 10**-2,atol = 0)
+
 def test_findHydroBoundaries():
     res1,res2 = np.zeros((N,5)),np.zeros((N,5))
     psiN = 1-0.5*rng.random(N)

tests/test_Hydrodynamics.py

@@ -218,6 +218,23 @@ def test_findMatching():
     res = hydrodynamics.findMatching(0.9)
     np.testing.assert_allclose(res,(0.9, 0.894957,0.9,0.918446),rtol = 10**-2,atol = 0)
 
+
+# Test efficiency factor in two-step model
+def test_efficiencyFactor():
+    model1 = TestModel2Step(0.2,0.1,0.4,0.7)
+    hydrodynamics = WallGo.Hydrodynamics(model1, tmax, tmin, rtol, atol)
+    res = hydrodynamics.efficiencyFactor(0.4)
+    assert res == pytest.approx(0.140972, rel = 0.01)
+    res = hydrodynamics.efficiencyFactor(0.6)
+    assert res == pytest.approx(0.334666, rel = 0.01)
+    model1 = TestModel2Step(0.2,0.1,0.4,0.9)
+    hydrodynamics = WallGo.Hydrodynamics(model1, tmax, tmin, rtol, atol)
+    res = hydrodynamics.efficiencyFactor(0.4)
+    assert res == pytest.approx(0.0399354, rel = 0.01)
+    res = hydrodynamics.efficiencyFactor(0.6)
+    assert res == pytest.approx(0.197849, rel = 0.01)
+
+
 ### Test local thermal equilibrium solution in bag model