modified: source/io_data_process.py source/solifluction.py source/vof.py

number of layers calculation is revised to avoid a layer with zero thikness at surface, rhs is set not to be negative

Author: saeb-faraji-gargari

source/io_data_process.py

@@ -356,7 +356,9 @@ 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 + 1
+        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
     )
 
     h_total_initial = lfr.from_gdal(

source/solifluction.py

@@ -340,9 +340,9 @@ def solifluction_simulate(
 
             print("local_heat_transfer_iteration :", local_heat_transfer_iteration)
 
-        # --------------- End: compute temperatures in internal layers -------------------------
+        # --------------- End: compute temperatures in internal layers -----------------
 
-        # --------------- compute momentum u_x in internal layers ------------------------------
+        # --------------- compute momentum u_x in internal layers ----------------------
 
         for local_momentum_iteration in range(1, momentum_iteration_threshold + 1):
             # while (
@@ -370,6 +370,13 @@ def solifluction_simulate(
                     )
                 )
 
+                # NOTE: rhs cannot be negative as the flow cannot be to the upstream
+                rhs = lfr.where(
+                    (rhs > 0),
+                    rhs,
+                    0,
+                )
+
                 # print("type of layer_list[layer_id].T:", type(layer_list[layer_id].T))
 
                 phase_state = phase_detect_from_temperature(layer_list[layer_id].T)
@@ -518,6 +525,16 @@ def solifluction_simulate(
                     layer_id,
                 )
 
+            for layer_id in range(0, num_layers):
+
+                layer_h_mesh_numpy = lfr.to_numpy(layer_list[layer_id].h_mesh)
+                print(
+                    "layer_h_mesh_numpy[10,10]",
+                    layer_h_mesh_numpy[10, 10],
+                    "layer_id",
+                    layer_id,
+                )
+
             # local_momentum_time = local_momentum_time + dt_momentum
             # local_momentum_iteration = local_momentum_iteration + 1
 
@@ -529,10 +546,10 @@ def solifluction_simulate(
             #         time - (momentum_iteration_threshold * dt_momentum)
             #     ) + dt_mass_conservation
 
-            print("time inside momentum u_x", time)
-            # input("enter to continue ...")
+        print("time inside momentum u_x", time)
+        input("enter to continue ...")
 
-        # --------------- End: compute momentum u_x in internal layers -------------------------------
+        # --------------- End: compute momentum u_x in internal layers -----------------
 
         # --------------- compute VOF -------------------------------
 

source/vof.py

@@ -268,6 +268,8 @@ def h_mesh_assign(
     h_total: Any, num_layers: Any, prespecified_uniform_h_mesh_value: float
 ) -> list[Any]:
 
+    epsilon = 1e-5
+
     h_total_remain = h_total
     h_mesh_list = []
 
@@ -287,6 +289,7 @@ def h_mesh_assign(
             h_total_remain,
             prespecified_uniform_h_mesh_value,
         )
+        h_mesh = lfr.where(h_mesh < epsilon, 0.0, h_mesh)
         h_mesh_list.append(h_mesh)
 
         h_total_remain = h_total_remain - prespecified_uniform_h_mesh_value