From 54bda855e5adedc4a6ba020970887f357db0ff14 Mon Sep 17 00:00:00 2001 From: zomseffen Date: Mon, 20 Dec 2021 17:21:46 +0100 Subject: [PATCH] travking fluid particle --- Client/Client.py | 92 ++++++++++--- FluidSim/FluidSimParameters.py | 2 +- FluidSim/FluidSimulator.py | 185 ++++++++++++++++++++++++++ FluidSim/LatticeBoltzmann.py | 71 +++++++--- FluidSim/StaggeredArray.py | 108 ++++++++++++++++ FluidSim/__init__.py | 0 Objects/World.py | 229 +++++++++++++++++++++++++++++++++ WorldProvider/WorldProvider.py | 1 + tests/test_FluidSimulator.py | 29 +++++ tests/test_Staggered_Array.py | 22 ++++ 10 files changed, 702 insertions(+), 37 deletions(-) create mode 100644 FluidSim/FluidSimulator.py create mode 100644 FluidSim/StaggeredArray.py create mode 100644 FluidSim/__init__.py create mode 100644 tests/test_FluidSimulator.py create mode 100644 tests/test_Staggered_Array.py diff --git a/Client/Client.py b/Client/Client.py index 9baea65..d87efc3 100644 --- a/Client/Client.py +++ b/Client/Client.py @@ -103,6 +103,53 @@ class Client: self.world_provider.world.put_object(x_pos, y_pos, z_pos, Cuboid().setColor( random.randint(0, 100) / 100.0, random.randint(0, 100) / 100.0, random.randint(0, 100) / 100.0)) + colors = {} + for plate in range(int(np.max(self.world_provider.world.plates))): + colors[plate + 1] = (random.randint(0, 100) / 100.0, + random.randint(0, 100) / 100.0, + random.randint(0, 100) / 100.0) + + for x_pos in range(0, 100): + for y_pos in range(0, 100): + for z_pos in range(0, 1): + if self.world_provider.world.plates[x_pos, y_pos] == -1: + r, g, b, = 0, 0, 1 #0.5, 0.5, 0.5 + else: + r, g, b = colors[int(self.world_provider.world.plates[x_pos, y_pos])] + self.world_provider.world.set_color(x_pos, y_pos, z_pos, r, g, b) + + total_x = self.world_provider.world.chunk_n_x * self.world_provider.world.chunk_size_x + total_y = self.world_provider.world.chunk_n_y * self.world_provider.world.chunk_size_y + for x_pos in range(0, 100): + for y_pos in range(0, 100): + if self.world_provider.world.faults[x_pos, y_pos] == -2: + self.world_provider.world.set_color(x_pos, y_pos, 0, 0, 0, 0) + + for line_index, line in enumerate(self.world_provider.world.fault_lines): + for x_pos in range(0, 100): + for y_pos in range(0, 100): + if self.world_provider.world.faults[x_pos, y_pos] == line_index: + if line_index != 9: + self.world_provider.world.set_color(x_pos, y_pos, 0, 0, 0, 1) + else: + self.world_provider.world.set_color(x_pos, y_pos, 0, 1, 1, 1) + + for x_pos in range(0, 100): + for y_pos in range(0, 100): + for z_pos in range(0, 1): + if [x_pos, y_pos] in self.world_provider.world.fault_nodes: + r, g, b = 1, 0, 0 + self.world_provider.world.set_color(x_pos, y_pos, z_pos, r, g, b) + + # # visualize direction lengths + # lengths = np.sqrt(np.sum(np.square(self.world_provider.world.directions), axis=2)) + # lengths = lengths / np.max(lengths) + # for x_pos in range(0, 100): + # for y_pos in range(0, 100): + # for z_pos in range(0, 1): + # r, g, b = lengths[x_pos, y_pos], lengths[x_pos, y_pos], lengths[x_pos, y_pos] + # self.world_provider.world.set_color(x_pos, y_pos, z_pos, r, g, b) + self.projMatrix = perspectiveMatrix(45.0, 400 / 400, 0.01, MAX_DISTANCE) self.rx = self.cx = self.cy = 0 @@ -222,28 +269,31 @@ class Client: max_value_vel = np.max(vel) # max_value_vel = np.sqrt(3) - print('round') - print('sum n: %f' % np.sum(self.n)) - print('max n: %f' % np.max(self.n)) - print('min n: %f' % np.min(self.n)) - print('sum vel: %f' % np.sum(vel)) - print('max vel: %f' % np.max(vel)) - print('min vel: %f' % np.min(vel)) + # print('round') + # print('sum n: %f' % np.sum(self.n)) + # print('max n: %f' % np.max(self.n)) + # print('min n: %f' % np.min(self.n)) + # print('sum vel: %f' % np.sum(vel)) + # print('max vel: %f' % np.max(vel)) + # print('min vel: %f' % np.min(vel)) - for x_pos in range(0, 100): - for y_pos in range(0, 100): - for z_pos in range(0, 1): - if self.state == 2: - r, g, b = value_to_color(int(self.gravity_applies[x_pos, y_pos, z_pos]), 0, 1) - if self.state == 1: - r, g, b = value_to_color(vel[x_pos, y_pos, z_pos], min_value, max_value_vel) - if self.state == 0: - r, g, b = value_to_color(self.n[x_pos, y_pos, z_pos], min_value, max_value_n) - - self.world_provider.world.set_color(x_pos, y_pos, z_pos, r, g, b) - self.world_provider.world.set_color(int(round(self.test_pixel[0])), - int(round(self.test_pixel[1])), - int(round(self.test_pixel[2])), 1.0, 1.0, 1.0) + # for x_pos in range(0, 100): + # for y_pos in range(0, 100): + # for z_pos in range(0, 1): + # # if self.state == 2: + # # r, g, b = value_to_color(int(self.gravity_applies[x_pos, y_pos, z_pos]), 0, 1) + # # if self.state == 1: + # # r, g, b = value_to_color(vel[x_pos, y_pos, z_pos], min_value, max_value_vel) + # # if self.state == 0: + # # r, g, b = value_to_color(self.n[x_pos, y_pos, z_pos], min_value, max_value_n) + # r, g, b, = 128, 128, 128 + # if [x_pos, y_pos] in self.world_provider.world.fault_nodes: + # r, g, b = 128, 0, 0 + # + # self.world_provider.world.set_color(x_pos, y_pos, z_pos, r, g, b) + # self.world_provider.world.set_color(int(round(self.test_pixel[0])), + # int(round(self.test_pixel[1])), + # int(round(self.test_pixel[2])), 1.0, 1.0, 1.0) # print(1.0 / (time.time() - self.time)) self.time = time.time() diff --git a/FluidSim/FluidSimParameters.py b/FluidSim/FluidSimParameters.py index 2170be4..0680ea5 100644 --- a/FluidSim/FluidSimParameters.py +++ b/FluidSim/FluidSimParameters.py @@ -1,7 +1,7 @@ class FluidSimParameter: viscosity = 0.1 / 3.0 # Pr = 1.0 - Pr = 100.0 + Pr = 1.0 # vc = 1.0 vc = 0.5 diff --git a/FluidSim/FluidSimulator.py b/FluidSim/FluidSimulator.py new file mode 100644 index 0000000..ab520e0 --- /dev/null +++ b/FluidSim/FluidSimulator.py @@ -0,0 +1,185 @@ +from FluidSim.StaggeredArray import StaggeredArray2D +import numpy as np +import scipy +import scipy.sparse +import scipy.sparse.linalg + +class FluidSimulator2D: + def __init__(self, x_n: int, y_n: int): + self.x_n = x_n + self.y_n = y_n + + self.array = StaggeredArray2D(self.x_n, self.y_n) + + self.coordinate_array = np.zeros((x_n, y_n, 2), dtype=np.int) + for x in range(x_n): + for y in range(y_n): + self.coordinate_array[x, y, :] = x, y + + def advect(self, field: np.ndarray, delta_t: float): + u_x, u_y = self.array.get_velocity_arrays() + u = np.stack([u_x, u_y], axis=2) + + def runge_kutta_layer(input, time_elapsed, border_handling='clamp'): + shifted_pos = np.round(self.coordinate_array - u * time_elapsed).astype(np.int) * (u < 0) +\ + (self.coordinate_array - u * time_elapsed).astype(np.int) * (u > 0) + \ + self.coordinate_array * (u == 0) + # border handling + if border_handling == 'clamp': + shifted_pos = np.maximum(0, shifted_pos) + shifted_pos[:, :, 0] = np.minimum(len(input), shifted_pos[:, :, 0]) + shifted_pos[:, :, 1] = np.minimum(len(input[0]), shifted_pos[:, :, 1]) + pass + layer = np.zeros(field.shape, dtype=field.dtype) + for x in range(self.x_n): + for y in range(self.y_n): + layer[x, y] = field[shifted_pos[x, y][0], shifted_pos[x, y][1]] - field[x, y] + return layer + + k1 = runge_kutta_layer(field, delta_t) + + k2 = runge_kutta_layer(field + 0.5 * delta_t * k1, 0.5 * delta_t) + + k3 = runge_kutta_layer(field + 0.75 * delta_t * k2, 0.75 * delta_t) # maybe 0.25 instead? + + # new_field = field + 2.0 / 9.0 * delta_t * k1 + 3.0 / 9.0 * delta_t * k2 + 4.0 / 9.0 * delta_t * k3 + new_field = field + k1 + + return new_field + + def get_timestep(self, f, h=1.0, k_cfl=1.0): + f_length = np.max(np.sqrt(np.sum(np.square(f), axis=-1))) + + u_x, u_y = self.array.get_velocity_arrays() + u_length = np.max(np.sqrt(np.square(u_x) + np.square(u_y))) + + # return k_cfl * h / (u_length + np.sqrt(h + f_length)) + return k_cfl * h / u_length + + def update_velocity(self, timestep, border_handling='constant'): + if border_handling == 'constant': + p_diff_x = np.pad(self.array.p, [(0, 1), (0, 0)], mode='constant', constant_values=0) -\ + np.pad(self.array.p, [(1, 0), (0, 0)], mode='constant', constant_values=0) + borders_fluid_x = np.pad(self.array.has_fluid, [(0, 1), (0, 0)], mode='constant', constant_values=False) +\ + np.pad(self.array.has_fluid, [(1, 0), (0, 0)], mode='constant', constant_values=False) + + p_diff_y = np.pad(self.array.p, [(0, 0), (0, 1)], mode='constant', constant_values=0) -\ + np.pad(self.array.p, [(0, 0), (1, 0)], mode='constant', constant_values=0) + borders_fluid_y = np.pad(self.array.has_fluid, [(0, 0), (0, 1)], mode='constant', constant_values=False) +\ + np.pad(self.array.has_fluid, [(0, 0), (1, 0)], mode='constant', constant_values=False) + else: + p_diff_x = 0 + p_diff_y = 0 + borders_fluid_x = False + borders_fluid_y = False + + u_x_new = self.array.u_x - (timestep * p_diff_x) * (1.0 * borders_fluid_x) + u_y_new = self.array.u_y - (timestep * p_diff_y) * (1.0 * borders_fluid_y) + + # clear all components that do not border the fluid + u_x_new *= (1.0 * borders_fluid_x) + u_y_new *= (1.0 * borders_fluid_y) + + return u_x_new, u_y_new + + def calculate_divergence(self, h=1.0): + dx_u = (self.array.u_x[1:, :] - self.array.u_x[:-1, :]) / h + dy_u = (self.array.u_y[:, 1:] - self.array.u_y[:, -1]) / h + + return dx_u + dy_u + + def pressure_solve(self, divergence): + new_p = np.zeros((self.array.x_n, self.array.y_n)) + connection_matrix = np.zeros((self.array.x_n * self.array.y_n, self.array.x_n * self.array.y_n)) + flat_divergence = np.zeros((self.array.x_n * self.array.y_n)) + + for x in range(self.array.x_n): + for y in range(self.array.y_n): + flat_divergence[x * self.array.y_n + y] = divergence[x, y] + + neighbors = 4 + if x == 0: + neighbors -= 1 + else: + connection_matrix[x * self.array.y_n + y, (x - 1) * self.array.y_n + y] = 1 + if y == 0: + neighbors -= 1 + else: + connection_matrix[x * self.array.y_n + y, x * self.array.y_n + y - 1] = 1 + if x == self.array.x_n - 1: + neighbors -= 1 + else: + connection_matrix[x * self.array.y_n + y, (x + 1) * self.array.y_n + y] = 1 + if y == self.array.y_n - 1: + neighbors -= 1 + else: + connection_matrix[x * self.array.y_n + y, x * self.array.y_n + y + 1] = 1 + + connection_matrix[x * self.array.y_n + y, x * self.array.y_n + y] = -neighbors + + p = scipy.sparse.linalg.spsolve(connection_matrix, -flat_divergence) + + for x in range(self.array.x_n): + for y in range(self.array.y_n): + new_p[x, y] = p[x * self.array.y_n + y] + return new_p + + def timestep(self, external_f, h=1.0, k_cfl=1.0): + """ + :param external_f: external forces to be applied + :param h: grid size + :param k_cfl: speed up multiplier (reduces accuracy if > 1.0. Does not make much sense if smaller than 1.0) + :return: + """ + delta_t = self.get_timestep(external_f, h, k_cfl) + + has_fluid = self.advect(self.array.has_fluid * 1.0, delta_t) + self.array.density = self.advect(self.array.density, delta_t) + + # temp_u_x = self.advect(self.array.u_x, delta_t) + # temp_u_y = self.advect(self.array.u_y, delta_t) + # self.array.u_x = temp_u_x + # self.array.u_y = temp_u_y + #TODO advect velocity + test_u = np.stack(self.array.get_velocity_arrays(), axis=2) + test = self.advect(np.stack(self.array.get_velocity_arrays(), axis=2), delta_t) + + # TODO maybe use dynamic threshold to conserve amount of cells containing fluid + self.array.has_fluid = has_fluid >= 0.5 + # add more stuff to advect here. For example temperature, density, and other things. Maybe advect velocity. + + # self.array.u_x, self.array.u_y = self.update_velocity(delta_t) + # TODO add forces (a = F / m) -> add a to u + + dx_u = (self.array.u_x[1:, :] - self.array.u_x[:-1, :]) / h + dy_u = (self.array.u_y[:, 1:] - self.array.u_y[:, -1]) / h + + dx_u = self.advect(dx_u, delta_t) + dy_u = self.advect(dy_u, delta_t) + + divergence = dx_u + dy_u + # divergence = self.calculate_divergence(h) + + self.array.p = self.pressure_solve(divergence) + + self.array.u_x, self.array.u_y = self.update_velocity(delta_t) + + +if __name__ == '__main__': + fs = FluidSimulator2D(50, 50) + + fs.array.has_fluid[10, 10] = True + + for i in range(100): + fs.timestep(0) + print(i) + + print(fs.array.has_fluid[10, 10]) + print(np.sum(fs.array.has_fluid * 1.0)) + # test = fs.advect(fs.array.has_fluid * 1.0, 1.0) + # + # test2 = fs.update_velocity(1.0) + # + # test3 = fs.calculate_divergence() + # + # test4 = fs.pressure_solve(test3) diff --git a/FluidSim/LatticeBoltzmann.py b/FluidSim/LatticeBoltzmann.py index b39a8b3..07d5acb 100644 --- a/FluidSim/LatticeBoltzmann.py +++ b/FluidSim/LatticeBoltzmann.py @@ -37,6 +37,10 @@ def main(): # Initial Conditions N = np.ones((Ny, Nx, NL)) # * rho0 / NL temperature = np.ones((Ny, Nx, NL), np.float) # * rho0 / NL + + tracked_fluid = np.zeros((Ny, Nx, NL), np.float) # * rho0 / NL + tracked_fluid[50:, :, 0] = 1 + has_fluid = np.ones((Ny, Nx), dtype=np.bool) has_fluid[int(Ny/2):, :] = False np.random.seed(42) @@ -56,13 +60,17 @@ def main(): # Cylinder boundary X, Y = np.meshgrid(range(Nx), range(Ny)) + cylinder = (X - Nx / 4) ** 2 + (Y - Ny / 2) ** 2 < (Ny / 4) ** 2 - inner_cylinder = (X - Nx / 4) ** 2 + (Y - Ny / 2) ** 2 < (Ny / 4 - 2) ** 2 + cylinder[:, :] = False + N[cylinder] = 0 N[0, :] = 0 N[Ny - 1, :] = 0 temperature[cylinder] = 0 + + tracked_fluid[cylinder] = 0 # N[int(Ny/2):, :] = 0 has_fluid[cylinder] = False @@ -98,14 +106,17 @@ def main(): temperature[:, :, i] = np.roll(temperature[:, :, i], cx, axis=1) temperature[:, :, i] = np.roll(temperature[:, :, i], cy, axis=0) + tracked_fluid[:, :, i] = np.roll(tracked_fluid[:, :, i], cx, axis=1) + tracked_fluid[:, :, i] = np.roll(tracked_fluid[:, :, i], cy, axis=0) + # has_fluid = new_has_fluid > 0.5 # new_has_fluid[F_sum == 0] += has_fluid[F_sum == 0] * 1.0 # new_has_fluid[(np.abs(F_sum) < 0.000000001)] = 0 - - fluid_sum = np.sum(has_fluid * 1.0) - has_fluid = (new_has_fluid / np.sum(new_has_fluid * 1.0)) * fluid_sum - - print('fluid_cells: %d' % np.sum(has_fluid * 1)) + # + # fluid_sum = np.sum(has_fluid * 1.0) + # has_fluid = (new_has_fluid / np.sum(new_has_fluid * 1.0)) * fluid_sum + # + # print('fluid_cells: %d' % np.sum(has_fluid * 1)) # for i in idxs: # N[:, :, i] *= has_fluid @@ -126,9 +137,13 @@ def main(): bndryTemp = temperature[bndry, :] bndryTemp = bndryTemp[:, reflection_mapping] + bndryTracked = tracked_fluid[bndry, :] + bndryTracked = bndryTracked[:, reflection_mapping] + sum_f = np.sum(N) print('Sum of Particles: %f' % sum_f) print('Sum of Temperature: %f' % np.sum(temperature)) + print('Sum of tacked particles: %f' % np.sum(tracked_fluid)) # sum_f_cyl = np.sum(N[cylinder]) # print('Sum of Forces in cylinder: %f' % sum_f_cyl) @@ -147,6 +162,9 @@ def main(): # Calculate fluid variables rho = np.sum(N, 2) temp_rho = np.sum(temperature, 2) + + tracked_rho = np.sum(tracked_fluid, 2) + ux = np.sum(N * cxs, 2) / rho uy = np.sum(N * cys, 2) / rho @@ -154,16 +172,28 @@ def main(): uy[(np.abs(rho) < epsilon)] = 0 g = -params.g * (temp_rho - yy / Ny) + + uy1 = (uy + g * params.t1 * (tracked_rho * 0.9 + 0.1)) + uy2 = (uy + g * params.t2 * (tracked_rho * 0.9 + 0.1)) + # uy[np.abs(rho) >= epsilon] += g[np.abs(rho) >= epsilon] / 2.0 - uy += g / 2.0 + # uy += g / 2.0 # u_length = np.maximum(np.abs(ux), np.abs(uy)) - u_length = np.sqrt(np.square(ux) + np.square(uy)) + u_length1 = np.sqrt(np.square(ux) + np.square(uy1)) + u_length2 = np.sqrt(np.square(ux) + np.square(uy2)) + + u_max_length1 = np.max(u_length1) + u_max_length2 = np.max(u_length2) + u_max_length = max(u_max_length1, u_max_length2) + + if u_max_length > 2: + test = 1 - u_max_length = np.max(u_length) if u_max_length > np.sqrt(2): ux = (ux / u_max_length) * np.sqrt(2) - uy = (uy / u_max_length) * np.sqrt(2) + uy1 = (uy1 / u_max_length) * np.sqrt(2) + uy2 = (uy2 / u_max_length) * np.sqrt(2) print('max vector part: %f' % u_max_length) # ux /= u_max_length @@ -193,13 +223,17 @@ def main(): # Apply Collision temperature_eq = np.zeros(temperature.shape) + tracked_eq = np.zeros(temperature.shape) Neq = np.zeros(N.shape) for i, cx, cy, w in zip(idxs, cxs, cys, weights): Neq[:, :, i] = rho * w * ( - 1 + 3 * (cx * ux + cy * uy) + 9 * (cx * ux + cy * uy) ** 2 / 2 - 3 * (ux ** 2 + uy ** 2) / 2) + 1 + 3 * (cx * ux + cy * uy1) + 9 * (cx * ux + cy * uy1) ** 2 / 2 - 3 * (ux ** 2 + uy1 ** 2) / 2) temperature_eq[:, :, i] = temp_rho * w * ( - 1 + 3 * (cx * ux + cy * uy) + 9 * (cx * ux + cy * uy) ** 2 / 2 - 3 * (ux ** 2 + uy ** 2) / 2) + 1 + 3 * (cx * ux + cy * uy2) + 9 * (cx * ux + cy * uy2) ** 2 / 2 - 3 * (ux ** 2 + uy2 ** 2) / 2) + + tracked_eq[:, :, i] = tracked_rho * w * ( + 1 + 3 * (cx * ux + cy * uy1) + 9 * (cx * ux + cy * uy1) ** 2 / 2 - 3 * (ux ** 2 + uy1 ** 2) / 2) # test1 = np.sum(Neq) test2 = np.sum(N-Neq) if abs(test2) > 0.0001: @@ -212,10 +246,12 @@ def main(): N += -(1.0 / params.t1) * (N - Neq) temperature += -(1.0 / params.t2) * (temperature - temperature_eq) + tracked_fluid += -(1.0 / params.t1) * (tracked_fluid - tracked_eq) # Apply boundary N[bndry, :] = bndryN temperature[bndry, :] = bndryTemp + tracked_fluid[bndry, :] = bndryTracked # temperature[0, :, 0] = 0 # temperature[1, :, 0] = 0 @@ -223,8 +259,9 @@ def main(): temperature[0, :, 0] /= 2 temperature[1, :, 0] /= 2 - temperature[Ny - 1, :, 0] = 1 - temperature[Ny - 2, :, 0] = 1 + if it <= 3000: + temperature[Ny - 1, :, 0] = 1 + temperature[Ny - 2, :, 0] = 1 # n_sum = np.sum(N, 2) # n_sum_min = np.min(n_sum) @@ -280,8 +317,12 @@ def main(): plt.subplot(2, 2, 3) max_N = np.max(np.sum(N, 2)) plt.imshow(np.sum(N, 2) / max_N * 2.0 - 1.0, cmap='bwr') - plt.clim(-.1, .1) + + plt.subplot(2, 2, 4) + plt.imshow(np.sum(tracked_fluid, 2) * 2.0 - 1.0, cmap='bwr') + plt.clim(-.1, .1) + # ax = plt.gca() # ax.invert_yaxis() # ax.get_xaxis().set_visible(False) diff --git a/FluidSim/StaggeredArray.py b/FluidSim/StaggeredArray.py new file mode 100644 index 0000000..ca17a93 --- /dev/null +++ b/FluidSim/StaggeredArray.py @@ -0,0 +1,108 @@ +import numpy as np +from scipy.signal import convolve + + +class StaggeredArray2D: + def __init__(self, x_n, y_n): + """ + creates a staggered array + :param x_n: x size of the array + :param y_n: y size of the array + """ + self.x_n = x_n + self.y_n = y_n + self.p = np.zeros((x_n, y_n), dtype=np.float) + + self.u_x = np.zeros((x_n + 1, y_n), dtype=np.float) + self.u_y = np.zeros((x_n, y_n + 1), dtype=np.float) + + self.has_fluid = np.zeros((x_n, y_n), dtype=np.bool) + self.density = np.zeros((x_n, y_n), dtype=np.float) + + def get_velocity(self, x, y): + """ + get mid point value for the coordinates + :param x: x coordinate + :param y: y coordinate + :return: + """ + assert 0 <= x < self.x_n, 'x value out of bounds!' + assert 0 <= y < self.y_n, 'y value out of bounds!' + + lower_x = self.u_x[x, y] + upper_x = self.u_x[x + 1, y] + + lower_y = self.u_y[x, y] + upper_y = self.u_y[x, y + 1] + + return (lower_x + upper_x) / 2.0, (lower_y + upper_y) / 2.0 + + def get_velocity_arrays(self): + c_x = np.array([[0.5], [0.5]]) + u_x = convolve(self.u_x, c_x, mode='valid') + + c_y = np.array([[0.5, 0.5]]) + u_y = convolve(self.u_y, c_y, mode='valid') + + return u_x, u_y + + +class StaggeredArray3D: + def __init__(self, x_n, y_n, z_n): + """ + creates a staggered array + :param x_n: x size of the array + :param y_n: y size of the array + :param z_n: z size of the array + """ + + self.x_n = x_n + self.y_n = y_n + self.z_n = z_n + + self.p = np.zeros((x_n, y_n, z_n), dtype=np.float) + + self.u_x = np.zeros((x_n + 1, y_n, z_n), dtype=np.float) + self.u_y = np.zeros((x_n, y_n + 1, z_n), dtype=np.float) + self.u_z = np.zeros((x_n, y_n, z_n + 1), dtype=np.float) + + self.has_fluid = np.zeros((x_n, y_n, z_n), dtype=np.bool) + + def get_velocity(self, x, y, z): + """ + get mid point value for the coordinates + :param x: x coordinate + :param y: y coordinate + :param z: z coordinate + :return: + """ + assert 0 <= x < self.x_n, 'x value out of bounds!' + assert 0 <= y < self.y_n, 'y value out of bounds!' + assert 0 <= z < self.z_n, 'y value out of bounds!' + + lower_x = self.u_x[x, y, z] + upper_x = self.u_x[x + 1, y, z] + + lower_y = self.u_y[x, y, z] + upper_y = self.u_y[x, y + 1, z] + + lower_z = self.u_z[x, y, z] + upper_z = self.u_z[x, y, z + 1] + + return (lower_x + upper_x) / 2.0, (lower_y + upper_y) / 2.0, (lower_z + upper_z) / 2.0 + + +if __name__ == '__main__': + sa = StaggeredArray2D(10, 10) + + for x in range(11): + for y in range(10): + sa.u_x[x, y] = y + + for x in range(10): + for y in range(11): + sa.u_y[x, y] = x + + ux, uy = sa.get_velocity_arrays() + + ux2, uy2 = sa.get_velocity(0, 0) \ No newline at end of file diff --git a/FluidSim/__init__.py b/FluidSim/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/Objects/World.py b/Objects/World.py index 4443f01..e21bdfe 100644 --- a/Objects/World.py +++ b/Objects/World.py @@ -7,6 +7,8 @@ from OpenGL.GLU import * from OpenGL.GL import * import math import numpy as np +import random +import sys class WorldChunk(Structure): def __init__(self, width: int, length: int, height: int, programs: dict): @@ -200,6 +202,11 @@ class World(Renderable): self.chunk_n_z = chunk_n_z self.programs = programs + self.fault_nodes = [] + self.fault_lines = [] + self.plates = None + self.directions = None + self.chunks: [[[WorldChunk]]] = [] for x in range(chunk_n_x): self.chunks.append([]) @@ -208,6 +215,228 @@ class World(Renderable): for z in range(chunk_n_z): self.chunks[x][y].append(None) + def generate(self, seed=None, sea_height=50, continental_height=200): + if seed is None: + seed = random.randrange(2**32) + seed = 229805811 + print('Generation seed is %i' % seed) + random.seed(seed) + np.random.seed(seed) + node_n = self.chunk_n_x + self.chunk_n_y + total_x = self.chunk_n_x * self.chunk_size_x + total_y = self.chunk_n_y * self.chunk_size_y + nodes = [] + for _ in range(node_n): + nodes.append([random.randint(0, total_x - 1), random.randint(0, total_y - 1)]) + + # connections = np.random.randint(2, 5, len(nodes)) #np.zeros(len(nodes)) + 3 + connections = (np.abs(np.random.normal(0, 5, len(nodes))) + 2).astype(np.int) + + def calc_min_vector(start, end): + dx = end[0] - start[0] + wrapped_dx = dx % total_x + + dy = end[1] - start[1] + wrapped_dy = dy % total_y + + vector = np.array([dx, dy]) + if wrapped_dx < abs(dx): + vector[0] = wrapped_dx + if wrapped_dy < abs(dy): + vector[1] = wrapped_dy + return vector + + def is_intersecting_any(start, end, edges): + vec1 = calc_min_vector(start, end) + + for (start2_index, end2_index) in edges: + start2 = nodes[start2_index] + end2 = nodes[end2_index] + + vec2 = calc_min_vector(start2, end2) + + norm1 = vec1 / np.sqrt(np.sum(np.square(vec1))) + norm2 = vec2 / np.sqrt(np.sum(np.square(vec2))) + + # parrallel + parallel_threshold = 0.0001 + if np.sqrt(np.sum(np.square(norm1 - norm2))) < parallel_threshold or np.sqrt(np.sum(np.square(norm1 + norm2))) < parallel_threshold: + t = (start[0] - start2[0]) / vec2[0] + t2 = (end[0] - start2[0]) / vec2[0] + s = (start2[0] - start[0]) / vec1[0] + s2 = (end2[0] - start[0]) / vec1[0] + if (start2[1] + t * vec2[1]) - start[1] < parallel_threshold: + if (0 <= t <= 1.0 and 0 <= t2 <= 1.0) or (0 <= s <= 1.0 and 0 <= s2 <= 1.0): + return True + else: + if start != start2 and end != end2 and end != start2 and start != end2: + t = (vec1[0] * start[1] + vec1[1] * start2[0] - vec1[1] * start[0] - start2[1] * vec1[0]) / (vec2[1] * vec1[0] - vec2[0] * vec1[1]) + if 0 <= t <= 1.0: + intersection = np.array(start2) + vec2 * t + s = (intersection[0] - start[0]) / vec1[0] + if 0 <= s <= 1.0: + return True + + return False + + + + for index, node in enumerate(nodes): + distances = [] + for other_index, other_node in enumerate(nodes): + if node != other_node and (index, other_index) not in self.fault_lines and\ + (other_index, index) not in self.fault_lines: + if (not is_intersecting_any(node, other_node, self.fault_lines)) and (not is_intersecting_any(other_node, node, self.fault_lines)): + distances.append((other_index, np.sqrt(np.sum(np.square(calc_min_vector(node, other_node)))))) + + distances.sort(key=lambda element: element[1]) + while connections[index] > 0 and len(distances) > 0: + self.fault_lines.append((index, distances[0][0])) + connections[distances[0][0]] -= 1 + connections[index] -= 1 + distances.pop(0) + + self.fault_nodes = nodes + + plates = np.zeros((total_x, total_y)) + faults = np.zeros((total_x, total_y)) - 1 + plate_bordering_fault = {} + # draw fault lines + for fault_index, fault_line in enumerate(self.fault_lines): + start = self.fault_nodes[fault_line[0]] + end = self.fault_nodes[fault_line[1]] + vector = calc_min_vector(start, end) + vector = vector / np.sqrt(np.sum(np.square(vector))) + + point = np.array(start, dtype=np.float) + plate_bordering_fault[fault_index] = [] + + while np.sqrt(np.sum(np.square(point - np.array(end)))) > 0.5: + plates[int(point[0]), int(point[1])] = -1 + if faults[int(point[0]), int(point[1])] == -1: + faults[int(point[0]), int(point[1])] = fault_index + elif faults[int(point[0]), int(point[1])] != fault_index: + faults[int(point[0]), int(point[1])] = -2 + point += 0.5 * vector + point[0] %= total_x + point[1] %= total_y + self.faults = faults + + plate = 1 + while np.any(plates == 0): + start = np.where(plates == 0) + start = (start[0][0], start[1][0]) + plates[start] = plate + work_list = [start] + + while len(work_list) > 0: + work = work_list.pop() + + up = (work[0], (work[1] + 1) % total_y) + down = (work[0], (work[1] - 1) % total_y) + left = ((work[0] - 1) % total_x, work[1]) + right = ((work[0] + 1) % total_x, work[1]) + + if plates[up] == -1 and plates[down] == -1 and plates[left] == -1 and plates[right] == -1: + plates[work] = -1 + continue + + if plates[up] <= 0: + if plates[up] == 0: + work_list.append(up) + plates[up] = plate + if plates[down] <= 0: + if plates[down] == 0: + work_list.append(down) + plates[down] = plate + if plates[left] <= 0: + if plates[left] == 0: + work_list.append(left) + plates[left] = plate + if plates[right] <= 0: + if plates[right] == 0: + work_list.append(right) + plates[right] = plate + + if faults[up] > -1: + if plate not in plate_bordering_fault[faults[up]]: + plate_bordering_fault[faults[up]].append(plate) + if faults[down] > -1: + if plate not in plate_bordering_fault[faults[down]]: + plate_bordering_fault[faults[down]].append(plate) + if faults[left] > -1: + if plate not in plate_bordering_fault[faults[left]]: + plate_bordering_fault[faults[left]].append(plate) + if faults[right] > -1: + if plate not in plate_bordering_fault[faults[right]]: + plate_bordering_fault[faults[right]].append(plate) + plate += 1 + + plate_num = plate + for plate in range(1, plate_num): + if np.sum(plates == plate) < 20: + plates[plates == plate] = -1 + for key, item in plate_bordering_fault.items(): + if plate in item: + item.remove(plate) + + directions = np.zeros((total_x, total_y, 3)) + heights = np.zeros((total_x, total_y)) + for plate in range(1, plate_num): + if random.randint(1, 2) == 1: + heights[plates == plate] = sea_height + else: + heights[plates == plate] = continental_height + + coords = np.zeros((total_x, total_y, 2)) + for x in range(total_x): + for y in range(total_y): + coords[x, y, 0] = x + coords[x, y, 1] = y + + for fault_index, fault_line in enumerate(self.fault_lines): + start = self.fault_nodes[fault_line[0]] + end = self.fault_nodes[fault_line[1]] + vector = calc_min_vector(start, end) + vector = vector / np.sqrt(np.sum(np.square(vector))) + + perpendicular = np.array([vector[1], -vector[0]]) + + if len(plate_bordering_fault[fault_index]) == 2: + for plate in plate_bordering_fault[fault_index]: + vecs = coords - np.array(start) + lengths = np.sqrt(np.sum(np.square(vecs), axis=2, keepdims=True)) + norm_vecs = vecs / lengths + scalars = np.sum(norm_vecs * perpendicular, axis=2, keepdims=True) + scalars[lengths == 0] = 0 + + end_vecs = coords - np.array(end) + end_lengths = np.sqrt(np.sum(np.square(end_vecs), axis=2, keepdims=True)) + end_min_length = np.min(end_lengths[np.logical_and(plates == plate, end_lengths[:, :, 0] > 0)]) + end_min_length_scalar = scalars[np.logical_and(plates == plate, end_lengths[:, :, 0] == end_min_length)][0, 0] + + min_length = np.min(lengths[np.logical_and(plates == plate, lengths[:, :, 0] > 0)]) + min_length_scalar = scalars[np.logical_and(plates == plate, lengths[:, :, 0] == min_length)][0, 0] + + mean_scalar = np.mean(scalars[plates == plate]) + + if (min_length_scalar / abs(min_length_scalar)) == (end_min_length_scalar / abs(end_min_length_scalar)): + scalar = min_length_scalar + else: + if (min_length_scalar / abs(min_length_scalar)) == (mean_scalar / abs(mean_scalar)): + scalar = min_length_scalar + else: + scalar = end_min_length_scalar + + directions[plates == plate, :2] += perpendicular * (scalar / abs(scalar)) + pass + + + + + self.plates = plates + self.directions = directions + def set_color(self, x: int, y: int, z: int, r: float, g: float, b: float): x = x % (self.chunk_size_x * self.chunk_n_x) y = y % (self.chunk_size_y * self.chunk_n_y) diff --git a/WorldProvider/WorldProvider.py b/WorldProvider/WorldProvider.py index 0401eca..a9629b7 100644 --- a/WorldProvider/WorldProvider.py +++ b/WorldProvider/WorldProvider.py @@ -4,6 +4,7 @@ from Objects.World import World class WorldProvider: def __init__(self, programs): self.world: World = World(10, 10, 10, 10, 10, 10, programs) + self.world.generate() def update(self): pass diff --git a/tests/test_FluidSimulator.py b/tests/test_FluidSimulator.py new file mode 100644 index 0000000..bc6b31b --- /dev/null +++ b/tests/test_FluidSimulator.py @@ -0,0 +1,29 @@ +from FluidSim.FluidSimulator import FluidSimulator2D +import numpy as np + + +def test_stand_still(): + fs = FluidSimulator2D(50, 50) + + fs.array.has_fluid[10, 10] = True + + for i in range(100): + fs.timestep(0) + + assert fs.array.has_fluid[10, 10], "Fluid not on the same spot anymore" + assert np.sum(fs.array.has_fluid * 1.0) == 1.0, "Fluid amount changed" + + +def test_move_positive_x(): + fs = FluidSimulator2D(50, 50) + + fs.array.has_fluid[10, 10] = True + fs.array.u_x[10, 10] = 1.0 + fs.array.u_x[11, 10] = 1.0 + # fs.array.u_x[9, 10] = -1.0 + + for i in range(10): + fs.timestep(0) + assert np.sum(fs.array.has_fluid * 1.0) == 1.0, "Fluid amount changed" + + assert fs.array.has_fluid[0, 10], "Fluid not on the right border" diff --git a/tests/test_Staggered_Array.py b/tests/test_Staggered_Array.py new file mode 100644 index 0000000..ea61960 --- /dev/null +++ b/tests/test_Staggered_Array.py @@ -0,0 +1,22 @@ +from FluidSim.StaggeredArray import StaggeredArray2D + + +def test_staggered_array_2D(): + sa = StaggeredArray2D(10, 10) + + for x in range(11): + for y in range(10): + sa.u_x[x, y] = y + + for x in range(10): + for y in range(11): + sa.u_y[x, y] = x + + ux, uy = sa.get_velocity_arrays() + + for x in range(10): + for y in range(10): + ux2, uy2 = sa.get_velocity(x, y) + + assert ux[x, y] == ux2, 'value output should be consistent!' + assert uy[x, y] == uy2, 'value output should be consistent!'