Merge pull request #34 from softwareengineerprogrammer/fix-list-param…

…eters

Fix list parameters

src/geophires_x/GeoPHIRESUtils.py

@@ -554,7 +554,7 @@ def read_input_file(return_dict_1, logger=None, input_file_name=None):
                     comment = comment + elements[i]
 
             # done with parsing, now create the object and add to the dictionary
-            p_entry = ParameterEntry(description, s_val, comment)
+            p_entry = ParameterEntry(description, s_val, comment, line)
             return_dict_1[description] = p_entry  # make the dictionary element
 
     else:

src/geophires_x/Parameter.py

@@ -43,6 +43,7 @@ class ParameterEntry:
     Name: str
     sValue: str
     Comment: Optional[str] = None
+    raw_entry: Optional[str] = None
 
 
 @dataclass
@@ -253,10 +254,10 @@ def ReadParameter(ParameterReadIn: ParameterEntry, ParamToModify, model):
     """
     ReadParameter: A method to take a single ParameterEntry object and use it to update the associated Parameter.
     Does validation as well as Unit and Currency conversion
-    :param ParameterEntry: The value the user wants to change and the value they want to change it to (as a string)
+    :param ParameterReadIn: The value the user wants to change and the value they want to change it to (as a string)
      and  any comment they provided with it (as a string) - all in one object (ParameterEntry) that is passed in
       to this method as a parameter itself (ParameterReadIn) - see ParameterEntry class for details on the fields in it
-    :type ParameterEntry: :class:`~geophires_x.Parameter.ParameterEntry`
+    :type ParameterReadIn: :class:`~geophires_x.Parameter.ParameterEntry`
     :param ParamToModify: The Parameter that will be modified (assuming it passes validation and conversion) - this is
       the object that will be modified by this method - see Parameter class for details on the fields in it
     :type ParamToModify: :class:`~geophires_x.Parameter.Parameter`
@@ -293,9 +294,9 @@ def ReadParameter(ParameterReadIn: ParameterEntry, ParamToModify, model):
 
     def default_parameter_value_message(new_val: Any, param_to_modify_name: str, default_value: Any) -> str:
         return (
-            f'Parameter given ({str(New_val)}) for {ParamToModify.Name} is the same as the default value. '
-            f'Consider removing {ParamToModify.Name} from the input file unless you wish '
-            f'to change it from the default value of ({str(ParamToModify.DefaultValue)})'
+            f'Parameter given ({str(new_val)}) for {param_to_modify_name} is the same as the default value. '
+            f'Consider removing {param_to_modify_name} from the input file unless you wish '
+            f'to change it from the default value of ({str(default_value)})'
         )
 
     if isinstance(ParamToModify, intParameter):
@@ -370,16 +371,23 @@ def default_parameter_value_message(new_val: Any, param_to_modify_name: str, def
                 model.logger.warning(msg)
             model.logger.info(f'Complete {str(__name__)}: {sys._getframe().f_code.co_name}')
             return
-        # All is good.  With a list, we have to use the last character of the Description to get the position.
-        # I.e., "Gradient 1" should yield a position = 0 ("1" - 1)
         else:
-            parts = ParameterReadIn.Name.split(' ')
-            position = int(parts[1]) - 1
-            if position >= len(ParamToModify.value):
-                ParamToModify.value.append(New_val)  # we are adding to the list, so use append
-            else:  # we are replacing a value, so pop the value we want to replace, then insert a new one
-                ParamToModify.value.pop(position)
-                ParamToModify.value.insert(position, New_val)
+            if ' ' in ParamToModify.Name:
+                # Some list parameters are read in with enumerated parameter names;  in these cases we use the last
+                # character of the description to get the position i.e., "Gradient 1" is position 0.
+                parts = ParameterReadIn.Name.split(' ')
+                position = int(parts[1]) - 1
+                if position >= len(ParamToModify.value):
+                    ParamToModify.value.append(New_val)  # we are adding to the list, so use append
+                else:  # we are replacing a value, so pop the value we want to replace, then insert a new one
+                    ParamToModify.value.pop(position)
+                    ParamToModify.value.insert(position, New_val)
+            else:
+                # In an ideal world this would be handled in ParameterEntry such that its sValue and Comment are
+                # correct; however that would only be practical if ParameterEntry had typing information to know
+                # whether to treat text after a second comma as a comment or list entry.
+                ParamToModify.value = [float(x.strip()) for x in ParameterReadIn.raw_entry.split('--')[0].split(',')[1:]
+                                       if x.strip() != '']
     elif isinstance(ParamToModify, boolParameter):
         if ParameterReadIn.sValue == "0":
             New_val = False

src/geophires_x/Reservoir.py

@@ -593,30 +593,15 @@ def read_parameters(self, model: Model) -> None:
                             # fracshape = 4  Rectangular fracture
                             ParameterToModify.value = FractureShape.RECTANGULAR
 
-                    elif ParameterToModify.Name.startswith('Gradient'):
+                    elif ParameterToModify.Name.startswith('Gradient '):
                         parts = ParameterReadIn.Name.split(' ')
                         position = int(parts[1]) - 1
                         model.reserv.gradient.value[position] = ParameterToModify.value
-                        if model.reserv.gradient.value[position] > 1.0:
-                            # TODO refactor to avoid heuristic-based unit conversions
-                            model.reserv.gradient.value[position] = model.reserv.gradient.value[
-                                                                        position] / 1000.0  # convert C/m
-                            model.reserv.gradient.CurrentUnits = TemperatureGradientUnit.DEGREESCPERM
 
-                        if model.reserv.gradient.value[position] < 1e-6:
-                            # convert 0 C/m gradients to very small number, avoids divide by zero errors later
-                            model.reserv.gradient.value[position] = 1e-6
-
-                    elif ParameterToModify.Name.startswith('Thickness'):
+                    elif ParameterToModify.Name.startswith('Thickness '):
                         parts = ParameterReadIn.Name.split(' ')
                         position = int(parts[1]) - 1
                         model.reserv.layerthickness.value[position] = ParameterToModify.value
-                        if model.reserv.layerthickness.value[position] < 100.0:
-                            model.reserv.layerthickness.value[position] = model.reserv.layerthickness.value[
-                                                                              position] * 1000.0  # convert m
-                            model.reserv.layerthickness.CurrentUnits = LengthUnit.METERS
-                        # set thickness of bottom segment to large number to override lower, unused segments
-                        model.reserv.layerthickness.value[position + 1] = 100_000.0
 
                     elif ParameterToModify.Name.startswith("Fracture Separation"):
                         self.fracsepcalc.value = self.fracsep.value
@@ -635,6 +620,31 @@ def read_parameters(self, model: Model) -> None:
 
         coerce_int_params_to_enum_values(self.ParameterDict)
 
+        for position in range(len(model.reserv.gradient.value)):
+            if model.reserv.gradient.value[position] > 1.0:
+                # TODO refactor to avoid heuristic-based unit conversions
+                model.reserv.gradient.value[position] = model.reserv.gradient.value[
+                                                            position] / 1000.0  # convert to C/m
+                model.reserv.gradient.CurrentUnits = TemperatureGradientUnit.DEGREESCPERM
+
+            if model.reserv.gradient.value[position] < 1e-6:
+                # convert 0 C/m gradients to very small number, avoids divide by zero errors later
+                model.reserv.gradient.value[position] = 1e-6
+
+        for position in range(len(model.reserv.layerthickness.value)):
+            if model.reserv.layerthickness.value[position] < 100.0:
+                # TODO refactor to avoid heuristic-based unit conversions
+                model.reserv.layerthickness.value[position] = model.reserv.layerthickness.value[
+                                                                  position] * 1000.0  # convert to m
+                model.reserv.layerthickness.CurrentUnits = LengthUnit.METERS
+
+        # set thickness of bottom segment to large number to override lower, unused segments
+        while len(model.reserv.layerthickness.value) < model.reserv.numseg.value:
+            model.reserv.layerthickness.value.append(100_000.0)
+
+        model.reserv.layerthickness.value[model.reserv.numseg.value-1] = 100_000.0
+
+
         model.logger.info(f'complete {str(__class__)}: {sys._getframe().f_code.co_name}')
 
     @lru_cache(maxsize=1024)

src/geophires_x/SBTReservoir.py

@@ -5,6 +5,7 @@
 import pandas as pd
 from scipy.special import erf, erfc, jv, yv, exp1
 from scipy.interpolate import interp1d
+import scipy.io as sio
 import matplotlib.pyplot as plt
 
 import geophires_x.Model as Model
@@ -451,7 +452,7 @@ def Calculate_Coaxial(self, model):
         g = 9.81  # Gravitational acceleration [m/s²]
         gamma = 0.577215665  # Euler's constant
         alpha_m = self.krock.value / (self.rhorock.value * self.cprock.value)  # Rock thermal diffusivity [m²/s]
-        alpha_m_boiler = self.krock.value_boiler / (self.rhorock.value * self.cprock.value)  # Boiler rock thermal diffusivity [m²/s]
+        alpha_m_boiler = self.krock.value / (self.rhorock.value * self.cprock.value)  # Boiler rock thermal diffusivity [m²/s]
 
         outerradiuscenterpipe = radiuscenterpipe + thicknesscenterpipe  # Outer radius of inner pipe [m]
         A_flow_annulus = np.pi * (radius ** 2 - outerradiuscenterpipe ** 2)  # Flow area of annulus pipe [m²]