diff --git a/opengl/glfw_compute_shader_heat_equation.c b/opengl/glfw_compute_shader_heat_equation.c index 743fa00..fd221d6 100644 --- a/opengl/glfw_compute_shader_heat_equation.c +++ b/opengl/glfw_compute_shader_heat_equation.c @@ -1,22 +1,21 @@ /* -Compute shader that increments an SSBO, and copies it to a texture that is then displayed. +Heat equation. TODO CPU only for now. Usage: - ./prog [width [work_group_width [cpu]]]] + ./prog [width [window_width [work_group_width [cpu [conduction_coeff [select_boundary]]]]] + +Example: + + ./prog 512 512 16 1 0.5 2 - width: window width - work_group_width. Must divide width. - cpu: if '1', use CPU, else GPU. If set, work_group_width is ignored. -TODO The GPU was much slower than CPU on my tests. Am I doing something wrong (global memory access?) -or is it just a GPU limitation of reading / writting the temperature vector? +It is hard to view function values with colors like this, because the range is too limited. -## Work groups - -Work group size is important in this example. - -Work groups of size 1 are very slow (60FPS). 16 by 16 was 1400FPS! +A solution would be to us a 3D plot, but that might require moving the camera around to see certain spots... */ #include "common.h" @@ -33,6 +32,7 @@ static const GLuint indices[] = { 0, 1, 2, 0, 2, 3, }; +static const float conduction_coeff = 0.5; static const GLchar *vertex_shader_source = "#version 330 core\n" @@ -77,9 +77,76 @@ void error_callback(int error, const char* description) { puts(description); } +void init_boundary(GLfloat *temperatures, size_t width, size_t height, int which) { + const float PI2 = 2.0 * acos(-1.0); + const unsigned int sin_periods = 2; + switch (which) { + case 1: + /* Linear decrease. + * + * 1.0 | 0.5 + * --------- + * 0.5 | 0.0 + * */ + for (size_t i = 0; i < height; ++i) { + temperatures[i * width + 0 ] = 0.5 + (0.5 * (i / (float)height)); + temperatures[i * width + width - 1] = 0.5 * (i / (float)height); + } + for (size_t j = 1; j < width - 1; ++j) { + temperatures[ j] = 0.5 * (1.0 - (j / (float)width)); + temperatures[(height - 1) * width + j] = 0.99 - (0.5 * (j / (float)width)); + } + break; + + case 2: + /* sin: + * + * 0.5 | 0.5 + * --------- + * 0.5 | 0.5 + * + * Does N periods in the middle. + * */ + for (size_t i = 0; i < height; ++i) { + float f = 0.5 * (0.99 + sin(sin_periods * PI2 * (i / (float)height))); + temperatures[i * width + 0 ] = f; + temperatures[i * width + width - 1] = f; + } + for (size_t j = 1; j < width - 1; ++j) { + float f = 0.5 * (0.99 + sin(sin_periods * PI2 * (j / (float)width))); + temperatures[ j] = f; + temperatures[(height - 1) * width + j] = f; + } + break; + + case 0: + default: + /* Const at 1.0. + * + * 1.0 | 0.5 + * --------- + * 1.0 | 0.5 + * + * */ + for (size_t i = 0; i < height; ++i) { + temperatures[i * width + 0 ] = 0.99; + temperatures[i * width + width - 1] = 0.99; + } + for (size_t j = 1; j < width - 1; ++j) { + temperatures[ j] = 0.99; + temperatures[(height - 1) * width + j] = 0.99; + } + break; + } +} + int main(int argc, char **argv) { GLFWwindow *window; - GLfloat *temperatures; + GLfloat + *temperatures = NULL, + *temperatures2 = NULL, + *temperature_buf = NULL + ; GLint coord2d_location, textureSampler_location, @@ -90,15 +157,18 @@ int main(int argc, char **argv) { compute_program, ebo, height, + window_height, program, ssbo, texture, width, + window_width, work_group_width, vao, vbo ; - int cpu; + int cpu, which_boundary; + float conduction_coeff; size_t n_temperatures; char *compute_shader_source, *work_group_width_str; @@ -110,21 +180,50 @@ int main(int argc, char **argv) { } height = width; if (argc > 2) { - work_group_width = strtol(argv[2], NULL, 10); + window_width = strtol(argv[2], NULL, 10); + } else { + window_width = WIDTH; + } + window_height = window_width; + if (argc > 3) { + work_group_width = strtol(argv[3], NULL, 10); } else { work_group_width = WORK_GROUP_WIDTH; } - if (argc > 3) { - cpu = (argv[3][0] == '1'); + if (argc > 4) { + cpu = (argv[4][0] == '1'); } else { cpu = 0; } + /* Must be between 0.0 and 1.0. + * + * Physics allows it to be in 0 / infinity. + * + * Anything greater than 1.0 leads to numeric instabilities + * for our simplistic method, as equilibrium is reached in less + * time than one time step, and goes over. + * + * Negative values make temperatures unbounded and breaks energy conservation. + * + * But you obviously will try out "bad" values in the simulation to see what happens. + * The behaviour of this value around 1.99, 2.0, 2.01, 3.0 is specially interesting. + * */ + if (argc > 5) { + conduction_coeff = strtod(argv[5], NULL); + } else { + conduction_coeff = conduction_coeff; + } + if (argc > 6) { + which_boundary = strtol(argv[6], NULL, 10); + } else { + which_boundary = 0; + } /* Window. */ glfwInit(); glfwSetErrorCallback(error_callback); glfwWindowHint(GLFW_RESIZABLE, GL_FALSE); - window = glfwCreateWindow(width, height, __FILE__, NULL, NULL); + window = glfwCreateWindow(window_width, window_height, __FILE__, NULL, NULL); glfwMakeContextCurrent(window); glfwSwapInterval(1); glewInit(); @@ -179,18 +278,22 @@ int main(int argc, char **argv) { srand(time(NULL)); n_temperatures = width * height; temperatures = malloc(n_temperatures * sizeof(temperatures[0])); - for (size_t i = 0; i < n_temperatures; ++i) { - /* Doable with just a simple fragment shader + time uniform, but easier to see what is going on. */ - temperatures[i] = 0.0; - /* Not doable with simple fragment shader, so closer to real applications. */ - /*temperatures[i] = rand() / (float)RAND_MAX;*/ + /* Initial condition. TODO: make continuous with boundary conditions. */ + for (size_t i = 1; i < height - 1; ++i) { + for (size_t j = 1; j < width - 1; ++j) { + temperatures[i * width + j] = 0.0; + } } - if (!cpu) { + init_boundary(temperatures, width, height, which_boundary); + if (cpu) { + temperatures2 = malloc(n_temperatures * sizeof(temperatures[0])); + /* Boundary must also be initialized for this buffer, + * since the boundary is never touched after the beginning. */ + init_boundary(temperatures2, width, height, which_boundary); + } else { glGenBuffers(1, &ssbo); glBindBuffer(GL_SHADER_STORAGE_BUFFER, ssbo); - /* GL_DYNAMIC_COPY because we are copying to GPU memory several times. */ glBufferData(GL_SHADER_STORAGE_BUFFER, n_temperatures * sizeof(temperatures[0]), temperatures, GL_DYNAMIC_COPY); - /* 0 corresponds to the 0 of the shader at: "layout (std430, binding=0)". */ glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, ssbo); glBindBuffer(GL_SHADER_STORAGE_BUFFER, 0); free(temperatures); @@ -209,25 +312,48 @@ int main(int argc, char **argv) { } /* Constant state. */ - glViewport(0, 0, width, height); + glViewport(0, 0, window_width, window_height); glClearColor(1.0f, 1.0f, 1.0f, 1.0f); /* Main loop. */ common_fps_init(); while (!glfwWindowShouldClose(window)) { if (cpu) { - for (unsigned int i = 0; i < height; ++i) { - for (unsigned int j = 0; j < width; ++j) { - GLfloat *t = &temperatures[i * width + j]; - *t = fmod(*t + 0.01, 1.0); + for (unsigned int i = 1; i < height - 1; ++i) { + for (unsigned int j = 1; j < width - 1; ++j) { + /* TODO: find the name of this naive integration. + * + * TODO how imprecise it it? + * + * One limitation of this method is that temperature changes can + * only move one square at a time. + * + * In the real word, a temperature change on one side of the room + * immediately affects the other side (very little). + * + * And a denser grid implies slower propagation which makes no sense. + * + * The decent method will likely require matrix multiplications? + * */ + size_t idx = i * width + j; + temperatures2[idx] = + (1.0 - conduction_coeff) * temperatures2[idx] + + conduction_coeff * ( + temperatures[idx - 1] + + temperatures[idx + 1] + + temperatures[idx - width] + + temperatures[idx + width] + ) / 4.0; } } - /* GL_RED because there is no better way to do grayscale. - * http://stackoverflow.com/questions/680125/can-i-use-a-grayscale-image-with-the-opengl-glteximage2d-function */ glTexImage2D( GL_TEXTURE_2D, 0, GL_RED, width, height, - 0, GL_RED, GL_FLOAT, temperatures + 0, GL_RED, GL_FLOAT, temperatures2 ); + /* Swap old and new. */ + temperature_buf = temperatures; + temperatures = temperatures2; + temperatures2 = temperature_buf; } else { /* Compute. */ glUseProgram(compute_program); @@ -251,10 +377,10 @@ int main(int argc, char **argv) { /* Cleanup. */ glDeleteBuffers(1, &ebo); - if (!cpu) { - glDeleteBuffers(1, &ssbo); - } else { + if (cpu) { free(temperatures); + } else { + glDeleteBuffers(1, &ssbo); } glDeleteBuffers(1, &vbo); glDeleteVertexArrays(1, &vao); diff --git a/opengl/glfw_compute_shader_ssbo_inc.c b/opengl/glfw_compute_shader_ssbo_inc.c index b283c71..d19fe43 100644 --- a/opengl/glfw_compute_shader_ssbo_inc.c +++ b/opengl/glfw_compute_shader_ssbo_inc.c @@ -249,10 +249,10 @@ int main(int argc, char **argv) { /* Cleanup. */ glDeleteBuffers(1, &ebo); - if (!cpu) { - glDeleteBuffers(1, &ssbo); - } else { + if (cpu) { free(temperatures); + } else { + glDeleteBuffers(1, &ssbo); } glDeleteBuffers(1, &vbo); glDeleteVertexArrays(1, &vao); diff --git a/opengl/offscreen.c b/opengl/offscreen.c index cbaf8ff..aba103b 100644 --- a/opengl/offscreen.c +++ b/opengl/offscreen.c @@ -112,6 +112,7 @@ static void screenshot_png(const char *filename, unsigned int width, unsigned in png_write_info(png, info); png_write_image(png, *png_rows); png_write_end(png, NULL); + png_destroy_write_struct(&png, &info); fclose(f); } #endif diff --git a/png/open_manipulate_write.c b/png/open_manipulate_write.c index 65b9f51..9c1edac 100644 --- a/png/open_manipulate_write.c +++ b/png/open_manipulate_write.c @@ -85,6 +85,7 @@ static void write_png_file(char *filename) { /*png_set_filler(png, 0, PNG_FILLER_AFTER);*/ png_write_image(png, row_pointers); png_write_end(png, NULL); + png_destroy_write_struct(&png, &info); for (y = 0; y < height; y++) { free(row_pointers[y]); }