-- pde_nl13a.adb  solve non linear second order, third degree PDE
--                automatic generation of nonlinear system of equations
--                later, merge boundaries so not in equations
--
-- solve  u^2 u'' + 2 u u' + 3u = f(x)
--        f(x) = (41/3)-(2116/9)x+(4780/3)x^2-(48976/9)x^3+
--               (48976.0/9.0)*x*x*x + 9984.0*x*x*x*x -
--               (88064.0/9.0)*x*x*x*x*x + (14336.0/3.0)*x*x*x*x*x*x -
--               9984x^4-(88064/9)x^5+(14336/3)x^6-(8192/9)x^7
-- with boundary  u(1)=1, u(2)=1  0 < x < 1
-- set up system of equations, use specialized Jacobian, simeq_newton5
-- solution u(x) = 1.0 -(20.0/3.0)*x +12.0*x*x -(16.0/3.0)*x*x*x
--
--

with Ada.Text_IO;
use  Ada.Text_IO;
with Ada.Numerics.Long_Elementary_Functions;
use  Ada.Numerics.Long_Elementary_Functions;
with Real_Arrays; -- types real, real_vector, real_matrix
use  Real_Arrays;
with Integer_Arrays; -- types integer_vector
use  Integer_Arrays;
with deriv;
with rderiv;
with Simeq_Newton5;

