liquidsph3d.java

JAVA program for Liquid SPH 3D models

    double sgm;
	
    sgm = 0.25;
    i=0;

    for (ix=0; ix<2; ix++) {
      for (iy=0; iy<2; iy++) {
        for (iz=0; iz<2; iz++) {
          wkx[i] = -0.5*sgm+(xMax+sgm)*(double)(ix);
          wky[i] = -0.5*sgm+(yMax+sgm)*(double)(iy);
          wkz[i] = -0.5*sgm+(zMax+sgm)*(double)(iz);
          i += 1;
        }
      }
    }
  }


/*-------------------------- set initial condition --------- */

  void setInitialCondition() {
    int i,j,ix,iy,iz;
    double s;

    t = 0.0;

    i = 0;
    for (ix=0; ix<ncube; ix++) {
      for (iy=0; iy<ncube; iy++) {
        for (iz=0; iz<ncube; iz++) {
          kind[i] = 0;
          mass[i] = 1.0;
          xx[i] = 0.4*(ix)+0.5;
          yy[i] = 0.4*(iy)+1.5;
          zz[i] = 0.4*(iz)+0.5;
          vx[i] = 0.5*(Math.random()-0.5);
          vy[i] = 0.5*(Math.random()-0.5);
          vz[i] = 0.5*(Math.random()-0.5);
          ax[i] = 0.0;
          ay[i] = 0.0;
          az[i] = 0.0;
          hh[i] = hh0;
          i += 1;
        }
      }
    }
/*
    for (ix=0; ix<10; ix++) {
      for (iy=0; iy<10; iy++) {
        for (iz=0; iz<1; iz++) {
          kind[i] = 1;
          mass[i] = 1.0;
          xx[i] = 2.0+0.866*rBond*ix;
          yy[i] = 2.0+rBond*iy+0.5*rBond*(ix%2);
          zz[i] = 7.0;
          vx[i] = 0.0;
          vy[i] = 0.0;
          vz[i] = 0.0;
          ax[i] = 0.0;
          ay[i] = 0.0;
          az[i] = 0.0;
          hh[i] = hh0;
          i += 1;
        }
      }
    }
*/

    NNp = i;
    bonding();

    for (i=0; i<NNp; i++) {
      energy[i] = 1.0;
      power[i] = 0.0;
    }

    dividingSection();
    setDensity();
    dens0 = density[ncube*ncube*ncube/2+ncube*ncube/2+ncube/2]*0.85;

    rotate(theta,fai);
    zSortInit();
  }

/* --------------------------- bonding manager ---------- */

  void bonding() {
    int i,j;
    double r2,rb,rb2;

    rb = rBond*1.2;
    rb2 = rb*rb;

    bondInit();

    for (i=0; i<NNp; i++) {
      if (kind[i]!=0) {
        for (j=i+1; j<NNp; j++) {
          if (kind[j]!=0) {
            r2 = rSq(i,j);
            if (r2<=rb2) {
              bondReg(i, j); bondReg(j, i);
            }
          }
        }
      }
    }
  }

  void bondInit() {
    int i,j;

    for (i=0; i<NNp; i++) {
      for (j=0; j<20; j++) {
        bond[i][j]=-1;
      }
    }
  }

  void bondReg(int ii, int jj) {
    int j;

    for (j=0; j<20; j++) {
      if (bond[ii][j]==-1) {
        bond[ii][j] = jj;
        break;
      }
    }
  }

  void bondUnreg(int ii, int jj) {
    int j;

    for (j=0; j<20; j++) {
      if (bond[ii][j]==jj) {
        bond[ii][j] = -1;
        break;
      }
    }
  }

  int bondQ(int ii, int jj) {
    int j,q;
    double r2,rb,rb2;

    rb = rBond*1.2;
    rb2 = rb*rb;

    q = 0;
    for (j=0; j<20; j++) {
      if (bond[ii][j]==jj) {
        q = 1;
        r2 = rSq(ii,jj);
        if (r2>rb2) {
          bondReg(ii,jj); bondReg(jj,ii);
          q = 0;
        }
        break;
      }
    }
    return (q);
  }

  double rSq(int i, int j) {
    double xij,yij,zij;

    xij = xx[i]-xx[j]; yij = yy[i]-yy[j]; zij = zz[i]-zz[j];
    return(xij*xij+yij*yij+zij*zij);
  }

