Md/stabilize solid1d

Project.toml

@@ -1,28 +1,37 @@
 name = "Obliqua"
 uuid = "70fbf7c6-2f8c-431c-b73b-6bd747db0edf"
-authors = ["Marijn van Dijk <m.r.van.dijk.3@student.rug.nl>", "Tim Lichtenberg <tim.lichtenberg@rug.nl>", "Harrison Nicholls <harrison.nicholls@physics.ox.ac.uk>", "Hamish Hay <hamish@tides.rocks>"]
+authors = ["Marijn van Dijk <m.r.van.dijk.3@student.rug.nl>", "Tim Lichtenberg <tim.lichtenberg@rug.nl>", "Harrison Nicholls <harrison.nicholls@physics.ox.ac.uk>"]
 version = "0.1.0"
 
 [deps]
 AssociatedLegendrePolynomials = "2119f1ac-fb78-50f5-8cc0-dda848ebdb19"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
+NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 
 [compat]
 FFTW = "1.10.0"
+GenericLinearAlgebra = "0"
 Interpolations = "0.16.2"
 JSON = "1.3.0"
 LoggingExtras = "1.2.0"
+NCDatasets = "0.14.15"
 Plots = "1.41.2"
 Printf = "1.11.0"
+SparseArrays = "1.11.0"
+SpecialFunctions = "2.7"
 Statistics = "1.11.1"
 TOML = "1.0.3"

README.md