procedure pde_nl13a is
  nx : Integer := 5;    -- number of coefficients for discrete derivative
  neqn : Integer := nx; -- number of equations and number of linear unknowns
  nmax : Integer := (nx-2)*(nx-2)*(nx-2); -- number of nonlinear terms
                                          -- more than enough for this PDE
  nlin : Integer := 0;                    -- computed actual nonlinear terms
  ntot : Integer := neqn+nmax; -- maximum total number of terms
  k  : Real_Matrix(1..neqn,1..ntot); -- nonlinear system of equations
  xg : Real_Vector(1..nx); -- x grid
  f  : Real_Vector(1..neqn); -- RHS of equations
  u  : Real_Vector(1..ntot); -- solution being computed
  Ua : Real_Vector(1..neqn);
  cx  : Real_Vector(1..nx);    -- numerical first derivative
  cxx : Real_Vector(1..Nx);   -- numerical second derivative
  hx : real;
  err, avgerr, Maxerr, x : real;
  eps : real := 1.0e-6;
  xmin : Real := 0.0;
  xmax : Real := 1.0;
  j : Integer;

  -- linear, then nonlinear terms
  Var1 : Integer_vector(1..Ntot); -- to be generated
         -- first neqn are 1, 2, ..., neqn
  Var2 : Integer_Vector(1..Ntot); -- zero if not second order term
  Var3 : Integer_vector(1..Ntot); -- zero if not third order term
  Vari1 : Integer_vector(1..Ntot); -- zero if not third order term
  Vari2 : Integer_vector(1..Ntot); -- zero if not third order term

  -- definition of PDE to be solved
  function FF(x:real) return real is -- forcing function
  begin
    return (41.0/3.0)-(2116.0/9.0)*x +(4780.0/3.0)*x*x -
           (48976.0/9.0)*x*x*x + 9984.0*x*x*x*x -
           (88064.0/9.0)*x*x*x*x*x + (14336.0/3.0)*x*x*x*x*x*x -
           (8192.0/9.0)*x*x*x*x*x*x*x;
  end FF;

  function uana(x:real) return real is -- analytic solution and boundary
  begin
    return 1.0 -(20.0/3.0)*x +12.0*x*x -(16.0/3.0)*x*x*x;
  end uana;

  procedure check_soln(U : Real_Vector) is -- check when solution unknown
    -- u^2 u'' + 2 u u' + 3u = f(x) check for each xg inside boundary
    xi, err, Smaxerr : real;
  begin
    Put_line("check_soln");
    smaxerr := 0.0;
    for I in 2..neqn-1 loop
      xi := xg(i);
      rderiv(1, nx, i, hx, cx);
      rderiv(2, nx, i, hx, cxx);
      err := 3.0*U(i)-FF(xi);
      for J in 1..nx loop
        err := err + U(i)*U(i)*cxx(j)*U(j);
        err := err + 2.0*U(i)*cx(j)*U(j);
      end loop; -- j
      if abs(err)>smaxerr then Smaxerr := abs(err); end if;
      Put_line("i="&Integer'Image(I)&", err="&Real'Image(err));
    end loop; -- i
    Put_line("check_soln max error="&Real'Image(smaxerr));
  end Check_Soln;

  function find_term(t1, t2, t3 : integer) return integer is -- term index
    index : Integer := 1;
    V1 : Integer := t1;
    V2 : Integer := t2;
    V3 : Integer := t3;
    V  : Integer;
  begin
    -- must have largest first for uniqueness
    if V2>V1 then V:=V1; V1:=V2; V2:=V; end if;
    if V3>V1 then V:=V1; V1:=V3; V3:=V; end if;
    if V3>V2 then V:=V2; V2:=V3; V3:=V; end if;

    while V1/=Var1(index) or V2/=Var2(index) or V3/=Var3(index) loop
      if index>neqn+nlin then
        Var1(index) := V1;
        Var2(index) := V2;
        Var3(index) := V3;
        nlin := nlin+1;
        exit;
      end if;
      index := index+1;
    end loop;
    return index;
  end find_term;

  package Real_IO is new Float_IO(Real);
  use Real_IO;
begin
  Put_Line("pde_nl13a.adb running ");
  Put_Line("solve  u^2 u'' + 2 u u' + 3u = f(x)");
  Put_Line("       f(x) = (41/3)-(2116/9)x+(4780/3)x^2-(48976/9)x^3+");
  Put_Line("              (48976.0/9.0)*x*x*x + 9984.0*x*x*x*x -");
  Put_Line("              (88064.0/9.0)*x*x*x*x*x + (14336.0/3.0)*x*x*x*x*x*x -");
  Put_Line("              9984x^4-(88064/9)x^5+(14336/3)x^6-(8192/9)x^7");
  Put_Line(" ");
  Put_Line("   0<x<1  Boundary u(0) = 1, u(1) = 1");
  Put_Line("Analytic solution ");
  Put_Line("U(x) = 1.0 -(20.0/3.0)*x +12.0*x^2 -(16.0/3.0)*x^3");
  Put_Line("xmin="&Real'Image(Xmin)&", xmax="&Real'Image(Xmax));

  Put_Line("nx="&Integer'Image(nx)&" points for numeric derivative");
  hx := (xmax-xmin)/real(nx-1);
  Put_Line("x grid and analytic solution ");
  for i in 1..neqn loop
    xg(i) := xmin+Real(i-1)*hx;
    Ua(i) := uana(xg(i));
    Put("i="&Integer'Image(i)&", Ua(");
    Put(xg(i),3,5,0);
    Put(")=");
    Put(Ua(i),3,5,0);
    New_Line;
  end loop;
  New_Line;

  -- initialize k and f
  for i in 1..neqn loop
    for j in 1..ntot loop
      k(i,j) := 0.0;
    end loop;
    f(i) := 0.0;
  end loop;
  -- set boundary conditions
  k(1,1) := 1.0;
  put_line("k("&Integer'Image(1)&","&Integer'Image(1)&")="&Real'Image(k(1,1)));
  u(1) := uana(xmin);
  f(1) := u(1);
  put_line("f("&real'Image(xmin)&")="&Real'Image(f(1))&
           ", for f at i="&Integer'Image(1));
  k(neqn,neqn) := 1.0;
  put_line("k("&Integer'Image(neqn)&","&Integer'Image(neqn)&
           ")="&Real'Image(k(neqn,neqn)));
  u(neqn) := uana(xmax);
  f(neqn) := u(neqn);
  put_line("f("&real'Image(xmax)&")="&Real'Image(f(neqn))&
           ", for f at i="&Integer'Image(neqn));

  -- initialize terms
  for i in 1..Neqn loop
     Var1(i) := i;
     Var2(i) := 0;
     Var3(i) := 0;
     Vari1(i) := 0;
     Vari2(i) := 0;
  end loop;
  for i in neqn+1..neqn+nmax loop
     Var1(i) := 0;
     Var2(i) := 0;
     Var3(i) := 0;
     Vari1(i) := 0;
     Vari2(i) := 0;
  end loop;

  -- compute entries in stiffness matrix
  -- method
  --  u^2 u'' + 2 u u' + 3u = f(x)  becomes system of non linear equations
  --  at xg(1), xg(2), ... , xg(nx-1)
  --  for example: denoted x1, x2, ..., x5
  --  example case x1 and x5 are boundary and thus value is known
  --  thus U(x1) and U(x5) are known, U(x2) U(x3) and U(x4) are linear DOF
  --
  --U(x2)*U(x2)*cxx(1)*U(x1)+U(x2)*U(x2)*cxx(2)*U(x2)+U(x2)*U(x2)*cxx(3)*U(x3)+
  --U(x2)*U(x2)*cxx(4)*U(x4)+U(x2)*U(x2)*cxx(5)*U(x5)+
  --2*U(x2)*cx(1)*U(x1)+2*U(x2)*cx(2)*U(x2)+2*U(x2)*cx(3)*U(x3)+
  --2*U(x2)*cx(4)*U(x4)+2*U(x2)*cx(5)*U(x5)+
  --3*U(x2)=f(x2)
  --
  --U(x3)*U(x3)*cxx(1)*U(x1)+U(x3)*U(x3)*cxx(2)*U(x2)+U(x3)*U(x3)*cxx(3)*U(x3)+
  --U(x3)*U(x3)*cxx(4)*U(x4)+U(x3)*U(x3)*cxx(5)*U(x5)+
  --2*U(x3)*cx(1)*U(x1)+2*U(x3)*cx(2)*U(x2)+2*U(x3)*cx(3)*U(x3)+
  --2*U(x3)*cx(4)*U(x4)+2*U(x3)*cx(5)*U(x5)+
  --3*U(x3)=f(x3)
  --
  --U(x4)*U(x4)*cxx(1)*U(x1)+U(x4)*U(x4)*cxx(2)*U(x2)+U(x4)*U(x4)*cxx(3)*U(x3)+
  --U(x4)*U(x4)*cxx(4)*U(x4)+U(x4)*U(x4)*cxx(5)*U(x5)+
  --2*U(x4)*cx(1)*U(x1)+2*U(x4)*cx(2)*U(x2)+2*U(x4)*cx(3)*U(x3)+
  --2*U(x4)*cx(4)*U(x4)+2*U(x4)*cx(5)*U(x5)+
  --3*U(x4)=f(x4)
  --

  Put_line("compute non linear stiffness matrix");
  for i in 2..neqn-1 loop -- first=1 and last=neqn are boundary
    x := xg(i);
    rderiv(1, nx, i, hx, cx); -- discrete derivative coefficients
    rderiv(2, nx, i, hx, cxx);

    -- U^2 U'' term
    for jj in 1..nx loop
      if jj=1 or jj=nx then -- involves boundary
        -- compute coefficient contribution for k(i,j)
        j := find_term(i,i,0);
        k(i,j) := k(i,j) + cxx(jj)*u(jj);
        put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&
                 Real'Image(k(i,j)));
      else
        -- compute coefficient contribution for k(i,j)
        j := find_term(i,i,jj);
        k(i,j) := k(i,j) + cxx(jj);
        put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&
                 Real'Image(k(i,j)));
      end if;
    end loop; -- jj

    -- 2 U U' term
    for jj in 1..nx loop
      if jj=1 or jj=nx then -- involves boundary
        -- compute coefficient contribution for k(i,j)
        j := find_term(i,0,0);
        k(i,j) := k(i,j) + 2.0*cx(jj)*u(jj);
        put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&
                 Real'Image(k(i,j)));
      else
        -- compute coefficient contribution for k(i,j)
        j := find_term(i,jj,0);
        k(i,j) := k(i,j) + 2.0*cx(jj);
        put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&
                 Real'Image(k(i,j)));
      end if;
    end loop; -- jj

    -- 3 U term
    k(i,i) := k(i,i) + 3.0;

    -- RHS
    f(i) := FF(x);
  end loop; --i

  Put_Line("k computed stiffness matrix");
  ntot := neqn+nlin; -- now know actual total
  for i in 1..neqn loop
    for j in 1..ntot loop
      Put_Line("i="&Integer'Image(I)&", j="&Integer'Image(J)&", k(i,j)="&
                Real'Image(k(i,j)));
    end loop;
  end loop;
  New_Line;
  Put_Line("f computed forcing function ");
  for i in 1..neqn loop
    Put_Line("f("&Integer'Image(i)&")="&Real'Image(f(i)));
  end loop;
  New_Line;

  -- guess at initial conditions, boundary known
  u(2) := uana(xmin); -- also set up previously
  u(neqn) := uana(xmax);
  for I in 2..neqn-1 loop
    u(i) := 1.0;
  end loop;

  Put_Line("finds a solution, not the one expected");
  Simeq_Newton5(neqn, nlin, k, f, Var1, Var2, Var3, Vari1, Vari2, U);

  avgerr := 0.0;
  maxerr := 0.0;
  Put_Line("u computed from nonlinear equations, Ua analytic, error ");
  for i in 1..neqn loop
    err := abs(u(i)-Ua(i));
    if err>maxerr then  maxerr := err; end if;
    avgerr := avgerr + err;
    Put("u("&Integer'Image(I)&")=");
    Put(u(i),3,5,0);
    Put(", Ua="); Put(Ua(i),3,5,0);
    Put(", err="); Put(u(i)-Ua(i),3,5,0);
    New_Line;
  end loop;
  Put_Line("                   maxerr="&
           Real'Image(maxerr)&", avgerr="&
           Real'image(avgerr/Real(neqn)));
  New_Line;

  Put_line("check_soln on computed solution against PDE");
  check_soln(u);  -- for use when analytic solution is not known
  Put_line("just checking check_soln when given correct solution");
  check_soln(Ua); -- checking on checker
  New_Line;

  -- try giving the correct solution, to check on solver
  Put_line("just checking simeq_newton5 when given correct solution \n");
  for I in 1..Neqn loop
    u(i) := uana(xg(i));
  end loop;
  simeq_newton5(neqn, nlin, k, f, var1, var2, var3, Vari1, Vari2, U);
  Put_Line("pde_nl13a.adb finished");
end pde_nl13a;