-- fem_check44_la.adb Galerkin FEM -- solve uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) + -- 4 utttt(x,y,z,t) + 5 uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) + -- 7 uztt(x,y,z,t) + 8 uxxx(x,y,z,t) + 9 uyyt(x,y,z,t) + -- 10 uzt(x,y,z,t) + 11 ut(x,y,z,t) + 12 u(x,y,z,t) = FF(x,y,z,t) -- FF(x,y,z,t) 706.0 + 120.0*t*t + 140.0*y*y*Z + 768.0*x + 180.0*z*z + -- 200.0*y*y*z*t + 44.0*t*t*t + 110.0*y*y*z*z*t + 66.0*x*y + -- 12.0*t*t*t*t + 24.0*z*z*z*z + 36.0*Y*Y*Y*Y + 48.0*x*x*x*x + -- 60.0*y*y*z*z*t*t + 72.0*x*y*t + 84.0*x*z + 96.0*t ; -- boundary conditions computed using u(x,y,z,t) -- analytic solution may be given by u(x,y,z,t) = t^4 + 2 z^4 + 3 y^4 + -- 4 x^4 + 5 y^2 z^2 t^2 + 6 x y t + 7 x z + 8 t + 9 -- Gauss-Legendre integration used -- solve for Uijkl numerical approximation U(x_i,y_j,z_k,t_l) -- kg * ug = fg, solve simultaneous equations for U with Ada.Text_IO; use Ada.Text_IO; with Ada.Calendar; use Ada.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_check44_la is nx : Integer := 5; ny : Integer := 5; nz : Integer := 5; nt : Integer := 5; kg : Real_Matrix(1..nx*ny*nz*nt,1..nx*ny*nz*nt); --((i-1)*ny*nz*nt+(ii-1)*nz*nt+(iii-1)*nt+iiii, -- (j-1)*ny*nz*nt+(jj-1)*nz*nt+(jjj-1)*nt+jjjj) xg : Real_Vector(1..nx); yg : Real_Vector(1..ny); zg : Real_Vector(1..nz); tg : Real_Vector(1..nt); -- (i-1)*ny*nz*nt+(ii-1)*nz*nt+(iii-1)*nt+iiii fg : Real_Vector(1..nx*ny*nz*nt); ug : Real_Vector(1..nx*ny*nz*nt); Ua : Real_Vector(1..nx*ny*nz*nt); xmin : real := 0.0; xmax : real := 1.0; ymin : real := 0.0; ymax : real := 1.0; zmin : real := 0.0; zmax : real := 1.0; tmin : real := 0.0; tmax : real := 1.0; x, y, z, t, hx, hy, hz, ht : real; npx : Integer := 12; npy : Integer := 12; npz : Integer := 12; npt : 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); tt : Real_Vector(1..npt); wt : Real_Vector(1..npt); val, err, avgerr, maxerr : real; tstart, tnow, tlast : Duration; -- in seconds tclock : Time; function s(i:Integer; ii:Integer; iii:Integer; iiii:Integer) return Integer is -- general subscript begin return (i-1)*ny*nz*nt+(ii-1)*nz*nt+(iii-1)*nt+iiii; end s; function FF(x,y,z,t:real) return real is -- forcing function begin return 706.0 + 120.0*t*t + 140.0*y*y*z + 768.0*x + 180.0*z*z*t + 200.0*y*y*z*t + 44.0*t*t*t + 110.0*y*y*z*z*t + 66.0*x*y + 12.0*t*t*t*t + 24.0*z*z*z*z + 36.0*Y*Y*Y*Y + 48.0*x*x*x*x + 60.0*y*y*z*z*t*t + 72.0*x*y*t + 84.0*x*z + 96.0*t ; end FF; function uana(x,y,z,t:real) return real is -- analytic solution and boundary begin return t**4 + 2.0*z**4 + 3.0*y**4 + 4.0*x**4 + 5.0*y*y*z*z*t*t + 6.0*x*y*t + 7.0*x*z + 8.0*t + 9.0; end uana; function galk(x,y,z,t:Real; i,j,ii,jj,iii,jjj,iiii,jjjj:integer) return real is -- Galerkin k stiffness function begin -- for this specific problem return (phipppp(x,j,nx,xg)* phi(y,jj,ny,yg)* phi(z,jjj,nz,zg)* phi(t,jjjj,nt,tg)+ 2.0* phi(x,j,nx,xg)*phipppp(y,jj,ny,yg)* phi(z,jjj,nz,zg)* phi(t,jjjj,nt,tg)+ 3.0* phi(x,j,nx,xg)* phi(y,jj,ny,yg)*phipppp(z,jjj,nz,zg)* phi(t,jjjj,nt,tg)+ 4.0* phi(x,j,nx,xg)* phi(y,jj,ny,yg)* phi(z,jjj,nz,zg)*phipppp(t,jjjj,nt,tg)+ 5.0* phip(x,j,nx,xg)* phip(y,jj,ny,yg)* phi(z,jjj,nz,zg)* phip(t,jjjj,nt,tg)+ 6.0* phi(x,j,nx,xg)* phipp(y,jj,ny,yg)* phipp(z,jjj,nz,zg)* phi(t,jjjj,nt,tg)+ 7.0* phi(x,j,nx,xg)* phi(y,jj,ny,yg)* phip(z,jjj,nz,zg)* phipp(t,jjjj,nt,tg)+ 8.0* phippp(x,j,nx,xg)* phi(y,jj,ny,yg)* phi(z,jjj,nz,zg)* phi(t,jjjj,nt,tg)+ 9.0* phi(x,j,nx,xg)* phipp(y,jj,ny,yg)* phi(z,jjj,nz,zg)* phip(t,jjjj,nt,tg)+ 10.0* phi(x,j,nx,xg)* phi(y,jj,ny,yg)* phip(z,jjj,nz,zg)* phip(t,jjjj,nt,tg)+ 11.0* phi(x,j,nx,xg)* phi(y,jj,ny,yg)* phi(z,jjj,nz,zg)* phip(t,jjjj,nt,tg)+ 12.0* phi(x,j,nx,xg)* phi(y,jj,ny,yg)* phi(z,jjj,nz,zg)* phi(t,jjjj,nt,tg))* (phi(x,i,nx,xg)* phi(y,ii,ny,yg)* phi(z,iii,nz,zg)* phi(t,iiii,nt,tg)); end galk; function galf(x,y,z,t:Real; i,ii,iii,iiii:integer) return real is -- Galerkin f function begin return FF(x,y,z,t)*phi(x,i,nx,xg)*phi(y,ii,ny,yg)*phi(z,iii,nz,zg)* phi(t,iiii,nt,tg); end galf; package Real_IO is new Float_IO(Real); use Real_IO; procedure pbound(i,ii,iii,iiii:Integer; x,y,z,t,v:Real) is begin 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(", iiii="&Integer'Image(iiii)); Put(", t="); Put(t,2,3,0); Put(" is "); Put(v); New_Line; end pbound; begin Put_Line("fem_check44_la.adb running "); tclock := Clock; tstart := Seconds(tclock); Put_Line("Solve uxxxx(x,y,z,t) +2 uyyyy(x,y,z,t) +3 uzzzz(x,y,z,t) +"); Put_Line(" 4 utttt(x,y,z,t) +5 uxyt(x,y,z,t) +6 uyyzz(x,y,z,t) +"); Put_Line(" 7 uztt(x,y,z,t) +8 uxxx(x,y,z,t) +9 uyyt(x,y,z,t) +"); Put_Line(" 10 uzt(x,y,z,t) +11 ut(x,y,z,t) +12 u(x,y,z,t) ="); Put_Line(" FF(x,y,z,t)"); Put_Line("FF(x,y,z,t) 706.0 + 120.0*t*t + 140.0*y*y*Z + 768.0*x + 180.0*z*z +"); Put_Line(" 200.0*y*y*z*t + 44.0*t*t*t + 110.0*y*y*z*z*t + 66.0*x*y +"); Put_Line(" 12.0*t*t*t*t + 24.0*z*z*z*z + 36.0*Y*Y*Y*Y + 48.0*x*x*x*x +"); Put_Line(" 60.0*y*y*z*z*t*t + 72.0*x*y*t + 84.0*x*z + 96.0*t "); Put_Line("boundary conditions computed using u(x,y,z,t)"); Put_Line("analytic solution given by u(x,y,z,t) = t^4 + 2 z^4 + 3 y^4 + "); Put_Line(" 4 x^4 + 5 y^2 z^2 t^2 + 6 x y t + 7 x z + 8 t + 9 "); 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); Put(", tmin="); Put(tmin,2,2,0); Put(", tmax="); Put(tmax,2,2,0); New_Line; Put_Line("nx="&Integer'Image(nx)&", ny="&Integer'Image(ny)& ", nz="&Integer'Image(nz)&", nt="&Integer'Image(nt)); -- 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,tmin "); for i in 1..nx loop xg(i) := xmin+Real(i-1)*hx; Ua(i) := uana(xg(i),ymin,Zmin,tmin); 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(tmin,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,tmin"); for ii in 1..ny loop yg(ii) := ymin+Real(ii-1)*hy; Ua(ii) := uana(xmin,yg(ii),zmin,tmin); 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(tmin,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, tmin"); for iii in 1..nz loop zg(iii) := zmin+Real(iii-1)*hz; Ua(iii) := uana(xmin,ymin,zg(iii),tmin); 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(tmin,2,5,0); Put(")="); Put(Ua(iii),3,5,0); New_Line; end loop; -- iii New_Line; -- initialize tg array. does not have to be uniform spacing ht := (tmax-tmin)/real(nt-1); Put_Line("t grid and analytic solution at xmin, ymin, zmin"); for iiii in 1..nt loop tg(iiii) := tmin+Real(iiii-1)*ht; Ua(iiii) := uana(xmin,ymin,zmin,tg(iiii)); Put("iiii="&Integer'Image(iiii)&", Ua("); Put(xmin,2,5,0); Put(","); Put(ymin,2,5,0); Put(","); Put(zmin,2,5,0); Put(","); Put(tg(iiii),2,5,0); Put(")="); Put(Ua(iiii),3,5,0); New_Line; end loop; -- 1iii 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; for iiii in 1..nt loop t := tmin+Real(iiii-1)*ht; Ua(s(i,ii,iii,iiii)) := uana(x,y,z,t); -- 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(", iiii="&Integer'Image(iiii)); -- Put(", t="); Put(t,2,3,0); -- Put(" is "); Put(Ua(s(i,ii,iii,iiii)),3,3,0); -- New_Line; end loop; -- iiii 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 for iiii in 1..nt loop for jjjj in 1..nt loop kg(s(i,ii,iii,iiii),s(j,jj,jjj,jjjj)) := 0.0; end loop; -- jjjj end loop; -- iiii end loop; -- jjj end loop; -- iii end loop; -- jj end loop; -- ii end loop; -- j end loop; -- i -- set boundary conditions, kg diagonal and fg sides, overkill for i in 1..nx loop for ii in 1..ny loop for iii in 1..nz loop for iiii in 1..nt loop x := xg(i); y := yg(ii); z := zg(iii); t := tg(iiii); kg(s(i,ii,iii,iiii),s(i,ii,iii,iiii)) := 1.0; val := uana(x,y,z,t); fg(s(i,ii,iii,iiii)) := val; end loop; -- iiii end loop; -- iii end loop; -- ii pbound(i,ny,nz,nt,x,y,z,t,val); end loop; -- i -- 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 Put_Line("calling gaulegf npt="&Integer'Image(npt)); gaulegf(tmin, tmax, tt, wt, npt); -- tt and wt set up for integrations Put_Line("xx(1)="&Real'Image(xx(1))&", xx(2)="&Real'Image(xx(2))); -- same index as C Put_Line("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))); Put_Line("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))); Put_Line("wz(1)="&Real'Image(wz(1))&", wz(2)="&Real'Image(wz(2))); Put_Line("tt(1)="&Real'Image(tt(1))&", tt(2)="&Real'Image(tt(2))); Put_Line("wt(1)="&Real'Image(wt(1))&", wt(2)="&Real'Image(wt(2))); Val := galk(xx(2),yy(2),zz(2),tt(2),2,2,2,2,2,2,2,2);-- i,j,etc. one higher than C Put_Line("galk=(xx(2),yy(2),zz(2),tt(2),2,2,2,2,2,2,2,2)="&Real'Image(val)); Val := galf(xx(2),yy(2),zz(2),tt(2),2,2,2,2); Put_Line("galf(xx(2),yy(2),zz(2),tt(2),2,2,2,2)="&Real'Image(val)); tclock := Clock; tnow := Seconds(tclock); Put_Line("initialization took "&Duration'Image(tnow-tstart)&" seconds"); -- 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 for iiii in 2..nt-1 loop for jjjj in 1..nt loop val := 0.0; for px in 1..npx loop for py in 1..npy loop for pz in 1..npz loop for pt in 1..npt loop val := val + wx(px)*wy(py)*wz(pz)*wt(pt)* galk(xx(px),yy(py),zz(pz),tt(pt), i,j,ii,jj,iii,jjj,iiii,jjjj); end loop; --pt end loop; --pz end loop; --py end loop; --px kg(s(i,ii,iii,iiii),s(j,jj,jjj,jjjj)) := val; end loop; -- jjjj end loop; -- iiii end loop; -- jjj end loop; -- iii end loop; -- jj end loop; -- ii Put_Line("Legendre integration="&Real'Image(val)& ", at i="&Integer'Image(i)&", j="&Integer'Image(j)& ", ii="&Integer'Image(ny-1)&", jj="&Integer'Image(ny)& ", iii="&Integer'Image(nz-1)&", jjj="&Integer'Image(nz)& ", iiii="&Integer'Image(nt-1)&", jjjj="&Integer'Image(nt)); end loop; -- j end loop; -- i -- compute integral for f(i), not changing boundary for i in 2..nx-1 loop for ii in 2..ny-1 loop for iii in 2..nz-1 loop for iiii in 2..nt-1 loop val := 0.0; for px in 1..npx loop for py in 1..npy loop for pz in 1..npx loop for pt in 1..npt loop val := val + wx(px)*wy(py)*wz(pz)*wt(pt)* galf(xx(px),yy(py),zz(pz),tt(pt),i,ii,iii,iiii); end loop; -- pt end loop; -- pz end loop; -- py end loop; -- px fg(s(i,ii,iii,iiii)) := val; end loop; -- iiii end loop; -- iii end loop; -- ii Put_Line("Legendre integration="&Real'Image(val)& ", F at i="&Integer'Image(i)&", ii="&Integer'Image(ny-1)& ", iii="&Integer'Image(nz-1)&", iiii="&Integer'Image(nt-1)); 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*nt)&" in "& Integer'Image(nx*ny*nz*nt)&"unknowns"); tclock := Clock; tlast := tnow; tnow := Seconds(tclock); Put_Line("integration took "&Duration'Image(tnow-tlast)&" seconds"); ug := simeq(kg, fg); tclock := Clock; tlast := tnow; tnow := Seconds(tclock); Put_Line("simeq took "&Duration'Image(tnow-tlast)&" seconds"); -- extra check, not part of solution declare nxyzt : Integer := nx*ny*nz*nt; fgck : Real_Vector(1..nxyzt); begin maxerr := 0.0; for i in 1..nxyzt loop fgck(i) := 0.0; for j in 1..nxyzt loop fgck(i) := fgck(i) + kg(i,j)*ug(j); end loop; err := abs(fgck(i)-fg(i)); if err>maxerr then maxerr := err; end if; end loop; Put_Line("simeq check maxerr="&Real'Image(Maxerr)); tclock := Clock; tlast := tnow; tnow := Seconds(tclock); Put_Line("extra fgck took "&Duration'Image(tnow-tlast)&" seconds"); end; -- extra check 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 for iiii in 1..nt loop err := abs(ug(s(i,ii,iii,iiii))-Ua(s(i,ii,iii,iiii))); if err>maxerr then maxerr := err; end if; avgerr := avgerr + err; Put("ug("&Integer'Image(i)&","&Integer'Image(ii)&","& Integer'Image(iii)&","&Integer'Image(iiii)&")="); Put(ug(s(i,ii,iii,iiii)),3,5,0); Put(", Ua="); Put(Ua(s(i,ii,iii,iiii)),3,5,0); Put(", err="); Put(err,3,5,0); New_Line; end loop; -- iiii end loop; -- iii end loop; -- ii end loop; -- i Put_Line(" maxerr="& Real'Image(maxerr)&", avgerr="& Real'image(avgerr/Real(nx*ny*nz*nt))); tclock := Clock; tnow := Seconds(tclock); Put_Line("total run took "&Duration'Image(tnow-tstart)&" seconds"); Put_Line("end fem_check44_la.adb"); end fem_check44_la;