png deallocate properly (credits to mgmalheiros), improve heat equation

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);

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);

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

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]);
     }