From 09c2cba92da445785c4876c285e2fd3ce7b9960f Mon Sep 17 00:00:00 2001 From: Marijn Date: Tue, 8 Sep 2026 15:02:15 +0000 Subject: [PATCH 01/11] Added small bug fix, removes ~10% discrepancy between shooting and relax solvers. --- src/solid1d.jl | 3 ++- src/solid1d_mush.jl | 5 +++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/solid1d.jl b/src/solid1d.jl index a17df88..4f2e017 100644 --- a/src/solid1d.jl +++ b/src/solid1d.jl @@ -321,7 +321,8 @@ module solid1d Ks = K ./ μ0 ωs = ω / ω0 - y_start = get_Ic(ωs, rs[end,1], ρ_core/ρ0, gs[end,1], μ_core/μ0, κ_core/μ0, core, n; G0=G0, Y=[1,2,3,4,5,6]) + # Core basis must be evaluated at the core-mantle boundary itself (rs[1,1]) + y_start = get_Ic(ωs, rs[1,1], ρ_core/ρ0, gs[1,1], μ_core/μ0, κ_core/μ0, core, n; G0=G0, Y=[1,2,3,4,5,6]) y1_4 = zeros(precc, 6, 3, nsublayers-1, nlayers) # Three linearly independent y solutions diff --git a/src/solid1d_mush.jl b/src/solid1d_mush.jl index 57b9640..69afaee 100644 --- a/src/solid1d_mush.jl +++ b/src/solid1d_mush.jl @@ -482,11 +482,12 @@ module solid1d_mush ks = k./R0^2 # Define starting vector as the core solution matrix, Y_r_C (Eq. S5.15) + # Core basis must be evaluated at the core-mantle boundary itself (rs[1,1]) if porous_layer[1] - y_start = get_Ic(ωs, rs[end,1], ρ_core/ρ0, gs[end,1], μ_core/μ0, κ_core/μ0, core, n; G0=G0, Y=[1,2,3,4,5,6,7,8]) + y_start = get_Ic(ωs, rs[1,1], ρ_core/ρ0, gs[1,1], μ_core/μ0, κ_core/μ0, core, n; G0=G0, Y=[1,2,3,4,5,6,7,8]) else y_start = zeros(precc, 8, 4) - y_start[1:6, 1:3] .= get_Ic(ωs, rs[end,1], ρ_core/ρ0, gs[end,1], μ_core/μ0, κ_core/μ0, core, n; G0=G0, Y=[1,2,3,4,5,6]) + y_start[1:6, 1:3] .= get_Ic(ωs, rs[1,1], ρ_core/ρ0, gs[1,1], μ_core/μ0, κ_core/μ0, core, n; G0=G0, Y=[1,2,3,4,5,6]) end y1_4 = zeros(precc, 8, 4, nsublayers-1, nlayers) # Four linearly independent y solutions From dfa33ce3d739221e980e99a17e1408927c6696a2 Mon Sep 17 00:00:00 2001 From: Marijn Date: Tue, 8 Sep 2026 18:12:40 +0000 Subject: [PATCH 02/11] Updated tests for solid1d and solid1d_mush. --- test/test_obliqua.jl | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/test/test_obliqua.jl b/test/test_obliqua.jl index 5abc655..fa9e47d 100644 --- a/test/test_obliqua.jl +++ b/test/test_obliqua.jl @@ -44,9 +44,9 @@ using Obliqua.solid1d_relax.common perm, phi = Obliqua.interior.limit_porosity(perm, phi, cfg) bulkd = Obliqua.interior.get_drained_bulk(bulk, phi, cfg) - power_prf_expt = [5.947109168352234e-17, 3.1933710995124006e-18, 3.163189030611521e-18, 3.0313968498306403e-18, 2.957234184019153e-18, 2.901052165774541e-18, 7.139262688501286e-20, 8.116278641807716e-21, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - power_blk_expt = 1.116704658520055e6 - imag_k2_expt = [0.0015300931674779095] + power_prf_expt = [5.811903654413909e-17, 3.1331643176818168e-18, 3.1025899880099643e-18, 2.9724177811069162e-18, 2.898779016027967e-18, 2.8427679069919094e-18, 6.972115536420127e-20, 9.094049487001497e-21, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + power_blk_expt = 1.093766208671846e6 + imag_k2_expt = [0.0014986632230270696] power_prf, power_blk, _, _, LNk = Obliqua.run_tides( omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, bulkd, phi, perm, cfg @@ -72,9 +72,9 @@ using Obliqua.solid1d_relax.common perm, phi = Obliqua.interior.limit_porosity(perm, phi, cfg) bulkd = Obliqua.interior.get_drained_bulk(bulk, phi, cfg) - power_prf_expt = [7.093015825910385e-17, 7.866205023233557e-22, 7.828755302672636e-22, 7.1630669966287245e-22, 6.930940156995041e-22, 6.796575518404947e-22, 2.067944736376056e-23, 2.6792800694886212e-24, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - power_blk_expt = 1.043240911950073e6 - imag_k2_expt = [0.0014294341652731294] + power_prf_expt = [6.930763584609395e-17, 7.717943205767067e-22, 7.678856966385867e-22, 7.023802241566042e-22, 6.794063964493952e-22, 6.660180949962994e-22, 2.0194178851454583e-23, 2.999024557829818e-24, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + power_blk_expt = 1.0210408255188541e6 + imag_k2_expt = [0.0013990159160908932] power_prf, power_blk, _, _, LNk = Obliqua.run_tides( omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, bulkd, phi, perm, cfg @@ -175,9 +175,9 @@ using Obliqua.solid1d_relax.common perm, phi = Obliqua.interior.limit_porosity(perm, phi, cfg) bulkd = Obliqua.interior.get_drained_bulk(bulk, phi, cfg) - power_prf_expt = [5.523694183459497e-17, 2.9026418454483742e-18, 2.8798176332637832e-18, 2.7648964196486443e-18, 2.7020842662437843e-18, 2.655345585757709e-18, 6.840499967001757e-20, 1.5011935200727582e-20, 1.5435493282770822e-19, 3.7231798866713286e-19, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - power_blk_expt = 1.0317994367738063e6 - imag_k2_expt = [0.0014137572153656449] + power_prf_expt = [5.3950076375562126e-17, 2.8407437126296814e-18, 2.818017456639662e-18, 2.705284895036029e-18, 2.6435353304508365e-18, 2.5975027640527154e-18, 6.681457041216717e-20, 1.6896966982997684e-20, 1.7546748749581933e-19, 4.269886897967008e-19, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + power_blk_expt = 1.0100820769031402e6 + imag_k2_expt = [0.0013840003913923272] power_prf, power_blk, _, _, LNk = Obliqua.run_tides( omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, bulkd, phi, perm, cfg @@ -203,9 +203,9 @@ using Obliqua.solid1d_relax.common perm, phi = Obliqua.interior.limit_porosity(perm, phi, cfg) bulkd = Obliqua.interior.get_drained_bulk(bulk, phi, cfg) - power_prf_expt = [6.59329217956111e-17, 7.150225033540763e-22, 7.127417705214227e-22, 6.5332044304041115e-22, 6.332667084217917e-22, 6.220528058460307e-22, 1.9820333213539122e-23, 4.993102056808919e-24, 5.802505797821157e-23, 1.1587864466255244e-19, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - power_blk_expt = 962513.5224966361 - imag_k2_expt = [0.0013188226207715324] + power_prf_expt = [6.439313496369688e-17, 6.997885032837296e-22, 6.974614792513135e-22, 6.392492308651947e-22, 6.195600031186184e-22, 6.085179544465906e-22, 1.9359067755994564e-23, 5.617084352918099e-24, 6.593993133185726e-23, 1.3445953880034122e-19, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + power_blk_expt = 941057.7310122688 + imag_k2_expt = [0.00128942419415749] power_prf, power_blk, _, _, LNk = Obliqua.run_tides( omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, bulkd, phi, perm, cfg @@ -376,9 +376,9 @@ using Obliqua.solid1d_relax.common perm, phi = Obliqua.interior.limit_porosity(perm, phi, cfg) bulkd = Obliqua.interior.get_drained_bulk(bulk, phi, cfg) - power_prf_expt = [2.305906768268315e-14, 1.3218894513174959e-15, 1.3094017433075674e-15, 1.254865626054167e-15, 1.2241749182735129e-15, 1.2009254964452314e-15, 2.9551938004902827e-17, 3.360326218245033e-18, 3.4664482750907783e-18, 1.1922878409360185e-18, 4.083645111622018e-19, 1.3928661517229027e-19, 4.731391426312715e-20, 1.600697002079666e-20, 5.393770965404668e-21, 1.8103406061654148e-21, 6.052478834142824e-22, 2.015726504449689e-22, 6.687673376567115e-23, 2.2104589027693444e-23, 7.279010165466779e-24, 2.388158711559352e-24, 7.806786714731434e-25, 2.5428252010576644e-25, 8.253008730303911e-26, 2.669166668912046e-26, 8.602451848695603e-27, 2.7629137952499543e-27, 8.843539405574548e-28, 2.8210607331075036e-28, 8.968920682266818e-29, 2.841990858390102e-29, 8.975793708031064e-30, 2.8255545395877367e-30, 8.865994977508642e-31, 2.7730440198682795e-31, 8.645752829854655e-32, 2.687056276502569e-32, 8.325100636241318e-33, 2.5712953420809534e-33, 7.917252866312166e-34, 2.430341490098225e-34, 7.437738796963028e-35, 2.2693777370451804e-35, 6.903977026653167e-36, 2.0955861618165803e-36, 6.391194744399865e-37, 2.1077576736530604e-37, 1.2393351176789043e-37, 2.4075842948919993e-37, 7.584848660826297e-37, 2.5477707351847457e-36, 8.632137324650337e-36, 2.934797275755495e-35, 1.0007743921631831e-34, 3.4227100678540193e-34, 1.1740163363934534e-33, 4.0387232629503395e-33, 1.3933977862639133e-32, 4.821293511853955e-32, 1.673038470279199e-31, 5.822358436804119e-31, 2.03207513562298e-30, 7.112555260102203e-30, 2.49663627497486e-29, 8.788700196651061e-29, 3.102652862245179e-28, 1.098446762706548e-27, 3.899968789874166e-27, 1.3886054980816842e-26, 4.9582915455296047e-26, 1.7754998382537867e-25, 6.3759584027734245e-25, 2.2961899059708784e-24, 8.292933678414047e-24, 3.003644687871403e-23, 1.091016957195519e-22, 3.9743010232373613e-22, 1.4519162678640861e-21, 5.319619756058623e-21, 1.9547227580361233e-20, 7.2038495216176e-20, 2.66274826773369e-19, 9.87177345285578e-19, 3.671157510452736e-18, 1.3694762279275003e-17, 5.124745890875036e-17, 1.9239126084706928e-16, 7.246566861091503e-16, 2.7388294901300264e-15, 1.0388384225614868e-14, 1.0532040117949126e-14, 1.069254481915923e-14, 1.0873852449000426e-14, 1.108088550763046e-14, 9.18918563012912e-15, 7.202303249026087e-15, 5.674876263134549e-15, 4.509006923637211e-15] - power_blk_expt = 2.8408034975037627e9 - imag_k2_expt = [0.01161526252837546, 0.011119738496597285, 0.010623311307753108, 0.010125985043568942, 0.009627770742160606, 0.00912868799587736, 0.008628767019528763, 0.008128051376163703, 0.007626601642311584] + power_prf_expt = [2.253416932569662e-14, 1.2969687152089551e-15, 1.284317849973822e-15, 1.2304513290501027e-15, 1.1999767161471653e-15, 1.1767973140210735e-15, 2.885987911789033e-17, 3.7652262670483354e-18, 3.884135423480617e-18, 1.3359516918923412e-18, 4.575700941205927e-19, 1.5606985394675695e-19, 5.3014969741074366e-20, 1.793571816061222e-20, 6.0436894510768036e-21, 2.0284762728179813e-21, 6.78176784245978e-22, 2.2586099946282244e-22, 7.493499686480194e-23, 2.476806530791136e-23, 8.156089169057235e-24, 2.6759181480123464e-24, 8.747459767433619e-25, 2.849221062988066e-25, 9.247448978256785e-26, 2.9907859535632453e-26, 9.638997989499756e-27, 3.095829071291302e-27, 9.909135215132926e-28, 3.1609823818779192e-28, 1.0049624217014523e-28, 3.184434467291166e-29, 1.0057325402989003e-29, 3.166017666281719e-30, 9.93429655477572e-31, 3.1071799299123705e-31, 9.68751649069829e-32, 3.010831156328364e-32, 9.328227545596274e-33, 2.881121692652761e-33, 8.871236451606357e-34, 2.723183691518359e-34, 8.333943096732691e-35, 2.542822786970244e-35, 7.73580183189588e-36, 2.3478770145423524e-36, 7.154129916797755e-37, 2.337795304623464e-37, 1.308493285083435e-37, 2.4283157772817694e-37, 7.591045451357273e-37, 2.54795543236148e-36, 8.632192217828748e-36, 2.934798902609151e-35, 1.0007744402424495e-34, 3.422710082023477e-34, 1.1740163368098833e-33, 4.038723263072389e-33, 1.393397786267481e-32, 4.8212935118549946e-32, 1.6730384702792294e-31, 5.822358436804129e-31, 2.0320751356229802e-30, 7.112555260102203e-30, 2.49663627497486e-29, 8.788700196651061e-29, 3.102652862245179e-28, 1.098446762706548e-27, 3.899968789874167e-27, 1.3886054980816842e-26, 4.9582915455296047e-26, 1.7754998382537867e-25, 6.3759584027734245e-25, 2.2961899059708784e-24, 8.292933678414047e-24, 3.003644687871403e-23, 1.091016957195519e-22, 3.974301023237362e-22, 1.4519162678640861e-21, 5.319619756058623e-21, 1.9547227580361233e-20, 7.2038495216176e-20, 2.66274826773369e-19, 9.871773452855781e-19, 3.671157510452736e-18, 1.3694762279275003e-17, 5.1247458908750365e-17, 1.9239126084706933e-16, 7.246566861091503e-16, 2.7388294901300264e-15, 1.0388384225614868e-14, 1.0532040117949126e-14, 1.069254481915923e-14, 1.0873852449000426e-14, 1.108088550763046e-14, 9.18918563012912e-15, 7.202303249026087e-15, 5.674876263134549e-15, 4.509006923637211e-15] + power_blk_expt = 2.8317958612079725e9 + imag_k2_expt = [0.011581176958188387, 0.01108652173121107, 0.010590988673882713, 0.010094581843049406, 0.009597312117327252, 0.009099198763898173, 0.008600271466978374, 0.008100573001742285, 0.007600162830519118] power_prf, power_blk, _, _, LNk = Obliqua.run_tides( omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, bulkd, phi, perm, cfg From 76060a48f8102a9633deeeab51ba938ea8b0f75a Mon Sep 17 00:00:00 2001 From: Marijn Date: Wed, 9 Sep 2026 20:00:54 +0000 Subject: [PATCH 03/11] Updated docs/reference. --- docs/src/how-to-guides/config_file.md | 31 +- docs/src/reference/configuration-file.md | 295 ++++++++---- docs/src/reference/forcing-frequency.md | 19 +- docs/src/reference/index.md | 67 ++- docs/src/reference/liquid-phase.md | 36 +- docs/src/reference/rheology.md | 16 +- docs/src/reference/solid-phase.md | 334 +++++++------- docs/src/reference/solid/solid0d.md | 12 +- docs/src/reference/solid/solid1d.md | 28 +- .../reference/solid/solid1d_equil_relax.md | 12 +- docs/src/reference/solid/solid1d_mush.md | 16 +- .../src/reference/solid/solid1d_mush_relax.md | 424 +++--------------- docs/src/reference/solid/solid1d_relax.md | 28 +- docs/src/reference/surface-loading.md | 16 +- docs/src/reference/tidal-potentials.md | 12 +- 15 files changed, 654 insertions(+), 692 deletions(-) diff --git a/docs/src/how-to-guides/config_file.md b/docs/src/how-to-guides/config_file.md index 3376e27..b3da7cd 100644 --- a/docs/src/how-to-guides/config_file.md +++ b/docs/src/how-to-guides/config_file.md @@ -53,28 +53,6 @@ This block controls output and logging. --- -### Stellar Parameters - -Defines the host star. - -```@raw html -

config [star]

-``` - -```@raw html -
-``` - -| NAME | TYPE | DESCRIPTION | -| :--- | :--- | :--- | -| `mass` | float | Stellar mass in solar masses ($M_\odot$). | - -```@raw html -
-``` - ---- - ### Tidal Model Parameters Controls the tidal response model. @@ -94,7 +72,6 @@ Controls the tidal response model. | `optimize_scales` | bool | Boolean flag to optimize scaling factors for numerical stability. | | `solid_shell` | bool | Boolean flag to add an infinitesimal solid shell around the core to couple y2 and y4 in fluid mantles. | | `min_frac` | float | Minimum segment fraction of total mantle before it is considered. | -| `max_frac` | float | Maximum segment fraction of total mantle before it is considered. | | `visc_l` | float | Liquid viscosity. | | `visc_lus` | float | Liquid-Mush handoff viscosity. | | `visc_s` | float | Solid viscosity. | @@ -105,12 +82,12 @@ Controls the tidal response model. | `N_sigma` | int | Number of sampled forcing frequencies. | | `p_min` | float | Minimum period ($\log_{10}$ kyr). | | `p_max` | float | Maximum period ($\log_{10}$ kyr). | -| `s_min` | int | Minimum Fourier mode. | -| `s_max` | int | Maximum Fourier mode. | +| `s_min` | int or `"none"` | Minimum Fourier mode. `"none"` derives it from the eccentricity relation. | +| `s_max` | int or `"none"` | Maximum Fourier mode. `"none"` derives it from the eccentricity relation. | | `material_mu` | str | Rheological model for shear modulus (`"andrade"`, `"maxwell"`, or `"elastic"`). | | `material_k` | str | Rheological model for bulk modulus (`"andrade"`, `"maxwell"`, or `"elastic"`). | | `alpha` | float | Andrade power-law exponent. | -| `module_solid` | str | Solid interior model (`"solid0d"`, `"solid1d"`, `"solid1d-relax"`, `"solid1d-mush"`, `"solid1d-mush-relax"`, or `"solid1d-equil-relax"`). | +| `module_solid` | str | Solid interior model (`"none"`, `"solid0d"`, `"solid1d"`, `"solid1d-relax"`, `"solid1d-mush"`, `"solid1d-mush-relax"`, or `"solid1d-equil-relax"`). | | `module_mushy` | str | Mushy layer model (`"none"` or `"interp"`). | | `module_fluid` | str | Fluid layer model (`"none"`, `"fluid0d"`, or `"fluid1d"`). | @@ -135,7 +112,7 @@ Controls the tidal response model. | `ncalc` | int | Number of radial layers (shooting method). | | `dr_min` | int | Minimum grid spacing for relaxation solver [m]. | | `dr_max` | int | Maximum grid spacing for relaxation solver [m]. | -| `core` | str | Core boundary condition (`"liquid"`, `"solid"`, `"inertial"`). | +| `core` | str | Core boundary condition (`"liquid"`, `"solid"`, `"inertial-liquid"`, `"inertial"`). | | `core_props` | str | Core properties (shear modulus, bulk modulus) to use for CMB boundary condition (`"core"`, `"mantle"`). | | `inertial_terms` | bool | Boolean flag to include inertial terms in the motion matrix. | | `bulk_l` | float | Liquid bulk modulus [Pa]. | diff --git a/docs/src/reference/configuration-file.md b/docs/src/reference/configuration-file.md index f6e6e62..b0e5804 100644 --- a/docs/src/reference/configuration-file.md +++ b/docs/src/reference/configuration-file.md @@ -6,118 +6,259 @@ CollapsedDocStrings = true # Configuration -The configuration files follow the conventions used within PROTEUS. The default `all_options.toml` file contains all available parameters. +The configuration files follow the conventions used within PROTEUS. The default `all_options.toml` file (`res/config/all_options.toml`) contains all available parameters together with their defaults; the full parameter table is also available on the [Configuration file](@ref) how-to-guide page. This page instead walks through the `[orbit.obliqua]` block topic by topic and links each group of parameters to the reference page that derives the underlying model. ### Globals -These define metadata and output behavior: -* **`title`**: Identifier for the simulation setup. -* **`version`**: Configuration file version for reproducibility. ---- +```@raw html +
+``` -### Execution Parameters (`[params]`) -This block controls output and logging: -* **`path`**: Directory where output files are stored. -* **`time`**: Current time of the simulation run in years (used for output file naming). -* **`logging`**: Logging level (e.g., `INFO`, `DEBUG`). -* **`plot_fmt`**: Output format for generated plots. +| NAME | TYPE | DESCRIPTION | +| :--- | :--- | :--- | +| `title` | str | Identifier for the simulation setup. | +| `version` | str | Configuration file version for reproducibility. | + +```@raw html +
+``` --- -### Stellar Parameters (`[star]`) -Defines the host star: -* **`mass`**: Stellar mass in solar masses ($M_\odot$). +### Execution Parameters ---- +```@raw html +

