Commit 05b6dd0a authored by Murali's avatar Murali
Browse files

Initial commit

parent 48fa3ca6

README.md

+382 −90

File changed.

Preview size limit exceeded, changes collapsed.

func_2DCylinder.py

0 → 100644
+317 −0
Original line number Diff line number Diff line 
import numpy as np

from skimage import draw

from skimage.draw import polygon



def make_grid(Lx, Ly, nx, ny):

    """Make a rectangular grid, with coordinates centered on (0, 0)."""

    x = np.linspace(0, Lx, num=nx) - Lx/2

    y = np.linspace(0, Ly, num=ny) - Ly/2

    X, Y = np.meshgrid(x, y)

    return X, Y





# ## Find the discrete boundary on the grid we have defined

def transform_grid_coords_to_xy(grid_coords, grid):

    """Transform rc coords to xy coords from the grid."""

    X, Y = grid

    x_disc = np.array([X[rc[0], rc[1]] for rc in grid_coords])

    y_disc = np.array([Y[rc[0], rc[1]] for rc in grid_coords])

    return x_disc, y_disc



def rasterize_points_to_grid_coords(x, y, grid):

    """Rasterizes a discrete set of (x, y) points to a grid. Returns (r, c) coordinates."""

    X, Y = grid

    X_flat, Y_flat = X.flatten(), Y.flatten()



    r = np.arange(X.shape[0])

    c = np.arange(X.shape[1])

    C, R = np.meshgrid(c, r)

    R_flat, C_flat = R.flatten(), C.flatten()



    coords = np.c_[X_flat, Y_flat]

    coords_rc = np.c_[R_flat, C_flat]

    

    # first part: discretize each point from x and y while avoiding duplicates

    rcs = []

    discrete_coords = []

    already_seen_points = set()

    for xx, yy in zip(x, y):

        dist = ((coords - np.array([xx, yy]).reshape(1, -1))**2).sum(axis=1)

        argmin = dist.argmin()

        if argmin not in already_seen_points:            

            rc = coords_rc[argmin]

            rcs.append(rc)

            discrete_coords.append(coords[argmin])

            already_seen_points.add(argmin)



    discrete_coords = np.array(discrete_coords)

    x_disc, y_disc = transform_grid_coords_to_xy(rcs, grid)

    

    # sanity check

    assert np.allclose(np.c_[x_disc, y_disc], discrete_coords)



    # second part: use the Bresenham algorithm to discretize each segment

    N = len(rcs)

    grid_coords = []

    for start, stop in zip(range(N), range(1, N+1)):

        if stop >= N:

            stop = stop % N

        r0, c0 = rcs[start]

        r1, c1 = rcs[stop]

        rr, cc = draw.line(r0, c0, r1, c1)

        grid_coords.extend([r,c] for r, c in zip(rr[:-1], cc[:-1]))

        

    return np.array(grid_coords)





# # The normal vectors on the interior boundary

# 

# Since we have now establised a clear link between the parametrized curve and the grid, we can move on to the computation of normal vectors at each grid point.

def rasterize_vectors_to_grid(x, y, nx, ny, grid):

    """Rasterizes a discrete set of (x, y) points and vectors (nx, ny) to a grid. 

    Returns (nx_disc, ny_disc) vector field."""



    X, Y = grid

    X_flat, Y_flat = X.flatten(), Y.flatten()



    r = np.arange(X.shape[0])

    c = np.arange(X.shape[1])

    C, R = np.meshgrid(c, r)

    R_flat, C_flat = R.flatten(), C.flatten()



    coords = np.c_[X_flat, Y_flat]

    coords_rc = np.c_[R_flat, C_flat]

    

    # first part: discretize each point from x and y while avoiding duplicates

    rcs = []

    discrete_coords = []

    discrete_vectors = []

    already_seen_points = set()

    for xx, yy, nxx, nyy in zip(x, y, nx, ny):

        dist = ((coords - np.array([xx, yy]).reshape(1, -1))**2).sum(axis=1)

        argmin = dist.argmin()

        if argmin not in already_seen_points:            

            rc = coords_rc[argmin]

            rcs.append(rc)

            discrete_coords.append(coords[argmin])

            discrete_vectors.append([nxx, nyy])

            already_seen_points.add(argmin)



    discrete_coords = np.array(discrete_coords)

    x_disc, y_disc = transform_grid_coords_to_xy(rcs, grid)

    

    # sanity check

    assert np.allclose(np.c_[x_disc, y_disc], discrete_coords)



    # second part: use the Bresenham algorithm to discretize each segment

    N = len(rcs)

    grid_coords = []

    rasterized_vector_field = []

    for start, stop in zip(range(N), range(1, N+1)):

        if stop >= N:

            stop = stop % N

        r0, c0 = rcs[start]

        r1, c1 = rcs[stop]

        nx0, ny0 = discrete_vectors[start]

        nx1, ny1 = discrete_vectors[stop]

        

        rr, cc = draw.line(r0, c0, r1, c1)

        N_interp = len(rr) 

        for i in range(N_interp - 1):

            alpha = i/N_interp

            grid_coords.append([rr[i], cc[i]])

            rasterized_vector_field.append([(1 - alpha) * nx0 + alpha * nx1,

                                            (1 - alpha) * ny0 + alpha * ny1])

    return np.array(rasterized_vector_field)





