Gravity

.gitignore

@@ -4,3 +4,5 @@
 *.pdf
 *.swp
 __pycache__/
+.DS_Store
+

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,16 +29,18 @@ class SphStability:
 
     chebyshev_degree: int
     gamma: float
-    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:
@@ -86,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)
@@ -140,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:
+                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")
@@ -182,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())
 
@@ -204,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
@@ -241,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"""
 
@@ -294,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