diff --git a/.gitignore b/.gitignore index d64fca2..1d1838c 100644 --- a/.gitignore +++ b/.gitignore @@ -8,4 +8,5 @@ output/ build/ .vscode/ venv/ -*.so \ No newline at end of file +venv_dry/ +*.mp4 diff --git a/CMakeLists.txt b/CMakeLists.txt index d33101a..72b2172 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,6 +6,11 @@ find_package(pybind11 REQUIRED) # Create platform-agnostic core library add_library(core STATIC HGD/motion/core.cpp HGD/motion/stress.cpp) + +# === 【關鍵修改】加入這一行來解決 -fPIC 錯誤 === +set_target_properties(core PROPERTIES POSITION_INDEPENDENT_CODE ON) +# ============================================ + target_include_directories(core PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/HGD/motion) target_compile_features(core PUBLIC cxx_std_14) @@ -20,4 +25,4 @@ target_compile_options(d2q4_cpp PRIVATE -O3 -march=native) # Ensure the module is placed inside HGD/motion/ set_target_properties(d2q4_cpp PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${CMAKE_SOURCE_DIR}/HGD/motion -) +) \ No newline at end of file diff --git a/HGD/motion/d2q4_array_v3.py b/HGD/motion/d2q4_array_v3.py new file mode 100644 index 0000000..2692d89 --- /dev/null +++ b/HGD/motion/d2q4_array_v3.py @@ -0,0 +1,449 @@ +import sys +import numpy as np +from numpy.typing import ArrayLike +from HGD import operators +from HGD import stress + + +# @njit + +import numpy as np + + +def particle_potential_velocity( + u: np.ndarray, + v: np.ndarray, # [nx, ny, nm] + s: np.ndarray, # [nx, ny, nm] + alpha_p_cell2d: np.ndarray, # [nx, ny] - Solid volume fraction + p, # Parameter object + *, + dt: float | None = None, + gx: float = 0.0, + gy: float | None = None, +): + """ + Simplified Inertial Velocity Calculation: + 1. Assumes Static Gas (Ug = 0, Vg = 0). + 2. Uses a Blend Model (Wen-Yu + Ergun) for drag force calculation. + 3. Integrates gravity and drag to calculate potential velocity after time dt. + """ + + # --- 1. Basic Parameter Setup --- + if dt is None: + dt = p.dt + if gy is None: + gy = -abs(getattr(p, "g", 9.81)) + + # Physical constants + rho_g = p.gas_density # Gas density + mu_g = p.gas_viscosity # Gas viscosity + rho_p = p.solid_density # Particle density + eps = 1e-30 # Epsilon to prevent division by zero + + # Expand dimensions to match [nx, ny, nm] for broadcasting + # alpha_p_cell2d is originally [nx, ny], needs to be [nx, ny, 1] + alpha_p3 = alpha_p_cell2d[:, :, np.newaxis] + + # Determine occupancy (where particles exist) + occ = np.isfinite(s) + + # Calculate fluid volume fraction (alpha_f) + # alpha_f = 1 - alpha_s + alpha_f3 = np.clip(1.0 - alpha_p3, 1e-6, 1.0) + + # Safe particle diameter (handle NaNs by replacing with 1.0, though masked later) + d_safe = np.where(occ, np.maximum(s, eps), 1.0) + + # --- 2. Calculate Relative Velocity & Reynolds Number --- + Ur = np.sqrt(u * u + v * v) + + # Reynolds Number (Re_p) + Re = rho_g * Ur * d_safe / (mu_g + eps) + Re = np.maximum(Re, 1e-12) + + # Standard Drag Coefficient Cd (Schiller-Naumann correlation) + Cd = 24.0 / Re * (1.0 + 0.15 * Re**0.687) + + # --- 3. Calculate Beta (Drag Coefficient Term) --- + + # A. Wen-Yu Model (Suitable for sparse flow) + beta_wy = 0.75 * Cd * rho_g * (Ur / (d_safe + eps)) * np.power(alpha_f3, -2.65) + + # B. Ergun Model (Suitable for dense flow / packed bed) + beta_er = 150.0 * mu_g * (alpha_p3**2) / (alpha_f3 * (d_safe**2) + eps) + 1.75 * rho_g * alpha_p3 * ( + Ur / (d_safe + eps) + ) + mask_wy = alpha_f3 >= 0.8 + beta = np.where(mask_wy, beta_wy, beta_er) + + # Ensure drag is only applied where particles actually exist + beta = np.where(occ, beta, 0.0) + + # --- 4. Time Integration --- + # We solve the equation: dv/dt = g_eff - a*v + # where 'a' is the drag relaxation coefficient: a = beta / rho_p + + # Calculate relaxation coefficient 'a' + a = beta / (rho_p + eps) + + # Clip 'a' to prevent numerical instability + a = np.clip(a, 1e-12, None) + + # Calculate exponential decay factor + expf = np.exp(-a * dt) + + # Buoyancy-corrected gravity + buoy = (rho_p - rho_g) / (rho_p + eps) + bx, by = gx * buoy, gy * buoy + + # Analytic Solution for Linear Drag over time step dt: + # v_new = v_terminal + (v_old - v_terminal) * exp(-a*dt) + # where v_terminal = force / a + u_pot = (bx / a) + (u - bx / a) * expf + v_pot = (by / a) + (v - by / a) * expf + + return np.nan_to_num(u_pot), np.nan_to_num(v_pot) + + +def move_voids( + u: ArrayLike, + v: ArrayLike, + s: ArrayLike, + p, + diag: int = 0, + c: None | ArrayLike = None, + T: None | ArrayLike = None, + chi: None | ArrayLike = None, + last_swap: None | ArrayLike = None, +) -> tuple[ArrayLike, ArrayLike, ArrayLike, None | ArrayLike, None | ArrayLike]: + """ + Function to move voids each timestep. + + Args: + u: Storage container for counting how many voids moved horizontally + v: Storage container for counting how many voids moved vertically + s: 3D array containing the local sizes everywhere. `NaN`s represent voids. Other values represent the grain size. The first two dimensions represent real space, the third dimension represents the micro-structural coordinate. + diag: Should the voids swap horizontally (von neumnann neighbourhood, `diag=0`) or diagonally upwards (moore neighbourhood, `diag=1`). Default value `0`. + c: If ArrayLike, a storage container for tracking motion of differently labelled particles. If `None`, do nothing. + T: If ArrayLike, the temperature field. If `None`, do nothing. + boundary: If ArrayLike, a descriptor of cells which voids cannot move into (i.e. boundaries). If `internal_boundary` is defined in the params file, allow for reduced movement rates rather than zero. If `None`, do nothing. + + Returns: + u: The updated horizontal velocity + v: The updated vertical velocity + s: The new locations of the grains + c: The updated concentration field + T: The updated temperature field + """ + + options = [(1, -1)] # , (0, -1), (0, 1)] # up, left, right + for d in range(1, p.max_diff_swap_length + 1): + options.extend([(0, -d), (0, d)]) + + P_diff_weighting = ( + p.max_diff_swap_length * (p.max_diff_swap_length + 1) * (2 * p.max_diff_swap_length + 1) / 6 + ) + # np.random.shuffle(options) # oh boy, this is a massive hack + N_swap = np.zeros([p.nx, p.ny], dtype=int) + + nu = operators.get_solid_fraction(s) + + solid = nu >= p.nu_cs + Solid = np.repeat(solid[:, :, np.newaxis], p.nm, axis=2) + + unstable = np.isnan(s) * ~Solid # & ~swapped # & ~Skip + + s_bar = operators.get_average(s) + S_bar = np.repeat(s_bar[:, :, np.newaxis], p.nm, axis=2) + + if p.slope_stability_model == "stress": + potential_free_surface = None + else: + potential_free_surface = operators.empty_up(nu, p) + + swap_indices = [] + dest_indices = [] + + u_new = np.zeros_like(u) + v_new = np.zeros_like(v) + if p.inertia is not False: + u, v = particle_potential_velocity(u, v, s, nu, p, dt=p.dt, gx=0, gy=-9.81) + + for axis, d in options: + s_dest = np.roll(s, d, axis=axis) + S_bar_dest = np.roll(S_bar, d, axis=axis) + + if p.advection_model == "average_size": + u_here = np.sqrt(p.g * s_bar) + U_dest = np.sqrt(p.g * S_bar_dest) + elif p.advection_model == "freefall": + u_here = U_dest = np.sqrt(p.g * p.dy) + elif p.advection_model == "stress": + sigma = stress.calculate_stress(s, last_swap, p) + pressure = np.abs( + stress.get_pressure(sigma, p) + ) # HACK: PRESSURE SHOULD BE POSITIVE BUT I HAVE ISSUES WITH THE STRESS MODEL + u_here = np.sqrt(2 * pressure / p.solid_density) + U = np.repeat(u_here[:, :, np.newaxis], p.nm, axis=2) + U_dest = np.roll( + U, d, axis=axis + ) # NEED TO TAKE DESTINATION VALUE BECAUSE PRESSURE IS ZERO AT OUTLET!!! + + with np.errstate(divide="ignore", invalid="ignore"): + beta = np.exp(-p.P_stab * p.dt / (p.dx / u_here)) + # print(beta) + + if axis == 1: # vertical + s_inv_bar = operators.get_hyperbolic_average(s) + S_inv_bar = np.repeat(s_inv_bar[:, :, np.newaxis], p.nm, axis=2) + S_inv_bar_dest = np.roll(S_inv_bar, d, axis=axis) + + P = (p.dt / p.dy) * U_dest * (S_inv_bar_dest / s_dest) ** 0.1 + + if p.inertia is not False: + v_new_swapper = np.roll(v, d, axis=axis) + P = (np.abs(v_new_swapper)) * p.dt / p.dy + if d == -1: + P[v_new_swapper >= 0] = 0 # up + else: + P[v_new_swapper <= 0] = 0 # down + P[:, -1, :] = 0 # no swapping up from top row + + elif axis == 0: # horizontal + P = p.alpha * U_dest * s_dest * (p.dt / p.dy / p.dy) * P_diff_weighting + + if p.inertia is not False: + if d < 0: # left + U_valid = np.where(u < 0, u, 0) + else: # right + U_valid = np.where(u > 0, u, 0) + P += (p.dt / p.dy) * np.roll(U_valid, d, axis=axis) + + if d > 0: # left + P[:d, :, :] = 0 # no swapping left from leftmost column + else: # right + P[d:, :, :] = 0 # no swapping right from rightmost column + + if p.slope_stability_model == "gradient": + slope_stable = operators.stable_slope_fast(s, d, p, chi, potential_free_surface) + elif p.slope_stability_model == "stress": + slope_stable = operators.stable_slope_stress(s, p, last_swap) + else: + sys.exit(f"Invalid slope stability model: {p.slope_stability_model}") + + # Prevent swaps OUT FROM stable slope cells + unstable = np.logical_or(unstable, ~slope_stable) + # P[slope_stable] = 0 + + filled_dest = ~np.isnan(s_dest) + swap_possible = np.logical_and(unstable, filled_dest) + + if p.inertia: + if axis == 0: + u_new += np.where(swap_possible, p.alpha * U_dest * s_dest / p.dy * d, 0) + u_new[0, :, :] = 0 # no swapping left from leftmost column + u_new[-1, :, :] = 0 # no swapping right from rightmost column + elif axis == 1: + v_new += np.where(swap_possible, U_dest, 0) + v_new[:, -1, :] = 0 # no swapping up from top row + + P = np.where(swap_possible, P, 0) + swap = np.random.rand(*P.shape) < P + + this_swap_indices = np.argwhere(swap) + this_dest_indices = this_swap_indices.copy() + this_dest_indices[:, axis] -= d + + swap_indices.extend(this_swap_indices) + dest_indices.extend(this_dest_indices) + + # Prevent conflicts by filtering out swaps that would cause two voids to swap into the same cell + swap_indices_conflict_free, dest_indices_conflict_free = prevent_conflicts(swap_indices, dest_indices) + + if len(swap_indices_conflict_free) > 0: + # Prevent overfilling by limiting swaps into cells based on maximum allowed nu + swap_indices_filtered, dest_indices_filtered = prevent_overfilling( + swap_indices_conflict_free, dest_indices_conflict_free, nu, potential_free_surface, p + ) + + ( + s[swap_indices_filtered[:, 0], swap_indices_filtered[:, 1], swap_indices_filtered[:, 2]], + s[dest_indices_filtered[:, 0], dest_indices_filtered[:, 1], dest_indices_filtered[:, 2]], + ) = ( + s[dest_indices_filtered[:, 0], dest_indices_filtered[:, 1], dest_indices_filtered[:, 2]], + s[swap_indices_filtered[:, 0], swap_indices_filtered[:, 1], swap_indices_filtered[:, 2]], + ) + + last_swap[swap_indices_filtered[:, 0], swap_indices_filtered[:, 1], swap_indices_filtered[:, 2]] = ( + 2 * axis - 1 + ) # 1 for up, -1 for left or right + + # N_swap[swap_indices_filtered[:, 0], swap_indices_filtered[:, 1]] += 1 + np.add.at(N_swap, (swap_indices_filtered[:, 0], swap_indices_filtered[:, 1]), 1) + + delta = dest_indices_filtered - swap_indices_filtered + + if p.inertia is False: + # Update the new velocities where the mass is + u_new[swap_indices_filtered[:, 0], swap_indices_filtered[:, 1], swap_indices_filtered[:, 2]] = ( + -delta[:, 0] * p.dx / p.dt + ) + v_new[swap_indices_filtered[:, 0], swap_indices_filtered[:, 1], swap_indices_filtered[:, 2]] = ( + -delta[:, 1] * p.dy / p.dt + ) + else: + ( + v[swap_indices_filtered[:, 0], swap_indices_filtered[:, 1], swap_indices_filtered[:, 2]], + v[dest_indices_filtered[:, 0], dest_indices_filtered[:, 1], dest_indices_filtered[:, 2]], + ) = ( + v[dest_indices_filtered[:, 0], dest_indices_filtered[:, 1], dest_indices_filtered[:, 2]], + v[swap_indices_filtered[:, 0], swap_indices_filtered[:, 1], swap_indices_filtered[:, 2]], + ) + + ( + u[swap_indices_filtered[:, 0], swap_indices_filtered[:, 1], swap_indices_filtered[:, 2]], + u[dest_indices_filtered[:, 0], dest_indices_filtered[:, 1], dest_indices_filtered[:, 2]], + ) = ( + u[dest_indices_filtered[:, 0], dest_indices_filtered[:, 1], dest_indices_filtered[:, 2]], + u[swap_indices_filtered[:, 0], swap_indices_filtered[:, 1], swap_indices_filtered[:, 2]], + ) + # Zero out the velocities for the swapped voids + u[dest_indices_filtered[:, 0], dest_indices_filtered[:, 1], dest_indices_filtered[:, 2]] = 0 + v[dest_indices_filtered[:, 0], dest_indices_filtered[:, 1], dest_indices_filtered[:, 2]] = 0 + + if p.inertia is False: + u = 0 + v = 0 + else: + u[Solid] = u[Solid] + v[Solid] = v[Solid] + last_swap[np.isnan(s)] = np.nan + + # chi_new = ( + chi = ( + N_swap / p.nm / p.P_stab + ) # NOTE: Should this be P_stab??? Should be one when everything swaps as fast as possible... + # chi = beta * chi + (1 - beta) * chi_new + + return u, v, s, c, T, chi, last_swap + + +def prevent_conflicts(swap_indices, dest_indices): + """ + Identifies and removes conflicts in swap and destination indices. + This function takes two lists or arrays of indices, `swap_indices` and + `dest_indices`, and identifies any conflicts where the same index appears + more than once in either list. It then filters out any swaps involved in + these conflicts, returning only the conflict-free swaps. + Parameters: + swap_indices (list or np.ndarray): Indices representing the source positions + for swaps. + dest_indices (list or np.ndarray): Indices representing the destination + positions for swaps. + Returns: + tuple: A tuple containing two numpy arrays: + - conflict-free swap indices + - conflict-free destination indices + """ + + swap_indices = np.array(swap_indices) + dest_indices = np.array(dest_indices) + + # Identify duplicates in swap_indices (reading conflicts) + _, swap_inv, swap_counts = np.unique(swap_indices, axis=0, return_inverse=True, return_counts=True) + swap_conflict_mask = swap_counts[swap_inv] > 1 + + # Identify duplicates in dest_indices (writing conflicts) + _, dest_inv, dest_counts = np.unique(dest_indices, axis=0, return_inverse=True, return_counts=True) + dest_conflict_mask = dest_counts[dest_inv] > 1 + + # Combine masks to filter out any swaps involved in conflicts + conflict_mask = swap_conflict_mask | dest_conflict_mask + conflict_free_mask = ~conflict_mask + + if len(conflict_free_mask) > 0: + # Extract conflict-free swaps + swap_indices = swap_indices[conflict_free_mask] + dest_indices = dest_indices[conflict_free_mask] + + return swap_indices, dest_indices + + +def prevent_overfilling(swap_indices, dest_indices, nu, potential_free_surface, p): + """ + Prevents overfilling of destination locations by limiting the number of swaps. + + Parameters: + swap_indices (ndarray): Array of indices representing the source locations for swaps. + dest_indices (ndarray): Array of indices representing the destination locations for swaps. + nu (ndarray): Array representing the current state of the system. + potential_free_surface (ndarray or None): Array indicating potential free surface points, or None if not applicable. + p (object): Parameter object containing the following attributes: + - nm (int): Scaling factor for the number of swaps. + - nu_cs (float): Critical value for nu. + - delta_limit (float or ndarray): Limit for the change in nu during swaps. + + Returns: + tuple: A tuple containing the filtered swap_indices and dest_indices after applying the overfilling prevention logic. + """ + # Randomly permute the swaps + perm = np.random.permutation(len(swap_indices)) + swap_indices = swap_indices[perm] + dest_indices = dest_indices[perm] + + # Compute the change in nu for each swap + nu_source = nu[swap_indices[:, 0], swap_indices[:, 1]] + nu_dest = nu[dest_indices[:, 0], dest_indices[:, 1]] + + # swap_sideways = swap_indices[:, 0] != dest_indices[:, 0] + swap_vertical = swap_indices[:, 1] != dest_indices[:, 1] + delta_nu = nu_dest - nu_source + + # Allow for vertical swaps to exceed the delta limit + delta_nu[swap_vertical] = 1 + + # Count the number of swaps into each source location + unique_locs, inverse_indices = np.unique(swap_indices[:, :2], axis=0, return_inverse=True) + # counts = np.bincount(inverse_indices) + + # For each unique location, compute the allowed number of swaps + if np.isscalar(p.nu_cs): + max_swaps_bulk = p.nm * (p.nu_cs - nu[unique_locs[:, 0], unique_locs[:, 1]]) + else: + max_swaps_bulk = p.nm * ( + p.nu_cs[unique_locs[:, 0], unique_locs[:, 1]] - nu[unique_locs[:, 0], unique_locs[:, 1]] + ) + + if potential_free_surface is None: + max_swaps = max_swaps_bulk[inverse_indices] + else: + if np.isscalar(p.delta_limit): + delta_limit = p.delta_limit + else: + delta_limit = p.delta_limit[unique_locs[:, 0], unique_locs[:, 1]][inverse_indices] + max_swaps_slope = ((delta_nu - delta_limit) * p.nm).astype(int) + max_swaps_slope = np.maximum(max_swaps_slope, 0) + + # Check for potential free surface points + unique_potential = potential_free_surface[unique_locs[:, 0], unique_locs[:, 1]] + + max_swaps = np.where( + unique_potential[inverse_indices], max_swaps_slope, max_swaps_bulk[inverse_indices] + ) + + max_swaps = np.minimum(max_swaps, max_swaps_bulk[inverse_indices]) + + # Find the first occurrence of each unique source location + _, first_indices = np.unique(inverse_indices, return_index=True) + + # Compute position in group + position_in_group = np.arange(len(swap_indices)) - first_indices[inverse_indices] + + # Build keep mask + keep_mask = position_in_group < max_swaps + + # Apply the mask to swap_indices and dest_indices + swap_indices = swap_indices[keep_mask] + dest_indices = dest_indices[keep_mask] + + return swap_indices, dest_indices diff --git a/HGD/params.py b/HGD/params.py index 218aec1..86abbd5 100644 --- a/HGD/params.py +++ b/HGD/params.py @@ -184,28 +184,34 @@ def update_before_time_march(self, cycles): P_ratio = self.alpha * self.s_M * P_n / self.dy print(f"n = {n}, n*P_diff/P_adv = {n*P_ratio}") - while not safe: - self.P_adv_ref = stability - self.dt = self.P_adv_ref * self.dy / self.free_fall_velocity - - # self.P_diff_ref = self.alpha * self.P_adv_ref - self.P_diff_max = ( - self.alpha - * self.s_M - * self.free_fall_velocity - * self.dt - / self.dy - / self.dy - * self.P_diff_weighting - ) - - self.P_adv_max = self.P_adv_ref * (self.s_M / self.s_m) - # self.P_diff_max = self.P_diff_ref * (self.s_M / self.s_m) - - if self.P_adv_max + (2 * self.max_diff_swap_length) * self.P_diff_max <= self.P_stab: - safe = True - else: - stability *= 0.95 + if self.inertia: + # === Inertia Mode: Use fixed time step === + self.dt = self.defined_time_step_size + else: + # === Non-Inertia Mode: Original automatic stability loop === + safe = False # Initialize loop condition + while not safe: + self.P_adv_ref = stability + self.dt = self.P_adv_ref * self.dy / self.free_fall_velocity + + # self.P_diff_ref = self.alpha * self.P_adv_ref + self.P_diff_max = ( + self.alpha + * self.s_M + * self.free_fall_velocity + * self.dt + / self.dy + / self.dy + * self.P_diff_weighting + ) + + self.P_adv_max = self.P_adv_ref * (self.s_M / self.s_m) + # self.P_diff_max = self.P_diff_ref * (self.s_M / self.s_m) + + if self.P_adv_max + (2 * self.max_diff_swap_length) * self.P_diff_max <= self.P_stab: + safe = True + else: + stability *= 0.95 if self.t_f is None: self.nt = 0 diff --git a/HGD/plotter.py b/HGD/plotter.py index 5f64c65..0ea16ef 100644 --- a/HGD/plotter.py +++ b/HGD/plotter.py @@ -6,7 +6,7 @@ import matplotlib.pyplot as plt from matplotlib.colors import LinearSegmentedColormap import matplotlib.cm as cm - +import matplotlib.colors as mcolors from PIL import Image from HGD import operators from HGD import stress @@ -271,6 +271,8 @@ def update(p, state, *args): plot_s(s, p) if "s_blue_pink" in p.plot: plot_s_blue_pink(s, p) + if "combined_plot" in p.plot: + plot_combined(s, p, p.tstep, p.t_f / p.nt, v, u) if "nu" in p.plot: plot_nu(s, p) if "rel_nu" in p.plot: @@ -1172,3 +1174,134 @@ def stack_videos(paths, name, p): # buffer = generate_colormap_image(256, 256) # with open('colormap.png', 'wb') as f: # f.write(buffer.read()) + + +def plot_combined(s, p, t, dt, v, u): + """ + Plots the Solid Fraction (left) and Particle Size (right) side-by-side. + Boundary masks have been removed. + """ + + # Create a figure with 1 row and 2 columns + fig, axs = plt.subplots(1, 2, figsize=(14, 6)) + + # Set figure and axes backgrounds to white + fig.patch.set_facecolor("white") + axs[0].set_facecolor("white") + axs[1].set_facecolor("white") + + # Scale factors for display units + scale_factor = 1_000_000 # Convert meters to micrometers (µm) + scale_factor_nu = 1 # Solid fraction is dimensionless + + # ---- Left Plot: Solid Fraction (nu) ---- + + # Calculate solid fraction: 1 minus the average void fraction + nu_raw = (1 - np.mean(np.isnan(s), axis=2).T) * scale_factor_nu # shape (ny, nx) + + # Mask cells where solid fraction is exactly 0 so the background remains white + nu = np.ma.masked_where(nu_raw == 0.0, nu_raw) + + # Define grid coordinates + x = np.linspace(0, p.W, p.nx) + y = np.linspace(0, p.H, p.ny) + + # Plot solid fraction heatmap + c1 = axs[0].pcolormesh( + x, + y, + nu, + cmap="inferno_r", + vmin=0, + vmax=scale_factor_nu, + shading="auto", + ) + + # Set labels and aspect ratio for the left plot + axs[0].set_xlabel("X (m)") + axs[0].set_ylabel("Y (m)") + axs[0].set_aspect("equal", adjustable="box") + + # ---- Right Plot: Particle Size (s_plot) ---- + + # Calculate average particle size, ignoring NaNs (voids) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) + s_plot = np.nanmean(s, axis=2).T * scale_factor + + # Prepare data for plotting (mask NaNs so they show as background) + s_plot = np.ma.masked_where(np.isnan(s_plot), s_plot) + + # Create a custom colormap (Pink to Blue) + # pink_blue_cmap = mcolors.LinearSegmentedColormap.from_list( + # "pink_blue", ["#ff69b4", "#0000ff"], N=256 + # ) + + # # Select colormap based on parameters + # if getattr(p, "orange_blue_cmap", False): + # # Assuming orange_blue_cmap is defined elsewhere globally or in imports + # # If not, you might need to define it or pass it in. + # # For safety, falling back to pink_blue if variable not found is safer, + # # but adhering to original logic here. + # try: + # cmap_to_use = orange_blue_cmap + # except NameError: + # cmap_to_use = pink_blue_cmap + # else: + # cmap_to_use = pink_blue_cmap + + # Plot particle size heatmap + c2 = axs[1].pcolormesh( + x, + y, + s_plot, + cmap=orange_blue_cmap, + vmin=np.nanmin(s) * 1e6, + vmax=np.nanmax(s) * 1e6, + shading="auto", + ) + + # Set labels and aspect ratio for the right plot + axs[1].set_xlabel("X (m)") + axs[1].set_ylabel("Y (m)") + axs[1].set_aspect("equal", adjustable="box") + + # ---- Titles and Colorbars ---- + + # Set global title with simulation time + time = t * dt + fig.suptitle(f"Time = {time:.2f} s", fontsize=12, y=1.05) + + # Add colorbar for Solid Fraction (Top Left) + cbar_ax1 = fig.add_axes([0.1, 0.92, 0.35, 0.02]) + cbar1 = fig.colorbar(c1, cax=cbar_ax1, orientation="horizontal") + cbar1.ax.xaxis.set_ticks_position("top") + cbar1.ax.xaxis.set_label_position("top") + cbar1.set_label("Solid fraction (-)", labelpad=10) + + # Add colorbar for Particle Size (Top Right) + cbar_ax2 = fig.add_axes([0.55, 0.92, 0.35, 0.02]) + cbar2 = fig.colorbar(c2, cax=cbar_ax2, orientation="horizontal") + cbar2.ax.xaxis.set_ticks_position("top") + cbar2.ax.xaxis.set_label_position("top") + cbar2.set_label("Particle size (µm)", labelpad=10) + + # Adjust layout spacing + plt.subplots_adjust(left=0.1, right=0.9, bottom=0.1, top=0.88, wspace=0.4) + + # ---- Save Output ---- + + # Define output path + output_folder = os.path.join(p.folderName, "combined") + os.makedirs(output_folder, exist_ok=True) + + # Save figure + plt.savefig( + os.path.join(output_folder, "combined" + str(t).zfill(8) + ".png"), + bbox_inches="tight", + dpi=100, + ) + + # Close plots to free memory + plt.close(fig) + plt.close() diff --git a/json/collapse.json5 b/json/collapse.json5 index ef6c2a1..df5bc55 100644 --- a/json/collapse.json5 +++ b/json/collapse.json5 @@ -28,4 +28,6 @@ // saving saves : 20, t_f : 5, // s + motion_model : 'd2q4_array_v3', // 'd2q4_slow', 'd2q4_array' or 'd2q9_slow' + } diff --git a/json/collapse_bi.json5 b/json/collapse_bi.json5 index 1f50e56..ac63ceb 100644 --- a/json/collapse_bi.json5 +++ b/json/collapse_bi.json5 @@ -13,10 +13,12 @@ repose_angle : 30, // (degrees) gsd_mode : 'bi', large_concentration: 0.5, - s_m : [0.0001,0.0003,0.0005], // minimum particle size (m) + s_m : [0.0001], // minimum particle size (m) s_M : 0.001, // maximum particle size (m) nu_cs : 0.5, alpha : 1, + inertia: [true,false], + // initial condition: fill_ratio : 0.25, @@ -27,6 +29,14 @@ save_inc : 100, t_f : null, // s videos : ['nu','s'], - save : ['nu','s'], - mask_s_bar : true + // save : ['nu','s'], + plot : ['combined_plot'], + + mask_s_bar : true, + motion_model : 'd2q4_array_v3', // 'd2q4_slow', 'd2q4_array' or 'd2q9_slow' + t_f:100, + + + defined_time_step_size:1e-3, + } diff --git a/json/defaults.json5 b/json/defaults.json5 index 1578cb7..99c1f15 100644 --- a/json/defaults.json5 +++ b/json/defaults.json5 @@ -53,12 +53,16 @@ stress_mode : 'active', // 'isotropic', 'active', 'passive', 'K_0', 'anisotropic' or 'harr' solid_density : 2700, // density of the solid particles (kg/m^3/m) wall_friction_angle : 30, // friction coefficient between the walls and the particles + + // gas + gas_density:1.225, + gas_viscosity:1.8e-5, // saving folderName : "output/", // name of the folder to save the data in save_inc : 1, // how often to save the state of the system (timesteps) t_f : 5, // final time to end the simulation (s) - plot : ['nu', 's'], // which variables to plot. Should at least include the ones to make videos. + plot : ['nu', 's', 'combined_plot'], // which variables to plot. Should at least include the ones to make videos. videos : ['nu','s'], // which variables to save videos of. save : [], // which variables to save. By default none. mask_s_bar : false, // should the images of s_bar be masked to hide areas of low density