search

jrenaud90 / TidalPy

Software suite to perform solid-body tidal dissipation calculations for rocky and icy worlds
Other
15 stars 3 forks source link

Degree-1 Love Number Issues #53

Open nlwagner opened 7 months ago

nlwagner commented 7 months ago

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.

jrenaud90 commented 7 months ago

Awesome, thanks for working on this.

Here is where the boundary condition is defined: https://github.com/jrenaud90/TidalPy/blob/edab615f5ade69b8923536615d7cb9942ac5dce5/TidalPy/RadialSolver/solver.pyx#L306

And here is where it is applied: https://github.com/jrenaud90/TidalPy/blob/edab615f5ade69b8923536615d7cb9942ac5dce5/TidalPy/RadialSolver/boundaries/boundaries.pyx#L47

You could make modifications to the first link if it is as simple as a factor off. The 3 BCs that are defined there are applied to y1, y3, and y6. If you need to apply a BC to a different y (maybe y5 as in LoadDef's table 5.2) then you will need to hack the code in the 2nd link. Notice the [0 * max_num_y + 1], [0 * max_num_y + 3], [0 * max_num_y + 5], etc? Those are pulling out the values of y2, y4, y6 respectively. So applying a BC to the other ys would require changing that index.

Note any changes you make will need to be recompiled by doing a python -m pip install -v . on your local, modified version of TidalPy.

jrenaud90 commented 7 months ago

Putting a note here to remind us that if the trick of subtracting off k_1 from the other love numbers is indeed the correct approach that it should be easy to implement as a post-process. So that is a ToDo if we go that route.

nlwagner commented 7 months ago

truthfully through reading LoadDef's manual more and testing I think the k_1 trick is enough to get more accurate h1s. However! Ive found something kind of fascinating. In the attached plot I have h1 as a function of forcing frequency. (note these are all Load Love numbers and elastic) Figure_1

I figured this sort of insight might be useful for you. The slight decrease in h1 after about 100 hours is strange to me though, do you have any insight? k2 and h2 don't have this decrease around the same value: Figure_2

nlwagner commented 7 months ago

Note: I was getting stable h1s with a simpler interior model. Increasing complexity makes h1 want to misbehave. I'll work on quantifying this more in the future.

nlwagner commented 7 months ago

Here's a more complete plot before I sign off for the day: Figure_4

jrenaud90 commented 7 months ago

That is really interesting. Thanks for sharing - would you be willing to share your inputs to radial solver? I can poke around to see if I can learn anything more.

My one guess is that you are hitting some sort of resonance with the rheology. A good check is if you change the viscosity (not sure if you are using a constant model or some $\eta(r)$ function. But either way if you take it and bump it up or down $eta(r)_n = 1000 \times \eta(r)_o$ (or divide by 1000). And see if that changes anything? I am surprised that would change l=1 and not l>1 though...

nlwagner commented 7 months ago

Oops sorry just saw this. Definitely. Once I push my workspace to my forked repository I'll let you know.

Additionally, I've noticed singularities occuring at really long periods (>1000 hours) that appear/disappear with varying integration tolerances. Seen below. Figure_8

If we tighten the tolerances we see: Figure_9

So probably no surprise we see numerical issues at really long periods I'm just trying to double check degree 1 isn't actually dependent on frequency in an elastic case :)

nlwagner commented 7 months ago

Within my forked repository "mars", the file that generated the above plots is deg1deg2_stability_complexmodel.py

Note I had to add 1500 slices for each layer so it takes a while to run.