@@ -16,17 +16,6 @@ Forked from the [original repository](https://github.com/hamishHay/Love.jl) of H
 ### Documentation
 https://proteus-framework.org/Obliqua
 
-## Contributors
-
-| Name                    | Email address                         |
-| -                       | -                                     |
-| Marijn van Dijk         | m.r.van.dijk.3[at]student.rug.nl      |
-| Tim Lichtenberg         | tim.lichtenberg[at]rug.nl             |
-| Harrison Nicholls       | harrison.nicholls[at]physics.ox.ac.uk |
-| Mohammad Farhat         | farhat[at]berkeley.edu                |
-| Hamish Hay              | hamish[at]tides.rocks                 |
-
-
 ### Repository structure
 * `README.md`           - This file
 * `LICENSE.txt`         - License for modification, distribution, etc.

docs/make.jl

@@ -11,6 +11,7 @@ format = Documenter.HTML(
 # Build the docs
 makedocs(
     sitename="Obliqua",
+    checkdocs = :none,
     format=format,
     modules = [Obliqua],
     pages = [
@@ -19,17 +20,25 @@ makedocs(
         "Compass" => "compass.md",
         PageNode("Tutorials" => "tutorials/index.md", [
             "1 - Loading Data" => "tutorials/loading-data.md", 
-            "2 - Running Model" => "tutorials/running-model.md", 
-            "3 - Plotting" => "tutorials/plotting.md"
+            "2 - Configuring Model" => "tutorials/configuration-file.md", 
+            "3 - Running Model" => "tutorials/running-model.md", 
+            "4 - Plotting" => "tutorials/plotting.md"
             ]
         ),
         "How-to guides" => "how-to-guides/index.md",
         PageNode("Reference" => "reference/index.md", [
-            "1 - Rheology" => "reference/rheology.md", 
-            "2 - Solid-phase" => "reference/solid-phase.md", 
-            "3 - Mush layer" => "reference/mush-layer.md", 
-            "4 - Liquid-phase" => "reference/liquid-phase.md",
-            "5 - Forcing Frequency" => "reference/forcing-frequency.md",
+            "1 - Forcing Frequency" => "reference/forcing-frequency.md",
+            "2 - Rheology" => "reference/rheology.md", 
+            PageNode("3 -Solid-phase" => "reference/solid-phase.md", [
+                "Solid0d" => "reference/solid/solid0d.md",
+                "Solid1d" => "reference/solid/solid1d.md", 
+                "Solid1d-mush" => "reference/solid/solid1d_mush.md", 
+                "Solid1d-relax" => "reference/solid/solid1d_relax.md", 
+                "Solid1d-mush-relax" => "reference/solid/solid1d_mush_relax.md"
+                ]
+            ),
+            "4 - Mush layer" => "reference/mush-layer.md", 
+            "5 - Liquid-phase" => "reference/liquid-phase.md",
             "6 - Surface Loading" => "reference/surface-loading.md",
             "7 - Tidal Potentials" => "reference/tidal-potentials.md"
             ]

docs/src/development.md

@@ -23,6 +23,8 @@ Modules = [
     Obliqua.solid0d, 
     Obliqua.solid1d, 
     Obliqua.solid1d_mush, 
+    Obliqua.solid1d_relax, 
+    Obliqua.solid1d_mush_relax, 
     Obliqua.fluid0d,
     Obliqua.Hansen
 ]

docs/src/index.md

@@ -2,9 +2,6 @@
     <a href="https://opensource.org/license/mit">
         <img src="https://img.shields.io/badge/License-MIT-blue.svg">
     </a>
-    <a href="https://arxiv.org/abs/2507.11266">
-        <img src="https://img.shields.io/badge/arXiv-2507.11266-b31b1b">
-    </a>
 ```
 
 # Obliqua (Tidal heating model)
@@ -14,7 +11,7 @@ A Julia package to calculate the tidal deformation (i.e., tidal Love numbers) of
 Forked from the [original repository](https://github.com/hamishHay/Love.jl) of Hamish Hay. Distributed under the MIT License.
 
 ### Documentation
-https://fwl-proteus.readthedocs.io
+PROTEUS: [https://fwl-proteus.readthedocs.io](https://fwl-proteus.readthedocs.io)
 
 ### Repository structure
 * `README.md`           - This file

docs/src/reference/forcing-frequency.md

@@ -1,25 +1,27 @@
 
-### Reference (5)
+### Reference (1)
 
 # Forcing Frequency
 
-Given the fact that the tidal forcing magnitude decreases exponentially with harmonnic order, we limit the calculation to only the lowest harmonnic frequency ``n = 2``. Generally, in the limit of ``R_p = a``, it suffices to only consider the quadrupolar harmonic (``n = 2``). Moreover, as we are considering a coplanar geometry, terms with ``\ell = 1`` associated with obliquity tides vanish. The Fourier modes of the second harmonnic have frequencies given by
+Given the fact that the tidal forcing magnitude decreases exponentially with harmonic order, we may limit the calculation to only the lowest harmonic degree ``n = 2``. Generally, in the limit of ``R_p \ll a``, it suffices to only consider the quadrupolar harmonic (``n = 2``). Nevertheless, Obliqua can also determine higher degree contributions. As of now, we are considering a coplanar geometry. Hence, terms with ``m = 1`` associated with obliquity tides vanish. The Fourier modes of the second harmonic have frequencies given by
 
 ```math
-\sigma = 2\Omega - k n_{\mathrm{orb}},
+\sigma = m\Omega - k n_{\mathrm{orb}},
 ```
 
-where ``\Omega`` is spin rate and ``n_{\mathrm{orb}}`` orbital mean motion, and for integer values of ``k``. In our formalism tides are occuring over a large time interval, a time step ``\Delta t``. As such, we must account for tidal excitations that occur over a wide range of frequencies. We calculate the imaginary part of the second harmonnic degree (``k_2``) Love number (``\Im[k_{2}(\sigma)]``) for a handful of frequencies. The resulting descrete profile can subsequently be interpolated across positive and negative frequencies to yield a continous profile. 
-
-Note, since we are only interested in the the Fourier modes of the second harmonnic, we don't actually use the entire spectrum. We can determine which frequencies are relevant using the Hansen Coefficients. 
+where ``\Omega`` is spin rate and ``n_{\mathrm{orb}}`` orbital mean motion, and for integer values of order ``-2 \leq m \leq 2`` and harmonic ``\infty \leq k \leq \infty``. In our formalism tides are occuring over a large time interval, a time step ``\Delta t``. As such, we must account for tidal excitations that occur over a wide range of frequencies. We calculate the imaginary part of the ``n``th harmonic degree (``k_n``) Love number (``\Im[k_{n}(\sigma)]``) for all relevant harmonnic frequencies for which the Hansen coefficient 
 
 ```math
-X_{22k}(e) =
+X^{-(n+1), m}_k(e) =
 \frac{1}{2\pi} \int_0^{2\pi}
-\left(\frac{r}{a}\right)^2
-e^{2i v - ikM}\,dM,
+\left(\frac{r}{a}\right)^n
+e^{im\Omega - ikn_{\mathrm{orb}}}\,dn_{\mathrm{orb}} \geq 0.01,
 ```
 
-where ``v`` is the true anomoly. 
+(i.e. ~1% corrections). This implies that we are considering the following set $K$ of harmonnic frequencies $k$:
+
+```math
+\{k \in K \, \forall \, k : X^{-(n+1), m}_k \geq 0.01\ | k \in Z\}
+```
 
 ---

docs/src/reference/index.md

@@ -4,7 +4,7 @@ In this component you will be guided through the employed formalism and theoreti
 
 - [Rheology](@ref)
 - [Solid-Phase](@ref)
-- [Mush layer](@ref)
+- [Mush layer - interp](@ref)
 - [Liquid-Phase](@ref)
 - [Forcing Frequency](@ref)
 - [Surface Loading](@ref)

docs/src/reference/liquid-phase.md

@@ -1,5 +1,5 @@
 
-### Reference (4)
+### Reference (5)
 
 # Liquid-Phase
 
@@ -53,9 +53,9 @@ This corresponds to the simplest assumption that turbulence or small-scale mixin
 
 The exponential profile concentrates dissipation near the lower boundary and decreases with height:
 
-[
+```math
 D(z) \propto e^{-z/H_R}
-]
+```
 
 This can approximate scenarios where dissipation is strongest near the interface with the underlying solid layer, for example when tidal flows interact with boundary roughness or generate shear-driven turbulence.
 
@@ -69,19 +69,19 @@ This provides a simple way to model dissipation that is concentrated toward the
 
 The quadratic profile falls off more steeply with height:
 
-[
+```math
 D(z) \propto (1 - z/H_R)^2
-]
+```
 
 This concentrates dissipation even more strongly toward the base of the layer. It can approximate cases where turbulence is primarily generated near the boundary and decays rapidly away from it.
 
 ### Dynamic dissipation
 
 The dynamic profile introduces a depth-dependent mixing length
 
-[
+```math
 \ell_{\text{mix}} = \min(z, H_R)
-]
+```
 
 so that the dissipation scale adjusts with distance from the boundary.
 

docs/src/reference/mush-layer.md

@@ -1,4 +1,44 @@
 
-### Reference (3)
+### Reference (4)
 
-# Mush layer
+# Mush layer - interp
+
+The `interp` model is a simplified approach to modeling dissipation within a planetary region where active tidal equations (like those in `solid1d`) are not explicitly solved, such as a thin mushy layer or a transition zone. Instead of solving the full poro-viscoelastic system, it approximates the heating profile and subsequent Love numbers using exponential decay from the layer boundaries.
+
+### Methodology
+
+The model assumes that tidal dissipation peaks at the interfaces (upper and lower) and decays exponentially into the interior of the segment.
+
+#### 1. Heating Profile Generation
+The heating profile $P(r)$ is constructed as a superposition of two exponential decay functions originating from the top ($r_{top}$) and bottom ($r_{bot}$) interfaces:
+
+$$P(r) = P_t \exp\left( -\frac{|r - r_{top}|}{l_t} \right) + P_b \exp\left( -\frac{|r - r_{bot}|}{l_b} \right)$$
+
+where:
+*   **$P_t, P_b$**: The heating intensities at the top and bottom interfaces.
+*   **$l_t, l_b$**: The characteristic decay lengths, defined as a fraction of the total segment thickness ($l = \text{width} \cdot \Delta R$).
+
+
+
+#### 2. Inferring Love Numbers
+To remain consistent with the energy dissipation within the global model, the imaginary part of the Tidal Love number $k_2^T$ is inferred directly from the generated power profile. For a shell with volume $V_s$, the contribution to the complex Love number is:
+
+$$\text{Im}(k_{2,s}^T) = -\frac{5 R \omega}{8\pi G} (P_s V_s)$$
+
+The total $k_2^T$ for the segment is the sum of these contributions across all radial shells. Note that in this interpolation mode, the load Love number $k_2^L$ is typically returned as zero (or ignored) as the model focuses purely on energy dissipation rather than structural loading response.
+
+---
+
+### Function Documentation
+
+```@docs
+Obliqua.run_interp
+```
+
+### Configuration Parameters
+When using this model via the TOML configuration, the following parameters in `[orbit.obliqua.mushy]` are relevant:
+
+*   `t_width`: Controls the decay length from the top interface ($l_t$).
+*   `b_width`: Controls the decay length from the bottom interface ($l_b$).
+
+This model is particularly useful for representing "mushy" regions where the physical properties are highly uncertain, but the dissipation is expected to be concentrated near the boundaries of solid or liquid layers.