// pde_cylinderical_eq.java  cylinderical coordinates
// solve del^2 U(radius,theta,z) + U(radius,theta,z) = f(r,t,z)
// given Dirchlet boundary values Ub(r,t,z), and f(r,t,z)
// Much checking in this code, analytic solution is
// U(r,t,z) = r^3 + t^3/216 + z^3
//
// f(r,t,z) := 12*r + 1/36*t/(r^2) + 6*z
//             r^3 + 1/216*t^3 + z^3
//
// x = r cos theta    0 <= theta < 2 Pi
// y = r sin theta    r < 0  restriction for derivative
// z = z              z double
// r = sqrt(x^2+y^2+z^2)
// theta = arctan(y,x)  if theta<0  theta = 2 Pi + arctan(y,x)
// z = z
//
// use existing partial derivative computation in Cartesian coordinates to
// compute partial derivatives in cylinderical coordinates
//
// compute del U(radius,theta,z) r=radius, t=theta angle, z=height
// del U  =  (d/dr U, 1/r d/dt U, d/dz U)
//
// compute del^2 U(radius,theta,z) r=radius, t=theta angle, z=height
// del^2 U  =  d/dr(r d/dr U) + 1/r^2 d^2/dt^2 U + d^2/dz^2 U
//          =  2/r d/dr U + d^2/dr^2 U + 1/r^2 d^2/dt^2 U + d^2/dz^2 U

// needs nuderiv.java
// needs simeq.java

public class pde_cylinderical_eq
{
  final int Nr = 5; // number radius
  final int Nt = 7; // number theta angle
  final int Nz = 5; // number Z coordinates
  final int Nrtz = (Nr-2)*(Nt-2)*(Nz-2); // matrix, equations
  double A[][] = new double[Nrtz][Nrtz+1]; // no bound, has RHS
  double X[] = new double[Nrtz]; // computed solution to U
  double Rg[] = new double[Nr]; // r grid values
  double Tg[] = new double[Nt]; // theta grid values
  double Zg[] = new double[Nz]; // z grid values
  double Ua[] = new double[Nr*Nt*Nz]; // analytic solution for checking
  double Cr[] = new double[Nr]; // first deriv wrt r
  double Crr[] = new double[Nr]; // second deriv wrt r
  double Ctt[] = new double[Nt]; // second deriv wrt t
  double Czz[] = new double[Nz]; // second deriv wrt z

