Simple simulation of Cahn Hilliard equation

cahn_hilliard.c

// Compile with:
// $ gcc -Ofast -march=native -fopenmp ./cahn_hilliard.c -o cahn_hilliard -lm
// View the simulation using mpv:
// $ ./cahn_hilliard | mpv - --scale=ewa_lanczossoft --fs

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <stdbool.h>

#define W (1920 / 4)
#define H (1080 / 4)

// Convenience defines for OpenMP
#define FOREACH_POINT \
    _Pragma("omp parallel for collapse(2) schedule(guided)")\
        for(int i = 0; i < H; i ++) for(int j = 0; j < W; j ++)

#define SECTIONS \
    _Pragma("omp sections")
#define SECTION \
    _Pragma("omp section")

// Helper macro or accessing the grid, taking boundary wrapping into account
#define F(__dx, __dy) (grid[(x + H + __dx) % H][(y + W + __dy) % W])

// Discretization parameters
static const int SKIP_FRAMES = 10; // Steps between each output frame
static const double TIME_STEP = 0.000001;
static const double GRID_STEP = 0.1;

// Fairly arbitrary function that turns a double value into a color on a gradient
static void color(int x, int y, const double grid[H][W], char * c_val){
    double c = grid[x][y];

    // Sharpness of gradient
    const double SHARP = 10;
    // Adjusts the tint of negative and positive values
    const double TINT_A = 0.7;
    const double TINT_B = 0.9;

    // Use tanh to squish the value into (-1, 1) and fabs to avoid negatives
    double q = fabs(tanh(c * SHARP));

    c_val[0] = 0.5 + 255.0 * q;
    c_val[1] = 0.5 + 255.0 * (q * (TINT_A * (c < 0) + (1-TINT_A) * (c > 0)));
    c_val[2] = 0.5 + 255.0 * (q * (TINT_B * (c > 0) + (1-TINT_B) * (c < 0)));
}

static double step(int x, int y, const double grid[H][W]){
    // Square root of width of transition regions
    const double B = 3;

    double u = F(0,0);
    double dx, dy, abs_grad_sq;
    double laplacian, biharmonic;

    // Orthogonally adjacent grid points
    double orth1 = F(1, 0) + F(-1, 0) + F(0, 1) + F(0, -1);
    // Diagonally adjacent grid points
    double diag = F(1, 1) + F(1, -1) + F(-1, 1) + F(-1, -1);
    // Distance 2 orthogonal points
    double orth2 = F(2, 0) + F(0, 2) + F(-2, 0) + F(0, -2);

    // First derivatives approximated by finite differences (second order)
    dx = (-F(2, 0) + 8 * F(1, 0) - 8 * F(-1, 0) + F(-2, 0)) / (12.0 * GRID_STEP);
    dy = (-F(0, 2) + 8 * F(0, 1) - 8 * F(0, -1) + F(0, -2)) / (12.0 * GRID_STEP);

    // Five point stencil
    laplacian = (orth1 - 4 * u) / pow(GRID_STEP, 2);

    // Stencil for biharmonic operator
    biharmonic = (orth2 + 2*diag - 8*orth1 + 20*u) / pow(GRID_STEP, 4);

    // Squared absolute value of gradient
    abs_grad_sq = pow(dx, 2) + pow(dy, 2);

    // Forward Euler timestepping for Cahn-Hilliard equation
    return u + TIME_STEP * (3*u*(2*abs_grad_sq + laplacian) - laplacian - B * biharmonic);
}

int main(void){
    static double grid[2][H][W];

    // Maximum possible length of PPM header
    const int MAX_HEADER_LEN = 128;

    // Entire PPM output
    char frame[W * H * 3 + MAX_HEADER_LEN];

    // Points to the actual bitmap part of output
    char *image;

    size_t header_len = snprintf(frame, MAX_HEADER_LEN, "P6\n%d %d\n255\n", W, H);
    size_t frame_len = 3 * W * H + header_len;

    // Image starts right after the header, so we just shift the pointer
    image = frame + header_len;


    srand48(time(NULL));

    // Randomize the initial conditions
    FOREACH_POINT
        grid[0][i][j] = 2*drand48() - 1;

    while(true) SECTIONS {
        // Output the grid
        SECTION {
            FOREACH_POINT
                color(i, j, grid[0], &image[3 * (i * W + j)]);

            fwrite(frame, 1, frame_len, stdout);
        }

        SECTION {
            // Run time-stepping
            for(int k = 0; k < SKIP_FRAMES; k ++) for(int flag = 0; flag < 2; flag ++)
                FOREACH_POINT
                    grid[!flag][i][j] = step(i, j, grid[flag]);
        }
    }

    return 0;
}