-- pde_nl13.adb  solve non linear second order, third degree PDE
--               improvement: compute nlin and Var's based on nx
--               this finds one of multiple solutions, not the expected
--
-- 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_nl13 is
  nx : Integer := 5;    -- number of coefficients for discrete derivative
  neqn : Integer := nx; -- number of equations and number of linear unknowns
  nlin : Integer := 15; -- number of nonlinear terms, specific to this PDE, nx
  ntot : Integer := neqn+nlin; -- 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;
  i, j : Integer;

  -- linear, then nonlinear terms
  Var1 : Integer_vector(1..Ntot) :=
           ( 1, 2, 3, 4, 5, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4);
  Var2 : Integer_Vector(1..Ntot) :=
           ( 0, 0, 0, 0, 0, 2, 3, 4, 3, 4, 2, 2, 3, 4, 2, 3, 4, 2, 3, 4);
  Var3 : Integer_vector(1..Ntot) :=
           ( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 3, 3, 3, 4, 4, 4);
  Vari1 : Integer_vector(1..Ntot) :=
           ( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
  Vari2 : Integer_vector(1..Ntot) :=
           ( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);

  -- 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 Print_Equations(n : integer; nlin : integer;
                            Var1 : integer_vector; Var2 : integer_vector;
                            Var3 : integer_vector; Vari1 : integer_vector;
                            Vari2 : Integer_Vector) is
    Xnames : array(1..5) of String(1..2) := ("X1", "X2", "X3", "X4", "X5");
  begin
    Put_Line("system of equations to be solved, i=1.."&Integer'Image(n));
    for j in 1..n loop
      Put_line("k(i,"&Integer'Image(j)&")*"&Xnames(j)&"+");
    end loop;
    for j in n+1..n+nlin loop
      Put("k(i,"&Integer'Image(J)&")");
      if Var1(J)>0 then
         Put("*"&Xnames(Var1(J)));
      end if;
      if Var2(J)>0 then
         Put("*"&Xnames(Var2(J)));
      end if;
      if Var3(J)>0 then
         Put("*"&Xnames(Var3(J)));
      end if;
      if Vari1(J)>0  then
         Put("/("&Xnames(Vari1(J)));
      end if;
      if Vari2(J)>0  then
         Put("*"&Xnames(Vari2(J)));
      end if;
      if Vari1(J)>0 or Vari2(J)>0 then
         Put(")");
      end if;
      if J /= n+nlin then
        Put_Line("+");
      else
        New_Line;
      end if;
    end loop;
    Put_Line(" = f(i)");
    New_Line;
  end Print_Equations;

  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;

  package Real_IO is new Float_IO(Real);
  use Real_IO;
begin
  Put_Line("pde_nl13.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("one possible 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));

  -- 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)  denoted x2, x3, ..., x5
  --  special 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)
  --
  -- The linear and non linear unknown terms are:
  -- U(x2)              U(x3)              U(x4)
  -- U(x2)*U(x2)        U(x3)*U(x3)        U(x4)*U(x4)
  -- U(x2)*U(x3)        U(x3)*U(x4)        U(x4)*U(x2)
  -- U(x2)*U(x2)*U(x2)  U(x3)*U(x3)*U(x2)  U(x4)*U(x4)*U(x2)
  -- U(x2)*U(x2)*U(x3)  U(x3)*U(x3)*U(x3)  U(x4)*U(x4)*U(x3)
  -- U(x2)*U(x2)*U(x4)  U(x3)*U(x3)*U(x4)  U(x4)*U(x4)*U(x4)
  --
  -- The coefficients of the above terms are:
  -- 3+                 3+                 3+
  -- 2*cx(1)*U(x1)+
  -- 2*cx(5)*U(x5)
  -- cxx(1)*U(x1)+      cxx(1)*U(x1)+      cxx(1)*U(x1)+
  -- cxx(5)*U(x5)+      cxx(5)*U(x5)       cxx(5)*U(x5)
  -- 2*cx(2)            2*cx(3)            2*cx(4)
  -- 2*cx(3)+2*cx(2)    2*cx(4)+2*cx(3)    2*cx(2)+2*cx(4)
  -- cxx(2)             cxx(2)             cxx(2)
  -- cxx(3)             cxx(3)             cxx(3)
  -- cxx(4)             cxx(4)             cxx(4)
  --
  -- Var1() = ( 1, 2, 3, 4, 5, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4);
  -- Var2() = ( 0, 0, 0, 0, 0, 2, 3, 4, 3, 4, 2, 2, 3, 4, 2, 3, 4, 2, 3, 4);
  -- Var3() = ( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 3, 3, 3, 4, 4, 4);
  -- j index    1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20


  Put_line("compute non linear stiffness matrix");
  i := 2; -- 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);
  -- compute val for k(i,j)
  j := 2;  -- U(x2) term
  k(i,j) := k(i,j) + 3.0 + 2.0*cx(1)*U(1) + 2.0*cx(neqn)*U(neqn);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 6;  -- U(x2)*U(x2) term
  k(i,j) := k(i,j) + k(i,j) + cxx(1)*U(1) + cxx(neqn)*U(neqn) + 2.0*cx(2);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 9;  -- U(x2)*U(x3) term
  k(i,j) := k(i,j) + 2.0*cx(3);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 11; -- U(x4)*U(x2) term
  k(i,j) := k(i,j) + 2.0*cx(4);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 12;  -- U(x2)*U(x2)*U(x2) term
  k(i,j) := k(i,j) + cxx(2);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 15; -- U(x2)*U(x2)*U(x3) term
  k(i,j) := k(i,j) + cxx(3);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 18; -- U(x2)*U(x2)*U(x4) term
  k(i,j) := k(i,j) + cxx(4);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  f(i) := FF(x); -- fill in RHS for f(i)
  put_line("f("&real'Image(X)&")="&Real'Image(f(i))&
           ", for f at i="&Integer'Image(i));


  I := 3; -- actually the third equation
  x := xg(i);
  rderiv(1, nx, i, hx, cx); -- discrete derivative coefficients
  rderiv(2, nx, i, hx, cxx);
  -- compute val for k(i,j)
  j := 3;  -- U(x3) term
  k(i,j) := k(i,j) + 3.0 + 2.0*cx(1)*U(1) + 2.0*cx(neqn)*U(neqn);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 7;  -- U(x3)*U(x3) term
  k(i,j) := k(i,j) + cxx(1)*U(1) + cxx(neqn)*U(neqn) + 2.0*cx(3);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 9;  -- U(x3)*U(x2) term
  k(i,j) := k(i,j) + 2.0*cx(2);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 10;  -- U(x3)*U(x4) term
  k(i,j) := k(i,j) + 2.0*cx(4);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 13; -- U(x3)*U(x3)*U(x2) term
  k(i,j) := k(i,j) + cxx(2);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 16; -- U(x3)*U(x3)*U(x3) term
  k(i,j) := k(i,j) + cxx(3);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 19; -- U(x3)*U(x3)*U(x4) term
  k(i,j) := k(i,j) + cxx(4);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  f(i) := FF(x); -- fill in RHS for f(i)
  put_line("f("&real'Image(X)&")="&Real'Image(f(i))&
           ", for f at i="&Integer'Image(i));


  I := 4; -- third equation
  x := xg(i);
  rderiv(1, nx, i, hx, cx); -- discrete derivative coefficients
  rderiv(2, nx, i, hx, cxx);
  -- compute val for k(i,j)
  j := 4;  -- U(x4) term
  k(i,j) := k(i,j) + 3.0 + 2.0*cx(1)*U(1) + 2.0*cx(neqn)*U(neqn);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 8;  -- U(x4)*U(x4) term
  k(i,j) := k(i,j) + cxx(1)*U(1) + cxx(neqn)*U(neqn) + 2.0*cx(4);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 10;  -- U(x4)*U(x3) term
  k(i,j) := k(i,j) + 2.0*cx(3);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 11; -- U(x4)*U(x2) term
  k(i,j) := k(i,j) + 2.0*cx(2);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 14; -- U(x4)*U(x4)*U(x2) term
  k(i,j) := k(i,j) + cxx(2);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 17; -- U(x4)*U(x4)*U(x3) term
  k(i,j) := k(i,j) + cxx(3);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  j := 20; -- U(x4)*U(x4)*U(x4) term
  k(i,j) := k(i,j) + cxx(4);
  put_line("k("&Integer'Image(i)&","&Integer'Image(j)&")="&Real'Image(k(i,j)));
  f(i) := FF(x); -- fill in RHS for f(i)
  put_line("f("&real'Image(X)&")="&Real'Image(f(i))&
           ", for f at i="&Integer'Image(i));


  Put_Line("k computed stiffness matrix");
  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;

  Print_Equations(neqn, nlin, Var1, Var2, Var3, Vari1, Vari2);

  -- guess at initial conditions, boundary known
  U(1) := 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("one possible solution found, not unique, not expected");
  Simeq_Newton5(neqn, nlin, k, f, Var1, Var2, Var3, Vari1, Vari2,
                U, 10, 0.001, 1);

  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, 6, 0.00001, 1);
  Put_Line("pde_nl13.adb finished");
end pde_nl13;