config [params.out]

+``` + +```@raw html +
+``` + +| NAME | TYPE | DESCRIPTION | +| :--- | :--- | :--- | +| `path` | str | Directory where output files are stored. | +| `time` | float | Current time of the simulation run in years (used for output file naming). | +| `logging` | str | Logging level (e.g. `"INFO"`, `"DEBUG"`). | +| `plot_fmt` | str | Output format for generated plots (`"png"` or `"pdf"` recommended). | -### Planetary Parameters (`[planet]`) -Defines the planet's physical properties: -* **`mass_tot`**: Total planetary mass in Earth masses ($M_\oplus$). +```@raw html +
+``` --- -### Tidal Model Parameters (`[orbit.obliqua]`) +### Tidal Model Parameters + Controls the tidal response model. -* **`store_3D`**: Boolean flag to store 3D tidal response data. -* **`enforce_ec`**: Boolean flag to enforce energy conservation in tidal response calculations. -* **`optimize_scales`**: Boolean flag to optimize scaling factors for numerical stability. -* **`solid_shell`**: Boolean flag to add an infinitesimal solid shell around the core to couple y2 and y4 in fluid mantles. + +```@raw html +

config [orbit.obliqua]

+``` + +```@raw html +
+``` + +| NAME | TYPE | DESCRIPTION | +| :--- | :--- | :--- | +| `store_3D` | bool | Store 3D tidal response (displacement, stress, strain) for each layer. Generates large output files; if `false`, radial profiles are stored instead. | +| `enforce_ec` | bool | Enforce energy conservation in tidal response calculations. Improves stability in fluid-mush cases at low forcing frequencies (< 1e-7 Hz); does not affect the Love numbers. | +| `optimize_scales` | bool | Optimize non-dimensionalization scales for the relaxation method, for numerical stability at low forcing frequencies. Do not combine with BigFloat precision. | +| `solid_shell` | bool | Insert an infinitesimal solid shell around the core to patch a $y_2$/$y_4$ decoupling instability in fluid layers. Only relevant for `solid1d-relax` or `solid1d-mush-relax`. | #### Rheology and Viscosity -* **`min_frac`**, **`max_frac`**: Minimum segment fraction of total mantle before it is considered. -* **`visc_l`**, **`visc_lus`**: Liquid and liquidus viscosities. -* **`visc_s`**, **`visc_sus`**: Solid and solidus viscosities. + +See [Rheology](@ref) for the underlying model. + +| NAME | TYPE | DESCRIPTION | +| :--- | :--- | :--- | +| `min_frac` | float | Minimal segment radius fraction before smoothing [dimensionless]. | +| `visc_l` | float | Pure liquid viscosity [Pa s]. | +| `visc_lus` | float | Liquidus viscosity to use [Pa s]. | +| `visc_s` | float | Pure solid viscosity [Pa s]. | +| `visc_sus` | float | Solidus viscosity to use [Pa s]. | +| `material_mu` | str | Rheological model for the complex shear modulus (`"andrade"`, `"maxwell"`, or `"elastic"`). | +| `material_k` | str | Rheological model for the complex bulk modulus (`"andrade"`, `"maxwell"`, or `"elastic"`). | +| `alpha` | float | Andrade power-law exponent (free parameter), only used for Andrade rheology. | #### Spectral and Forcing Parameters -* **`n`**: Radial dependence exponent in $(r/a)^n$. -* **`m`**: Tidal harmonic (e.g., $m=2$ for semidiurnal tides). -* **`spectrum`**: Frequency sampling strategy (`"full"` or `"adaptive"`). -* **`N_sigma`**: Number of sampled forcing frequencies. -* **`p_min`**, **`p_max`**: Period range ($\log_{10}$ kyr). -* **`s_min`**, **`s_max`**: Fourier mode range. - -#### Material Model -* **`material_mu`**: Rheological model for shear modulus (`"andrade"`, `"maxwell"`, or `"elastic"`). -* **`material_k`**: Rheological model for bulk modulus (`"andrade"`, `"maxwell"`, or `"elastic"`). -* **`alpha`**: Andrade power-law exponent. + +See [Forcing Frequency](@ref) for the underlying model. + +| NAME | TYPE | DESCRIPTION | +| :--- | :--- | :--- | +| `n` | array | Radial dependence exponent(s) in $(r/a)^n$; since $r \ll a$, only $n=2$ contributes significantly. | +| `m` | array | Tidal harmonic(s) of the true anomaly (e.g. $m=2$ semidiurnal, $m=1$ diurnal). | +| `spectrum` | str | Frequency sampling strategy: `"full"` samples the whole k2 spectrum, `"adaptive"` samples only the region of interest. | +| `N_sigma` | int | Number of probe frequencies to evaluate k2 at (used when `spectrum = "full"`). | +| `p_min` | float | Minimum period for orbital and axial frequencies [$\log_{10}$ kyr]. | +| `p_max` | float | Maximum period for orbital and axial frequencies [$\log_{10}$ kyr]. | +| `s_min` | int or `"none"` | Minimum tidal mode (Fourier index in mean anomaly). `"none"` derives it from the eccentricity relation. | +| `s_max` | int or `"none"` | Maximum tidal mode (Fourier index in mean anomaly). `"none"` derives it from the eccentricity relation. | + +#### Phase specific Tidal Models + +See [Tidal Models](@ref) for all included models. + +| NAME | TYPE | DESCRIPTION | +| :--- | :--- | :--- | +| `module_solid` | str | Solid-phase tidal model, see below. | +| `module_mushy` | str | Mushy-phase tidal model, see below. | +| `module_fluid` | str | Fluid-phase tidal model, see below. | + +```@raw html +
+``` --- -### Interior Structure Models +### Tidal Models -#### Solid Model (`module_solid`) -* **`solid0d`**: Homogeneous solid approximation. -* **`solid1d`**: Radially resolved structure. -* **`solid1d-relax`**: Relaxation-based solver (more stable). -* **`solid1d-mush`**: Includes partially molten regions. -* **`solid1d-mush-relax`**: Relaxation-based solver (more stable), includes partially molten regions. -* **`solid1d-equil-relax`**: Relaxation-based solver (more stable), specifically for equilibrium tides. +```@raw html +

models [orbit.obliqua.module_solid]

+``` -#### Mushy Layer (`module_mushy`) -* **`none`**: No explicit treatment. -* **`interp`**: Smooth transition between solid and liquid. +See [Solid-Phase](@ref) for the underlying theory. -#### Fluid Model (`module_fluid`) -* **`none`**: No fluid layer. -* **`fluid0d`**: Bulk fluid approximation. -* **`fluid1d`**: Radially structured fluid with heating. +* **`"none"`**: No solid-phase tidal model. +* **`"solid0d"`**: Homogeneous solid approximation. The solid region is treated as a single effective layer with averaged mechanical properties; fast, but cannot resolve radial structure. +* **`"solid1d"`**: Radially resolved structure, solved with the shooting method. Allows realistic rigidity, density, and rheology profiles. +* **`"solid1d-relax"`**: Same as `solid1d`, but uses the relaxation method instead of shooting. Generally more stable, but slightly slower. +* **`"solid1d-mush"`**: Same as `solid1d`, but accounts for a partially molten/porous ("mushy") (interface) layer. +* **`"solid1d-mush-relax"`**: Same as `solid1d-relax`, but accounts for partially molten/porous ("mushy") regions. Uses the relaxation method instead of shooting. +* **`"solid1d-equil-relax"`**: Same as `solid1d-relax`, but limited to the fully fluid interior case (equilibrium tide). Automatically selected by the other `solid1d*` models when the forcing frequency is below the inverse Hubble time. + +```@raw html +

models [orbit.obliqua.module_mushy]

+``` + +See [Mush layer - interp](@ref) for the underlying model. + +* **`"none"`**: No explicit mushy layer treatment. +* **`"interp"`**: Tidal dissipation transitions smoothly between solid and fluid regimes using interpolation across the melt fraction range. Accounts for imaginary part of the Lovenumber, but not the real part. + +```@raw html +

models [orbit.obliqua.module_fluid]

+``` + +See [Liquid-Phase](@ref) for the underlying model. + +* **`"none"`**: No fluid-phase tidal model. +* **`"fluid0d"`**: Bulk fluid approximation. The fluid region is treated as a single effective layer with averaged mechanical properties; fast, but cannot resolve radial structure. +* **`"fluid1d"`**: Same as `fluid0d`, but allows a user-specified radial heating distribution (see `[orbit.obliqua.fluid].sigma_R_prf`). --- -### Solid Interior Parameters (`[orbit.obliqua.solid]`) -* **`ncalc`**: Number of radial layers (shooting method). -* **`dr_min`**, **`dr_max`**: Grid spacing for relaxation solver [m]. -* **`core`**: Core boundary condition (`"liquid"`, `"solid"`, `"inertial"`). -* **`core_props`**: Core properties (shear modulus, bulk modulus) to use for CMB boundary condition (`"core"`, `"mantle"`). -* **`inertial_terms`**: Boolean flag to include inertial terms in the motion matrix. -* **`bulk_l`**: Liquid bulk modulus [Pa]. -* **`dbulk_power`**: Drained bulk modulus powerlaw scaling exponent. -* **`porosity_thresh`**: Threshold below which no mush is formed. +### Solid Interior Parameters + +```@raw html +

config [orbit.obliqua.solid]

+``` + +```@raw html +
+``` + +| NAME | TYPE | DESCRIPTION | +| :--- | :--- | :--- | +| `ncalc` | int | Number of sublayers for `solid1d` solvers using the shooting method. | +| `dr_min` | float | Minimum spacing between grid points, relaxation method [m]. | +| `dr_max` | float | Maximum spacing between grid points, relaxation method [m]. | +| `core` | str | Core solution used as the CMB boundary condition: `"liquid"`, `"solid"`, `"inertial-liquid"`, or `"inertial"`. See [Solid-Phase](@ref) ("The four `core` options"). | +| `core_props` | str | Core properties (shear modulus, bulk modulus) to use for the CMB boundary condition: `"core"` or `"mantle"`. | +| `inertial_terms` | bool | Include inertial terms in the solid tidal response equations. | +| `bulk_l` | float | Liquid bulk modulus [Pa]. | +| `dbulk_power` | float | Drained bulk modulus power-law scaling exponent. | +| `porosity_thresh` | float | Percolation threshold, at which there is a first-order transition from fully connected to fully isolated pore space [dimensionless]. | + +```@raw html +
+``` --- -### Fluid Parameters (`[orbit.obliqua.fluid]`) -* **`sigma_R`**: Rayleigh drag at the interface. -* **`sigma_R_inf`**: Drag in the bulk fluid. -* **`sigma_R_prf`**: Vertical drag profile. -* **`H_R`**: Scale height [m]. -* **`efficiency`**: Drag efficiency factor. +### Fluid Parameters + +```@raw html +

config [orbit.obliqua.fluid]

+``` + +```@raw html +
+``` + +| NAME | TYPE | DESCRIPTION | +| :--- | :--- | :--- | +| `sigma_R` | float | Rayleigh drag coefficient at the interface [dimensionless]. | +| `sigma_R_inf` | float | Rayleigh drag coefficient in the bulk (pure) fluid [dimensionless]. | +| `sigma_R_prf` | str | Vertical drag profile: `"uniform"`, `"exp"`, `"linear"`, `"quadratic"`, `"dynamic"`, or `"dynamic_interp"` (default). See [Liquid-Phase](@ref) for the physical meaning of each. | +| `H_R` | float | Rayleigh drag scale height [m]. | +| `efficiency` | float | Rayleigh drag efficiency at the core interface [dimensionless]. | + +```@raw html +
+``` --- -### Mushy Layer Parameters (`[orbit.obliqua.mushy]`) -* **`b_width`**, **`t_width`**: Width of dissipation peaks as fraction of layer thickness. +### Mushy Layer Parameters + +```@raw html +

config [orbit.obliqua.mushy]

+``` + +```@raw html +
+``` + +| NAME | TYPE | DESCRIPTION | +| :--- | :--- | :--- | +| `b_width` | float | Width of the bottom dissipation peak, as a fraction of layer thickness. | +| `t_width` | float | Width of the top dissipation peak, as a fraction of layer thickness. | + +```@raw html +
+``` --- -### Planetary Structure (`[struct]`) -Defines bulk planetary properties: -* **`core_density`**: Core density [kg m$^{-3}$]. -* **`core_shear`**: Core shear [Pa s]. -* **`core_bulk`**: Core bulk modulus [Pa ]. +### Planetary Structure + +```@raw html +

config [struct]

+``` + +```@raw html +
+``` + +| NAME | TYPE | DESCRIPTION | +| :--- | :--- | :--- | +| `core_density` | float | Core density [kg m$^{-3}$]. | +| `core_shear` | float | Core shear modulus [Pa]. | +| `core_bulk` | float | Core bulk modulus [Pa]. | + +```@raw html +
+``` --- -### Interior Energetics (`[interior_energetics]`) -Defines interior energy transfer dynamics related properties: -* **`grain_size`**: Grain size [m]. \ No newline at end of file +### Interior Energetics + +```@raw html +

config [interior_energetics]

