Md/stabilize solid1d

[MarijnJ0]

Description

This PR introduces a new ODE solver for the solid1d module based on the relaxation method. The new solver is much more stable than the original shooting method based solver. This allows for the inclusion of inertia effects. Additionally, the new method has been added to the documentation.

Besides this, some small issues have been addressed regarding the treatment of the core in the solid tides calculations.

The get_k_range function has also been updated, it now grabs the relevant ranges from a table that is generated by the user before starting Obliqua. A notebook is also provided so users can generate their own table based on specific requirements.

The fluid1d model has been updated, heat is now distributed radially after the global heating has been calculated.

Closes #26

I have to update the tests, when I do this I will add the relevant logs to the discussion below.

Validation of changes

The new solver obtains a solution for the y functions, to confirm the overall stability I have plotted for many cases these solutions. Below is one such plot

The oscilating bahavior persists after increasing the resolution, suggesting that it represents the physical medium response at some intermediate viscosity. Overall, all cases are internally consistent, and produce overall reasonable heating rates that are in in agreement with the shooting method results in the pure solid limit.

Tested on linux, ubuntu, with Julia Version 1.11.6.

Checklist

• I have followed the contributing guidelines
• My code follows the style guidelines of this project
• I have performed a self-review of my code
• My changes generate no new warnings or errors
• I have checked that the tests still pass on my computer
• I have updated the docs, as appropriate
• I have added tests for these changes, as appropriate
• I have checked that all dependencies have been updated, as required

Relevant people

@nichollsh @timlichtenberg

[Copilot]

Pull request overview

This PR adds a new, more stable solid-tides solver option for solid1d based on a relaxation method, expands configuration/documentation to support it, and updates related solid/fluid modeling utilities (core mass handling, Hansen k-range lookup, fluid1d radial heating distribution).

Changes:

• Added solid1d_relax module and integrated it into Obliqua via a new module_solid = "solid1d-relax" option with new config parameters (dr_min, dr_max, core).
• Updated solid gravity calculations to include explicit core mass (m_core) and adjusted Hansen get_k_range to use a precomputed NetCDF lookup table (+ notebook to generate it).
• Updated fluid1d model to distribute heating radially after computing global heating; documentation refreshed accordingly.

Reviewed changes

Copilot reviewed 15 out of 16 changed files in this pull request and generated 12 comments.

Show a summary per file

