starry-phoenix · starry-phoenix · Mar 10, 2026 · Mar 10, 2026 · Mar 10, 2026 · Mar 10, 2026
diff --git a/docs/source/quick_start/example/jupyter_run.ipynb b/docs/source/quick_start/example/jupyter_run.ipynb
diff --git a/docs/source/quick_start/example/jupyter_run_algae.ipynb b/docs/source/quick_start/example/jupyter_run_algae.ipynb
diff --git a/docs/source/quick_start/example/main_original.ipynb b/docs/source/quick_start/example/main_original.ipynb
diff --git a/docs/source/user_guide/main_original.ipynb b/docs/source/user_guide/main_original.ipynb
@@ -119,7 +119,7 @@
    "source": [
     "### Model settings\n",
     "\n",
-    "for more info on all the model parameters have a look at the documentation: [userinput](http://[::1]:8000/api/spyice.parameters.user_input.html)"
+    "for more info on all the model parameters have a look at the documentation: [userinput](https://starry-phoenix.github.io/spyice-thesis/build/html/api/spyice.parameters.user_input.html)"
    ]
   },
   {
@@ -278,7 +278,7 @@
    "metadata": {},
    "source": [
     "### Visualization of Model:\n",
-    "for more info on other visualization options look at:  [visualize model](http://[::1]:8000/api/spyice.postprocess.visualise_model.html)"
+    "for more info on other visualization options look at:  [visualize model](https://starry-phoenix.github.io/spyice-thesis/build/html/api/spyice.postprocess.visualise_model.html)"
    ]
   },
   {

diff --git a/src/spyice/models/advection_diffusion.py b/src/spyice/models/advection_diffusion.py
@@ -37,6 +37,7 @@ def __init__(
         Buffo=False,
         Voller=False,
         bc_neumann=None,
+        top_temperature = None,
     ):
         """
         Args:
@@ -82,7 +83,10 @@ def __init__(
         self.t_passed = t_passed
         self.S_IC = S_IC
         self.S_bc_top = "Dirichlet"
-        self.top_temp = top_temp  # OR "T_const_250" OR "T_const_260" OR "T_const_265" OR "T_W3" OR "Stefan"
+        if top_temperature:
+            self.top_temp = top_temperature
+        else:
+            AssertionError(f"{self.argument} boundary value missing!")
         ### Compute a, b, c, d
         self.Delta_W = np.zeros(nz)
         [self.a, self.b, self.c, self.d] = update_coefficients(
@@ -218,11 +222,9 @@ def set_up_tridiagonal(self):
                     if i > 0:
                         self.lower_A[i - 1] = -self.factor1_plus[i]
                 else:
-                    self.main_A[i] = (
-                        1
-                        + self.factor1_plus[i]
-                        + self.factor1_minus[i]
-                    )
+                    self.main_A[i] = 1
+                        # + self.factor1_plus[i]
+                        # + self.factor1_minus[i]
             else:
                 self.main_A[i] = 2 * self.factor1[i] + 1.0 - self.factor2[i]
                 if i < self.nz - 1:
@@ -327,7 +329,6 @@ def unknowns_matrix(self, temperature_melt, non_constant_physical_properties=Fal
             self.temp_grad,
         )
 
-        # TODO: is X_wind necessary? since advective term is considered in the assemble_tridiagonal() matrix
         X_wind = correct_for_brine_movement(
             self.argument,
             self.X_initial,

diff --git a/src/spyice/models/sea_ice_model.py b/src/spyice/models/sea_ice_model.py
@@ -28,11 +28,6 @@
 from spyice.utils.spyice_exceptions import ConvergenceError, InvalidPhaseError
 
 plt.style.use("spyice.utils.custom")
-# plt.rcParams.update(
-#     {
-#         "text.usetex": True,
-#     }
-# )
 plt.rcParams["text.latex.preamble"].join(
     [
         r"\usepackage{dashbox}",
@@ -173,7 +168,7 @@ def t_running(self, fig, ax1, t_stefan, t_k, t_k_buffo=None, count=0):
         depth_arr = np.linspace(0, -self.preprocess_data.Z, self.preprocess_data.nz)
         ax1.plot(t_k, depth_arr, "r--")
         ax1.plot(t_stefan, depth_arr, "k")
-        if t_k_buffo is not None:
+        if self.preprocess_data.is_buffo:
             ax1.plot(t_k_buffo, depth_arr, "b-.", alpha=0.5)
         ax1.set_ylabel("Depth [m]")
         ax1.set_xlabel("Temperature [K]")
@@ -326,7 +321,6 @@ def convergence_loop_iteration(
             source_term_temperature,
             source_term_salinity,
         )
-        # TODO: truncate field variables until boundary (thickness_index)
         (
             t_km1,
             s_km1,
@@ -484,8 +478,6 @@ def run_while_convergence_iteration(
         )
         residual_voller = 1
         # Run the while loop until convergence is reached
-        # FIXME: tolerance check and change for new enthalpy update for inhomognoeus physical values
-        # TODO: Document Voller method on overleaf
         while (
             ((residual_voller > self.preprocess_data.temperature_tolerance) & (stefan))
             or (t_err > self.preprocess_data.temperature_tolerance)
@@ -816,23 +808,6 @@ def record_iteration_data(self):
             self.preprocess_data.s_k_iter,
         ) = [], [], [], [], [], [], []
 
-        # parameter_to_record_and_reset_list = [
-        #     [self.preprocess_data.t_k_iter, self.preprocess_data.t_k_iter_all],
-        #     [self.preprocess_data.phi_k_iter, self.preprocess_data.phi_k_iter_all],
-        #     [self.preprocess_data.all_phi_iter, self.preprocess_data.all_phi_iter_all],
-        #     [
-        #         self.preprocess_data.t_k_before_convergence,
-        #         self.preprocess_data.t_k_before_convergence_all,
-        #     ],
-        #     [
-        #         self.preprocess_data.mush_indx_list,
-        #         self.preprocess_data.mush_indx_list_all,
-        #     ],
-        # ]
-
-        # for parameter_list, parameter_all_list in parameter_to_record_and_reset_list:
-        #     self.check_and_reset_any_iteration_data(parameter_list, parameter_all_list)
-
     def initialize_state_variables(self, t, t_km1, s_km1, phi_km1):
         """Initializes the state variables for the sea ice model.
 
@@ -989,9 +964,7 @@ def run_sea_ice_model(self):
             np.array([]),
             self.preprocess_data.upwind_velocity,
         )