/* ----------------------------- timeEvolution -------------- */

  void timeEvolution() {
    int i;

    if (resetSW==1) {
      setInitialCondition();
	  resetSW = 0;
    }
    if (started==1) {
      for (i=0; i<incRate; i++) {
        timeStep();
      }
    }
  }

/*-------------------------------  particle motion   ---------*/

  void timeStep() {
    int i,j;
	double dtv2,s,rr;

    t = t + dt;
    dtv2 = dt/2.0;

    dividingSection();

    setDensity();

    for (i=0; i<NNp; i++) {
      pressure[i] = 10.0*(Math.pow(density[i]/dens0,7.0)-1);
    }

    for (i=0; i<NNp; i++) {
	  if (kind[i]>=0) {
        vx[i] += dtv2*ax[i];
        vy[i] += dtv2*ay[i];
		vz[i] += dtv2*az[i];
        xx[i] += vx[i]*dt;
        yy[i] += vy[i]*dt;
		zz[i] += vz[i]*dt;
	  }
    }
    changeCalc();
    for (i=0; i<NNp; i++) {
      if (kind[i]>=0) {
        vx[i] += dtv2*ax[i];
        vy[i] += dtv2*ay[i];
		vz[i] += dtv2*az[i];
        energy[i] += power[i]*dt;
	  }
    }

    rr = 1.0;
    for (i=0; i<NNp; i++) {
      if (xx[i] < 0.0) { 
        xx[i] = 0.0; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i];  vz[i] = rr*vz[i]; 
      }
      if (xx[i] > xMax) {
        xx[i] = xMax; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = rr*vz[i];  
      }
      if (yy[i] < 0.0) { 
        yy[i] = 0.0; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; vz[i] = rr*vz[i]; 
      }
      if (yy[i] > yMax) {
        yy[i] = yMax; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; vz[i] = rr*vz[i];
      }
      if (zz[i] < 0.0) { 
        zz[i] = 0.0; vx[i] = rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = -rr*vz[i]; 
      }
      if (zz[i] > zMax) {
        zz[i] = zMax; vx[i] = rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = -rr*vz[i];
      }
    }
  }

  void changeCalc() {
    int i,j,k,i0,i1,j0,j1,k0,k1,ii,jj,kk,iq;
    double ai,aj,b,rij,gradW,gradWx,gradWy,gradWz,a;
    double rcc,rcc2,r2,vr,mu,pij;

    rcc = hh0*2.0;
    rcc2 = rcc*rcc;

    for (i=0; i<NNp; i++) {
	  ax[i] = 0.0; ay[i] = gravity; az[i] = 0.0;
      boundaryForce(i, rcc, 100.0);
	  power[i] = 0.0;
	  ai = pressure[i]/(density[i]*density[i]);

	  i0 = (int)(Nsx*(xx[i]-rcc)/xMax); if (i0<0) i0 = 0;
      i1 = (int)(Nsx*(xx[i]+rcc)/xMax); if (i1>=Nsx) i1 = Nsx-1;
      j0 = (int)(Nsy*(yy[i]-rcc)/yMax); if (j0<0) j0 = 0;
      j1 = (int)(Nsy*(yy[i]+rcc)/yMax); if (j1>=Nsy) j1 = Nsy-1;
      k0 = (int)(Nsz*(zz[i]-rcc)/zMax); if (k0<0) k0 = 0;
      k1 = (int)(Nsz*(zz[i]+rcc)/zMax); if (k1>=Nsz) k1 = Nsy-1;
      for(ii=i0;ii<=i1;ii++) {
        for(jj=j0;jj<=j1;jj++) {
		  for(kk=k0;kk<=k1;kk++) {
            for(iq=1;iq<=section[ii][jj][kk][0];iq++) {
              j = section[ii][jj][kk][iq];
			  r2=(xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j])+(zz[i]-zz[j])*(zz[i]-zz[j]);
			  if (r2<rcc2) {
		        rij = Math.sqrt(r2);
		        if (rij>0) {
                  aj = pressure[j]/(density[j]*density[j]);
				  vr = (xx[i]-xx[j])*(vx[i]-vx[j])+(yy[i]-yy[j])*(vy[i]-vy[j])+(zz[i]-zz[j])*(vz[i]-vz[j]);
				  pij = 0;
				  if (vr<0.0) {
				    mu = 0.5*vr/(r2+0.001*0.5*0.5);
				    pij = (-alpha*0.5*mu+0.0*mu*mu)/((density[i]+density[j])/2.0);
				  }
		          b = mass[j]*(ai+aj+pij);
		          gradW = dwwvdr(rij, (hh[i]+hh[j])/2.0);
		          gradWx = gradW*(xx[i]-xx[j])/rij;
		          gradWy = gradW*(yy[i]-yy[j])/rij;
		          gradWz = gradW*(zz[i]-zz[j])/rij;
		          ax[i] += -b*gradWx;
                  ay[i] += -b*gradWy;
			      az[i] += -b*gradWz;
		          power[i] += 0.5*b*((vx[i]-vx[j])*gradWx+(vy[i]-vy[j])*gradWy+(vz[i]-vz[j])*gradWz);
				  if (kind[j]!=0 || kind[i]!=0) {
				    a = force(rij,i,j)/mass[i];
				    ax[i] += a*(xx[i]-xx[j])/rij;
				    ay[i] += a*(yy[i]-yy[j])/rij;
				    az[i] += a*(zz[i]-zz[j])/rij;
				  }
				}
			  }
			}
		  }
		}
	  }
    }
  }

  double force(double r, int ii, int jj) {
    int i,j;
    double ret,rc,ra,ra1,ra2,ram,fa,f,r1,r3,kc;


    if (kind[ii]==0 || kind[jj]==0) {

      rc = hh0*1.2;
      ret = 0.0;
      if (r<rc) {
	    ra = (rc-r)/rc;
	    ret = 1000.0*ra*ra;
		return ( ret );
	  }

	}

    if( kind[ii]>0 && kind[jj]>0 ) {

      kc = 2000.0;
      ra = rBond;
      ra1 = 1.1*ra;
      ra2 = 1.2*ra;
      ram = (ra1+ra2)/2.0;
      fa = 2.0/((ra2-ra1)*(ra2-ra1));

      if (ii<jj) {
        i= ii; j=jj;
      } else {
        i= jj; j=ii;
      }
      if (r>ra2) {
        bondUnreg(i,j); bondUnreg(j,i);
      } 

      f = 0.0;
      if (bondQ(i,j)==0) {
        if (r<ra) {
          f = kc*(r-ra);
        }
      } else {
        if (r<=ra1) {
          f = kc*(r-ra);
        } else if (r>=ra2) {
          f = 0.0;
        } else if (ra1<=r && r<ram) {
          f = kc*(r-ra)*(1.0-fa*(r-ra1)*(r-ra1));
        } else if (ram<=r && r<ra2) {
          f = kc*(r-ra)*fa*(r-ra2)*(r-ra2);
        }
      }
      return ( -f );

	}

    return (0.0);
  }

  void boundaryForce(int i, double rc, double k) {
    double r;

    if (xx[i] < rc) {
	  r = rc-xx[i];
      ax[i] += k*r*r; 
    }
    if (xx[i] > xMax-rc) {
	  r = xx[i]-(xMax-rc);
      ax[i] += -k*r*r; 
    }
    if (yy[i] < rc) {
	  r = rc-yy[i];
      ay[i] += k*r*r; 
    }
    if (yy[i] > yMax-rc) {
	  r = yy[i]-(yMax-rc);
      ay[i] += -k*r*r; 
    }
    if (zz[i] < rc) {
	  r = rc-zz[i];
      az[i] += k*r*r; 
    }
    if (zz[i] > zMax-rc) {
	  r = zz[i]-(zMax-rc);
      az[i] += -k*r*r; 
    }
  }


  void setDensity() {
    int i,j,k,i0,i1,j0,j1,k0,k1,ii,jj,kk,iq;
    double s,rij;
    double rcc,rcc2,r2;

    rcc = hh0*2.0;
    rcc2 = rcc*rcc;

    for (i=0; i<NNp; i++) {
	  s = 0.0;
	  i0 = (int)(Nsx*(xx[i]-rcc)/xMax); if (i0<0) i0 = 0;
      i1 = (int)(Nsx*(xx[i]+rcc)/xMax); if (i1>=Nsx) i1 = Nsx-1;
      j0 = (int)(Nsy*(yy[i]-rcc)/yMax); if (j0<0) j0 = 0;
      j1 = (int)(Nsy*(yy[i]+rcc)/yMax); if (j1>=Nsy) j1 = Nsy-1;
      k0 = (int)(Nsz*(zz[i]-rcc)/zMax); if (k0<0) k0 = 0;
      k1 = (int)(Nsz*(zz[i]+rcc)/zMax); if (k1>=Nsz) k1 = Nsz-1;
      for(ii=i0;ii<=i1;ii++) {
        for(jj=j0;jj<=j1;jj++) {
		  for(kk=k0;kk<=k1;kk++) {
            for(iq=1;iq<=section[ii][jj][kk][0];iq++) {
              j = section[ii][jj][kk][iq];
			  r2=(xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j])+(zz[i]-zz[j])*(zz[i]-zz[j]);
			  if (r2<rcc2) {
		        rij = Math.sqrt(r2);
                s += mass[j]*kernel(rij, (hh[i]+hh[j])/2.0);
			  }
			}
		  }
		}
	  }
	  density[i] = s;
    }
  }

