Fixed issues with non-dim

.gitignore

@@ -16,6 +16,7 @@ TidalPy-Run*
 !TidalPy/RadialSolver/rs_solution_.hpp
 !TidalPy/Material/eos/eos_solution_.cpp
 !TidalPy/Material/eos/eos_solution_.hpp
+!TidalPy/utilities/dimensions/nondimensional_.hpp
 
 # End TidalPy
 

MANIFEST.in

@@ -20,3 +20,4 @@ include TidalPy/RadialSolver/rs_solution_.cpp
 include TidalPy/RadialSolver/rs_solution_.hpp
 include TidalPy/Material/eos/eos_solution_.cpp
 include TidalPy/Material/eos/eos_solution_.hpp
+include TidalPy/utilities/dimensions/nondimensional_.hpp

Tests/Test_Utilities/Test_Dimensions/test_nondimensional.py

@@ -5,62 +5,52 @@
 import TidalPy
 TidalPy.test_mode = True
 
-from TidalPy.constants import G
-from TidalPy.utilities.dimensions.nondimensional import non_dimensionalize_physicals, redimensionalize_physicals
+from math import isnan, isclose, sqrt
 
+from TidalPy.constants import G, PI_DBL
+from TidalPy.utilities.dimensions.nondimensional import NonDimensionalScalesClass, build_nondimensional_scales
 
-def test_non_dimensionalize_physicals():
 
-    frequency     = 1.0e-3
-    mean_radius   = 1.0e6
-    bulk_density  = 5500.
-    radius_array  = np.linspace(0., mean_radius, 10)
-    density_array = np.linspace(7000, 3500, 10)
-    bulk_array    = 100.0e9 * np.ones_like(radius_array)
-    gravity_array = np.linspace(0.1, 7.0, 10)
-    shear_array   = (50.0e9 + 20.0e3j) * np.ones(radius_array.size, dtype=np.complex128)
+frequency     = 1.0e-3
+mean_radius   = 1.0e6
+bulk_density  = 5500.
 
-    # Make copies of the arrays so we can keep the original values for comparison.
-    radius_array_out  = radius_array.copy()
-    density_array_out = density_array.copy()
-    bulk_array_out    = bulk_array.copy()
-    gravity_array_out = gravity_array.copy()
-    shear_array_out   = shear_array.copy()    
+def test_non_dimensionalize_structure():
+    """ Test that the non-dimensionalize structure initializes correctly. """
 
-    # non-dimensionalize the arrays and get some non-dim constants.
-    output_freq_nondim, G_nondim = non_dimensionalize_physicals(
-        frequency, mean_radius, bulk_density,
-        radius_array_out, density_array_out, gravity_array_out, bulk_array_out, shear_array_out)
+    test_struct = NonDimensionalScalesClass()
 
