Untitled

//Demonstration of the shallow water equations in a channel.  This demo has
//periodic boundary conditions in the East/West and walls in the North/South.

tributary.trace = false;
//rectangle area parameters
var xT = 30;
var yT = 30;

var W = 1024;
var H = 388;

//graphic parameters
var backgroundColor = "#F1F1F1";
var numberOfColors = 64

//number of points to use in solution line
var m = 40
var n = 40

//domain sizes, in km
var length = 6000.
var width = 4400.

//initial condition parameters
var h0 = 2000
var h1 = -220.
var h2 = 133.

// Gravity, m/s^2
var gp = 9.81
// Coriolis factor
var f = 1.0e-4
// Beta-plane 1/(m s)
var beta = 1.5e-11

//time stretching
var timeStretch = 10
var colorMap = null;

//var vizMode = "phi";
//var vizMode = "u";
//var vizMode = "v";
var vizMode = "uv";
var oldt, u_old, v_old, phi_old, dx, dy, realT, k, u_new,v_new,phi_new;
tributary.init = function(g) {

  oldt  = null
  u_old = null
  v_old = null
  phi_old = null
  dx = 2*length*1000/ (m-2)
  dy = 2*width*1000/ (n-2)
  realT = 0.
  k = 0
  u_new = create2DArray(m, n)
  v_new = create2DArray(m, n)
  phi_new = create2DArray(m, n)
  //initialConditions(u_new, v_new, phi_new)
  //getColorMap();

  //initialize svg stuff
  var inds = d3.range(n*m);
  var cells = g.selectAll("g.cell")
    .data(inds)
    .enter()
    .append("g").classed("cell", true)
   //cells.attr("transform",
   cells
    .append("rect");

};

tributary.run = function(g,t) {

  if (u_old != null) {
    dt = timeStretch*10000*(t - oldt)
    //advanceModel(u_new, v_new, phi_new, u_old, v_old, phi_old, k, dt)
    k++;
  }
  else {
    dt = 0.;
  }

  //drawShallowWater(ctx, u_new, v_new, phi_new);
  //updateValues();

  realT = realT + dt;
  oldt  = t;
};

//drawing functions
function drawShallowWater(ctx, u_new, v_new, phi_new) {
  if (vizMode == "phi") {
    vmin = 14000;
    vmax = 24000;
    vsrc = phi_new
  }
  else if (vizMode == "u") {
    vmin = -30;
    vmax = 10;
    vsrc = u_new
  }
  else if (vizMode == "v") {
    vmin = -20;
    vmax = 20;
    vsrc = v_new;
  }
  else {
    vmin = 0;
    vmax = 30;
    // vsrc specified below
  }

  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++) {
      if (vizMode == "uv") {
        lv = Math.sqrt(Math.pow(u_new[i][j],2) + Math.pow(v_new[i][j],2)) - vmin;
      }
      else {
        lv = vsrc[i][j] - vmin;
      }
      hv = lv/(vmax-vmin)
      ci = Math.min(Math.max(Math.round(numberOfColors*hv),0),numberOfColors-1)
      if (!isNaN(ci)) {
        color = colorMap[ci];
        ctx.fillStyle = RGB2Color(color)
      }
      else {
        ctx.fillStyle = "#000000"
      }

      drawRectT(ctx, (i + 0.)/m, (n - j - 1.)/n, 1./m, 1./n)
    }
  }
}

function getColorMap() {
  colorMap = create2DArray(numberOfColors, 3);

  nc4 = numberOfColors /4;
  // red
  colorTriangle(0, 3*nc4, nc4, numberOfColors);
  // green
  colorTriangle(1, 2*nc4, nc4, numberOfColors);
  // blue
  colorTriangle(2, 1*nc4, nc4, numberOfColors);
}

function colorTriangle(colorIndex, start, width, length) {
  lb  = start - 3*width/2 - 1;
  mb1 = start - width/2 - 1;
  mb2 = start + width/2 - 1;
  rb  = start + 3*width/2 - 1;

  for (i = 0; i < length; i++) {
    // rising regime
    if (i <= mb1 && i > lb) {
      colorMap[i][colorIndex] = lininterp(i, lb, mb1, 0, 1);
    }
    // one regime
    else if (i <= mb2 && i > mb1) {
      colorMap[i][colorIndex] = 1.;
    }
    // falling regime
    else if (i <= rb && i > mb2) {
      colorMap[i][colorIndex] = lininterp(i, mb2, rb, 1, 0);
    }
    // zero regime
    else {
      colorMap[i][colorIndex] = 0.;
    }
  }
}

