jrenaud90/tidalpy

Software suite to perform solid-body tidal dissipation calculations for rocky and icy worlds

tidalpy's Issues

Add solidus and liquidus as a function of pressure

This will be easy to do for silicates (see the equations 21 and 22 in Shoji & Kurita 2014).

Run these equations at layer initialization (TidalPy does not allow for a changing pressure at the moment, so these do not need to be recalculated after the planet is built).

Now phase change for other materials, namely ice, will be a much bigger challenge and deserves its own issue.

FLoating Point Issues in Takeuchi Starting Conditions

An earlier version of TidalPy used the full definition of Takeuchi & Saito (1972)'s starting conditions in the radial function.
However, testing revealed that there were issues matching the limiting values (as z -> 0) since the full solution has terms that are like: (Very Big) / (Very Big) ~ near 1. Errors would compound leaving the final answer larger or smaller than 1.

Current Solution:
Have reverted back to using the limiting version of these equations.

Add ability to calculate effective temperature based on rheology and strength

Currently, TidalPy supports calculations of the form:
viscosity = viscosity(T)
shear = shear(T)

Where the functional forms of viscosity and shear modulus are determined by the rheology.

It would be nice to invert these equations such that:
T = T(viscosity, shear)

There will be some degeneracy at higher temperatures when many temperatures can produce the same viscosity and shear (that being those of liquids). But, it may still be useful in a restricted temperature band, or in the full temperature domain but with a user warning if the effective viscosity is too low for an accurate temperature calculation.

See method set_strength in the Layer class for more information.

RadialSolver solution gives unallocated memory for Love number results

Users have reported that successful runs of RadialSolver.radial_solver give unallocated values when accessing the solutions.k, .h, .l etc.

# Problem
def run():
    solution = radial_solver(*input)
    k = solution.k
    return k

k2 = run()
print(k2)  # Random value. Usually something super small like ~1e-300

# Workaround
import numpy as np
def run():
    solution = radial_solver(*input)
    k = solution.k[0]
    # or 
    k = np.copy(solution.k)
    return k
    # or return the solution instance: `return solution`

k2 = run()
print(k2)  # Correct value!

This has occurred with, at least, TidalPy v0.5.4 but as discussed below, it is not a TidalPy problem per-say.

The problem is that the RadialSolverSolution class performs a heap allocation to memory in order to store the Love number results. It is responsible for freeing this memory when the class instance is destroyed. This particular problem comes up when the solution is destroyed too early.

To Reproduce
The common situation where this problem arises is when an outer python function calls radial solver and gets a solution but does not return the solution. This can be reproduced with pure Cython and no use of TidalPy:

To reproduce make a cython class and functions like...

## Cython
from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free

import numpy as np
cimport numpy as np

cdef class Solution:
    
    cdef size_t n
    cdef double complex* result_ptr
    cdef double complex[::1] result_view
    
    def __init__(self, size_t n):
        
        self.n = n
        
        self.result_ptr = <double complex*> PyMem_Malloc(self.n * sizeof(double complex))
        if not (self.result_ptr is NULL):
            self.result_view = <double complex[:self.n]>self.result_ptr
    
    @property
    def result(self):
        return np.ascontiguousarray(self.result_view, dtype=np.complex128)
    
    def __dealloc__(self):
        if not (self.result_ptr is NULL):
            PyMem_Free(self.result_ptr)
    

cdef Solution cy_run():
    
    cdef size_t n = 3
    cdef size_t i
    cdef Solution solution = Solution(n)
    
    for i in range(n):
        solution.result_ptr[i] = i * (10. + 1.j)
    
    return solution

def run():
    return cy_run()

Call this function from another python wrapping function...

## Python
import numpy as np
def run_outer():
    
    solution = run()
    k = solution.result
    
    return k

print(run_outer())

You will see unallocated values are printed.

Workaround
We can get the correct values by grabbing a specific value of the array, or by making a copy of the array, then the correct values are printed...

import numpy as np
def run_outer():
    
    result = run()
    k = np.copy(result.result)
    # or
    k = result.result[0]
    
    return k

print(run_outer())

Alternatively, and probably the better solution, is to return the solution instance itself so that it is no longer tied to the scope of run_outer

