-- pde_cylinderical_eq.adb  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 real
-- 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

with Ada.Text_IO;
use  Ada.Text_IO;
with Ada.Numerics.Long_Elementary_Functions;
use  Ada.Numerics.Long_Elementary_Functions;
with Real_Arrays;
use  Real_Arrays;
with Nuderiv; -- procedure
with Simeqb; -- procedure
procedure pde_cylinderical_eq is

  package Real_IO is new Ada.Text_IO.Float_IO(Real);
  use Real_IO;

  Pi : Real := 3.1415926535897932384626433832795028841971;
  Nr : Integer := 5;
  Rmin : Real := 0.2;
  Rmax : Real := 1.0;
  Hr : Real := (Rmax-Rmin)/Real(Nr-1);
  Nt : Integer := 7;
  Tmin : Real := 0.1;
  Tmax : Real := 300.0*Pi/180.0;
  -- Ht not applicable
  Nz : Integer := 5;
  Zmin : Real := -1.0;
  Zmax : Real :=  1.0;
  Hz : Real := (Zmax-Zmin)/Real(Nz-1);
  Rg : Real_Vector(1..Nr); -- r grid values
  Tg : Real_Vector(1..Nt); -- theta grid values
  Zg : Real_Vector(1..Nz); -- z grid values
  Ua : Real_Vector(1..Nr*Nt*Nz);
  Cr : Real_Vector(1..Nr); -- first deriv wrt r
  Crr : Real_Vector(1..Nr); -- second deriv wrt r
  Ctt : Real_Vector(1..Nt); -- second deriv wrt t
  Czz : Real_Vector(1..Nz); -- second deriv wrt z
  R, T, Z, Ana, Cmp, Err : Real;
  Urv, Urrv, Uttv, Uzzv : Real; -- computed derivatives
  Maxerrur, Maxerrurr, Maxerrutt, Maxerruzz, Maxerr : Real := 0.0;
  Xt, Yt, Tmp : Real;
  Nrtz : Integer := (Nr-2)*(Nt-2)*(Nz-2);
  A : Real_Matrix(1..Nrtz,1..Nrtz+1); -- no bound, has RHS
  X : Real_Vector(1..Nrtz); -- computer solution to U
  
  function Ub(R:Real; T:Real; Z:Real) return Real is
  begin
    return R**3 + T**3/216.0 + Z**3;
  end Ub;

  function Ur(R:Real; T:Real; Z:Real) return Real is
  begin
    return 3.0*R**2;
  end Ur;

  function Urr(R:Real; T:Real; Z:Real) return Real is
  begin
    return 6.0*R;
  end Urr;

  function Ut(R:Real; T:Real; Z:Real) return Real is
  begin
    return 3.0*T**2/216.0;
  end Ut;

  function Utt(R:Real; T:Real; Z:Real) return Real is
  begin
    return 6.0*T/216.0;
  end Utt;

  function Uz(R:Real; T:Real; Z:Real) return Real is
  begin
    return 3.0*Z**2;
  end Uz;

  function Uzz(R:Real; T:Real; Z:Real) return Real is
  begin
    return 6.0*Z;
  end Uzz;

  function Uuu(R:Real; T:Real; Z:Real) return Real is
  begin
    return 2.0*Ur(R,T,Z)/R + Urr(R,T,Z) + Utt(R,T,Z)/(R*R) + Uzz(R,T,Z);
  end Uuu;

  function D2(Urv:Real; Urrv:Real; Uttv:Real; Uzzv:Real; R:Real) return Real is
  begin
    return 2.0*Urv/R + Urrv + Uttv/(R*R)+ Uzzv;
  end D2;

  function S(I:Integer; J:Integer; K:Integer) return Integer is
  begin -- index of each grid point
    return (I-1)*Nt*Nz + (J-1)*Nz + K;
  end S;

  function Sd(I:Integer; J:Integer; K:Integer) return Integer is
  begin -- index of each grid point inside boundary, A and X
    return (I-2)*(Nt-2)*(Nz-2) + (J-2)*(Nz-2) + (K-1);
  end Sd;

  function F(R:Real; T:Real; Z:Real) return Real is -- RHS
  begin
    return 12.0*R + (1.0/36.0)*T/(R**2) + 6.0*Z +
           R*R*R + (1.0/216.0)*T*T*T + Z*Z*Z;
  end F;  -- RHS

  procedure build_A is -- A is discrete system of equations, including Y
                    -- solution is X, from A*X=Y
    Row, Idof, Cs : Integer;
    Urrs, Utts, Uzzs, Ud2ck, Err: Real;
    Err1, Err2, Err3 : Real;
  begin
    -- 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
    Put_Line("building discrete system of "&Integer'Image(Nrtz)&" equations");
    Cs := Nrtz+1; -- RHS
    -- initialize A and X
    for I in 2..Nr-1 loop
      for J in 2..Nt-1 loop
        for K in 2..Nz-1 loop
          Row := Sd(I, J, K); -- equation index
          for Ii in 2..Nr-1 loop
            for Jj in 2..Nt-1 loop
              for Kk in 2..Nz-1 loop
                A(Row,Sd(Ii,Jj,Kk)) := 0.0;
              end loop; -- Kk
            end loop; -- Jj
          end loop; -- Ii
          A(Row,Cs) := F(Rg(I),Tg(J),Zg(K)); -- RHS
          X(Row) := 0.0;
        end loop; -- K
      end loop; -- J
    end loop; -- I
    -- now build A inside boundary (could use sparse)
    for I in 2..Nr-1 loop
      R := Rg(I);
      for J in 2..Nt-1 loop
        T := Tg(J);
        for K in 2..Nz-1 loop
          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 in 1..Nr loop
            Urrs := 2.0*Cr(Ii)/R + Crr(Ii);
            Ud2ck := Ud2ck + Urrs;

            if Ii=1 or Ii=Nr then -- 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;
            end if;
          end loop; --Ii

          -- Utts term
          nuderiv(2, Nt, J, Tg, Ctt);
          for Jj in 1..Nt loop
            Utts := Ctt(Jj)/(R*R);
            Ud2ck := Ud2ck + Utts;
            if Jj=1 or Jj=Nt then -- 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;
            end if;
          end loop; -- Jj

          -- Uzzs term
          nuderiv(2, Nz, K, Zg, Czz);
          for Kk in 1..Nz loop
            Uzzs := Czz(Kk);
            Ud2ck := Ud2ck + Uzzs;
            if Kk=1 or Kk=Nz then -- 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;
            end if;
          end loop; -- Kk

          -- U term
          A(Row,Row) := A(Row,Row) + 1.0;

          -- check
          Err := abs(Ud2ck);
          if Err> 0.001 then
             Put_line("("&Integer'Image(I)&","&Integer'Image(J)&","&
                        Integer'Image(K)&") Ud2ck="&
                        Real'Image(Ud2ck)&", err="&Real'Image(Err));
          end if;
        end loop; -- K
        if Row=3 then -- debug print
          Err1 := 0.0;
          Err2 := 0.0;
          Err3 := 0.0;
          for Ii in 2..Nr-1 loop
            for Jj in 2..Nt-1 loop
              for Kk in 2..Nz-1 loop
                Put_line("("&Integer'Image(Ii)&","&Integer'Image(Jj)&","&
                        Integer'Image(Kk)&") A(1)="&
                        Real'Image(A(1,Sd(Ii,Jj,Kk)))&", A(2)="&
                        Real'Image(A(2,Sd(Ii,Jj,Kk)))&", A(3)="&
                        Real'Image(A(3,Sd(Ii,Jj,Kk))));
                Err1 := Err1 + A(1,Sd(Ii,Jj,Kk))*Ub(Rg(Ii),Tg(Jj),Zg(Kk));
                Err2 := Err2 + A(2,Sd(Ii,Jj,Kk))*Ub(Rg(Ii),Tg(Jj),Zg(Kk));
                Err3 := Err3 + A(3,Sd(Ii,Jj,Kk))*Ub(Rg(Ii),Tg(Jj),Zg(Kk));
              end loop; -- Kk
            end loop; -- Jj
          end loop; -- Ii
          Put_line("RHS Cs="&Integer'Image(Cs)&" A(1)="&
                   Real'Image(A(1,Cs))&", A(2)="&
                   Real'Image(A(2,Cs))&", A(3)="&
                   Real'Image(A(3,Cs)));
          Err1 := Err1 - A(1,Cs);
          Err2 := Err2 - A(2,Cs);
          Err3 := Err3 - A(3,Cs);
          Put_Line("row check");
          Put_Line("Err1="&Real'Image(Err1)&
                   ",  Err2="&Real'Image(Err2)&
                   ",  Err3="&Real'Image(Err3));
        end if; -- debug print
      end loop; -- J
    end loop; -- I
  end Build_A;

begin
  Put_Line("pde_cylinderical_eq.adb starting");
  Put_Line("r gird values, not necessarily equally spaced");
  for I in 1..Nr loop
     Rg(I) := Rmin + Hr*Real(I-1);
     Put_Line("Rg("&Integer'Image(I)&")="&Real'Image(Rg(I)));
  end loop;
  Put_Line("theta gird values, not equally spaced");
  Tg(1) := 0.1;
  Tg(2) := 40.0*Pi/180.0;
  Tg(3) := 85.0*Pi/180.0;
  Tg(4) := 130.0*Pi/180.0;
  Tg(5) := 180.0*Pi/180.0;
  Tg(6) := 240.0*Pi/180.0;
  Tg(7) := 300.0*Pi/180.0;
  for J in 1..Nt loop
     Put_Line("Tg("&Integer'Image(J)&")="&Real'Image(Tg(J)));
  end loop;
  Put_Line("z gird values, not necessarily equally spaced");
  for K in 1..Nz loop
     Zg(K) := Zmin+ Hz*Real(K-1);
     Put_Line("Tg("&Integer'Image(K)&")="&Real'Image(Zg(K)));
  end loop;

  Put_Line("analytic partial second derivative");
  for I in 1..Nr loop
    R := Rg(I);
    for J in 1..Nt loop
      T := Tg(J);
      for K in 1..Nz loop
      Z := Zg(K);
        Ua(S(I,J,K)) := Ub(R,T,Z);
        Put_Line("Ua("&Integer'Image(I)&
                  ","&Integer'Image(J)&
                  ","&Integer'Image(K)&
                  ")="&Real'Image(Ua(S(I,J,K)))&
                   ", D2="&Real'Image(D2(Ur(R,T,Z),Urr(R,T,Z),Utt(R,T,Z),
				      Uzz(R,T,Z),R)));
	Put_Line("Uuu+U="&Real'Image(Uuu(R,T,Z)+Ub(R,T,Z))&
		 ", F="&Real'Image(F(R,T,Z)));
        Xt := R*Cos(T);
        Yt := R*Sin(T);
        if abs(Sqrt(Xt**2+Yt**2)-R)>0.0001 then
	   Put_Line("R error");
        end if;
        Tmp := Arctan(Yt,Xt);
        if Tmp<0.0 then Tmp := 2.0*Pi + Tmp; end if;
        if abs(Tmp-T)>0.0001 then
	   Put_Line("Theta error");
        end if;	
      end loop;
    end loop;
  end loop;

  for I in 1..Nr loop
    R := Rg(I);
    for J in 1..Nt loop
      T := Tg(J);
      for K in 1..Nz loop
        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 in 1..Nr loop
           Urv := Urv + Cr(Ii)*Ub(Rg(Ii),T,Z);
           Urrv := Urrv + Crr(Ii)*Ub(Rg(Ii),T,Z);
        end loop;
        Uttv := 0.0;
        Nuderiv(2, Nt, J, Tg, Ctt);
        for Jj in 1..Nt loop
           Uttv := Uttv + Ctt(Jj)*Ub(R,Tg(Jj),Z);
        end loop;
        Uzzv := 0.0;
        Nuderiv(2, Nz, K, Zg, Czz);
        for Kk in 1..Nz loop
           Uzzv := Uzzv + Czz(Kk)*Ub(R,T,Zg(Kk));
        end loop;

        Cmp := D2(Urv, Urrv, Uttv, Uzzv, R);
        Ana := Uuu(R,T,Z);
        Err := Ana - Cmp;
        Maxerr := Real'Max(Maxerr, Err);
        Put_Line("i="&Integer'Image(I)&
                  ", j="&Integer'Image(J)&
                  ", k="&Integer'Image(K)&
                  ", Ana="&Real'Image(Ana)&
                  ", Cmp="&Real'Image(Cmp)&
                  ", err="&Real'Image(Err));
        Err := Ur(R,T,Z)-Urv;
        Maxerrur := Real'Max(Maxerrur, Err);
        Put_Line("i="&Integer'Image(I)&
                  ", j="&Integer'Image(J)&
                  ", k="&Integer'Image(K)&
                  ", Ur="&Real'Image(Ur(R,T,Z))&
                  ", Urv="&Real'Image(Urv)&
                  ", err="&Real'Image(Err));
        Err := Urr(R,T,Z)-Urrv;
        Maxerrurr := Real'Max(Maxerrurr, Err);
        Put_Line("i="&Integer'Image(I)&
                  ", j="&Integer'Image(J)&
                  ", k="&Integer'Image(K)&
                  ", Urr="&Real'Image(Urr(R,T,Z))&
                  ", Urrv="&Real'Image(Urrv)&
                  ", err="&Real'Image(Err));
        Err := Utt(R,T,Z)-Uttv;
        Maxerrutt := Real'Max(Maxerrutt, Err);
        Put_Line("i="&Integer'Image(I)&
                  ", j="&Integer'Image(J)&
                  ", k="&Integer'Image(K)&
                  ", Utt="&Real'Image(Utt(R,T,Z))&
                  ", Uttv="&Real'Image(Uttv)&
                  ", err="&Real'Image(Err));
        Err := Uzz(R,T,Z)-Uzzv;
        Maxerruzz := Real'Max(Maxerruzz, Err);
        Put_Line("i="&Integer'Image(I)&
                  ", j="&Integer'Image(J)&
                  ", k="&Integer'Image(K)&
                  ", Utt="&Real'Image(Uzz(R,T,Z))&
                  ", Uttv="&Real'Image(Uzzv)&
                  ", err="&Real'Image(Err));
        Put_Line(" ");
      end loop;
    end loop;
  end loop;
  Put_Line("Maxerr Ur="&Real'Image(Maxerrur));
  Put_Line("Maxerr Urr="&Real'Image(Maxerrurr));
  Put_Line("Maxerr Utt="&Real'Image(Maxerrutt));
  Put_Line("Maxerr Uzz="&Real'Image(Maxerruzz));
  Put_Line("Maxerr deriv="&Real'Image(Maxerr));

  build_A;
  simeqb(Nrtz, A, X); -- X is solution
  Put_Line("check against analytic solution");
  Maxerr := 0.0;
  for I in 2..Nr-1 loop
    for J in 2..Nt-1 loop
      for K in 2..Nz-1 loop
        Err := X(Sd(I, J, K))-Ua(S(I,J,K));
        Maxerr := Real'Max(Maxerr, Abs(Err));
        Put_line("("&Integer'Image(I)&","&
                   Integer'Image(J)&","&
                   Integer'Image(K)&") "&
                   Real'Image(X(Sd(I, J, K)))&", ana="&
                   Real'Image(Ua(S(I,J,K)))&", err="&
                   Real'Image(Err));
      end loop; -- K
    end loop; -- J
  end loop; -- I
  Put_Line("Maxerr soln="&Real'Image(Maxerr));

  Put_line("pde_cylinderical_eq.adb finishing");
end pde_cylinderical_eq;