  pde_cylinderical_eq()
  {
  double Pi = 3.1415926535897932384626433832795028841971;
  double Rmin = 0.2;
  double Rmax = 1.0;
  double Hr = (Rmax-Rmin)/(double)(Nr-1);
  double Tmin = 0.1;
  double Tmax = 300.0*Pi/180.0;
  // Ht not applicable
  double Zmin = -1.0;
  double Zmax =  1.0;
  double Hz = (Zmax-Zmin)/(double)(Nz-1);
  double R, T, Z, Ana, Cmp, Err;
  double Urv, Urrv, Uttv, Uzzv; // computed derivatives
  double Maxerrur=0.0, Maxerrurr=0.0, Maxerrutt=0.0;
  double Maxerruzz=0.0, Maxerr=0.0;
  double Xt, Yt, Tmp;
  double time_wall_start, time_wall_end;

  System.out.println("pde_cylinderical_eq.java starting");
  time_wall_start = System.currentTimeMillis()/1000.0;

  System.out.println("r gird values, not necessarily equally spaced");
  for(int i=0; i<Nr; i++)
  {
    Rg[i] = Rmin + Hr*(double)(i);
    System.out.println("Rg["+i+"]= "+Rg[i]);
  }
  System.out.println("theta gird values, not equally spaced");
  Tg[0] = 0.1;
  Tg[1] = 40.0*Pi/180.0;
  Tg[2] = 85.0*Pi/180.0;
  Tg[3] = 130.0*Pi/180.0;
  Tg[4] = 180.0*Pi/180.0;
  Tg[5] = 240.0*Pi/180.0;
  Tg[6] = 300.0*Pi/180.0;
  for(int j=0; j<Nt; j++)
  {
    System.out.println("Tg["+j+"]= "+Tg[j]);
  }
  System.out.println("z gird values, not necessarily equally spaced");
  for(int k=0; k<Nz; k++)
  {
    Zg[k] = Zmin+ Hz*(double)(k);
    System.out.println("Tg["+k+"]= "+ Zg[k]);
  }

  System.out.println("analytic partial second derivative");
  for(int i=0; i<Nr; i++)
  {
    R = Rg[i];
    for(int j=0; j<Nt; j++)
    {
      T = Tg[j];
      for(int k=0; k<Nz; k++)
      {
        Z = Zg[k];
        Ua[S(i,j,k)] = Ub(R,T,Z);
        System.out.println("["+i+"]["+j+"]["+k+"]");
        System.out.println("Ua="+Ua[S(i,j,k)]+", D2="+
                D2(Ur(R,T,Z),Urr(R,T,Z),Utt(R,T,Z),Uzz(R,T,Z),R)); 
        System.out.println("Urr="+Urr(R,T,Z)+", Utt="+
                Utt(R,T,Z)+", Uzz="+Uzz(R,T,Z));
	System.out.println("Uuu+Ub="+(Uuu(R,T,Z)+Ub(R,T,Z))+
                ", F="+F(R,T,Z));
        Xt = R*Math.cos(T);
        Yt = R*Math.sin(T);
        if(Math.abs(Math.sqrt(Xt*Xt+Yt*Yt)-R)>0.0001) 
          System.out.println("R error");

        Tmp = Math.atan2(Yt,Xt);
        if(Tmp<0.0) Tmp = 2.0*Pi + Tmp;
        if(Math.abs(Tmp-T)>0.0001) System.out.println("Theta error");
      } // k
    } // j
  } // i

  for(int i=0; i<Nr; i++)
  {
    R = Rg[i];
    for(int j=0; j<Nt; j++)
    {
      T = Tg[j];
      for(int k=0; k<Nz; k++)
      {
        Z = Zg[k];
        // compute Ur, Urr, Utt, Uzz and call D2
        new nuderiv(1, Nr, i, Rg, Cr);
        new nuderiv(2, Nr, i, Rg, Crr);
        Urv = 0.0;
        Urrv = 0.0;
        for(int ii=0; ii<Nr; ii++)
	{
           Urv = Urv + Cr[ii]*Ub(Rg[ii],T,Z);
           Urrv = Urrv + Crr[ii]*Ub(Rg[ii],T,Z);
        }
        Uttv = 0.0;
        new nuderiv(2, Nt, j, Tg, Ctt);
        for(int jj=0; jj<Nt; jj++)
        {
          Uttv = Uttv + Ctt[jj]*Ub(R,Tg[jj],Z);
        }
        Uzzv = 0.0;
        new nuderiv(2, Nz, k, Zg, Czz);
        for(int kk=0; kk<Nz; kk++)
        {
          Uzzv = Uzzv + Czz[kk]*Ub(R,T,Zg[kk]);
        }

        Cmp = D2(Urv, Urrv, Uttv, Uzzv, R);
        Ana = Uuu(R,T,Z);
        Err = Ana - Cmp;
        Maxerr = Math.max(Maxerr, Math.abs(Err));
        System.out.println("["+i+"]["+j+"]["+k+"]");
        System.out.println("Ana="+Ana+", Cmp="+Cmp+", Err="+Err);
        Err = Ur(R,T,Z)-Urv;
        Maxerrur = Math.max(Maxerrur, Math.abs(Err));
        System.out.println("Ur="+Ur(R,T,Z)+", Urv="+Urv+", Err="+Err);
        Err = Urr(R,T,Z)-Urrv;
        Maxerrurr = Math.max(Maxerrurr, Math.abs(Err));
        System.out.println("Urr="+Urr(R,T,Z)+", Urrv="+Urrv+", Err="+Err);
        Err = Utt(R,T,Z)-Uttv;
        Maxerrutt = Math.max(Maxerrutt, Math.abs(Err));
        System.out.println("Utt="+Utt(R,T,Z)+", Uttv="+Uttv+", Err="+Err);
        Err = Uzz(R,T,Z)-Uzzv;
        Maxerruzz = Math.max(Maxerruzz, Math.abs(Err));
        System.out.println("Uzz="+Uzz(R,T,Z)+", Uzzv="+Uzzv+", Err="+Err);
        System.out.println("");
      } // k
    } // j
  } // i
  System.out.println("Maxerr Ur = "+ Maxerrur);
  System.out.println("Maxerr Urr= "+ Maxerrurr);
  System.out.println("Maxerr Utt= "+ Maxerrutt);
  System.out.println("Maxerr Uzz= "+ Maxerruzz);
  System.out.println("Maxerr deriv= "+ Maxerr);

  build_A();
  new simeq(Nrtz, A, X); // X is solution
  System.out.println("check against analytic solution");
  Maxerr = 0.0;
  for(int i=1; i<Nr-1; i++)
  {
    for(int j=1; j<Nt-1; j++)
    {
      for(int k=1; k<Nz-1; k++)
      {
        Err = X[Sd(i, j, k)]-Ua[S(i,j,k)];
        Maxerr = Math.max(Maxerr, Math.abs(Err));
        System.out.println("["+i+"]["+j+"]["+k+"] X="+
               X[Sd(i, j, k)]+", Ua="+Ua[S(i,j,k)]+", Err="+Err);
      } // k
    } // j
  } // i
  time_wall_end = System.currentTimeMillis()/1000.0;
  System.out.println("wall time="+(time_wall_end-time_wall_start)+
                     ",  seconds");
  System.out.println("Maxerr soln="+ Maxerr);
  System.out.println("pde_cylinderical_eq.java finishing");
  } // end pde_cylinderical_eq

