// pde_cylinderical_eq.c  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

#include <stdio.h>
#include <math.h>
#include <time.h>
#include "nuderiv.h"
#include "simeq.h"

#undef  abs
#define abs(a) ((a)<0.0?(-(a)):(a))
#undef  min
#define min(a,b) ((a)<(b)?(a):(b))
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))

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

  // 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;
    int i, j, k, ii, jj, kk;

    // 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
    printf("building discrete system of %d equations\n", Nrtz);
    Cs = Nrtz; // RHS
    // initialize A and X
    for(i=1; i<Nr-1; i++)
    {
      for(j=1; j<Nt-1; j++)
      {
        for(k=1; k<Nz-1; k++)
	{
          Row = Sd(i, j, k); // equation index
          for(ii=1; ii<Nr-1; ii++)
	  {
            for(jj=1; jj<Nt-1; jj++)
	    {
              for(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(i=1; i<Nr-1; i++)
    {
      R = Rg[i];
      for(j=1; j<Nt-1; j++)
      {
        T = Tg[j];
        for(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
          nuderiv(1, Nr, i, Rg, Cr);
          nuderiv(2, Nr, i, Rg, Crr);
          for(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
          nuderiv(2, Nt, j, Tg, Ctt);
          for(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
          nuderiv(2, Nz, k, Zg, Czz);
          for(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 = abs(Ud2ck);
          if(Err> 0.001)
	  {
	    printf("[%1d][%1d][%1d] Ud2ck=%e, Err=%e\n",
                   i, j, k, Ud2ck, Err);
	  }
        } // k
        if(2==2)  // debug print
	{
          Err1 = 0.0;
          Err2 = 0.0;
          Err3 = 0.0;
          for(ii=1; ii<Nr-1; ii++)
	  {
            for(jj=1; jj<Nt-1; jj++)
	    {
              for(kk=1; kk<Nz-1; kk++)
	      {
                printf("[%1d][%1d][%1d] A[][]=%e, %e, %e\n", ii, jj, kk,
		       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
          printf("RHS Cs=%d, A[][Cs]=%e, %e, %e\n",
		 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];
          printf("row check\n");
          printf("Err= %e, %e, %e\n", Err1, Err2, Err3);
        } // debug print
      } // j
    } // i
  } // end Build_A

int main(int argc, char * argv[]) // 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;
  int i, j, k, ii, jj, kk;  
  double time_cpu_start, time_cpu_end, time_wall_start, time_wall_end;

  printf("pde_cylinderical_eq.c starting\n");
  time_cpu_start = (double)clock()/(double)CLOCKS_PER_SEC;
  time_wall_start = (double)time(NULL);

  printf("r gird values, not necessarily equally spaced\n");
  for(i=0; i<Nr; i++)
  {
    Rg[i] = Rmin + Hr*(double)(i);
    printf("Rg[%1d]=%f\n", i, Rg[i]);
  }
  printf("theta gird values, not equally spaced\n");
  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(j=0; j<Nt; j++)
  {
    printf("Tg[%1d]=%f\n", j, Tg[j]);
  }
  printf("z gird values, not necessarily equally spaced\n");
  for(k=0; k<Nz; k++)
  {
    Zg[k] = Zmin+ Hz*(double)(k);
    printf("Tg[%1d]=%f\n", k, Zg[k]);
  }

  printf("analytic partial second derivative\n");
  for(i=0; i<Nr; i++)
  {
    R = Rg[i];
    for(j=0; j<Nt; j++)
    {
      T = Tg[j];
      for(k=0; k<Nz; k++)
      {
        Z = Zg[k];
        Ua[S(i,j,k)] = Ub(R,T,Z);
        printf("[%1d][%1d][%1d] \n", i, j, k);
        printf("Ua=%f, D2=%f\n", Ua[S(i,j,k)],
                D2(Ur(R,T,Z),Urr(R,T,Z),Utt(R,T,Z),Uzz(R,T,Z),R)); 
        printf("Urr=%f, Utt=%f, Uzz=%f \n", Urr(R,T,Z),
                Utt(R,T,Z), Uzz(R,T,Z));
	printf("Uuu+Ub=%f, F=%f\n", Uuu(R,T,Z)+Ub(R,T,Z), F(R,T,Z));
        Xt = R*cos(T);
        Yt = R*sin(T);
        if(abs(sqrt(Xt*Xt+Yt*Yt)-R)>0.0001) printf("R error\n");

        Tmp = atan2(Yt,Xt);
        if(Tmp<0.0) Tmp = 2.0*Pi + Tmp;
        if(abs(Tmp-T)>0.0001) printf("Theta error\n");
      } // k
    } // j
  } // i

  for(i=0; i<Nr; i++)
  {
    R = Rg[i];
    for(j=0; j<Nt; j++)
    {
      T = Tg[j];
      for(k=0; k<Nz; k++)
      {
        Z = Zg[k];
        // compute Ur, Urr, Utt, Uzz and call D2
        nuderiv(1, Nr, i, Rg, Cr);
        nuderiv(2, Nr, i, Rg, Crr);
        Urv = 0.0;
        Urrv = 0.0;
        for(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;
        nuderiv(2, Nt, j, Tg, Ctt);
        for(jj=0; jj<Nt; jj++)
        {
          Uttv = Uttv + Ctt[jj]*Ub(R,Tg[jj],Z);
        }
        Uzzv = 0.0;
        nuderiv(2, Nz, k, Zg, Czz);
        for(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 = max(Maxerr, abs(Err));
        printf("[%1d][%1d][%1d] \n", i, j, k);
        printf("Ana=%f, Cmp=%f, Err=%e\n", Ana, Cmp, Err);
        Err = Ur(R,T,Z)-Urv;
        Maxerrur = max(Maxerrur, abs(Err));
        printf("Ur=%f, Urv=%f, Err=%e \n", Ur(R,T,Z), Urv, Err);
        Err = Urr(R,T,Z)-Urrv;
        Maxerrurr = max(Maxerrurr, abs(Err));
        printf("Urr=%f, Urrv=%f, Err=%e\n", Urr(R,T,Z), Urrv, Err);
        Err = Utt(R,T,Z)-Uttv;
        Maxerrutt = max(Maxerrutt, abs(Err));
        printf("Utt=%f, Uttv=%f, Err=%e\n", Utt(R,T,Z), (Uttv), Err);
        Err = Uzz(R,T,Z)-Uzzv;
        Maxerruzz = max(Maxerruzz, abs(Err));
        printf("Uzz=%f, Uzzv=%f, Err=%e\n", Uzz(R,T,Z), Uzzv, Err);
        printf("\n");
      } // k
    } // j
  } // i
  printf("Maxerr Ur =%e\n", Maxerrur);
  printf("Maxerr Urr=%e\n", Maxerrurr);
  printf("Maxerr Utt=%e\n", Maxerrutt);
  printf("Maxerr Uzz=%e\n", Maxerruzz);
  printf("Maxerr deriv=%e\n", Maxerr);

  build_A();
  simeq2b(Nrtz, A, X); // X is solution
  printf("check against analytic solution\n");
  Maxerr = 0.0;
  for(i=1; i<Nr-1; i++)
  {
    for(j=1; j<Nt-1; j++)
    {
      for(k=1; k<Nz-1; k++)
      {
        Err = X[Sd(i, j, k)]-Ua[S(i,j,k)];
        Maxerr = max(Maxerr, abs(Err));
        printf("[%1d][%1d][%1d] X=%f, Ua=%f, Err=%e\n",
               i, j, k, X[Sd(i, j, k)], Ua[S(i,j,k)], Err);
      } // k
    } // j
  } // i
  time_cpu_end = (double)clock()/(double)CLOCKS_PER_SEC;
  time_wall_end = (double)time(NULL);
  printf("wall time=%f seconds, cpu time=%f seconds\n",
         time_wall_end-time_wall_start, 
         time_cpu_end -time_cpu_start);
  printf("Maxerr soln=%e\n", Maxerr);
  printf("pde_cylinderical_eq.c finishing\n");
  return 0;
} // pde_cylinderical_eq.c