Future Fix
It is hard to know the solution to this problem that is not one of the workarounds mentioned above. We want RadialSolverSolution to own the heap memory otherwise its owner may never free it properly. The downside is a situation like this where it is deallocated too early.

Certain open to suggestions but for now this probably just requires updates to the documentation.

RadialSolver: Dynamic Liquid Layer Stability Issues

Describe the bug
This is a thread to track improvements / issues related to integration stability when using dynamic liquid layers in RadialSolver.

• Causes
  • Stability problems appear when a dynamic liquid layer is present in a world with a moderately large forcing period.
  • The value of the forcing period depends on the density of the liquid. A low density liquid tends to be more stable at higher forcing periods than a high density layer.

Note that "forcing period" refers to the tidal mode that is being calculated. A planet with a low orbital period may still require larger forcing period modes depending on its eccentricity, obliquity, and rotation rate.

Should Stefan number be changed to a function of melt fraction

Right now it is set as a configurable constant.

Looking at Equation 3 in Shoji + Kurita (2014, ApJ vol 789) they state:
St = derivative_vol_melt_fraction * heat_fusion * volume_viscoelastic / (specific_heat * volume_layer)

Is this an accurate model? Where is this derived?

Add a check so that andrade alpha stays between 0.0 and 1.0

The custom njit factorial function does not work outside this range. Physically that is fine since we do not expect alpha to be outside this range. However, the current code allows it to fall outside this range.

Multi-Layer Tidal Calculator Tests

Implement Tests for Multi-Layer Code

Requires #1 to be completed first.

Ensure that dill loading and saving is working as expected

As of version 0.1.0 the dill/pickle interface is only half-implemented.

Implement a more realistic lower layer temperature for OOP

Currently, the base temperature of a layer is equal to its average temperature (layer.temperature_lower = layer.temperature).

In reality, this temperature will be higher. Generally, this is found via an assumed viscosity drop across the layer (See Hussmann & Spohn 2004).

• Look into seeing if this is a decent enough model to use.
• Determine if the viscosity drop is going to be robust enough to handle different viscosity equations (Reference vs. Arrhenius).
• Implement it into the layer class.

Add Docstrings

Please read the style.md before making new doc strings.

Docstrings should be concise

Sub-Packages

• burnman_interface
• cooling
• dynamics
• integration
• radiogenics
• rheology
• stellar
• structures
• tides
• tools
• utilities

Other

• Loose /TidalPy/ Files
• Cookbooks and Examples
• Documentation Folder and/or Wiki
  • OOP vs. Functional
  • How-to add new rheologies
  • How-to add new planets
  • How-to Burnman Basics

Make sure all function calls to njit'd functions which have default arguments have all of their arguments spelled out

There is a problem with the current version of numba (see bug report here) where if a function has default args but they are not used in a function call then the call is significantly slower.

The current work around is to provide all default args whenever possible.

Performance Improvements

The current implementation of the functions in TidalPy relies on Numba's njit method which complies code at run-time. The resulting code is quite fast in many circumstances, especially when used alongside np.ndarrays. However, each time these functions are called they must be recompiled which can be quite slow. This tends not to be a problem for time studies which only need to compile once, but for static exploration (most of the Jupyter Notebooks) the compile times can be quite annoying.

In TidalPy v0.2.0, cacheable njit functions were introduced. These are only compiled once per unique call (where unique refers to the signature of the inputs, e.g., all floats, all arrays, some mix, etc.). They are then loaded from disk for subsequent calls. This can greatly increase the initial performance but does come at a small hit to run-time performance. There are also other issues with caching functions as mentioned in the "gotchas" documentation. Lastly, many of these functions should not change and do not need all of Python's overhead. There is no obvious reason they can not be pre-compiled using something like Cython.

Todo:

• Look into Cython vs. raw C vs. Fortran (f2py) as options to write some lower level functions. What are the performance vs. compatibility benefits? What are the benefits vs. Numba (which is getting faster each day)?
• If one of these approaches seems like a better option: how will they integrate in the OOP and functional scheme of TidalPy?
• List the functions that need this attention:
  • Eccentricity and Inclination functions
  • Mode collapse functions