//model functions
function initialConditions(u_new, v_new, phi_new) {
  hx = 1./(4.*dx)
  hy = 1./(4.*dy)
  hhx = 2.*hx
  hhy = 2.*hy
  pi = Math.PI
  mlength = length*1000
  mwidth = width*1000

  for (i=0; i < m; i++) {
    for (j=0; j < n; j++) {
      xix = (i - 1)*dx/2.
      triarg = 2.*pi*xix/mlength
      xtrig = Math.sin(triarg)
      dtrig = (2.*pi/mlength)*Math.cos(triarg)

      yiy = (j - 1)*dy/2.
      fr = f + beta*(yiy-mwidth/2.)

      co = 9.*((mwidth/2.) - yiy)/(2.*mwidth)

      ab1 = cosh(co)
      bb = 1./(ab1*ab1)
      dco = -4.5/mwidth
      derth = bb*dco

      aa = tanh(co)
      dersec = -2.*aa*derth

      u_new[i][j] = -gp*(h1*derth + h2*xtrig*dersec )/fr
      v_new[i][j] = gp*(h2*bb*dtrig)/fr
      phi_new[i][j] = gp*(h0 + h1*aa + h2*bb*xtrig)
    }
  }

  setBoundaryConditions(u_new, v_new, phi_new, m, n)
}

function advanceModel(u_new, v_new, phi_new, u_old, v_old, phi_old, k, dt) {

  if (u_old == null || v_old == null || phi_old == null)
    return;
  if (k == 0) {
    dtv = dt;
    kv = 1;
  }
  else {
    kv = 0;
    dtv = 2*dt;
  }

  hx = 1./(4.*dx)
  hy = 1./(4.*dy)
  hhx = 2.*hx
  hhy = 2.*hy

  for (i=1; i < m-1; i++) {
    for (j=1; j < n-1; j++) {
      ip1=i+1
      im1=i-1
      jp1=j+1
      jm1=j-1

      aa = (u_old[ip1][j][1] + u_old[im1][j][1])*
           (u_old[ip1][j][1] - u_old[im1][j][1])*hx

      bb = (v_old[i][jp1][1] + v_old[i][jm1][1])*
           (u_old[i][jp1][1] - u_old[i][jm1][1])*hy

      cc = (u_old[ip1][j][1] + u_old[im1][j][1])*
           (v_old[ip1][j][1] - v_old[im1][j][1])*hx

      dd = (v_old[i][jp1][1] + v_old[i][jm1][1])*
           (v_old[i][jp1][1] - v_old[i][jm1][1])*hy

      phix = (phi_old[ip1][j][1] - phi_old[im1][j][1])*hhx
      phiy = (phi_old[i][jp1][1] - phi_old[i][jm1][1])*hhy

      fu = f*u_old[i][j][1]
      fv = f*v_old[i][j][1]

      phiu = (phi_old[ip1][j][1]*u_old[ip1][j][1] -
              phi_old[im1][j][1]*u_old[im1][j][1])*hhx
      phiv = (phi_old[i][jp1][1]*v_old[i][jp1][1] -
              phi_old[i][jm1][1]*v_old[i][jm1][1])*hhy

      u_new[i][j] = u_old[i][j][kv] - dtv*(aa + bb + phix - fv)
      v_new[i][j] = v_old[i][j][kv] - dtv*(cc + dd + phiy + fu)
      phi_new[i][j] = phi_old[i][j][kv] - dtv*(phiu + phiv)
    }
  }

  setBoundaryConditions(u_new, v_new, phi_new, m, n)
}

