// nav1.java
// Navier Stokes test data.
// Velocity vector  u,v,w  in  x,y,z direction. Pressure p
// Space is x,y,z and time is t.
// u(x,y,z,t) is velocity in x direction at that point.
// ux(x,y,z,t) is derivative of u(x,y,z,t) wrt x.
// uyy(x,y,z,t) is second derivative of u(x,y,z,t) wrt y. etc.
// gx(x,y,z,t) is x direction forcing function, computed here for testing.
// test of funcions at end.

class nav1
{
  double mu = 1.4;
  double rho = 1.2;

  nav1()
  {
    System.out.println("nav1 instantiated");
  }
  double  u(double x, double y, double z, double t)
  {
    return 1.0 - 0.1*x + 0.1*y*y + 0.2*z*t;
  }
  double  v(double x, double y, double z, double t)
  {
    return 1.0 + 0.1*z*z - 0.1*y - 0.2*x*t;
  }
  double w(double x, double y, double z, double t)
  {
    return 1.2 + 0.1*y*y - 0.2*x*t + 0.2*z;
  }
  double ux(double x, double y, double z, double t)
  {
    return -0.1;
  }
  double vy(double x, double y, double z, double t)
  {
    return  -0.1;
  }
  double wz(double x, double y, double z, double t)
  {
    return 0.2;
  }
  double continuity(double x, double y, double z, double t)
  { // should return zero for all x,y,z,t
    return ux(x,y,z,t)+vy(x,y,z,t)+wz(x,y,z,t);
  }


  double uy(double x, double y, double z, double t)
  {
    return 0.2*y;
  }
  double uz(double x, double y, double z, double t)
  {
    return 0.2*t;
  }
  double ut(double x, double y, double z, double t)
  {
    return 0.2*z;
  }
  double vx(double x, double y, double z, double t)
  {
    return -0.2*t;
  }
  double vz(double x, double y, double z, double t)
  {
    return 0.2*z;
  }
  double vt(double x, double y, double z, double t)
  {
    return  -0.2*x;
  }
  double wx(double x, double y, double z, double t)
  {
    return -0.2*t;
  }
  double wy(double x, double y, double z, double t)
  {
    return 0.2*y;
  }
  double wt(double x, double y, double z, double t)
  {
    return -0.2*x;
  }
  double uxx(double x, double y, double z, double t)
  {
    return 0.0;
  }
  double uyy(double x, double y, double z, double t)
  {
    return 0.2;
  }
  double uzz(double x, double y, double z, double t)
  {
    return 0.2;
  }
  double vxx(double x, double y, double z, double t)
  {
    return 0.0;
  }
  double vyy(double x, double y, double z, double t)
  {
    return 0.0;
  }
  double vzz(double x, double y, double z, double t)
  {
    return 0.2;
  }
  double wxx(double x, double y, double z, double t)
  {
      return 0.0;
  }
  double wyy(double x, double y, double z, double t) 
  {
    return 0.2;
  }
  double wzz(double x, double y, double z, double t) 
  {
    return 0.0;
  }

  // Now pressure and derivatives.
  double p(double x, double y, double z, double t)
  {
    return 0.5 + 0.15*x*x  - 0.2*y + 0.25*z - 0.05*t;
  }
  double px(double x, double y, double z, double t)
  {
    return 0.3*x;
  }
  double py(double x, double y, double z, double t)
  {
    return -0.2;
  }
  double pz(double x, double y, double z, double t)
  {
    return 0.25; 
  }
  double gx(double x, double y, double z, double t)
  {
    return (ut(x,y,z,t)+u(x,y,z,t)*ux(x,y,z,t)+v(x,y,z,t)*uy(x,y,z,t)+
           w(x,y,z,t)*uz(x,y,z,t))+(1.0/rho)*px(x,y,z,t)-
           (mu/rho)*(uxx(x,y,z,t)+uyy(x,y,z,t)+uzz(x,y,z,t));
 
  }
  double gy(double x, double y, double z, double t)
  {
    return (vt(x,y,z,t)+u(x,y,z,t)*vx(x,y,z,t)+v(x,y,z,t)*vy(x,y,z,t)+
           w(x,y,z,t)*vz(x,y,z,t))+(1.0/rho)*py(x,y,z,t)-
           (mu/rho)*(vxx(x,y,z,t)+vyy(x,y,z,t)+vzz(x,y,z,t));
  }
  double gz(double x, double y, double z, double t)
  {
    return (wt(x,y,z,t)+u(x,y,z,t)*wx(x,y,z,t)+v(x,y,z,t)*wy(x,y,z,t)+
           w(x,y,z,t)*wz(x,y,z,t))+(1.0/rho)*pz(x,y,z,t)-
           (mu/rho)*(wxx(x,y,z,t)+wyy(x,y,z,t)+wzz(x,y,z,t));

  }
  double gxf(double x, double y, double z, double t)
  {
    return 0.2*z - 0.5666666667 + 0.26*x - 0.03*y*y  + 0.02*z*t + 0.2*y
         + 0.02*y*z*z  - 0.04*y*x*t + 0.24*t + 0.02*t*y*y  - 0.04*x*t*t;
  }
  double gyf(double x, double y, double z, double t)
  {
    return -0.2*x - 0.2*t + 0.04*x*t - 0.02*t*y*y  - 0.04*z*t*t
	   - 0.5 + 0.03*z*z  + 0.01*y + 0.24*z + 0.02*z*y*y  - 0.04*z*x*t;
  }
  double gzf(double x, double y, double z, double t)
  {
    return - 0.2*x - 0.2*t - 0.02*x*t - 0.02*t*y*y  - 0.04*z*t*t
           + 0.2*y + 0.02*y*z*z  - 0.04*y*x*t + 0.215 + 0.04*z;
  }
  void test()
  {
    double xmin = 0.0;
    double xmax = 1.0;
    double ymin = 0.0;
    double ymax = 1.0;
    double zmin = 0.0;
    double zmax = 1.0;
    double tmin = 0.0;
    double tmax = 1.0;
    double x, y, z, t;
    double hx, hy, hz, ht;
    int nx = 5;
    int ny = 5;
    int nz = 5;
    int nt = 5;

    hx = (xmax-xmin)/(nx-1);
    hy = (ymax-ymin)/(ny-1);
    hz = (zmax-zmin)/(nz-1);
    ht = (tmax-tmin)/(nt-1);
    for(int i=0; i<nx; i++)
    {
      x = xmin+hx*i;
      for(int j=0; j<nx; j++)
      {
        y = ymin+hy*j;
        for(int k=0; k<nx; k++)
        {
          z = zmin+hz*k;
          for(int m=0; m<nx; m++)
          {
            t = tmin+ht*m;
            // various equality tests
            if(Math.abs(continuity(x, y, z, t))>1.0e-6)
	    {
	      System.out.println("continuity error ="+continuity(x, y, z, t)
				 +" at x="+x+", y="+y+", z="+z+", t="+t);
	    }

            if(Math.abs(gxf(x,y,z,t)-gx(x,y,z,t))>1.0e-6)
	    {
		System.out.println("gxf error ="+(gxf(x,y,z,t)-gx(x,y,z,t))
				 +" at x="+x+", y="+y+", z="+z+", t="+t);
	    }

	  } // end t
	} // end z
      } // end y
    } // end x
  } // end test
} // end class nav1