Simulation of hydrodynamics using processing

BenGuy

Hello everyone here is a little program to simulate the movement of a liquid surface using processing. It is based on the shallow water equations (Saint Venant). I am interested in any comments and ideas on how to optimize the code and accelerate running, so the space resolution can be improved (currently 50x50 elements) and still keep the appearance of real time simulation. Thanks

/**
 * Solution of Saint Venant equations
 * Runs with Processing 1.5.1
 *
 * Simulation of drop impacts on water surface
 * Right click to initialize, left clik to drop
 *
 * By Ben Guy - May 2014
 */

int i, j, n, nm, np, nmm, npp, ii, c;
float[] Ho = new float[2601]; float[] HUo = new float[2601]; float[] HVo = new float[2601]; float[] Zo = new float[2601]; 
float[] H1o = new float[2601]; float[] HU1o = new float[2601]; float[] HV1o = new float[2601];
float amax;
float H, Him, Hip, Hjm, Hjp, H1;
float HU, HUim, HUip, HUjm, HUjp, HU1;
float HV, HVim, HVip, HVjm, HVjp, HV1;
float Z, Zim, Zip, Zjm, Zjp;
float U, Uim, Uip, Ujm, Ujp;
float V, Vim, Vip, Vjm, Vjp;
float Cip2g, Cim2d, Cjp2g, Cjm2d;
float Hip2g, Hip2d, Him2g, Him2d;
float Hjp2g, Hjp2d, Hjm2g, Hjm2d;
float FHip2g, FHim2d, FHjp2g, FHjm2d;
float HUip2g, HUip2d, HUim2g, HUim2d;
float HUjp2g, HUjp2d, HUjm2g, HUjm2d;
float HVip2g, HVip2d, HVim2g, HVim2d;
float HVjp2g, HVjp2d, HVjm2g, HVjm2d;
float FHUip2g, FHUim2d, FHUjp2g, FHUjm2d, Si;
float FHVip2g, FHVim2d, FHVjp2g, FHVjm2d, Sj;
//initialisation des constantes et plages
float g = 9.81; float dtx = 0.1; float dty = 0.1; float r = 0;
int k=10; int kz=100;
int ki, kip, kj, kjp;
float kzn; float kz1; float kz51; float kz52; 

void setup() {
  size(850, 650, P3D);
  noStroke();
}

void mousePressed() {
//initialisation de H, HU, HV, Z
  if (mouseButton == RIGHT) {
    for (n=0; n<2600; n+=1) {Ho[n]=2; HUo[n]=0; HVo[n]=0; Zo[n]=0;}
//obstacle
    for (j=20; j<30; j+=1) {
      for (i=24; i<26; i+=1) {
        n = i + 51*j;
        Zo[n]=2.5; Ho[n]=0;
      }
    }
  } else {
//goutte
    i=min(49, mouseX*50/850); j=min(49, mouseY*50/650);
    n = i + 51*j; Ho[n]=0.2; Ho[n+1]=0.2; Ho[n+51]=0.2; Ho[n+52]=0.2;
  }
}

