Add Giese solver option to the fluid solver

pttools/bubble/bubble.py

@@ -46,9 +46,14 @@ def __init__(
             t_end: float = const.T_END_DEFAULT,
             n_xi: int = const.N_XI_DEFAULT,
             thin_shell_t_points_min: int = const.THIN_SHELL_T_POINTS_MIN,
+            use_bag_solver: bool = False,
+            use_giese_solver: bool = False,
             log_success: bool = False,
             allow_invalid: bool = False,
             log_invalid: bool = True):
+        if use_bag_solver and use_giese_solver:
+            raise ValueError("Both bag and Giese solvers cannot be used at the same time.")
+
         if v_wall < 0 or v_wall > 1:
             raise ValueError(f"Invalid v_wall={v_wall}")
 
@@ -123,6 +128,8 @@ def __init__(
         self.negative_net_entropy_change = False
         self.numerical_error = False
         self.unphysical_alpha_plus = False
+        self.use_bag_solver = use_bag_solver
+        self.use_giese_solver = use_giese_solver
 
         # LaTeX labels are not supported in Plotly 3D plots.
         # https://github.com/plotly/plotly.js/issues/608
@@ -224,13 +231,19 @@ def solve(
             sum_rtol_error: float = 5e-2,
             error_prec: str = ".4f",
             use_bag_solver: bool = False,
+            use_giese_solver: bool = False,
             log_high_alpha_n_failures: bool = True,
             log_negative_entropy: bool = True):
         if self.solved:
             msg = "Re-solving an already solved bubble! Already computed quantities will not be updated due to caching."
             logger.warning(msg)
             self.add_note(msg)
 
+        use_bag_solver = self.use_bag_solver or use_bag_solver
+        use_giese_solver = self.use_giese_solver or use_giese_solver
+        if use_bag_solver and use_giese_solver:
+            raise ValueError("Both bag and Giese solvers cannot be used at the same time.")
+
         alpha_n_max_bag = alpha_n_max_deflagration_bag(self.v_wall)
         high_alpha_n = alpha_n_max_bag - self.alpha_n < 0.05
 
@@ -246,7 +259,7 @@ def solve(
                     wn=self.wn,
                     alpha_n_max_bag=alpha_n_max_bag,
                     high_alpha_n=high_alpha_n, t_end=self.t_end, n_xi=self.n_xi, thin_shell_limit=self.thin_shell_t_points_min,
-                    use_bag_solver=use_bag_solver,
+                    use_bag_solver=use_bag_solver, use_giese_solver=use_giese_solver,
                     log_success=self.log_success, log_high_alpha_n_failures=log_high_alpha_n_failures
                 )
             self.sn = self.model.s(self.wn, Phase.SYMMETRIC)

pttools/bubble/fluid.py

@@ -26,6 +26,11 @@
 
 logger = logging.getLogger(__name__)
 
+try:
+    from giese.lisa import kappaNuMuModel
+except ImportError:
+    kappaNuMuModel = None
+
 # The output consists of:
 # v, w, xi
 # vp, vm, vp_tilde, vm_tilde, v_sh, vm_sh, vm_tilde_sh, wp, wm, wm_sh
@@ -863,6 +868,7 @@ def sound_shell_generic(
             reverse: bool = False,
             allow_failure: bool = False,
             use_bag_solver: bool = False,
+            use_giese_solver: bool = False,
             log_success: bool = True,
             log_high_alpha_n_failures: bool = False
         ) -> tp.Tuple[
@@ -872,13 +878,17 @@ def sound_shell_generic(
 
     In most cases you should not have to call this directly. Create a Bubble instead.
     """
+    if use_giese_solver:
+        return sound_shell_giese(
+            model=model, v_wall=v_wall, alpha_n=alpha_n, wn=wn, wn_guess=wn_guess, wm_guess=wm_guess)
+
     start_time = time.perf_counter()
     if alpha_n_max_bag is None:
         alpha_n_max_bag = alpha.alpha_n_max_deflagration_bag(v_wall)
     if high_alpha_n is None:
         high_alpha_n = alpha_n > alpha_n_max_bag
 
-    if wn is None:
+    if wn is None or np.isnan(wn):
         wn = model.w_n(alpha_n, wn_guess=wn_guess)
     # The shock curve hits v=0 here
     cs_n = np.sqrt(model.cs2(wn, Phase.SYMMETRIC))
@@ -1047,3 +1057,67 @@ def sound_shell_generic(
 
 
 fluid_shell_generic = sound_shell_generic
+
+
+def sound_shell_giese(
+            model: "Model",
+            v_wall: float,
+            alpha_n: float,
+            wn: float = None,
+            wn_guess: float = None,
+            wm_guess: float = None,
+        ) -> tp.Tuple[
+            np.ndarray, np.ndarray, np.ndarray, SolutionType,
+            float, float, float, float, float, float, float, float, float, float, float, bool, float]:
+    if kappaNuMuModel is None:
+        raise ImportError("The Giese code has to be installed to use this solver.")
+
+    start_time = time.perf_counter()
+
+    if wn is None or np.isnan(wn):
+        wn = model.w_n(alpha_n, wn_guess=wn_guess)
+    if wm_guess is None or np.isnan(wm_guess):
+        wm_guess = 1.
+
+    try:
+        kappa_theta_bar_n, v, wow, xi, mode, vp, vm = kappaNuMuModel(
+            cs2b=model.cs2(wm_guess, Phase.BROKEN),
+            cs2s=model.cs2(wn, Phase.SYMMETRIC),
+            al=model.alpha_theta_bar_n_from_alpha_n(alpha_n=alpha_n, wn=wn),
+            vw=v_wall
+           )
+    except ValueError:
+        return const.nan_arr, const.nan_arr, const.nan_arr, SolutionType.ERROR, \
+            np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, True, time.perf_counter() - start_time
+
+    if mode == 0:
+        sol_type = SolutionType.SUB_DEF
+    elif mode == 1:
+        sol_type = SolutionType.HYBRID
+    elif mode == 2:
+        sol_type = SolutionType.DETON
+    else:
+        raise ValueError("Got invalid mode from Giese solver:", mode)
+    w: np.ndarray = wow * wn
+
+    # Velocities in the wall frame
+    vp_tilde: float = relativity.lorentz(xi=v_wall, v=vp)
+    vm_tilde: float = relativity.lorentz(xi=v_wall, v=vm)
+
+    # Shock
+    v_sh: float = xi[-3]
+    vm_sh: float = v[-3]
+    vm_tilde_sh: float = relativity.lorentz(xi=v_sh, v=vm_sh)
+    wm_sh: float = w[-3]
+
+    # Enthalpies
+    i_wall = np.argmax(v)
+    wp: float = w[i_wall]
+    wm: float = w[i_wall - 1]
+
+    # Other
+    v_cj = chapman_jouguet.v_chapman_jouguet(model, alpha_n=alpha_n, wn=wn, wn_guess=wn_guess)
+    solution_found = True
+    elapsed = time.perf_counter() - start_time
+
+    return v, w, xi, sol_type, vp, vm, vp_tilde, vm_tilde, v_sh, vm_sh, vm_tilde_sh, wp, wm, wm_sh, v_cj, not solution_found, elapsed