+``` + +```@raw html +
+``` + +| NAME | TYPE | DESCRIPTION | +| :--- | :--- | :--- | +| `grain_size` | float | Grain size [m]. | + +```@raw html +
+``` diff --git a/docs/src/reference/forcing-frequency.md b/docs/src/reference/forcing-frequency.md index 777a3b4..c7c5f04 100644 --- a/docs/src/reference/forcing-frequency.md +++ b/docs/src/reference/forcing-frequency.md @@ -9,19 +9,26 @@ Given the fact that the tidal forcing magnitude decreases exponentially with har \sigma = m\Omega - k n_{\mathrm{orb}}, ``` -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 +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^{-(n+1), m}_k(e) = +|X^{-(n+1), m}_k(e) |= \bigg| \frac{1}{2\pi} \int_0^{2\pi} \left(\frac{r}{a}\right)^n -e^{im\Omega - ikn_{\mathrm{orb}}}\,dn_{\mathrm{orb}} \geq 0.01, +e^{im\Omega - ikn_{\mathrm{orb}}}\,dn_{\mathrm{orb}} \bigg| \geq 0.001, ``` -(i.e. ~1% corrections). This implies that we are considering the following set $K$ of harmonnic frequencies $k$: +(i.e. ~0.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\} +K = \{k \in \mathbb{Z} \, \big| \, |X^{-(n+1), m}_k| \geq 0.001 \} ``` ---- \ No newline at end of file +--- + +### Function Documentation + +```@docs +Obliqua.Hansen.get_k_range +Obliqua.Hansen.get_hansen +``` \ No newline at end of file diff --git a/docs/src/reference/index.md b/docs/src/reference/index.md index 4c214da..c022a7c 100644 --- a/docs/src/reference/index.md +++ b/docs/src/reference/index.md @@ -1,16 +1,67 @@ ### Reference -In this component you will be guided through the employed formalism and theoretical background of the `Obliqua` package. +This is the theory behind `Obliqua`: where the forcing comes from, how a +material responds to it, and how that response is actually solved for +across a planet's interior. Read it start to finish for the full story, +or jump straight to the model you're configuring — each page stands on +its own and links back to the ones it builds on. -- [Rheology](@ref) -- [Solid-Phase](@ref) -- [Mush layer - interp](@ref) -- [Liquid-Phase](@ref) -- [Forcing Frequency](@ref) -- [Surface Loading](@ref) -- [Tidal potentials](@ref) +#### 1. Setting the stage: what's forcing the tides? +A planet in an eccentric orbit feels tidal forcing at more than just one +frequency. Before any material response can be computed, we first need +to know *which* frequencies matter and *how strongly* each one drives +the response. +- [Forcing Frequency](@ref) — which harmonics and Fourier modes actually + contribute, and why we can usually stop at $n=2$. +- [Tidal potentials](@ref) — normalizing the forcing itself, and turning + a Love-number spectrum into a heating rate. +#### 2. How a material responds +- [Rheology](@ref) — Maxwell, Andrade, and the purely elastic limit: the + three ways `Obliqua` lets shear and bulk moduli become + frequency-dependent (and dissipative). +#### 3. The solid interior + +The bulk of `Obliqua`'s machinery lives here: propagating the tidal +response equations through a radially-structured solid mantle, from a +homogeneous zeroth-order estimate up to a fully poro-viscoelastic, +core-inertia-aware solver. + +- [Solid-Phase](@ref) — the shared theory: $y$-functions, the motion + matrix $\pmb{A}_n(r)$, the core boundary matrix $\pmb{I}_C$ and its + four `core` options, and the shooting vs. relaxation split. +- [Solid-Phase - solid0d](@ref) — the zero-dimensional, single-layer + approximation. Fast, and a good sanity check for everything else. +- [Solid-Phase - solid1d](@ref) — radially resolved, shooting method. +- [Solid-Phase - solid1d-relax](@ref) — the same physics, but solved + with a Henyey-style relaxation scheme instead — generally far more + stable. +- [Solid-Phase - solid1d-mush](@ref) and + [Solid-Phase - solid1d-mush-relax](@ref) — their poro-viscoelastic + extensions, for a mantle with a partially molten layer. +- [Solid-Phase (Equilibrium) - solid1d-equil-relax](@ref) — the reduced, + two-component equilibrium-tide limit used automatically at very low + forcing frequencies. + +#### 4. The mushy transition + +- [Mush layer - interp](@ref) — a lightweight, purely interpolated + treatment of dissipation across a thin mushy or transitional region, + for when the full poro-viscoelastic machinery above is overkill. + +#### 5. The fluid response + +- [Liquid-Phase](@ref) — the Laplace tidal equations for a magma ocean + or other fully fluid layer, and the family of radial dissipation + profiles (including the energy-conserving `dynamic_interp` default) + used to distribute that heating with depth. + +#### 6. Reading out the answer + +- [Surface Loading](@ref) — turning surface boundary conditions into + tidal and load Love numbers, and how the solid and fluid contributions + are currently combined into one global $k_n$. diff --git a/docs/src/reference/liquid-phase.md b/docs/src/reference/liquid-phase.md index 5a94b89..4084a87 100644 --- a/docs/src/reference/liquid-phase.md +++ b/docs/src/reference/liquid-phase.md @@ -87,19 +87,39 @@ so that the dissipation scale adjusts with distance from the boundary. This mimics mixing-length arguments commonly used in geophysical and astrophysical fluid dynamics, where turbulence intensity depends on the available eddy size. +### Dynamic-interpolated dissipation (`dynamic_interp`) + +The `dynamic_interp` profile (the default `sigma_R_prf`) replaces the single depth-dependent shape function above with a composite of three physically distinct contributions, each assigned its own share of the total dissipated power so that the three shares sum exactly back to the total: + +$$P(z) = P_{\mathrm{sbd}}(z) + P_{\mathrm{drag}}(z) + P_{\mathrm{fric}}(z), \qquad +\int P\,dV = E_{\mathrm{sbd}} + E_{\mathrm{drag}} + E_{\mathrm{fric}} = E_{\mathrm{total}}.$$ + +* **Shear/bulk/Darcy component** ($E_{\mathrm{sbd}}$): the shear-, bulk-, and Darcy-heating already computed at the interface (e.g. from a neighbouring mush segment) is carried into the fluid layer with a Gaussian-shoulder radial decay from the top interface, $\propto \exp[-(z/H_{\mathrm{decay}})^{1.5}]$. The decay length starts at the configured scale height $H_R$, but is reduced (found by bisection) if needed so that this component never exceeds 25% of the frequency-independent baseline energy $E_{\infty}$ (the energy dissipated in the ``\sigma_R \to \sigma_{R,\infty}`` bulk-fluid limit). +* **Drag component** ($E_{\mathrm{drag}} = \max(E_{\infty} - E_{\mathrm{sbd}}, 0)$): the remaining bulk Rayleigh-drag dissipation, distributed with a smooth (``\tanh``) sigmoid that switches on around the depth $z_{\mathrm{visc}}$ where the viscosity profile first drops to the pure-liquid viscosity `visc_l`, i.e. at the base of the mush-to-liquid transition. +* **Friction component** ($E_{\mathrm{fric}} = \max(E_{\mathrm{total}} - E_{\infty}, 0)$): the excess dissipation associated with the interfacial (frequency-dependent) drag term, represented as a Gaussian pulse centred below the interface (at $z = 2H_R$, width $H_R/4$) and forced to vanish at $z=0$. + +This construction is designed to interpolate smoothly between the mush-transition heating handed off from the solid side and the bulk fluid response, without requiring the user to hand-pick a single idealized shape. + ### Choosing a dissipation profile Since the physical location of tidal energy dissipation is uncertain in many systems, these profiles should be viewed as parameterized hypotheses. Comparing results across multiple profiles is often more informative than adopting a single choice. In practice: -| Profile | Physical assumption | -| ----------- | --------------------------------------------- | -| Uniform | Energy dissipated evenly throughout the layer | -| Exponential | Strong boundary-layer dissipation | -| Linear | Mild bottom-enhanced dissipation | -| Quadratic | Strongly localized bottom dissipation | -| Dynamic | Mixing-length controlled turbulence | +| Profile | Physical assumption | +| ---------------- | ----------------------------------------------------------------- | +| Uniform | Energy dissipated evenly throughout the layer | +| Exponential | Strong boundary-layer dissipation | +| Linear | Mild bottom-enhanced dissipation | +| Quadratic | Strongly localized bottom dissipation | +| Dynamic | Mixing-length controlled turbulence | +| Dynamic-interpolated (`dynamic_interp`, default) | Energy-budget-conserving split between mush hand-off, bulk drag, and interfacial friction | + +--- +### Function Documentation ---- \ No newline at end of file +```@docs +Obliqua.run_fluid1d +Obliqua.fluid0d.compute_fluid_lovenumbers +``` \ No newline at end of file diff --git a/docs/src/reference/rheology.md b/docs/src/reference/rheology.md index 8d88091..19524fe 100644 --- a/docs/src/reference/rheology.md +++ b/docs/src/reference/rheology.md @@ -27,4 +27,18 @@ The same equations hold for the (drained) bulk modulus, just replace ``\mu`` wit \zeta (\phi) \approx \frac{\eta (\phi)}{\phi}. ``` ---- \ No newline at end of file +A third, purely elastic option is also available (`"elastic"`, equivalent to `"none"`), which switches off dissipation entirely by returning the real, frequency-independent modulus + +```math +\tilde{\mu}(\omega) = \mu. +``` + +This is useful as a dissipation-free limit for testing and for isolating the elastic part of the tidal response. + +--- + +### Function Documentation + +```@docs +Obliqua.complex_modulus +``` \ No newline at end of file diff --git a/docs/src/reference/solid-phase.md b/docs/src/reference/solid-phase.md index 4c55be0..ef7f17a 100644 --- a/docs/src/reference/solid-phase.md +++ b/docs/src/reference/solid-phase.md @@ -80,22 +80,34 @@ $$\begin{aligned} & U & \zeta_n & \tau & P \\ \hline y_{3}(a) & 0 & -g_e \zeta_n & 0 & -P_n \\ y_{4}(a) & 0 & 0 & \tau_n & 0 \\ -\dfrac{n+1}{a} y_{5}(a) + y_{6}(a) +y_{6}(a) & \dfrac{2n+1}{a} U_n & 4\pi G \zeta_n & 0 & 0 \end{array} \end{aligned}$$ -where $g_e$ is the norm of surface gravity. The surface -mass load can also be written as an external potential $U'$ such that -$\zeta_n = [(2n + 1)/4 \pi G a] U'_n$. +where $g_e$ is the norm of surface gravity. No $y_5$ term appears in the +third row because $y_6$ is already Takeuchi & Saito's (1972) combined +"potential stress" variable: its own radial equation, $dy_6/dr = +(n-1)y_6/r + \ldots$, has no coupling to $y_5$ (unlike the raw +potential-gradient definition used in some other conventions). The +surface mass load can also be written as an external potential $U'$ such +that $\zeta_n = [(2n + 1)/4 \pi G a] U'_n$. $$\begin{aligned} y_{3}(R) &= - \frac{(2n + 1)g_e}{4 \pi G R} U'_n - P_n \\ y_{4}(R) &= \tau_n \\ -\dfrac{n+1}{R} y_{5}(R) + y_{6}(R) +y_{6}(R) &= \frac{2n+1}{R} (U_n + U'_n) \end{aligned}$$ +Note that `get_surface_bc!` in `src/common.jl` does not literally apply +this $\zeta_n \leftrightarrow U'_n$ conversion; it instead sets +$(U,U',\tau,P)$ directly as dimensionless $0$/$1$ selector flags (see +below), which for the load case numerically works out to $y_3(a) = +-(2n+1)g(a)/(4\pi a^2)$ and $y_6(a) = (2n+1)G/a^2$ — see +[Solid-Phase - solid1d](@ref) for the concrete tidal/load values the +code actually produces. + For now we only use the following: the tidal Love number corresponding to the case of an external potential perturbation $U$ and the load Love number computed for a surface mass load perturbation @@ -171,14 +183,14 @@ $$\small \dfrac{4}{r}\!\left( \dfrac{3\kappa\mu}{r\beta} - \rho_{0}g \right) { - \rho_0 \omega^2} & \dfrac{\ell(\ell+1)}{r}\!\left(\rho_{0}g - \dfrac{6\kappa\mu}{r\beta}\right) & --\dfrac{4\mu}{r\beta} & \dfrac{\ell(\ell+1)}{r} & - \dfrac{\rho_{0}(\ell+1)}{r} & -\rho_{0} \\[1.2em] +-\dfrac{4\mu}{r\beta} & \dfrac{\ell(\ell+1)}{r} & \dfrac{\rho_{0}(\ell+1)}{r} & +-\rho_{0} \\[1.2em] \dfrac{1}{r}\!\left(\rho_{0}g - \dfrac{6\mu\kappa}{r\beta}\right) & \dfrac{2\mu}{r^{2}}\!\left[\ell(\ell+1)\!\left(1+\dfrac{\lambda}{\beta}\right)-1\right] { - \rho_0 \omega^2} & -\dfrac{\lambda}{r\beta} & - \dfrac{3}{r} & -\dfrac{\rho_{0}}{r} & 0 \\[1.2em] --4\pi G \rho_{0} & 0 & 0 & 0 & -\dfrac{\ell+1}{r} & 1 \\[1.2em] --\dfrac{4\pi G \rho_{0} (\ell+1)}{r} & \dfrac{4\pi G \rho_{0} \ell(\ell+1)}{r} & 0 & 0 & 0 & \dfrac{\ell-1}{r} +-\dfrac{\rho_{0}}{r} & 0 \\[1.2em] +4\pi G \rho_{0} & 0 & 0 & 0 & -\dfrac{\ell+1}{r} & 1 \\[1.2em] +\dfrac{4\pi G \rho_{0} (\ell+1)}{r} & -\dfrac{4\pi G \rho_{0} \ell(\ell+1)}{r} & 0 & 0 & 0 & \dfrac{\ell-1}{r} \end{pmatrix}$$ with ($\ell = n$, same thing, different notation) and @@ -195,37 +207,41 @@ As you can see, this matrix has size $6\times 6$, implying that it couples six $y$-functions. Specifically these are $U_{n,m} , V_{n,m} , X_{n,m} , Y_{n,m}, \Phi_{n,m} , \Psi_{n,m}$. We may extend the motion matrix to also include corrections from porosity, the -corresponding a matrix will have size $8\times 8$ +corresponding a matrix will have size $8\times 8$. A dot ($\cdot$) marks +an entry that is carried over unchanged from the $6\times 6$ matrix +above (whether or not that entry happens to be zero); an explicit "$0$" +means the entry is genuinely zero, or genuinely new and zero, in the +poro-elastic system. -$$\tiny +$$\scriptsize A = \begin{pmatrix} -\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \alpha \beta^{-1} & 0 \\ -\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 0 & 0 \\ -1i k \rho_\ell^2 g^2 n(n+1)\! \dfrac{r^{-2}}{\omega \eta_\ell} & \cdot & \cdot & \cdot & --\! (n+1) r^{-1} 1i k \rho_\ell^2 g n \dfrac{r^{-1}}{\omega \eta_\ell} & -\cdot & 1i k \rho_\ell g n(n+1)\!\dfrac{r^{-2}}{\omega \eta_\ell} - 4\mu \alpha \beta^{-1} r^{-1} & 1i k \rho_\ell^2 g^2 n(n+1)\!\dfrac{r^{-2}}{\omega \eta_\ell} - 4 \phi \rho_\ell g r^{-1} \\ -0 & 0 & 0 & 0 & 0 & 0 & 2\alpha\mu r^{-1}\beta^{-1} & \phi \rho_\ell g r^{-1} \\ -0 & 0 & 0 & 0 & 0 & 0 & 0 & 4\pi G \rho_\ell \phi \\ --1i 4\pi G n(n+1) r^{-1} k\rho_\ell^2 g\dfrac{r^{-1}}{\omega\eta_\ell} & 0 & 0 & 0 & -1i 4\pi n(n+1)G\rho_\ell^2 k\dfrac{r^{-2}}{\omega\eta_\ell} & 0 & --1i 4\pi n(n+1)G\rho_\ell k\dfrac{r^{-2}}{\omega\eta_\ell} & -4\pi G (n+1)r^{-1}\!\left(\phi\rho_\ell - 1i n k\rho_\ell^2 g\dfrac{r^{-1}}{\omega\eta_\ell}\right) \\ -\rho_\ell g r^{-1}\!\left(4 - 1i k\rho_\ell g n(n+1)\dfrac{r^{-1}}{\omega \phi \eta_\ell}\right) & +\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \alpha \beta^{-1} & 0 \\[0.6em] +\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 0 & 0 \\[0.6em] +i\, k \rho_\ell^2 g^2 n(n+1) \dfrac{r^{-2}}{\omega \eta_\ell} & \cdot & \cdot & \cdot & +-(n+1) r^{-1} \cdot i\, k \rho_\ell^2 g n \dfrac{r^{-1}}{\omega \eta_\ell} & +\cdot & i\, k \rho_\ell g n(n+1) \dfrac{r^{-2}}{\omega \eta_\ell} - 4\mu \alpha \beta^{-1} r^{-1} & i\, k \rho_\ell^2 g^2 n(n+1) \dfrac{r^{-2}}{\omega \eta_\ell} - 4 \phi \rho_\ell g r^{-1} \\[0.6em] +\cdot & \cdot & \cdot & \cdot & \cdot & 0 & 2\alpha\mu r^{-1}\beta^{-1} & \phi \rho_\ell g r^{-1} \\[0.6em] +\cdot & 0 & 0 & 0 & \cdot & \cdot & 0 & 4\pi G \rho_\ell \phi \\[0.6em] +-i\, 4\pi G n(n+1) r^{-1} k\rho_\ell^2 g\dfrac{r^{-1}}{\omega\eta_\ell} & \cdot & 0 & 0 & +i\, 4\pi n(n+1)G\rho_\ell^2 k\dfrac{r^{-2}}{\omega\eta_\ell} & \cdot & +-i\, 4\pi n(n+1)G\rho_\ell k\dfrac{r^{-2}}{\omega\eta_\ell} & +4\pi G (n+1)r^{-1}\left(\phi\rho_\ell - i\, n k\rho_\ell^2 g\dfrac{r^{-1}}{\omega\eta_\ell}\right) \\[0.6em] +\rho_\ell g r^{-1}\left(4 - i\, k\rho_\ell g n(n+1)\dfrac{r^{-1}}{\omega \phi \eta_\ell}\right) & -\rho_\ell n(n+1) g r^{-1} & 0 & 0 & --\rho_\ell (n+1)r^{-1}\!\left(1 - 1i k\rho_\ell g n\dfrac{r^{-1}}{\omega \phi \eta_\ell}\right) & +-\rho_\ell (n+1)r^{-1}\left(1 - i\, k\rho_\ell g n\dfrac{r^{-1}}{\omega \phi \eta_\ell}\right) & \rho_\ell & --1i k\rho_\ell g n(n+1)\!\dfrac{r^{-2}}{\omega \phi \eta_\ell} & --1i \omega \phi \eta_\ell / k - 4\pi G(\rho - \phi \rho_\ell)\rho_\ell + \rho_\ell g r^{-1}\!\left(4 - 1i k\rho_\ell g n(n+1)\dfrac{r^{-1}}{\omega \phi \eta_\ell}\right) \\ -r^{-1}\!\left(1i k\rho_\ell g n(n+1)\dfrac{r^{-1}}{\omega \phi \eta_\ell} - \dfrac{\alpha}{\phi} 4\mu\beta^{-1}\right) & +-i\, k\rho_\ell g n(n+1) \dfrac{r^{-2}}{\omega \phi \eta_\ell} & +-i\, \omega \phi \eta_\ell / k - 4\pi G(\rho - \phi \rho_\ell)\rho_\ell + \rho_\ell g r^{-1}\left(4 - i\, k\rho_\ell g n(n+1)\dfrac{r^{-1}}{\omega \phi \eta_\ell}\right) \\[0.6em] +r^{-1}\left(i\, k\rho_\ell g n(n+1)\dfrac{r^{-1}}{\omega \phi \eta_\ell} - \dfrac{\alpha}{\phi} 4\mu\beta^{-1}\right) & \dfrac{\alpha}{\phi} 2n(n+1)\mu \beta^{-1} r^{-1} & -\dfrac{\alpha}{\phi}\beta^{-1} & 0 & --1i k \rho_\ell n(n+1)\!\dfrac{r^{-2}}{\omega \phi \eta_\ell} & +-i\, k \rho_\ell n(n+1) \dfrac{r^{-2}}{\omega \phi \eta_\ell} & 0 & -1i k n(n+1)\!\dfrac{r^{-2}}{\omega \phi \eta_\ell} - \dfrac{1}{\phi}(S + \alpha^2 \beta^{-1}) & -1i k \rho_\ell g n(n+1)\!\dfrac{r^{-2}}{\omega \phi \eta_\ell} - 2r^{-1} +i\, k n(n+1) \dfrac{r^{-2}}{\omega \phi \eta_\ell} - \dfrac{1}{\phi}(S + \alpha^2 \beta^{-1}) & +i\, k \rho_\ell g n(n+1) \dfrac{r^{-2}}{\omega \phi \eta_\ell} - 2r^{-1} \end{pmatrix}$$ where @@ -243,10 +259,14 @@ formulation. It is important to note that the ordering of the elements in $\pmb{A}_n(r)$ and $\pmb{y}_{n,m}(r)$ must be consistent. Different papers will order the $y$-functions differently (specifically $y_2$ and -$y_3$ are often flipped). We will also perform a reorganization, where -we order $\pmb{y}_{n,m}(r)$ as $[\text{lower}, \text{upper}]^T$. -However, for now we can omit this as there is currently no benefit to -doing this yet. +$y_3$ are often flipped). For the shooting-method solvers (`solid1d`, +`solid1d-mush`) the standard ordering above is used directly, with no +benefit to reordering. The relaxation-method solvers do perform the +$[\text{lower}, \text{upper}]^T$ reorganization described here: this is +exactly the `Y6 = [1,2,4,5,3,6]` and `Y8 = [1,2,5,6,3,7,4,8]` ordering +used in `src/solid1d_mush_relax.jl` to embed the $6\times6$ elastic +system as a view into the $8\times8$ poro-elastic one — see +[Solid-Phase - solid1d-mush-relax](@ref). This puts us in a bit of an awkward spot, since we have a coupled system of ODEs, but we can only constrain half of the $y$-functions at the @@ -270,7 +290,10 @@ introducing three new unknowns: we are saved the system is no longer overestimated! Let us define $\pmb{C}$ the vector of constants of integration (weights) -$$\pmb{C} = (a_1, a_2, a_3),$$ these constants of integration must be + +$$\pmb{C} = (a_1, a_2, a_3),$$ + +these constants of integration must be determined using the boundary condition at the Earth's surface for the stress components of the spheroidal vector solution. The (true) modal solution vector can be written as a superposition of the elementary @@ -303,10 +326,10 @@ $$\small \begin{pmatrix} -r^n/g & 0 & 1 \\[1.2em] 0 & 1 & 0 \\[1.2em] -0 & 0 & g/\rho \\[1.2em] +0 & 0 & g \rho \\[1.2em] 0 & 0 & 0 \\[1.2em] -r^n & 0 & 0 \\[1.2em] -2(n-1)r^{(n-1)} & 0 & 4 \pi G \rho +-r^n & 0 & 0 \\[1.2em] +-2(n-1)r^{(n-1)} & 0 & -4 \pi G \rho \end{pmatrix}$$ The three columns correspond to the elementary @@ -320,145 +343,81 @@ inclusion of inertia effects will be essential for the Earth-Moon system, this does not change the method so we can proceed like this for now. -Given the above you may skip this subsection, here I will briefly -introduce $\pmb{I}_C$ for the compressible fluid core with inertial -effects (assuming shear is zero.) We start with the aforementioned -elementary solutions from Equations from Takeuchi & Saito (1972) - -$$\begin{aligned} - y_1 (r) &= n r^{n-1} \\ - y_2 (r) &= r^{n-1} \\ - y_3 (r) &= 0 \\ - y_4 (r) &= 0 \\ - y_5 (r) &= -(n\gamma - \omega^2 ) r^n \\ - y_6 (r) &= -[2(n-1)n\gamma - (2n+1)\omega^2] r^{n-1} -\end{aligned}$$ - -where $\gamma = 4\pi G\rho /3$. The second elementary -solution is more troublesome - -$$\begin{aligned} - y_1 (r) &= - \frac{r^{l+1}}{2n+3} \left[ \frac12 n h \psi_n(x) + f \phi_{n+1} (x) \right] \\ - y_2 (r) &= - \frac{r^{l+1}}{2n+3} \left[ \frac12 h \psi_n(x) - \phi_{n+1} (x) \right] \\ - y_3 (r) &= - \lambda r^n f \phi_n (x) \\ - y_4 (r) &= 0 \\ - y_5 (r) &= - r^{n+2} \left[ \frac{\alpha^2 f}{r^2} - \frac{3\gamma f}{2(2n+3)} \psi_n (x) \right] \\ - y_6 (r) &= - r^{n+1} \left[ \frac{(2n+1) \alpha^2 f}{r^2} - \frac{3\gamma [ (2n+1)f - nh]}{2(2n+3)}\psi_n(x)\right] -\end{aligned}$$ - -where - -$$\begin{aligned} - x &= kr \\ - k^2 &= \frac{1}{\alpha^2} \left[ \omega^2 + 4\gamma - \frac{n(n+1)\gamma^2}{\omega^2} \right] \\ - f &= -\frac{\omega^2}{\gamma} \\ - h &= f - (n+1) \\ - \phi_n (x) &= \frac{(2n+1)!!}{x^n} j_n (x) \\ - \psi_n (x) &= \frac{2(2n+3)}{x^2} [1-\phi_n (x)]. -\end{aligned}$$ - -Here $\alpha$ is the compressional wave speed and -$j_n (x)$ is the spherical Bessel function of the first kind. The issue -is that as we decrease the forcing frequency $\omega$, $k^2$ becomes -negative, and $x$ becomes large and imaginary. The spherical Bessel -function of the first kind grows exponentially, and the elementary -solution rapidly grows beyond numerical precision. One potential fix -would be to use non-dimensional scaling, however this is currently not -working properly. Another way to resolve this is to note that the -elementary solution may be scaled entirely by some arbitrary constant -($C_2$), as noted earlier in the text. As such, we may simply divide by -$j_n(x)$ to remain around unity. Note that without this the solution -diverges and the system cannot be resolved. With this in mind we need -the recursion relations for the spherical Bessel function of the first -kind, define - -$$z_{n} (x) = x j_{n+1} (x) / j_n (x)$$ - -such that - -$$z_{n-1} (x) = \frac{x^2}{(2n+1) - z_n(x)}$$ - -This recursion calculation -should be carried out in order of decreasing $n$ starting from a -sufficiently large wavenumber with an initial value -$z_n(x) = x^2/(2n+3)$. We can now do some algebraic manipulation of the -elementary solution to obtain a form with only terms $\propto 1/j_n(x)$ -or unity. Since, as mentioned, $j_n(x)$ will diverge for large imaginary -$x$, the newly obtained solution will be well-behaved around unity. Let -us define - -$$\begin{aligned} - \bar{\phi}_n (x) &\equiv \frac{\phi_n (x)}{j_n (x)} = \frac{(2n+1)!!}{x^n} \\ - \bar{\phi}_{n+1} (x) &\equiv \frac{\phi_{n+1} (x)}{j_n (x)} = \frac{(2n+3)!!}{x^{n+1}} \frac{j_{n+1}(x)}{j_n (x)} = \frac{(2n+3)!!}{x^{n+2}}z_n(x) \\ - \bar{\psi}_n (x) &= \frac{\psi_n(x)}{j_n(x)} = \frac{2(2n+3)}{x^2} \left[ \frac{1}{j_n (x)} - \bar{\phi}_n (x)\right] -\end{aligned}$$ - -The elementary solution is then written: - -$$\begin{aligned} - \bar{y}_1 (r) &= - \frac{r^{n+1}}{2n+3} \left[ \frac{1}{2} n h \bar{\psi}_n(x) + f \frac{(2n+3)!!}{x^{n+2}} z_n(x) \right] \\ - \bar{y}_2 (r) &= - \frac{r^{n+1}}{2n+3} \left[ \frac{1}{2} h \bar{\psi}_n(x) - \frac{(2n+3)!!}{x^{n+2}} z_n(x) \right] \\ - \bar{y}_3 (r) &= - \lambda r^n f \bar{\phi}_n (x) \\ - \bar{y}_4 (r) &= 0 \\ - \bar{y}_5 (r) &= - r^{n+2} \left[ \frac{\alpha^2 f}{r^2} \frac{1}{j_n(x)} - \frac{3\gamma f}{2(2n+3)} \bar{\psi}_n (x) \right] \\ - \bar{y}_6 (r) &= - r^{n+1} \left[ \frac{(2n+1) \alpha^2 f}{r^2} \frac{1}{j_n(x)} - \frac{3\gamma [ (2n+1)f - nh]}{2(2n+3)}\bar{\psi}_n(x)\right] -\end{aligned}$$ - -Now we note that we can simplify the system of equations -when $k^2 \ll -1$, which implies - -$$\begin{aligned} - \omega^2 + 4\gamma &\ll \frac{n(n+1)\gamma^2}{\omega^2} \\ - \omega^4 + 4\gamma \omega^2 &\ll n(n+1)\gamma^2 -\end{aligned}$$ - -Noting that the LHS is dominated by $4\gamma \omega^2$ as $\omega \to 0$, we have - -$$\begin{aligned} - 4\gamma \omega^2 &\ll n(n+1)\gamma^2 \\ - 4\omega^2 &\ll n(n+1)\gamma \qquad \wedge \qquad \gamma \neq 0 -\end{aligned}$$ - -We can now find a typical forcing frequency for the -Earth's iron core, below which we may further simplify the system. Let's -assume for the outer liquid core -$\rho \approx 10000 \text{ to } 12000 \text{ kg/m}^3$. Then, plugging in -to obtain $\gamma$, we have - -$$\gamma \approx \frac{4}{3} \pi (6.67 \times 10^{-11}) (11000) \approx 3 \times 10^{-6} \text{ s}^{-2}$$ - -For a degree $n=2$, that becomes - -$$\frac{n(n+1)}{4} \gamma = \frac{6}{4} (3 \times 10^{-6}) = 4.5 \times 10^{-6} \text{ s}^{-1}$$ - -So that would imply that $\omega \ll 10^{-3} \text{ s}^{-2}$. In this -regime we should use - -$$\begin{aligned} - \bar{y}_1 (r) &= - \frac{r^{n+1}}{2n+3} \left[ \frac{n h (2n+3)}{x^2} (-\bar{\phi}_n) + f \frac{(2n+3)!!}{x^{n+2}} z_n(x) \right] \\ - \bar{y}_2 (r) &= - \frac{r^{n+1}}{2n+3} \left[ \frac{h (2n+3)}{x^2} (-\bar{\phi}_n) - \frac{(2n+3)!!}{x^{n+2}} z_n(x) \right] \\ - \bar{y}_3 (r) &= - \lambda r^n f \bar{\phi}_n (x) \\ - \bar{y}_4 (r) &= 0 \\ - \bar{y}_5 (r) &= - r^{n+2} \left[ - \frac{3\gamma f}{2(2n+3)} \bar{\psi}_n (x) \right] \\ - \bar{y}_6 (r) &= - r^{n+1} \left[ - \frac{3\gamma [ (2n+1)f - nh]}{2(2n+3)}\bar{\psi}_n(x)\right] -\end{aligned}$$ - -which further simplify to - -$$\begin{aligned} - \bar{y}_1 (r) &= - \frac{r^{n+1}(2n+1)!!}{x^{n+2}} \left[ f z_n(x) - nh \right] \\ - \bar{y}_2 (r) &= - \frac{r^{n+1}(2n+1)!!}{x^{n+2}} \left[ -z_n(x) -h \right] \\ - \bar{y}_3 (r) &= - \frac{\lambda f r^n (2n+1)!!}{x^n} \\ - \bar{y}_4 (r) &= 0 \\ - \bar{y}_5 (r) &= - \frac{3\gamma f r^{n+2} (2n+1)!!}{x^{n+2}} \\ - \bar{y}_6 (r) &= - \frac{3\gamma [(2n+1)f - nh] r^{n+1} (2n+1)!!}{x^{n+2}} -\end{aligned}$$ - -Every term now shares the common factor -$\mathcal{K} = \frac{(2n+1)!!}{x^{n+2}}$. If you wish, you can even -divide the entire solution by $\mathcal{K}$ (since the solution is -scalable by an arbitrary constant $C_2$), which leaves us with a -remarkably clean system: +#### The four `core` options + +The `[orbit.obliqua.solid].core` option selects which $\pmb{I}_C$ is used. There +are two static (frequency-independent) limits and two dynamic +(frequency-dependent, $\omega \neq 0$) generalizations of them: + +* **`"liquid"`**: the static ($\omega = 0$) fluid core matrix $\pmb{I}_C$ given + above (Takeuchi & Saito 1972). Two columns carry the physical degrees of + freedom; the third is a pure tangential-slip vector ($y_2 = 1$, else $0$), + representing the discontinuity in horizontal displacement that a fluid core + permits at the fluid-solid interface. +* **`"solid"`**: the static, incompressible, elastic solid core of Love + (1911). All three columns are genuine elementary solutions, regular at + $r=0$; the incompressible limit ($K \to \infty$) means $K$ does not enter + the formula. There is no fluid-solid interface, so no column is a slip + vector. +* **`"inertial-liquid"`**: the general, oscillating ($\omega \neq 0$), + compressible fluid ($\mu = 0$) core (Takeuchi & Saito 1972; Korenaga 2025). + This is the case derived below. It reduces to `"liquid"` as $\omega \to 0$: + a slowly oscillating fluid approaches hydrostatic equilibrium. +* **`"inertial"`**: the general, oscillating, compressible, finite-$\mu,K$ + solid core (Takeuchi & Saito 1972; Kervazo et al. 2021). It reduces to + `"inertial-liquid"` as $\mu \to 0$ (an oscillating solid without shear + rigidity is an oscillating fluid), and to `"solid"` as $\omega \to 0$ and + $K \to \infty$ (the static, incompressible limit of the general oscillating + solid recovers the classical elastic solution). All three columns are + again genuine elementary solutions with no slip vector. + +Use `"solid"`/`"liquid"` for the static limit, and +`"inertial"`/`"inertial-liquid"` when frequency-dependent (dynamic) effects +at the core boundary matter. The full matrix entries for `"solid"` and +`"inertial"` are implemented in `get_Ic` (`src/common.jl`) and are not +reproduced here; the derivation below covers the `"inertial-liquid"` case in +detail, since it is algebraically the most tractable of the two dynamic +options. + +### Nullspace vs. direct use: shooting vs. relaxation + +The shooting- and relaxation-method solvers use $\pmb{I}_C$ in two +different ways, and it is worth being explicit about the distinction +since it is easy to conflate them. + +* **Shooting method** (`solid1d`, `solid1d-mush`): the three columns of + $\pmb{I}_C$ are used *directly* as the three starting vectors for the + numerical integration, $\pmb{y}^{(i)}(r_C^+) = \pmb{I}_C \pmb{e}_i$. + Since $\pmb{I}_C$ is constructed to already be regular at $r=0$ and + satisfy the core physics by definition, the CMB boundary condition is + automatically satisfied by construction — there is nothing further to + solve for at the core; only the weights $\pmb{C}$ are later fixed by + the surface condition. +* **Relaxation method** (`solid1d-relax`, `solid1d-mush-relax`, + `solid1d-equil-relax`): here the solver never propagates elementary + solutions — it solves for the full state vector $\pmb{y}_n(r)$ directly + at every grid point simultaneously, so the CMB condition must instead + be expressed as a *linear constraint* $\pmb{B}_1 \pmb{y}(r_C^+) = 0$ + that any admissible solution vector satisfies. This $\pmb{B}_1$ is + obtained numerically as the left nullspace of $\pmb{I}_C$ (`get_core_bc!` + in `src/common.jl`: `nullspace(transpose(Ic))`), i.e. the row vectors + orthogonal to every column of $\pmb{I}_C$. Since $\pmb{I}_C$'s three + columns span the admissible 3-dimensional subspace of the 6-dimensional + state space at the CMB, this 3-dimensional orthogonal complement is + exactly equivalent to using $\pmb{I}_C$ as starting vectors — it is the + same boundary condition, just expressed as a constraint rather than a + basis, because the relaxation method never explicitly constructs + $\pmb{C}$. + +The `"inertial-liquid"` and `"inertial"` core options are implemented +exactly this way in `get_Ic` (Takeuchi & Saito 1972; Korenaga 2025; +Kervazo et al. 2021 — see the citations above), using the full, +frequency-dependent Bessel-function solution rather than any low-frequency +approximation. In the low-forcing-frequency limit relevant to most solid +cores ($k^2 \ll -1$, i.e. roughly $\omega \ll 10^{-3}\,\text{s}^{-1}$ for +Earth's core), that general solution admits a considerably simpler, +purely algebraic closed form (no Bessel functions or their recursion +relations required): $$\begin{aligned} \bar{y}_1 (r) &= -r^{n+1} [f z_n(x) - nh] \\ @@ -469,6 +428,12 @@ $$\begin{aligned} \bar{y}_6 (r) &= -3\gamma [(2n+1)f - nh] r^{n+1} \end{aligned}$$ +with $\gamma = 4\pi G\rho/3$, $f = -\omega^2/\gamma$, $h = f-(n+1)$, and +$z_n(x)$ the Bessel ratio $x\,j_{n+1}(x)/j_n(x)$. This low-frequency +simplification is **not currently implemented** — `get_Ic` always +evaluates the general form — and is left as a possible future +optimization rather than derived in full here. + Finally, we should note that the porosity related $y$-functions need also be bounded from below and partially from above. Moreover, since a porous layer is likely to form at some point midway through the mantle, @@ -493,8 +458,17 @@ stable then the shooting method. For now we will finish this section by stating the elementary boundary conditions to be imposed on the poro-viscoelastic solution vector. At the lower boundary of the porous layer the pore pressure is nonzero -$b_{(7,4)} = 1$ and we have zero radial Darcy flux $b_{(8,4)} = 1$. At +$b_{(7,4)} = 1$ and we have zero radial Darcy flux $b_{(8,4)} = 0$. At the upper part of the porous layer we impose again zero radial Darcy -flux $b_{(8)} = 1$, without any constraint on the pore pressure. +flux $b_{(8)} = 0$, without any constraint on the pore pressure. We will now discuss the different solver schemes. + +--- + +### Function Documentation + +```@docs +Obliqua.solid1d.common.get_A +Obliqua.solid1d.common.get_Ic +``` diff --git a/docs/src/reference/solid/solid0d.md b/docs/src/reference/solid/solid0d.md index fcdad4c..562bce9 100644 --- a/docs/src/reference/solid/solid0d.md +++ b/docs/src/reference/solid/solid0d.md @@ -40,4 +40,14 @@ $$k_n^T = \frac{1}{1 + \mu^*_n} \left( \frac{3}{2(n - 1)} \right)$$ #### Load Love Number ($k_n^L$) The load Love number represents the response to a surface mass load: -$$k_n^L = -\frac{1}{1 + \mu^*_n}$$ \ No newline at end of file +$$k_n^L = -\frac{1}{1 + \mu^*_n}$$ + +--- + +### Function Documentation + +```@docs +Obliqua.solid0d.mean_cmu +Obliqua.solid0d.compute_solid_lovenumbers +Obliqua.run_solid0d +``` \ No newline at end of file diff --git a/docs/src/reference/solid/solid1d.md b/docs/src/reference/solid/solid1d.md index 90ed08a..55a90d4 100644 --- a/docs/src/reference/solid/solid1d.md +++ b/docs/src/reference/solid/solid1d.md @@ -8,7 +8,9 @@ The `solid1d` model uses the shooting method approach and is in many regards ide The structure of the model is as follows: 1. We impose the Core-Mantle Boundary (CMB) condition (three elementary solutions). 2. We solve the system at each radius: + $$\frac{d\pmb{y}_{n,m}(r)}{dr} = \pmb{A}_n (r) \pmb{y}_{n,m}(r) - \pmb{f}_{n,m}(r)$$ + stepping (shooting) from the CMB to the surface. For this model, we assume $\pmb{f}_{n,m}(r) = \pmb{0}$, as porosity effects are not included. @@ -20,9 +22,10 @@ $$\pmb{y}_{n,m}(r_{i+1}) = \pmb\Pi_n(r_{i}) \pmb{y}_{n,m}(r_{i})$$ Using the **RK4 (Runge-Kutta 4th Order)** method, we have: -$$\pmb\Pi_n(r_{i}) = \pmb{1} + \frac{1}{6} \left(\pmb{K}_1 + \pmb{K}_2 + \pmb{K}_3 + \pmb{K}_4 \right)$$ +$$\pmb\Pi_n(r_{i}) = \pmb{1} + \frac{1}{6} \left(\pmb{K}_1 + 2\pmb{K}_2 + 2\pmb{K}_3 + \pmb{K}_4 \right)$$ where: + $$\begin{aligned} \pmb{K}_1 &= \Delta r_i \pmb{A}_n (r_i) \\ \pmb{K}_2 &= \Delta r_i \pmb{A}_n (r_i + \Delta r_i / 2) \left[ \pmb{1} + \frac{1}{2} \pmb{K}_1 \right] \\ @@ -42,12 +45,15 @@ It is worth highlighting that this method purely propagates the solution vectors ### Numerical Integration Imposing continuity of the propagator (implying the solution itself is continuous): + $$\pmb \Pi_n (r_i^+,r') = \pmb \Pi_n (r_i^-,r')$$ And imposing CMB conditions in the general solution: + $$\pmb y_{n,m}(r_C^+) = \pmb y_0 = \pmb I_C \pmb C$$ We find: + $$\pmb y_{n,m}(r) = \pmb \Pi_n (r,r_C^+) \pmb I_C \pmb C$$ We iterate through the mantle, repeatedly multiplying the elementary solutions with the propagator matrix until we reach the surface ($R = a^-$). @@ -57,10 +63,22 @@ At the surface, we impose the boundary conditions on the upper components of the $$\pmb P_1 \pmb y (a^-) = \begin{pmatrix} -y_3(a^-) \\[1.2em] -y_4(a^-) \\[1.2em] -\frac{n+1}{a^-} y_5(a^-) + y_6(a^-) +y_3(a^-) \\ +y_4(a^-) \\ +y_6(a^-) \end{pmatrix} = \pmb P_1 \left( \pmb\Pi_n (a^-, r_C^+) \pmb I_C \pmb C \right) = \pmb b$$ -where $\pmb b$ is the 3-vector composed of the RHS of the surface boundary conditions. This allows for the determination of the coefficient vector $\pmb C$, which is used to combine the three elementary solutions into the true modal solution. \ No newline at end of file +where $\pmb b$ is the 3-vector composed of the RHS of the surface boundary conditions. No explicit $y_5$ term is needed in the third row because $y_6$ is already the combined "potential stress" variable of Takeuchi & Saito (1972) — its own radial equation, $dy_6/dr = (n-1)y_6/r + \ldots$, carries no coupling to $y_5$, unlike the raw potential-gradient definition. For unit-amplitude forcing, $\pmb b$ takes the classical free-surface values + +$$\pmb b_{\text{tidal}} = \begin{pmatrix} 0 \\ 0 \\ \dfrac{2n+1}{a} \end{pmatrix}, \qquad \pmb b_{\text{load}} = \begin{pmatrix} -\dfrac{(2n+1)g(a)}{4\pi a^2} \\ 0 \\ \dfrac{(2n+1)G}{a^2} \end{pmatrix}$$ + +for the tidal and load Love-number problems respectively. This allows for the determination of the coefficient vector $\pmb C$, which is used to combine the three elementary solutions into the true modal solution. + +--- + +### Function Documentation + +```@docs +Obliqua.run_solid1d +``` \ No newline at end of file diff --git a/docs/src/reference/solid/solid1d_equil_relax.md b/docs/src/reference/solid/solid1d_equil_relax.md index 7da7939..a7c7b1f 100644 --- a/docs/src/reference/solid/solid1d_equil_relax.md +++ b/docs/src/reference/solid/solid1d_equil_relax.md @@ -55,9 +55,19 @@ The assembled system again has the block-tridiagonal Henyey structure, but with * **Core step** (`core_boundary`): combines $B_1$ with the upper half of $C_1$ to form $S_1$, and initializes $R_1 = -S_1^{-1}Q_1$. * **Propagation step** (`propagate_solid`): for each interior layer, carries forward the "stored" lower half-rows of $C_n$ and $D_{n+1}$ from the previous step to build $P_n$, $S_n$, $Q_n$, then updates -$$X_n = P_n R_{n-1} + S_n, \qquad R_n = -X_n^{-1}Q_n$$ + + $$X_n = P_n R_{n-1} + S_n, \qquad R_n = -X_n^{-1}Q_n$$ + * **Surface step** (`surface_boundary`): applies $B_N$ (separately for the tidal and load cases) in place of the interior $C_N$ upper half, solves $X_N y = b$ for the surface potential pair, and back-fills $y_1$ and the shear-free components. The recursion logic is otherwise identical to the elastic $6\times6$ solver — only the block dimension changes, since the physics being relaxed is restricted to the potential equation rather than the full stress-displacement-potential system. Both a tidal ($y_t$) and a load ($y_l$) surface solution are produced from a single forward sweep, since the two cases share every $C_n$, $D_{n+1}$ block and differ only in the surface right-hand side $b$. +--- + +### Function Documentation + +```@docs +Obliqua.run_solid1d_equil_relax +``` + --- \ No newline at end of file diff --git a/docs/src/reference/solid/solid1d_mush.md b/docs/src/reference/solid/solid1d_mush.md index af08d09..2afa9e2 100644 --- a/docs/src/reference/solid/solid1d_mush.md +++ b/docs/src/reference/solid/solid1d_mush.md @@ -11,10 +11,12 @@ The governing equations are an extension of the `solid1d` model. The additional The solver initially treats the system as a $6 \times 6$ problem. When it encounters a porous layer, it updates the elementary solutions by embedding the $3 \times 1$ coefficient set into a $4 \times 1$ matrix and adding a fourth elementary solution vector, $\pmb{y}^{(4)}_n(r_i)$. At the initial solid-to-porous transition ($r_i$), the new elementary vector is defined as: + $$\pmb{y}^{(4)}_{n,m}(r_i) = (0, 0, 0, 0, 0, 0, 1, 0)^T$$ ### Continuity and Source Terms By turning back to the governing differential equation: + $$\frac{d\pmb{y}_{n,m}(r)}{dr} = \pmb{A}_n (r) \pmb{y}_{n,m}(r) - \pmb{f}_{n,m}(r)$$ We impose continuity across an infinitesimally thin layer ($\Delta r = 0$) where $\pmb\Pi_n(r_{i}) = \pmb{1}$. We allow $\pmb{y}_{n,m}(r_i^+)$ to be an $8 \times 1$ vector and $\pmb{y}_{n,m}(r_i^-)$ to be $6 \times 1$. To preserve continuity while introducing pore pressure, we include a source term $\pmb{f}_{n,m}(r_i^+)$: @@ -22,10 +24,13 @@ We impose continuity across an infinitesimally thin layer ($\Delta r = 0$) where $$\begin{pmatrix} \pmb{1}_{(8\times8)} \end{pmatrix} \pmb{y}_{n,m}(r_i^+) - \pmb{f}_{n,m}(r_i^+) = \begin{pmatrix} \pmb{1}_{(8\times6)} \\ \pmb{0}_{(2\times6)} \end{pmatrix} \pmb{y}_{n,m}(r_i^-)$$ Since the elementary solutions are linearly independent, the modal solution at the interface is: + $$\pmb{y}_n(r_i^+) = a_1 \mathbf{y}^{(1)}(r_i^+) + a_2 \mathbf{y}^{(2)}(r_i^+) + a_3 \mathbf{y}^{(3)}(r_i^+) + a_4 \mathbf{y}^{(4)}(r_i^+)$$ + $$\pmb{y}_{n,m}(r_i^+) = \{ \text{solution of } 6\times 6 \text{ system} \} + (0, 0, 0, 0, 0, 0, a_4, 0)^T$$ To preserve continuity, the source term must exactly cancel the introduced pressure: + $$\pmb{f}_{n,m}(r_i^+) = (0, 0, 0, 0, 0, 0, a_4, 0)^T$$ Understanding these source/sink steps is key for including effects like porosity or seismic activity in the relaxation-based solver. @@ -38,6 +43,7 @@ If the medium becomes solid again before the surface, we reverse the process. We $$\begin{pmatrix} \pmb{1}_{(8\times6)} \\ \pmb{0}_{(2\times6)} \end{pmatrix}^T \pmb{y}_{n,m}(r_i^+) = \begin{pmatrix} \pmb{1}_{(8\times8)} \end{pmatrix} \pmb{y}_{n,m}(r_i^-) - \pmb{f}_{n,m}(r_i^-)$$ We must constrain the **Darcy flux to be zero** at the interface (upper boundary condition). While a sink term can remove excess pore pressure, it cannot "force" a boundary condition like zero Darcy flux if the propagation has already diverged. Instead, we find: + $$\pmb{f}_{n,m}(r_i^-) = (0, 0, 0, 0, 0, 0, y_7(r_i^-), 0)^T$$ A simpler numerical approach is to manually set $y_7 = 0$ for all $r > r_i^-$, as it no longer interacts with the $6 \times 6$ system. @@ -48,4 +54,12 @@ A simpler numerical approach is to manually set $y_7 = 0$ for all $r > r_i^-$, a The remaining steps are identical to the `solid1d` model: 1. Propagate the elementary solution vectors ($8 \times 1$ or $6 \times 1$) to the surface. 2. Determine the coefficients $\pmb{C}$ (now including $a_4$ for the porous component). -3. Apply $\pmb{C}$ to obtain the final modal solution. \ No newline at end of file +3. Apply $\pmb{C}$ to obtain the final modal solution. + +--- + +### Function Documentation + +```@docs +Obliqua.run_solid1d_mush +``` \ No newline at end of file diff --git a/docs/src/reference/solid/solid1d_mush_relax.md b/docs/src/reference/solid/solid1d_mush_relax.md index 479ee59..5c6b73e 100644 --- a/docs/src/reference/solid/solid1d_mush_relax.md +++ b/docs/src/reference/solid/solid1d_mush_relax.md @@ -70,380 +70,66 @@ $$ \pmb{C}_n \pmb{y}_n + \pmb{D}_{n+1} \pmb{y}_{n+1} = \pmb{0}$$ with -$\pmb{A}_n$ the $6 \times 6$ motion matrix. At some point we encounter a -porous layer, at this point we introduce the porous $y$-functions. We -will introduce them in a decoupled fashion and allow them to couple to -the $6 \times 6$ system in a minute. First we need to introduce an -infinitesimal layer over which we impose continuity - -$$\begin{pmatrix} - 1 & & & & & & & \\ - & 1 & & & & & & \\ - & & 1 & & & & & \\ - & & & 1 & & & & \\ - & & & & 1 & & & \\ - & & & & & 1 & & \\ - & & & & & & 1 & \\ - & & & & & & & 1 - \end{pmatrix}_{(8\times8)} - \pmb{y}_{n,m}(r_i^+) + b(r_i^+) = - \begin{pmatrix} - 1 & & & & & \\ - & 1 & & & \\ - & & 1 & & \\ - & & & & & \\ - & & & 1 & \\ - & & & & 1 \\ - & & & & & 1 \\ - & & & & & - \end{pmatrix}_{(8\times6)} - \pmb{y}_{n,m}(r_i^-)$$ - -Evidently, for the LHS to have a zero on the -4th and 8th rows we must have - -$$\begin{aligned} - b_4(r_i^+) &= - y_4(r_i^+) \\ - b_8(r_i^+) &= - y_8(r_i^+) -\end{aligned}$$ - -while continuity of the $6 \times 6$ system implies that -the other components of $b(r_i^+)$ are all zero. One may be inclined to -put $b_4(r_i^+) = -1$ and $b_8(r_i^+) = 0$, i.e. none-zero pore pressure -and zero Darcy flux, but there is a catch! These conditions only hold -for the elementary solution, and not for the modal solution. Instead we -would have to put $b_7(r_i^+) =-a_4$ and $b_8(r_i^+) = 0$, but what is -$a_4$? Since we eliminated all weights already at the core this boundary -is not very useful. Actually, let us try the following - -$$\begin{aligned} - b_4(r_i^+) &= -y_4(r_i^+) \\ - b_8(r_i^+) &= 0 -\end{aligned}$$ - -and treat $y_4(r_i^+)$ as some undetermined quantity. - -This will enter into the equations as follows - -$$\begin{pmatrix} -B_1 \\ -C_1 & D_2 \\ - & C_2 & D_3 \\ - &&\ddots & \ddots \\ - &&&C_{i-1} & D_i \\ - &&&&C_{i}^- & D_i^+ \\ - &&&&& C_{i} & D_{i+1} \\ - &&&&&&\ddots & \ddots \\ - &&&&&&& C_{N-1} & D_N \\ -&&&&&&&&B_N -\end{pmatrix} -\begin{pmatrix} -\pmb{y} (r_C^+) \\ -\pmb{y} (r_2) \\ -\pmb{y} (r_3) \\ -\vdots \\ -\pmb{y} (r_i^-) \\ -\pmb{y} (r_i^+) \\ -\pmb{y} (r_{i+1}) \\ -\vdots \\ -\pmb{y}(r_{N-1}) \\ -\pmb{y}(r_N) -\end{pmatrix} -= \begin{pmatrix} -0 \\ -0 \\ -0 \\ -\vdots \\ -0 \\ -b(r_i^+) \\ -0 \\ -\vdots \\ -0 \\ -b -\end{pmatrix},$$ - -where - -$$C_{i}^- = \begin{pmatrix} - 1 & & & & & \\ - & 1 & & & \\ - & & 1 & & \\ - & & & & & \\ - & & & 1 & \\ - & & & & 1 \\ - & & & & & 1 \\ - & & & & & - \end{pmatrix}_{(8\times6)} - \qquad \text{and} \qquad - D_i^+ = -\begin{pmatrix} - 1 & & & & & & & \\ - & 1 & & & & & & \\ - & & 1 & & & & & \\ - & & & 1 & & & & \\ - & & & & 1 & & & \\ - & & & & & 1 & & \\ - & & & & & & 1 & \\ - & & & & & & & 1 - \end{pmatrix}_{(8\times8)}$$ - -At this point we need to make sure also -include the inhomogeneous part ($b_4(r_i^+)$ and $b_8(r_i^+)$) in our -forwarding scheme. The governing equations change slightly, and we need -to be careful about the dimensions of the matrices involved in the -calculation around the interface. In order to not accidentally remove -information from the systems, we shall embed the $6\times6$ system in -the $8\times8$ space. For some matrix $M$ this implies - -$$M_{(3\times6)} = \begin{bmatrix} - a_1 & b_1 & c_1 & d_1 & e_1 & f_1 \\ - a_2 & b_2 & c_2 & d_2 & e_2 & f_2 \\ - a_3 & b_3 & c_3 & d_3 & e_3 & f_3 \\ - \end{bmatrix}_{(3\times6)} \Rightarrow M_{(4\times8)} = - \begin{bmatrix} - a_1 & b_1 & c_1 & 0 & d_1 & e_1 & f_1 & 0 \\ - a_2 & b_2 & c_2 & 0 & d_2 & e_2 & f_2 & 0 \\ - a_3 & b_3 & c_3 & 0 & d_3 & e_3 & f_3 & 0 \\ - 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ - \end{bmatrix}_{(4\times8)}$$ - -Applying this logic to $C^{l}_{i-1}$ and ${C_i^-}^{u}$ - -$$P_i^- \equiv -\begin{bmatrix} -C^{l}_{i-1} \\ -0 -\end{bmatrix}_{(8\times8)}$$ - -$$S_i^- \equiv -\begin{bmatrix} -D^{l}_i \\ -{C_i^-}^{u} -\end{bmatrix}_{(8\times8)}$$ - -$$Q_i^- \equiv -\begin{bmatrix} -0 \\ -{D_i^+}^{u} -\end{bmatrix}_{(8\times8)}$$ - -We now obtain an $R_i^-$ matrix with dimensions $8 \times 8$. Now we -need to figure out how the $8 \times 8$ system departs from the -interface. Lets start by considering the governing equation for a system -with internal boundary - -$$\left(P_n R_{n-1} + S_n \right) y_n + Q_n y_{n+1} = b_n.$$ - -We can -reformulate this specifically for the interface layer as - -$$\left(P_i^+ R_{i}^- + S_i^+ \right) y_i^+ + Q_i^+ y_{i+1} = b(r_i^+)$$ - -where - -$$P_i^+ \equiv -\begin{bmatrix} -{C_{i}^-}^{l} \\ -0 -\end{bmatrix}$$ - -Again, we need to be careful here given the odd -dimensions of ${C_{i}^-}^{l}$. The continuity of the porous -$y$-functions into their own $2\times2$ system implies that we expand -${C_{i}^-}$ into the $8 \times 8$ identity matrix (as seen from the -porous layer). We thus have $P_i^+$ with dimensions $8 \times 8$. Next - -$$S_i^+ \equiv -\begin{bmatrix} -{D_i^+}^{l} \\ -{C_{i+1}}^{u} -\end{bmatrix}$$ - -is trivial, and similarly - -$$Q_i^+ \equiv -\begin{bmatrix} -0 \\ -{D_{i+1}}^{u} -\end{bmatrix}$$ - -is trivial, both are $8 \times 8$ matrices. - -We can now plug everything in and simplify - -$$\begin{aligned} - \left(P_i^+ R_{i}^- + S_i^+ \right) y_i^+ + Q_i^+ y_{i+1} &= b(r_i^+) -\end{aligned}$$ - -Recall that $b(r_i^+)$ depends on $y_i^+$, we may write $b(r_i^+) = K y_i^+$ where $K$ extracts the 7th and 8th components of $y_i^+$ and changes the sign. We thus write - -$$\begin{aligned} - \left(P_i^+ R_{i}^- + S_i^+ + K\right) y_i^+ &= - Q_i^+ y_{i+1} \\ - y_i^+ &= - \left(P_i^+ R_{i}^- + S_i^+ + K\right)^{-1} Q_i^+ y_{i+1} \\ - &= R_i^+ y_{i+1} -\end{aligned}$$ - -With these conditions on $\pmb{1}$ and $K$ we can -determine the interface form of the equations. - -However, we should be careful not to accidentally store too much -information in $R_i^-$ and $R_i^+$. Realize that we embedded the -$6 \times 6$ system in an $8 \times 8$ space to construct a square and -invertible matrix $X_i^-$. This came at the cost of padding the matrices -with zeros. If one instead combined the non-square matrices, they would -obtain $R_i^-$ with size $6 \times 8$ - -$$\begin{aligned} - P_i^- \cdot R_{i-1} + S_i^- &= X_i^- \\ - (7 \times 6) \cdot (6 \times 6) + (7 \times 6) &\Rightarrow (7 \times 6) \\ - -(X_i^- )^{-1} Q_i^- &= R_i^- \\ - (6 \times 7) \cdot (7 \times 8) &\Rightarrow (6 \times 8) -\end{aligned}$$ - -Hence, at this step the rows corresponding to $y_7$ and -$y_8$ to zero. Next we repeat these steps and include the internal -boundary condition on $y_7$ - -$$\begin{aligned} - P_i^+ \cdot R_i^- + S_i^+ + K &= X_i^+ \\ - (8 \times 6) \cdot (6 \times 8) + (8 \times 8) + (8 \times 8) &\Rightarrow (8 \times 8) \\ - -(X_i^+ )^{-1} Q_i^+ &= R_i^+ \\ - (8 \times 8) \cdot (8 \times 8) &\Rightarrow (8 \times 8) -\end{aligned}$$ - -The resulting system contains now information on $y_8$ -from below, hence we must set the row corresponding to $y_8$ to zero in -$R$. This closes the system and accounts for all dof. - -Now we proceed constructing $R_n$ through the porous layers +$\pmb{A}_n$ the $6 \times 6$ motion matrix, exactly as in +[Solid-Phase - solid1d-relax](@ref). At some point we may encounter a +porous layer, where we must couple the $6\times6$ elastic system to the +$8\times8$ poro-elastic one that carries the two extra porous +$y$-functions (pore pressure $P_{n,m}$ and Darcy flux $R_{n,m}$). + +### Solid-to-mush and mush-to-solid interfaces + +At a transition, the response matrices are still built from the same +$\pmb{C}_n, \pmb{D}_{n+1}$ recursion, but now in the padded $8\times8$ +space so that both sides of the interface can be expressed in a common +$8$-component vector. Two asymmetric adjustments are made to $\pmb C_n$ +or $\pmb D_{n+1}$ depending on the direction of the transition, since one +side of the interface has no porous degrees of freedom to match against: + +* **Solid $\to$ mush** (entering a porous layer from below): the incoming + $3\times6$ "stored" lower half-block from the solid side is scattered + into the $8$-slot ordering at the six non-porous column positions + (i.e. everywhere except the pore-pressure and Darcy-flux slots), and + the pore-pressure row of $\pmb C_n$ is zeroed, since the solid side has + no pore pressure to enforce continuity on. +* **Mush $\to$ solid** (leaving a porous layer from below): both the + pore-pressure and Darcy-flux rows of $\pmb D_{n+1}$ are zeroed, since + the solid side above has neither degree of freedom. + +In both cases a small regularizing term $\pmb K_n$ — a single $1$ placed +on the Darcy-flux diagonal entry — is added to $\pmb X_n = \pmb P_n \pmb +R_{n-1} + \pmb S_n$ before it is inverted, since dropping a row/column +pair can otherwise leave $\pmb X_n$ singular. Because this padding and +regularization can leave $\pmb X_n$ ill-conditioned, the response matrix +is obtained with the Moore-Penrose pseudo-inverse, +$\pmb R_n = -\pmb X_n^{+} \pmb Q_n$, rather than a direct solve. The +physical boundary condition actually enforced at the interface — zero +Darcy flux crossing it — is then applied directly and unconditionally +afterwards, by overwriting the Darcy-flux column (mush $\to$ solid) or +row (solid $\to$ mush) of the resulting $\pmb R_n$ with zeros. This +guarantees the no-flux condition regardless of what the padded, +pseudo-inverted solve produced for that entry. + +Once inside the porous layer, we continue constructing $R_n$ using the +full $8\times8$ motion matrix, $$ - \pmb{C}_n \pmb{y}_n + \pmb{D}_{n+1} \pmb{y}_{n+1} = \pmb{0}$$ with -$\pmb{A}_n$ the $8 \times 8$ motion matrix. - -When we reach the top of the porous layers, potentially at the surface, -but maybe sooner, we need to decouple $y_7$ and $y_8$ from the -$6 \times 6$ system. Impose again an infinitesimal layer - -$$\begin{pmatrix} - 1 & & & & & \\ - & 1 & & & \\ - & & 1 & & \\ - & & & & & \\ - & & & 1 & \\ - & & & & 1 \\ - & & & & & 1 \\ - & & & & & - \end{pmatrix}_{(8\times6)} - \pmb{y}_{n,m}(r_i^+) + b(r_i^+) = - \begin{pmatrix} - 1 & & & & & & & \\ - & 1 & & & & & & \\ - & & 1 & & & & & \\ - & & & 1 & & & & \\ - & & & & 1 & & & \\ - & & & & & 1 & & \\ - & & & & & & 1 & \\ - & & & & & & & 1 - \end{pmatrix}_{(8\times8)} - \pmb{y}_{n,m}(r_i^-)$$ - -Here we apply the constraint that the Darcy -flux vanishes at the interface. The pore pressure remains unconstrained. -Evidently, for the RHS to have a zero on the 8th row we must have - -$$\begin{aligned} - b_7(r_i^+) &= y_7(r_i^-) \\ - b_8(r_i^+) &= 0 -\end{aligned}$$ - -Similar to before, this interface will enter into the equations as -follows - -$$\begin{pmatrix} -B_1 \\ -C_1 & D_2 \\ - & C_2 & D_3 \\ - &&\ddots & \ddots \\ - &&&C_{i-1} & D_i \\ - &&&&C_{i}^- & D_i^+ \\ - &&&&& C_{i} & D_{i+1} \\ - &&&&&&\ddots & \ddots \\ - &&&&&&& C_{N-1} & D_N \\ -&&&&&&&&B_N -\end{pmatrix} -\begin{pmatrix} -\pmb{y} (r_C^+) \\ -\pmb{y} (r_2) \\ -\pmb{y} (r_3) \\ -\vdots \\ -\pmb{y} (r_i^-) \\ -\pmb{y} (r_i^+) \\ -\pmb{y} (r_{i+1}) \\ -\vdots \\ -\pmb{y}(r_{N-1}) \\ -\pmb{y}(r_N) -\end{pmatrix} -= \begin{pmatrix} -0 \\ -0 \\ -0 \\ -\vdots \\ -0 \\ -b(r_i^+) \\ -0 \\ -\vdots \\ -0 \\ -b -\end{pmatrix},$$ - -where - -$$C_{i}^- = \begin{pmatrix} - 1 & & & & & & & \\ - & 1 & & & & & & \\ - & & 1 & & & & & \\ - & & & 1 & & & & \\ - & & & & 1 & & & \\ - & & & & & 1 & & \\ - & & & & & & 1 & \\ - & & & & & & & 1 - \end{pmatrix}_{(8\times8)} - \qquad \text{and} \qquad - D_i^+ = -\begin{pmatrix} - 1 & & & & & \\ - & 1 & & & \\ - & & 1 & & \\ - & & & 1 & \\ - & & & & 1 \\ - & & & & & 1 \\ - & & & & & \\ - & & & & & - \end{pmatrix}_{(8\times6)}$$ + \pmb{C}_n \pmb{y}_n + \pmb{D}_{n+1} \pmb{y}_{n+1} = \pmb{0}, \qquad + \pmb{A}_n \text{ the } 8\times8 \text{ motion matrix},$$ -Again we will pad some of the matrices where necessary +with the same forward/backward Henyey sweep used for the solid segments +above, until the mush layer ends (another interface of the opposite +kind) or the surface is reached. -$$P_i^- \equiv -\begin{bmatrix} -C^{l}_{i-1} \\ -0 -\end{bmatrix}$$ +A simple practical workaround is to assign a porosity slightly above the +percolation threshold (`porosity_thresh`) everywhere a porous layer is +present, so that the solver is built as a single, interface-free +$8\times8$ propagator instead of exercising this coupling logic at all. -$$S_i^- \equiv -\begin{bmatrix} -D^{l}_i \\ -{C_i^-}^{u} -\end{bmatrix}$$ +--- -$$Q_i^- \equiv -\begin{bmatrix} -0 \\ -{D_i^+}^{u} -\end{bmatrix}$$ +### Function Documentation -and apply the conditions on $\pmb{1}$ and $K$ to obtain -the interface form of the equations. (This is still a partial work in -progress, however it seems to work in the numerical implementation.) A -simply work around these interfaces is to assign a porosity slightly -greater than the porosity threshold in the code, this way the solver -will be constructed as purely mush propagator without interfaces. +```@docs +Obliqua.run_solid1d_mush_relax +``` diff --git a/docs/src/reference/solid/solid1d_relax.md b/docs/src/reference/solid/solid1d_relax.md index b37298a..7c39b6a 100644 --- a/docs/src/reference/solid/solid1d_relax.md +++ b/docs/src/reference/solid/solid1d_relax.md @@ -7,13 +7,17 @@ A significant advantage of the `solid1d-relax` model is that it eliminates the n ### Numerical Formulation We begin with the standard problem statement: + $$\frac{d\pmb{y}_{n,m}(r)}{dr} = \pmb{A}_n (r) \pmb{y}_{n,m}(r) - \pmb{f}_{n,m}(r)$$ Assuming no internal sources ($\pmb{f}_{n,m}(r) = \pmb{0}$), we approximate the system using second-order finite differences: + $$\pmb{y}_{n+1} − \pmb{y}_{n} = \frac{\Delta r}{2} (\pmb{A}_{n+1} \pmb{y}_{n+1} + \pmb{A}_n \pmb{y}_n)$$ Rearranging terms yields the fundamental relaxation equation: + $$\pmb{C}_n \pmb{y}_n + \pmb{D}_{n+1} \pmb{y}_{n+1} = \pmb{0}$$ + where: * $\pmb{C}_n = \pmb{I} + \frac{\Delta r}{2} \pmb{A}_n$ * $\pmb{D}_{n+1} = −\pmb{I} + \frac{\Delta r}{2} \pmb{A}_{n+1}$ @@ -25,18 +29,23 @@ Unlike the shooting method, where we track elementary solutions, the relaxation #### Lower Boundary (CMB) At the CMB ($r = r_C^+$), we eliminate unknown coefficients to write the boundary condition as: + $$B_1 \mathbf{y}(r_C^+) = 0$$ + where $B_1$ is a $3 \times 6$ matrix that enforces continuity. By setting the RHS to zero, we ensure the solution remains physically consistent at the interface without needing to manually restart the integration. #### Upper Boundary (Surface) Similarly, at the surface ($r = a^-$), we define the $3 \times 6$ matrix $B_N$: -$$B_N = \begin{pmatrix} 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & (n+1)/a & 1 \end{pmatrix}$$ -This satisfies $B_N \mathbf{y} = b$, where $b$ contains the surface forcing terms. + +$$B_N = \begin{pmatrix} 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix}$$ + +This satisfies $B_N \mathbf{y} = b$, where $b$ contains the surface forcing terms (see [Solid-Phase - solid1d](@ref) for the explicit tidal/load values of $b$). No $y_5$ coefficient appears in the third row: $y_6$ is already the combined potential-stress variable, so matching to the exterior vacuum solution is a direct condition on $y_6(a^-)$ alone. --- ### The Global Matrix and Henyey Relaxation Combining the difference equations and boundary conditions results in a large, banded global system: + $$\mathcal{H} \pmb{y} = \pmb{b}$$ Because $\mathcal{H}$ is typically $6N \times 6N$ (with $N \approx 2000$), direct inversion is too expensive. We use a **Henyey-type relaxation** to solve the system recursively. We define square $6 \times 6$ submatrices $P_n, S_n,$ and $Q_n$: @@ -49,13 +58,18 @@ The system is then solved in two passes: 1. **Forward Sweep (Recursive Relations):** Calculate the response matrices $R_n$ from the CMB to the surface: + $$R_n = -X_n^{-1} Q_n, \quad X_n = P_n R_{n-1} + S_n$$ + starting with $R_1 = -S_1^{-1} Q_1$. 2. **Backward Sweep (Solution):** Determine the surface solution: + $$y_N = X_N^{-1} b$$ + Then propagate the solution back down to the CMB: + $$y_n = R_n y_{n+1}$$ > **Conceptual Analogy:** Imagine a string in a river tied to two endpoints. The river's flow (the physics/differential equations) dictates the string's curve, while the endpoints (boundaries) fix its absolute position. If a rock (internal boundary) is in the river, the string naturally bends around it. @@ -67,4 +81,12 @@ To resolve high gradients near the surface, we use an uneven exponential grid. T $$\alpha = \ln\left(\frac{\Delta r_{\max}}{\Delta r_{\min}}\right)$$ -$$r_i = R_E + (R_c - R_E) \left( \frac{\exp\left( \alpha \frac{N - i}{N - 1} \right) - 1}{\exp(\alpha) - 1} \right), \quad i = 1, \ldots, N$$ \ No newline at end of file +$$r_i = R_E + (R_c - R_E) \left( \frac{\exp\left( \alpha \frac{N - i}{N - 1} \right) - 1}{\exp(\alpha) - 1} \right), \quad i = 1, \ldots, N$$ + +--- + +### Function Documentation + +```@docs +Obliqua.run_solid1d_relax +``` diff --git a/docs/src/reference/surface-loading.md b/docs/src/reference/surface-loading.md index 396cc7c..752448d 100644 --- a/docs/src/reference/surface-loading.md +++ b/docs/src/reference/surface-loading.md @@ -14,16 +14,20 @@ For a given degree $n$, the boundary conditions for the state vector components | :--- | :--- | :--- | :--- | :--- | | **$y_3(R)$** (Normal Stress) | $0$ | $-g_e \zeta_n$ | $0$ | $-P_n$ | | **$y_4(R)$** (Tangential Stress) | $0$ | $0$ | $\tau_n$ | $0$ | -| **$\frac{n+1}{R} y_5(R) + y_6(R)$** | $\frac{2n+1}{R} U_n$ | $4\pi G \zeta_n$ | $0$ | $0$ | +| **$y_6(R)$** (Potential Stress) | $\frac{2n+1}{R} U_n$ | $4\pi G \zeta_n$ | $0$ | $0$ | + +No $y_5$ term appears in the third row: in Obliqua's convention $y_6$ is already Takeuchi & Saito's (1972) combined "potential stress" variable, whose own radial equation has no coupling back to $y_5$, so the surface condition is a direct statement about $y_6(R)$ alone (see [Solid-Phase](@ref) for the underlying motion matrix). By expressing a surface mass load $\zeta_n$ as an equivalent external potential $U'$, where $\zeta_n = \frac{2n + 1}{4 \pi G R} U'_n$, the system simplifies to: $$\begin{aligned} y_{3}(R) &= - \frac{(2n + 1)g_e}{4 \pi G R} U'_n - P_n \\ y_{4}(R) &= \tau_n \\ -\frac{n+1}{R} y_5(R) + y_6(R) &= \frac{2n+1}{R} (U_n + U'_n) +y_6(R) &= \frac{2n+1}{R} (U_n + U'_n) \end{aligned}$$ +Note that Obliqua's actual `get_surface_bc!` (`src/common.jl`) does not literally apply this $\zeta_n \leftrightarrow U'_n$ conversion; it sets $(U,U',\tau,P)$ directly as dimensionless $0$/$1$ selector flags, which for the load case numerically works out to $y_3(R) = -(2n+1)g(R)/(4\pi R^2)$ and $y_6(R) = (2n+1)G/R^2$ — see [Solid-Phase - solid1d](@ref) for the concrete tidal/load values the code actually produces. + ### Calculation of Love Numbers In `Obliqua`, Love numbers are non-dimensionalized by setting the forcing terms to either $1$ (present) or $0$ (absent). * **Tidal Love Number (TLN):** Calculated by setting $(U, U', \tau, P) = (1, 0, 0, 0)$. @@ -54,4 +58,10 @@ In this formulation: While additional corrections for atmospheric pressure or specialized crustal rheologies (Andrade/Maxwell) are not yet active, the framework is designed to incorporate these by calculating the respective pressure and loading Love numbers which are already implemented in the solver. ---- \ No newline at end of file +--- + +### Function Documentation + +```@docs +Obliqua.solid1d.common.get_surface_bc! +``` \ No newline at end of file diff --git a/docs/src/reference/tidal-potentials.md b/docs/src/reference/tidal-potentials.md index 223bbc6..b517717 100644 --- a/docs/src/reference/tidal-potentials.md +++ b/docs/src/reference/tidal-potentials.md @@ -10,7 +10,9 @@ To determine the total tidal heating, `Obliqua` loops over all $(n,m,k)$ pairs a ### Tidal Potential and Normalization For every triplet $(n, m, k)$, we calculate the Hansen coefficient $X_k^{-(n+1),m}(e)$ and the normalization factor $A_{n,m,k}$: -$$A_{n,m,k} = (2 - \delta_{m,0}\delta_{k,0}) (1 - \delta_{m,0}\delta_{k<0}) \sqrt{\frac{4\pi}{2n+1}\frac{(n-m)!}{(n+m)!}} P_n^m(0) X_k^{-(n+1),m}(e)$$ +$$A_{n,m,k} = 2\sqrt{\frac{4\pi}{2n+1}\frac{(n-m)!}{(n+m)!}} P_n^m(0) X_k^{-(n+1),m}(e)$$ + +(When `Obliqua`'s adaptive mode enumeration restricts the loop to $m \geq 0$ to exploit the $\pm m$ symmetry of the problem, it applies additional $m=0$ correction factors to the implementation of this sum so that the $m=0$ terms are not double-counted. Those factors are a bookkeeping detail of that loop-restriction optimization, not part of the definition of $A_{n,m,k}$ itself, and are omitted here.) The associated tidal potential $U_{n,m,k}$ is defined as: @@ -44,4 +46,10 @@ Since the tidal forcing consists of a discrete set of frequencies, the total hea ### Outputs for Orbital Dynamics In addition to the heating rates, `Obliqua` returns the complex tidal Love numbers $k_n(\sigma)$ and their corresponding forcing frequencies $\sigma$. these are required for calculating tidal torques and the long-term orbital evolution of the system. ---- \ No newline at end of file +--- + +### Function Documentation + +```@docs +Obliqua.run_tides +``` \ No newline at end of file From 1e5d55be3be2bafcd4f33ac433ad73824a544d37 Mon Sep 17 00:00:00 2001 From: Marijn Date: Wed, 9 Sep 2026 20:01:24 +0000 Subject: [PATCH 04/11] Added new core solution to config. --- res/config/all_options.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/res/config/all_options.toml b/res/config/all_options.toml index 6a13f90..690a69a 100644 --- a/res/config/all_options.toml +++ b/res/config/all_options.toml @@ -101,7 +101,7 @@ version = "1.0" # Version of this configuration file ncalc = 100 # Number of sublayers to use for dim=1 solid tides (for shooting method) dr_min = 10 # Minimum spacing between grid points (for relaxation method) [m] dr_max = 300 # Maximum spacing between grid points (for relaxation method) [m] - core = "liquid" # Core solution to use as CMB boundary condition, options: "liquid", "solid", "inertial" + core = "liquid" # Core solution to use as CMB boundary condition, options: "liquid", "solid", "inertial-liquid", "inertial" core_props = "core" # Core properties to use for CMB boundary condition, options: "core", "mantle" inertial_terms = true # Include inertial terms in the solid tidal response equations bulk_l = 1e9 # Liquid bulk modulus [Pa] From 226075bd53a25e3ed121cdc666d63b0f45521966 Mon Sep 17 00:00:00 2001 From: Marijn Date: Sat, 12 Sep 2026 20:17:49 +0000 Subject: [PATCH 05/11] Added spectrum=legacy for broader compatibility within PROTEUS. --- docs/src/reference/configuration-file.md | 2 +- res/config/all_options.toml | 6 ++-- src/Obliqua.jl | 43 ++++++++++++++++++++++-- 3 files changed, 46 insertions(+), 5 deletions(-) diff --git a/docs/src/reference/configuration-file.md b/docs/src/reference/configuration-file.md index b0e5804..810c530 100644 --- a/docs/src/reference/configuration-file.md +++ b/docs/src/reference/configuration-file.md @@ -90,7 +90,7 @@ See [Forcing Frequency](@ref) for the underlying model. | :--- | :--- | :--- | | `n` | array | Radial dependence exponent(s) in $(r/a)^n$; since $r \ll a$, only $n=2$ contributes significantly. | | `m` | array | Tidal harmonic(s) of the true anomaly (e.g. $m=2$ semidiurnal, $m=1$ diurnal). | -| `spectrum` | str | Frequency sampling strategy: `"full"` samples the whole k2 spectrum, `"adaptive"` samples only the region of interest. | +| `spectrum` | str | Frequency sampling strategy: `"full"` samples the whole k2 spectrum, `"adaptive"` samples only the region of interest, `"legacy"` reproduces the original LovePy module (hardcoded low-eccentricity, spin-synchronous $(n,m,k) = (2,0,1),(2,2,1),(2,2,3)$ triplet evaluated at a single forcing frequency $\omega$). `"legacy"` overrides `n`, `m`, `s_min`, and `s_max`. | | `N_sigma` | int | Number of probe frequencies to evaluate k2 at (used when `spectrum = "full"`). | | `p_min` | float | Minimum period for orbital and axial frequencies [$\log_{10}$ kyr]. | | `p_max` | float | Maximum period for orbital and axial frequencies [$\log_{10}$ kyr]. | diff --git a/res/config/all_options.toml b/res/config/all_options.toml index 690a69a..30e715c 100644 --- a/res/config/all_options.toml +++ b/res/config/all_options.toml @@ -37,7 +37,7 @@ version = "1.0" # Version of this configuration file eccentricity = 0.0 # initial eccentricity of planet's orbit [dimensionless] [orbit.obliqua] - store_3D = false # Store 3D tidal response (displacement, stress, strain) for each layer. Note: this will generate large output files. If false, radial profiles of tidal responses will be stored instead. + store_3D = false # Store 3D tidal response for each layer. Note: this will generate large output files. If false, radial profiles of tidal responses will be stored instead. enforce_ec = true # Enforce energy conservation in tidal response calculations, improves stability in fluid-mush cases at low forcing frequencies (< 1e-7 Hz). This does not affect the Lovenumbers. optimize_scales = false # Optimize non-dimensionalization scales for the relaxation method. This can improve convergence and stability, especially for low forcing frequencies (< 1e-7 Hz). However, do not use with Bigfloat precision. solid_shell = true # Insert an infinitesimal solid shell around the core. This patches an issue where y2 and y4 become decoupled and cause the solution to diverge in fluid layers. Only use with solid1d_relax or solid1d_mush_relax. @@ -49,7 +49,9 @@ version = "1.0" # Version of this configuration file visc_sus = 5e4 # Solidus viscosity to use [Pa s] n = [2] # Power of the radial factor (goes with (r/a)^{n}, since r<) or region of interest (<"adaptive">) + spectrum = "adaptive"# Obtain k2 spectrum for whole spectrum (<"full">), region of interest (<"adaptive">), or the + # hardcoded low-eccentricity, spin-synchronous (n,m,k) = (2,0,1),(2,2,1),(2,2,3) triplet + # evaluated at a single forcing frequency ω (<"legacy">, matches the original LovePy module) N_sigma = 20 # Number probe frequencies to evaluate k2 at p_min = -10 # Minimum period for orbital and axial frequencies [log(kyr)] i.e. forcing occurion over minimal 1e-17 yr interval p_max = 8 # Maximum period for orbital and axial frequencies [log(kyr)] i.e. forcing occurion over maximal 1e9 yr interval diff --git a/src/Obliqua.jl b/src/Obliqua.jl index 12b6abb..57e6cc2 100644 --- a/src/Obliqua.jl +++ b/src/Obliqua.jl @@ -335,6 +335,11 @@ module Obliqua N_σ = cfg["orbit"]["obliqua"]["N_sigma"] p_min = cfg["orbit"]["obliqua"]["p_min"] p_max = cfg["orbit"]["obliqua"]["p_max"] + elseif spectrum == "legacy" + # No further config keys required: (n, m, k) modes and forcing + # frequency are hardcoded to reproduce the original LovePy module. + else + throw("Invalid spectrum value: $spectrum. Must be 'adaptive', 'full', or 'legacy'.") end material_μ = cfg["orbit"]["obliqua"]["material_mu"] @@ -501,6 +506,38 @@ module Obliqua push!(σ_range, σ_range_i[i]) end @info "Using (n, m, k) = ($(nm[1][1]), $(nm[1][2]), 1) for full spectrum." + + elseif spectrum == "legacy" + # Reproduce the original LovePy module: hardcode the three dominant + # low-eccentricity (n, m, k) modes and force an identical forcing + # frequency across all three, so a single k2 Love number spectrum + # is evaluated (at ω) instead of the full/adaptive mode expansion. + # This assumes spin-orbit synchronisation and e << 1, matching the + # simplifications baked into LovePy; it is not a general-purpose + # replacement for "adaptive". + nmk = [(2, 0, 1), (2, 2, 1), (2, 2, 3)] + + # LovePy hardcodes forcing frequency = orbital mean motion (omega), + # which is only physically exact for spin-orbit synchronous rotation + # (axial == omega): warn loudly if that assumption is violated, + # since this mode silently ignores `axial` otherwise. + if !isapprox(axial, omega; rtol=1e-3) + @warn "Legacy spectrum assumes spin-orbit synchronisation (axial == omega), but axial=$axial rad/s and omega=$omega rad/s differ by more than 0.1%. Forcing frequency is still hardcoded to omega; results will not match a self-consistent tidal calculation for this rotation state." + end + + _, X_02 = Hansen.get_hansen(ecc, 2, 0, 1, 1) + _, X_22 = Hansen.get_hansen(ecc, 2, 2, 1, 3) + X_hansen = [X_02[1], X_22[1], X_22[3]] + + # Note we consistently drop the sign of the forcing frequency, as + # LovePy does, since only the imaginary part of k2 is used. + σ_range = fill(Float64(omega), 3) + + N_σ = length(σ_range) + + @info "Using legacy (LovePy-compatible) spectrum: (n, m, k) = (2, 0, 1), (2, 2, 1), (2, 2, 3), all evaluated at ω = $omega." + else + throw("Invalid spectrum value: $spectrum. Must be 'adaptive', 'full', or 'legacy'.") end # get frequency dependent complex shear modulus per mode @@ -1022,8 +1059,10 @@ module Obliqua # specify mode n_i, m_i, s_i = nmk[iss] - # calculate physical forcing frequency - σ = m_i*axial - s_i*omega + # forcing frequency for this mode (matches `σ_range[iss]` exactly for + # "adaptive", since it was built there with the same m_i*axial - s_i*omega + # formula; for "legacy" this is instead the hardcoded ω shared by all modes) + σ = σ_range[iss] # if forcing frequency is zero, then skip to next frequency (no heating) iszero(σ) && continue From cd9d4f6789df41e5f33b9e01326f01ffa7a0bed7 Mon Sep 17 00:00:00 2001 From: Marijn Date: Sun, 13 Sep 2026 13:41:21 +0000 Subject: [PATCH 06/11] Added LN cap for PROTEUS runs that integrate over extreme resonances. --- res/config/all_options.toml | 1 + src/Obliqua.jl | 65 ++++++++++++++++++++++++++++++++++--- test/test.toml | 1 + 3 files changed, 63 insertions(+), 4 deletions(-) diff --git a/res/config/all_options.toml b/res/config/all_options.toml index 30e715c..816fdfc 100644 --- a/res/config/all_options.toml +++ b/res/config/all_options.toml @@ -41,6 +41,7 @@ version = "1.0" # Version of this configuration file enforce_ec = true # Enforce energy conservation in tidal response calculations, improves stability in fluid-mush cases at low forcing frequencies (< 1e-7 Hz). This does not affect the Lovenumbers. optimize_scales = false # Optimize non-dimensionalization scales for the relaxation method. This can improve convergence and stability, especially for low forcing frequencies (< 1e-7 Hz). However, do not use with Bigfloat precision. solid_shell = true # Insert an infinitesimal solid shell around the core. This patches an issue where y2 and y4 become decoupled and cause the solution to diverge in fluid layers. Only use with solid1d_relax or solid1d_mush_relax. + cap_LN = false # Clamp each mode's Re(k2)/Im(k2) to 3x/2x the fluid Love-number limit for its degree n (3/(2(n-1)), heating is rescaled to match via enforce_ec. min_frac = 0.02 # Minimal segment radius fraction before smoothing [dimensionless] visc_l = 1e2 # Pure liquid viscosity [Pa s] diff --git a/src/Obliqua.jl b/src/Obliqua.jl index 57e6cc2..3535147 100644 --- a/src/Obliqua.jl +++ b/src/Obliqua.jl @@ -315,6 +315,14 @@ module Obliqua optimize_scales = cfg["orbit"]["obliqua"]["optimize_scales"] solid_shell = cfg["orbit"]["obliqua"]["solid_shell"] + # Optional: cap each mode's tidal/load Love number at a fixed, + # physically-motivated magnitude (Re: 1.5, the homogeneous- + # incompressible-body elastic k2 limit; Im: 1.0, mirroring the + # "potentially unbound" thresholds used by PROTEUS's own + # plot_lovenumber diagnostic). Not in req_keys: absent in older + # configs/callers, defaults to off. + cap_LN = get(cfg["orbit"]["obliqua"], "cap_LN", false) + min_frac = cfg["orbit"]["obliqua"]["min_frac"] visc_l = cfg["orbit"]["obliqua"]["visc_l"] @@ -373,6 +381,7 @@ module Obliqua enforce_ec = true_if_true(enforce_ec) optimize_scales = true_if_true(optimize_scales) solid_shell = true_if_true(solid_shell) + cap_LN = true_if_true(cap_LN) # convert "none" to nothing module_solid = nothing_if_none(module_solid) @@ -827,9 +836,21 @@ module Obliqua # if segment is water elseif seg == "water" - # calculate water tides in water region + # calculate water tides in water region knms_T[iss, iseg], knms_L[iss, iseg] = 0., 0. # no expression for this yet - @warn "Water layers are currently not supported. Skipping this segment..." + @warn "Water layers are currently not supported. Skipping this segment..." + end + + # Cap this mode's tidal/load Love number before the + # enforce_ec block below uses knms_T to rescale prf_total, + # so the heating profile stays consistent with whatever + # Love number is actually reported (rather than the heating + # reflecting an uncapped, potentially resonant, value while + # the reported Love number is capped). See cap_LN's + # docstring for why this exists. + if cap_LN + knms_T[iss, iseg] = cap_lovenumber(knms_T[iss, iseg], n_i) + knms_L[iss, iseg] = cap_lovenumber(knms_L[iss, iseg], n_i) end if interp_previous @@ -1171,12 +1192,48 @@ module Obliqua function nothing_if_none(val) if val == "none" return nothing - else - return val + else + return val end end + """ + cap_lovenumber(k_val::precc, n::Int) + + Clamp a single mode's Love number to a fixed multiple of the classical + fluid (zero-rigidity) Love number limit for degree `n`, as a stop-gap + against dynamic-tide (`inertial_terms=true`) normal-mode resonances + producing values far outside the range a real, damped solid body can + sustain for the duration PROTEUS then treats as constant (there is + currently no timestepping fine enough to resolve a resonance crossing + narrower than a macro-step; see `cap_LN`'s call site). + + The baseline is `k_n^fluid = 3 / (2(n-1))`, the μ -> 0 limit of the + homogeneous, incompressible, self-gravitating elastic sphere Love + number (classical result, see e.g. Munk & MacDonald, *The Rotation of + the Earth*, 1960 -- verify the exact equation before citing elsewhere; + not independently checked against a primary source here). For n=2 this + is the textbook k2 = 3/2. + + The Re/Im multipliers (3x, 2x) are an empirical choice, not a derived + physical bound -- there is no theorem bounding a real damped body's + *dynamic*-tide Love number at a fixed multiple of the fluid limit (a + resonance peak is bounded only by whatever damping is actually + present). They are chosen to give genuine, lightly-damped resonance + peaks room to appear while still bounding the clearly-unphysical + macro-step energy injections seen in practice; adjust them if that + balance turns out wrong for a given case. + """ + function cap_lovenumber(k_val::precc, n::Int)::precc + n = max(2, n) + fluid_limit = 3.0 / (2.0 * (n - 1)) + re_max = 3.0 * fluid_limit + im_max = 2.0 * fluid_limit + return precc(clamp(real(k_val), -re_max, re_max), clamp(imag(k_val), -im_max, im_max)) + end + + """ get_layers(r, η, η_l, η_s; min_frac=0.02) diff --git a/test/test.toml b/test/test.toml index 91183e0..72137a8 100644 --- a/test/test.toml +++ b/test/test.toml @@ -40,6 +40,7 @@ version = "1.0" # Version of this configuration file enforce_ec = false # Enforce energy conservation in tidal response calculations, improves stability in fluid-mush cases at low forcing frequencies (< 1e-7 Hz). This does not affect the Lovenumbers. optimize_scales = false # Optimize non-dimensionalization scales for the relaxation method. This can improve convergence and stability, especially for low forcing frequencies (< 1e-7 Hz). However, do not use with Bigfloat precision. solid_shell = false # Insert an infinitesimal solid shell around the core. This patches an issue where y2 and y4 become decoupled and cause the solution to diverge in fluid layers. Only use with solid1d_relax or solid1d_mush_relax. + cap_LN = false # Clamp each mode's Re(k2)/Im(k2) to 3x/2x the fluid Love-number limit for its degree n (3/(2(n-1)); for n=2, Re<=4.5, Im<=3.0). A stop-gap against dynamic-tide (inertial_terms=true) normal-mode resonances producing values a real damped body could not sustain for a full macro-step; heating is rescaled to match via enforce_ec. min_frac = 0.02 # Minimal segment radius fraction before smoothing [dimensionless] visc_l = 1e2 # Pure liquid viscosity [Pa s] From 4f8b12cd840f7a3c02b0869374fce6e85a858a61 Mon Sep 17 00:00:00 2001 From: Marijn Date: Sun, 13 Sep 2026 22:17:45 +0000 Subject: [PATCH 07/11] Config option. Test. --- res/config/all_options.toml | 2 +- test/test_obliqua.jl | 48 +++++++++++++++++++++++++++++++++++++ 2 files changed, 49 insertions(+), 1 deletion(-) diff --git a/res/config/all_options.toml b/res/config/all_options.toml index 816fdfc..36e6acd 100644 --- a/res/config/all_options.toml +++ b/res/config/all_options.toml @@ -114,7 +114,7 @@ version = "1.0" # Version of this configuration file [orbit.obliqua.fluid] sigma_R = 2e-4 # Rayleigh drag coefficient at interface [dimensionless] sigma_R_inf = 1e-4 # Rayleigh drag coefficient in pure fluid [dimensionless] - sigma_R_prf = "dynamic_interp" # Rayleigh drag profile ("uniform", "exp", "linear", "quadratic", "dynamic") + sigma_R_prf = "dynamic_interp" # Rayleigh drag profile ("uniform", "exp", "linear", "quadratic", "dynamic", "dynamic_interp") H_R = 1e5 # Rayleigh drag scale height [m] efficiency = 1e-2 # Rayleigh drag efficiency at core interface [dimensionless] diff --git a/test/test_obliqua.jl b/test/test_obliqua.jl index fa9e47d..1df1a08 100644 --- a/test/test_obliqua.jl +++ b/test/test_obliqua.jl @@ -25,6 +25,54 @@ using Obliqua.solid1d_relax.common @test Obliqua.open_config(config_path) !== nothing end + # ========================================================================= + # cap_lovenumber: n-dependent Love-number clamp used by cap_LN + # ========================================================================= + @testset "cap_lovenumber" begin + # n=2: fluid limit 3/(2*(2-1)) = 1.5, so re_max=4.5, im_max=3.0 + @testset "n=2 bounds" begin + # Within bounds: passed through unchanged (both signs) + k_small = precc(1.0, -0.5) + @test Obliqua.cap_lovenumber(k_small, 2) == k_small + + # Re exceeds +4.5: clamped to exactly the bound, Im untouched + k_re_hi = precc(14.19, -0.60) + capped = Obliqua.cap_lovenumber(k_re_hi, 2) + @test real(capped) == 4.5 + @test imag(capped) == -0.60 + + # Re exceeds -4.5 (negative side): clamped to exactly -4.5 + k_re_lo = precc(-14.19, 0.0) + @test real(Obliqua.cap_lovenumber(k_re_lo, 2)) == -4.5 + + # Im exceeds +3.0 and -3.0: clamped to exactly the bound + @test imag(Obliqua.cap_lovenumber(precc(0.0, 15.5), 2)) == 3.0 + @test imag(Obliqua.cap_lovenumber(precc(0.0, -15.5), 2)) == -3.0 + + # Exactly at the boundary is not altered + @test Obliqua.cap_lovenumber(precc(4.5, 3.0), 2) == precc(4.5, 3.0) + end + + # Discrimination guard: n=3 must give a DIFFERENT (smaller) bound than + # n=2, not the same fixed number for every degree -- this is what + # distinguishes the n-dependent formula from the earlier fixed-value + # (1.5/1.0) implementation. + @testset "degree dependence (n=3)" begin + fluid_limit_n3 = 3.0 / (2.0 * (3 - 1)) # = 0.75 + @test Obliqua.cap_lovenumber(precc(10.0, 10.0), 3) == + precc(3.0 * fluid_limit_n3, 2.0 * fluid_limit_n3) + @test real(Obliqua.cap_lovenumber(precc(10.0, 0.0), 3)) < + real(Obliqua.cap_lovenumber(precc(10.0, 0.0), 2)) + end + + # n=1 (translation, undefined fluid limit) is guarded up to n=2 + # rather than dividing by zero or clamping to zero. + @testset "n=1 guarded to n=2 bound" begin + @test Obliqua.cap_lovenumber(precc(10.0, 10.0), 1) == + Obliqua.cap_lovenumber(precc(10.0, 10.0), 2) + end + end + # ========================================================================= # 1. Solid1D Module Tests # ========================================================================= From 93f729ffe06a14ed35161c14e8516dc945e4d12c Mon Sep 17 00:00:00 2001 From: Marijn Date: Mon, 14 Sep 2026 09:47:01 +0000 Subject: [PATCH 08/11] Cleaned up comments. --- src/Obliqua.jl | 31 ++----------------------------- 1 file changed, 2 insertions(+), 29 deletions(-) diff --git a/src/Obliqua.jl b/src/Obliqua.jl index 3535147..69fa145 100644 --- a/src/Obliqua.jl +++ b/src/Obliqua.jl @@ -841,13 +841,7 @@ module Obliqua @warn "Water layers are currently not supported. Skipping this segment..." end - # Cap this mode's tidal/load Love number before the - # enforce_ec block below uses knms_T to rescale prf_total, - # so the heating profile stays consistent with whatever - # Love number is actually reported (rather than the heating - # reflecting an uncapped, potentially resonant, value while - # the reported Love number is capped). See cap_LN's - # docstring for why this exists. + # Cap this mode's tidal/load Love number before the enforce_ec block if cap_LN knms_T[iss, iseg] = cap_lovenumber(knms_T[iss, iseg], n_i) knms_L[iss, iseg] = cap_lovenumber(knms_L[iss, iseg], n_i) @@ -1202,28 +1196,7 @@ module Obliqua cap_lovenumber(k_val::precc, n::Int) Clamp a single mode's Love number to a fixed multiple of the classical - fluid (zero-rigidity) Love number limit for degree `n`, as a stop-gap - against dynamic-tide (`inertial_terms=true`) normal-mode resonances - producing values far outside the range a real, damped solid body can - sustain for the duration PROTEUS then treats as constant (there is - currently no timestepping fine enough to resolve a resonance crossing - narrower than a macro-step; see `cap_LN`'s call site). - - The baseline is `k_n^fluid = 3 / (2(n-1))`, the μ -> 0 limit of the - homogeneous, incompressible, self-gravitating elastic sphere Love - number (classical result, see e.g. Munk & MacDonald, *The Rotation of - the Earth*, 1960 -- verify the exact equation before citing elsewhere; - not independently checked against a primary source here). For n=2 this - is the textbook k2 = 3/2. - - The Re/Im multipliers (3x, 2x) are an empirical choice, not a derived - physical bound -- there is no theorem bounding a real damped body's - *dynamic*-tide Love number at a fixed multiple of the fluid limit (a - resonance peak is bounded only by whatever damping is actually - present). They are chosen to give genuine, lightly-damped resonance - peaks room to appear while still bounding the clearly-unphysical - macro-step energy injections seen in practice; adjust them if that - balance turns out wrong for a given case. + fluid (zero-rigidity) Love number limit for degree `n`. """ function cap_lovenumber(k_val::precc, n::Int)::precc n = max(2, n) From 1d1840e051c4e6d9283ac23f174f647106bbcfbd Mon Sep 17 00:00:00 2001 From: Marijn Date: Mon, 14 Sep 2026 13:37:46 +0000 Subject: [PATCH 09/11] Added new config items to docs. --- docs/src/how-to-guides/config_file.md | 3 ++- docs/src/reference/configuration-file.md | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/src/how-to-guides/config_file.md b/docs/src/how-to-guides/config_file.md index b3da7cd..9c06d1b 100644 --- a/docs/src/how-to-guides/config_file.md +++ b/docs/src/how-to-guides/config_file.md @@ -72,13 +72,14 @@ Controls the tidal response model. | `optimize_scales` | bool | Boolean flag to optimize scaling factors for numerical stability. | | `solid_shell` | bool | Boolean flag to add an infinitesimal solid shell around the core to couple y2 and y4 in fluid mantles. | | `min_frac` | float | Minimum segment fraction of total mantle before it is considered. | +| `cap_LN` | bool | Boolean flag to cap the Love number response to avoid divergences. | | `visc_l` | float | Liquid viscosity. | | `visc_lus` | float | Liquid-Mush handoff viscosity. | | `visc_s` | float | Solid viscosity. | | `visc_sus` | float | Solid-Mush handoff viscosity. | | `n` | array | Radial dependence exponent in $(r/a)^n$. | | `m` | array | Tidal harmonic (e.g., $m=2$ for semidiurnal tides). | -| `spectrum` | str | Frequency sampling strategy (`"full"` or `"adaptive"`). | +| `spectrum` | str | Frequency sampling strategy (`"full"`, `"adaptive"`, or `"legacy"`). | | `N_sigma` | int | Number of sampled forcing frequencies. | | `p_min` | float | Minimum period ($\log_{10}$ kyr). | | `p_max` | float | Maximum period ($\log_{10}$ kyr). | diff --git a/docs/src/reference/configuration-file.md b/docs/src/reference/configuration-file.md index 810c530..72a1086 100644 --- a/docs/src/reference/configuration-file.md +++ b/docs/src/reference/configuration-file.md @@ -66,6 +66,7 @@ Controls the tidal response model. | `enforce_ec` | bool | Enforce energy conservation in tidal response calculations. Improves stability in fluid-mush cases at low forcing frequencies (< 1e-7 Hz); does not affect the Love numbers. | | `optimize_scales` | bool | Optimize non-dimensionalization scales for the relaxation method, for numerical stability at low forcing frequencies. Do not combine with BigFloat precision. | | `solid_shell` | bool | Insert an infinitesimal solid shell around the core to patch a $y_2$/$y_4$ decoupling instability in fluid layers. Only relevant for `solid1d-relax` or `solid1d-mush-relax`. | +| `cap_LN` | bool | Clamp each mode's Re(k2)/Im(k2) to 3x/2x the fluid Love-number limit for its degree n, rescaling heating to match via `enforce_ec`. | #### Rheology and Viscosity @@ -90,7 +91,7 @@ See [Forcing Frequency](@ref) for the underlying model. | :--- | :--- | :--- | | `n` | array | Radial dependence exponent(s) in $(r/a)^n$; since $r \ll a$, only $n=2$ contributes significantly. | | `m` | array | Tidal harmonic(s) of the true anomaly (e.g. $m=2$ semidiurnal, $m=1$ diurnal). | -| `spectrum` | str | Frequency sampling strategy: `"full"` samples the whole k2 spectrum, `"adaptive"` samples only the region of interest, `"legacy"` reproduces the original LovePy module (hardcoded low-eccentricity, spin-synchronous $(n,m,k) = (2,0,1),(2,2,1),(2,2,3)$ triplet evaluated at a single forcing frequency $\omega$). `"legacy"` overrides `n`, `m`, `s_min`, and `s_max`. | +| `spectrum` | str | Frequency sampling strategy: `"full"` samples the whole k2 spectrum, `"adaptive"` samples only the region of interest, `"legacy"` reproduces the original LovePy module (hardcoded low-eccentricity, spin-synchronous $(n,m,k) = (2,0,1),(2,2,1),(2,2,3)$ triplet evaluated at a single forcing frequency. | | `N_sigma` | int | Number of probe frequencies to evaluate k2 at (used when `spectrum = "full"`). | | `p_min` | float | Minimum period for orbital and axial frequencies [$\log_{10}$ kyr]. | | `p_max` | float | Maximum period for orbital and axial frequencies [$\log_{10}$ kyr]. | From b8c4d3001986ccbb2465fb0a8bf459801b9fe2bd Mon Sep 17 00:00:00 2001 From: Marijn Date: Mon, 14 Sep 2026 14:44:11 +0000 Subject: [PATCH 10/11] Fix copilot comments. --- docs/src/reference/configuration-file.md | 2 +- res/config/all_options.toml | 2 +- src/Obliqua.jl | 32 ++--- test/test_obliqua.jl | 151 +++++++++++++++++++++++ test/test_solid1d_mush.jl | 26 ++++ 5 files changed, 191 insertions(+), 22 deletions(-) diff --git a/docs/src/reference/configuration-file.md b/docs/src/reference/configuration-file.md index 72a1086..7be3b5e 100644 --- a/docs/src/reference/configuration-file.md +++ b/docs/src/reference/configuration-file.md @@ -91,7 +91,7 @@ See [Forcing Frequency](@ref) for the underlying model. | :--- | :--- | :--- | | `n` | array | Radial dependence exponent(s) in $(r/a)^n$; since $r \ll a$, only $n=2$ contributes significantly. | | `m` | array | Tidal harmonic(s) of the true anomaly (e.g. $m=2$ semidiurnal, $m=1$ diurnal). | -| `spectrum` | str | Frequency sampling strategy: `"full"` samples the whole k2 spectrum, `"adaptive"` samples only the region of interest, `"legacy"` reproduces the original LovePy module (hardcoded low-eccentricity, spin-synchronous $(n,m,k) = (2,0,1),(2,2,1),(2,2,3)$ triplet evaluated at a single forcing frequency. | +| `spectrum` | str | Frequency sampling strategy: `"full"` samples the whole k2 spectrum, `"adaptive"` samples only the region of interest, `"legacy"` reproduces the original LovePy module (hardcoded low-eccentricity, spin-synchronous $(n,m,k) = (2,0,1),(2,2,1),(2,2,3)$ triplet evaluated at a single forcing frequency $\omega$). `"legacy"` overrides `n`, `m`, `s_min`, and `s_max`. | | `N_sigma` | int | Number of probe frequencies to evaluate k2 at (used when `spectrum = "full"`). | | `p_min` | float | Minimum period for orbital and axial frequencies [$\log_{10}$ kyr]. | | `p_max` | float | Maximum period for orbital and axial frequencies [$\log_{10}$ kyr]. | diff --git a/res/config/all_options.toml b/res/config/all_options.toml index 36e6acd..28405c0 100644 --- a/res/config/all_options.toml +++ b/res/config/all_options.toml @@ -41,7 +41,7 @@ version = "1.0" # Version of this configuration file enforce_ec = true # Enforce energy conservation in tidal response calculations, improves stability in fluid-mush cases at low forcing frequencies (< 1e-7 Hz). This does not affect the Lovenumbers. optimize_scales = false # Optimize non-dimensionalization scales for the relaxation method. This can improve convergence and stability, especially for low forcing frequencies (< 1e-7 Hz). However, do not use with Bigfloat precision. solid_shell = true # Insert an infinitesimal solid shell around the core. This patches an issue where y2 and y4 become decoupled and cause the solution to diverge in fluid layers. Only use with solid1d_relax or solid1d_mush_relax. - cap_LN = false # Clamp each mode's Re(k2)/Im(k2) to 3x/2x the fluid Love-number limit for its degree n (3/(2(n-1)), heating is rescaled to match via enforce_ec. + cap_LN = false # Clamp each mode's Re(k2)/Im(k2) to 3x/2x the fluid Love-number limit for its degree n (3/(2(n-1)); for n=2, Re<=4.5, Im<=3.0). Heating is rescaled to match via enforce_ec. min_frac = 0.02 # Minimal segment radius fraction before smoothing [dimensionless] visc_l = 1e2 # Pure liquid viscosity [Pa s] diff --git a/src/Obliqua.jl b/src/Obliqua.jl index 69fa145..4695473 100644 --- a/src/Obliqua.jl +++ b/src/Obliqua.jl @@ -314,13 +314,6 @@ module Obliqua enforce_ec = cfg["orbit"]["obliqua"]["enforce_ec"] optimize_scales = cfg["orbit"]["obliqua"]["optimize_scales"] solid_shell = cfg["orbit"]["obliqua"]["solid_shell"] - - # Optional: cap each mode's tidal/load Love number at a fixed, - # physically-motivated magnitude (Re: 1.5, the homogeneous- - # incompressible-body elastic k2 limit; Im: 1.0, mirroring the - # "potentially unbound" thresholds used by PROTEUS's own - # plot_lovenumber diagnostic). Not in req_keys: absent in older - # configs/callers, defaults to off. cap_LN = get(cfg["orbit"]["obliqua"], "cap_LN", false) min_frac = cfg["orbit"]["obliqua"]["min_frac"] @@ -518,18 +511,11 @@ module Obliqua elseif spectrum == "legacy" # Reproduce the original LovePy module: hardcode the three dominant - # low-eccentricity (n, m, k) modes and force an identical forcing - # frequency across all three, so a single k2 Love number spectrum - # is evaluated (at ω) instead of the full/adaptive mode expansion. - # This assumes spin-orbit synchronisation and e << 1, matching the - # simplifications baked into LovePy; it is not a general-purpose - # replacement for "adaptive". + # low-eccentricity (n, m, k) modes mathcing LovePy. nmk = [(2, 0, 1), (2, 2, 1), (2, 2, 3)] # LovePy hardcodes forcing frequency = orbital mean motion (omega), - # which is only physically exact for spin-orbit synchronous rotation - # (axial == omega): warn loudly if that assumption is violated, - # since this mode silently ignores `axial` otherwise. + # warn loudly if that assumption is violated. if !isapprox(axial, omega; rtol=1e-3) @warn "Legacy spectrum assumes spin-orbit synchronisation (axial == omega), but axial=$axial rad/s and omega=$omega rad/s differ by more than 0.1%. Forcing frequency is still hardcoded to omega; results will not match a self-consistent tidal calculation for this rotation state." end @@ -545,8 +531,6 @@ module Obliqua N_σ = length(σ_range) @info "Using legacy (LovePy-compatible) spectrum: (n, m, k) = (2, 0, 1), (2, 2, 1), (2, 2, 3), all evaluated at ω = $omega." - else - throw("Invalid spectrum value: $spectrum. Must be 'adaptive', 'full', or 'legacy'.") end # get frequency dependent complex shear modulus per mode @@ -836,9 +820,9 @@ module Obliqua # if segment is water elseif seg == "water" - # calculate water tides in water region + # calculate water tides in water region knms_T[iss, iseg], knms_L[iss, iseg] = 0., 0. # no expression for this yet - @warn "Water layers are currently not supported. Skipping this segment..." + @warn "Water layers are currently not supported. Skipping this segment..." end # Cap this mode's tidal/load Love number before the enforce_ec block @@ -858,6 +842,14 @@ module Obliqua prf_total[iss, i_sp:i_ep] .+= Δprf knms_T[iss, iseg-1] += ΔkT knms_L[iss, iseg-1] += ΔkL + + # Re-apply the cap: the previous segment's Love number was + # already clamped above, but this interpolation increment + # is added afterwards and can push it back out of bounds. + if cap_LN + knms_T[iss, iseg-1] = cap_lovenumber(knms_T[iss, iseg-1], n_i) + knms_L[iss, iseg-1] = cap_lovenumber(knms_L[iss, iseg-1], n_i) + end end # repeat for all probe forcing frequencies diff --git a/test/test_obliqua.jl b/test/test_obliqua.jl index 1df1a08..df60131 100644 --- a/test/test_obliqua.jl +++ b/test/test_obliqua.jl @@ -73,6 +73,157 @@ using Obliqua.solid1d_relax.common end end + # ========================================================================= + # spectrum = "legacy": hardcoded LovePy-compatible (n,m,k) triplet + # ========================================================================= + @testset "legacy spectrum" begin + # runtests_mantle.json has axial != omega, so this run also exercises + # the spin-synchronisation mismatch @warn path (Obliqua.jl ~L534). + cfg["orbit"]["obliqua"]["module_solid"] = "solid1d" + cfg["orbit"]["obliqua"]["module_fluid"] = "none" + cfg["orbit"]["obliqua"]["module_mushy"] = "none" + cfg["orbit"]["obliqua"]["material_mu"] = "andrade" + cfg["orbit"]["obliqua"]["spectrum"] = "legacy" + + omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, phi, ncalc = + load.load_interior_mush_full(interior_json_path, false) + + perm = Obliqua.interior.get_permeability(phi, cfg) + perm, phi = Obliqua.interior.limit_porosity(perm, phi, cfg) + bulkd = Obliqua.interior.get_drained_bulk(bulk, phi, cfg) + + power_blk_expt = 1.728319570165342e6 + imag_k2_expt = 0.0008308957062851288 + + power_prf, power_blk, nmk, sigma_range, LNk = Obliqua.run_tides( + omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, bulkd, phi, perm, cfg + ) + imag_k2 = -imag.(LNk) + + # Mode selection matches the hardcoded LovePy (n,m,k) triplet exactly. + @test nmk == [(2, 0, 1), (2, 2, 1), (2, 2, 3)] + + # All three modes share one forcing frequency (dropping the sign, as + # LovePy does), this is the whole point of the "legacy" mode. + @test length(unique(sigma_range)) == 1 + @test sigma_range[1] ≈ omega + + # ...and therefore also collapse onto a single shared Im(k2), not + # three independently-evaluated Love numbers. + @test length(unique(round.(imag_k2, sigdigits=12))) == 1 + @test all(isapprox.(imag_k2, imag_k2_expt; rtol=rtol)) + + @test isapprox(power_blk, power_blk_expt; rtol=rtol) + @test power_blk > 0.0 + + # Restore spectrum for subsequent testsets. + cfg["orbit"]["obliqua"]["spectrum"] = "adaptive" + end + + @testset "invalid spectrum value" begin + cfg["orbit"]["obliqua"]["module_solid"] = "solid1d" + cfg["orbit"]["obliqua"]["module_fluid"] = "none" + cfg["orbit"]["obliqua"]["module_mushy"] = "none" + cfg["orbit"]["obliqua"]["material_mu"] = "andrade" + cfg["orbit"]["obliqua"]["spectrum"] = "not_a_real_spectrum_mode" + + omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, phi, ncalc = + load.load_interior_mush_full(interior_json_path, false) + + perm = Obliqua.interior.get_permeability(phi, cfg) + perm, phi = Obliqua.interior.limit_porosity(perm, phi, cfg) + bulkd = Obliqua.interior.get_drained_bulk(bulk, phi, cfg) + + @test_throws String Obliqua.run_tides( + omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, bulkd, phi, perm, cfg + ) + + # Restore spectrum for subsequent testsets. + cfg["orbit"]["obliqua"]["spectrum"] = "adaptive" + end + + # ========================================================================= + # cap_LN pipeline wiring: cap_lovenumber is applied inside run_tides + # ========================================================================= + @testset "cap_LN pipeline wiring" begin + cfg["orbit"]["obliqua"]["module_solid"] = "solid1d" + cfg["orbit"]["obliqua"]["module_fluid"] = "none" + cfg["orbit"]["obliqua"]["module_mushy"] = "none" + cfg["orbit"]["obliqua"]["material_mu"] = "andrade" + cfg["orbit"]["obliqua"]["spectrum"] = "adaptive" + cfg["orbit"]["obliqua"]["cap_LN"] = true + + omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, phi, ncalc = + load.load_interior_mush_full(interior_json_path, false) + + perm = Obliqua.interior.get_permeability(phi, cfg) + perm, phi = Obliqua.interior.limit_porosity(perm, phi, cfg) + bulkd = Obliqua.interior.get_drained_bulk(bulk, phi, cfg) + + # For this config, Im(k2)/Re(k2) sit far below the fluid-limit + # cap, so clamping must be a no-op: cap_LN=true must reproduce the + # exact same result as the "solid1d module / Andrade Rheology" + # cap_LN=false test above. + power_blk_expt = 1.093766208671846e6 + imag_k2_expt = 0.0014986632230270696 + + _, power_blk, _, _, LNk = Obliqua.run_tides( + omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, bulkd, phi, perm, cfg + ) + imag_k2 = -imag.(LNk) + + @test isapprox(power_blk, power_blk_expt; rtol=rtol) + @test all(isapprox.(imag_k2, imag_k2_expt; rtol=rtol)) + + # Physical bound implied by cap_LN for n=2 (im_max = 2 * 3/(2*(2-1))). + @test all(imag_k2 .<= 3.0) + + # Restore configuration for subsequent testsets. + cfg["orbit"]["obliqua"]["cap_LN"] = false + end + + # The `module_mushy == "interp"` path adds a correction (ΔkT/ΔkL) onto the + # *previous* segment's Love number strictly after that segment's own + # cap_LN clamp already ran, which could silently push it back out of + # bounds. Check that the re-clamp after interpolation is applied correctly. + @testset "cap_LN with module_mushy=interp (re-clamp after interpolation)" begin + cfg["orbit"]["obliqua"]["module_solid"] = "solid1d" + cfg["orbit"]["obliqua"]["module_fluid"] = "fluid1d" + cfg["orbit"]["obliqua"]["module_mushy"] = "interp" + cfg["orbit"]["obliqua"]["s_min"] = -2 + cfg["orbit"]["obliqua"]["s_max"] = 6 + cfg["orbit"]["obliqua"]["material_mu"] = "andrade" + cfg["orbit"]["obliqua"]["cap_LN"] = true + + omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, phi, ncalc = + load.load_interior_mush_full(interior_json_path, false) + + perm = Obliqua.interior.get_permeability(phi, cfg) + perm, phi = Obliqua.interior.limit_porosity(perm, phi, cfg) + bulkd = Obliqua.interior.get_drained_bulk(bulk, phi, cfg) + + # Same config as "Complete model" (cap_LN=false there). + power_blk_expt = 2.8317958612079725e9 + imag_k2_expt = [0.011581176958188387, 0.01108652173121107, 0.010590988673882713, 0.010094581843049406, 0.009597312117327252, 0.009099198763898173, 0.008600271466978374, 0.008100573001742285, 0.007600162830519118] + + _, power_blk, _, _, LNk = Obliqua.run_tides( + omega, axial, ecc, sma, S_mass, rho, radius, visc, shear, bulk, bulkd, phi, perm, cfg + ) + imag_k2 = -imag.(LNk) + + @test isapprox(power_blk, power_blk_expt; rtol=rtol) + @test all(isapprox.(imag_k2, imag_k2_expt; rtol=rtol, atol=atol)) + @test all(imag_k2 .<= 3.0) + + # Restore configuration (in particular s_min/s_max) to the test.toml + # defaults for subsequent testsets. + cfg["orbit"]["obliqua"]["cap_LN"] = false + cfg["orbit"]["obliqua"]["module_fluid"] = "none" + cfg["orbit"]["obliqua"]["module_mushy"] = "none" + cfg["orbit"]["obliqua"]["s_min"] = 1 + cfg["orbit"]["obliqua"]["s_max"] = 1 + end + # ========================================================================= # 1. Solid1D Module Tests # ========================================================================= diff --git a/test/test_solid1d_mush.jl b/test/test_solid1d_mush.jl index 1c6545e..10a3197 100644 --- a/test/test_solid1d_mush.jl +++ b/test/test_solid1d_mush.jl @@ -148,6 +148,32 @@ using Obliqua.constants @test size(y_sol) == (8, size(r, 1)-1, size(r, 2)) end + # ------------------------------------------------------------------------- + # Test 8b: compute_M with a porous bottom layer (porous_layer[1] == true) + # ------------------------------------------------------------------------- + # Regression test for the core-basis-radius fix: get_Ic must be evaluated + # at the core-mantle boundary itself (rs[1,1]), not the top of the bottom + # sublayer (rs[end,1]). Test 8 above only exercises the `else` branch of + # `if porous_layer[1]`; a porous bottom layer takes the `if` branch, which + # calls get_Ic with the full 8-component core basis instead. + @testset "Interior Boundary Matrix with porous core-adjacent layer" begin + r = [0.1 10.0; 5.0 20.0; 10.0 30.0] + g = [0.01 1.0; 0.5 1.5; 1.0 2.0] + ρ = [3000.0, 3200.0]; μ = [1e10+0im, 1.2e10+0im]; K = [2e10+0im, 2.5e10+0im] + ρₗ = [1000.0, 1000.0]; Kl = [2e9, 2e9]; Kd = [1e10+0im, 1.1e10+0im] + α = [0.5+0im, 0.6+0im]; ηₗ = [1.0, 1.0] + ϕ = [0.1, 0.0] # porous_layer = [true, false]: bottom (core-adjacent) layer is porous + k = [1e-12, 1e-12] + ω = 1e-3; n = 2; ρ_core = 5000.0; μ_core = complex(0.0); κ_core = complex(1e11) + + M_mat, y1_4 = solid1d_mush.compute_M(ω, r, ρ, g, μ, K, ρₗ, Kl, Kd, α, ηₗ, ϕ, k, n, ρ_core, μ_core, κ_core, ones(prec, 3); core="liquid") + + @test size(M_mat) == (4, 4) + @test size(y1_4) == (8, 4, size(r, 1)-1, size(r, 2)) + # The core basis must actually be finite/nonzero once evaluated at the CMB. + @test !all(iszero, y1_4[:, :, 1, 1]) + end + # ------------------------------------------------------------------------- # Test 9: Physical Property Tensors (Strain, Displacement, Pore Pressure) # ------------------------------------------------------------------------- From e5a6a229a094879f770666c728ec8643573b5bdd Mon Sep 17 00:00:00 2001 From: Marijn Date: Tue, 15 Sep 2026 15:13:30 +0000 Subject: [PATCH 11/11] Fixed typos. --- docs/src/explanation/iter_proc.md | 12 ++++++------ docs/src/explanation/main_loop.md | 6 +++--- docs/src/explanation/post_proc.md | 2 +- docs/src/how-to-guides/config_file.md | 2 +- docs/src/how-to-guides/usage.md | 4 +++- docs/src/reference/solid-phase.md | 5 ++--- docs/src/reference/surface-loading.md | 8 ++++---- 7 files changed, 20 insertions(+), 19 deletions(-) diff --git a/docs/src/explanation/iter_proc.md b/docs/src/explanation/iter_proc.md index ffb67d8..207f268 100644 --- a/docs/src/explanation/iter_proc.md +++ b/docs/src/explanation/iter_proc.md @@ -19,11 +19,11 @@ deformation, a.k.a. the $n$th degree Lovenumber $k_n(\sigma)$; the planet-wide loading Lovenumber $k'_n(\sigma)$; and the normalized heating profile in the segment. All the details regarding these models will be given in the corresponding sections below. The currently -available tidal models are "solid0d", "solid1d", "solid1d_relax", -"solid1d_mush", "solid1d_mush_relax"; "fluid0d", "fluid1d"; "interp", -"none". For details see the Reference documantation. +available tidal models are `"solid0d"`, `"solid1d"`, `"solid1d_relax"`, +`"solid1d_mush"`, `"solid1d_mush_relax"`, `"solid1d_equil_relax"`; `"fluid0d"`, `"fluid1d"`; `"interp"`, +`"none"`. For details see the Reference documantation. -The "interp" model requires knowledge of heating at both interfaces, as +The `"interp"` model requires knowledge of heating at both interfaces, as such an additional code block is included to update the heating in the -"interp" region during the tidal calculation in the next segment after -the "interp" region. +`"interp"` region during the tidal calculation in the next segment after +the `"interp"` region. diff --git a/docs/src/explanation/main_loop.md b/docs/src/explanation/main_loop.md index 2d122a9..4a7c271 100644 --- a/docs/src/explanation/main_loop.md +++ b/docs/src/explanation/main_loop.md @@ -40,15 +40,15 @@ determined by the code. The desired range increases with eccentricity, for our purposes the desired range $K$ contains ```math -\{k \in K \, \forall \, k : X^{-(n+1), m}_k \geq 0.01\ | k \in Z\} +K = \{k \in \mathbb{Z} \, \big| \, |X^{-(n+1), m}_k| \geq 0.001 \} ``` where $X^{-(n+1), m}_k$ is the Hansen coefficient. Basically, we only include $k$ in the range $K$ if the corresponding Hansen coefficient -signifies a contribution greater than 1% to the complete tidal response. +signifies a contribution greater than 0.1% to the complete tidal response. One may specify a different criterion and generate their $K$ using the included Notebook on the Obliqua Github repository -"/examples/hansen_k_table.ipynb". +`"/examples/hansen_k_table.ipynb"`. For testing convenience Obliqua includes two modes: "full" and "adaptive". In principle, one should use "full" to test the tidal diff --git a/docs/src/explanation/post_proc.md b/docs/src/explanation/post_proc.md index 52a4ce9..775e7da 100644 --- a/docs/src/explanation/post_proc.md +++ b/docs/src/explanation/post_proc.md @@ -48,7 +48,7 @@ The tidal potential is $$U_{n,m,1} = \frac{GM}{a} \left(\frac{R}{a}\right)^n A_{n,m,1}(e)$$ -The prefactor is $$\text{prefactor} \, = \frac{(2n + 1)R}{8πG} \sigma$$ The +The prefactor is $$\text{prefactor} \, = \frac{(2n + 1)R}{8πG} \sigma$$. The normalized heating profile is then simply $$H(r, \sigma) = \tilde{H}(r, \sigma) \times |U_{n,m,1}|^2$$ diff --git a/docs/src/how-to-guides/config_file.md b/docs/src/how-to-guides/config_file.md index 9c06d1b..99dd9e8 100644 --- a/docs/src/how-to-guides/config_file.md +++ b/docs/src/how-to-guides/config_file.md @@ -71,8 +71,8 @@ Controls the tidal response model. | `enforce_ec` | bool | Boolean flag to enforce energy conservation in tidal response calculations. | | `optimize_scales` | bool | Boolean flag to optimize scaling factors for numerical stability. | | `solid_shell` | bool | Boolean flag to add an infinitesimal solid shell around the core to couple y2 and y4 in fluid mantles. | -| `min_frac` | float | Minimum segment fraction of total mantle before it is considered. | | `cap_LN` | bool | Boolean flag to cap the Love number response to avoid divergences. | +| `min_frac` | float | Minimum segment fraction of total mantle before it is considered. | | `visc_l` | float | Liquid viscosity. | | `visc_lus` | float | Liquid-Mush handoff viscosity. | | `visc_s` | float | Solid viscosity. | diff --git a/docs/src/how-to-guides/usage.md b/docs/src/how-to-guides/usage.md index 822bf60..97771ad 100644 --- a/docs/src/how-to-guides/usage.md +++ b/docs/src/how-to-guides/usage.md @@ -1,7 +1,7 @@ # Usage -This section describes how to use Obliqua. The module can be run three ways: +This section describes how to use Obliqua. The module can be run four ways: - **Full spectrum (Standalone)**: Compute the tidal ``k``-Love number response for a full spectrum of forcing frequencies. This mode is agnostic to the orbital parameters that feed into the Hansen mode weights and dissipative response. The generated spectrum can be used as a lookup table, or to probe the quantative dissipative response of the tidal model to a wide range of forcing frequencies. @@ -9,6 +9,8 @@ This section describes how to use Obliqua. The module can be run three ways: - **Adaptive (PROTEUS)**: By extension of the previous mode, the adaptive mode can be used in conjunction with the PROTEUS framework. This allows Obliqua to interact with both dynamically evolving orbital parameters and interior properties. The tidal ``k``-Love number response is computed on-the-fly and is used to update the orbital evolution while the dissipative response feedsback into the interior. +- **Legacy (Standalone)**: Compute the tidal ``k``-Love number response for the subset of dominant forcing frequencies at small eccentricities. Given that the other modes explicitely reduce to this mode, this mode is provided for legacy purposes and is not recommended for use in new applications. Specifically, the `legacy` mode forces the returned Love numbers to be the same across the included set of modes, this makes the output directly compatible with simplified orbital dynamics models that do not resolve the full tidal spectrum. + Naturally, one can also use the full spectrum mode in conjunction with PROTEUS in post-processing. This can forexample be used to validate the adaptive mode or to study the impact of different forcing frequencies on the tidal response. Moreover, it may be used to test model convergence. Below, we provide here an example run of the full spectrum mode in conjunction with PROTEUS computed in post-processing. In all cases you configure the model through a configuration file, described in the [configuration guide](@ref "Configuration file"). If you run into problems, see the [troubleshooting](@ref "Troubleshooting") page. diff --git a/docs/src/reference/solid-phase.md b/docs/src/reference/solid-phase.md index ef7f17a..8592b5c 100644 --- a/docs/src/reference/solid-phase.md +++ b/docs/src/reference/solid-phase.md @@ -94,10 +94,9 @@ surface mass load can also be written as an external potential $U'$ such that $\zeta_n = [(2n + 1)/4 \pi G a] U'_n$. $$\begin{aligned} -y_{3}(R) &= - \frac{(2n + 1)g_e}{4 \pi G R} U'_n - P_n \\ +y_{3}(R) &= - \frac{(2n + 1)g_e}{4 \pi G R} \left[\frac{G}{R} U'_n \right] - P_n \\ y_{4}(R) &= \tau_n \\ -y_{6}(R) - &= \frac{2n+1}{R} (U_n + U'_n) +y_{6}(R) &= \frac{2n+1}{R} \left(U_n + \left[\frac{G}{R} U'_n \right] \right) \end{aligned}$$ Note that `get_surface_bc!` in `src/common.jl` does not literally apply diff --git a/docs/src/reference/surface-loading.md b/docs/src/reference/surface-loading.md index 752448d..1ea5053 100644 --- a/docs/src/reference/surface-loading.md +++ b/docs/src/reference/surface-loading.md @@ -21,10 +21,10 @@ No $y_5$ term appears in the third row: in Obliqua's convention $y_6$ is already By expressing a surface mass load $\zeta_n$ as an equivalent external potential $U'$, where $\zeta_n = \frac{2n + 1}{4 \pi G R} U'_n$, the system simplifies to: $$\begin{aligned} -y_{3}(R) &= - \frac{(2n + 1)g_e}{4 \pi G R} U'_n - P_n \\ -y_{4}(R) &= \tau_n \\ -y_6(R) &= \frac{2n+1}{R} (U_n + U'_n) -\end{aligned}$$ +y_{3}(R) &= - \frac{(2n + 1)g_e}{4 \pi G R} \left[\frac{G}{R} U'_n \right] - P_n \\ +y_{4}(R) &= \tau_n \\ +y_{6}(R) &= \frac{2n+1}{R} \left(U_n + \left[\frac{G}{R} U'_n \right] \right) +\end{aligned}$$ Note that Obliqua's actual `get_surface_bc!` (`src/common.jl`) does not literally apply this $\zeta_n \leftrightarrow U'_n$ conversion; it sets $(U,U',\tau,P)$ directly as dimensionless $0$/$1$ selector flags, which for the load case numerically works out to $y_3(R) = -(2n+1)g(R)/(4\pi R^2)$ and $y_6(R) = (2n+1)G/R^2$ — see [Solid-Phase - solid1d](@ref) for the concrete tidal/load values the code actually produces.