-- fem_check33_la.adb Galerkin FEM -- solve uxxx(x,y,z) + 2 uyyy(x,y,z) + 3 uzzz(x,y,z) + 4 uxxy(x,y,z) + -- 5 uxxz(x,y,z) + 6 uyyx(x,y,z) + 7 uyyz(x,y,z) + 8 uzzx(x,y,z) + -- 9 uzzy(x,y,z) + 10 uxx(x,y,z) + 11 uyy(x,y,z) + 12 uzz(x,y,z) + -- 13 uxy(x,y,z) + 14 uxz(x,y,z) + 15 uyz(x,y,z) + 16 ux(x,y,z) + -- 17 uy(x,y,z) + 18 uz(x,y,z) + 19 u(x,y,z) = FF(x,y,z) -- FF(x,y,z) = 406.0 + 208.0*x*y*z*z + 240.0*x*x*y*z + 88.0*x*x*z*z + -- 136.0*x*x*y*z*z + 128.0*x*y*y*z*z + 217.0*x + 302.0*y + -- 80.0*y*y*z*z + 144.0*x*x*y*y*z + 96.0*x*x*y*y + 224.0*x*y*y*z + -- 76.0*x*x*y*y*z*z +144.0*x*x*y + 80.0*y*y*z + 64.0*y*z*z + -- 133.0*y*z + 128.0*x*y*y + 112.0*x*x*z + 96.0*x*z*z + 114.0*x*z + -- 95.0*x*y + 431.0*z + 19.0*x*x*x*x + 38.0*y*y*y*y + 57.0*z*z*z*z + -- 64.0*x*x*x + 136.0*y*y*y + 216.0*z*z*z + 120.0*x*x + 264.0*y*y + -- 432.0*z*z; -- boundary conditions computed using u(x,y,z) -- analytic solution may be given by u(x,y,z) = x^4 + 2 y^4 + 3 z^4 + -- 4 x^2 y^2 z^2 + 5 x y + 6 x z + 7 y z + 8 -- Gauss-Legendre integration used -- solve for Uijk numerical approximation U(x_i,y_j,z_k) of u(x_i,y_j,z_k) -- kg * ug = fg, solve simultaneous equations for U with Ada.Text_IO; use Ada.Text_IO; with Calendar; with Real_Arrays; -- types real, real_vector, real_matrix use Real_Arrays; with laphi; -- phi, phip, phipp, phippp, phipppp use laphi; with simeq; -- function simeq with gaulegf; -- function gaulegf procedure fem_check33_la is nx : Integer := 5; ny : Integer := 5; nz : Integer := 5; kg : Real_Matrix(1..nx*ny*nz,1..nx*ny*nz); --((i-1)*ny*nz+(ii-1)*nz+iii,(j-1)*ny*nz+(jj-1)*nz+jjj) xg : Real_Vector(1..nx); yg : Real_Vector(1..ny); zg : Real_Vector(1..nz); fg : Real_Vector(1..nx*ny*nz); -- (i-1)*ny*nz+(ii-1)*nz+iii ug : Real_Vector(1..nx*ny*nz); Ua : Real_Vector(1..nx*ny*nz); xmin : real := 0.0; xmax : real := 1.0; ymin : real := 0.0; ymax : real := 1.0; zmin : real := 0.0; zmax : real := 1.0; x, y, z, hx, hy, hz : real; npx : Integer := 12; npy : Integer := 12; npz : Integer := 12; xx : Real_Vector(1..npx); -- for Gauss-Legendre wx : Real_Vector(1..npx); yy : Real_Vector(1..npy); wy : Real_Vector(1..npy); zz : Real_Vector(1..npz); wz : Real_Vector(1..npz); i, ii , iii : Integer; val, err, avgerr, maxerr : real; time_start, time_now, time_last : Duration; function s(i:Integer; ii:Integer; iii:integer) return Integer is -- general subscript begin return (i-1)*ny*nz+(ii-1)*nz+iii; end s; function FF(x,y,z:real) return real is -- forcing function begin return 406.0 + 208.0*x*y*z*z + 240.0*x*x*y*z + 88.0*x*x*z*z + 136.0*x*x*y*z*z + 128.0*x*y*y*z*z + 217.0*x + 302.0*y + 80.0*y*y*z*z + 144.0*x*x*y*y*z + 96.0*x*x*y*y + 224.0*x*y*y*z + 76.0*x*x*y*y*z*z +144.0*x*x*y + 80.0*y*y*z + 64.0*y*z*z + 133.0*y*z + 128.0*x*y*y + 112.0*x*x*z + 96.0*x*z*z + 114.0*x*z + 95.0*x*y + 431.0*z + 19.0*x*x*x*x + 38.0*y*y*y*y + 57.0*z*z*z*z + 64.0*x*x*x + 136.0*y*y*y + 216.0*z*z*z + 120.0*x*x + 264.0*y*y + 432.0*z*z; end FF; function uana(x,y,z:real) return real is -- analytic solution and boundary begin return x**4 + 2.0*y**4 + 3.0*z**4 + 4.0*x*x*y*y*z*z + 5.0*x*y + 6.0*x*z + 7.0*y*z + 8.0; end uana; function galk(x,y,z:Real; i,j,ii,jj,iii,jjj:integer) return real is -- Galerkin k stiffness function begin -- for this specific problem return(1.0 * phippp(x,j,nx,xg)* phi(y,jj,ny,yg)* phi(z,jjj,nz,zg) + 2.0 * phi(x,j,nx,xg)*phippp(y,jj,ny,yg)* phi(z,jjj,nz,zg) + 3.0 * phi(x,j,nx,xg)* phi(y,jj,ny,yg)*phippp(z,jjj,nz,zg) + 4.0 * phipp(x,j,nx,xg)* phip(y,jj,ny,yg)* phi(z,jjj,nz,zg) + 5.0 * phipp(x,j,nx,xg)* phi(y,jj,ny,yg)* phip(z,jjj,nz,zg) + 6.0 * phip(x,j,nx,xg)* phipp(y,jj,ny,yg)* phi(z,jjj,nz,zg) + 7.0 * phi(x,j,nx,xg)* phipp(y,jj,ny,yg)* phip(z,jjj,nz,zg) + 8.0 * phip(x,j,nx,xg)* phi(y,jj,ny,yg)* phipp(z,jjj,nz,zg) + 9.0 * phi(x,j,nx,xg)* phip(y,jj,ny,yg)* phipp(z,jjj,nz,zg) + 10.0 * phipp(x,j,nx,xg)* phi(y,jj,ny,yg)* phi(z,jjj,nz,zg) + 11.0 * phi(x,j,nx,xg)* phipp(y,jj,ny,yg)* phi(z,jjj,nz,zg) + 12.0 * phi(x,j,nx,xg)* phi(y,jj,ny,yg)* phipp(z,jjj,nz,zg) + 13.0 * phip(x,j,nx,xg)* phip(y,jj,ny,yg)* phi(z,jjj,nz,zg) + 14.0 * phip(x,j,nx,xg)* phi(y,jj,ny,yg)* phip(z,jjj,nz,zg) + 15.0 * phi(x,j,nx,xg)* phip(y,jj,ny,yg)* phip(z,jjj,nz,zg) + 16.0 * phip(x,j,nx,xg)* phi(y,jj,ny,yg)* phi(z,jjj,nz,zg) + 17.0 * phi(x,j,nx,xg)* phip(y,jj,ny,yg)* phi(z,jjj,nz,zg) + 18.0 * phi(x,j,nx,xg)* phi(y,jj,ny,yg)* phip(z,jjj,nz,zg) + 19.0 * phi(x,j,nx,xg)* phi(y,jj,ny,yg)* phi(z,jjj,nz,zg) )* (phi(x,i,nx,xg)* phi(y,ii,ny,yg)* phi(z,iii,nz,zg)); end galk; function galf(x,y,z:Real; i,ii,iii:integer) return real is -- Galerkin f function begin return FF(x,y,z)*phi(x,i,nx,xg)*phi(y,ii,ny,yg)*phi(z,iii,nz,zg); end galf; package Real_IO is new Float_IO(Real); use Real_IO; package Duration_IO is new Fixed_IO(Duration); use Duration_IO; begin Put_Line("fem_check33_la.adb running "); Put_Line("Solve uxxx(x,y,z) +2 uyyy(x,y,z) +3 uzzz(x,y,z) +4 uxxy(x,y,z)+"); Put_Line(" 5 uxxz(x,y,z) +6 uyyx(x,y,z) +7 uyyz(x,y,z) +8 uzzx(x,y,z)+"); Put_Line(" 9 uzzy(x,y,z) +10 uxx(x,y,z) +11 uyy(x,y,z) +12 uzz(x,y,z)+"); Put_Line(" 13 uxy(x,y,z) +14 uxz(x,y,z) +15 uyz(x,y,z) +16 ux(x,y,z)+"); Put_Line(" 17 uy(x,y,z) +18 uz(x,y,z) +19 u(x,y,z) = FF(x,y,z)"); Put_Line("FF(x,y,z)="); Put_Line(" 406.0 + 208.0*x*y*z*z + 240.0*x*x*y*z + 88.0*x*x*z*z +"); Put_Line(" 136.0*x*x*y*z*z + 128.0*x*y*y*z*z + 217.0*x + 302.0*y +"); Put_Line(" 80.0*y*y*z*z + 144.0*x*x*y*y*z + 96.0*x*x*y*y + 224.0*x*y*y*z +"); Put_Line(" 76.0*x*x*y*y*z*z +144.0*x*x*y + 80.0*y*y*z + 64.0*y*z*z +"); Put_Line(" 133.0*y*z + 128.0*x*y*y + 112.0*x*x*z + 96.0*x*z*z + 114.0*x*z +"); Put_Line(" 95.0*x*y + 431.0*z + 19.0*x*x*x*x + 38.0*y*y*y*y + 57.0*z*z*z*z +"); Put_Line(" 64.0*x*x*x + 136.0*y*y*y + 216.0*z*z*z + 120.0*x*x + 264.0*y*y +"); Put_Line(" 432.0*z*z;"); Put_Line("boundary conditions computed using u(x,y,z)"); Put_Line("analytic solution given by u(x,y,z)=x^4 + 2 y^4 + 3 z^4 + "); Put_Line(" 4 x^2 y^2 z^2 + 5 x y + 6 x z + 7 y z + 8"); Put_Line("Gauss-Legendre integration used"); New_Line; Put("xmin="); Put(xmin,2,2,0); Put(", xmax="); Put(xmax,2,2,0); Put(", ymin="); Put(ymin,2,2,0); Put(", ymax="); Put(ymax,2,2,0); Put(", zmin="); Put(zmin,2,2,0); Put(", zmax="); Put(zmax,2,2,0); New_Line; Put_Line("nx="&Integer'Image(nx)&", ny="&Integer'Image(ny)& ", nz="&Integer'Image(nz)); time_start := Calendar.Seconds(Calendar.Clock); Put(time_start); Put_Line("= start time in seconds"); time_last := time_start; -- initialize xg array. does not have to be uniform spacing hx := (xmax-xmin)/real(nx-1); Put_Line("x grid and analytic solution at ymin,zmin "); for i in 1..nx loop xg(i) := xmin+Real(i-1)*hx; Ua(i) := uana(xg(i),ymin,zmin); Put("i="&Integer'Image(i)&", Ua("); Put(xg(i),2,5,0); Put(","); Put(ymin,2,5,0); Put(","); Put(zmin,2,5,0); Put(")="); Put(Ua(i),3,5,0); New_Line; end loop; New_Line; -- i -- initialize yg array. does not have to be uniform spacing hy := (ymax-ymin)/real(ny-1); Put_Line("y grid and analytic solution at xmin,zmin"); for ii in 1..ny loop yg(ii) := ymin+Real(ii-1)*hy; Ua(ii) := uana(xmin,yg(ii),zmin); Put("ii="&Integer'Image(ii)&", Ua("); Put(xmin,2,5,0); Put(","); Put(yg(ii),2,5,0); Put(","); Put(zmin,2,5,0); Put(")="); Put(Ua(ii),3,5,0); New_Line; end loop; New_Line; -- ii -- initialize zg array. does not have to be uniform spacing hz := (zmax-zmin)/real(nz-1); Put_Line("z grid and analytic solution at xmin, ymin"); for iii in 1..nz loop zg(iii) := zmin+Real(iii-1)*hz; Ua(iii) := uana(xmin,ymin,zg(iii)); Put("iii="&Integer'Image(iii)&", Ua("); Put(xmin,2,5,0); Put(","); Put(ymin,2,5,0); Put(","); Put(zg(iii),2,5,0); Put(")="); Put(Ua(iii),3,5,0); New_Line; end loop; -- iii New_Line; -- put solution in Ua vector for i in 1..nx loop x := xmin+Real(i-1)*hx; for ii in 1..ny loop y := ymin+Real(ii-1)*hy; for iii in 1..nz loop z := zmin+Real(iii-1)*hz; Ua(s(i,ii,iii)) := uana(x,y,z); Put("solution at i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(Ua(s(i,ii,iii)),3,3,0); New_Line; end loop; -- iii end loop; -- ii end loop; -- i -- initialize k and f for i in 1..nx loop for j in 1..nx loop for ii in 1..ny loop for jj in 1..ny loop for iii in 1..nz loop for jjj in 1..nz loop kg(s(i,ii,iii),s(j,jj,jjj)) := 0.0; end loop; -- jjj end loop; -- iii end loop; -- jj end loop; -- ii end loop; -- j end loop; -- i for i in 1..nx loop for ii in 1..ny loop for iii in 1..nz loop fg(s(i,ii,iii)) := 0.0; end loop; -- iii end loop; -- ii end loop; -- i -- set boundary conditions, six sides for i in 1..nx loop x := xg(i); for ii in 1..ny loop y := yg(ii); iii := 1; z := zg(iii); kg(s(i,ii,iii),s(i,ii,iii)) := 1.0; fg(s(i,ii,iii)) := uana(x,y,z); Put("boundary i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(fg(s(i,ii,iii))); New_Line; end loop; -- ii for ii in 1..ny loop y := yg(ii); iii := nz; z := zg(iii); kg(s(i,ii,iii),s(i,ii,iii)) := 1.0; fg(s(i,ii,iii)) := uana(x,y,z); Put("boundary i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(fg(s(i,ii,iii))); New_Line; end loop; -- ii for iii in 1..nz loop z := zg(iii); ii := 1; y := yg(ii); kg(s(i,ii,iii),s(i,ii,iii)) := 1.0; fg(s(i,ii,iii)) := uana(x,y,z); Put("boundary i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(fg(s(i,ii,iii))); New_Line; end loop; -- iii for iii in 1..nz loop z := zg(iii); ii := ny; y := yg(ii); kg(s(i,ii,iii),s(i,ii,iii)) := 1.0; fg(s(i,ii,iii)) := uana(x,y,z); Put("boundary i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(fg(s(i,ii,iii))); New_Line; end loop; -- iii end loop; -- nx for ii in 1..ny loop y := yg(ii); for i in 1..nx loop x := xg(i); iii := 1; z := zg(iii); kg(s(i,ii,iii),s(i,ii,iii)) := 1.0; fg(s(i,ii,iii)) := uana(x,y,z); Put("boundary i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(fg(s(i,ii,iii))); New_Line; end loop; -- i for i in 1..nx loop x := xg(i); iii := nz; z := zg(iii); kg(s(i,ii,iii),s(i,ii,iii)) := 1.0; fg(s(i,ii,iii)) := uana(x,y,z); Put("boundary i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(fg(s(i,ii,iii))); New_Line; end loop; -- i for iii in 1..nz loop z := zg(iii); i := 1; x := xg(i); kg(s(i,ii,iii),s(i,ii,iii)) := 1.0; fg(s(i,ii,iii)) := uana(x,y,z); Put("boundary i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(fg(s(i,ii,iii))); New_Line; end loop; -- iii for iii in 1..nz loop z := zg(iii); i := nx; x := xg(i); kg(s(i,ii,iii),s(i,ii,iii)) := 1.0; fg(s(i,ii,iii)) := uana(x,y,z); Put("boundary i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(fg(s(i,ii,iii))); New_Line; end loop; -- iii end loop; -- ny for iii in 1..nz loop z := zg(iii); for i in 1..nx loop x := xg(i); ii := 1; y := yg(ii); kg(s(i,ii,iii),s(i,ii,iii)) := 1.0; fg(s(i,ii,iii)) := uana(x,y,z); Put("boundary i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(fg(s(i,ii,iii))); New_Line; end loop; -- i for i in 1..nx loop x := xg(i); ii := ny; y := yg(ii); kg(s(i,ii,iii),s(i,ii,iii)) := 1.0; fg(s(i,ii,iii)) := uana(x,y,z); Put("boundary i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(fg(s(i,ii,iii))); New_Line; end loop; -- i for ii in 1..ny loop y := yg(ii); i := 1; x := xg(i); kg(s(i,ii,iii),s(i,ii,iii)) := 1.0; fg(s(i,ii,iii)) := uana(x,y,z); Put("boundary i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(fg(s(i,ii,iii))); New_Line; end loop; -- ii for ii in 1..ny loop y := yg(ii); i := nx; x := xg(i); kg(s(i,ii,iii),s(i,ii,iii)) := 1.0; fg(s(i,ii,iii)) := uana(x,y,z); Put("boundary i="&Integer'Image(i)); Put(", x="); Put(x,2,3,0); Put(", ii="&Integer'Image(ii)); Put(", y="); Put(y,2,3,0); Put(", iii="&Integer'Image(iii)); Put(", z="); Put(z,2,3,0); Put(" is "); Put(fg(s(i,ii,iii))); New_Line; end loop; -- ii end loop; -- nz -- integrate over complete domain Put_Line("calling gaulegf npx="&Integer'Image(npx)); gaulegf(xmin, xmax, xx, wx, npx); -- xx and wx set up for integrations Put_Line("calling gaulegf npy="&Integer'Image(npy)); gaulegf(ymin, ymax, yy, wy, npy); -- yy and wy set up for integrations Put_Line("calling gaulegf npz="&Integer'Image(npz)); gaulegf(zmin, zmax, zz, wz, npz); -- zz and wz set up for integrations -- debug Put_Line("xx(1)="&Real'Image(xx(1))&" xx(2)="&Real'Image(xx(2))& " wx(1)="&Real'Image(wx(1))&" wx(2)="&Real'Image(wx(2))); Put_Line("yy(1)="&Real'Image(yy(1))&" yy(2)="&Real'Image(yy(2))& " wy(1)="&Real'Image(wy(1))&" wy(2)="&Real'Image(wy(2))); Put_Line("zz(1)="&Real'Image(zz(1))&" zz(2)="&Real'Image(zz(2))& " wz(1)="&Real'Image(wz(1))&" wz(2)="&Real'Image(wz(2))); val := galk(xx(2),yy(2),zz(2),2,2,2,2,2,2); -- six 1's in C Put_Line("galk(xx(2),yy(2),zz(2))="&Real'Image(val)& ", at i="&Integer'Image(2)&", j="&Integer'Image(2)& ", ii="&Integer'Image(2)&", jj="&Integer'Image(2)& ", iii="&Integer'Image(2)&", jjj="&Integer'Image(2)); val := galf(xx(2),yy(2),zz(2),2,2,2); -- three 1's in C Put_Line("galf(xx(2),yy(2),zz(2))="&Real'Image(val)); -- compute entries in stiffness matrix Put_Line("compute stiffness matrix"); for i in 2..nx-1 loop for j in 1..nx loop for ii in 2..ny-1 loop for jj in 1..ny loop for iii in 2..nz-1 loop for jjj in 1..nz loop -- compute integral for ks(i,j) val := 0.0; for px in 1..npx loop for py in 1..npy loop for pz in 1..npz loop val := val + wx(px)*wy(py)*wz(pz)* galk(xx(px),yy(py),zz(pz),i,j,ii,jj,iii,jjj); end loop; --pz end loop; --py end loop; --px Put_Line("Legendre integration="&Real'Image(val)& ", at i="&Integer'Image(i)&", j="&Integer'Image(j)& ", ii="&Integer'Image(ii)&", jj="&Integer'Image(jj)& ", iii="&Integer'Image(iii)&", jjj="&Integer'Image(jjj)); kg(s(i,ii,iii),s(j,jj,jjj)) := val; end loop; -- jjj end loop; -- iii end loop; -- jj end loop; -- ii end loop; -- j end loop; -- i -- compute integral for f(i) for i in 2..nx-1 loop for ii in 2..ny-1 loop for iii in 2..nz-1 loop val := 0.0; for px in 1..npx loop for py in 1..npy loop for pz in 1..npx loop val := val + wx(px)*wy(py)*wz(pz)* galf(xx(px),yy(py),zz(pz),i,ii,iii); end loop; -- pz end loop; -- py end loop; -- px Put_Line("Legendre integration="&Real'Image(val)& ", i="&Integer'Image(i)&", ii="&Integer'Image(ii)& ", iii="&Integer'Image(iii)); fg(s(i,ii,iii)) := val; end loop; -- iii end loop; -- ii end loop; -- i Put_Line("k computed stiffness matrix, see above"); Put_Line("f computed forcing function, see above "); Put_Line("solving "&Integer'Image(nx*ny*nz)&" in "&Integer'Image(nx*ny*nz)& "unknowns"); ug := simeq(kg, fg); Put_Line("equations solved"); time_now := Calendar.Seconds(Calendar.Clock); Put(time_now); Put("= time now, "); Put(time_now-time_last); Put_Line("= seconds for previous section"); time_last := time_now; avgerr := 0.0; maxerr := 0.0; Put_Line("ug computed Galerkin, Ua analytic, error "); for i in 1..nx loop for ii in 1..ny loop for iii in 1..nz loop err := abs(ug(s(i,ii,iii))-Ua(s(i,ii,iii))); if err>maxerr then maxerr := err; end if; avgerr := avgerr + err; Put("ug("&Integer'Image(i)&","&Integer'Image(ii)&")="); Put(ug(s(i,ii,iii)),3,5,0); Put(", Ua="); Put(Ua(s(i,ii,iii)),3,5,0); Put(", err="); Put(err,3,5,0); New_Line; end loop; end loop; end loop; time_now := Calendar.Seconds(Calendar.Clock); Put(time_now-time_start); Put_Line("= total CPU time in seconds"); Put_Line(" maxerr="& Real'Image(maxerr)&", avgerr="& Real'image(avgerr/Real(nx*ny*nz))); Put_Line("end fem_check33_la.adb"); end fem_check33_la;