/*-------------------------------  registration   -------------*/

  void dividingSection() {
    int i,j,k,ip,iq;

    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        for(k=0;k<Nsz;k++) {
          section[i][j][k][0] = 0;
        }
      }
    }

    for(ip=0;ip<NNp;ip++) {
      i = (int)(Nsx*xx[ip]/xMax); if (i>=Nsx) i = Nsx-1;
      j = (int)(Nsy*yy[ip]/yMax); if (j>=Nsy) j = Nsy-1;
      k = (int)(Nsz*zz[ip]/zMax); if (k>=Nsz) k = Nsz-1;
      iq = section[i][j][k][0]+1;
      section[i][j][k][0] = iq;
      section[i][j][k][iq] = ip;
    }
  }

  int maxSection() {
    int i,j,k,m;

    m = 0;
    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        for(k=0;k<Nsz;k++) {
          if (section[i][j][k][0]>m) m = section[i][j][k][0];
		}
      }
    }
    return ( m );
  }

/*----------------------------  smoothed particle ---------*/

  double kernel(double r, double h) {
    double a,q,ret;

    ret = 0.0;
	q = r/h;
	a = 1.0/3.141592/(h*h*h);
	if (q<1.0) {
      ret = a*(1.0-1.5*q*q+0.75*q*q*q);
    } else if (q<2.0) {
      ret = a*0.25*(2.0-q)*(2.0-q)*(2.0-q);
	}
    return ( ret );
  }

  double dwwvdr(double r, double h) {
    double a,q,ret;

    ret = 0.0;
	q = r/h;
	a = 1.0/3.141592/(h*h*h*h);
	if (q<1.0) {
      ret = a*(-3.0*q+2.25*q*q);
    } else if (q<=2.0) {
      ret = -a*0.75*(2.0-q)*(2.0-q);
	}
    return ( ret );
  }

/* ----------------------------- end of applet -------------- */

}