Inconsistent definition of initial conditions for multilayer tidal shooting method

There is a mismatch in definitions between TS72 Eq. 99 (Takeuchi, H., and M. Saito (1972), Seismic surface waves, Methods Comput. Phys., 11, 217–295.) and KMN15 Eq. B13 (Kamata+ (2015; JGR-P; DOI: 10.1002/2015JE004821)) for the initial guess used in the shooting method for solid layers. TS72 has a "minus/plus" for k2 whereas KMN has a "plus/minus" basically flipping solution 1 and 2.

• Check if this makes a difference to begin with (try flipping solution 1 and solution 2)
• If it does then try to determine which is correct
• Update code and TODOs at:
  • Around line 90 and 314 in tides.multilayer.numerical_int.initial_solution_dynamic.py
  • Around line 90 and 216 in tides.multilayer.numerical_int.initial_solution_static.py

Finalize setup.py

Ensure all dependencies are in setup.py and all information is up to date.

cibuildwheel fails for x86 and musl-linux

Having issues with cibuildwheel building wheels for x86 architectures (on windows and linux) as well as musl-linux. Have disabled these for now in pyproject.toml.

Add to PyPI

Add TidalPy to PyPI

Add method to world class to allow for more elegant config changes

Currently, to change a configuration on a planet you need to A). know a potentially complex series of configuration keys and B). make a separate call to the planet's reinit(). For example to change the rheology of the earth's lower mantle:
earth.config['layers']['Lower_Mantle']['rheology']['compliance_model'] = 'Andrade'
earth.reinit()

All of those keys are case sensitive and any mistake will cause a crash due to KeyError. This can be confusing for someone new to TidalPy.

It would be nice to have a method that:

