sandbox/Antoonvh/rt.c
The Rayleigh-Taylor instability and the ‘Black body’ colourbar
Inspired by the reference for the so-called ‘cool-warm’ colourbar, I visited the website of Kenneth Moreland and found that the ‘black body’ colour bar is simply stunning. We examplify it with the Rayleigh-Taylor instability:
This looks an awfull lot like this \mu-HH based video, which puts the artistic ‘freedom’ argument into its perspective.
#include "navier-stokes/centered.h"
#include "tracer.h"
#include "diffusion.h"
scalar b[];
scalar * tracers = {b};
face vector av[];
int maxlevel = 11;
long unsigned int n;
Implementation
We follow Stephane’s implementation of the ‘cool-warm’ diverging colour map and use the colour codes from Kenneth’s website.
void black_body (double cmap[NCMAP][3])
{
/* black body color map from:
* http://www.kennethmoreland.com/color-advice/
*/
static double basemap[33][3] = {
{0.0,0.0,0.0},
{0.0857913205762,0.0309874526184,0.0173328711915},
{0.133174636606,0.0588688899571,0.0346802666087},
{0.180001956037,0.0730689545154,0.0515393237212},
{0.22981556179,0.0840603593119,0.0647813713857},
{0.281397607223,0.093912584278,0.075408501413},
{0.334521638801,0.102639499627,0.0842454688083},
{0.388957802186,0.110254429637,0.0927990674821},
{0.444611925648,0.116732501721,0.101402659637},
{0.501422312285,0.122025816585,0.110058408122},
{0.559331322331,0.126067584009,0.118767796491},
{0.618285970576,0.128767919785,0.127531801155},
{0.678237857955,0.130007052818,0.136351016263},
{0.712849583079,0.181721849923,0.13081678256},
{0.743632057947,0.232649759358,0.120991817028},
{0.774324938583,0.279315911516,0.108089917959},
{0.804936242903,0.323627020047,0.0907961686083},
{0.835473266757,0.366524681419,0.0662363460741},
{0.865942668698,0.408541395043,0.026029485466},
{0.876634426153,0.46401951695,0.0173065426095},
{0.883455346031,0.518983528803,0.0149628730405},
{0.88905246237,0.572164381169,0.013499801006},
{0.893375939063,0.624108797455,0.0130334871745},
{0.89637036663,0.675180034619,0.013680092215},
{0.897973818846,0.725630730259,0.015555776796},
{0.898116710502,0.775642817733,0.0187767015864},
{0.896720396485,0.825350944866,0.023459027255},
{0.927670131094,0.859991226192,0.319086199143},
{0.956158602738,0.893933112845,0.503316730316},
{0.97827065392,0.92856476667,0.671307024002},
{0.993196411712,0.963913323002,0.83560909192},
{1.0,1.0,1.0},
};
for (int i = 0; i < NCMAP; i++) {
double x = i*(31 - 1e-10)/(NCMAP - 1);
int j = x; x -= j;
for (int k = 0; k < 3; k++)
cmap[i][k] = (1. - x)*basemap[j][k] + x*basemap[j+1][k];
}
}
Simulating the Rayleigh-Taylor instability
All steps are rather straight forward.
int main(){
periodic(left);
a = av; //Link gravity acceleration
init_grid (32);
b.gradient = minmod2; //Sharp interface advection
const face vector muc[] = {1E-6, 1E-6};
mu = muc;
run();
}
event init (t = 0){
DT = 0.01;
do{ // Initialize an unstable stratification with a sharp interface.
foreach(){
b[] = (y < Y0 + L0/2.);
u.y[] = 0.001*noise(); // Add noise to kick-off the growth of the instability
}
}while (adapt_wavelet({b}, (double[]){0.01}, maxlevel).nf);
}
event acceleration (i++){ //The effect of Gravity
foreach_face(y)
av.y[] = (b[] + b[0,-1])/2.;
}
event tracer_diffusion (i++)
diffusion (b, dt, mu);
event adapt (i++)
adapt_wavelet ((scalar*){b,u}, (double[]){0.01, 0.01, 0.01}, maxlevel);
event output (t += 0.01; t <= 1.001){ // Output a .mp4 movie and a .png image, then stop the simulation
boundary ({b});
scalar db[];
foreach(){
db[] = 0.;
foreach_dimension()
db[] += sq((b[1] - b[-1])/(2*Delta));
if (db[] > 0.)
db[] = log (sqrt (db[]) + 1.);
}
boundary({db});
output_ppm(db, file = "rt.mp4", n = (1 << (maxlevel)), map = black_body,
linear = true, box = {{0., 0.45}, {1., 0.55}}, min = 0., max = 6.);
if (t > 0.999){
output_ppm(db, file = "rt.png", n = (1 << (maxlevel)), map = black_body,
linear = true, box = {{0., 0.45}, {1., 0.55}}, min = 0., max = 6.);
foreach()
db[] *= - 1;
boundary({db});
output_ppm(db, file = "rt_banner.png", n = (1 << (maxlevel)), map = black_body,
linear = true, box = {{0., 0.47}, {1., 0.53}}, min = -6.5, max = -1);
}
}
The banner with inverted colors looks like this:
The movie reveals how the solution got into the state depicted above. The original rendering is 2048 pixels wide, we display it with half (or 1/4-th) of the pixels.