-- pde3b.adb  3D non equal h boundary value PDE, iterative solution
--  u(x,y,z) = x^3  + y^3  +z^3  + 2x^2y + 3xy^2 + 4z^2x + 5y + 6
--     du/dx = 3x^2 + 0    + 0   + 4xy   + 3y^2  + 4z^2  + 0  + 0
-- d^2u/dx^2 = 6x   + 0    + 0   + 4y    + 0     + 0     + 0  + 0
--     du/dy = 0    + 3y^2 + 0   + 2x^2  + 6xy   + 0     + 5  + 0
-- d^2u/dy^2 = 0    + 6y   + 0   + 0     + 6x    + 0     + 0  + 0
--     du/dz = 0    + 0    +3z^2 + 0     + 0     + 8zx   + 0  + 0
-- d^2u/dz^2 = 0    + 0    +6z   + 0     + 0     + 8x    + 0  + 0
--
-- solve d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 = c(x,y,z)
-- by second order method on cube -1,-1,-1 to 1.2,1,1

-- second order numerical difference
-- d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 =
-- 1/hx^2(u(x-hx,y,z)+u(x+hx,y,z)-2u(x,y,z)) +
-- 1/hy^2(u(x,y-hy,z)+u(x,y+hy,z)-2u(x,y,z))
-- 1/hz^2(u(x,y,z-hz)+u(x,y,z+hz)-2u(x,y,z)) + O(h^2)
--  = c(x,y,z) = 20x + 10y + 6z
-- solve for new u(x,y,z) =
-- ((u(x-hx,y,z)+u(x+hx,y,z))/hx^2 +
--  (u(x,y-hy,z)+u(x,y+hy,z))/hy^2 +
--  (u(x,y,z-hz)+u(x,y,z+hz))/hz^2 - c(x,y,z)) /
--  (2/hx^2 + 2/hy^2 + 2/hz^2)

with ada.text_io; use ada.text_io;