void draw() {
  for (j=1; j<50; j+=1) {
    for (i=1; i<50; i+=1) {
      n = i + 51*j; nm=n-1; np=n+1; nmm=n-51; npp=n+51;
//U vitesse positive vers ligne+, V vitesse positive vers colonne+
//initialisation de H, HU, HV, Z, pour (i,j), (i+1,j), (i-1,j), (i,j+1), (i,j-1)
      H = Ho[n]; Him = Ho[nm]; Hip = Ho[np]; Hjm = Ho[nmm]; Hjp = Ho[npp];
      HU = HUo[n]; HUim = HUo[nm]; HUip = HUo[np]; HUjm = HUo[nmm]; HUjp = HUo[npp];
      HV = HVo[n]; HVim = HVo[nm]; HVip = HVo[np]; HVjm = HVo[nmm]; HVjp = HVo[npp];
      if (H > 0) {U = HU / H; V = HV / H;} else {V = 0; U = 0;}
      Z = Zo[n]; Zim = Zo[nm]; Zip = Zo[np]; Zjm = Zo[nmm]; Zjp = Zo[npp];
//variables en i+1/2 et i-1/2
//  hauteurs et vitesses
      if (Him > 0) {Uim = HUim / Him; Vim = HVim / Him;} else {Uim = 0; Vim = 0;}
      if (Hip > 0) {Uip = HUip / Hip; Vip = HVip / Hip;} else {Uip = 0; Vip = 0;}
      amax=max(Z, Zip); Hip2g=max(0, H+Z-amax); Hip2d=max(0, Hip+Zip-amax);
      amax=max(Zim, Z); Him2g=max(0, Him+Zim-amax); Him2d=max(0, H+Z-amax);
      HUip2g = Hip2g * U; HUip2d = Hip2d * Uip;
      HUim2g = Him2g * Uim; HUim2d = Him2d * U;
      HVip2g = Hip2g * V; HVip2d = Hip2d * Vip;
      HVim2g = Him2g * Vim; HVim2d = Him2d * V;
//  termes de flux
      Cip2g = max(abs(U) + sqrt(g * Hip2g), abs(Uip) + sqrt(g * Hip2d));
      Cim2d = max(abs(Uim) + sqrt(g * Him2g), abs(U) + sqrt(g * Him2d));
      FHip2g = (HUip2d + HUip2g - Cip2g * (Hip2d - Hip2g)) / 2;
      FHim2d = (HUim2d + HUim2g - Cim2d * (Him2d - Him2g)) / 2;
      FHUip2g = (HUip2d * Uip + g * Hip2d * Hip2d / 2 + HUip2g * U + g * Hip2g * Hip2g / 2) / 2 - Cip2g * (HUip2d - HUip2g) / 2;
      FHUim2d = (HUim2d * U + g * Him2d * Him2d / 2 + HUim2g * Uim + g * Him2g * Him2g / 2) / 2 - Cim2d * (HUim2d - HUim2g) / 2;
      FHVip2g = (HVip2d * Uip + HVip2g * U) / 2 - Cip2g * (HVip2d - HVip2g) / 2;
      FHVim2d = (HVim2d * U + HVim2g * Uim) / 2 - Cim2d * (HVim2d - HVim2g) / 2;
      Si = g * (Him2d * Him2d - Hip2g * Hip2g) / 2;
//variables en j+1/2 et j-1/2
//  hauteurs et vitesses
      if (Hjm > 0) {Ujm = HUjm / Hjm; Vjm = HVjm / Hjm;} else {Ujm = 0; Vjm = 0;}
      if (Hjp > 0) {Ujp = HUjp / Hjp; Vjp = HVjp / Hjp;} else {Ujp = 0; Vjp = 0;}
      amax=max(Z, Zjp); Hjp2g=max(0, H+Z-amax); Hjp2d=max(0, Hjp+Zjp-amax);
      amax=max(Zjm, Z); Hjm2g=max(0, Hjm+Zjm-amax); Hjm2d=max(0, H+Z-amax);
      HUjp2g = Hjp2g * U; HUjp2d = Hjp2d * Ujp;
      HUjm2g = Hjm2g * Ujm; HUjm2d = Hjm2d * U;
      HVjp2g = Hjp2g * V; HVjp2d = Hjp2d * Vjp;
      HVjm2g = Hjm2g * Vjm; HVjm2d = Hjm2d * V;
//  termes de flux
      Cjp2g = max(abs(V) + sqrt(g * Hjp2g), abs(Vjp) + sqrt(g * Hjp2d));
      Cjm2d = max(abs(Vjm) + sqrt(g * Hjm2g), abs(V) + sqrt(g * Hjm2d));
      FHjp2g = (HVjp2d + HVjp2g - Cjp2g * (Hjp2d - Hjp2g)) / 2;
      FHjm2d = (HVjm2d + HVjm2g - Cjm2d * (Hjm2d - Hjm2g)) / 2;
      FHUjp2g = (HUjp2d * Vjp + HUjp2g * V) / 2 - Cjp2g * (HUjp2d - HUjp2g) / 2;
      FHUjm2d = (HUjm2d * V + HUjm2g * Vjm) / 2 - Cjm2d * (HUjm2d - HUjm2g) / 2;
      FHVjp2g = (HVjp2d * Vjp + g * Hjp2d * Hjp2d / 2 + HVjp2g * V + g * Hjp2g * Hjp2g / 2) / 2 - Cjp2g * (HVjp2d - HVjp2g) / 2;
      FHVjm2d = (HVjm2d * V + g * Hjm2d * Hjm2d / 2 + HVjm2g * Vjm + g * Hjm2g * Hjm2g / 2) / 2 - Cjm2d * (HVjm2d - HVjm2g) / 2;
      Sj = g * (Hjm2d * Hjm2d - Hjp2g * Hjp2g) / 2;
//calcul des nouvelles valeurs
      H1 = H - dtx * (FHip2g - FHim2d) - dty * (FHjp2g - FHjm2d);
      HU1 = HU - dtx * (FHUip2g - FHUim2d + Si) - dty * (FHUjp2g - FHUjm2d);
      HV1 = HV - dtx * (FHVip2g - FHVim2d) - dty * (FHVjp2g - FHVjm2d + Sj);
//mémorisation des nouvelles valeurs de H, HU, HV
      H1o[n] = H1; HU1o[n] = HU1; HV1o[n] = HV1;
    }
//calcul et mémorisation des nouvelles valeurs aux limites en i=0 et i=50
    H1o[51*j] = H1o[1+51*j]; H1o[50+51*j] = H1o[49+51*j];
    HU1o[51*j] = r * HU1o[1+51*j]; HU1o[50+51*j] = r * HU1o[49+51*j]; //réflexion si r=0, libre si r=1 - Débit libre
    HV1o[51*j] = HV1o[1+51*j]; HV1o[50+51*j] = HV1o[49+51*j];
  }
//calcul et mémorisation des nouvelles valeurs aux limites en j=0 et j=50
  for (i=0; i<51; i+=1) {
    H1o[i] = H1o[i+51]; H1o[i+51*50] = H1o[i+51*49];
    HU1o[i] = HU1o[i+51]; HU1o[i+51*50] = HU1o[i+51*49];
    HV1o[i] = r * HV1o[i+51]; HV1o[i+51*50] = r * HV1o[i+51*49];      //réflexion si r=0, libre si r=1 - Débit libre
//    HV1o[i] = HVo[i+51]; HV1o[i+51*50] = HV1o[i+51*49];               //Débit entrant en haut fixé, libre en bas
  }
//écrasement des anciennes valeurs de H, HU, HV, et affichage
  for (n=0; n<2601; n+=1) {
    Ho[n] = H1o[n]; HUo[n] = HU1o[n]; HVo[n] = HV1o[n];
  }

  translate(80, 350, 0); rotateX(PI/3.5); rotateZ(-PI/5);
  background(51);
  directionalLight(255, 255, 255, 1, -1, -0.5);
  for (j=1; j<49; j+=1) {
    for (i=1; i<49; i+=1) {
      n = i + 51*j; ki=k*i; kip=k*(i+1); kj=k*j; kjp=k*(j+1);
      kzn=kz*(H1o[n]+Zo[n]); kz1=kz*(H1o[n+1]+Zo[n+1]); kz52=kz*(H1o[n+52]+Zo[n+52]); kz51=kz*(H1o[n+51]+Zo[n+51]);
      beginShape();
        vertex(ki, kj, kzn);
        vertex(kip, kj, kz1);
        vertex(kip, kjp, kz52);
        vertex(ki, kjp, kz51);
      endShape(CLOSE);
    }
  }
}