  // Ub for computing boundaries and checking computed solution
  double Ub(double R, double T, double Z)
  {
    return R*R*R + T*T*T/216.0 + Z*Z*Z;
  } // end Ub

  // analytic derivatives for checking
  double Ur(double R, double T, double Z)
  {
    return 3.0*R*R;
  } // end Ur

  double Urr(double R, double T, double Z)
  {
    return 6.0*R;
  } // end Urr

  double Ut(double R, double T, double Z)
  {
    return 3.0*T*T/216.0;
  } // end Ut

  double Utt(double R, double T, double Z)
  {
    return 6.0*T/216.0;
  } // end Utt

  double Uz(double R, double T, double Z)
  {
    return 3.0*Z*Z;
  } // end Uz;

  double Uzz(double R, double T, double Z)
  {
    return 6.0*Z;
  } // end Uzz

  double Uuu(double R, double T, double Z)
  {
    return 2.0*Ur(R,T,Z)/R + Urr(R,T,Z) + Utt(R,T,Z)/(R*R) + Uzz(R,T,Z);
  } // end Uuu;

  double D2(double Urv, double Urrv, double Uttv, double Uzzv, double R)
  {
    return 2.0*Urv/R + Urrv + Uttv/(R*R)+ Uzzv;
  } // end D2

  int S(int i, int j, int k)
  { // index of each grid point
    return i*Nt*Nz + j*Nz + k;
  } // end S

  int Sd(int i, int j, int k)
  { // index of each grid point inside boundary, A and X
    return (i-1)*(Nt-2)*(Nz-2) + (j-1)*(Nz-2) + (k-1);
  } // end Sd

  double F(double R, double T, double Z)
  {
    return 12.0*R + (1.0/36.0)*T/(R*R) + 6.0*Z +
           R*R*R + (1.0/216.0)*T*T*T + Z*Z*Z;
  } // end F   RHS