# # Setting up the finite difference problem geometry

# 

# Let's first distinguish the three regions we need to set up the problem:

# 

# * the outer boundary

# * the inner boundary

# * the rest of the region

def make_tags(grid, x, y):

    """Create a grid with tags defining each cell's type."""

    X, Y = grid

    r = np.arange(X.shape[0])

    c = np.arange(X.shape[1])

    C, R = np.meshgrid(c, r)



    xmin, xmax = X.min(), X.max()

    ymin, ymax = Y.min(), Y.max()

    EPS = 1e-5



    grid_coords = rasterize_points_to_grid_coords(x, y, grid)

    grid_coords_set = set((r,c) for (r,c) in grid_coords)

    mask = polygon(np.array(grid_coords)[:, 0], np.array(grid_coords)[:, 1], X.shape)

    mask_set = set((r, c) for r, c in zip(*mask))



    tags = np.zeros_like(R, dtype=int) * np.nan



    for r, c in zip(R.flatten(), C.flatten()):

        # left boundary

        if abs(X[r, c] - xmin) < EPS:

            tags[r, c] = 1

        # right boundary

        elif abs(X[r, c] - xmax) < EPS:

            tags[r, c] = 2

        # bottom boundary

        elif abs(Y[r, c] - ymin) < EPS:

            tags[r, c] = 3

        # top boundary

        elif abs(Y[r, c] - ymax) < EPS:

            tags[r, c] = 4

        # interior

        elif (r, c) in mask_set:

            # boundary

            if (r, c) in  grid_coords_set:

                tags[r, c] = 5

            # volume inside

            else:

                tags[r, c] = 6

        else:

            tags[r, c] = 7



    for r,c in grid_coords:

        tags[r, c] = 5

    return tags





# # Doing the equations: constant flow case

# 

# As a first step, we will model the potential value at each grid point and do as if there was no interior region.

