diff --git a/.github/workflows/build_docs.yml b/.github/workflows/build_docs.yml
index 547cbc1..77685f1 100644
--- a/.github/workflows/build_docs.yml
+++ b/.github/workflows/build_docs.yml
@@ -49,7 +49,7 @@ jobs:
       
       - name: Upgrade pip and setuptools
         run: |
-          python3 -m pip install --break-system-packages pip setuptools --upgrade
+          python3 -m pip install pip setuptools --upgrade
           
       # See: https://github.com/marketplace/actions/setup-miniconda
       - name: Setup Conda environment
diff --git a/src/vasp/simulations/simulation_common.py b/src/vasp/simulations/simulation_common.py
index e65ee16..3768c83 100644
--- a/src/vasp/simulations/simulation_common.py
+++ b/src/vasp/simulations/simulation_common.py
@@ -272,7 +272,6 @@ def calculate_and_print_flow_properties(dt: float, mesh: Mesh, v: Function, inle
     # Calculate the DG vector of velocity magnitudes
     DG = FunctionSpace(mesh, "DG", 0)
     V_vector = local_project(sqrt(inner(v, v)), DG, local_rhs).vector().get_local()
-    h = mesh.hmin()
 
     # Calculate flow rate at the inlet
     flow_rate_inlet = abs(assemble(inner(v, n) * dsi))
@@ -288,6 +287,10 @@ def calculate_and_print_flow_properties(dt: float, mesh: Mesh, v: Function, inle
     V_min = comm.gather(local_V_min, 0)
     V_max = comm.gather(local_V_max, 0)
 
+    # compute the minimum cell diameter in the mesh
+    h = mesh.hmin()
+    h_min = MPI.min(MPI.comm_world, h)
+
     if MPI.rank(comm) == 0:
         # Calculate mean, min, and max velocities
         v_mean = np.mean(V_mean)
@@ -301,9 +304,9 @@ def calculate_and_print_flow_properties(dt: float, mesh: Mesh, v: Function, inle
         Re_max = rho_f * v_max * diam_inlet / mu_f
 
         # Calculate CFL numbers
-        CFL_mean = v_mean * dt / h
-        CFL_min = v_min * dt / h
-        CFL_max = v_max * dt / h
+        CFL_mean = v_mean * dt / h_min * v.ufl_element().degree()
+        CFL_min = v_min * dt / h_min * v.ufl_element().degree()
+        CFL_max = v_max * dt / h_min * v.ufl_element().degree()
 
         # Print the flow properties
         print("Flow Properties:")
diff --git a/tests/test_simulations.py b/tests/test_simulations.py
index 679fa83..1fe0c6f 100644
--- a/tests/test_simulations.py
+++ b/tests/test_simulations.py
@@ -121,7 +121,7 @@ def test_cylinder_problem(input_mesh, tmpdir):
     output_match_cfl_number = re.findall(output_cfl_number, str(result))
     assert output_match_cfl_number is not None, "Regular expression did not match the output."
 
-    expected_cfl_number = [0.008465020210929876, 1.5793871333017697e-05, 0.022455733137609308]
+    expected_cfl_number = [0.016930040421859752, 3.1587742666035394e-05, 0.044911466275218616]
     cfl_number_mean_min_max = [float(output_match_cfl_number[-1][0]), float(output_match_cfl_number[-1][1]),
                                float(output_match_cfl_number[-1][2])]
 
@@ -187,7 +187,7 @@ def test_aneurysm_problem(input_mesh, tmpdir):
     output_match_cfl_number = re.findall(output_cfl_number, str(result))
     assert output_match_cfl_number is not None, "Regular expression did not match the output."
 
-    expected_cfl_number = [0.002380256687906308, 2.217345370176909e-17, 0.009231963373337334]
+    expected_cfl_number = [0.004760513375812616, 4.434690740353818e-17, 0.01846392674667467]
     cfl_number_mean_min_max = [float(output_match_cfl_number[-1][0]), float(output_match_cfl_number[-1][1]),
                                float(output_match_cfl_number[-1][2])]