feat: remove pressure boundaries, improve params

1 files changed, 22 insertions, 38 deletions

diff --git a/src/main.py b/src/main.py
index 5904fbb..ed363a0 100644
--- a/src/main.py
+++ b/src/main.py
@@ -8,20 +8,19 @@ cb = None
 plot = None
 
 
-Re = 200
+Re = 170
 nu = 1 / Re
 domain_size = (1, 2)
 step = 0.05
 N = int(domain_size[0] / step)
 M = int(domain_size[1] / step)
-
-
 shape = (N, M)
+
 alpha = 0.8  # Pressure under-relaxation coefficient
 
 t_m = 1
-h_c = 0.3
-l_c = 0.6
+h_c = 0.8
+l_c = 0.5
 
 # Staggered vars
 u = np.zeros(shape=shape, dtype=float)
@@ -38,11 +37,6 @@ d_n = np.zeros(shape=shape, dtype=float)
 b = np.zeros(shape=shape, dtype=float)
 
 
-def u_boundary(t, y):
-    def f(t):
-        return 1 if t_m < t else 0.5 * (np.sin(0.5 * PI * (2 * t / t_m - 1)) + 1)
-    return 6 * f(t) * (y - h_c) * (1 - y) / (1 - h_c)**2
-
 def assert_rule2(value):
     assert value > -0.01, f'Coefficient must be positive: {value}' # > 0
 
@@ -50,7 +44,7 @@ def assert_rule2(value):
 error = 1
 precision = 10 ** -7
 
-t = 0.4
+t = 0.8
 
 iteration = 0
 while error > precision:
@@ -58,23 +52,9 @@ while error > precision:
     # Inflow boundary condition
     for i in range(N):
         y = i * step
-        u_star[i][0] = u_boundary(t, y)
+        u_star[i][0] = 1
         v_star[i][0] = 0
 
-    # Sides boundary conditions
-    for j in range(M):
-        v_star[0][j] = 0
-        v_star[N - 1][j] = 0
-        u_star[0][j] = 0
-        u_star[N - 1][j] = 0
-
-    # Backwards-facing step boundary conditions (same as sides)
-    for i in range(int(h_c / step)):
-        for j in range(int(l_c / step)):
-            u_star[i][j] = 0
-            v_star[i][j] = 0
-            p_star[i][j] = 0
-
     # x-momentum
     for i in range(1, N - 1):
         for j in range(1, M - 1):
@@ -125,6 +105,19 @@ while error > precision:
             v_star[i][j] = (a_E * v[i][j + 1] + a_W * v[i][j - 1] + a_N * v[i - 1][j] + a_S * v[i + 1][j]) / a_n
             + d_n[i][j] * (p_star[i][j] - p_star[i + 1][j])
 
+    # Sides boundary conditions
+    for j in range(M):
+        v_star[0][j] = 0
+        v_star[N - 1][j] = 0
+        u_star[0][j] = 0
+        u_star[N - 1][j] = 0
+
+    # Backwards-facing step boundary conditions (same as sides)
+    for i in range(int(h_c / step)):
+        for j in range(int(l_c / step)):
+            u_star[i][j] = 0
+            v_star[i][j] = 0
+
     # Pressure correction
     p_prime = np.zeros(shape=shape, dtype=float)
     for i in range(1, N - 1):
@@ -144,14 +137,6 @@ while error > precision:
     p = p_star + p_prime * alpha
     p_star = p
 
-    # Pressure boundaries
-    for i in range(N - 1):
-        p[i][0] = p[i][1]
-        p[i][M - 1] = p[i][M - 2]
-    for j in range(M - 1):
-        p[0][j] = p[1][j]
-        p[N - 1][j] = p[N - 2][j]
-
     # Velocity correction
     for i in range(1, N - 1):
         for j in range(1, M - 1):
@@ -163,7 +148,6 @@ while error > precision:
         for j in range(int(l_c / step)):
             u[i][j] = 0
             v[i][j] = 0
-            p[i][j] = 0
 
     # Continuity residual as error measure
     new_error = 0
@@ -171,7 +155,7 @@ while error > precision:
         for j in range(M):
             new_error += abs(b[i][j])
     new_error /= N * M
-    if abs(new_error - error) > 10 ** -12:
+    if abs(new_error - error) > 10 ** -20:
         error = new_error
     else:
         print('Error decrease too small, you probably wont do better')
@@ -194,7 +178,7 @@ while error > precision:
         u_normalized = u / factor
         v_normalized = v / factor
 
-        density = 2
+        density = 1
         plot = plt.quiver(
             x[::density, ::density],
             y[::density, ::density],
@@ -202,7 +186,7 @@ while error > precision:
             v_normalized[::density, ::density],
             p[::density, ::density],
             scale=30,
-            cmap='inferno'
+            cmap='plasma'
         )
         cb = plt.colorbar()
         plt.pause(0.0001)