feat: randomize initial pressure, add assertions

1 files changed, 32 insertions, 13 deletions

diff --git a/src/main.py b/src/main.py
index 5e0d682..fffb570 100644
--- a/src/main.py
+++ b/src/main.py
@@ -8,7 +8,7 @@ figure, axes = plt.subplots()
 cb = None
 
 
-Re = 100
+Re = 400
 nu = 1 / Re
 domain_size = (1, 2)
 step = 0.05
@@ -32,9 +32,7 @@ v = np.zeros(shape=shape, dtype=float)
 v_star = np.zeros(shape=shape, dtype=float)
 v_new = np.zeros(shape=shape, dtype=float)
 
-p = np.zeros(shape=shape, dtype=float)
-p_star = np.zeros(shape=shape, dtype=float)
-p_new = np.zeros(shape=shape, dtype=float)
+p = np.random.rand(shape[0], shape[1])
 
 d_e = np.zeros(shape=shape, dtype=float)
 d_n = np.zeros(shape=shape, dtype=float)
@@ -46,19 +44,21 @@ def u_boundary(t, y):
         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) # > 0
 
 # Loop
 error = 1
 precision = 10 ** -7
 
-t = 0.35
+t = 0.4
 
 iteration = 0
 while error > precision:
     iteration += 1
     # Inflow boundary condition
     for i in range(N):
-        y = domain_size[0] - i * step
+        y = i * step
         u_star[i][0] = u_boundary(t, y)
         v_star[i][0] = 0
 
@@ -69,6 +69,13 @@ while error > precision:
         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[i][j] = 0
+
     # x-momentum
     for i in range(1, N - 1):
         for j in range(1, M - 1):
@@ -81,6 +88,10 @@ while error > precision:
             a_W = +0.5 * u_W * step + nu
             a_N = -0.5 * v_N * step + nu
             a_S = +0.5 * v_S * step + nu
+            assert_rule2(a_E)
+            assert_rule2(a_W)
+            assert_rule2(a_N)
+            assert_rule2(a_S)
 
             a_e = 0.5 * step * (u_E - u_W + v_N - v_S) + 4 * nu
             A_e = -step
@@ -102,6 +113,10 @@ while error > precision:
             a_W = +0.5 * u_W * step + nu
             a_N = -0.5 * v_N * step + nu
             a_S = +0.5 * v_S * step + nu
+            assert_rule2(a_E)
+            assert_rule2(a_W)
+            assert_rule2(a_N)
+            assert_rule2(a_S)
 
             a_n = 0.5 * step * (u_E - u_W + v_N - v_S) + 4 * nu
             A_n = -step
@@ -111,13 +126,6 @@ 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[i][j] - p[i + 1][j])
 
-    # 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_new[i][j] = 100
-
     # Pressure correction
     p_c = np.zeros(shape=shape, dtype=float)
     for i in range(1, N - 1):
@@ -126,6 +134,10 @@ while error > precision:
             a_W = -d_e[i][j-1] * step
             a_N = -d_n[i-1][j] * step
             a_S = -d_n[i][j] * step
+            assert_rule2(a_E)
+            assert_rule2(a_W)
+            assert_rule2(a_N)
+            assert_rule2(a_S)
             a_P = a_E + a_W + a_N + a_S
             b[i][j] = step * (-(u_star[i][j] - u_star[i][j-1]) + (v_star[i][j] - v_star[i-1][j]))
 
@@ -146,6 +158,13 @@ while error > precision:
             u_new[i][j] = u_star[i][j] + alpha * d_e[i][j] * (p_c[i + 1][j] - p_c[i][j])
             v_new[i][j] = v_star[i][j] + alpha * d_n[i][j] * (p_c[i][j] - p_c[i + 1][j])
 
+    # Backwards-facing step boundary conditions enforce
+    for i in range(int(h_c / step)):
+        for j in range(int(l_c / step)):
+            u_new[i][j] = 0
+            v_new[i][j] = 0
+            p_new[i][j] = 0
+
     # Continuity residual as error measure
     error = 0
     for i in range(N):