Pseudo-evolution for all SO definitions

SOAP/compute_halo_properties.py

@@ -540,6 +540,12 @@ def compute_halo_properties():
                 "BoundSubhalo/EncloseRadius is not enabled. This means apertures with r > r_enclose will be calculated explicitly, rather than copying over values from smaller apertures"
             )
         category_filter.print_filters()
+        cosmology_errors = parameter_file.check_cosmology(cellgrid.cosmology)
+        if len(cosmology_errors):
+            print("The snapshot cosmology is incompatible with some properties:")
+            for error in cosmology_errors:
+                print(f"  {error}", flush=True)
+            comm_world.Abort(1)
 
         # Properties enabled in the parameter file must be computed, so abort
         # if the input files do not contain the datasets they require

SOAP/core/parameter_file.py

@@ -588,6 +588,30 @@ def print_invalid_properties(self, halo_prop_list) -> None:
             for base_halo_type, prop in invalid_properties:
                 print(f"  {base_halo_type}  {prop}")
 
+    def check_cosmology(self, cosmology: Dict) -> List[str]:
+        """
+        Check that the cosmology is compatible with the properties that will
+        be calculated. This must be called after all the halo types have been
+        created, since that is when the property filters are set.
+
+        Parameters:
+         - cosmology: Dict
+           Cosmology attributes read from the snapshot.
+
+        Returns a list of error messages, which is empty if there are no problems.
+        """
+        errors = []
+
+        # The pseudo-evolution correction for the flow rates assumes flat LCDM
+        SO_filters = self.property_filters.get("SOProperties", {})
+        if any(f for name, f in SO_filters.items() if name.endswith("FlowRate")):
+            if abs(cosmology["Omega_k"]) > 1e-6:
+                errors.append("SO flow rates can only be computed if Omega_k=0")
+            if (cosmology["w_0"] != -1) or (cosmology["w_a"] != 0):
+                errors.append("SO flow rates can only be computed if w_0=-1, w_a=0")
+
+        return errors
+
     def has_enabled_properties(self, base_halo_type: str) -> bool:
         """
         Return True if the parameter file enables at least one property for the

SOAP/core/swift_cells.py

@@ -142,6 +142,51 @@ def identify_datasets(filename, nr_files, ptypes, registry):
     return metadata
 
 
+def compute_virBN98(cosmology, a):
+    """
+    Compute the Bryan & Norman (1998) critical density multiple at the
+    given scale factor.
+
+    Parameters:
+     - cosmology: dict
+       Cosmology attributes read from the snapshot.
+     - a: float
+       Scale factor.
+
+    Returns the critical density multiple.
+    """
+    Omega_k = cosmology["Omega_k"]
+    Omega_Lambda = cosmology["Omega_lambda"]
+    Omega_m = cosmology["Omega_m"]
+    bnx = -(Omega_k / a**2 + Omega_Lambda) / (
+        Omega_k / a**2 + Omega_m / a**3 + Omega_Lambda
+    )
+    return 18.0 * np.pi**2 + 82.0 * bnx - 39.0 * bnx**2
+
+
+def compute_dlog_virBN98_dloga(cosmology, a, eps=1e-5):
+    """
+    Compute the logarithmic derivative of the Bryan & Norman (1998) critical
+    density multiple with respect to the scale factor. We difference the
+    expression used to set the multiple itself, so that the two cannot
+    become inconsistent.
+
+    Parameters:
+     - cosmology: dict
+       Cosmology attributes read from the snapshot.
+     - a: float
+       Scale factor.
+     - eps: float
+       Step size in log(a) used for the central difference.
+
+    Returns dlog(virBN98)/dlog(a).
+    """
+    return (
+        np.log(compute_virBN98(cosmology, a * np.exp(eps)))
+        - np.log(compute_virBN98(cosmology, a * np.exp(-eps)))
+    ) / (2 * eps)
+
+
 class SWIFTCellGrid:
     def get_unit(self, name):
         return unyt.Unit(name, registry=self.snap_unit_registry)
@@ -270,16 +315,14 @@ def __init__(
             )
 
             # Compute the BN98 critical density multiple
-            Omega_k = self.cosmology["Omega_k"]
-            Omega_Lambda = self.cosmology["Omega_lambda"]
-            Omega_m = self.cosmology["Omega_m"]
-            bnx = -(Omega_k / self.a**2 + Omega_Lambda) / (
-                Omega_k / self.a**2 + Omega_m / self.a**3 + Omega_Lambda
-            )
-            self.virBN98 = 18.0 * np.pi**2 + 82.0 * bnx - 39.0 * bnx**2
+            self.virBN98 = compute_virBN98(self.cosmology, self.a)
             if self.virBN98 < 50.0 or self.virBN98 > 1000.0:
                 raise RuntimeError("Invalid value for virBN98!")
 
+            # The BN98 density multiple is time dependent, so its logarithmic
+            # derivative is required to calculate the pseudo-evolution of R_BN98.
+            self.dlog_virBN98_dloga = compute_dlog_virBN98_dloga(self.cosmology, self.a)
+
             # Get the box size. Assume it's comoving with no h factors.
             comoving_length_unit = self.get_unit("snap_length") * self.a_unit
             self.boxsize = unyt.unyt_quantity(

SOAP/particle_selection/SO_properties.py

@@ -249,6 +249,7 @@ def __init__(
         observer_position: unyt.unyt_array,
         core_excision_fraction: float,
         virial_definition: bool,
+        compute_flow_rates: bool,
         search_radius: unyt.unyt_quantity,
         cosmology: dict,
         boxsize: unyt.unyt_quantity,
@@ -273,6 +274,9 @@ def __init__(
          - virial_definition: bool
            Whether to calculate the properties that are only valid for virial SO
            definitions
+         - compute_flow_rates: bool
+           Whether to calculate the flow rates. These are not valid for SO
+           definitions with a fixed physical radius.
          - search_radius: unyt.unyt_quantity
            Current search radius. Particles are guaranteed to be included up to
            this radius.
@@ -311,6 +315,7 @@ def __init__(
         self.observer_position = observer_position
         self.core_excision_fraction = core_excision_fraction
         self.virial_definition = virial_definition
+        self.compute_flow_rates = compute_flow_rates
         self.search_radius = search_radius
 
     def get_dataset(self, name: str) -> unyt.unyt_array:
@@ -2796,17 +2801,12 @@ def calculate_flow_rate(
             # Adding Hubble flow term
             if hubble:
                 v_r += radii[r_mask] * self.cosmology["H"]
-            # Account for expansion of R_SO
+            # Account for expansion of R_SO. The coefficient depends on the
+            # SO definition, and is set when this calculation is constructed.
             if pseudo_evolve:
-                G = unyt.Unit("newton_G", registry=masses.units.registry)
-                R_dot = (2 / 3) * (G * self.SO_mass * self.cosmology["H"] / 100) ** (
-                    1 / 3
+                v_r -= (
+                    R * self.cosmology["H"] * self.cosmology["pseudo_evolution_coeff"]
                 )
-                R_dot *= (
-                    2 * self.cosmology["Omega_g"] + (3 / 2) * self.cosmology["Omega_m"]
-                )
-                R_dot *= R_frac
-                v_r -= R_dot
 
             # Calculate different flow types
             # We want both the inflow and outflow rates to be positive values
@@ -2845,7 +2845,7 @@ def DarkMatterMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of dark matter through 3 spherical shells
         """
