Merge pull request #354 from RocketPy-Team/enh/liquid-motors-optimiza…

…tion

Enh/liquid motors optimization

rocketpy/Function.py

@@ -188,7 +188,7 @@ def source(x):
         if callable(source):
             # Set source
             self.source = source
-            # Set geValueOpt2
+            # Set getValueOpt
             self.getValueOpt = source
             # Set arguments name and domain dimensions
             parameters = signature(source).parameters
@@ -543,7 +543,7 @@ def setDiscrete(
             self.__interpolation__ = "shepard"
         return self
 
-    def setDiscreteBasedOnModel(self, modelFunction, oneByOne=True):
+    def setDiscreteBasedOnModel(self, modelFunction, oneByOne=True, keepSelf=True):
         """This method transforms the domain of Function instance into a list of
         discrete points based on the domain of a model Function instance. It does so by
         retrieving the domain, domain name, interpolation method and extrapolation
@@ -564,6 +564,11 @@ def setDiscreteBasedOnModel(self, modelFunction, oneByOne=True):
             If True, evaluate Function in each sample point separately. If
             False, evaluates Function in vectorized form. Default is True.
 
+        keepSelf : boolean, optional
+            If True, the original Function interpolation and extrapolation methods
+            will be kept. If False, those are substituted by the ones from the model
+            Function. Default is True.
+
         Returns
         -------
         self : Function
@@ -635,8 +640,12 @@ def setDiscreteBasedOnModel(self, modelFunction, oneByOne=True):
             Zs = np.array(self.getValue(mesh))
             self.setSource(np.concatenate(([Xs], [Ys], [Zs])).transpose())
 
-        self.setInterpolation(self.__interpolation__)
-        self.setExtrapolation(self.__extrapolation__)
+        interp = self.__interpolation__ if keepSelf else modelFunction.__interpolation__
+        extrap = self.__extrapolation__ if keepSelf else modelFunction.__extrapolation__
+
+        self.setInterpolation(interp)
+        self.setExtrapolation(extrap)
+
         return self
 
     def reset(
@@ -1714,7 +1723,7 @@ def __interpolateAkima__(self):
         self.__akimaCoefficients__ = coeffs
 
     def __neg__(self):
-        """Negates the Function objetive. The result has the same effect as
+        """Negates the Function object. The result has the same effect as
         multiplying the Function by -1.
 
         Returns
@@ -1770,15 +1779,18 @@ def __ge__(self, other):
                     # Other is lambda based Function
                     return self.yArray >= other(self.xArray)
                 except ValueError:
-                    raise ValueError("Operands should have the same discretization.")
+                    raise ValueError(
+                        "Comparison not supported between instances of the "
+                        "Function class with different domain discretization."
+                    )
             else:
                 # Other is not a Function
                 try:
                     return self.yArray >= other
                 except TypeError:
                     raise TypeError(
                         "Comparison not supported between instances of "
-                        f"'Function' and '{type(other)}'"
+                        f"'Function' and '{type(other)}'."
                     )
         else:
             # self is lambda based Function
@@ -1787,8 +1799,8 @@ def __ge__(self, other):
                     return self(other.xArray) >= other.yArray
                 except AttributeError:
                     raise TypeError(
-                        "Cannot compare lambda based Function with "
-                        "lambda based Function."
+                        "Comparison not supported between two instances of "
+                        "the Function class with callable sources."
                     )
 
     def __le__(self, other):
@@ -1829,7 +1841,7 @@ def __le__(self, other):
                 except TypeError:
                     raise TypeError(
                         "Comparison not supported between instances of "
-                        f"'Function' and '{type(other)}'"
+                        f"'Function' and '{type(other)}'."
                     )
         else:
             # self is lambda based Function
@@ -1838,8 +1850,8 @@ def __le__(self, other):
                     return self(other.xArray) <= other.yArray
                 except AttributeError:
                     raise TypeError(
-                        "Cannot compare lambda based Function with "
-                        "lambda based Function."
+                        "Comparison not supported between two instances of "
+                        "the Function class with callable sources."
                     )
 
     def __gt__(self, other):
@@ -1911,7 +1923,7 @@ def __add__(self, other):
         # If other is Function try...
         try:
             # Check if Function objects source is array or callable
-            # Check if Function objects have same interpolation and domain
+            # Check if Function objects have the same domain discretization
             if (
                 isinstance(other.source, np.ndarray)
                 and isinstance(self.source, np.ndarray)
@@ -2037,7 +2049,7 @@ def __mul__(self, other):
         # If other is Function try...
         try:
             # Check if Function objects source is array or callable
-            # Check if Function objects have same interpolation and domain
+            # Check if Function objects have the same domain discretization
             if (
                 isinstance(other.source, np.ndarray)
                 and isinstance(self.source, np.ndarray)
@@ -2120,15 +2132,15 @@ def __truediv__(self, other):
         # If other is Function try...
         try:
             # Check if Function objects source is array or callable
-            # Check if Function objects have same interpolation and domain
+            # Check if Function objects have the same domain discretization
             if (
                 isinstance(other.source, np.ndarray)
                 and isinstance(self.source, np.ndarray)
                 and self.__domDim__ == other.__domDim__
                 and np.array_equal(self.xArray, other.xArray)
             ):
                 # Operate on grid values
-                with np.errstate(divide="ignore"):
+                with np.errstate(divide="ignore", invalid="ignore"):
                     Ys = self.source[:, 1] / other.source[:, 1]
                     Ys = np.nan_to_num(Ys)
                 Xs = self.source[:, 0]
@@ -2141,7 +2153,7 @@ def __truediv__(self, other):
                 # Create new Function object
                 return Function(source, inputs, outputs, interpolation)
             else:
-                return Function(lambda x: (self.getValueOpt2(x) / other(x)))
+                return Function(lambda x: (self.getValueOpt(x) / other(x)))
         # If other is Float except...
         except AttributeError:
             if isinstance(other, (float, int, complex)):
@@ -2159,10 +2171,10 @@ def __truediv__(self, other):
                     # Create new Function object
                     return Function(source, inputs, outputs, interpolation)
                 else:
-                    return Function(lambda x: (self.getValueOpt2(x) / other))
+                    return Function(lambda x: (self.getValueOpt(x) / other))
             # Or if it is just a callable
             elif callable(other):
-                return Function(lambda x: (self.getValueOpt2(x) / other(x)))
+                return Function(lambda x: (self.getValueOpt(x) / other(x)))
 
     def __rtruediv__(self, other):
         """Divides 'other' by a Function object and returns a new Function
@@ -2194,10 +2206,10 @@ def __rtruediv__(self, other):
                 # Create new Function object
                 return Function(source, inputs, outputs, interpolation)
             else:
-                return Function(lambda x: (other / self.getValueOpt2(x)))
+                return Function(lambda x: (other / self.getValueOpt(x)))
         # Or if it is just a callable
         elif callable(other):
-            return Function(lambda x: (other(x) / self.getValueOpt2(x)))
+            return Function(lambda x: (other(x) / self.getValueOpt(x)))
 
     def __pow__(self, other):
         """Raises a Function object to the power of 'other' and
@@ -2223,7 +2235,7 @@ def __pow__(self, other):
         # If other is Function try...
         try:
             # Check if Function objects source is array or callable
-            # Check if Function objects have same interpolation and domain
+            # Check if Function objects have the same domain discretization
             if (
                 isinstance(other.source, np.ndarray)
                 and isinstance(self.source, np.ndarray)
@@ -2243,7 +2255,7 @@ def __pow__(self, other):
                 # Create new Function object
                 return Function(source, inputs, outputs, interpolation)
             else:
-                return Function(lambda x: (self.getValueOpt2(x) ** other(x)))
+                return Function(lambda x: (self.getValueOpt(x) ** other(x)))
         # If other is Float except...
         except AttributeError:
             if isinstance(other, (float, int, complex)):
@@ -2484,7 +2496,7 @@ def identityFunction(self):
         follows the same discretization, and has linear interpolation and
         extrapolation.
         If the Function is defined by a lambda, the identity Function is the
-        indentity map 'lambda x: x'.
+        identity map 'lambda x: x'.
 
         Returns
         -------
@@ -2494,15 +2506,13 @@ def identityFunction(self):
 
         # Check if Function object source is array
         if isinstance(self.source, np.ndarray):
-            identity = Function(
-                [(-1, -1), (1, 1)],
+            return Function(
+                np.column_stack((self.xArray, self.xArray)),
                 inputs=self.__inputs__,
                 outputs=f"identity of {self.__outputs__}",
                 interpolation="linear",
                 extrapolation="natural",
             )
-            return identity.setDiscreteBasedOnModel(self)
-
         else:
             return Function(
                 lambda x: x,
@@ -2840,9 +2850,9 @@ def __new__(
         source,
         inputs=["Scalar"],
         outputs=["Scalar"],
-        interpolation="akima",
+        interpolation="spline",
         extrapolation=None,
-        datapoints=50,
+        datapoints=100,
     ):
         """
         Creates a piecewise function from a dictionary of functions. The keys of the dictionary

rocketpy/motors/HybridMotor.py

@@ -330,14 +330,11 @@ def inertiaTensor(self):
             self.liquid.mass * (self.liquid.centerOfMass - self.centerOfMass) ** 2
         )
 