function setBoundaryConditions(u, v, phi, m, n) {

  for (i=1; i < n-1; i++) {
    // left boundary, periodic
    u[  0][i] = u[  m-2][i]
    v[  0][i] = v[  m-2][i]
    phi[0][i] = phi[m-2][i]
  }

  for (i=0; i < m; i++) {
    // top boundary, wall
    u[  i][n-1] = u[  i][n-2]
    v[  i][n-1] = 0.
    phi[i][n-1] = phi[i][n-2]
  }

  for (i = 1; i < n-1; i++) {
    // right boundary, periodic
    u[  m-1][i] = u[  1][i]
    v[  m-1][i] = v[  1][i]
    phi[m-1][i] = phi[1][i]
  }

  for (i=0; i < m; i++) {
    // bottom boundary, wall
    u[i][0]   = u[i][1]
    v[i][0]   = 0.
    phi[i][0] = phi[i][1]
  }

  // top-right corner
  u[m-1][n-1] = (u[m-2][n-1] + u[m-1][n-2] + u[m-2][n-2])/3.
  v[m-1][n-1] = (v[m-2][n-1] + v[m-1][n-2] + v[m-2][n-2])/3.
  phi[m-1][n-1] = (phi[m-2][n-1] + phi[m-1][n-2] + phi[m-2][n-2])/3.

  // bottom-right corner
  u[m-1][0] = (u[m-2][0] + u[m-1][1] + u[m-2][1])/3.
  v[m-1][0] = (v[m-2][0] + v[m-1][1] + v[m-2][1])/3.
  phi[m-1][0] = (phi[m-2][0] + phi[m-1][1] + phi[m-2][1])/3.

  // bottom-left corner
  u[0][0] = (u[0][1] + u[1][0] + u[1][1])/3.
  v[0][0] = (v[0][1] + v[1][0] + v[1][1])/3.
  phi[0][0] = (phi[0][1] + phi[1][0] + phi[1][1])/3.

  // top-left corner
  u[0][n-1] = (u[0][n-2] + u[1][n-1] + u[1][n-2])/3.
  v[0][n-1] = (v[0][n-2] + v[1][n-1] + v[1][n-2])/3.
  phi[0][n-1] = (phi[0][n-2] + phi[1][n-1] + phi[1][n-2])/3.
}

function updateValues() {
  if (u_old == null) {
    u_old = create3DArray(m, n, 2);
    v_old = create3DArray(m, n, 2);
    phi_old = create3DArray(m, n, 2);
    type = 0;
  }
  else {
    type = 1;
  }

  if (type != 0) {
    for (i=0; i < m; i++) {
      for (j=0; j < n; j++) {
        u_old  [i][j][0] = u_old  [i][j][1]
        v_old  [i][j][0] = v_old  [i][j][1]
        phi_old[i][j][0] = phi_old[i][j][1]
      }
    }
  }

  for (i=0; i < m; i++) {
    for (j=0; j < n; j++) {
      u_old  [i][j][1] = u_new  [i][j]
      v_old  [i][j][1] = v_new  [i][j]
      phi_old[i][j][1] = phi_new[i][j]
    }
  }
}

//drawing primitves
function drawRectT(ctx, x, y, w, h){
  ctx.fillRect(x*W+xT, y*H+yT, w*W, h*H);
}

function lineT(ctx, x0, y0, x1, y1){
  ctx.beginPath();
  moveToT(ctx, x0, y0);
  lineToT(ctx, x1, y1);
  ctx.stroke();
}

function moveToT(ctx, x, y){
  ctx.moveTo(x*W+xT, (1-y)*H+yT);
}

function lineToT(ctx, x, y){
  ctx.lineTo(x*W+xT, (1-y)*H+yT);
}

//helper functions
function create2DArray(m, n) {
  var a = new Array(m);

  for (i = 0; i < m; i++) {
      a[i] = new Array(n);
  }

  return a;
}

function create3DArray(m, n, o) {
  var a = new Array(m);

  for (i = 0; i < m; i++) {
    a[i] = new Array(n);
    for (j = 0; j < n; j++) {
      a[i][j] = new Array(o);
    }
  }

  return a;
}


// Hyperbolic cosine of x
function cosh(x) {
  return (Math.exp(x) + Math.exp(-x)) / 2;
}

// Hyperbolic tangent of x
function tanh(x) {
  return (Math.exp(x) - Math.exp(-x)) / (Math.exp(x) + Math.exp(-x));
}

// Convert r, g, b to hex color
function RGB2Color(color)
{
  return '#' + byte2Hex(255*color[0] % 256) + byte2Hex(255*color[1] % 256) + byte2Hex(255*color[2] % 256);
}

function byte2Hex(n)
{
  // from http://krazydad.com/tutorials/makecolors.php
  var nybHexString = "0123456789ABCDEF";
  return String(nybHexString.substr((n >> 4) & 0x0F,1)) + nybHexString.substr(n & 0x0F,1);
}

function lininterp(x, x0, x1, y0, y1) {
    delta = x - x0;
    return y0 + (delta*y1 - delta*y0)/(x1 - x0);
}