Merge pull request #51 from jrenaud90/ver-0.5.2

Ver 0.5.2

CHANGES.md

@@ -4,6 +4,16 @@
 * Other Major Changes
   * Remove support for the older non-cythonized `radial_solver` module.
 
+### Version 0.5.2 (2024-02-22)
+
+Documentation
+* Improved RadialSolver documentation regarding higher degree-l.
+* Added info about issues that can arise from using non C-continguous arrays in cythonized functions.
+
+Fixes:
+* Added error message to `RadialSolver.radial_solver` if length of provided assumption tuples is not the same.
+* Fixed issue where non C-continguous arrays were allowed in cythonized functions that required them.
+
 ### Version 0.5.1 (2024-02-14)
 * Removed Python 3.8 support due to issues with building SciPy.
 

Documentation/Rheology.md

@@ -50,6 +50,9 @@ print(complex_shear)
 TidalPy provides helper methods to efficiently parse over arrays. These functions use multithreading when possible to 
 quickly calculate results over large arrays.
 
+Note that all arrays must be [C-contiguous](https://stackoverflow.com/questions/26998223/what-is-the-difference-between-contiguous-and-non-contiguous-arrays).
+If you suspect that an array may not be C-contiguous you can use the numpy function `arr = np.ascontiguousarray(arr)` to ensure that they are before being passed to the rheology methods.
+
 ```python
 from TidalPy.models import Andrade
 

Documentation/Using RadialSolver.md

@@ -14,6 +14,9 @@ To learn more about this method please review the literature cited in the refere
 ## `TidalPy.RadialSolver.radial_solver` Function
 The `radial_solver` function, contained in the `TidalPy.RadialSolver` module can be used with the following arguments.
 
+Note that all arrays must be [C-contiguous](https://stackoverflow.com/questions/26998223/what-is-the-difference-between-contiguous-and-non-contiguous-arrays).
+If you suspect that an array may not be C-contiguous you can use the numpy function `arr = np.ascontiguousarray(arr)` to ensure that they are before being passed to the rheology methods.
+
 ```python
 from TidalPy.RadialSolver import radial_solver
 radial_solver_solution = radial_solver(
@@ -66,7 +69,9 @@ radial_solver_solution = radial_solver(
     # Harmonic degree used. 
     # Note that the stability of the solution gets worse with higher-l. You may need to increase the relative tolerance
     # or use simpler layer assumptions to get a successful solution. 
-    # It particularly starts to break down for l > 10.
+    # It particularly starts to break down for l > 10 but stable solutions exist up to l=40.
+    # For higher degrees, it is recommended to start integration higher in the planet. I.,e. above the
+    #  core-mantle-boundary (if applicable).
 
     solve_for = None,
     # What to solve for (type: tuple[str, ...])

README.md

@@ -11,7 +11,7 @@
 
 ---
 
-<a href="https://github.com/jrenaud90/TidalPy/releases"><img src="https://img.shields.io/badge/TidalPy-0.5.1 Alpha-orange" alt="TidalPy Version 0.5.1 Alpha" /></a>
+<a href="https://github.com/jrenaud90/TidalPy/releases"><img src="https://img.shields.io/badge/TidalPy-0.5.2 Alpha-orange" alt="TidalPy Version 0.5.1 Alpha" /></a>
 
 **Tidal Dynamics and Thermal-Orbital Evolution Software Suite Implemented in Cython and Python**
 

TidalPy/RadialSolver/derivatives/odes.pyx

@@ -752,7 +752,7 @@ cdef RadialSolverBase cf_build_solver(
     cdef RadialSolverBase solver
 
     # Convert the y0 pointer to a memoryview in order to work with CyRK's CySolver __init__
-    cdef double[:] y0_view = <double[:num_ys]> y0_ptr
+    cdef double[::1] y0_view = <double[:num_ys]> y0_ptr
 
     if (layer_type == 0):
         # Solid layer

TidalPy/RadialSolver/love.pyx

@@ -30,18 +30,18 @@ cdef void find_love_cf(
 
 
 def find_love(
-    double complex[:] complex_love_numbers_view,
-    double complex[:] surface_solutions_view,
+    double complex[::1] complex_love_numbers_view,
+    double complex[::1] surface_solutions_view,
     double surface_gravity
     ):
     """
     Find the complex Love and Shida numbers given the surface radial solutions for a planet.
 
     Parameters
     ----------
-    complex_love_numbers_view : double complex[:], array, output
+    complex_love_numbers_view : double complex[::1], array, output
         Array to store complex Love numbers. There must be space for 3 double complex numbers.
-    surface_solutions_view : double complex[:], array, input
+    surface_solutions_view : double complex[::1], array, input
         Array of radial solutions (y_i) values at the surface of a planet.
     surface_gravity : double, input
         Acceleration due to gravity at the planet's surface [m s-2].

TidalPy/RadialSolver/solver.pyx

@@ -1045,11 +1045,11 @@ cdef RadialSolverSolution cf_radial_solver(
 
 
 def radial_solver(
-        double[:] radius_array,
-        double[:] density_array,
-        double[:] gravity_array,
-        double[:] bulk_modulus_array,
-        double complex[:] complex_shear_modulus_array,
+        double[::1] radius_array,
+        double[::1] density_array,
+        double[::1] gravity_array,
+        double[::1] bulk_modulus_array,
+        double complex[::1] complex_shear_modulus_array,
         double frequency,
         double planet_bulk_density,
         tuple layer_types,
@@ -1184,6 +1184,17 @@ def radial_solver(
     # Unpack inefficient user-provided tuples into bool arrays and pass by pointer
     cdef size_t num_layers
     num_layers = len(layer_types)
+
+    # Check that number of assumptions match.
+    if len(is_static_by_layer) != num_layers:
+        raise AttributeError('Number of `is_static_by_layer` must match number of `layer_types`.')
+    if len(is_incompressible_by_layer) != num_layers:
+        raise AttributeError('Number of `is_incompressible_by_layer` must match number of `layer_types`.')
+    if len(upper_radius_by_layer) != num_layers:
+        raise AttributeError('Number of `upper_radius_by_layer` must match number of `layer_types`.')
+
+    # Build array of assumptions
+    # OPT: Perhaps set a maximum number of layers then we can put these on the stack rather than heap allocating them.
     cdef int* layer_assumptions_ptr = <int *> allocate_mem(
         3 * num_layers * sizeof(int),
         'layer_assumptions_ptr (radial_solver; init)'

TidalPy/rheology/base.pyx

@@ -93,10 +93,10 @@ cdef class RheologyModelBase(TidalPyBaseExtensionClass):
 
     def vectorize_frequency(
             self,
-            double[:] frequency_view,
+            double[::1] frequency_view,
             double modulus,
             double viscosity,
-            double complex[:] output_view,
+            double complex[::1] output_view,
             ):
 
         cdef Py_ssize_t n, n2
@@ -111,9 +111,9 @@ cdef class RheologyModelBase(TidalPyBaseExtensionClass):
     def vectorize_modulus_viscosity(
             self,
             double frequency,
-            double[:] modulus_view,
-            double[:] viscosity_view,
-            double complex[:] output_view,
+            double[::1] modulus_view,
+            double[::1] viscosity_view,
+            double complex[::1] output_view,
             ):
 
         cdef Py_ssize_t n, n2, n3

TidalPy/utilities/dimensions/nondimensional.pyx

@@ -151,11 +151,11 @@ def non_dimensionalize_physicals(
         double frequency,
         double mean_radius,
         double bulk_density,
-        double[:] radius_array_view,
-        double[:] density_array_view,
-        double[:] gravity_array_view,
-        double[:] bulk_array_view,
-        double_numeric[:] shear_array_view,
+        double[::1] radius_array_view,
+        double[::1] density_array_view,
+        double[::1] gravity_array_view,
+        double[::1] bulk_array_view,
+        double_numeric[::1] shear_array_view,
         ):
 
     cdef size_t num_radius = radius_array_view.size
@@ -174,11 +174,11 @@ def redimensionalize_physicals(
         double frequency,
         double mean_radius,
         double bulk_density,
-        double[:] radius_array_view,
-        double[:] density_array_view,
-        double[:] gravity_array_view,
-        double[:] bulk_array_view,
-        double_numeric[:] shear_array_view,
+        double[::1] radius_array_view,
+        double[::1] density_array_view,
+        double[::1] gravity_array_view,
+        double[::1] bulk_array_view,
+        double_numeric[::1] shear_array_view,
         ):
 
     cdef size_t num_radius = radius_array_view.size
@@ -201,7 +201,7 @@ def redimensionalize_radial_functions(
 
     Parameters
     ----------
-    radial_function_ptr : complex128*
+    radial_function_view : double complex[:, ::1]
         Non-dimensionalized radial solutions as a function of radius.
     mean_radius : float64
         Mean radius of the planet, used in scaling [m]

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name='TidalPy'
-version = '0.5.1'
+version = '0.5.2'
 description='Tidal Dynamics and Thermal-Orbital Evolution Software Suite Implemented in Cython and Python'
 authors= [ 
     {name = 'Joe P. Renaud', email = 'TidalPy@gmail.com'}