Simulating dynamic systems and postprocessing

Inspired by this question, I modified the script to run a simple simulation of Taylor-Green vortex using LBM and post-processes it using P5*js (an official Javascript port of the Processing API). Unfortunately, it is not performing well, but as I am 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.

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

Edit 1: Running the profiler in Chrome results in:

This shows that most of the CPU is being used by P5*JS rather than the code.