From 2c86dff6f09a24f7aa7bb9f40be7f5e63567c857 Mon Sep 17 00:00:00 2001 From: Line Colin Date: Fri, 16 May 2025 16:48:20 +0200 Subject: [PATCH 1/3] Add a gravitational acceleration that depends on the radius (linearly, g(r) = -g_0 r e_r). Add an option in the definition of the problem and modify the matrix --- stablinrb/spherical.py | 1 + 1 file changed, 1 insertion(+) diff --git a/stablinrb/spherical.py b/stablinrb/spherical.py index 12fc67e..967ad86 100644 --- a/stablinrb/spherical.py +++ b/stablinrb/spherical.py @@ -29,6 +29,7 @@ class SphStability: chebyshev_degree: int gamma: float + uniform_gravity: bool temperature: phy.AdvDiffEq | None = phy.AdvDiffEq( bc_top=phy.Zero(), bc_bot=phy.Zero(), From d4bf9d97e706e25066a51e8a3c23e44e17efc156 Mon Sep 17 00:00:00 2001 From: Line Colin Date: Fri, 16 May 2025 17:08:47 +0200 Subject: [PATCH 2/3] Squashed commit of the following: commit 5cbbdcedca625852772a7e31fa0092b8466fb0ab Author: Line Colin Date: Fri May 16 17:05:14 2025 +0200 Add the option for uniform or radius dependant gravity acceleration This update introduces the option of choosing between uniform and radius-dependent gravity acceleration in the SphStability class. The gravity acceleration that depends of the radius is defined as: g(r) = - g_0 * r * e_r The 'uniform_gravity' boolean flag has been added to the dataclass to control this option. When set to 'False', gravity acceleration varies linearly with radius. This change required modifications in 'lmat'. --- .gitignore | 2 ++ stablinrb/spherical.py | 27 ++++++++++++++++----------- 2 files changed, 18 insertions(+), 11 deletions(-) diff --git a/.gitignore b/.gitignore index af5bdfd..6b3990f 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,5 @@ *.pdf *.swp __pycache__/ +.DS_Store + diff --git a/stablinrb/spherical.py b/stablinrb/spherical.py index 967ad86..20c8d5c 100644 --- a/stablinrb/spherical.py +++ b/stablinrb/spherical.py @@ -14,7 +14,7 @@ from .rheology import Isoviscous if typing.TYPE_CHECKING: - from typing import Callable + from typing import Callable, Optional from numpy.typing import NDArray @@ -29,17 +29,18 @@ class SphStability: chebyshev_degree: int gamma: float - uniform_gravity: bool - temperature: phy.AdvDiffEq | None = phy.AdvDiffEq( + # NOTE: add uniform gravity to choose uniform gravity or radius dependant gravity + uniform_gravity: bool = True + temperature: Optional[phy.AdvDiffEq] = phy.AdvDiffEq( bc_top=phy.Zero(), bc_bot=phy.Zero(), ref_prof=DiffusiveProf(bcs_top=Dirichlet(0.0), bcs_bot=Dirichlet(1.0)), ) - composition: phy.AdvDiffEq | None = None + composition: Optional[phy.AdvDiffEq] = None bc_mom_top: phy.BCMomentum = phy.FreeSlip() bc_mom_bot: phy.BCMomentum = phy.FreeSlip() rheology: Rheology = Isoviscous() - cooling_smo: tuple[Callable, Callable] | None = None + cooling_smo: Optional[tuple[Callable, Callable]] = None frozen_time: bool = False def name(self) -> str: @@ -87,7 +88,7 @@ def slices(self) -> Slices: return Slices(var_specs=var_specs, nnodes=self.nodes.size) def eigen_problem( - self, l_harm: int, ra_num: float | None, ra_comp: float | None = None + self, l_harm: int, ra_num: Optional[float], ra_comp: Optional[float] = None ) -> EigenvalueProblem: ops = self.operators(l_harm) dr1, dr2 = ops.diff_r(1), ops.diff_r(2) @@ -141,7 +142,11 @@ def eigen_problem( lmat.add_term(Bulk("q"), ops.lapl, "q") if temp_terms: assert ra_num is not None - lmat.add_term(Bulk("q"), -ra_num * orl1, "T") + if self.uniform_gravity==True: + lmat.add_term(Bulk("q"), -ra_num * orl1, "T") + else: + # NOTE: Modifiaction for a radius dependant gravity g(r) = -g_O r e_r and multiplcation by 1/R+ to keep g(r = R+) = - g_O + lmat.add_term(Bulk("q"), -ra_num * one * (1 - self.gamma), "T") if self.composition is not None: assert ra_comp is not None lmat.add_term(Bulk("q"), -ra_comp * orl1, "c") @@ -183,7 +188,7 @@ def eigen_problem( return EigenvalueProblem(lmat, rmat) def growth_rate( - self, harm: int, ra_num: float | None, ra_comp: float | None = None + self, harm: int, ra_num: Optional[float], ra_comp: Optional[float] = None ) -> np.floating: return np.real(self.eigen_problem(harm, ra_num, ra_comp).max_eigval()) @@ -205,7 +210,7 @@ def neutral_ra( self, harm: int, ra_guess: float = 600, - ra_comp: float | None = None, + ra_comp: Optional[float] = None, eps: float = 1.0e-8, ) -> float: """Find Ra which gives neutral stability of a given harmonic @@ -242,7 +247,7 @@ def neutral_ra( ) def fastest_mode( - self, ra_num: float, ra_comp: float | None = None, harm: int = 2 + self, ra_num: float, ra_comp: Optional[float] = None, harm: int = 2 ) -> tuple[float, int]: """Find the fastest growing mode at a given Ra""" @@ -295,7 +300,7 @@ def ran_l_mins(self) -> tuple[tuple[int, float], tuple[int, float]]: return ((1, ran_mod1), (l_mod2, ran_mod2)) def critical_ra( - self, harm: int = 2, ra_guess: float = 600, ra_comp: float | None = None + self, harm: int = 2, ra_guess: float = 600, ra_comp: Optional[float] = None ) -> tuple[float, int]: """Find the harmonic with the lowest neutral Ra From c7c90f3aa6136cc4091c1db532bcab05c767a444 Mon Sep 17 00:00:00 2001 From: Line Colin Date: Fri, 16 May 2025 17:44:03 +0200 Subject: [PATCH 3/3] uses 'if self.uniform_gravity:' for truth checks --- stablinrb/spherical.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stablinrb/spherical.py b/stablinrb/spherical.py index 20c8d5c..dd6c6c1 100644 --- a/stablinrb/spherical.py +++ b/stablinrb/spherical.py @@ -142,7 +142,7 @@ def eigen_problem( lmat.add_term(Bulk("q"), ops.lapl, "q") if temp_terms: assert ra_num is not None - if self.uniform_gravity==True: + if self.uniform_gravity: lmat.add_term(Bulk("q"), -ra_num * orl1, "T") else: # NOTE: Modifiaction for a radius dependant gravity g(r) = -g_O r e_r and multiplcation by 1/R+ to keep g(r = R+) = - g_O