-- pde3.adb  just checking second order approximation         
--      z(x,y) = x^3  + y^3  + 2x^2y + 3xy^2 + 4x + 5y + 6  
--       dz/dx = 3x^2 + 0    + 4xy   + 3y^2  + 4  + 0  + 0  
--   d^2z/dx^2 = 6x   + 0    + 4y    + 0     + 0  + 0  + 0  
--       dz/dy = 0    + 3y^2 + 2x^2  + 6xy   + 0  + 5  + 0  
--   d^2z/dy^2 = 0    + 6y   + 0     + 6x    + 0  + 0  + 0  
--                                                          
-- solve d^2z/dx^2 + d^2z/dy^2 = 12x + 10y = c(x,y)         
-- by second order method on square -1,-1 to 1,1            

-- second order numerical difference
-- d^2z/dx^2 + d^2z/dy^2 = 
-- 1/h^2(z(x-h,y)-2z(x,y)+z(x+h,y))+
-- 1/h^2(z(x,y-h)-2z(x,y)+z(x,y+h)) + O(h^2)
--  = c(x,y)
-- solve for new z(x,y) = 
-- (z(x-h,y)+z(x+h,y)+z(x,y-h)+z(x,y+h)-c(x,y)*h^2)/4

with ada.text_io; use ada.text_io;

procedure pde3 is
  type real_matrix is array(integer range <>, integer range <>) of long_float;
  type real_vector is array(integer range <>) of long_float;
  zs : real_matrix(0..100, 0..100); -- solution being iterated
  zt : real_matrix(0..100, 0..100); -- temporary for partial results
  nx : integer := 10; -- one less than 'C' when Ada subscript starts at zero
  ny : integer := 10; -- nx-1 by ny-1 internal, non boundary cells 
  xmin : long_float := -1.0;
  ymin : long_float := -1.0;
  xmax : long_float :=  1.0;
  ymax : long_float :=  1.0;
  h    : long_float;         -- x  and  y step = (xmax-xmin)/(n-1) = (ymax-ymin) 
  diff : long_float := 99.0; -- sum abs difference form last iteration 
  error : long_float;        -- sum of absolute exact-computed         
  avg_error : long_float;    -- average error                          
  max_error : long_float;    -- maximum cell error
  iter : integer := 0;       -- iteration count
  
  function z(x:long_float; y:long_float) return long_float is -- solution for boundary and check 
  begin
    return x*x*x + y*y*y + 2.0*x*x*y + 3.0*x*y*y + 4.0*x + 5.0*y +6.0;
  end z;

  function c(x:long_float; y:long_float) return long_float is -- from solve
  begin
    return 12.0*x + 10.0*y;
  end c;
 
  function z00(i:integer; j:integer) return long_float is -- new iterated value 
  begin
   return ( zs(i-1,j) + zs(i+1,j) + zs(i,j-1) + zs(i,j+1) - 
           c(xmin+long_float(i)*h,ymin+long_float(j)*h)*h*h )/4.0;
  end z00;

  procedure printexact is
  begin
    new_line;
    put_line("exact solution");
    for i in 0..nx loop
      for j in 0..ny loop
        put_line("z("&integer'image(i)&","&integer'image(j)&")="&
                 long_float'image(z(xmin+long_float(i)*h,ymin+long_float(j)*h)));
      end loop;
    end loop;
  end printexact;

  procedure printzs is
  begin
    new_line;
    put_line("computed solution");
    for i in 0..nx loop
      for j in 0..ny loop
        put_line("zs("&integer'image(i)&","&integer'image(j)&")="&
                 long_float'image(zs(i,j)));
      end loop;
    end loop;
  end printzs;

  procedure printdiff is
  begin
    new_line;
    put_line("difference exact-computed solution"); 
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
        put_line("diff("&integer'image(i)&","&integer'image(j)&")="&
        long_float'image(z(xmin+long_float(i)*h,ymin+long_float(j)*h)-zs(i,j)));
      end loop;
    end loop; 
  end printdiff;

  procedure iterate is
  begin
    diff := 0.0;
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
        zt(i,j) := z00(i,j);
        diff := diff + abs(zt(i,j)-zs(i,j));
      end loop;
    end loop;
    -- copy temp to solution 
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
        zs(i,j) := zt(i,j);
      end loop;
    end loop;
  end iterate;

  procedure check is
    zexact, err : long_float;
  begin
    error := 0.0;
    max_error := 0.0;
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
        zexact := z(xmin+long_float(i)*h,ymin+long_float(j)*h);
        err := abs(zs(i,j)-zexact);
        error := error + err;
        if err>max_error then
          max_error := err;
        end if;
      end loop;
    end loop;
  end check;

  procedure init_boundary is
  begin
    -- set boundary on four sides 
    for i in 0..nx loop
      zs(i,0)  := z(xmin+long_float(i)*h, ymin);
      zs(i,ny) := z(xmin+long_float(i)*h, ymax);
    end loop;
    for j in 0..ny loop
      zs(0,j)  := z(xmin, ymin+long_float(j)*h);
      zs(nx,j) := z(xmax, ymin+long_float(j)*h);
    end loop;
    -- zero internal cells 
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
        zs(i,j) := 0.0;
      end loop;
    end loop;
  end init_boundary;

begin  -- pde3
  
  put_line("pde3.adb running");
  h  := (xmax-xmin)/long_float(nx); -- step size in x 
  -- := (ymax-ymin)/long_float(ny); -- step size in y 
  put_line("xmin="&long_float'image(xmin)&", xmax="&long_float'image(xmax)&
           ", nx="&integer'image(nx)&", h="&long_float'image(h));
  put_line("ymin="&long_float'image(ymin)&", ymax="&long_float'image(ymax)&
           ", ny="&integer'image(ny)&", h="&long_float'image(h));
  
  printexact;
  init_boundary;
  new_line;
  put_line("initial condition");
  printzs;
  
  while diff>0.0001 and iter<200 loop
    iterate;
    check;
    avg_error := error/(long_float(nx-1)*long_float(ny-1));
    put_line("iter="&integer'image(iter)&", diff="&long_float'image(diff)&
             ", error="&long_float'image(error)&" avg="&long_float'image(avg_error)&
             ", max="&long_float'image(max_error));
    iter := iter+1;
    if iter mod 50 =0  then
      printzs;
      printdiff;
    end if;
  end loop;
end pde3;