Merge pull request #422 from RocketPy-Team/enh/solid-ode-jacobian

FIX: Motors UnitTesting Warnings

Author: MateusStano

rocketpy/mathutils/function.py

@@ -2242,7 +2242,7 @@ def integral(self, a, b, numerical=False):
             ans = np.trapz(y_integration_data, x_integration_data)
         else:
             # Integrate numerically
-            ans, _ = integrate.quad(self, a, b, epsabs=0.001, limit=10000)
+            ans, _ = integrate.quad(self, a, b, epsabs=1e-4, epsrel=1e-3, limit=1000)
         return integration_sign * ans
 
     def differentiate(self, x, dx=1e-6, order=1):

rocketpy/motors/solid_motor.py

@@ -469,18 +469,35 @@ def evaluate_geometry(self):
 
         density = self.grain_density
         rO = self.grain_outer_radius
+        n_grain = self.grain_number
 
         # Define system of differential equations
         def geometry_dot(t, y):
-            grain_mass_dot = self.mass_flow_rate(t) / self.grain_number
+            # Store physical parameters
+            volume_diff = self.mass_flow_rate(t) / (n_grain * density)
+
+            # Compute state vector derivative
             rI, h = y
-            rIDot = (
-                -0.5 * grain_mass_dot / (density * np.pi * (rO**2 - rI**2 + rI * h))
-            )
-            hDot = (
-                1.0 * grain_mass_dot / (density * np.pi * (rO**2 - rI**2 + rI * h))
-            )
-            return [rIDot, hDot]
+            burn_area = 2 * np.pi * (rO**2 - rI**2 + rI * h)
+            rI_dot = -volume_diff / burn_area
+            h_dot = -2 * rI_dot
+
+            return [rI_dot, h_dot]
+
+        # Define jacobian of the system of differential equations
+        def geometry_jacobian(t, y):
+            # Store physical parameters
+            volume_diff = self.mass_flow_rate(t) / (n_grain * density)
+
+            # Compute jacobian
+            rI, h = y
+            factor = volume_diff / (2 * np.pi * (rO**2 - rI**2 + rI * h) ** 2)
+            drI_dot_drI = factor * (h - 2 * rI)
+            drI_dot_dh = factor * rI
+            dh_dot_drI = -2 * drI_dot_drI
+            dh_dot_dh = -2 * drI_dot_dh
+
+            return [[drI_dot_drI, drI_dot_dh], [dh_dot_drI, dh_dot_dh]]
 
         def terminate_burn(t, y):
             end_function = (self.grain_outer_radius - y[0]) * y[1]
@@ -494,6 +511,7 @@ def terminate_burn(t, y):
             geometry_dot,
             t_span,
             y0,
+            jac=geometry_jacobian,
             events=terminate_burn,
             atol=1e-12,
             rtol=1e-11,

tests/acceptance/test_ndrt_2020_rocket.py

@@ -27,7 +27,7 @@ def test_ndrt_2020_rocket_data_asserts_acceptance():
         "rocket_mass": (23.321 - 2.475 - 1, 0.010),
         # propulsion details
         "impulse": (4895.050, 0.033 * 4895.050),
-        "burn_time": (3.51, 0.1),
+        "burn_time": (3.45, 0.1),
         "nozzle_radius": (49.5 / 2000, 0.001),
         "throat_radius": (21.5 / 2000, 0.001),
         "grain_separation": (3 / 1000, 0.001),

tests/test_flight.py

@@ -585,6 +585,7 @@ def test_liquid_motor_flight(mock_show, calisto_liquid_modded):
         rail_length=5,
         inclination=85,
         heading=0,
+        max_time_step=0.25,
     )
 
     assert test_flight.all_info() == None
@@ -609,6 +610,7 @@ def test_hybrid_motor_flight(mock_show, calisto_hybrid_modded):
         rail_length=5,
         inclination=85,
         heading=0,
+        max_time_step=0.25,
     )
 
     assert test_flight.all_info() == None