From 4f99e7d9019b8b84eb6ef48b02e13fa7436b3d5d Mon Sep 17 00:00:00 2001 From: Benoit Date: Thu, 24 Oct 2024 14:44:45 -0400 Subject: [PATCH 1/7] Coded a function to compute the efficiency factor in Hydrodynamics --- src/WallGo/hydrodynamics.py | 81 ++++++++++++++++++++++++++++++++++++- 1 file changed, 79 insertions(+), 2 deletions(-) diff --git a/src/WallGo/hydrodynamics.py b/src/WallGo/hydrodynamics.py index cfc77320..37e710d6 100644 --- a/src/WallGo/hydrodynamics.py +++ b/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 @@ -558,7 +558,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 +903,83 @@ def shock( ) return float(sol.root) + def efficiencyFactor(self, vw: float) -> float: + r""" + Computes the efficiency factor :math:`\kappa`. + + 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-8], + 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 = self.thermodynamics.wHighT(T) + + # Integrate the solution to get kappa + kappaSW = 4 * simpson( + xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, + xi, + ) / (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-8], + xi0T0, + rtol=self.rtol, + atol=0, + ) # solve differential equation all the way from v = v- to v = 0 + vPlasma = solRarefaction.t + xi = solRarefaction.y[0] + T = solRarefaction.y[1] + enthalpy = self.thermodynamics.wLowT(T) + + # Integrate the solution to get kappa + kappaRW = 4 * simpson( + xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, + xi, + ) / (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 From c06ccb7b5fded44918298647d1f9ff78f189e817 Mon Sep 17 00:00:00 2001 From: Benoit Date: Thu, 24 Oct 2024 15:06:28 -0400 Subject: [PATCH 2/7] Coded it for the template model as well --- src/WallGo/hydrodynamics.py | 4 +- src/WallGo/hydrodynamicsTemplateModel.py | 58 +++++++++++++++++++++++- 2 files changed, 58 insertions(+), 4 deletions(-) diff --git a/src/WallGo/hydrodynamics.py b/src/WallGo/hydrodynamics.py index 37e710d6..8d3420d9 100644 --- a/src/WallGo/hydrodynamics.py +++ b/src/WallGo/hydrodynamics.py @@ -953,7 +953,7 @@ def shock(v: float, xiAndT: np.ndarray | list) -> float: kappaSW = 4 * simpson( xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, xi, - ) / (self.thermodynamics.wHighT(self.Tnucl)*self.template.alN) + ) / (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): @@ -976,7 +976,7 @@ def shock(v: float, xiAndT: np.ndarray | list) -> float: kappaRW = 4 * simpson( xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, xi, - ) / (self.thermodynamics.wHighT(self.Tnucl)*self.template.alN) + ) / (vw**3*self.thermodynamics.wHighT(self.Tnucl)*self.template.alN) return kappaSW + kappaRW diff --git a/src/WallGo/hydrodynamicsTemplateModel.py b/src/WallGo/hydrodynamicsTemplateModel.py index 3f3515f3..ddf9c930 100644 --- a/src/WallGo/hydrodynamicsTemplateModel.py +++ b/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 @@ -673,3 +673,57 @@ 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`. + + 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( + xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, + 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) + vPlasma = solRarefaction.t + xi = solRarefaction.y[0] + enthalpy = solRarefaction.y[1] + + # Integrate the solution to get kappa + kappaRW = 4 * simpson( + xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, + xi, + ) / (vw**3 * self.alN) + + return kappaSW + kappaRW From 681cd6bcb359155380c540bd2476f9fb567093d8 Mon Sep 17 00:00:00 2001 From: Benoit Date: Thu, 24 Oct 2024 17:24:08 -0400 Subject: [PATCH 3/7] Made a test for the efficiency factor --- src/WallGo/hydrodynamics.py | 33 ++++++++++++----------- src/WallGo/hydrodynamicsTemplateModel.py | 34 ++++++++++++++++-------- tests/test_HydroTemplateModel.py | 16 +++++++++++ 3 files changed, 57 insertions(+), 26 deletions(-) diff --git a/src/WallGo/hydrodynamics.py b/src/WallGo/hydrodynamics.py index 8d3420d9..f1d332b7 100644 --- a/src/WallGo/hydrodynamics.py +++ b/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, simpson +from scipy.integrate import solve_ivp, simps 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: @@ -938,7 +940,7 @@ def shock(v: float, xiAndT: np.ndarray | list) -> float: # Integrate the shock wave solShock = solve_ivp( self.shockDE, - [vpcent, 1e-8], + [vpcent, 1e-10], xi0T0, events=shock, rtol=self.rtol, @@ -947,10 +949,10 @@ def shock(v: float, xiAndT: np.ndarray | list) -> float: vPlasma = solShock.t xi = solShock.y[0] T = solShock.y[1] - enthalpy = self.thermodynamics.wHighT(T) + enthalpy = np.array([self.thermodynamics.wHighT(t) for t in T]) # Integrate the solution to get kappa - kappaSW = 4 * simpson( + kappaSW = 4 * simps( xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, xi, ) / (vw**3*self.thermodynamics.wHighT(self.Tnucl)*self.template.alN) @@ -962,18 +964,19 @@ def shock(v: float, xiAndT: np.ndarray | list) -> float: # Integrate the rarefaction wave solRarefaction = solve_ivp( self.shockDE, - [vmcent, 1e-8], + [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 = self.thermodynamics.wLowT(T) + enthalpy = np.array([self.thermodynamics.wLowT(t) for t in T]) # Integrate the solution to get kappa - kappaRW = 4 * simpson( + kappaRW = -4 * simps( xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, xi, ) / (vw**3*self.thermodynamics.wHighT(self.Tnucl)*self.template.alN) diff --git a/src/WallGo/hydrodynamicsTemplateModel.py b/src/WallGo/hydrodynamicsTemplateModel.py index ddf9c930..0f20a317 100644 --- a/src/WallGo/hydrodynamicsTemplateModel.py +++ b/src/WallGo/hydrodynamicsTemplateModel.py @@ -6,7 +6,7 @@ import warnings import logging import numpy as np -from scipy.integrate import solve_ivp, simpson +from scipy.integrate import solve_ivp, simps from scipy.optimize import root_scalar, minimize_scalar, OptimizeResult from .exceptions import WallGoError from .helpers import gammaSq, boostVelocity @@ -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 @@ -708,20 +720,20 @@ def efficiencyFactor(self, vw: float) -> float: enthalpy = solShock.y[1] # Integrate the solution to get kappa - kappaSW = 4 * simpson( + kappaSW = 4 * simps( xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, 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) + 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( + kappaRW = -4 * simps( xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, xi, ) / (vw**3 * self.alN) diff --git a/tests/test_HydroTemplateModel.py b/tests/test_HydroTemplateModel.py index 7da9a0eb..b4102d68 100644 --- a/tests/test_HydroTemplateModel.py +++ b/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) From 19e666a7fee229328d730c35053f8a488e943713 Mon Sep 17 00:00:00 2001 From: Benoit Laurent Date: Thu, 24 Oct 2024 17:37:24 -0400 Subject: [PATCH 4/7] Replaced simps by simpson --- src/WallGo/hydrodynamics.py | 6 +++--- src/WallGo/hydrodynamicsTemplateModel.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/WallGo/hydrodynamics.py b/src/WallGo/hydrodynamics.py index f1d332b7..7dd3e5a3 100644 --- a/src/WallGo/hydrodynamics.py +++ b/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, simps +from scipy.integrate import solve_ivp, simpson from .exceptions import WallGoError from .thermodynamics import Thermodynamics from .hydrodynamicsTemplateModel import HydrodynamicsTemplateModel @@ -952,7 +952,7 @@ def shock(v: float, xiAndT: np.ndarray | list) -> float: enthalpy = np.array([self.thermodynamics.wHighT(t) for t in T]) # Integrate the solution to get kappa - kappaSW = 4 * simps( + kappaSW = 4 * simpson( xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, xi, ) / (vw**3*self.thermodynamics.wHighT(self.Tnucl)*self.template.alN) @@ -976,7 +976,7 @@ def shock(v: float, xiAndT: np.ndarray | list) -> float: enthalpy = np.array([self.thermodynamics.wLowT(t) for t in T]) # Integrate the solution to get kappa - kappaRW = -4 * simps( + kappaRW = -4 * simpson( xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, xi, ) / (vw**3*self.thermodynamics.wHighT(self.Tnucl)*self.template.alN) diff --git a/src/WallGo/hydrodynamicsTemplateModel.py b/src/WallGo/hydrodynamicsTemplateModel.py index 0f20a317..2a1ea708 100644 --- a/src/WallGo/hydrodynamicsTemplateModel.py +++ b/src/WallGo/hydrodynamicsTemplateModel.py @@ -6,7 +6,7 @@ import warnings import logging import numpy as np -from scipy.integrate import solve_ivp, simps +from scipy.integrate import solve_ivp, simpson from scipy.optimize import root_scalar, minimize_scalar, OptimizeResult from .exceptions import WallGoError from .helpers import gammaSq, boostVelocity @@ -720,7 +720,7 @@ def efficiencyFactor(self, vw: float) -> float: enthalpy = solShock.y[1] # Integrate the solution to get kappa - kappaSW = 4 * simps( + kappaSW = 4 * simpson( xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, xi, ) / (vw**3 * self.alN) @@ -733,7 +733,7 @@ def efficiencyFactor(self, vw: float) -> float: enthalpy = solRarefaction.y[1] # Integrate the solution to get kappa - kappaRW = -4 * simps( + kappaRW = -4 * simpson( xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, xi, ) / (vw**3 * self.alN) From 5d8cfb7e166b22c9ffa2610a604c12df808b5ce3 Mon Sep 17 00:00:00 2001 From: Benoit Laurent Date: Thu, 24 Oct 2024 19:34:24 -0400 Subject: [PATCH 5/7] Try something to make the tests work --- src/WallGo/hydrodynamics.py | 8 ++++---- src/WallGo/hydrodynamicsTemplateModel.py | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/WallGo/hydrodynamics.py b/src/WallGo/hydrodynamics.py index 7dd3e5a3..0e4f8307 100644 --- a/src/WallGo/hydrodynamics.py +++ b/src/WallGo/hydrodynamics.py @@ -953,8 +953,8 @@ def shock(v: float, xiAndT: np.ndarray | list) -> float: # Integrate the solution to get kappa kappaSW = 4 * simpson( - xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, - xi, + 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 @@ -977,8 +977,8 @@ def shock(v: float, xiAndT: np.ndarray | list) -> float: # Integrate the solution to get kappa kappaRW = -4 * simpson( - xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, - xi, + y=xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, + x=xi ) / (vw**3*self.thermodynamics.wHighT(self.Tnucl)*self.template.alN) return kappaSW + kappaRW diff --git a/src/WallGo/hydrodynamicsTemplateModel.py b/src/WallGo/hydrodynamicsTemplateModel.py index 2a1ea708..8df8d02f 100644 --- a/src/WallGo/hydrodynamicsTemplateModel.py +++ b/src/WallGo/hydrodynamicsTemplateModel.py @@ -721,8 +721,8 @@ def efficiencyFactor(self, vw: float) -> float: # Integrate the solution to get kappa kappaSW = 4 * simpson( - xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, - xi, + y=xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, + x=xi ) / (vw**3 * self.alN) # If hybrid or detonation, computes the rarefaction wave contribution @@ -734,8 +734,8 @@ def efficiencyFactor(self, vw: float) -> float: # Integrate the solution to get kappa kappaRW = -4 * simpson( - xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, - xi, + y=xi**2*vPlasma**2*gammaSq(vPlasma)*enthalpy, + x=xi ) / (vw**3 * self.alN) return kappaSW + kappaRW From adcd03ab1be403f04d3d14bb61bf4f0a693463d3 Mon Sep 17 00:00:00 2001 From: Jorinde van de Vis Date: Mon, 28 Oct 2024 14:32:08 +0100 Subject: [PATCH 6/7] added a test for the efficiency factor in the 2-step --- tests/test_Hydrodynamics.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/tests/test_Hydrodynamics.py b/tests/test_Hydrodynamics.py index a7623adb..55e7b163 100644 --- a/tests/test_Hydrodynamics.py +++ b/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 From 32cb4b85d9a7d34005a33b706497dca029c2ef71 Mon Sep 17 00:00:00 2001 From: Benoit Laurent Date: Mon, 28 Oct 2024 14:44:02 -0400 Subject: [PATCH 7/7] Added the formula in the docstring --- src/WallGo/hydrodynamics.py | 3 ++- src/WallGo/hydrodynamicsTemplateModel.py | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/WallGo/hydrodynamics.py b/src/WallGo/hydrodynamics.py index 0e4f8307..8aeb44ce 100644 --- a/src/WallGo/hydrodynamics.py +++ b/src/WallGo/hydrodynamics.py @@ -907,7 +907,8 @@ def shock( def efficiencyFactor(self, vw: float) -> float: r""" - Computes the efficiency factor :math:`\kappa`. + 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 ---------- diff --git a/src/WallGo/hydrodynamicsTemplateModel.py b/src/WallGo/hydrodynamicsTemplateModel.py index 8df8d02f..307faa34 100644 --- a/src/WallGo/hydrodynamicsTemplateModel.py +++ b/src/WallGo/hydrodynamicsTemplateModel.py @@ -688,7 +688,8 @@ def detonationVAndT(self, vw: float) -> tuple[float, ...]: def efficiencyFactor(self, vw: float) -> float: r""" - Computes the efficiency factor :math:`\kappa`. + 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 ----------