Add a new file parser option for non-homologous expansion

tardis/io/config_reader.py

@@ -10,6 +10,7 @@
 from tardis.io import config_validator
 from tardis.io.util import YAMLLoader, yaml_load_file, HDFWriterMixin
 from tardis.io.parsers.csvy import load_yaml_from_csvy
+from tardis.montecarlo import montecarlo_configuration
 
 pp = pprint.PrettyPrinter(indent=4)
 
@@ -96,6 +97,18 @@ def __init__(self, value=None):
             if "v_outer_boundary" in csvy_yml:
                 model["v_outer_boundary"] = csvy_yml["v_outer_boundary"]
 
+            if NONHOMOLOGOUS_EXPANSION_ENABLED is True and (
+                "r_inner_boundary" not in csvy_yml
+                or "r_outer_boundary" not in csvy_yml
+            ):
+                raise ValueError(
+                    "r_inner_boundary and/or r_outer_boundary not in csvy_model file."
+                )
+            if "r_inner_boundary" in csvy_yml:
+                model["r_inner_boundary"] = csvy_yml["r_inner_boundary"]
+            if "r_outer_boundary" in csvy_yml:
+                model["r_outer_boundary"] = csvy_yml["r_outer_boundary"]
+
             self.__setitem__("model", model)
             for key in self.model:
                 self.model.__setitem__(key, self.model[key])
@@ -257,6 +270,33 @@ def from_config_dict(cls, config_dict, validate=True, config_dirname=""):
         validated_config_dict["config_dirname"] = config_dirname
 
         montecarlo_section = validated_config_dict["montecarlo"]
+        if montecarlo_section["convergence_strategy"]["type"] == "damped":
+            montecarlo_section[
+                "convergence_strategy"
+            ] = parse_convergence_section(
+                montecarlo_section["convergence_strategy"]
+            )
+        elif montecarlo_section["convergence_strategy"]["type"] == "custom":
+            raise NotImplementedError(
+                'convergence_strategy is set to "custom"; '
+                "you need to implement your specific convergence treatment"
+            )
+        else:
+            raise ValueError(
+                'convergence_strategy is not "damped" ' 'or "custom"'
+            )
+
+        enable_full_relativity = montecarlo_section["enable_full_relativity"]
+        spectrum_integrated = (
+            validated_config_dict["spectrum"]["method"] == "integrated"
+        )
+        if enable_full_relativity and spectrum_integrated:
+            raise NotImplementedError(
+                "The spectrum method is set to 'integrated' and "
+                "enable_full_rela tivity to 'True'.\n"
+                "The FormalIntegrator is not yet implemented for the full "
+                "relativity mode. "
+            )
         Configuration.validate_montecarlo_section(montecarlo_section)
 
         if "csvy_model" in validated_config_dict.keys():

tardis/io/model_reader.py

@@ -51,6 +51,7 @@ def read_density_file(filename, filetype):
         "artis": read_artis_density,
         "simple_ascii": read_simple_ascii_density,
         "cmfgen_model": read_cmfgen_density,
+        "non_homologous": read_csv_nonhomologous_density,
     }
 
     electron_densities = None
@@ -86,7 +87,6 @@ def read_density_file(filename, filetype):
         raise ConfigurationError(
             "Invalid volume of following cell(s):\n" f"{message:s}"
         )
-
     return (
         time_of_model,
         velocity,
@@ -232,6 +232,42 @@ def read_simple_ascii_density(fname):
     return time_of_model, velocity, mean_density
 
 
+def read_csv_nonhomologous_density(fname):
+    """
+    For non-homologous expansion. The same format as read_simple_ascii_density, but assumes "radius" is in the file.
+    #index velocity [km/s] density [g/cm^3]
+
+
+    Parameters
+    ----------
+    fname : str
+        filename or path with filename
+
+    Returns
+    -------
+    time_of_model : astropy.units.Quantity
+        time at which the model is valid
+    data : pandas.DataFrame
+        data frame containing index, velocity (in km/s) and density
+    """
+
+    with open(fname) as fh:
+        time_of_model_string = fh.readline().strip()
+        time_of_model = parse_quantity(time_of_model_string)
+
+    data = recfromtxt(
+        fname,
+        skip_header=1,
+        names=("index", "radius", "velocity", "density"),
+        dtype=(int, float, float, float),
+    )
+    radius = (data["radius"] * u.km).to("cm")
+    velocity = (data["velocity"] * u.km / u.s).to("cm/s")
+    mean_density = (data["density"] * u.Unit("g/cm^3"))[1:]
+
+    return time_of_model, radius, velocity, mean_density
+
+
 def read_artis_density(fname):
     """
     Reading a density file of the following structure (example; lines starting with a hash will be ignored):

tardis/io/schemas/csvy_model.yml

@@ -39,6 +39,14 @@ properties:
     type: quantity
     description: velocity of the outer boundary of the last shell
 
+  r_inner_boundary:
+    type: quantity
+    description: radius of the inner boundary
+
+  r_outer_boundary:
+    type: quantity
+    description: radius of the outer boundary of the last shell
+
   velocity:
     type: object
     properties:
@@ -49,6 +57,17 @@ properties:
       num:
         type: number
         multipleOf: 1.0
+
+  radius:
+    type: object
+    properties:
+      start:
+        type: quantity
+      stop:
+        type: quantity
+      num:
+        type: number
+        multipleOf: 1.0
 
   density:
     oneOf:

tardis/montecarlo/montecarlo_configuration.py

@@ -16,3 +16,4 @@
 VPACKET_LOGGING = False
 RPACKET_TRACKING = False
 CONTINUUM_PROCESSES_ENABLED = False
+NONHOMOLOGOUS_EXPANSION_ENABLED = False

tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py

@@ -34,7 +34,7 @@ def velocity_dvdr(r_packet, numba_model):
 @njit(**njit_dict_no_parallel)
 def tau_sobolev_factor(r_packet, numba_model):
     """
-    The angle and velocity dependent Tau Sobolev factor component. Is called when ENABLE_NONHOMOLOGOUS_EXPANSION is set to True.
+    The angle and velocity dependent Tau Sobolev factor component. Is called when NONHOMOLOGOUS_EXPANSION_ENABLED is set to True.
 
     Note: to get Tau Sobolev, this needs to be multiplied by tau_sobolevs found from plasma
     Parameters