modified: run_simulation.py source/interpolation.py source/io_data_process.py source/momentum.py source/solifluction.py

heat_transfer does not work correct in the previous version as the surface temperature does not pass corectly from the file in time.
It is revised in the current version.

Author: saeb-faraji-gargari

run_simulation.py

@@ -57,15 +57,15 @@ def main() -> None:
         days_temperature_file,
         temps_temperature_file,
         results_pathname,
-        g_sin,
+        slope_radian,
     ) = read_config_file(param_path)
 
     # ---------------------  initial information --------------------
 
     dx, dz, array_shape, max_h_total = read_tif_info_from_gdal(
         h_total_initial_file_name
     )
-    num_rows, num_cols = array_shape
+    # num_rows, num_cols = array_shape
 
     partition_shape: tuple[int, int] = 2 * (partition_shape_size,)
 
@@ -81,36 +81,6 @@ def main() -> None:
 
     # ---------------------  initial information --------------------
 
-    # solifluction_simulate(
-    #     dx,
-    #     dz,
-    #     num_cols,
-    #     num_rows,
-    #     max_h_total,
-    #     bed_depth_elevation,
-    #     h_total_initial_file,
-    #     mu_value,
-    #     density_value,
-    #     k_conductivity_value,
-    #     rho_c_heat_value,
-    #     dt_momentum,
-    #     dt_mass_conservation,
-    #     dt_heat_transfer,
-    #     write_intervals_time,
-    #     momentum_iteration_threshold,
-    #     time_end_simulation,
-    #     heat_transfer_warmup,
-    #     heat_transfer_warmup_iteration,
-    #     h_mesh_step_value,
-    #     g_sin,
-    #     nu_x,
-    #     nu_z,
-    #     days_temperature_file,
-    #     temps_temperature_file,
-    #     partition_shape,
-    #     results_pathname,
-    # )
-
     solifluction(
         array_shape=array_shape,
         partition_shape=partition_shape,
@@ -133,7 +103,7 @@ def main() -> None:
         density_value=density_value,
         k_conductivity_value=k_conductivity_value,
         rho_c_heat_value=rho_c_heat_value,
-        g_sin=g_sin,
+        slope_radian=slope_radian,
         nu_x=nu_x,
         nu_z=nu_z,
         days_temperature_file=days_temperature_file,

source/interpolation.py

@@ -4,6 +4,8 @@ def interpolate_temperature(
     seconds_per_day = 86400
     t_days = t_seconds / seconds_per_day
 
+    # print("t_days: ", t_days)
+
     # handel the values out of input data boundaries
     if t_days <= days[0]:
         return temps[0]
@@ -12,7 +14,7 @@ def interpolate_temperature(
 
     # linear interpolation
     for i in range(len(days) - 1):
-        if days[i] <= t_days < days[i + 1]:
+        if days[i] <= t_days <= days[i + 1]:
             temp_lower = temps[i]
             temp_upper = temps[i + 1]
             frac = (t_days - days[i]) / (days[i + 1] - days[i])

source/io_data_process.py

@@ -221,7 +221,7 @@ def read_config_file(param_path: Path) -> tuple[
     list[float],  # days_temperature_file
     list[float],  # temps_temperature_file
     str,  # results_pathname
-    float,  # g_sin
+    float,  # slope_radian
 ]:
 
     input_variables: ConfigParser = load_config(param_path)
@@ -292,7 +292,7 @@ def clean_int(val: str) -> int:
         raise ValueError(f"Missing results path: {err}") from err
 
     slope_radian = clean_float(input_variables["simulation"]["alfa_slope"])
-    g_sin = np.sin(slope_radian) * 9.81
+    # g_sin = np.sin(slope_radian) * 9.81
 
     return (
         number_of_iterations,
@@ -315,7 +315,7 @@ def clean_int(val: str) -> int:
         days_temperature_file,
         temps_temperature_file,
         results_pathname,
-        g_sin,
+        slope_radian,
     )
 
 
@@ -371,11 +371,16 @@ def initiate_layers_variables(
     print("initiate_initial_layers_variables is in run")
 
     num_layers: int = int(
-        np.ceil(max_h_total - bed_depth_elevation)
-        / h_mesh_step_value
+        np.ceil((max_h_total - bed_depth_elevation) / h_mesh_step_value)
         # np.ceil(max_h_total - bed_depth_elevation) / h_mesh_step_value + 1
     )
 
+    print("max_h_total: ", max_h_total)
+    print("h_mesh_step_value: ", h_mesh_step_value)
+    print("bed_depth_elevation: ", bed_depth_elevation)
+
+    print("num_layers: ", num_layers)
+
     h_total_initial = lfr.from_gdal(
         h_total_initial_file, partition_shape=partition_shape
     )
@@ -526,15 +531,39 @@ def write_tif_file(
     lfr.to_gdal(array, output_path)
 
 
-def save_tif_file(
-    array: Any, output_file_name: str, time_label: float, output_folder_path: Path
-) -> None:
-    """Write a lue array to a .tif file using lfr.to_gdal."""
+# def save_tif_file(
+#     array: Any, output_file_name: str, time_label: float, output_folder_path: Path
+# ) -> None:
+#     """Write a lue array to a .tif file using lfr.to_gdal."""
 
-    os.makedirs(output_folder_path, exist_ok=True)
+#     os.makedirs(output_folder_path, exist_ok=True)
 
-    output_path: str = os.path.join(
-        output_folder_path, f"{output_file_name}-{time_label}.tif"
-    )
+#     output_path: str = os.path.join(
+#         output_folder_path, f"{output_file_name}-{time_label}.tif"
+#     )
 
-    lfr.to_gdal(array, output_path)
+#     lfr.to_gdal(array, output_path)
+
+
+def save_u_x_tem_time(u_x_tem_time, filename="u_x_tem_time.csv") -> None:
+    """
+    Save simulation results surface u_x and temperature in time to CSV file.
+
+    Parameters
+    ----------
+    filename : str
+        Output CSV file name.
+    """
+    data = np.array(u_x_tem_time)
+
+    print("Data shape:", data.shape)
+    print("Data type:", data.dtype)
+
+    np.savetxt(
+        filename,
+        data,
+        delimiter=",",
+        header="time,surface_temp,u_x_50_50",
+        comments="",
+        fmt="%.10e",
+    )

source/momentum.py

@@ -399,3 +399,201 @@ def momentum_ux(
     # )
 
     return phi
+
+
+# def momentum_ux_steady_state(
+#     phi: Any,
+#     phase_state: Any,
+#     dx: float,
+#     dz: float,
+#     lue_u_x: Any,
+#     lue_u_z: Any,
+#     nu_x: float,
+#     nu_z: float,
+#     rhs: Any,
+#     h_mesh: Any,
+#     boundary_loc: Any,
+#     boundary_type: Any,
+#     Dirichlet_boundary_value: Any,  # noqa: N803
+#     Neumann_boundary_value: Any,  # noqa: N803
+# ) -> Any:
+
+#     # lue_u_x : phi
+#     # lue_u_y : 0
+#     # nu_x : +(mu_mesh/gama_soil_mesh)
+#     # rhs_g_sin : g*sin(alfa)
+#     # gama_prime_mesh: (gama_soil/cos(alfa))-(gama_water*cos(alfa))
+#     # rhs : g*sin(alfa) - ((gama_soil/cos(alfa))-(gama_water*cos(alfa)))*dh/dx
+
+#     # kernel_i_j = np.array(
+#     #     [
+#     #         [0, 0, 0],
+#     #         [0, 1, 0],
+#     #         [0, 0, 0],
+#     #     ],
+#     #     dtype=np.uint8,
+#     # )
+
+#     dt: float = 1
+
+#     # kernel_im1_j   i-1, j
+#     kernel_im1_j: NDArray[np.uint8] = np.array(
+#         [
+#             [0, 0, 0],
+#             [1, 0, 0],
+#             [0, 0, 0],
+#         ],
+#         dtype=np.uint8,
+#     )
+
+#     """
+#     # kernel_i_jm1   i, j-1
+#     kernel_i_jm1 = np.array(
+#         [
+#             [0, 0, 0],
+#             [0, 0, 0],
+#             [0, 1, 0],
+#         ],
+#         dtype=np.uint8,
+#     )
+#     """
+
+#     # kernel_i_jm1   i, j-1
+#     kernel_i_jm1: NDArray[np.uint8] = np.array(
+#         [
+#             [0, 1, 0],
+#             [0, 0, 0],
+#             [0, 0, 0],
+#         ],
+#         dtype=np.uint8,
+#     )
+
+#     # kernel_ip1_j   i+1, j
+#     kernel_ip1_j: NDArray[np.uint8] = np.array(
+#         [
+#             [0, 0, 0],
+#             [0, 0, 1],
+#             [0, 0, 0],
+#         ],
+#         dtype=np.uint8,
+#     )
+
+#     """
+#     # kernel_i_jp1   i, j+1
+#     kernel_i_jp1 = np.array(
+#         [
+#             [0, 1, 0],
+#             [0, 0, 0],
+#             [0, 0, 0],
+#         ],
+#         dtype=np.uint8,
+#     )
+#     """
+
+#     # kernel_i_jp1   i, j+1
+#     kernel_i_jp1: NDArray[np.uint8] = np.array(
+#         [
+#             [0, 0, 0],
+#             [0, 0, 0],
+#             [0, 1, 0],
+#         ],
+#         dtype=np.uint8,
+#     )
+
+#     # coeff_map_i_j = (
+#     #     1
+#     #     + (-(dt / dx) * lfr.abs(lue_u_x))
+#     #     + (-(dt / dy) * lfr.abs(lue_u_y))
+#     #     + (2 * (dt / (dx**2)) * nu_x)
+#     #     + (2 * (dt / (dy**2)) * nu_y)
+#     # )
+
+#     # coeff_map_i_j: Any = (
+#     #     # 1  # check this   # This term is omitted since ∂u/∂t = 0 in the steady-state case
+#     #     +(-(dt / dx) * lfr.abs(lue_u_x))
+#     #     + (-(dt / dz) * lfr.abs(lue_u_z))
+#     #     + (2 * (dt / (dx**2)) * nu_x)
+#     #     + (2 * (dt / (dz**2)) * nu_z)
+#     # )
+
+#     coeff_map_i_j: Any = (
+#         (-(dt / dx) * lfr.abs(lue_u_x))
+#         + (-(dt / dz) * lfr.abs(lue_u_z))
+#         + (2 * (dt / (dx**2)) * nu_x)
+#         + (2 * (dt / (dz**2)) * nu_z)
+#     )
+
+#     coeff_map_im1_j: Any = lfr.where(
+#         lue_u_x >= 0,
+#         ((dt / dx) * lfr.abs(lue_u_x)) + (-(dt / (dx**2)) * nu_x),
+#         (-(dt / (dx**2)) * nu_x),
+#     )
+
+#     coeff_map_ip1_j: Any = lfr.where(
+#         lue_u_x < 0,
+#         ((dt / dx) * lfr.abs(lue_u_x)) + (-(dt / (dx**2)) * nu_x),
+#         (-(dt / (dx**2)) * nu_x),
+#     )
+
+#     coeff_map_i_jm1: Any = lfr.where(
+#         lue_u_z >= 0,
+#         ((dt / dz) * lfr.abs(lue_u_z)) + (-(dt / (dz**2)) * nu_z),
+#         (-(dt / (dz**2)) * nu_z),
+#     )
+
+#     coeff_map_i_jp1: Any = lfr.where(
+#         lue_u_z < 0,
+#         ((dt / dz) * lfr.abs(lue_u_z)) + (-(dt / (dz**2)) * nu_z),
+#         (-(dt / (dz**2)) * nu_z),
+#     )
+
+#     # NOTE: coeff_<.> should be implemented on boundaries, for instance on boundary_tpe
+#     # 1 (left boundary) phi_i-1,j is located outside of domain and we need
+#     # forward discretization
+#     # For now this implementation is ignored as we impose certain boundary conditions
+#     # on the boundaries which overwrite phi on the boundaries but in the future
+#     # this should be considered and for each boundary use exclusive discretization
+
+#     phi_all_internal_domain_i_j: Any = -(
+#         (coeff_map_i_j * phi)  # (coeff_map_i_j * lfr.focal_sum(solution, kernel_i_j))
+#         + (coeff_map_im1_j * lfr.focal_sum(phi, kernel_im1_j))
+#         + (coeff_map_i_jm1 * lfr.focal_sum(phi, kernel_i_jm1))
+#         + (coeff_map_ip1_j * lfr.focal_sum(phi, kernel_ip1_j))
+#         + (coeff_map_i_jp1 * lfr.focal_sum(phi, kernel_i_jp1))
+#         + (dt * rhs)
+#     )
+
+#     # phase_state: 0 solid  --> (frozen soil), 1 --> (fluid or unfrozen), now vegetation
+#     # is ignored in phase_state but it is considered in vegetation_vol_fraction
+
+#     phi_internal: Any = lfr.where(
+#         (phase_state != 0) & (h_mesh > 0),  # fluid or unfrozen
+#         phi_all_internal_domain_i_j,
+#         0,
+#     )
+
+#     phi = boundary_set(
+#         phi_internal,
+#         boundary_loc,
+#         boundary_type,
+#         Dirichlet_boundary_value,
+#         Neumann_boundary_value,
+#         dx,
+#         dz,
+#     )
+
+#     # ----------------
+
+#     # phi_u_x_numpy = lfr.to_numpy(phi)
+#     # print(
+#     #     "inside momentum_ux function phi_u_x_numpy[10,10]",
+#     #     phi_u_x_numpy[10, 10],
+#     # )
+
+#     # rhs_numpy = lfr.to_numpy(rhs)
+#     # print(
+#     #     "inside momentum_ux function rhs_numpy[10,10]",
+#     #     rhs_numpy[10, 10],
+#     # )
+
+#     return phi