  void build_A()    // A is discrete system of equations, including Y
  {                 // solution is X, from A*X=Y
    int Row, Idof, Cs;
    double Urrs, Utts, Uzzs, Ud2ck, Err;
    double Err1, Err2, Err3 ;
    double R, T, Z;

    // equation Urrs(r,t,z)+Utts(r,t,z)+Uzzs(r,t,z)+U(r,t,z) = f(r,t,z)
    // for every (r,t,z) inside boundary
    System.out.println("building discrete system of "+Nrtz+" equations");
    Cs = Nrtz; // RHS
    // initialize A and X
    for(int i=1; i<Nr-1; i++)
    {
      for(int j=1; j<Nt-1; j++)
      {
        for(int k=1; k<Nz-1; k++)
	{
          Row = Sd(i, j, k); // equation index
          for(int ii=1; ii<Nr-1; ii++)
	  {
            for(int jj=1; jj<Nt-1; jj++)
	    {
              for(int kk=1; kk<Nz-1; kk++)
	      {
                A[Row][Sd(ii,jj,kk)] = 0.0;
	      } // kk
            } // jj
          } // ii
          A[Row][Cs] = F(Rg[i],Tg[j],Zg[k]); // RHS
          X[Row] = 0.0;
        } // k
      } // j
    } // i
    // now build A inside boundary (could use sparse)
    for(int i=1; i<Nr-1; i++)
    {
      R = Rg[i];
      for(int j=1; j<Nt-1; j++)
      {
        T = Tg[j];
        for(int k=1; k<Nz-1; k++)
	{
          Z = Zg[k];
          Row = Sd(i, j, k); // equation index
          Ud2ck = 0.0; // only for development checking

          // equation Urrs(r,t,z)+Utts(r,t,z)+Uzzs(r,t,z)+U(r,t,z) = f(r,t,z)
          // Urrs term
          new nuderiv(1, Nr, i, Rg, Cr);
          new nuderiv(2, Nr, i, Rg, Crr);
          for(int ii=0; ii<Nr; ii++)
	  {
            Urrs = 2.0*Cr[ii]/R + Crr[ii];
            Ud2ck = Ud2ck + Urrs;

            if(ii==0 || ii==Nr-1) // boundary, negative on RHS
	    {
              A[Row][Cs] = A[Row][Cs] - Urrs*Ub(Rg[ii],T,Z);
	    }
            else
	    {
              Idof = Sd(ii,j,k);
              A[Row][Idof] = A[Row][Idof] + Urrs;
	    }
          } // ii

          // Utts term
          new nuderiv(2, Nt, j, Tg, Ctt);
          for(int jj=0; jj<Nt; jj++)
	  {
            Utts = Ctt[jj]/(R*R);
            Ud2ck = Ud2ck + Utts;
            if(jj==0 || jj==Nt-1) // boundary, negative on RHS
	    {
               A[Row][Cs] = A[Row][Cs] - Utts*Ub(R,Tg[jj],Z);
	    }
            else
	    {
              Idof = Sd(i,jj,k);
              A[Row][Idof] = A[Row][Idof] + Utts;
	    }
	  } // jj

          // Uzzs term
          new nuderiv(2, Nz, k, Zg, Czz);
          for(int kk=0; kk<Nz; kk++)
	  {
            Uzzs = Czz[kk];
            Ud2ck = Ud2ck + Uzzs;
            if(kk==0 || kk==Nz-1) // boundary, negative on RHS
	    {
               A[Row][Cs] = A[Row][Cs] - Uzzs*Ub(R,T,Zg[kk]);
	    }
            else
	    {
              Idof = Sd(i,j,kk);
              A[Row][Idof] = A[Row][Idof] + Uzzs;
            }
          } // kk

          // U term
          A[Row][Row] = A[Row][Row] + 1.0;

          // check
          Err = Math.abs(Ud2ck);
          if(Err> 0.001)
	  {
	    System.out.println("["+i+"]["+j+"]["+k+"] Ud2ck="+Ud2ck+
                               ", Err="+Err);
	  }
        } // k
        if(2==2)  // debug print
	{
          Err1 = 0.0;
          Err2 = 0.0;
          Err3 = 0.0;
          for(int ii=1; ii<Nr-1; ii++)
	  {
            for(int jj=1; jj<Nt-1; jj++)
	    {
              for(int kk=1; kk<Nz-1; kk++)
	      {
                System.out.println("["+ii+"]["+jj+"]["+kk+"] A[][]="+
		       A[0][Sd(ii,jj,kk)]+", "+
		       A[1][Sd(ii,jj,kk)]+", "+
                       A[2][Sd(ii,jj,kk)]);
                Err1 = Err1 + A[0][Sd(ii,jj,kk)]*Ub(Rg[ii],Tg[jj],Zg[kk]);
                Err2 = Err2 + A[1][Sd(ii,jj,kk)]*Ub(Rg[ii],Tg[jj],Zg[kk]);
                Err3 = Err3 + A[2][Sd(ii,jj,kk)]*Ub(Rg[ii],Tg[jj],Zg[kk]);
              } // kk
            } // jj
          } // ii
          System.out.println("RHS Cs="+Cs+", A[][Cs]="+
		             A[0][Cs]+", "+A[1][Cs]+", "+A[2][Cs]);
          Err1 = Err1 - A[0][Cs];
          Err2 = Err2 - A[1][Cs];
          Err3 = Err3 - A[2][Cs];
          System.out.println("row check");
          System.out.println("Err= "+Err1+", "+Err2+", "+Err3);
        } // debug print
      } // j
    } // i
  } // end Build_A

  public static void main (String[] args)
  {
     new pde_cylinderical_eq();
  }
} // end class and pde_cylinderical_eq.java