-    # Build conversions for comparison
-    second2_conversion = 1. / (np.pi * G * bulk_density)
-    second_conversion  = np.sqrt(second2_conversion)
-    length_conversion  = mean_radius
-    density_conversion = bulk_density
-    mass_conversion    = bulk_density * mean_radius**3
-    pascal_conversion  = mass_conversion / (length_conversion * second2_conversion)
-
-    # Check conversion is correct
-    assert np.allclose(radius_array / length_conversion, radius_array_out)
-    assert np.allclose(density_array / density_conversion, density_array_out)
-    assert np.allclose(gravity_array / (length_conversion / second2_conversion), gravity_array_out)
-    assert np.allclose(bulk_array / pascal_conversion, bulk_array_out)
-    assert np.allclose(shear_array / pascal_conversion, shear_array_out)
-
-    assert np.isclose(G_nondim, G / (length_conversion**3 / (mass_conversion * second2_conversion)))
-    assert np.isclose(output_freq_nondim, frequency / (1. / second_conversion))
-
-    # Check that the redimensionalize function works too.
-    # Use the output from the non-dim as inputs to the redim.
-    output_freq_redim, G_redim = redimensionalize_physicals(
-        frequency, mean_radius, bulk_density,
-        radius_array_out, density_array_out, gravity_array_out, bulk_array_out, shear_array_out)
-
-    # These should all match the original again.
-    assert np.allclose(radius_array, radius_array_out)
-    assert np.allclose(density_array, density_array_out)
-    assert np.allclose(gravity_array, gravity_array_out)
-    assert np.allclose(bulk_array, bulk_array_out)
-    assert np.allclose(shear_array, shear_array_out)
-    assert np.isclose(G_redim, G)
-    assert np.isclose(output_freq_redim, frequency)
+    # All conversions should initialize as nans
+    assert isnan(test_struct.second2_conversion)
+    assert isnan(test_struct.second_conversion)
+    assert isnan(test_struct.length_conversion)
+    assert isnan(test_struct.length3_conversion)
+    assert isnan(test_struct.density_conversion)
+    assert isnan(test_struct.mass_conversion)
+    assert isnan(test_struct.pascal_conversion)
+
+
+def test_build_nondimensional_scales():
+    """ Test building a non-dimensionalize structure with real inputs. """
+
+    # Build
+    test_struct = build_nondimensional_scales(frequency, mean_radius, bulk_density)
+
+    # Check values were set correctly
+    assert not isnan(test_struct.second2_conversion)
+    assert not isnan(test_struct.second_conversion)
+    assert not isnan(test_struct.length_conversion)
+    assert not isnan(test_struct.length3_conversion)
+    assert not isnan(test_struct.density_conversion)
+    assert not isnan(test_struct.mass_conversion)
+    assert not isnan(test_struct.pascal_conversion)
+
+    # Check that they are the correct values.
+    second2_conv = 1. / (PI_DBL * G * bulk_density)
+    assert isclose(test_struct.second2_conversion, second2_conv)
+    assert isclose(test_struct.second_conversion, second2_conv)
+    assert isclose(test_struct.length_conversion, mean_radius)
+    assert isclose(test_struct.length3_conversion, mean_radius**3)
+    assert isclose(test_struct.density_conversion, bulk_density)
+    assert isclose(test_struct.mass_conversion, bulk_density * mean_radius**3)
+    assert isclose(test_struct.pascal_conversion, (bulk_density * mean_radius**3) / (mean_radius * second2_conv))

TidalPy/Material/eos/eos_solution.pxd

@@ -7,6 +7,11 @@ cnp.import_array()
 
 from CyRK cimport CySolverResult
 
+from TidalPy.utilities.dimensions.nondimensional cimport NonDimensionalScalesCC
+
+cdef extern from "nondimensional_.hpp" nogil:
+    pass
+
 cdef extern from "eos_solution_.cpp" nogil:
 
     const size_t EOS_Y_VALUES
@@ -16,6 +21,8 @@ cdef extern from "eos_solution_.cpp" nogil:
     cdef cppclass EOSSolutionCC:
             int error_code
             int iterations
+            int nondim_status
+            int solution_nondim_status
             cpp_bool success
             cpp_bool max_iters_hit
             cpp_bool radius_array_set
@@ -35,6 +42,13 @@ cdef extern from "eos_solution_.cpp" nogil:
             double mass
             double moi
 
+            double redim_length_scale
+            double redim_gravity_scale
+            double redim_mass_scale
+            double redim_density_scale
+            double redim_moi_scale
+            double redim_pascal_scale
+
             vector[double] upper_radius_bylayer_vec
             vector[shared_ptr[CySolverResult]] cysolver_results_sptr_bylayer_vec
 
@@ -64,14 +78,13 @@ cdef extern from "eos_solution_.cpp" nogil:
                 const size_t layer_index,
                 const double radius,
                 double* y_interp_ptr)
-            void call_vectorize(
-                const size_t layer_index,
-                const double* radius_array_ptr,
-                size_t len_radius_array,
-                double* y_interp_ptr)
 
             void change_radius_array(
                 double* new_radius_ptr,
                 size_t new_radius_size)
 
             void interpolate_full_planet()
+
+            void dimensionalize_data(
+                NonDimensionalScalesCC* nondim_scales,
+                cpp_bool redimensionalize)

TidalPy/Material/eos/eos_solution_.cpp