-        # with alive_bar(self.preprocess_data.max_iterations, force_tty=True) as bar:
         for t in range(1, self.preprocess_data.max_iterations):
-            # time.sleep(0.005)
             if self.preprocess_data.is_buffo:
                 (
                     t_k_buffo,
@@ -1093,9 +1066,6 @@ def run_sea_ice_model(self):
             t_stefan, error_depth_t, thickness_stefan = (
                 self.choose_phase_type_iteration(t)
             )
-
-            # TODO: Add radiation and algae
-            # TODO: create a parameter script for algae model
             (
                 radiative_source_term,
                 salinity_source_term,
@@ -1163,7 +1133,6 @@ def run_sea_ice_model(self):
                 salinity_source_term,
                 buffo=self.preprocess_data.is_buffo,
             )
-            # bar()
 
             if t % 1000 == 0:
                 count += 1

diff --git a/src/spyice/models/stefan_problem.py b/src/spyice/models/stefan_problem.py
@@ -2,15 +2,13 @@
 
 import matplotlib.pyplot as plt
 
-# TODO: check if the concentration is dependent on previous iterations or vice-versa
-# TODO: implement mushy layer consideration with liquid fraction check paper
 import numpy as np
 import scipy.optimize as opt
 from scipy.special import erfc
 
 from spyice.parameters.user_input import UserInput
 
-ui = UserInput()
+# ui = UserInput()
 
 
 class StefanProblem:
@@ -175,7 +173,7 @@ def stefan_problem_twophase(t, ui: UserInput):
         return 2 * lambda_stefan * np.sqrt(ui.constants.D_s * t)
 
 
-    def _plot_stefan_temp_twophase(z_depth=0.5):
+    def _plot_stefan_temp_twophase(ui, z_depth=0.5):
         dt = ui.grid_timestep_dt
         t_passed = 0
         temperature_array = []
@@ -203,5 +201,6 @@ def _plot_stefan_temp_twophase(z_depth=0.5):
         plt.show()
         return temperature_array, salinity_arr
 
-if __name__ == "__main__":
-    StefanProblem._plot_stefan_temp_twophase()
+# if __name__ == "__main__":
+#     ui = UserInput()
+#     StefanProblem._plot_stefan_temp_twophase(ui = ui)
diff --git a/src/spyice/parameters/user_input.py b/src/spyice/parameters/user_input.py
@@ -204,7 +204,7 @@ class UserInput:
 
     # --- Model Switches ---
     is_stefan: bool = True
-    is_buffo: bool = True
+    is_buffo: bool = False
     is_voller: bool = False
 
     is_salinity_equation: bool = True
@@ -224,8 +224,8 @@ class UserInput:
     # --- Boundary and Geometry ---
     boundary_condition_type: BoundaryConditionType = BoundaryConditionType.DIRICHLET.value
     geometry_type: int = field(init=False)
-    boundary_salinity: float = field(init=False)
-    boundary_top_temperature: float = field(init=False)
+    boundary_salinity: float = 34.0
+    boundary_top_temperature: float = 265.0
     temperature_top_type: TopTemperatureType = TopTemperatureType.STEFAN.value
 
     # --- Tolerances ---
@@ -263,8 +263,6 @@ def __post_init__(self):
         _dt_stability_validator(self.grid_resolution_dz, self.grid_timestep_dt)
 
         if isinstance(self.constants, RealConstants):
-            self.boundary_salinity = 34
-            self.boundary_top_temperature = 265
 
             # melt temperature affects the liquid relation: Frezchem or Normal in src/update_physical_values.py script
             method = InitialMeltTemperature.ONEPHASE.value

diff --git a/src/spyice/postprocess/visualise_model.py b/src/spyice/postprocess/visualise_model.py
@@ -448,7 +448,7 @@ def plot_temperature(
             + ".csv"
         )
 
-    def plot_temperature_heatmap(self, savefig: bool = True, export_csv=False, show: bool = True):
+    def plot_temperature_heatmap(self, savefig: bool = True, export_csv=False, show: bool = False):
         """Plots the temperature heatmap."""
 
         print("Plotting Temperature heatmap...")
@@ -499,9 +499,9 @@ def plot_temperature_heatmap(self, savefig: bool = True, export_csv=False, show:
         plt.close(fig)
 
     def plot_salinity_heatmap(self, savefig: bool = True, show: bool = False):
-        """Plots the temperature heatmap."""
+        """Plots the Salinity heatmap."""
 
-        print("Plotting Temperature heatmap...")
+        print("Plotting Salinity heatmap...")
         x_axis_iter = np.arange(0, self.ui_object.max_iterations - 1, 1)
         heatmap_data = self.results_object.s_k_list[1:]
         fig, ax = plt.subplots()
@@ -543,10 +543,10 @@ def plot_salinity_heatmap(self, savefig: bool = True, show: bool = False):
             plt.show()
         plt.close(fig)
 
-    def plot_liquidfraction_heatmap(self, savefig: bool = True):
-        """Plots the temperature heatmap."""
+    def plot_liquidfraction_heatmap(self, savefig: bool = True, show=False):
+        """Plots the Liquidfraction heatmap."""
 
-        print("Plotting Temperature heatmap...")
+        print("Plotting Liquid-Fraction heatmap...")
         x_axis_iter = np.arange(0, self.ui_object.max_iterations - 1, 1)
         heatmap_data = self.results_object.phi_k_list[1:]
         fig, ax = plt.subplots()
@@ -583,6 +583,8 @@ def plot_liquidfraction_heatmap(self, savefig: bool = True):
                 self.ui_object.dir_output_name + "/Liquidfraction.pdf",
                 backend="pgf",
             )
+        if show:
+            plt.show()
         plt.close(fig)
 
     def plot_temperature_heatmap_as_gif(self):
@@ -1806,90 +1808,4 @@ def plot_temperature_3d_contours(self, savefig=True):
                 self.ui_object.dir_output_name + "/temperature_3d_contour_profile.pdf",
                 backend="pgf",
             )
