diff --git a/crs.py b/crs.py
index 5e2c14b512add932269dff487ecb728c1bce9c65..faf676419d5c7460d190cd8f09812f58e9508b10 100644
--- a/crs.py
+++ b/crs.py
@@ -22,9 +22,38 @@ def build_crs_matrix(n, row_callback):
- `row_ptr`: NumPy-array with `int` entries
"""
# TODO replace these
- val = np.empty((0,), dtype=float)
- col_idx = np.empty((0,), dtype=int)
- row_ptr = np.empty((0,), dtype=int)
+ elt = list()
+ nb_val = 0
+ row_ptr_lst = [0]
+ index = 1
+
+ # (*) change here if index starts with 1
+ last_id = 0
+
+ for i in range(n):
+ elt_l = row_callback(i)
+ if elt_l != list():
+ elt.append(elt_l)
+ nb_val += len(elt_l)
+
+ last_id += len(elt_l)
+
+ row_ptr_lst.append(last_id)
+ else:
+ row_ptr_lst.append(last_id)
+
+
+ col_idx = np.zeros((nb_val,), dtype=int)
+ row_ptr = np.array(row_ptr_lst, dtype=int)
+ val = np.zeros((nb_val,), dtype=float)
+
+
+ index = 0
+ for elt_i in elt:
+ for elt_j in elt_i:
+ val[index] = elt_j[0]
+ col_idx[index] = elt_j[1]
+ index +=1
return CrsMatrix(val, col_idx, row_ptr)
@@ -38,7 +67,49 @@ def crs_mvm(A, x):
Returns:
np.ndarray: The result vector y.
"""
- pass # TODO
+ # TODO
+ val, col_idx, row_ptr = A
+ n, = np.shape(x)
+ y = np.zeros((n,))
+ for i in range(n):
+ for j in range(row_ptr[i], row_ptr[i + 1]):
+ y[i] += val[j] * x[col_idx[j]]
+ return y
+
+def cg2(A, b, x0, epsilon=1e-6, max_iter=100):
+ """Runs the conjugate gradient algorithm on the system Ax = b, with starting vector x0, until the
+ norm of the residual falls below epsilon.
+
+ Args:
+ A (CrsMatrix): The system matrix, stored in CRS format
+ b (np.ndarray): Right-hand side of the system
+ x0 (np.ndarray): The starting vector of the iteration
+
+ Returns:
+ (np.ndarray, int): The approximate solution reached after `nmax` steps, or earlier if
+ a residual norm of <= epsilon was reached; and the number of steps the algorithm took.
+
+ """
+ # TODO
+ count = 0
+ x_n = x0
+ g_n = b - crs_mvm(A, x0)
+ s_n = -g_n
+ roh = np.dot(g_n, g_n) / (np.dot(s_n, crs_mvm(A, s_n)))
+ for i in range(max_iter):
+ x_n_plus = x_n + roh * s_n
+ g_n_plus = g_n + roh * crs_mvm(A, s_n)
+ beta = (np.dot(g_n_plus, g_n_plus)) / np.dot(g_n, g_n)
+ g_n = g_n_plus
+ s_n = -g_n - beta * s_n
+ roh = np.dot(g_n, g_n) / (np.dot(s_n, crs_mvm(A, s_n)))
+ if (np.linalg.norm(g_n) < epsilon):
+ return (g_n, count)
+ count += 1
+ return (g_n, count)
+
+
+
def cg(A, b, x0, epsilon=1e-6, max_iter=100):
"""Runs the conjugate gradient algorithm on the system Ax = b, with starting vector x0, until the
@@ -54,10 +125,36 @@ def cg(A, b, x0, epsilon=1e-6, max_iter=100):
a residual norm of <= epsilon was reached; and the number of steps the algorithm took.
"""
- pass # TODO
+ # TODO
+
+ n = len(b)
+ x = x0
+ g = b - crs_mvm(A, x)
+ s = g
+ gg_old = np.dot(g, g)
+ count = 0
+
+ for i in range(max_iter):
+ As = crs_mvm(A, s)
+ if (np.dot(s, As) < epsilon):
+ break
+ count += 1
+ alpha = gg_old / np.dot(s, As)
+ x = x + alpha * s
+ g = g - alpha * As
+ gg_new = np.dot(g, g)
+ if np.sqrt(gg_new) < epsilon:
+ break
+ s = g + (gg_new / gg_old) * s
+ gg_old = gg_new
+
+ return (x, count)
def main():
- pass
+ def a_cb1(j):
+ return [(1, j)]
+ val, col_idx, row_ptr = build_crs_matrix(4, a_cb1)
+ y = crs_mvm((val, col_idx, row_ptr), np.array([1, 2, 3, 4.0]))
if __name__ == "__main__":
main()
\ No newline at end of file
diff --git a/lbm.py b/lbm.py
index fcf4d8c1771fdbb0ab75f203e42b5ebcf11520fd..fac2b95f09d25ab0127f647886aed95686563979 100644
--- a/lbm.py
+++ b/lbm.py
@@ -28,7 +28,8 @@ def density(f):
A NumPy-array of dimension (W, H) that stores the density of
every cell.
"""
- pass # TODO
+ # TODO
+ return np.sum(f, axis=0)
def velocity(f):
@@ -44,7 +45,13 @@ def velocity(f):
Two NumPy-arrays of dimension (W, H) that represent the velocity
in x- and y-direciton.
"""
- pass # TODO
+ # TODO
+ d = density(f)
+ e_x = directions[:, 0]
+ e_y = directions[:, 1]
+ u_x = (np.sum(f * e_x[:, np.newaxis, np.newaxis], axis = 0)) / d
+ u_y = (np.sum(f * e_y[:, np.newaxis, np.newaxis], axis = 0)) / d
+ return (u_x, u_y)
def f_eq(f):
@@ -60,7 +67,33 @@ def f_eq(f):
A NumPy-array of dimension (9, W, H) that represents the
equilibrium distribution.
"""
- pass # TODO
+ # TODO
+ new_weights = weights[:,np.newaxis, np.newaxis]
+ d = density(f)
+ e_x = directions[:, 0]
+ e_x = e_x[:, np.newaxis, np.newaxis]
+ e_y = directions[:, 1]
+ e_y = e_y[:, np.newaxis, np.newaxis]
+ u_x = velocity(f)[0]
+ u_y = velocity(f)[1]
+ #expands d
+ d_exp = np.expand_dims(d, axis=0)
+ d_exp = np.repeat(d_exp, 9, axis=0)
+
+ #expand u_x
+ u_x_exp = np.expand_dims(u_x, axis=0)
+ u_x_exp = np.repeat(u_x_exp, 9, axis=0)
+
+ #expand u_y
+ u_y_exp = np.expand_dims(u_y, axis=0)
+ u_y_exp = np.repeat(u_y_exp, 9, axis=0)
+
+ tmp1 = 3 * (e_x * u_x_exp + e_y * u_y_exp)
+ tmp2 = (9 / 2) * ((tmp1 / 3) * (tmp1 / 3))
+ tmp3 = (3 / 2) * (u_x_exp * u_x_exp + u_y_exp * u_y_exp)
+ result = (new_weights * d_exp) * (1 + tmp1 + tmp2 - tmp3)
+ return result
+
def collide(f, omega):
@@ -79,7 +112,8 @@ def collide(f, omega):
A NumPy-array of dimension (9, W, H) that represents the state
after the collision operator has been applied.
"""
- pass # TODO
+ # TODO
+ return f - omega * (f - f_eq(f))
def stream(f):
@@ -95,7 +129,15 @@ def stream(f):
A NumPy-array of dimension (9, W, H) that represents the state
after the steam operator has been applied.
"""
- pass # TODO
+ # TODO
+ N, W, H = np.shape(f)
+ f_copy = np.copy(f)
+ for i in range(N):
+ for j in range(1, W-1):
+ for k in range(1, H-1):
+ f[i, j, k] = f_copy[i, j+directions[8-i, 0], k+directions[8-i, 1]]
+
+ return f
def noslip(f, masklist):
@@ -114,7 +156,13 @@ def noslip(f, masklist):
A NumPy-array of dimension (9, W, H) that represents the state
after the initialization of the noslip boundary conditions.
"""
- pass # TODO
+ # TODO
+ N, W, H = np.shape(f)
+ M, _ = np.shape(masklist)
+ f_copy = np.copy(f)
+ for i in range(N):
+ for l in range(M):
+ f[i, masklist[l, 0], masklist[l, 1]] = f_copy[inverse_dir[i], masklist[l, 0], masklist[l, 1]]
return f
diff --git a/poisson.py b/poisson.py
index 82f4dbeb269db833e846175a288e6813bb7512dc..071086504a0627a20724d5b06ba22b4d80a5e015 100644
--- a/poisson.py
+++ b/poisson.py
@@ -14,10 +14,29 @@ def generate_poisson_matrix(N):
Returns:
CrsMatrix: The system matrix in CRS format
"""
+ h = 1 / (N - 1) # Step size
+ num_inner_points = (N - 2) ** 2 # Number of inner points
+
def row_callback(row):
- pass # TODO
-
- pass # TODO
+ i = row // (N - 2) + 1
+ j = row % (N - 2) + 1
+ #index = i * N + j - 1 # -1 added
+ index = (i - 1) * (N - 2) + (j - 1)
+ #index = (i - 1) * (N - 2) + (j - 1)
+ neighbors = [(i, j - 1), (i, j + 1), (i - 1, j), (i + 1, j)]
+
+ entries = [(4 / h ** 2, index)] # Diagonal entry
+
+ for (ni, nj) in neighbors:
+ if 1 <= ni <= N - 2 and 1 <= nj <= N - 2:
+ #neighbor_index = (ni - 1) * (N - 2) + (nj - 1)
+ neighbor_index = (ni - 1) * (N - 2) + (nj - 1)
+ #neighbor_index = ni * N + nj - 1 # -1 added
+ entries.append((-1 / h ** 2, neighbor_index)) # Neighbor entries
+ entries.sort(key=lambda x: x[1])
+ return entries
+
+ return build_crs_matrix(num_inner_points, row_callback)
def generate_rhs(N, f, phi):
"""Generates the right-hand side vector for the Poisson problem.
@@ -31,7 +50,38 @@ def generate_rhs(N, f, phi):
Returns:
np.ndarray: Right-hand side vector as one-dimensional numpy array.
"""
- pass # TODO
+ # TODO
+ num_inner_points = (N - 2) ** 2 # Number of inner points
+ b = np.zeros(num_inner_points) # Initialize the right-hand side vector
+
+ h = 1 / (N - 1) # Step size
+
+ # Fill in the values for the inner points
+ for row in range(1, N - 1):
+ for col in range(1, N - 1):
+ x = col * h
+ y = row * h
+ index = (row - 1) * (N - 2) + (col - 1)
+ b[index] = f(x, y)
+
+ # Handle the boundary points
+ for row in range(1, N - 1):
+ x0 = 0
+ y = row * h
+ b[(row - 1) * (N - 2)] += phi(x0, y)
+
+ xN = 1
+ b[row * (N - 2) - 1] += phi(xN, y)
+
+ for col in range(1, N - 1):
+ x = col * h
+ y0 = 0
+ b[col - 1] += phi(x, y0)
+
+ yN = 1
+ b[(N - 3) * (N - 2) + col - 1] += phi(x, yN)
+
+ return b
def solve_poisson(N, f, phi):
"""Solves the discretized Poisson problem on an N x N square grid with Dirichlet boundary condition
@@ -47,7 +97,15 @@ def solve_poisson(N, f, phi):
np.ndarray: Two-dimensional array of shape (n, n) containing the values of the discrete solution u,
with u_ij stored at index [j, i].
"""
- pass # TODO
+ # TODO
+ A = generate_poisson_matrix(N)
+ b = generate_rhs(N, f, phi)
+ u = cg(A, b)
+
+ # Reshape the solution vector u to a grid
+ u_grid = u.reshape((N - 2, N - 2))
+
+ return u_grid
def main():
import matplotlib.pyplot as plt