procedure pde3b is
  type real_matrix is array(integer range <>, integer range <>, integer range <>)
                      of long_float;
  type real_vector is array(integer range <>) of long_float;
  us : real_matrix(0..40, 0..40, 0..40); -- solution being iterated
  ut : real_matrix(0..40, 0..40, 0..40); -- temporary for partial results
                         -- nx-1 by ny-1 by nz-1 internal, non boundary cells
  nx : integer := 10; -- one less than 'C' when Ada subscript starts at zero
  ny : integer := 10;
  nz : integer := 8;
  xmin : long_float := -1.0;
  ymin : long_float := -1.0;
  zmin : long_float := -1.0;
  xmax : long_float :=  1.2;
  ymax : long_float :=  1.0;
  zmax : long_float :=  1.0;
  hx   : long_float;         -- x step = (xmax-xmin)/nx
  hy   : long_float;         -- y step = (ymax-ymin)/ny
  hz   : long_float;         -- z step = (zmax-zmin)/nz
  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;

  function u(x:long_float; y:long_float; z:long_float) return long_float is
           -- solution for boundary and check
  begin
    return x*x*x + y*y*y + z*z*z + 2.0*x*x*y + 3.0*x*y*y +
           4.0*z*z*x + 5.0*y +6.0;
  end u;

  function c(x:long_float; y:long_float; z:long_float) return long_float is
           -- from solve
  begin
    return 20.0*x + 10.0*y + 6.0*z;
  end c;

  function u000(i:integer; j:integer; k:integer) return long_float is
           -- new iterated value
  begin
   return ((us(i-1,j,k) + us(i+1,j,k))/(hx*hx) +
           (us(i,j-1,k) + us(i,j+1,k))/(hy*hy) +
           (us(i,j,k-1) + us(i,j,k+1))/(hz*hz) -
           c(xmin+long_float(i)*hx,ymin+long_float(j)*hy,zmin+long_float(k)*hz)) /
           (2.0/(hx*hx) + 2.0/(hy*hy) + 2.0/(hz*hz));
  end u000;

  procedure printexact is
  begin
    new_line;
    put_line("exact solution");
    for i in 0..nx loop
      for j in 0..ny loop
        for k in 0..nz loop
          put_line("u("&integer'image(i)&","&integer'image(j)&
                   ","&integer'image(k)&")="&long_float'image(
                   u(xmin+long_float(i)*hx,ymin+long_float(j)*hy,zmin+long_float(k)*hz)));
        end loop;
      end loop;
    end loop;
  end printexact;

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

  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
        for k in 1..nz-1 loop
          put_line("diff("&integer'image(i)&","&integer'image(j)&
                        ","&integer'image(k)&")="&long_float'image(
          u(xmin+long_float(i)*hx,ymin+long_float(j)*hy,zmin+long_float(k)*hz)-
          us(i,j,k)));
        end loop;
      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
        for k in 1..nz-1 loop
          ut(i,j,k) := u000(i,j,k);
          diff := diff + abs(ut(i,j,k)-us(i,j,k));
        end loop;
      end loop;
    end loop;
    -- copy temp to solution
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
        for k in 1..nz-1 loop
          us(i,j,k) := ut(i,j,k);
        end loop;
      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
        for k in 1..nz-1 loop
          zexact := u(xmin+long_float(i)*hx,
                      ymin+long_float(j)*hy,
                      zmin+long_float(k)*hz);
          err := abs(us(i,j,k)-zexact);
          error := error + err;
          if err>max_error then
            max_error := err;
          end if;
        end loop;
      end loop;
    end loop;
  end check;

  procedure init_boundary is
  begin
    -- set boundary on six sides
    for i in 0..nx loop
      for j in 0..ny loop
        us(i,j,0)  := u(xmin+long_float(i)*hx, ymin+long_float(j)*hy, zmin);
        us(i,j,nz) := u(xmin+long_float(i)*hx, ymin+long_float(j)*hy, zmax);
      end loop;
    end loop;
    for j in 0..ny loop
      for k in 0..nz loop
        us(0,j,k)  := u(xmin, ymin+long_float(j)*hy, zmin+long_float(k)*hz);
        us(nx,j,k) := u(xmax, ymin+long_float(j)*hy, zmin+long_float(k)*hz);
      end loop;
    end loop;
    for i in 0..nx loop
      for k in 0..nz loop
        us(i,0,k)  := u(xmin+long_float(i)*hx,ymin, zmin+long_float(k)*hz);
        us(i,ny,k) := u(xmin+long_float(i)*hx,ymax, zmin+long_float(k)*hz);
      end loop;
    end loop;
    -- zero internal cells
    for i in 1..nx-1 loop
      for j in 1..ny-1 loop
        for k in 1..nz-1 loop
          us(i,j,k) := 0.0;
        end loop;
      end loop;
    end loop;
  end init_boundary;

begin  -- pde3b

  put_line("pde3b.adb running");
  hx := (xmax-xmin)/long_float(nx); -- step size in x
  hy := (ymax-ymin)/long_float(ny); -- step size in y
  hz := (zmax-zmin)/long_float(nz); -- step size in z
  put_line("xmin="&long_float'image(xmin)&", xmax="&long_float'image(xmax)&
           ", nx="&integer'image(nx)&", hx="&long_float'image(hx));
  put_line("ymin="&long_float'image(ymin)&", ymax="&long_float'image(ymax)&
           ", ny="&integer'image(ny)&", hy="&long_float'image(hy));
  put_line("zmin="&long_float'image(zmin)&", zmax="&long_float'image(zmax)&
           ", nz="&integer'image(nz)&", hz="&long_float'image(hz));

  init_boundary;
  put_line("initial conditions set");

  while diff>0.00001 and iter<149 loop
    iterate;
    check;
    avg_error := error/(long_float(nx-1)*long_float(ny-1)*long_float(nz-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
      printus;
      printdiff;
    end if;
  end loop;
  printexact;
  printus;
  printdiff;
end pde3b;