def assemble_only(grid, x, y, nx, ny, tags):

    """Assembles and solves the potential flow from a curve (x, y, nx, ny) and on the given grid."""

    grid_coords = rasterize_points_to_grid_coords(x, y, grid)

    rasterized_vector_field = rasterize_vectors_to_grid(x, y, nx, ny, grid)



    assert grid_coords.shape[0] == rasterized_vector_field.shape[0]



    rc2normal = {}

    for rc, n in zip(grid_coords, rasterized_vector_field):

        rc2normal[(rc[0], rc[1])] = n



    X, Y = grid

    r = np.arange(X.shape[0])

    c = np.arange(X.shape[1])

    C, R = np.meshgrid(c, r)



    has_unknown = (tags != 6)

    rcs = np.c_[R[has_unknown].flatten(), C[has_unknown].flatten()]



    # Let's create some mappings to simplify the loops.



    rc2id = {}

    id2rc = {}



    for ind, (r, c) in enumerate(rcs):

        rc2id[(r, c)] = ind

        id2rc[ind] = (r, c)



    # coef matrix

    A = np.zeros((rcs.shape[0], rcs.shape[0]))

    # rhs vector

    B = np.zeros((rcs.shape[0],))

    # v0

    v0 = np.array([1., 0])

    # grid params

    dx = X[0, 1] - X[0, 0]

    dy = Y[1, 0] - Y[0, 0]



    for (r, c) in rc2id:

        this_point = rc2id[(r, c)]

        # left edge

        if tags[r, c] == 1:

            id_right = rc2id[(r, c+1)]

            id_left = rc2id[(r, c)]

            A[this_point, id_left] += 1/dx

            A[this_point, id_right] += -1/dx

            B[this_point] = -v0[0]

        # right edge

        elif tags[r, c] == 2:

            id_right = rc2id[(r, c)]

            id_left = rc2id[(r, c-1)]

            A[this_point, id_right] += 1/dx

            A[this_point, id_left] += -1/dx

            B[this_point] = v0[0]

        # top edge

        elif tags[r, c] == 4:

            id_top = rc2id[(r, c)]

            id_bottom = rc2id[(r-1, c)]

            A[this_point, id_top] += 1/dy

            A[this_point, id_bottom] += -1/dy

            B[this_point] = 0

        # bottom edge

        elif tags[r, c] == 3:

            id_top = rc2id[(r, c)]

            id_bottom = rc2id[(r+1, c)]

            A[this_point, id_top] += 1/dy

            A[this_point, id_bottom] += -1/dy

            B[this_point] = 0

        # point on the interior boundary

        elif tags[r, c] == 5:

            nxx, nyy = rc2normal[(r, c)]

            id_center = rc2id[(r, c)]

            # vertical gradient

            if (r-1, c) in rc2id:

                id_bottom = rc2id[(r-1, c)]

                A[this_point, id_center] += 1/(dy) * nyy

                A[this_point, id_bottom] += -1/(dy) * nyy

            else:

                id_top = rc2id[(r+1, c)]

                A[this_point, id_center] += -1/(dy) * nyy

                A[this_point, id_top] += 1/(dy) * nyy



            # horizontal gradient

            if (r, c+1) in rc2id:

                id_right = rc2id[(r, c+1)]

                A[this_point, id_right] += 1/(dx) * nxx

                A[this_point, id_center] += -1/(dx) * nxx

            else:

                id_left = rc2id[(r, c-1)]

                A[this_point, id_center] += 1/(dx) * nxx

                A[this_point, id_left] += -1/(dx) * nxx



        # inside the volume

        elif tags[r, c] == 6:

            id_center = rc2id[(r, c)]

            A[this_point, id_center] = 1

            B[this_point] = 0

        else:

            id_bottom = rc2id[(r-1, c)]

            id_top = rc2id[(r+1, c)]

            id_right = rc2id[(r, c+1)]

            id_left = rc2id[(r, c-1)]

            id_center = rc2id[(r, c)]

            A[this_point, id_top] += -1/dy

            A[this_point, id_bottom] += -1/dy

            A[this_point, id_right] += -1/dx

            A[this_point, id_left] += -1/dx

            A[this_point, id_center] += (2/dx + 2/dy)



    print(f'Size of matrices: A = {A.shape}; B = {B.shape}')

    return A, B, id2rc, X, Y



def compute_analytical_solution(grid, U, R):

    """Applies analytical solution of cylinder flow to every point in the grid."""

    X, Y = grid

    r = np.sqrt(X**2 + Y**2)

    theta = np.arctan2(Y, X)

    phi = U * r * (1 + R**2/r**2) * np.cos(theta) * (r > R)

    phi -= phi.min() # normalizing by minimum value

    v_r = U * (1 - R**2 / r**2) * np.cos(theta) * (r > R)

    v_theta = -U * (1 + R**2 / r**2) * np.sin(theta) * (r > R)

    v_x = np.cos(theta) * v_r - np.sin(theta) * v_theta

    v_y = np.sin(theta) * v_r + np.cos(theta) * v_theta



    return phi, v_x, v_y



def compute_gradient(grid, phi):

    """Numerically computes velocity as gradient of phi on grid."""

    X, Y = grid

    dx = X[0, 1] - X[0, 0]

    dy = Y[1, 0] - Y[0, 0]

    v, u = np.gradient(phi, dx, dy)

    return u, v
 
 No newline at end of file

func_HeleShaw.py

0 → 100644
+174 −0
Original line number Diff line number Diff line 
import numpy as np