File	Description
src/solid1d_relax.jl	New relaxation-based ODE/FD solver implementation for solid1d.
src/Obliqua.jl	Wires in solid1d-relax, updates core mass handling, Hansen k-range usage, and fluid1d selection/behavior.
src/solid1d.jl	Updates solid ODE matrices (inertia terms), adds scaling/QR stabilization and core mass in gravity calculation.
src/solid1d_mush.jl	Updates gravity calculation to include m_core.
src/Hansen.jl + res/hansen_k_table.nc	Switches k-range estimation to a NetCDF lookup table and adds a default table to res/.
res/config/all_options.toml + test/test.toml	Adds new configuration keys and documents the new solver option.
examples/hansen_k_table.ipynb	Provides a workflow to generate custom Hansen k-range tables.
docs/src/reference/*.md	Documents normalization/scaling, relaxation method, and fluid dissipation profiles; minor formatting/spelling updates.
src/plotting.jl	Adds a helper plot for relaxation solution visualization.
docs/src/development.md, docs/src/index.md, README.md	Documentation maintenance updates (module listing, links, etc.).

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

solve_radial_system applies non-dimensional scaling (ωs/rs/ρs/gs/μs/Ks) when building A/C/D matrices, but the core and surface boundary conditions are built from unscaled inputs (get_core_bc!(ω, r[1], …; G0=1) and get_surface_bc!(R_planet, n)). This mixes dimensional and non-dimensional systems, so the recursion matrices and final solve are not consistent. Boundary conditions should be computed in the same scaled variables (use ωs/rs/… and the same G0, and apply the S/Sinv normalization consistently), or remove the scaling altogether.

Suggested change

	# boundary conditions
	B1 = get_core_bc!(ω, r[1], ρ_core, g[1], μ[1], K[1], core, n; G0=1)
	BN, b = get_surface_bc!(R_planet, n)
	ρ_cores = ρ_core / ρ0
	# boundary conditions must be expressed in the same scaled and
	# normalized state variables as the recursion matrices
	B1 = get_core_bc!(ωs, rs[1], ρ_cores, gs[1], μs[1], Ks[1], core, n; G0=G0) * Sinv
	BN, b = get_surface_bc!(R_planet / R0, n)
	BN = BN * Sinv

[Copilot]

Docstrings in this section still describe the old function signatures (e.g., compute_M(r, ρ, g, μ, K, n; ...) and get_B_product!(..., r, ρ, g, μ, K, n)), but the implementations now require ω (and accept G0). Updating the docstrings to reflect the new required arguments will prevent confusion for callers and for generated docs.

[MarijnJ0]

Hey, I have made the final updates. Additional to the solid1d_relax model there is now also the solid1d_mush_relax model. The former is working properly. The latter is mostly working; specifically, the boundary components / solid propagator component and interface components all appear to work. The mush propagator is still outputting some unphysical results, where the solution oscillates instantaneously between extrema. This may be resolved by using the proper non-dimensionalization. However, this is also still not working entirely as expected. I think that for now we can proceed with just the solid1d_relax model, and in due time also fix the last few bugs in the solid1d_mush_relax model. (The solid1d_relax model works fine without non-dimensional scaling.)

I have also updated the surface boundary condition, it is now given by the most general form that applies to both tides and surface loads. I made sure to verify this with other works.

The odd (r/R)^2 term in the heating profile scaling has been moved to the k_n Lovenumber. I realized that this term arises from the fact that the solid heating evaluated upto r < R, is combined with a tidal potential evaluated at R. Similarly k_n is combined with the same tidal potential as well as a prefactor. Overall this yields the observed scaling mismatch. I believe the correction might need to still be changed. I will think about this some more and update the code accordingly before we merge this PR.

For now, it would be ideal if we can review the solid1d_relax implementation. Please do not spend too much time on the solid1d_mush_relax model yet.

Lmk if there is anything I can do to help/clarify! :)

@nichollsh @timlichtenberg

[nichollsh]

Thanks, @MarijnJ0. I'll review the solid1d_relax implementation soon and then get back to the rest at a later stage. Please let me know later when you'd like us to review the remaining parts.

One initial comment: the model is now becoming very robust, and also very complicated. This means that I am only just keeping up with its inner workings and will need to think hard about the physics. This partly reflects the complexity of tides modelling (I think). So, in addition to the new explanations in the docs, could you perhaps consider adding a diagram to aid with the model description? Even a sketch similar to the one you showed in our last meeting.

[MarijnJ0]

Thanks @nichollsh, I will send a more polished overview of all the eqs and functions sometime this week! The current structure of the relaxation-based modules is exactly like in the diagram I showed during the meeting, I will make a nice latex version and add this as well. Regardless, the relevant equations should all be in the new docs and the PDF i sent last week over discord, so hopefully this can also help for now.

[nichollsh]

I tried running Obliqua as described in the documentation, and got an error about the calc_solid_tides function not existing. I think the documentation tutorials are out of date.

Anyway, I looked at the runtests.jl file and tried following the code layout from there in the Julia REPL. The solid1d_relax model seems to run and produce the plot below.

I also tried the solid1d_mush model and got some weird behaviour. The output says:

[ Info: Expected bulk heating: 7.283377429171197481236606466699505636155722518266908299731777758125457597271871e+08
[ Info: Obtained bulk heating: 2.799762844091837278996729935972998560593428872895119572803561123918216564506373e+17

and I get the plots below...

Why is the model dissipating so much heat at this near-surface layer? I wondered if this was because of how the viscosity profile is defined, but I plotted the viscosity profile and can see that it's not the case...

Not sure what's up with that, but let me know when you want us to look at this other mush module.

Also, I wonder if you could improve the performance of the model by using Float64 rather than BigFloat? Not sure if this will introduce instabilities, but BigFloat is 40 bytes while Float64 is 'only' 8 bytes. Definitely worth trying, given that it's already inconsistently using the two types throughout the code anyway.

[nichollsh]

@MarijnJ0 please let me know when this needs another review

[MarijnJ0]

@nichollsh, to answer your questions from before:

• I have updated the docs, indeed these functions have all been replaced by run_tides
• The output on the first plot looks good, here you are using a low solidus viscosity, hence the dissipation nearly reaches the surface. Here you really need the relaxation based solver.
• This shows you that the shooting method completely diverges at low viscosities. The result is indeed completely useless. However, the shooting method is very fast, so if you can model only solid tides (i.e. solidus visc > 1e15), than this model should still work for some intermediate forcing frequencies.
• To add to this, I have the new relaxation based mush solver working now, it is internally consistent for a wide range of forcing frequencies. It can be reviewed now.
• Regarding the numerical precision, I have not yet tested with lower precision, however given that we use prec and precc throughout the code, we can simply redefine them to be e.g. float64. We can definitely test this more now that the relaxation models are working. (If I used a different type than either its a typo, or there is a function that does not support bigfloat. In the latter case I need to check and find a better package. Specifically, this is the case for the inertial core.)

Besides these points, I have removed the global variables as requested in the previous PR. I have added (currently commented-out lines) the ability to plot the heating map at the surface. I have tested and compared the new relaxation models, here are the outputs for an Earth-like planet throughout solidification:
solid1d_relax

and solid1d_mush_relax

I have created the figure describing the general structure of the latter model:

And I have created the flowchart for the main run_tides function in Obliqua

For detailed description and analysis please refer to the reader I sent on Discord.

The tests are currently passing, so I think we can do the final review and hopefully merge these changes!

[nichollsh]

Some minor comments, but generally looks good to me.

Please can you also incorporate my suggestions about the Project.toml file dependencies? i.e.

GenericLinearAlgebra = "0"
SpecialFunctions = "2.7"

I have also run the tests on my laptop and they pass successfully (Julia 1.11.9 on Fedora 44). However, it's unclear what the code-coverage percentage is at the moment. We should add this at some point, potentially adopting the same testing framework as AGNI.

given that we use prec and precc throughout the code, we can simply redefine them to be e.g. float64

The code doesn't consistently use prec and precc. It should probably be updated to use these consistently, and then define them equal to Float64 (or another type) as appropriate.

Please request another review when you are ready. Thanks, again!

[timlichtenberg]

This looks great Marijn, amazing work! :) I believe there may be an inconsistency in how the updated BC is implemented. I think there is a disagreement between the tech document you sent us and some places in the code:

(a) Sec. 11 says we should not include the Korenaga-style ζ/4 terms in the boundary, which are the 4πGζ contributions to b[6] in solid1d_relax.jl. Do I understand this correctly?
(b) solid1d.jl:425 and solid1d_mush.jl:588 write b[3] = (2n+1)/r · U − 4πGζ.
(c) Relax solvers solid1d_relax.jl:508 and solid1d_mush_relax.jl:977 (M==8 branch) then: b[6] = (2n+1)/R · U + 4πGζ/G0.
(d) M==6 branch at solid1d_mush_relax.jl:989 does the same expression as (c) but with −.

ζ = (2n+1)/(4πGR)U_prime is thus the same throughout but with different signs apparently. So the new relax solvers are seemingly going in the opposite direction to the shooting reference everywhere except in the M==6 branch of solid1d_mush_relax.jl. The M==6 branch and the M==8 branch of the same file seem to enforce opposite signs at the same surface. Which sign is correct?

Once this is resolved you could add a test where solid1d-relax and solid1d are run on the same interior with m_core=0, no inertia, single-frequency forcing, and Im(k2) is required to agree to a few percent. This could help enforcing that the BC is correct.

[MarijnJ0]

Hey @nichollsh and @timlichtenberg, I have implemented the requested changes/suggestions.

I think that the package versions were already updated in the previous commits, but if the issue persists than please do let me know Harrison. Besides this, I have changed the y-function loop to use prec and precc consistently. The heating profile loop now uses float64 and complexf64 exclusively, which has improved the speed significantly. There is still one function that needs more work in the future, namely the inertial core solution, which relies on the spherical bessel function. To add to this, I have implemented a low frequency asymptotic approximation so that it can be used without issues, despite this the results are still incorrect. For now we can use the liquid core solution instead.

Regarding the surface boundary conditions, I have fixed the plus-minus inconsistency. This did not affect the previously generated output, since it was part of the loading term and not the tidal term. I am not sure if the zeta used by Korenaga is the same zeta we use to denote the surface load. Based on the derivation provided in the PDF I strongly doubt it, since it seems to correspond to the tidal potential (not load potential).

In order to compare the boundary conditions between the shooting and relaxation based solvers, we would need to model a fluid on top of a solid, without the interp model inbetween. This is not possible since the shooting models will diverge before reaching the fluid. Nevertheless, I have made some tweaks to the interior data file used by the tests, and the output between the different solid1d models is now within good agreement.

I agree with the comment on code-coverage percentage, it will be a good change to implement. Maybe this can be a part of the next PR? There are still some small tweaks that need to be made to for example the summation of the heating profiles as discussed in our previous meeting. As well as the (n,m,k) based loop. Since these are all relatively small changes they can likely be included in the next PR.

Hopefully this answers all remaining questions, if there are any remaining concerns do lmk!

[nichollsh]

Thanks, @MarijnJ0 - excellent work!

I absolutely agree that the code coverage task can be considered in a follow-up PR. Also, the package dependency issues have indeed been fixed. Tests pass on my laptop with Julia 1.12.

1) Precision

I am a bit confused about the usage of prec and precc. I see that the floating point types are now somewhat more consistent throughout the code, but they are still defined to be BigFloat, so it surprises me that there's any performance difference. Also, there's still several places where the types are hard-coded at Float64 (e.g. in solid1d_mush). Not really a problem, and it's great that this has already improved performance, but could you clarify this please?

2) Core solution

You say that "... despite this the results are still incorrect. For now we can use the liquid core solution instead."

Does this mean that the solid-core option does not work at all, at the moment?

3) Docstrings

In several places, the function docstrings are not consistent with the function parameters. In compute_solid_lovenumbers(), resample_profiles(), propagate_solid(), as three examples.

[MarijnJ0]

@nichollsh, I have made some final tweaks to the code. To answer your questions, the change from Bigfloat to Float64 sped up calculations, specifically when calculating the heating profile. This is now handled purely using Float64. Additionally, the types should now be consistent throughout (but maybe I missed something), so we can simply set prec and precc at the top of the .jl files to whatever precision is needed. I tested Float64 with the solid1d_relax model, which seems to work fine and is much faster.

The liquid and solid core solution vectors are working as expected. However, the new inertial core solution vector does not behave as expected. I have to admit that I have not had the time to properly derive it from first principles, and as such it maybe simply be invalid for the forcing frequency range im trying to apply it to. Alternatively, there may be an error in the equations, I will have to check this at some point.

I have further updated the docstrings.

Let me know if there are any final changes that need to be made before we can merge this!

[nichollsh]

Fantastic work! Thanks for this @MarijnJ0.

I tried the tutorials with the updated documentation, and it all seems to work well. The plots make sense to me, apart from the point you've noted regarding how to handle the core BC.

Great to hear that performance is improved by simply using Float64 throughout. I think there's probably still room for more performance improvements, so maybe using a profiler would be useful.

For example, running run_tides() with the @time macro reported:

22.583407 seconds (372.41 M allocations: 37.394 GiB, 22.78% gc time)

This is a lot of memory allocations (40 GB!) - maybe there are arrays being allocated in a loop, or data files being repeatedly read-in.

This PR is already large enough, and a great improvement, so am happy for this to be merged.

[MarijnJ0]

Great to see this get merged!

I was not aware that the code was using this much memory, maybe this is related to the fact that Bigfloat dumps some of its info to memory, in which case switching to Float64 should already help reduce this a lot. I will definitely look into this some more. One thing to note is that I tried using the multifloats package, but it did not work. It seems to not define arctan for its float types, maybe we can work around this, or potentially try one of the slower options. Ultimately, I think the aim should be to use Float64 for both relaxation based solvers, which will eliminate this issue completely.

[nichollsh]

I don't think it is using 40 GB of memory at once, but is doing 40 GB of memory allocations/frees during a single run. This is still expensive and is probably due to arrays being created inside loops.

No worries about MultiFloats. I agree that sticking with just the native Float64 is ideal anyway, since it also means fewer dependencies.