Inspired by this [question][1], I modified the script to run a simple simulation of [Taylor-Green vortex][2] using [LBM][3] and post-processes it using [P5*js][4] (an official js port of the [Processing][5] API). Unfortunately, it is not performing well, but as i am a new to Javascript programming, I lack the insight to really optimize the code. Getting some advice and suggestions for improvement would be greatly appreciated.
A working codepen version can be found [here][6].
// 2D vector class
function vec2(x,y){
this.x = x;
this.y = y;
}
vec2.prototype = {
scale: function(s){
return new vec2(this.x*s,this.y*s);
},
add: function(v){
return new vec2(this.x+v.x,this.y+v.y);
},
subtract: function(v){
return new vec2(this.x-v.x,this.y-v.y);
},
dot: function(v){
return this.x*v.x+this.y+v.y;
},
length: function(){
return Math.sqrt(this.x*this.x+this.y*this.y);
},
normalise: function(){
var l = 1/this.length();
return new vec2(this.x*l, this.y*l);
},
};
// domain class
function domain(nx, ny) {
this.nx = nx; // domain width
this.ny = ny; // domain height
this.f = []; // distribution array
this.ftmp = []; // temporary array
this.dens = []; // density array
this.vel = []; // velocity array
this.omega = 1; // relaxation frequency
this.e = [ // discrete velocity set
new vec2(0,0),
new vec2(0,1), new vec2(1,0), new vec2(0,-1), new vec2(-1,0),
new vec2(1,1), new vec2(1,-1), new vec2(-1,-1), new vec2(-1,1)
];
this.w = [ // weights
4/9,
1/9, 1/9, 1/9, 1/9,
1/36, 1/36, 1/36, 1/36
];
// Arrays initialization
for (var x=0; x<this.nx; x++) {
this.f[x] = [];
this.ftmp[x] = [];
this.dens[x] = [];
this.vel[x] = [];
for (var y=0; y<this.ny; y++) {
this.f[x][y] = [];
this.ftmp[x][y] = [];
}
}
}
domain.prototype = {
init: function(){
// Initializes Taylor-Green vortex
var kx = 2*Math.PI/this.nx;
var ky = 2*Math.PI/this.ny;
var kxkx = kx*kx;
var kyky = ky*ky;
var ksq = kxkx + kyky;
var k = Math.sqrt(ksq);
var dens0 = 1;
this.umax = 0.1;
var u0 = 4*this.umax;
this.densmax = dens0 + 3*dens0*u0*u0/4;
for (var x=0; x<this.nx; x++){
for (var y=0; y<this.ny; y++){
var u = u0*ky/k*Math.cos(kx*x)*Math.sin(ky*y);
var v = -u0*kx/k*Math.sin(kx*x)*Math.cos(ky*y);
this.dens[x][y] = dens0 + 3*dens0*u0*u0/4*(kyky/ksq*Math.cos(2*kx*x)+kxkx/ksq*Math.sin(2*ky*y));
this.vel[x][y] = new vec2(u,v);
for(var i=0; i<9; i++){
// Initialize using equilibrium distribution
var uu = this.vel[x][y].x*this.vel[x][y].x + this.vel[x][y].y*this.vel[x][y].y;
var eu = this.e[i].x*this.vel[x][y].x + this.e[i].y*this.vel[x][y].y;
this.f[x][y][i] = this.w[i]*this.dens[x][y]*(1+3*eu+4.5*eu*eu-1.5*uu);
}
}
}
},
collide: function(){
for(var x=1; x<this.nx-1; x++){
for(var y=1; y<this.ny-1; y++){
// calculate density
var rho = 0;
for(var i=0; i<9; i++){
rho += this.f[x][y][i];
}
this.dens[x][y] = rho;
// calculate velocity
var u = new vec2(0,0);
for(var i=1; i<9; i++){
u = u.add( this.e[i].scale( this.f[x][y][i] ) );
}
u = u.scale( 1/rho );
this.vel[x][y] = u;
// Perform collision step and save to temp array ftmp
var uu = u.x*u.x + u.y*u.y;
for(var i=0; i<9; i++){
var eu = u.x*this.e[i].x + u.y*this.e[i].y;
var fiEq = this.w[i]*rho*(1+3*eu+4.5*eu*eu-1.5*uu);
var fiCol = -this.omega*(this.f[x][y][i]-fiEq); // bgk
this.ftmp[x][y][i] = this.f[x][y][i] + fiCol;
}
}
}
},
periodic: function(){
// Apply periodic boundary conditions on ftmp
// x-periodic
for(var y=1; y<this.ny-1; y++){
for(var i=0; i<9; i++){
this.ftmp[0][y][i] = this.ftmp[this.nx-2][y][i];
this.ftmp[this.nx-1][y][i] = this.ftmp[1][y][i];
}
}
// y-periodic
for(var x=1; x<this.nx-1; x++){
for(var i=0; i<9; i++){
this.ftmp[x][0][i] = this.ftmp[x][this.ny-2][i];
this.ftmp[x][this.ny-1][i] = this.ftmp[x][1][i];
}
}
// corner treatment
for(var i=0; i<9; i++){
this.ftmp[0][0][i] = this.ftmp[this.nx-2][this.ny-2][i];
this.ftmp[this.nx-1][this.ny-1][i] = this.ftmp[1][1][i];
this.ftmp[this.nx-1][0][i] = this.ftmp[1][this.ny-2][i];
this.ftmp[0][this.ny-1][i] = this.ftmp[this.nx-2][1][i];
}
},
stream: function(){
// Perform streaming step ftmp -> f
for(var x=1; x<this.nx-1; x++){
for(var y=1; y<this.ny-1; y++){
for(var i=0; i<9; i++){
this.f[x][y][i] = this.ftmp[x-this.e[i].x][y-this.e[i].y][i];
}
}
}
}
}
function simulation(){
var sim = function(p) {
var nx = 200, ny = 200;
var myDomain = new domain(nx, ny);
p.setup = function() {
p.createCanvas(nx, ny)
.parent('sim');
//p.frameRate(30);
myDomain.init();
p.noStroke();
p.colorMode(p.RGB, 1);
}
p.draw = function() {
myDomain.collide();
myDomain.periodic();
myDomain.stream();
var dens = myDomain.dens;
var vel = myDomain.vel;
var densmax = myDomain.densmax;
var umax = myDomain.umax;
for (var x = 0; x < p.width; x++) {
for (var y = 0; y < p.height; y++ ) {
//var v = dens[x][y]/densmax;
var v = vel[x][y].length()/umax;
p.stroke(v, 0, 1-v);
p.point(x, y);
}
}
}
}
return new p5(sim);
}
$( document ).ready(simulation());
[1]: http://physics.stackexchange.com/questions/208924/problems-implementing-lattice-boltzmann-method
[2]: https://en.wikipedia.org/wiki/Taylor%E2%80%93Green_vortex
[3]: https://en.wikipedia.org/wiki/Lattice_Boltzmann_methods
[4]: http://p5js.org/
[5]: https://processing.org/
[6]: http://codepen.io/nlooije/pen/MaoEMJ