def ij2idx(row_i, column_j, column):

    return (row_i*column) + column_j



def Laplacian_pressure(P, P_left, P_right, dx, dy):

    ny, nx = P.shape

    LP = np.zeros((ny*nx, ny*nx))

    b  = np.zeros((ny*nx))

    for i in range(ny):

        for j in range(nx):

            ii = ij2idx(i, j, nx)

            if j==0: # left boundary middle cells

                # forward in x and central in y

                LP[ii,ij2idx(i  ,j  , nx)] += 1

                b[ii] = P_left

            elif j==nx-1: # right boundary middle cells

                # backward in x and central in y

                LP[ii,ij2idx(i  ,j  , nx)] += 1

                b[ii] = P_right

            elif i==0: # bottom boundary

                if j==0: # left-bottom corner

                    # forward in x and forward in y

                    LP[ii,ij2idx(i  ,j  , nx)] += 1

                    b[ii] = P_left

                elif j==nx-1: # right-bottom corner

                    # backward in x and forward in y

                    LP[ii,ij2idx(i  ,j  , nx)] += 1

                    b[ii] = P_right

                else: # bottom middle cells

                    # central in x and forward in y

                    LP[ii,ij2idx(i  ,j  , nx)] += -2/(dx**2)

                    LP[ii,ij2idx(i  ,j-1, nx)] += 1/(dx**2)

                    LP[ii,ij2idx(i  ,j+1, nx)] += 1/(dx**2)

                    LP[ii,ij2idx(i  ,j  , nx)] += 1/(dy**2)

                    LP[ii,ij2idx(i+1,j  , nx)] += -2/(dy**2)

                    LP[ii,ij2idx(i+2,j  , nx)] += 1/(dy**2)

                    b[ii] = 0

            elif i==ny-1: # top boundary

                if j==0: # left-top corner

                    # forward in x and backward in y

                    LP[ii,ij2idx(i  ,j  , nx)] += 1

                    b[ii] = P_left

                elif j==nx-1: # right-top corner

                    # backward in x and backward in y

                    LP[ii,ij2idx(i  ,j  , nx)] += 1

                    b[ii] = P_right

                else: # top middle cells

                    # central in x and backward in y

                    LP[ii,ij2idx(i  ,j  , nx)] += -2/(dx**2)

                    LP[ii,ij2idx(i  ,j-1, nx)] += 1/(dx**2)

                    LP[ii,ij2idx(i  ,j+1, nx)] += 1/(dx**2)

                    LP[ii,ij2idx(i  ,j  , nx)] += 1/(dy**2)

                    LP[ii,ij2idx(i-1,j  , nx)] += -2/(dy**2)

                    LP[ii,ij2idx(i-2,j  , nx)] += 1/(dy**2)

                    b[ii] = 0

            else: # middle cells

                # central in x and central in y

                LP[ii,ij2idx(i  ,j  , nx)] += -2/(dx**2)

                LP[ii,ij2idx(i  ,j-1, nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j+1, nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j  , nx)] += -2/(dy**2)

                LP[ii,ij2idx(i-1,j  , nx)] += 1/(dy**2)

                LP[ii,ij2idx(i+1,j  , nx)] += 1/(dy**2)

                b[ii] = 0

    return LP, b