-        plt.close(fig)
-
-
-    # def residual_plot(self):
-    #     res_0p01_1000iter_voller1 = np.load("outputs/2024-09-26/123923_real_47.0_1000_0.01/residuals.npy", allow_pickle=True)
-
-    # TODO: Add the following methods
-    # def plot_phi(self, timestep, savefig=True):
-    #     print(f"Plotting Liquid Fraction evolution at {timestep}H...")
-    #     x_axis_iter = np.arange(0,self.nz,1)
-    #     #x_axis_iter = np.arange(0,22970,1)
-    #     index = timestep*3600/self.ui_object.grid_timestep_dt
-    #     phi_k = self.phi_k_list[int(index)]
-    #     if self.is_buffo is True:
-    #         phi_buffo = self.phi_buffo_list[int(index)]
-
-    #     fig1,(ax1) = plt.subplots(figsize=(10,6))
-    #     plt.grid()
-    #     ax1.plot(x_axis_iter*self.ui_object.grid_timestep_dt/3600, phi_k, 'r--',label='Numerical Temperature')
-    #     if self.is_buffo is True:
-    #         ax1.plot(x_axis_iter*self.ui_object.grid_timestep_dt/3600, phi_buffo, 'b--',alpha=0.5,label='Buffo')
-    #     ax1.set_xlabel("Depth nodes")
-    #     ax1.set_ylabel("Phi")
-    #     #ax1.legend()
-    #     color = 'gray'
-    #     ax1.legend()
-    #     fig1.tight_layout()
-    #     if savefig:
-    #         fig1.savefig(self.dir_output_name + "/Liquid Fraction evolution at"+ str(timestep) +"s.png")
-
-    # def plot_salinity(self, z_depth, savefig=True):
-    #     x_axis_iter = np.arange(0, iter_max-1,1)
-    #     index = z_depth*self.nz
-    #     S_k = self.S_k_list[:,int(index)]
-    #     if self.is_buffo is True:
-    #         S_buffo = self.S_buffo_list[:,int(index)]
-    #     fig1,(ax1) = plt.subplots(figsize=(10,6))
-    #     plt.grid()
-    #     ax1.plot(x_axis_iter*self.ui_object.grid_timestep_dt/3600, S_k[x_axis_iter], 'r--',label='Numerical Temperature')
-    #     if self.is_buffo is True:
-    #         ax1.plot(x_axis_iter*self.ui_object.grid_timestep_dt/3600, S_buffo[x_axis_iter], 'b--',alpha=0.5,label='Buffo')
-    #     ax1.set_xlabel("t in hours")
-    #     ax1.set_ylabel("Salinity in ppt")
-    #     #ax1.legend()
-    #     color = 'gray'
-    #     ax1.legend()
-    #     fig1.tight_layout()
-    #     if savefig:
-    #         fig1.savefig(self.dir_output_name + "/Salinity evolution at"+ str(z_depth) +"m.png")
-
-    # def plot_enthalpy(self, timestep, savefig=False):
-    #     print(f"Plotting Liquid Fraction Vs Temperature at {timestep}H...")
-    #     x_axis_iter = np.arange(0, iter_max-1,1)
-    #     index = timestep*3600/self.ui_object.grid_timestep_dt
-    #     mush = self.phi_slope(int(index))
-    #     phi = self.phi_k_list[int(index)]
-    #     H_k = self.H_k_list[int(index)]
-    #     H_solid = self.H_solid_list[int(index)]
-    #     T_k = self.T_k_list[int(index)]
-    #     H = (H_k - H_solid)/334774
-    #     T_melt = Tm_w
-    #     T_interface = self.Temp_interface[int(index)]
-    #     fig1,(ax1) = plt.subplots(figsize=(10,6))
-    #     plt.grid()
-    #     ax1.plot(T_k,phi, 'r--',label='Phi')
-    #     ax1.fill_betweenx(H, T_k[mush][0], T_k[mush][-1], color='gray', alpha=0.2, label='Mushy Layer')
-    #     ax1.set_xlabel("Temperature in K")
-    #     ax1.set_ylabel(r"Liquid Fraction $\phi$")
-    #     color = 'gray'
-    #     ax3 = ax1.twinx()
-    #     ax3.axvline(
-    #         T_melt,
-    #         color='b',
-    #         linestyle='dashed',
-    #         label=f'T_melt:{str(round(T_melt, 2))}',
-    #     )
-    #     ax3.axvline(
-    #         T_interface,
-    #         color='g',
-    #         linestyle='dashed',
-    #         label=f'T_interface:{str(round(T_interface, 2))}',
-    #     )
-    #     fig1.tight_layout()
-    #     ax1.legend(loc=0)
-    #     ax3.legend(loc=1)
-    #     if savefig:
-    #         fig1.savefig(self.dir_output_name + "/LiqVsTemp evolution at"+ str(timestep) +"h.png")
+        plt.close(fig)
diff --git a/src/spyice/preprocess/initial_boundary_conditions.py b/src/spyice/preprocess/initial_boundary_conditions.py
@@ -203,10 +203,11 @@ def set_boundary_temperature(t_passed, temperature_bottom, **kwargs):
     """
 
     top_temp = kwargs.get("top_temp")
+    top_temp_value = kwargs.get("top_temp_value")
     if top_temp == "T_const_250":
         temperature_top = 250.0
     if top_temp == "Stefan":
-        temperature_top = boundary_top_temperature
+        temperature_top = top_temp_value
         temperature_bottom = temperature_bottom
     elif top_temp == "T_const_260":
         temperature_top = 260.0