@@ -29,7 +29,7 @@ EOSSolutionCC::EOSSolutionCC(
 
 EOSSolutionCC::~EOSSolutionCC( )
 {
-    printf("TidalPy::EOSSolutionCC Deconstructor Called.\n");
+    printf("TidalPy::EOSSolutionCC Deconstructor Called; addr = %p.\n", this);
     // Reset each shared pointer in the cysolver vector
     for (size_t i = 0; i < this->cysolver_results_sptr_bylayer_vec.size(); i++)
     {
@@ -63,28 +63,62 @@ void EOSSolutionCC::call(
         const double radius,
         double* y_interp_ptr) const
 {
-    if (layer_index < this->current_layers_saved) [[unlikely]]
+    if (layer_index < this->current_layers_saved) [[likely]]
     {
         this->cysolver_results_sptr_bylayer_vec[layer_index]->call(radius, y_interp_ptr);
-    }
-    else
-    {
-        // TODO: Better error handling
-        printf("Error! EOSSolution::call was asked to interpolate a layer that it has not saved.");
-        std::exception();
-    }
-}
 
+        if (this->nondim_status == 1)
+        {
+            // Acceleration due to Gravity
+            y_interp_ptr[0] *= this->redim_gravity_scale;
 
-void EOSSolutionCC::call_vectorize(
-        const size_t layer_index,
-        const double* radius_array_ptr,
-        size_t len_radius_array,
-        double* y_interp_ptr) const
-{
-    if (layer_index < this->current_layers_saved)
-    {
-        this->cysolver_results_sptr_bylayer_vec[layer_index]->call_vectorize(radius_array_ptr, len_radius_array, y_interp_ptr);
+            // Pressure
+            y_interp_ptr[1] *= this->redim_pascal_scale;
+
+            // These assume sphereical symmetry
+            // Total mass
+            y_interp_ptr[2] *= this->redim_mass_scale;
+
+            // Moment of inertia (r^2 multipled by dm which was found above)
+            y_interp_ptr[3] *= this->redim_moi_scale;
+
+            // Density
+            y_interp_ptr[4] *= this->redim_density_scale;
+
+            // Shear modulus (real and complex)
+            y_interp_ptr[5] *= this->redim_pascal_scale;
+            y_interp_ptr[6] *= this->redim_pascal_scale;
+
+            // Bulk modulus (real and complex)
+            y_interp_ptr[7] *= this->redim_pascal_scale;
+            y_interp_ptr[8] *= this->redim_pascal_scale;
+        }
+        else if (this->nondim_status == -1)
+        {
+            // Acceleration due to Gravity
+            y_interp_ptr[0] /= this->redim_gravity_scale;
+
+            // Pressure
+            y_interp_ptr[1] /= this->redim_pascal_scale;
+
+            // These assume sphereical symmetry
+            // Total mass
+            y_interp_ptr[2] /= this->redim_mass_scale;
+
+            // Moment of inertia (r^2 multipled by dm which was found above)
+            y_interp_ptr[3] /= this->redim_moi_scale;
+
+            // Density
+            y_interp_ptr[4] /= this->redim_density_scale;
+
+            // Shear modulus (real and complex)
+            y_interp_ptr[5] /= this->redim_pascal_scale;
+            y_interp_ptr[6] /= this->redim_pascal_scale;
+
+            // Bulk modulus (real and complex)
+            y_interp_ptr[7] /= this->redim_pascal_scale;
+            y_interp_ptr[8] /= this->redim_pascal_scale;
+        }
     }
     else
     {
@@ -147,6 +181,8 @@ void EOSSolutionCC::change_radius_array(
 void EOSSolutionCC::interpolate_full_planet()
 {
     printf("TidalPy::EOSSolutionCC.interpolate_full_planet called.\n");
+    this->solution_nondim_status = this->nondim_status;
+
     if (this->current_layers_saved == 0)
     {
         // No layers have been saved, we can't perform the interpolation.
@@ -208,4 +244,117 @@ void EOSSolutionCC::interpolate_full_planet()
 
     // Finished
     this->other_vecs_set = true;
+}
+
+#include <cstdio>
+
+void EOSSolutionCC::dimensionalize_data(NonDimensionalScalesCC* nondim_scales, bool redimensionalize)
+{   
+
+    printf("EOSSolutionCC::dimensionalize_data Called.\n");
+    printf("nondim scales (ptr = %p):\n", nondim_scales);
+    printf("\t length = %e\n", nondim_scales->length_conversion);
+    printf("\t leng3  = %e\n", nondim_scales->length3_conversion);
+    printf("\t densit = %e\n", nondim_scales->density_conversion);
+    printf("\t pascal = %e\n", nondim_scales->pascal_conversion);
+    printf("\t sec    = %e\n", nondim_scales->second_conversion);
+    printf("\t sec2   = %e\n", nondim_scales->second2_conversion);
+    printf("\t mass   = %e\n", nondim_scales->mass_conversion);
+
+    // Save scalers
+    this->redim_length_scale  = nondim_scales->length_conversion;
+    this->redim_gravity_scale = nondim_scales->length_conversion / nondim_scales->second2_conversion;
+    this->redim_mass_scale    = nondim_scales->mass_conversion;
+    this->redim_density_scale = nondim_scales->density_conversion;
+    this->redim_moi_scale     = nondim_scales->mass_conversion * nondim_scales->length_conversion * nondim_scales->length_conversion;
+    this->redim_pascal_scale  = nondim_scales->pascal_conversion;
+
+    // Figure out how to set flags used during eos solution calls
+    if (this->solution_nondim_status == 0)
+    {
+        // EOS was not explicitly non-dim'd or re-dim'd when solution was found. 
+        if (redimensionalize)
+        {
+            // User is now "redimensionalizing." Thus we assume that the solution was non-dim'd
+            this->nondim_status = 1;
+        }
+        else
+        {
+            // User is now "dimensionalizing." Thus we assume that the solution was dim'd
+            this->nondim_status = -1;
+        }
+    }
+    else
+    {
+        // TODO: Deal with this case! For now push the problem to the user when they try to call...
+    }
+
+    // Update other constants
+    if (redimensionalize)
+    {
+        this->pressure_error *= this->redim_pascal_scale;
+    }
+    else
+    {
+        this->pressure_error /= this->redim_pascal_scale;
+    }
+
+    // Update layer data
+    for (size_t layer_i = 0; layer_i < this->num_layers; layer_i++)
+    {
+        if (redimensionalize)
+        {
+            this->upper_radius_bylayer_vec[layer_i] *= this->redim_length_scale;
+        }
+        else
+        {
+            this->upper_radius_bylayer_vec[layer_i] /= this->redim_length_scale;
+        }
+    }
+
+    if (this->other_vecs_set)
+    {
+        // Work through arrays
+        for (size_t slice_i = 0; slice_i < this->radius_array_size; slice_i++)
+        {
+            if (redimensionalize)
+            {
+                this->radius_array_vec[slice_i]   *= this->redim_length_scale;
+                this->gravity_array_vec[slice_i]  *= this->redim_gravity_scale;
+                this->pressure_array_vec[slice_i] *= this->redim_pascal_scale;
+                this->mass_array_vec[slice_i]     *= this->redim_mass_scale;
+                this->moi_array_vec[slice_i]      *= this->redim_moi_scale;
+                this->density_array_vec[slice_i]  *= this->redim_density_scale;
+
+                // These complex arrays are stored as double arrays with twice the length (Cython and C++ don't play nicely with complex across all systems)
+                this->complex_shear_array_vec[2 * slice_i]     *= this->redim_pascal_scale;
+                this->complex_shear_array_vec[2 * slice_i + 1] *= this->redim_pascal_scale;
+                this->complex_bulk_array_vec[2 * slice_i]      *= this->redim_pascal_scale;
+                this->complex_bulk_array_vec[2 * slice_i + 1]  *= this->redim_pascal_scale;
+            }
+            else
+            {
+                this->radius_array_vec[slice_i]   /= this->redim_length_scale;
+                this->gravity_array_vec[slice_i]  /= this->redim_gravity_scale;
+                this->pressure_array_vec[slice_i] /= this->redim_pascal_scale;
+                this->mass_array_vec[slice_i]     /= this->redim_mass_scale;
+                this->moi_array_vec[slice_i]      /= this->redim_moi_scale;
+                this->density_array_vec[slice_i]  /= this->redim_density_scale;
+
+                // These complex arrays are stored as double arrays with twice the length (Cython and C++ don't play nicely with complex across all systems)
+                this->complex_shear_array_vec[2 * slice_i]     /= this->redim_pascal_scale;
+                this->complex_shear_array_vec[2 * slice_i + 1] /= this->redim_pascal_scale;
+                this->complex_bulk_array_vec[2 * slice_i]      /= this->redim_pascal_scale;
+                this->complex_bulk_array_vec[2 * slice_i + 1]  /= this->redim_pascal_scale;
+            }
+        }
+
+        // Update global constants
+        this->radius           = this->radius_array_vec[this->radius_array_size - 1];
+        this->surface_gravity  = this->gravity_array_vec[this->radius_array_size - 1];
+        this->surface_pressure = this->pressure_array_vec[this->radius_array_size - 1];
+        this->mass             = this->mass_array_vec[this->radius_array_size - 1];
+        this->moi              = this->moi_array_vec[this->radius_array_size - 1];
+        this->central_pressure = this->pressure_array_vec[0];
+    }
 }