Skip to content

Instantly share code, notes, and snippets.

@laserbat
Last active March 22, 2021 17:24
Show Gist options
  • Save laserbat/d6be8a145f246df77584ed17d8df8f3f to your computer and use it in GitHub Desktop.
Save laserbat/d6be8a145f246df77584ed17d8df8f3f to your computer and use it in GitHub Desktop.
Simple simulation of Cahn Hilliard equation
// 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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment