Nov-20-2023, 09:30 PM
(This post was last modified: Nov-21-2023, 10:18 AM by Larz60+.
Edit Reason: fixed error tags
)
Hi im getting this error:
Full code just in case:
Error:Traceback (most recent call last):
File "C:\Users\nitro-pc\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.11_qbz5n2kfra8p0\LocalCache\local-packages\Python311\site-packages\matplotlib\cbook\__init__.py", line 304, in process
func(*args, **kwargs)
File "C:\Users\nitro-pc\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.11_qbz5n2kfra8p0\LocalCache\local-packages\Python311\site-packages\matplotlib\animation.py", line 904, in _start
self._init_draw()
File "C:\Users\nitro-pc\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.11_qbz5n2kfra8p0\LocalCache\local-packages\Python311\site-packages\matplotlib\animation.py", line 1748, in _init_draw
self._draw_frame(frame_data)
File "C:\Users\nitro-pc\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.11_qbz5n2kfra8p0\LocalCache\local-packages\Python311\site-packages\matplotlib\animation.py", line 1767, in _draw_frame
self._drawn_artists = self._func(framedata, *self._args)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "c:\Users\nitro-pc\Desktop\fluid simulatio\test 5.py", line 309, in update
velocity_x, velocity_y = project(vx_updated, vy_updated, mask) # Ensure mask is used in projection if required
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "c:\Users\nitro-pc\Desktop\fluid simulatio\test 5.py", line 257, in project
div = filter(div, (width, width)) # Ensure width is passed as a tuple
^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "c:\Users\nitro-pc\Desktop\fluid simulatio\test 5.py", line 251, in filter
diffused_f = gaussian_filter(f, width)
^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\nitro-pc\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.11_qbz5n2kfra8p0\LocalCache\local-packages\Python311\site-packages\autograd\tracer.py", line 48, in f_wrapped
return f_raw(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^
File "c:\Users\nitro-pc\Desktop\fluid simulatio\test 5.py", line 209, in gaussian_filter
return scipy.ndimage.gaussian_filter(x, sigma, mode='reflect')
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\nitro-pc\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.11_qbz5n2kfra8p0\LocalCache\local-packages\Python311\site-packages\scipy\ndimage\_filters.py", line 375, in gaussian_filter
axes = [(axes[ii], sigmas[ii], orders[ii], modes[ii], radiuses[ii])
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\nitro-pc\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.11_qbz5n2kfra8p0\LocalCache\local-packages\Python311\site-packages\scipy\ndimage\_filters.py", line 376, in <listcomp>
for ii in range(num_axes) if sigmas[ii] > 1e-15]
^^^^^^^^^^^^^^^^^^@autograd.extend.primitive
def gaussian_filter(x, sigma):
# Ensure sigma is a tuple with length equal to x's number of dimensions
if np.isscalar(sigma):
sigma = (sigma,) * x.ndim # Make sigma a tuple with the same value for all dimensions
return scipy.ndimage.gaussian_filter(x, sigma, mode='reflect')
def _gaussian_filter_vjp(ans, x, sigma):
return lambda g: autograd.numpy.sum(g * ans)
autograd.extend.defvjp(gaussian_filter, _gaussian_filter_vjp)
def filter(f, width):
if np.isscalar(width):
width = (width,) * f.ndim # Ensure width is a tuple
diffused_f = gaussian_filter(f, width)
return diffused_fAny help in shedding light to this would be greatly appreciated.Full code just in case:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.gridspec import GridSpec
import scipy.ndimage
import autograd
import autograd.numpy as anp
import scipy.ndimage
# Parameters from your fluid dynamics simulation
nx, ny = 400, 400 # Grid points - making it a square grid for a circular domain
dx, dy = 2.0 / nx, 8192 / ny # Grid spacing
dt = 0.001 # Time step
viscosity = 0.1 # Fluid viscosity
alpha = 0.01 # Thermal diffusivity
g = 9.81 # Gravity
beta = 0.0034 # Thermal expansion
# Define the center and radius of the circle
center_x, center_y = nx // 2, ny // 2
radius = min(nx, ny) // 2
# Calculate radial distance from the center for each point
y, x = np.indices((ny, nx))
r = np.sqrt((x - center_x)**2 + (y - center_y)**2)
# Create a mask for the circular domain
mask = r <= radius
def smooth_mask_transition_circular(mask, radius, softness, sigma):
"""
Apply a smooth transition to the mask for a circular domain using a radial sigmoid function
and Gaussian smoothing for antialiasing.
Args:
mask (numpy.ndarray): The initial hard-edged mask.
radius (float): The radius of the circular domain.
softness (float): A value between 0 (hard edge) and 1 (soft edge) for the sigmoid transition.
sigma (float): The standard deviation for the Gaussian kernel, for antialiasing.
Returns:
numpy.ndarray: The smoothed mask.
"""
ny, nx = mask.shape
y, x = np.ogrid[-ny//2:ny//2, -nx//2:nx//2]
r = np.sqrt(x**2 + y**2)
transition_zone = np.logical_and(r >= radius * softness, r <= radius)
mask[transition_zone] = (1 - (r[transition_zone] - radius * softness) / (radius * (1 - softness)))
mask[r > radius] = 0
# Apply Gaussian smoothing to the mask
smoothed_mask = scipy.ndimage.gaussian_filter(mask.astype(float), sigma=sigma)
# Renormalize to keep the mask values between 0 and 1
smoothed_mask = np.clip(smoothed_mask, 0, 1)
return smoothed_mask
# Now create the antialiased smoothed mask for the circular domain
sigma = 2 # The sigma value for Gaussian smoothing, adjust as needed
smoothed_mask = smooth_mask_transition_circular(mask, radius, softness=0.1, sigma=sigma)
# Visualize the smoothed mask to check the transition
plt.imshow(smoothed_mask, cmap='gray')
plt.title('Smoothed Mask')
plt.colorbar()
plt.show()
# Initialize fields for fluid dynamics
velocity_x = np.zeros((ny, nx)) * mask # Apply mask to velocity fields
velocity_y = np.zeros((ny, nx)) * mask
pressure = np.zeros((ny, nx)) * mask
temperature = np.full((ny, nx), 2000.0) * mask
# Initialize the temperature field with a radial gradient
max_temp = 8000 # Maximum temperature at the center
temperature = np.zeros((ny, nx)) # Start with an array of zeros
max_allowed_temp = 8000 # Or another value based on your simulation's physical parameters
temperature = np.minimum(temperature, max_allowed_temp)
# Define a radius for the central heat source area
heat_source_radius = radius * 0.1 # for example, 10% of the overall radius
heat_source_rate = 100 # This is the rate at which the temperature increases in the central area
def apply_heat_source(temperature, dt, heat_source_rate, center_x, center_y, heat_source_radius):
ny, nx = temperature.shape
y, x = np.ogrid[-center_y:ny-center_y, -center_x:nx-center_x]
r = np.sqrt(x**2 + y**2)
central_area = r <= heat_source_radius
# Add heat source rate to the temperature within the central area
temperature[central_area] += heat_source_rate * dt
# Ensure the temperature does not exceed the maximum allowed temperature
temperature = np.minimum(temperature, max_allowed_temp)
return temperature
# Adjust heat_source_rate accordingly
temperature = apply_heat_source(temperature, heat_source_rate, dt, heat_source_radius, center_x, center_y)
# Apply a radial gradient from the edge of the heat source to the boundary of the circular domain
gradient_mask = (r > heat_source_radius) & (r <= radius)
temperature[gradient_mask] = max_temp * (1 - (r[gradient_mask] - heat_source_radius) / (radius - heat_source_radius))
# Enforce the heat sink condition at the boundary of the circular domain
temperature[~mask] = 0
# Define initial smoke density
initial_smoke_density = 0.9
occlusion = 0
red_smoke = np.zeros((ny, nx))
blue_smoke = np.zeros((ny, nx))
red_smoke[ny//4:ny//2, :] = initial_smoke_density # Initial red smoke band
blue_smoke[ny//2:3*ny//4, :] = initial_smoke_density # Initial blue smoke band
# Additional helper function to compute the strain rate tensor
def compute_strain_rate(velocity_x, velocity_y, dx, dy):
du_dx = np.gradient(velocity_x, axis=1) / dx
du_dy = np.gradient(velocity_x, axis=0) / dy
dv_dx = np.gradient(velocity_y, axis=1) / dx
dv_dy = np.gradient(velocity_y, axis=0) / dy
S_ij = np.sqrt(2 * (du_dx**2 + dv_dy**2 + 0.5 * (du_dy + dv_dx)**2))
return S_ij
# Make sure to define the Smagorinsky constant Cs somewhere in your code
Cs = 0.1 # or whatever value is appropriate for your simulation
# Smagorinsky subgrid-scale model function
def smagorinsky_sgs(velocity_x, velocity_y, dx, dy, Cs, mask):
# Calculate the strain rate tensor
S_ij = compute_strain_rate(velocity_x, velocity_y, dx, dy)
# Calculate the subgrid scale turbulent viscosity
delta = (dx * dy)**0.5 # Filter width, assumed to be proportional to the grid size
nu_t = (Cs * delta)**2 * S_ij # Turbulent viscosity
# Add the SGS turbulent viscosity to the effective viscosity
velocity_x += nu_t * mask
velocity_y += nu_t * mask
return velocity_x, velocity_y
# Fluid dynamics functions
def step(velocity_x, velocity_y, temperature, mask, dt, dx, dy, alpha, beta, g, center_x, center_y, heat_source_radius, max_temp, Cs):
ny, nx = temperature.shape # Assuming 'temperature' is a 2D numpy array
y, x = np.indices((ny, nx))
r = np.sqrt((x - center_x)**2 + (y - center_y)**2)
# Normalize radial distance to be between 0 and 1
r_normalized = r / np.max(r)
# Calculate buoyancy force
buoyancy_force = g * beta * (temperature - np.mean(temperature[mask])) * r_normalized
# Avoid division by zero at the center by using a mask
near_center = r < 0.5 # Adjust as necessary to avoid the centermost point(s)
r[near_center] = np.inf # Assign infinity to avoid division by zero
# Calculate the unit vector components in the radial direction
radial_unit_vector_x = np.zeros_like(x, dtype=np.float64)
radial_unit_vector_y = np.zeros_like(y, dtype=np.float64)
radial_unit_vector_x = (x - center_x) / r
radial_unit_vector_y = (y - center_y) / r
# Normalize the radial unit vectors to unit length, avoiding the centermost point(s)
radial_unit_vector_x[near_center] = 0
radial_unit_vector_y[near_center] = 0
# Update velocities with buoyancy
velocity_x += dt * buoyancy_force * radial_unit_vector_x
velocity_y += dt * buoyancy_force * radial_unit_vector_y
# Apply the mask to confine updates within the circular domain
velocity_x *= smoothed_mask
velocity_y *= smoothed_mask
# Include the effect of the SGS model
velocity_x, velocity_y = smagorinsky_sgs(velocity_x, velocity_y, dx, dy, Cs, mask)
# Diffusion calculation for temperature, ensuring it only applies within the mask
for i in range(1, ny - 1):
for j in range(1, nx - 1):
if mask[i, j]: # Apply updates only within the circular domain
temperature[i, j] = temperature[i, j] + alpha * dt * (
(temperature[i+1, j] - 2 * temperature[i, j] + temperature[i-1, j]) / (dx**2) +
(temperature[i, j+1] - 2 * temperature[i, j] + temperature[i, j-1]) / (dy**2)
)
# Set boundary condition for temperature at the edge of the mask to emulate a heat sink
temperature[~mask] = 0
# Apply the maximum temperature within a small radius at the center
y, x = np.ogrid[-center_y:ny-center_y, -center_x:nx-center_x]
center_mask = x**2 + y**2 <= heat_source_radius**2
temperature[center_mask] = max_temp
return velocity_x, velocity_y, temperature
# Example usage of the step function with SGS model
Cs = 0.1 # Smagorinsky constant, which you may need to adjust based on your simulation
velocity_x, velocity_y, temperature = step(velocity_x, velocity_y, temperature, mask, dt, dx, dy, alpha, beta, g, center_x, center_y, heat_source_radius, max_temp, Cs)
@autograd.extend.primitive
def gaussian_filter(x, sigma):
# Ensure sigma is a tuple with length equal to x's number of dimensions
if np.isscalar(sigma):
sigma = (sigma,) * x.ndim # Make sigma a tuple with the same value for all dimensions
return scipy.ndimage.gaussian_filter(x, sigma, mode='reflect')
def _gaussian_filter_vjp(ans, x, sigma):
return lambda g: autograd.numpy.sum(g * ans)
autograd.extend.defvjp(gaussian_filter, _gaussian_filter_vjp)
# A set of functions for simulating the wind tunnel
def occlude(x, occlusion):
return x
def sigmoid(x):
return 0.5*(np.tanh(x) + 1.0)
def advect(f, vx, vy):
"""Instead of moving the cell centers forward in time using the velocity fields,
we look for the particles which end up exactly at the cell centers by tracing
backwards in time from the cell centers. See 'implicit Euler integration.'"""
rows, cols = f.shape
cell_xs, cell_ys = np.meshgrid(np.arange(cols), np.arange(rows))
center_xs = (cell_xs - vx).ravel() # look backwards one timestep
center_ys = (cell_ys - vy).ravel()
left_ix = np.floor(center_ys).astype(int) # get locations of source cells.
top_ix = np.floor(center_xs).astype(int)
rw = center_ys - left_ix # relative weight of cells on the right
bw = center_xs - top_ix # same for cells on the bottom
left_ix = np.mod(left_ix, rows) # wrap around edges of simulation.
right_ix = np.mod(left_ix + 1, rows)
top_ix = np.mod(top_ix, cols)
bot_ix = np.mod(top_ix + 1, cols)
# a linearly-weighted sum of the 4 cells closest to the source of the cell center.
flat_f = (1 - rw) * ((1 - bw)*f[left_ix, top_ix] + bw*f[left_ix, bot_ix]) \
+ rw * ((1 - bw)*f[right_ix, top_ix] + bw*f[right_ix, bot_ix])
return np.reshape(flat_f, (rows, cols))
def filter(f, width):
if np.isscalar(width):
width = (width,) * f.ndim # Ensure width is a tuple
diffused_f = gaussian_filter(f, width)
return diffused_f
def project(vx, vy, width=0.4):
p = np.zeros(vx.shape)
div = -0.5 * (np.roll(vx, -1, axis=1) - np.roll(vx, 1, axis=1) + np.roll(vy, -1, axis=0) - np.roll(vy, 1, axis=0))
div = filter(div, (width, width)) # Ensure width is passed as a tuple
for k in range(50):
p = (div + np.roll(p, 1, axis=1) + np.roll(p, -1, axis=1) + np.roll(p, 1, axis=0) + np.roll(p, -1, axis=0)) / 4.0
p = filter(p, (width, width))
vx = vx - 0.5 * (np.roll(p, -1, axis=1) - np.roll(p, 1, axis=1))
vy = vy - 0.5 * (np.roll(p, -1, axis=0) - np.roll(p, 1, axis=0))
return vx, vy
def calculate_vorticity(vx, vy):
dvx_dy = np.gradient(vx, axis=0)
dvy_dx = np.gradient(vy, axis=1)
vorticity = (dvy_dx - dvx_dy) * mask
return vorticity
# Set up the figure for visualization
fig = plt.figure(figsize=(15, 5))
gs = GridSpec(1, 3, figure=fig)
# Temperature plot
ax1 = fig.add_subplot(gs[0, 0])
cax1 = ax1.imshow(temperature, cmap='hot', interpolation='nearest')
fig.colorbar(cax1, ax=ax1, label='Temperature (°C)')
ax1.title.set_text('Temperature Distribution')
# Velocity magnitude plot
ax2 = fig.add_subplot(gs[0, 1])
cax2 = ax2.imshow(np.sqrt(velocity_x**2 + velocity_y**2), cmap='viridis', interpolation='nearest')
fig.colorbar(cax2, ax=ax2, label='Velocity Magnitude')
ax2.title.set_text('Velocity Magnitude')
# Setup for vorticity visualization (cax3)
ax3 = fig.add_subplot(gs[0, 2])
cax3 = ax3.imshow(calculate_vorticity(velocity_x, velocity_y), cmap='RdBu', interpolation='nearest')
fig.colorbar(cax3, ax=ax3, label='Vorticity')
ax3.title.set_text('Vorticity')
# Function to update the state in each simulation step
def update(frame):
global velocity_x, velocity_y, temperature, mask, Cs # Include Cs as a global variable
# Pass Cs as the last argument to the step function
velocity_x, velocity_y, temperature = step(
velocity_x, velocity_y, temperature, mask, dt, dx, dy, alpha, beta, g,
center_x, center_y, heat_source_radius, max_temp, Cs
)
# Calculate vorticity and apply the mask
vorticity = calculate_vorticity(velocity_x, velocity_y)
# Wind tunnel smoke advection and projection
vx_updated = advect(velocity_x, velocity_x, velocity_y)
vy_updated = advect(velocity_y, velocity_x, velocity_y)
velocity_x, velocity_y = project(vx_updated, vy_updated, mask) # Ensure mask is used in projection if required
# Update visualization for temperature and velocity magnitude
cax1.set_data(temperature)
cax2.set_data(np.sqrt(velocity_x**2 + velocity_y**2))
cax3.set_data(vorticity)
# Set up the figure for visualization
fig = plt.figure(figsize=(15, 5))
gs = GridSpec(1, 3, figure=fig)
# Temperature plot
ax1 = fig.add_subplot(gs[0, 0])
cax1 = ax1.imshow(temperature, cmap='hot', interpolation='nearest')
fig.colorbar(cax1, ax=ax1, label='Temperature (°C)')
ax1.title.set_text('Temperature Distribution')
# Velocity magnitude plot
ax2 = fig.add_subplot(gs[0, 1])
cax2 = ax2.imshow(np.sqrt(velocity_x**2 + velocity_y**2), cmap='viridis', interpolation='nearest')
fig.colorbar(cax2, ax=ax2, label='Velocity Magnitude')
ax2.title.set_text('Velocity Magnitude')
# Add a subplot for vorticity visualization
ax3 = fig.add_subplot(gs[0, 2])
cax3 = ax3.imshow(calculate_vorticity(velocity_x, velocity_y), cmap='RdBu', interpolation='nearest')
fig.colorbar(cax3, ax=ax3, label='Vorticity')
ax3.title.set_text('Vorticity')
# Create and start the animation
ani = FuncAnimation(fig, update, frames=100, interval=20)
plt.show()
