-- corner0c.adb    uniform grid boundary value PDE
--                 known solution for testing method
--                 alpha = 1/2 on first quadrant
--                 r = sqrt(x^2+y^2)
--                 theta = arctan(y,x) 0..2Pi
--                 ub(x,y) = r^alpha * sin(alpha*theta) solution
--                 differential equation to solve
--                 d^2u/dx^2 + d^2u/dy^2 = 0
--
-- uniform grid on rectangle 0.01 to 1, 0.01 to 1
-- set up and solve system of linear equations

with Ada.Text_IO;
use  Ada.Text_IO;
with Ada.Calendar;
use  Ada.Calendar;
with Ada.Numerics;
use  Ada.Numerics;
with Ada.Numerics.Long_Elementary_Functions;
use  Ada.Numerics.Long_Elementary_Functions;
with Real_Arrays;
use  Real_Arrays;
with Nuderiv;

procedure corner0c is
  nx : constant :=11; -- including boundary at 1 and nx, 1..nx
  ny : constant :=11; -- including boundary at 1 and ny, 1..ny
  xmin : real := 0.01;
  xmax : real := 1.0;
  ymin : real := 0.01;
  ymax : real := 1.0;
  hx : real := (xmax-xmin)/real(nx-1);
  hy : real := (ymax-ymin)/real(ny-1);
  Start, Now : Duration;

  neqn : constant := (nx-2)*(ny-2);
  cs : integer := neqn+1;
  xg : real_vector(1..nx);
  yg : real_vector(1..ny);
  -- nx-2 by ny-2  internal, non boundary cells
  cxx : real_vector(1..nx);
  cyy : real_vector(1..ny);
  us : real_vector(1..neqn); -- solution vector
  ut : real_matrix(1..neqn, 1..neqn+1); -- includes RHS

  package fltio is new float_io(real);
  use fltio;
  package intio is new integer_io(integer);
  use intio;

  function Atan2(Y : Real; X : Real) return Real is
    Theta : Real; -- 0 .. 2Pi
  begin
    if abs(X)<1.0e-6 then
      if abs(Y)<1.0e-6 then
        return 0.0;
      elsif Y>=0.0 then
        Theta := Pi/2.0;
      else
        Theta := 3.0*Pi/2.0;
      end if;
    else
      Theta := arctan(Y,X);
    end if;
    if Theta<0.0 then
      Theta := 2.0*Pi+Theta;
    end if;
    return Theta;
  end Atan2;

  function U(X:real; Y:real) return real is
           -- solution for boundary and optional check
    Omega, Alpha, R, Theta : Real;
  begin
    Omega := 2.0*Pi;
    Alpha := Pi/Omega;
    R := sqrt(X*X+Y*Y);
    if R < 0.001 then return 0.0; end if;
    Theta := Atan2(Y,X);
    return R**Alpha * sin(Alpha*Theta);
  end U;

  function c(X:real; Y:real) return real is
     -- RHS of differential equation to solve
  begin
    return 0.0;
  end c;

  function S(i:integer; j:integer) return integer is
  begin
    return (i-2)*(ny-2)+(j-2)+1; -- index into us and ut
  end S;

  procedure Print_soln is
    error : real := 0.0;
    max_error : real := 0.0;
    avg_error : real := 0.0;
    val : real;
  begin
    new_line;
    put_line("exact solution u, computed us, error");
    for i in 1..nx loop
      for j in 1..ny loop
         put("xg("&integer'image(i)&")=");
         put(xg(i),3,3,0);
         put(", yg("&integer'image(j)&")=");
         put(yg(j),3,3,0);
         put(", u=");
         put(u(xg(i),yg(j)),3,5,0);
         if i=1 or i=nx or j=1 or j=ny then -- boundary
           val := u(xg(i),yg(j));
           error := 0.0;
         else
           val := us(S(i,j));
           error := abs(val-u(xg(i),yg(j)));
           avg_error := avg_error + error;
           if error > max_error then
              max_error := error;
           end if;
         end if;
         put(", us=");
         put(val,3,5,0);
         put(", err=");
         put_line(real'image(error));
      end loop;
    end loop;
    put_line("avg_error="&real'image(avg_error/real(neqn))&
             ", max_error="&real'image(max_error));
  end Print_soln;

  procedure init_matrix is
    ii, kj, ik : integer;
  begin
    -- cs is RHS, constant column subscript
    -- zero internal cells
    for i in 2..nx-1 loop
      for j in 2..ny-1 loop
        for ii in 2..nx-1 loop
          for jj in 2..ny-1 loop
            ut(S(i,j),S(ii,jj)) := 0.0;
          end loop;
        end loop;
        ut(S(i,j),cs) := 0.0;
      end loop;
    end loop;

    put_line("internal cells zeroed ");
    for i in 2..nx-1 loop -- inside boundary
      for j in 2..ny-1 loop
        ii := S(i,j); -- equation
          -- for each term

          -- d^2 u(x,y)/dx^2
          Nuderiv(2,nx,i,xg,cxx);
          for k in 1..nx loop
            if k=1 or k=nx then
              ut(ii,cs) := ut(ii,cs) - cxx(k)*u(xg(k),yg(j));
            else
              kj := S(k,j);
              ut(ii,kj) := ut(ii,kj) + cxx(k);
            end if;
          end loop; -- k
          -- d^2 u(x,y)/dy^2
          Nuderiv(2,ny,j,yg,cyy);
          for k in 1..ny loop
            if k=1 or k=ny then
              ut(ii,cs) := ut(ii,cs) - cyy(k)*u(xg(i),yg(k));
            else
              ik := S(i,k);
              ut(ii,ik) := ut(ii,ik) + cyy(k);
            end if;
          end loop; --k
          -- RHS constant
          ut(ii,cs) := ut(ii,cs) + c(xg(i),yg(j));
        -- end terms
      end loop; -- j
    end loop; -- i
  end init_matrix;

  procedure Print_B is
  begin
    put_line("matrix B includes RHS");
    for i in 2..nx-1 loop
      for j in 2..ny-1 loop
        for ii in 2..nx-1 loop
          for jj in 2..ny-1 loop
             put_line("i="&Integer'Image(i)&", j="&Integer'Image(j)&
                      ",  ii="&Integer'Image(ii)&", jj="&Integer'Image(jj)&
                       ", B(i,j)="&
                       Real'Image(ut(S(i,j),S(ii,jj))));
          end loop; -- jj
        end loop; -- ii
        put_line("i="&Integer'Image(i)&", j="&Integer'Image(j)&
                 ",                RHS="&
                       Real'Image(ut(S(i,j),cs)));
      end loop; -- j
    end loop; -- i
    New_line;
    New_line;
  end Print_B;

  procedure check_rows is -- only for checking
    -- row is in ut solution coordinates
    row, col : Integer;
    sum : Real;
  begin
    for i in 2..nx-1 loop
      for j in 2..ny-1 loop
        row := S(i,j);
        sum := 0.0;
        for ii in 2..nx-1 loop
          for jj in 2..ny-1 loop
            col := S(ii,jj);
            sum := sum + ut(row,col)*u(xg(ii),yg(jj)); -- u is check solution
          end loop; -- ii
        end loop; -- jj
        sum := sum - ut(row,cs);
        if abs(sum)>0.001 then
           put_line("i="&Integer'Image(i)&",j="&Integer'Image(j)&
                    " is bad err="&Real'Image(sum));
        end if;
      end loop; -- j
    end loop; --i
  end check_rows;

  procedure simeq(n : integer; b : in out real_matrix; x : out real_vector) is
    --  tailored version for this application, one based, in place
    row : array(1..n)of integer;   -- row interchange indices
    hold , i_pivot : integer;      -- pivot indices
    pivot : real;            -- pivot element value
    abs_pivot : real;        -- abs of pivot element
    norm1 : real := 1.0;     -- 1 norm of matrix
    matrix_data_error : exception;
  begin
    for k in 1..n loop
      row(k):= k;
    end loop;
    --                               begin main reduction loop
    for k in 1..n loop
      --                             find largest element for pivot
      pivot := b(row(k), k);
      abs_pivot := abs(pivot);
      i_pivot := k;
      for i in k..n loop
        if abs(b(row(i), k)) > abs_pivot then
          i_pivot := i;
          pivot := b(row(i), k);
          abs_pivot := abs(pivot);
        end if;
      end loop;
      --                             check for near singular
      if abs_pivot < real'epsilon * norm1 then
        raise matrix_data_error;
      end if;
      --                             have pivot, interchange row pointers
      hold := row(k);
      row(k):= row(i_pivot);
      row(i_pivot):= hold;
      --                           reduce about pivot
      for j in k+1..n+1 loop
        b(row(k), j) := b(row(k), j) / b(row(k), k);
      end loop;
      --                           inner reduction loop
      for i in 1..n loop
        if i /= k then
          for j in k+1..n+1 loop
            b(row(i), j) := b(row(i), j) - b(row(i), k) * b(row(k), j);
          end loop;
        end if;
      end loop;
      --                             finished inner reduction
    end loop;
    --                               build  x  for return, unscrambling rows
    for i in 1..n loop
      x(i):= b(row(i), n+1);
    end loop;
  end simeq;

  procedure check_soln is -- check when solution unknown
    smaxerr, err, x, y, uxx, uyy : Real;
  begin
    --  ux(x,y) + uy(x,y) = f(x,y)
    put_line("check_soln against PDE");
    smaxerr := 0.0;

    for i in 2..nx-1 loop -- through dof equations
      Nuderiv(2,nx,i,xg,cxx);
      x := xg(i);
      for ii in 2..ny-1 loop
        Nuderiv(2,ny,ii,yg,cyy);
        y := yg(ii);
        err := -c(x,y); -- add other side of equation terms
        uxx := 0.0;
        for j in 1..nx loop -- compute numerical derivative at (i,ii)
          if j=1 or j=nx then
            uxx := uxx + cxx(j)*u(xg(j),yg(ii)); -- boundary
          else
            uxx := uxx + cxx(j)*Us(S(j,ii));
          end if;
        end loop;
        uyy := 0.0;
        for jj in 1..ny loop
          if jj=1 or jj=ny then
            uyy := uyy + cyy(jj)*U(xg(i),yg(jj)); -- boundary
          else
            uyy := uyy + cyy(jj)*Us(S(i,jj));
          end if;
        end loop;
        err := err + uxx + uyy;
        if abs(err) > smaxerr then
          smaxerr := abs(err);
        end if;
      end loop; -- ii Y
    end loop; -- i X

    put_line("check soln against PDE max error="&Real'Image(smaxerr));
  end Check_Soln;

  -- write_soln2  for gnuplot
  procedure Write_Soln2(Ux : Real_Vector; File_Name : String) is
    Fout : File_Type;
  begin
    Create(Fout, Out_File, File_Name);

    Put_Line("writing "&File_name);
    for  i in 1..nx loop
      for ii in 1..ny loop
        Put(Fout," ");
        Put(Fout,xg(i),2,7,0);
        Put(Fout," ");
        Put(Fout,yg(ii),2,7,0);
        Put(Fout," ");
        if i=1 or i=nx or ii=1 or ii=ny then
          Put(Fout,U(xg(i),yg(ii)),4,7,0); -- boundary
        else
--          Put(Fout,U(xg(i),yg(ii)),4,7,0); -- check solution
          Put(Fout,Ux(S(i,ii)),4,7,0); -- computed solution
        end if;
        New_Line(Fout);
      end loop; -- ii
      New_Line(Fout);
    end loop; -- i
    Close(Fout);
    Put_Line("finished writing "&File_Name);
  end Write_Soln2;

  procedure CheckU(X : Real; Y : Real) is
  begin -- print a few boundary values
    Put("x="); Put(X,2,4,0);
    Put(", y="); Put(Y,2,4,0);
    Put(", th="); Put(Atan2(Y,X),2,4,0);
    put(", U="); put(U(X,Y),2,4,0);
    New_Line;
  end CheckU;

begin  -- pde21_eq

  put_line("corner0c.adb running");
  put_line("differential equation to solve");
  put_line("d^2u/dx^2 + d^2u/dy^2 = 0");
  put_line("uniform grid on rectangle -1,1 to -1,1");
  put_line("known Solution, for testing method");
  Put_Line("omega = 2 Pi");
  Put_Line("alpha = Pi/omega");
  Put_Line("theta = atan2(y,x)  0..2Pi");
  Put_Line("r = sqrt(x^2+y^2)");
  put_line("u(x,y) = r^alpha * sin(alpha*theta)");
  New_Line;

  Start := Seconds(Clock);
  for i in 1..nx loop
    xg(i) := xmin+real(i-1)*hx;
    put_line("xg("&integer'image(i)&")="&real'image(xg(i)));
  end loop;
  for j in 1..ny loop
    yg(j) := ymin+real(j-1)*hy;
    put_line("yg("&integer'image(j)&")="&real'image(yg(j)));
  end loop;
  put("xmin=");
  put(xmin,4,3,0);
  put(", xmax=");
  put(xmax,4,3,0);
  put(", hx=");
  put(hx,4,3,0);
  put_line(", nx="&integer'image(nx));
  put("ymin=");
  put(ymin,4,3,0);
  put(", ymax=");
  put(ymax,4,3,0);
  put(", hy=");
  put(hy,4,3,0);
  put_line(", ny="&integer'image(ny));
  CheckU( 1.0, 0.0);
  CheckU( 1.0, 1.0);
  CheckU( 0.0, 1.0);
  CheckU(-1.0, 1.0);
  CheckU(-1.0, 0.0);
  CheckU(-1.0,-1.0);
  CheckU( 0.0,-1.0);
  CheckU( 1.0,-1.0);
  CheckU( 0.0, 0.0);
  New_Line;

  init_matrix;
  put_line("matrix initialized");
  -- print_B;
  check_rows;
  simeq(neqn, ut, us); -- tailored version, returned us
  check_soln;
  Print_soln;
  Write_Soln2(us, "corner0c_ada.dat");
  Now := Seconds(Clock);
  put_line("finished corner0c.adb in "&Duration'Image(Now-Start)&" seconds");
end corner0c;