• Reduces the number of keys needed (by making some assumptions/guess about what the user really wants to change.
• Cleans up the remaining keys so that if the user did/didn't capitalize something it doesn't cause a crash.
• Allows the user to change multiple configs with one method call.
• Automatically calls reinit once all the configs have been changed (user togglable).

Move viscosity (at least pre-melt) calculations to the rheology package

Right now it is in the thermal package but uses defaults from rheology - and it fits in more with rheology.

This will break older versions so if it is done it must be done before 1.0

Recrusive formula for z-function in TS/KMN starting conditions does not match actual

The recrusision formula shown in TS72 and KMN15 (See TS72 Eq. 97) does not reproduce the full version of this function (TS72 Eq. 96). I am finding it to be off by about 30% at l = 2 and then the error improves to 10% at l=10. Perhaps thre is a problem with the implementation.

For this reason, for now, TidalPy does not use this recrusive formula (see cf_z_calc in "RadialSolver.starting.common.pyx").

Inclination is not shared between bodies - Ensure code takes this into account

Dual dissipating bodies will have their own inclinations (I_host and I_target). Right now it is treated more like eccentricity where it is shared between them.

This will require changing the orbit/mode/ package, as well as some current research scripts and test cases.

Use Sphinx to build documentation

The beginnings of this were implemented in version 0.2.1, but the sphinx documentation is sparse compared to the older markdown files.

• Rewrite old markdown files into new rst files that are used by sphinx (Documentation/source/index.rst).
• Look into how to host the html documentation files on github.

Commit Workflow

• Run sphinx-apidoc to make rst files for all TidalPy modules
• Run make html (make.bat html for windows) to build the html files (these make files should be in the Documentation directory)

Tests for Multi-Layer Orbital Coupling

Make tests for any improvements made by #3

setup.py file missing?

Title explains it. I've both cloned the repository and downloaded the repository as a zip file and I can't find the setup.py file. I've searched in every directory within TidalPy and made sure setuptools was installed. I think I might be missing something obvious, but until it's revealed to me I'm still very unsure where it might be.

I created a new conda environment to run this all under. A summary of commands run thus far are:

conda create -n tidalpy python=3.9
conda activate tidalpy
git clone https://github.com/jrenaud90/TidalPy.git
cd TidalPy

I feel embarrassed opening up this issue because it's so simple but any help is appreciated.

Implement Sabadini & Vermeersen Multi-Layer Tidal Love Calculator

As is used by Henning & Hurford 2014

Couple Multi-Layer to Orbital Mechanics

Ensure that Multi-Layer code is working efficiently with orbital evolution (time studies).

May end up using Cython to speed things up.

TidalPy.RadialSolver.solver missing after crash

I'll try to detail the conditions that led up to this crash as best I can.

I was testing radial_solver with much finer layers and more property variations than previously. I'm currently using Python 3.9 and Spyder version 5.5.1 through conda.

After rebooting Spyder, however, I'm met with this error:

runfile('/Users/nlwagner25/Desktop/PYTHON/TidalPy/testing_AK_model_elastic.py', wdir='/Users/nlwagner25/Desktop/PYTHON/TidalPy')
2024-02-16 14:33:05(+00:00:28::286448) - INFO     : TidalPy initializing...
2024-02-16 14:33:05(+00:00:28::287662) - INFO     : TidalPy initialization complete.
Traceback (most recent call last):

  File ~/opt/anaconda3/envs/tidalpy/lib/python3.9/site-packages/spyder_kernels/py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File ~/Desktop/PYTHON/TidalPy/testing_AK_model_elastic.py:15
    from TidalPy.RadialSolver import radial_solver

  File ~/Desktop/PYTHON/TidalPy/TidalPy/RadialSolver/__init__.py:1
    from TidalPy.RadialSolver.solver import radial_solver

ModuleNotFoundError: No module named 'TidalPy.RadialSolver.solver'

I tried reinstalling TidalPy through Pip, but no luck. I'm completely deleting and recreating my conda environment and will update with a comment once that finishes.

Question of starting codition order TS72 vs. KMN15

Kamata et al. (2015) has the two starting solutions for RadialSolver k^2 +/- (see their Eq. B13) implying that solution 1 uses the + version and solution 2 uses the minus.

Takeuchi & Saito (1972) have similar starting solution definitions but order them as -/+ (see their Eq. 99).

TidalPy has the option to use both starting solutions. For Kamata we use +/- and for Takeuchi we use -/+. See the respective equations in "RadialSolver.starting".

Underlying Object has Vanished

Issue with cached numba functions vanishing when called from multiple scripts at once.

Removing cache from some of the most affected functions for now.

Allow the OrbitBase class setters to accept lists of planet refs

So that one call can change multiple planet orbital parameters at once.

Add Time Lag model for Gas Giant Dissipation

Implement a simplified global rheology for gas giant dissipation.

Based on Hussmann & Spohn 2004:
-Im[k2] = k_J (t) / Q
k_J (t) = k_J (t_0) * (1 - n_io(t) / spin_J) / (1 - n_io(t_0) / spin_J)

n_io = Io's mean orbital motion
spin_J = Jupiter's spin rate (H&S assumed constant over time)
k_J(t_0) = Jupiter's modern k2 = .379
Q = Jupiter's fixed Q. H&S let this be a free parameter >= 6.6e4

• Is this still an accepted estimate for gas giant dissipation to first order?
• Can this be easily applied to other gas giants (Saturn) that have a dominant dissipating satellite (Titan?? Enceladus??)

Degree-1 Love Number Issues

Opening this bug report to document exploring why degree-1 Love numbers are often too large. Will add reports on test cases to see if this is an integration issue or not.

Turn "TODO"'s and "OPT"'s into github issues

Ctrl-f "TODO" and "OPT" throughout the project and make a new issue for each item. It is okay if you don't know what exactly to put (labels, comments, etc.). Just try to make educated guesses and use the context around the comment.

"OPT" refers to areas of the code that could be optimized in the future.

Add build_system function

Add a function to quickly build common systems (instantiate and return the orbit class for each system):

• Jovian system
• Earth-Moon system
• Pluto-Charon system
• TRAPPIST-1 system
• Generic super-Earth around m-star

Test what layer scales work well for RadialSolver

The RadialSolver method has the option of scaling relative and/or absolute tolerances based on the layer type.
See, for example, here.
The scales that are currently in place are just best guesses. Some thorough benchmarking and testing is needed to provide better scales.

Have dill look at the planet's tidalpy version

Have the dill loader check the TidalPy version that a planet was saved under and either raise an error or a warning if a user tries to load an incompatible version