This Jupyter Notebook explores voxel-based simulations for physical systems. It discusses optimizing voxel tessellation algorithms, including the considerations for using cubic or truncated octahedral voxel shapes. It also entails practical implementation of these concepts in Python, utilizing numpy and matplotlib for operations like defining 3D kernels for neighbor interactions in voxel grids, computing temperature differences for a modeled sphere, and visualizing these voxel grids. At its core, the notebook serves as a tutorial or example code for conducting physics simulations or studying material properties at the voxel level, including aspects like energy conservation, pressure, and temperature variations within the system.
- Should I consider changing the tesselation from cubic to truncated Octahedral - This would make some things easier (like all neighbors would share area with all other neighbors). However, some things would be harder, for example the area of each face is not constant so I have to keep track then of what kind of face makes contact. Further, I think the seperation between grid points becomes a bit harder to work with)
import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import Normalize from matplotlib.cm import ScalarMappable from scipy.ndimage import convolve, correlate from mplEasyAnimate import animation from matplotlib.cm import viridis from tqdm.notebook import tqdm plt.style.use(pubStyle)
def make_3d_kernels(): # Central kernel ck = np.zeros((3, 3, 3)) ck[1, 1, 1] = 1 # Initialize kernels array to hold all 26 kernels kernels = np.zeros((26, 3, 3, 3)) # Counter for filling the kernels array counter = 0 # Generate face neighbor kernels for axis in range(3): # Loop through x, y, z axes for delta in [-1, 1]: # Negative and positive directions k = np.copy(ck) k[1 + delta*(axis == 0), 1 + delta*(axis == 1), 1 + delta*(axis == 2)] = -1 kernels[counter, :, :, :] = k counter += 1 # Generate edge neighbor kernels for x in [-1, 1]: for y in [-1, 1]: k = np.copy(ck) k[1 + x, 1 + y, 0] = -1/2 kernels[counter, :, :, :] = k counter += 1 k = np.copy(ck) k[1 + x, 0, 1 + y] = -1/2 kernels[counter, :, :, :] = k counter += 1 k = np.copy(ck) k[0, 1 + x, 1 + y] = -1/2 kernels[counter, :, :, :] = k counter += 1 # Generate corner neighbor kernels for x in [-1, 1]: for y in [-1, 1]: for z in [-1, 1]: k = np.copy(ck) k[1 + x, 1 + y, 1 + z] = -1/3 kernels[counter, :, :, :] = k counter += 1 return kernels
def pretty_print_3d_array(array3D): # Determine the maximum number width for formatting maxNumWidth = max(len(str(abs(x))) for x in array3D.flatten() if x != 0) # Iterate through each layer for layerIndex in range(array3D.shape[0]): for rowIndex in range(array3D.shape[1]): for colIndex in range(array3D.shape[2]): element = array3D[layerIndex, rowIndex, colIndex] # Using format specifiers to align the output print(f"{element:{maxNumWidth}f}" if element != 0 else " " * maxNumWidth, end=" ") # Move to the next row print() # Add an empty line between layers for better readability print()
class EnergyConservationError(Exception): def __init__(self, msg): self.message = msg
class TimeEvolver: def __init__(self, sphere, times): self.sphere = sphere self._tracker = np.zeros_like(shape=(*self.sphere.temperature.shape, 4, len(times))) self.times = times self.dts = np.gradient(self.times) def evolve(self): for tID, dt in enumerate(self.dts): self.sphere.timestep(dt)
class CallbackNDArray(np.ndarray): def __new__(cls, input_array, callback=None, *args, **kwargs): obj = np.asarray(input_array).view(cls) obj.callback = callback return obj def __setitem__(self, key, value): super().__setitem__(key, value) if hasattr(self, 'callback') and self.callback: self.callback()
class VoxelSphere: mH = 1.00784 mHe = 4.002602 R = 8.3145e7 R = 1 _3DDiffKernels = make_3d_kernels() def __init__(self, radius, resolution, t0=0, X=0.75, Y=0.25, Z=0, D=1, k=1): self._radius = radius self._resolution = resolution self._t = t0 self._X, self._Y, self._Z = X, Y, Z self._effectiveMolarMass = self._X * self.mH + self._Y * self.mHe self._create_voxel_sphere() self._D = D self._k = k @property def radius(self): return self._radius @property def resolution(self): return self._resolution def _create_voxel_sphere(self): gridSize = np.ceil(self.radius * self.resolution) * 2 + 1 self.x = np.linspace(-self.radius, self.radius, int(gridSize)) self.y = np.linspace(-self.radius, self.radius, int(gridSize)) self.z = np.linspace(-self.radius, self.radius, int(gridSize)) xGrid, yGrid, zGrid = np.meshgrid(self.x, self.y, self.z, indexing='ij') self._dr = np.linalg.norm(np.array([xGrid[0,0,0] - xGrid[0,1,0], yGrid[0,0,0] - yGrid[1,0,0], zGrid[0,0,0] - zGrid[0,0,1]])) self._voxelVolume = np.abs((self.x[1] - self.x[0]) * (self.y[1] - self.y[0]) * (self.z[1] - self.z[0])) self._voxelFaceArea = np.abs((self.x[1] - self.x[0]) * (self.y[1] - self.y[0])) self._temperature = CallbackNDArray(np.ones_like(xGrid), callback=self._recompute_energy) self._density = CallbackNDArray(np.ones_like(xGrid), callback=self._recompute_energy) self._pressure = CallbackNDArray(np.ones_like(xGrid), callback=self._recompute_energy) distanceFromCenter = np.sqrt(xGrid**2 + yGrid**2 + zGrid**2) self._voxelSphere = distanceFromCenter <= self.radius self._recompute_energy() def _recompute_energy(self): self._n = (self._pressure * self._voxelVolume)/self._effectiveMolarMass self._energy = (3/2) * self._n * self.R * self._temperature self._energy[~self._voxelSphere] = 0 def _recompute_temperature(self): self._n = (self._pressure * self._voxelVolume)/self._effectiveMolarMass self._temperature = (2/3) * self._energy/(self._n*self.R) self._temperature[~self._voxelSphere] = 0 @property def temperature(self): return self._temperature def _temperature_differences(self): diffs = np.zeros(shape=(26, *self._temperature.shape)) for kernelID, kernel in enumerate(self._3DDiffKernels): correlation = correlate(self._temperature, kernel, mode='constant', cval=0.0) diffs[kernelID] = correlation return diffs def timestep(self, dt, ect = 0.01): totalEnergyInit = self._energy.sum() dTs = self._temperature_differences() Q = np.zeros_like(dTs) for orientationID, dT in enumerate(dTs): orientationQ = -self._k*self._voxelFaceArea*dT/(self._dr) self._energy = self._energy + dE pretty_print_3d_array(self._energy) totalEnergyFinal = self._energy.sum() fchange = abs((totalEnergyFinal - totalEnergyInit)/totalEnergyInit) if fchange > ect: raise EnergyConservationError(f"Total Energy varied more than allowed minimum fraction ({ect*100:0.3f}%) {fchange*100:0.3f}%") self._recompute_temperature() def plot_voxel_region(self, voxelSphere, voxelColors='red', **kwargs): fig, ax = plt.subplots(1, 1, figsize=(10, 7), subplot_kw=dict(projection='3d')) filled = np.zeros(voxelSphere.shape, dtype=bool) filled[np.where(voxelSphere)] = True ax.voxels(filled, facecolors=voxelColors, edgecolor='k', **kwargs) ax.set_xlabel('X Axis') ax.set_ylabel('Y Axis') ax.set_zlabel('Z Axis') return fig, ax def _mk_voxel_color_array(self, array): Norm = Normalize(vmin=array.min(), vmax=array.max()) colors = viridis(Norm(array).flatten()).reshape(*array.shape, 4) return colors def _pcparam_to_color(self, pcparam): if pcparam == None: colors = np.zeros(shape=(*self._voxelSphere.shape, 4)) colors[:,:,:] = [1,1,1,1] elif pcparam == 'temperature': colors = self._mk_voxel_color_array(self._temperature) elif pcparam == 'pressure': colors = self._mk_voxel_color_array(self._pressure) elif pcparam == 'energy': colors = self._mk_voxel_color_array(self._energy) return colors def plot_voxel_sphere(self, pcparam=None, cmap='viridis', **kwargs): colors = self._pcparam_to_color(pcparam) return self.plot_voxel_region(self._voxelSphere, voxelColors=colors, **kwargs)
sphere = VoxelSphere(1, 1) sphere.timestep(1)
0.853956
0.853956
0.853956 0.853956 0.853956
0.853956
0.853956
================================
4.269781
4.269781
4.269781 5.123737 4.269781
4.269781
4.269781
================================
5.123737
5.123737
5.123737 5.977693 5.123737
5.123737
5.123737
---------------------------------------------------------------------------
EnergyConservationError Traceback (most recent call last)
Cell In[236], line 2
1 sphere = VoxelSphere(1, 1)
----> 2 sphere.timestep(1)
Cell In[235], line 74, in VoxelSphere.timestep(self, dt, ect)
72 fchange = abs((totalEnergyFinal - totalEnergyInit)/totalEnergyInit)
73 if fchange > ect:
---> 74 raise EnergyConservationError(f"Total Energy varied more than allowed minimum fraction ({ect*100:0.3f}%) {fchange*100:0.3f}%")
75 self._recompute_temperature()
EnergyConservationError: Total Energy varied more than allowed minimum fraction (1.000%) 514.286%
a = np.zeros((10, 10)) a[5, :] = 5 a[6, :] = 6 k = np.zeros((3,3)) k[1,1] = 1 k[0,1] = -1 a
array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[5., 5., 5., 5., 5., 5., 5., 5., 5., 5.],
[6., 6., 6., 6., 6., 6., 6., 6., 6., 6.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])
array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[5., 5., 5., 5., 5., 5., 5., 5., 5., 5.],
[6., 6., 6., 6., 6., 6., 6., 6., 6., 6.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])
pkb = np.zeros((3,3)) pkb[1,1] = 1 pkb[0,1] = -1 dkb = np.zeros((3,3)) dkb[1,1] = 1 dkb[0,0] = -1/2 diffs = np.zeros(shape=(*a.shape, 8)) for rot in range(4): pk = np.rot90(pkb, k=rot) dk = np.rot90(dkb, k=rot) pc = correlate(a, pk, mode='constant', cval=0.0) dc = correlate(a, dk, mode='constant', cval=0.0) diffs[:, :, rot] = pc diffs[:, :, rot+4] = dc
def make_2d_kernels(): # for the peripheral elements pkb = np.zeros((3,3)) pkb[1,1] = 1 pkb[0,1] = -1 # for the diagonal elements dkb = np.zeros((3,3)) dkb[1,1] = 1 dkb[0,0] = -1/2 kernels = np.zeros(shape=(8,3,3)) for rot in range(4): pk = np.rot90(pkb, k=rot) dk = np.rot90(dkb, k=rot) kernels[rot,:,:] = pk kernels[rot+4,:,:] = dk return kernels
make_3d_kernels()[2]
array([[[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.]],
[[ 0., -1., 0.],
[ 0., 1., 0.],
[ 0., 0., 0.]],
[[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.]]])
np.rot90(dkb, k=0)
array([[-1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 0.]])
python