-        if (self.Ndm == 0) or (not self.virial_definition):
+        if (self.Ndm == 0) or (not self.compute_flow_rates):
             return None
 
         # Particles outside the SO radius are required to calculate the
@@ -2861,7 +2861,7 @@ def StellarMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of stars through 3 spherical shells
         """
-        if (self.Nstar == 0) or (not self.virial_definition):
+        if (self.Nstar == 0) or (not self.compute_flow_rates):
             return None
 
         # Particles outside the SO radius are required to calculate the
@@ -2877,7 +2877,7 @@ def HIMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of HI through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         # Particles outside the SO radius are required to calculate the
@@ -2903,7 +2903,7 @@ def H2MassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of H2 through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         # Particles outside the SO radius are required to calculate the
@@ -2931,7 +2931,7 @@ def MetalMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of metals through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         # Particles outside the SO radius are required to calculate the
@@ -2983,7 +2983,7 @@ def ColdGasMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of cold gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmax = 1.0e3 * unyt.K
@@ -2994,7 +2994,7 @@ def CoolGasMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of cool gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e3 * unyt.K
@@ -3008,7 +3008,7 @@ def WarmGasMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of warm gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e5 * unyt.K
@@ -3022,7 +3022,7 @@ def HotGasMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of hot gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e7 * unyt.K
@@ -3033,7 +3033,7 @@ def ColdGasEnergyFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the energy flow rate of cold gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmax = 1.0e3 * unyt.K
@@ -3046,7 +3046,7 @@ def CoolGasEnergyFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the energy flow rate of cool gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e3 * unyt.K
@@ -3060,7 +3060,7 @@ def WarmGasEnergyFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the energy flow rate of warm gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e5 * unyt.K
@@ -3074,7 +3074,7 @@ def HotGasEnergyFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the energy flow rate of hot gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e7 * unyt.K
@@ -3087,7 +3087,7 @@ def ColdGasMomentumFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the momentum flow rate of cold gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmax = 1.0e3 * unyt.K
@@ -3100,7 +3100,7 @@ def CoolGasMomentumFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the momentum flow rate of cool gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e3 * unyt.K
@@ -3114,7 +3114,7 @@ def WarmGasMomentumFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the momentum flow rate of warm gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e5 * unyt.K
@@ -3128,7 +3128,7 @@ def HotGasMomentumFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the momentum flow rate of hot gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e7 * unyt.K
@@ -3375,8 +3375,6 @@ def __init__(
         self.cosmology["H"] = cellgrid.cosmology[
             "H [internal units]"
         ] / cellgrid.get_unit("code_time")
-        self.cosmology["Omega_g"] = cellgrid.cosmology["Omega_g"]
-        self.cosmology["Omega_m"] = cellgrid.cosmology["Omega_m"]
 
         # This specifies how large a sphere is read in:
         # we use default values that are sufficiently small/large to avoid reading in too many particles
@@ -3397,6 +3395,31 @@ def __init__(
             self.virial_definition = True
         elif type == "physical":
             self.physical_radius_mpc = 0.001 * SOval
+        # Flow rates are not computed for a fixed physical radius, since it
+        # does not pseudo-evolve
+        self.compute_flow_rates = type != "physical"
+
+        # Coefficient used to correct the flow rates for the pseudo-evolution
+        # of the SO radius: Rdot = coeff * R * H, where
+        # coeff = -(1/3) dln(rho_ref)/dln(a) at fixed SO mass, and rho_ref is
+        # the reference density used to define the SO radius.
+        # Derivation is in documentation/pseudo_evolution.pdf
+        # The Omega values in the snapshot are z=0 values, so we scale them.
+        H0_over_H_sq = (
+            cellgrid.cosmology["H0 [internal units]"]
+            / cellgrid.cosmology["H [internal units]"]
+        ) ** 2
+        Omega_m = float(cellgrid.mean_density / cellgrid.critical_density)
+        Omega_r = cellgrid.cosmology["Omega_r"] * H0_over_H_sq / cellgrid.a**4
+        one_plus_q = 2 * Omega_r + 1.5 * Omega_m
+        if type == "mean":
+            self.cosmology["pseudo_evolution_coeff"] = 1.0
+        elif type == "crit":
+            self.cosmology["pseudo_evolution_coeff"] = (2 / 3) * one_plus_q
+        elif type == "BN98":
+            self.cosmology["pseudo_evolution_coeff"] = (
+                2 * one_plus_q - cellgrid.dlog_virBN98_dloga
+            ) / 3
 
         # Give this calculation a name so we can select it on the command line
         if type in ["mean", "crit"]:
@@ -3573,6 +3596,7 @@ def calculate(
                 self.observer_position,
                 self.core_excision_fraction,
                 self.virial_definition,
+                self.compute_flow_rates,
                 search_radius,
                 self.cosmology,
                 self.boxsize,