def Laplacian_velocity_x(U, U_top, U_bottom, P, dx, dy):

    ny, nx = U.shape

    LP = np.zeros((ny*nx, ny*nx))

    b  = np.zeros((ny*nx))

    for i in range(ny):

        for j in range(nx):

            ii = ij2idx(i, j, nx)

            if i==0: # bottom boundary

                LP[ii,ij2idx(i  ,j  , nx)] += 1

                b[ii] = U_bottom

            elif i==ny-1: # top boundary

                LP[ii,ij2idx(i  ,j  , nx)] += 1

                b[ii] = U_top

            elif j==0: # left boundary middle cells

                # forward in x and central in y

                LP[ii,ij2idx(i  ,j  , nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j+1, nx)] += -2/(dx**2)

                LP[ii,ij2idx(i  ,j+2, nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j  , nx)] += -2/(dy**2)

                LP[ii,ij2idx(i-1,j  , nx)] += 1/(dy**2)

                LP[ii,ij2idx(i+1,j  , nx)] += 1/(dy**2)

                px_ii = (P[i,j+1] - P[i,j]) / dx

                b[ii] = px_ii

            elif j==nx-1: # right boundary middle cells

                # backward in x and central in y

                LP[ii,ij2idx(i  ,j  , nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j-1, nx)] += -2/(dx**2)

                LP[ii,ij2idx(i  ,j-2, nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j  , nx)] += -2/(dy**2)

                LP[ii,ij2idx(i-1,j  , nx)] += 1/(dy**2)

                LP[ii,ij2idx(i+1,j  , nx)] += 1/(dy**2)

                px_ii = (P[i,j] - P[i,j-1]) / dx

                b[ii] = px_ii

            else: # middle cells

                # central in x and central in y

                LP[ii,ij2idx(i  ,j  , nx)] += -2/(dx**2)

                LP[ii,ij2idx(i  ,j-1, nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j+1, nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j  , nx)] += -2/(dy**2)

                LP[ii,ij2idx(i-1,j  , nx)] += 1/(dy**2)

                LP[ii,ij2idx(i+1,j  , nx)] += 1/(dy**2)

                px_ii = 0.5 * (P[i,j+1] - P[i,j-1]) / dx

                b[ii] = px_ii

    return LP, b



def Laplacian_velocity_y(V, V_top, V_bottom, P, dx, dy):

    ny, nx = V.shape

    LP = np.zeros((ny*nx, ny*nx))

    b  = np.zeros((ny*nx))

    for i in range(ny):

        for j in range(nx):

            ii = ij2idx(i, j, nx)

            if i==0: # bottom boundary

                LP[ii,ij2idx(i  ,j  , nx)] += 1

                b[ii] = V_bottom

            elif i==ny-1: # top boundary

                LP[ii,ij2idx(i  ,j  , nx)] += 1

                b[ii] = V_top

            elif j==0: # left boundary middle cells

                # forward in x and central in y

                LP[ii,ij2idx(i  ,j  , nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j+1, nx)] += -2/(dx**2)

                LP[ii,ij2idx(i  ,j+2, nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j  , nx)] += -2/(dy**2)

                LP[ii,ij2idx(i-1,j  , nx)] += 1/(dy**2)

                LP[ii,ij2idx(i+1,j  , nx)] += 1/(dy**2)

                py_ii = 0.5 * (P[i+1,j] - P[i-1,j])/(dy**2)

                b[ii] = py_ii

            elif j==nx-1: # right boundary middle cells

                # backward in x and central in y

                LP[ii,ij2idx(i  ,j  , nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j-1, nx)] += -2/(dx**2)

                LP[ii,ij2idx(i  ,j-2, nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j  , nx)] += -2/(dy**2)

                LP[ii,ij2idx(i-1,j  , nx)] += 1/(dy**2)

                LP[ii,ij2idx(i+1,j  , nx)] += 1/(dy**2)

                py_ii = 0.5 * (P[i+1,j] - P[i-1,j])/(dy**2)

                b[ii] = py_ii

            else: # middle cells

                # central in x and central in y

                LP[ii,ij2idx(i  ,j  , nx)] += -2/(dx**2)

                LP[ii,ij2idx(i  ,j-1, nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j+1, nx)] += 1/(dx**2)

                LP[ii,ij2idx(i  ,j  , nx)] += -2/(dy**2)

                LP[ii,ij2idx(i-1,j  , nx)] += 1/(dy**2)

                LP[ii,ij2idx(i+1,j  , nx)] += 1/(dy**2)

                py_ii = 0.5 * (P[i+1,j] - P[i-1,j])/(dy**2)

                b[ii] = py_ii

    return LP, b



# # Comparing with the analytical solution

def HeleShaw(x, y, pin, pout, L, D, mu):

    '''Analytical solution to the Hele-Shaw flow'''

    P = pin + ((pout-pin)*x/L)

    U = -(0.5/mu) * ( (pout-pin) * y * (D-y) / L)

    return P, U

func_matrix_vector.py

0 → 100644
+269 −0

File added.

Preview size limit exceeded, changes collapsed.

func_qc.py

0 → 100644
+248 −0

File added.

Preview size limit exceeded, changes collapsed.