Simulating dynamic systems and postprocessing

Question

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:

enter image description here

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

• TheCoffeeCup: Hi! Welcome to Code Review. We only review working code; I know your code is working, but we cannot really answer "Anyone any idea why dot method goes wrong?"

Answer 1

I found that this specific part of the processing code was performing badly:

for (var x = 0; x < p.width; x++) {
 for (var y = 0; y < p.height; y++ ) {
 // removed irrelevant code
 p.stroke(v, 0, 1-v);
 p.point(x, y);
 }
}

It turns out that this write to the screen at every point. It is much better to buffer the writes by set-ting the pixels to certain color and then writing the whole buffer using updatePixels() like:

for (var x = 0; x < p.width; x++) {
 for (var y = 0; y < p.height; y++ ) {
 // removed irrelevant code
 p.color(v, 0, 1-v);
 p.set(x, y, c);
 }
}
p.updatePixels();

Now the draw calls take up about as much time as the collide calls (about 35%). I'm still looking to reduce this further.