summaryrefslogtreecommitdiff
path: root/field_tests
diff options
context:
space:
mode:
authorBrett Weiland <brett_weiland@gmail.com>2024-06-11 14:50:14 -0500
committerBrett Weiland <brett_weiland@gmail.com>2024-06-11 14:50:14 -0500
commitcb69732f68c0bd46c1574de16ce1aee6f38e439b (patch)
treedef1daaec81a0d4cd7b3d44b2c26b9535e07579c /field_tests
restartingHEADmaster
Diffstat (limited to 'field_tests')
-rw-r--r--field_tests/backup62
-rw-r--r--field_tests/basic_field_test.py168
-rwxr-xr-xfield_tests/field.py191
-rwxr-xr-xfield_tests/gpu_migration.py137
-rw-r--r--field_tests/kernel.c26
-rw-r--r--field_tests/makefile4
6 files changed, 588 insertions, 0 deletions
diff --git a/field_tests/backup b/field_tests/backup
new file mode 100644
index 0000000..2e60664
--- /dev/null
+++ b/field_tests/backup
@@ -0,0 +1,62 @@
+#!/usr/bin/env python3
+import numpy as np
+import matplotlib.pyplot as plt
+from alive_progress import alive_bar
+
+img_res_x = 1000
+img_res_y = 1000
+total_pixels = img_res_x * img_res_y # so we don't gotta compute it every time
+
+periods = .25
+square_x = 1
+square_y = 1
+
+xmin = (-periods * np.pi) + (square_x * np.pi)
+xmax = (periods * np.pi) + (square_x * np.pi)
+ymin = (-periods * np.pi) + (square_y * np.pi)
+ymax = (periods * np.pi) + (square_y * np.pi)
+
+escape = 10000
+iterations = 255*3
+c_x = 2 * np.pi
+c_y = 2 * np.pi
+
+
+image = np.empty([img_res_y, img_res_x])
+
+with alive_bar(img_res_y, bar = 'filling', spinner = 'waves') as bar:
+ for pix_y, y in enumerate(np.linspace(ymin, ymax, img_res_y)):
+ for pix_x, x in enumerate(np.linspace(xmin, xmax, img_res_x)):
+ on_x = x
+ on_y = y
+ for i in range(iterations):
+ completed_ratio = (((pix_x * pix_y * 1)) / total_pixels)
+ next_x = (completed_ratio * (on_x/np.sin(on_y))) + ((1 - completed_ratio) * on_x/np.tan(on_y))
+ on_y = (completed_ratio * (on_y/np.sin(on_x))) + ((1 - completed_ratio) * on_y/np.tan(on_x))
+ on_x = next_x
+ if on_x**2 + on_y**2 > escape:
+ break
+ image[pix_y][pix_x] = i
+ bar()
+
+
+
+
+plt.style.use('dark_background')
+# fuck this shit
+fig = plt.figure(frameon=False)
+fig.set_size_inches(img_res_x/fig.dpi, img_res_y/fig.dpi)
+#fig.set_size_inches(width/height, 1, forward=False)
+
+ax = plt.Axes(fig, [0., 0., 1., 1.])
+ax.set_axis_off()
+fig.add_axes(ax)
+
+cmap = plt.cm.viridis
+cmap.set_bad((0,0,0))
+cmap.set_over((0,0,0))
+cmap.set_under((0,0,0))
+
+ax.imshow(image, norm="log", aspect="auto", cmap=cmap)
+fig.savefig("linear_transform_sin_tan_arnolds_tongue_hotspot.png")
+plt.show()
diff --git a/field_tests/basic_field_test.py b/field_tests/basic_field_test.py
new file mode 100644
index 0000000..faf3c2c
--- /dev/null
+++ b/field_tests/basic_field_test.py
@@ -0,0 +1,168 @@
+#!/usr/bin/env python3
+import numpy as np
+import matplotlib.pyplot as plt
+from alive_progress import alive_bar
+
+img_res_x = 100
+img_res_y = 100
+total_pixels = img_res_x * img_res_y # so we don't gotta compute it every time
+
+periods = 1
+square_x = 0
+square_y = 0
+
+#xmin = (-periods * np.pi) + (square_x * np.pi)
+#xmax = (periods * np.pi) + (square_x * np.pi)
+#ymin = (-periods * np.pi) + (square_y * np.pi)
+#ymax = (periods * np.pi) + (square_y * np.pi)
+
+xmin = -10
+xmax = 10
+ymin = -10
+ymax = 10
+
+escape = 10000
+iterations = 255*3
+c_x = 2 * np.pi
+c_y = 2 * np.pi
+
+
+image = np.empty([img_res_y, img_res_x])
+grid = np.meshgrid(np.linspace(ymin, ymax, img_res_y), np.linspace(xmin, xmax, img_res_x))
+print(grid[0].dtype)
+
+
+class point_charge():
+ def __init__(self, x, y, c, mod):
+ self.x = x
+ self.y = y
+ self.c = c
+ self.mod = mod
+ def get_field(self, to_x, to_y):
+ if(self.mod):
+ to_x = (to_x % self.mod)
+ to_y = (to_y % self.mod)
+ return (
+ ((self.c * (self.x - to_x)) / ((self.x - to_x)**2 + (self.y - to_y)**2)**1.5),
+ ((self.c * (self.y - to_y)) / ((self.x - to_x)**2 + (self.y - to_y)**2)**1.5))
+
+#will remove all the point charge code if it turns out to be good enough to be impliemnted into openCL
+#point_charges = [point_charge(-5, -5, 100), point_charge(-5, 5, -100), point_charge(5, 0, 100)]
+point_charges = [point_charge(5,5, 100, 10), point_charge(0,0,-100, 0)]
+
+
+plt.ion()
+ax = plt.gca()
+fig = plt.gcf()
+ax.set_autoscale_on(False)
+ax.set_xlim([xmin, xmax])
+ax.set_ylim([ymin, ymax])
+
+vector_arrows = None
+
+def show_field():
+ global vector_arrows
+ grid_f = np.zeros_like(grid)
+ for p in point_charges:
+ grid_f += p.get_field(grid[0], grid[1])
+ #plt.streamplot(grid[0], grid[1], grid_f[0], grid_f[1], density=5)
+ vector_arrows = plt.quiver(grid[0], grid[1], grid_f[0], grid_f[1])
+ plt.show(block=False)
+ plt.pause(.1)
+
+
+show_field()
+
+timestep = .1
+def test_sim():
+ particle_grid = np.meshgrid(np.linspace(ymin, ymax, 100), np.linspace(xmin, xmax, 100))
+ pos = particle_grid
+ acceleration = np.zeros_like(particle_grid)
+ velocity = np.zeros_like(particle_grid)
+ velocity = [np.ones_like(particle_grid[0]) * 1, np.ones_like(particle_grid[0]) * .5]
+ mass = 10
+ charge = 1
+ particle_plot = ax.plot(velocity[0], velocity[1], 'bo', animated=True)
+ #velocity += .1
+
+ background = fig.canvas.copy_from_bbox(ax.bbox)
+ ax.draw_artist(vector_arrows)
+ fig.canvas.blit(fig.bbox)
+
+ while True:
+ fig.canvas.restore_region(background)
+ field = np.zeros_like(particle_grid)
+ # TODO can make this quicker by skipping initilization
+ for p in point_charges:
+ field += p.get_field(pos[0], pos[1])
+ acceleration = ((charge * field) / mass) * timestep
+ #print(acceleration)
+ velocity += acceleration * timestep
+ pos += velocity * timestep
+
+ fig.canvas.restore_region(background)
+ particle_plot[0].set_data(pos[0],pos[1])
+ ax.draw_artist(particle_plot[0])
+ fig.canvas.blit(fig.bbox)
+ fig.canvas.flush_events()
+ plt.pause(1/60)
+
+ #fig.canvas.draw_idle()
+test_sim()
+
+exit(1)
+
+
+
+
+
+#with alive_bar(iterations, bar = 'filling', spinner = 'waves') as bar:
+# for i in range(iterations):
+# next_x = xx / np.sin(yy)
+# yy = yy / np.sin(xx)
+# xx = next_x
+# bar()
+#image = np.vstack([xx.ravel(), yy.ravel()])
+
+
+#meshgrid makes things slower as we can't test individual points for breaking to infinity
+fractal_test = False
+if fractal_test:
+ with alive_bar(img_res_y, bar = 'filling', spinner = 'waves') as bar:
+ for pix_y, y in enumerate(np.linspace(ymin, ymax, img_res_y)):
+ for pix_x, x in enumerate(np.linspace(xmin, xmax, img_res_x)):
+ on_x = x
+ on_y = y
+ for i in range(iterations):
+ completed_ratio = (((pix_x * pix_y * 1)) / total_pixels)
+ next_x = (completed_ratio * (on_x/np.sin(on_y))) + ((1 - completed_ratio) * on_x/np.tan(on_y))
+ on_y = (completed_ratio * (on_y/np.sin(on_x))) + ((1 - completed_ratio) * on_y/np.tan(on_x))
+ on_x = next_x
+ if on_x**2 + on_y**2 > escape:
+ break
+ image[pix_y][pix_x] = i
+ bar()
+else:
+ exit()
+
+
+exit(1)
+
+plt.style.use('dark_background')
+# fuck this shit
+fig = plt.figure(frameon=False)
+fig.set_size_inches(img_res_x/fig.dpi, img_res_y/fig.dpi)
+#fig.set_size_inches(width/height, 1, forward=False)
+
+ax = plt.Axes(fig, [0., 0., 1., 1.])
+ax.set_axis_off()
+fig.add_axes(ax)
+
+cmap = plt.cm.viridis
+cmap.set_bad((0,0,0))
+cmap.set_over((0,0,0))
+cmap.set_under((0,0,0))
+
+ax.imshow(image, norm="log", aspect="auto", cmap=cmap)
+fig.savefig("linear_transform_sin_tan_arnolds_tongue_hotspot.png")
+plt.show()
diff --git a/field_tests/field.py b/field_tests/field.py
new file mode 100755
index 0000000..af2abbf
--- /dev/null
+++ b/field_tests/field.py
@@ -0,0 +1,191 @@
+#!/usr/bin/env python3
+import numpy as np
+import matplotlib.pyplot as plt
+from alive_progress import alive_bar
+
+img_res_x = 250
+img_res_y = 250
+total_pixels = img_res_x * img_res_y # so we don't gotta compute it every time
+
+periods = 1
+square_x = 0
+square_y = 0
+
+xmin = (-periods * np.pi) + (square_x * np.pi)
+xmax = (periods * np.pi) + (square_x * np.pi)
+ymin = (-periods * np.pi) + (square_y * np.pi)
+ymax = (periods * np.pi) + (square_y * np.pi)
+
+#xmin = -10
+#xmax = 10
+#ymin = -10
+#ymax = 10
+
+escape = 10000
+iterations = 255*3
+c_x = 2 * np.pi
+c_y = 2 * np.pi
+
+
+grid = np.meshgrid(np.linspace(ymin, ymax, 200), np.linspace(xmin, xmax, 200))
+#image = np.meshgrid(np.linspace(ymin, ymax, img_res_y), np.linspace(xmin, xmax, img_res_x))
+image = np.zeros([img_res_y, img_res_x], dtype=np.double)
+
+
+class point_charge():
+ def __init__(self, x, y, c, mod):
+ self.x = x
+ self.y = y
+ self.c = c
+ self.mod = mod
+ def get_field(self, to_x, to_y):
+ if(self.mod):
+ to_x = (to_x % self.mod)
+ to_y = (to_y % self.mod)
+ return np.array([
+ ((self.c * (self.x - to_x)) / ((self.x - to_x)**2 + (self.y - to_y)**2)**1.5),
+ ((self.c * (self.y - to_y)) / ((self.x - to_x)**2 + (self.y - to_y)**2)**1.5)])
+
+#will remove all the point charge code if it turns out to be good enough to be impliemnted into openCL
+#point_charges = [point_charge(-5, -5, 100), point_charge(-5, 5, -100), point_charge(5, 0, 100)]
+#point_charges = [point_charge(-1,-1, 100, 10), point_charge(1,1,-100, 0)]
+point_charges = []
+
+
+#plt.ion()
+ax = plt.gca()
+fig = plt.gcf()
+#ax.set_autoscale_on(False)
+#ax.set_xlim([xmin, xmax])
+#ax.set_ylim([ymin, ymax])
+
+vector_arrows = None
+
+def show_field():
+ global vector_arrows
+ grid_f = np.zeros_like(grid)
+ for p in point_charges:
+ grid_f += p.get_field(grid[0], grid[1])
+ #plt.streamplot(grid[0], grid[1], grid_f[0], grid_f[1], density=5)
+ vector_arrows = plt.quiver(grid[0], grid[1], grid_f[0], grid_f[1])
+ plt.show(block=False)
+ plt.pause(.1)
+
+
+#show_field()
+
+timestep = .1
+def test_sim():
+ particle_grid = np.meshgrid(np.linspace(ymin, ymax, 100), np.linspace(xmin, xmax, 100))
+ pos = particle_grid
+ acceleration = np.zeros_like(particle_grid)
+ velocity = np.zeros_like(particle_grid)
+ velocity = [np.ones_like(particle_grid[0]) * 1, np.ones_like(particle_grid[0]) * .5]
+ mass = 10
+ charge = 1
+ particle_plot = ax.plot(velocity[0], velocity[1], 'bo', animated=True)
+ #velocity += .1
+
+ background = fig.canvas.copy_from_bbox(ax.bbox)
+ ax.draw_artist(vector_arrows)
+ fig.canvas.blit(fig.bbox)
+
+ while True:
+ fig.canvas.restore_region(background)
+ field = np.zeros_like(particle_grid)
+ # TODO can make this quicker by skipping initilization
+ for p in point_charges:
+ field += p.get_field(pos[0], pos[1])
+ acceleration = ((charge * field) / mass) * timestep
+ #print(acceleration)
+ velocity += acceleration * timestep
+ pos += velocity * timestep
+
+ fig.canvas.restore_region(background)
+ particle_plot[0].set_data(pos[0],pos[1])
+ ax.draw_artist(particle_plot[0])
+ fig.canvas.blit(fig.bbox)
+ fig.canvas.flush_events()
+ plt.pause(1/60)
+
+ #fig.canvas.draw_idle()
+#test_sim()
+
+#exit(1)
+
+
+max_timesteps = 10
+
+def get_fractal_iter(img):
+ next_x = img[1][y][x] / np.sin(img[0][y][x])
+ img[0][y][x] = img[0][y][x] / np.sin(img[1][y][x])
+ img[1][y][x] = next_x
+ if (np.square(img[0][y][x]) + np.square(img[1][y][x])) >= escape:
+ z[y][x] = i
+
+
+#meshgrid makes things slower as we can't test individual points for breaking to infinity;
+#however, I will fix that later.
+cmap = plt.cm.viridis
+cmap.set_bad((0,0,0))
+cmap.set_over((0,0,0))
+cmap.set_under((0,0,0))
+
+with alive_bar(img_res_y, bar = 'filling', spinner = 'waves') as bar:
+ for pix_y, y in enumerate(np.linspace(ymin, ymax, img_res_y)):
+ for pix_x, x in enumerate(np.linspace(xmin, xmax, img_res_x)):
+ on_x = x
+ on_y = y
+ for i in range(iterations):
+ completed_ratio = (((pix_x * pix_y * 1)) / total_pixels)
+ next_x = on_x/np.sin(on_y)
+ on_y = on_y/np.sin(on_x)
+ on_x = next_x
+
+
+ # do physics here - we could just use vectors but i keep rewriting things
+ timesteps = max_
+ acceleration = np.array([0, 0], dtype=np.double)
+ velocity = np.array([on_x, on_y], dtype=np.double) # maybe multiply by stuff
+ pos = np.array([0,0], dtype=np.double)
+ for t in range(timesteps):
+ for p in point_charges:
+ acceleration += p.get_field(on_x, on_y)
+ velocity += acceleration
+ pos += velocity * (timesteps)
+ on_x += pos[0]
+ on_y += pos[1]
+
+ if on_x**2 + on_y**2 > escape:
+ image[pix_x][pix_y] = i
+ break
+ bar()
+
+ax.imshow(image, norm="log", aspect="auto", cmap=cmap)
+plt.show()
+
+
+# yeah, I shouldn't have switched to a meshgrid, oh well
+#z = np.empty_like(image[0])
+exit(0)
+with alive_bar(img_res_x, bar = 'filling', spinner = 'waves') as bar:
+ for y in range(img_res_y):
+ for x in range(img_res_x):
+ for i in range(iterations):
+ if image[0][y][x] == np.NAN:
+ continue
+ next_x = image[1][y][x] / np.sin(image[0][y][x])
+ image[0][y][x] = image[0][y][x] / np.sin(image[1][y][x])
+ image[1][y][x] = next_x
+ if (np.square(image[0][y][x]) + np.square(image[1][y][x])) >= escape:
+ z[y][x] = i
+ image[0][y][x] = np.NAN
+ image[1][y][x] = np.NAN
+ break
+# #do physics here
+
+
+
+ bar()
+#image = np.clip(image, -escape, escape)
+
diff --git a/field_tests/gpu_migration.py b/field_tests/gpu_migration.py
new file mode 100755
index 0000000..718fe7d
--- /dev/null
+++ b/field_tests/gpu_migration.py
@@ -0,0 +1,137 @@
+#!/usr/bin/env python3
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+import pyopencl as cl
+from alive_progress import alive_bar
+from matplotlib.backend_bases import MouseButton
+
+img_res_x = 1000
+img_res_y = 1000
+total_pixels = img_res_x * img_res_y # so we don't gotta compute it every time
+
+periods = 1
+square_x = 0
+square_y = 0
+
+x_min = (-periods * np.pi) + (square_x * np.pi)
+x_max = (periods * np.pi) + (square_x * np.pi)
+y_min = (-periods * np.pi) + (square_y * np.pi)
+y_max = (periods * np.pi) + (square_y * np.pi)
+
+escape = 10000
+iterations = 255*3
+c_x = 2 * np.pi
+c_y = 2 * np.pi
+
+animation_progres_save = "ani1_less.mp4"
+single_frame_save = "asdf.png"
+
+opencl_context = cl.create_some_context(interactive=False)
+opencl_queue = cl.CommandQueue(opencl_context)
+
+kernel_src_path = "./kernel.c"
+
+frames = 1
+rendered_frames = []
+
+image = np.empty([img_res_x, img_res_y], np.uint32)
+image_buffer = cl.Buffer(opencl_context, cl.mem_flags.WRITE_ONLY, image.nbytes)
+
+plt.style.use('dark_background')
+# fuck this shit
+fig = plt.figure(frameon=False)
+fig.set_size_inches(img_res_x/fig.dpi, img_res_y/fig.dpi)
+#fig.set_size_inches(width/height, 1, forward=False)
+
+
+ax = plt.Axes(fig, [0., 0., 1., 1.])
+#ax.set_axis_off()
+fig.add_axes(ax)
+
+cmap = plt.cm.viridis
+cmap.set_bad((0,0,0))
+cmap.set_over((0,0,0))
+cmap.set_under((0,0,0))
+
+bruh = None
+
+def on_click(event):
+ global bruh
+ split_ratio = 1
+ if (event.button is MouseButton.MIDDLE) and (event.inaxes):
+ # there's probably a way to set global coordinates;
+ # I don't expect this to go anywhere so I don't really care
+ on_x = ((event.xdata / img_res_x) * abs(x_max - x_min)) + x_min
+ on_y = ((event.ydata / img_res_y) * abs(y_max - y_min)) + y_min
+
+ #I, uh, also don't know the best way to replicate the openCL code automaticly in python
+ #so ajust if nessesary
+ x_hops = []
+ y_hops = []
+ for i in range(iterations):
+ x_hops.append(((on_x - x_min) / abs(x_max - x_min)) * img_res_x)
+ y_hops.append(((on_y - y_min) / abs(y_max - y_min)) * img_res_y)
+ next_x = on_x/np.tan(on_y)
+ on_y = on_y/np.tan(on_x)
+ on_x = next_x
+ if on_x**2 + on_y**2 > escape:
+ break
+ print(y_hops[0])
+ print(on_y)
+ print("{} hops".format(len(x_hops)))
+ if bruh:
+ bruh.pop(0).remove()
+ bruh = ax.plot(x_hops, y_hops)
+ plt.draw()
+
+
+
+ print(on_x, on_y)
+
+
+print("compiling openCL kernel...")
+with open(kernel_src_path, 'r') as kernel_src:
+ compiled_kernel = cl.Program(opencl_context, kernel_src.read()).build()
+
+encoding_progress = alive_bar(frames, bar = 'filling', spinner = 'waves')
+
+def display_encoder_progress(current_frame: int, total_frames: int):
+ print("Encoding: frame {}/{}".format(current_frame, frames))
+
+
+
+print("Rendering {} frames...".format(frames))
+if frames > 1:
+ with alive_bar(frames, bar = 'filling', spinner = 'waves') as bar_total:
+ for frame_i in range(0, frames):
+ compiled_kernel.render_frame(opencl_queue, image.shape, None, image_buffer,
+ np.double(abs(x_max - x_min) / img_res_x),
+ np.double(abs(y_max - y_min) / img_res_y),
+ np.double(x_min), np.double(y_min),
+ np.uint32(iterations), np.uint32(escape),
+ np.double(frame_i / frames))
+
+
+ cl.enqueue_copy(opencl_queue, image, image_buffer).wait()
+ rendered_frame = ax.imshow(image, norm="log", aspect="auto", cmap=cmap, animated="True")
+ rendered_frames.append([rendered_frame])
+ bar_total()
+ print("Encoding/Saving...")
+ ani = animation.ArtistAnimation(fig, rendered_frames, interval=30, blit=True)
+ ani.save(animation_progres_save, extra_args=['-preset', 'lossless'], progress_callback=display_encoder_progress, codec="h264_nvenc")
+else:
+ compiled_kernel.render_frame(opencl_queue, image.shape, None, image_buffer,
+ np.double(abs(x_max - x_min) / img_res_x),
+ np.double(abs(y_max - y_min) / img_res_y),
+ np.double(x_min), np.double(y_min),
+ np.uint32(iterations), np.uint32(escape),
+ np.double(1))
+ cl.enqueue_copy(opencl_queue, image, image_buffer).wait()
+ ax.imshow(image, norm="log", aspect="auto", cmap=cmap)
+ fig.savefig(single_frame_save)
+ plt.autoscale(False)
+ plt.connect('motion_notify_event', on_click)
+ plt.show()
+
+
diff --git a/field_tests/kernel.c b/field_tests/kernel.c
new file mode 100644
index 0000000..08b9137
--- /dev/null
+++ b/field_tests/kernel.c
@@ -0,0 +1,26 @@
+__kernel void basic_test() {
+ printf("this cores ID: (%lu, %lu)\n", get_global_id(0), get_global_id(1));
+}
+
+__kernel void render_frame(__global unsigned int *frame_output,
+ double x_step, double y_step,
+ double x_start, double y_start,
+ unsigned int iterations, unsigned int escape, double ratio) {
+ unsigned int result;
+ double on_x = (get_global_id(0) * x_step) + x_start;
+ double on_y = (get_global_id(1) * y_step) + y_start;
+ double next_x;
+ unsigned int iter;
+
+ orig_x = on_x;
+ orig_y = on_y;
+ for(iter = 0; iter < iterations; iter++) {
+ if(orig_x == 123 && orig_y == 5) printf("%d, %d\n", orig_x,orig_y);
+ next_x = on_x / sin(on_y);
+ on_y = on_y / sin(on_x);
+ on_x = next_x;
+ if((pow(on_x, 2) + pow(on_y, 2)) >= escape) break;
+ }
+
+ frame_output[(get_global_id(1) * get_global_size(1)) + get_global_id(0)] = iter;
+}
diff --git a/field_tests/makefile b/field_tests/makefile
new file mode 100644
index 0000000..539deb1
--- /dev/null
+++ b/field_tests/makefile
@@ -0,0 +1,4 @@
+make:
+ gcc -Wall -fpic -c field.c
+ gcc -Wall -shared -o field.so field.o
+ python3 field.py