-        solidInertia = self.solid.inertiaTensor
-        liquidInertia = self.liquid.inertiaTensor
-
         self.InertiaI = (
-            solidInertia[0] + solidCorrection + liquidInertia[0] + liquidCorrection
+            self.solid.I_11 + solidCorrection + self.liquid.I_11 + liquidCorrection
         )
         self.InertiaZ = (
-            solidInertia[2] + solidCorrection + liquidInertia[2] + liquidCorrection
+            self.solid.I_33 + solidCorrection + self.liquid.I_33 + liquidCorrection
         )
 
         # Set naming convention

rocketpy/motors/LiquidMotor.py

@@ -183,7 +183,7 @@ def inertiaTensor(self):
         tuple (of Functions)
             Pricipal moment of inertia tensor of the motor, in kg*m^2.
         """
-        self.inertiaI = self.inertiaZ = Function(0)
+        self.inertiaI = self.inertiaZ = 0
         centerOfMass = self.centerOfMass
 
         for positioned_tank in self.positioned_tanks:
@@ -226,7 +226,7 @@ def I_11(self):
         ----------
         .. [1] https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor
         """
-        I_11 = Function(0)
+        I_11 = 0
         centerOfMass = self.centerOfMass
 
         for positioned_tank in self.positioned_tanks:

rocketpy/motors/SolidMotor.py

@@ -343,7 +343,7 @@ def centerOfMass(self):
         t : float
             Time in seconds.
 
-        Returns
+        Returnsg
         -------
         Function
             Position of the center of mass as a function
@@ -399,7 +399,15 @@ def terminateBurn(t, y):
         terminateBurn.terminal = True
 
         # Solve the system of differential equations
-        sol = integrate.solve_ivp(geometryDot, t_span, y0, events=terminateBurn)
+        sol = integrate.solve_ivp(
+            geometryDot,
+            t_span,
+            y0,
+            events=terminateBurn,
+            atol=1e-12,
+            rtol=1e-11,
+            method="LSODA",
+        )
 
         self.grainBurnOut = sol.t[-1]