Merge pull request #211 from NCAR/main

Version 2.3.1

src/acom_music_box/__init__.py

@@ -4,7 +4,7 @@
 This package contains modules for handling various aspects of a music box,
 including species, products, reactants, reactions, and more.
 """
-__version__ = "2.3.0"
+__version__ = "2.3.1"
 
 from .utils import convert_time, convert_pressure, convert_temperature, convert_concentration
 from .species import Species

src/acom_music_box/main.py

@@ -1,6 +1,6 @@
 import os
 import argparse
-from acom_music_box import MusicBox, Examples
+from acom_music_box import MusicBox, Examples, __version__ 
 import datetime
 import sys
 import logging
@@ -40,6 +40,11 @@ def parse_arguments():
         default=0,
         help='Increase logging verbosity. Use -v for info, -vv for debug.'
     )
+    parser.add_argument(
+        '--version',
+        action='version',
+        version=f'MusicBox {__version__}',
+    )
     parser.add_argument(
         '--color-output',
         action='store_true',

src/acom_music_box/music_box.py

@@ -199,7 +199,9 @@ def generateConfig(self, directory):
             elif reaction_rate.reaction.reaction_type == "LOSS":
                 name = "LOSS." + reaction_rate.reaction.name + ".s-1"
             elif reaction_rate.reaction.reaction_type == "EMISSION":
-                name = "EMISSION." + reaction_rate.reaction.name + ".s-1"
+                name = "EMIS." + reaction_rate.reaction.name + ".s-1"
+            elif reaction_rate.reaction.reaction_type == "USER_DEFINED":
+                name = "USER." + reaction_rate.reaction.name + ".s-1"
 
             reaction_names.append(name)
             reaction_rates.append(reaction_rate.rate)
@@ -425,7 +427,7 @@ def create_solver(
         # Create a solver object using the configuration file
         self.solver = musica.create_solver(
             path_to_config,
-            musica.micmsolver.rosenbrock,
+            solver_type,
             number_of_grid_cells)
 
     def solve(self, output_path=None):
@@ -469,6 +471,7 @@ def solve(self, output_path=None):
         headers.append("time")
         headers.append("ENV.temperature")
         headers.append("ENV.pressure")
+        headers.append("ENV.number_density_air")
 
         if (self.solver is None):
             raise Exception("Error: MusicBox object {} has no solver."
@@ -500,33 +503,20 @@ def solve(self, output_path=None):
 
         while (curr_time <= self.box_model_options.simulation_length):
 
-            # outputs to output_array if enough time has elapsed
-            if (next_output_time <= curr_time):
-                row = []
-                row.append(next_output_time)
-                row.append(curr_conditions.temperature)
-                row.append(curr_conditions.pressure)
-                for conc in ordered_concentrations:
-                    row.append(conc)
-                output_array.append(row)
-                next_output_time += self.box_model_options.output_step_time
-
             # iterates evolving  conditions if enough time has elapsed
             while (
                     next_conditions is not None and next_conditions_time <= curr_time):
 
                 curr_conditions.update_conditions(next_conditions)
+                ordered_rate_constants = self.order_reaction_rates(
+                    curr_conditions, rate_constant_ordering)
 
                 # iterates next_conditions if there are remaining evolving
                 # conditions
                 if (len(self.evolving_conditions) > next_conditions_index + 1):
                     next_conditions_index += 1
                     next_conditions = self.evolving_conditions.conditions[next_conditions_index]
                     next_conditions_time = self.evolving_conditions.times[next_conditions_index]
-
-                    ordered_rate_constants = self.order_reaction_rates(
-                        curr_conditions, rate_constant_ordering)
-
                 else:
                     next_conditions = None
 
@@ -537,20 +527,40 @@ def solve(self, output_path=None):
             air_density = curr_conditions.pressure / \
                 (GAS_CONSTANT * curr_conditions.temperature)
 
+            # outputs to output_array if enough time has elapsed
+            if (next_output_time <= curr_time):
+                row = []
+                row.append(next_output_time)
+                row.append(curr_conditions.temperature)
+                row.append(curr_conditions.pressure)
+                row.append(air_density)
+                for conc in ordered_concentrations:
+                    row.append(conc)
+                output_array.append(row)
+                next_output_time += self.box_model_options.output_step_time
+
+            # ensure the time step is not greater than the next update to the
+            # evolving conditions or the next output time
+            time_step = self.box_model_options.chem_step_time
+            if (next_conditions is not None and next_conditions_time > curr_time):
+                time_step = min(time_step, next_conditions_time - curr_time)
+            if (next_output_time > curr_time):
+                time_step = min(time_step, next_output_time - curr_time)
+
             # solves and updates concentration values in concentration array
             if (not ordered_concentrations):
                 logger.info("Warning: ordered_concentrations list is empty.")
             musica.micm_solve(
                 self.solver,
-                self.box_model_options.chem_step_time,
+                time_step,
                 curr_conditions.temperature,
                 curr_conditions.pressure,
                 air_density,
                 ordered_concentrations,
                 ordered_rate_constants)
 
             # increments time
-            curr_time += self.box_model_options.chem_step_time
+            curr_time += time_step
 
         df = pd.DataFrame(output_array[1:], columns=output_array[0])
         # outputs to file if output is present
@@ -687,16 +697,6 @@ def speciesOrdering(self):
         """
         return musica.species_ordering(self.solver)
 
-    def userDefinedReactionRates(self):
-        """
-        Retrieves the user-defined reaction rates from the solver.
-
-        This function calls the `user_defined_reaction_rates` function from the `musica` module,
-        passing the solver instance from the current object.
-
-        Returns:
-            dict: The dictionary of user-defined reaction rates used in the solver.
-        """
     @classmethod
     def order_reaction_rates(self, curr_conditions, rate_constant_ordering):
         """
@@ -721,6 +721,8 @@ def order_reaction_rates(self, curr_conditions, rate_constant_ordering):
                 key = "LOSS." + rate.reaction.name
             elif (rate.reaction.reaction_type == "EMISSION"):
                 key = "EMIS." + rate.reaction.name
+            elif (rate.reaction.reaction_type == "USER_DEFINED"):
+                key = "USER." + rate.reaction.name
             rate_constants[key] = rate.rate
 
         ordered_rate_constants = len(rate_constants.keys()) * [0.0]

src/acom_music_box/reaction.py

@@ -75,6 +75,7 @@ def short_type(self):
             "ARRHENIUS": "ARRH",
             "TUNNELING": "TUNN",
             "TROE_TERNARY": "TROE",
+            "USER_DEFINED": "USER",
         }
         return type_map.get(self.reaction_type, "UNKNOWN")
 

tests/test_chapman.py

@@ -58,7 +58,7 @@ def test_run(self):
                 assert math.isclose(
                     float(model_output_concs[i][j]),
                     float(test_output_concs[i][j]),
-                    rel_tol=1e-8,
+                    rel_tol=1e-7,
                     abs_tol=1e-15,
                 ), f"Arrays differ at index ({i}, {j}) for species {concs_to_test[j]}"