jrenaud90 / tidalpy Goto Github PK
View Code? Open in Web Editor NEWSoftware suite to perform solid-body tidal dissipation calculations for rocky and icy worlds
License: Other
Software suite to perform solid-body tidal dissipation calculations for rocky and icy worlds
License: Other
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.
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.
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.
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.
Describe the bug
This is a thread to track improvements / issues related to integration stability when using dynamic liquid layers in RadialSolver.
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.
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?
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.
Implement Tests for Multi-Layer Code
Requires #1 to be completed first.
As of version 0.1.0 the dill/pickle interface is only half-implemented.
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).
Please read the style.md before making new doc strings.
Docstrings should be concise
Sub-Packages
Other
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.
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:
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.
Ensure all dependencies are in setup.py and all information is up to date.
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 TidalPy to PyPI
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:
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
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").
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.
The beginnings of this were implemented in version 0.2.1, but the sphinx documentation is sparse compared to the older markdown files.
Commit Workflow
Make tests for any improvements made by #3
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.
As is used by Henning & Hurford 2014
Ensure that Multi-Layer code is working efficiently with orbital evolution (time studies).
May end up using Cython to speed things up.
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.
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".
Issue with cached numba functions vanishing when called from multiple scripts at once.
Removing cache from some of the most affected functions for now.
So that one call can change multiple planet orbital parameters at once.
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
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.
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 a function to quickly build common systems (instantiate and return the orbit class for each system):
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 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
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.