-- pde48h_eq.adb Biharmonic 4th order PDE in seven dimensions -- sample problem: -- solve Uxxxx(x,y,z,t,u,v,w,p) + Uyyyy(x,y,z,t,u,v,w,p) + Uzzzz(x,y,z,t,u,v,w,p) + -- Utttt(x,y,z,t,u,v,w,p) + Uuuuu(x,y,z,t,u,v,w,p) + Uvvvv(x,y,z,t,u,v,w,p) + -- Uwwww(x,y,z,t,u,v,w,p) + Upppp(x,y,z,t,u,v,w,p) + -- 2 Uxx(x,y,z,t,u,v,w,p) + 2 Uyy(x,y,z,t,u,v,w,p) + 2 Uzz(x,y,z,t,u,v,w,p) + -- 2 Utt(x,y,z,t,u,v,w,p) + 2 Uuu(x,y,z,t,u,v,w,p) + 2 Uvv(x,y,z,t,u,v,w,p) + -- 2 Uww(x,y,z,t,u,v,w,p) + 2 Upp(x,y,z,t,u,v,w,p) + -- 8 U(x,y,z,t,u,v,w,p) = f(x,y,z,t,u,v,w,p) -- -- boundary conditions computed using Ub(x,y,z,t,u,v,w,p) -- analytic solution is U(x,y,z,t,u,v,w,p) = sin(x+y+z+t+u+v+w+p) with ada.text_io; use ada.text_io; with ada.numerics; use ada.numerics; with ada.numerics.Long_elementary_Functions; use ada.numerics.Long_elementary_Functions; with real_arrays; use real_arrays; -- defines type real with simeq; with rderiv; with ada.calendar; use ada.Calendar; procedure pde48h_eq is nx : integer := 5; -- problem parameters ny : integer := 5; nz : integer := 5; nt : integer := 5; nu : integer := 5; nv : integer := 5; nw : integer := 5; np : integer := 5; nxyztuvwp : integer := (nx-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2)*(np-2); A: real_matrix(1..nxyztuvwp,1..nxyztuvwp); -- A * U = R R : real_vector(1..nxyztuvwp); -- RHS known UU : real_vector(1..nxyztuvwp); -- solution being computed - not boundary Ua : real_vector(1..nxyztuvwp); -- solution for checking - not boundary xg : real_vector(1..nx); yg : real_vector(1..ny); zg : real_vector(1..nz); tg : real_vector(1..nt); ug : real_vector(1..nu); vg : real_vector(1..nv); wg : real_vector(1..nw); pg : real_vector(1..np); xmin : real := 0.0; -- problem parameters xmax : real := 0.5; ymin : real := 0.0; ymax : real := 0.5; zmin : real := 0.0; zmax : real := 0.5; tmin : real := 0.0; tmax : real := 0.5; umin : real := 0.0; umax : real := 0.5; vmin : real := 0.0; vmax : real := 0.5; wmin : real := 0.0; wmax : real := 0.5; pmin : real := 0.0; pmax : real := 0.5; hx, hy, hz, ht, hu, hv, Hw, hp : real; package real_io is new float_io(real); use real_io; -- internal functions function S(i:integer; ii:integer; i3:integer; i4:integer; i5:integer; i6:integer; i7:integer; i8:integer) return integer is -- for internal x,y,z,t,u,v,w,p begin return (i- 2)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2)*(np-2) + (ii-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2)*(np-2) + (i3-2)*(Nt-2)*(nu-2)*(nv-2)*(nw-2)*(np-2) + (i4-2)*(nu-2)*(nv-2)*(nw-2)*(np-2) + (i5-2)*(nv-2)*(nw-2)*(np-2) + (i6-2)*(nw-2)*(np-2) + (i7-2)*(np-2) + (i8-1); end S; -- problem RHS, function f(x:real; y:real; z:real; t:real; u:real; v:real; w:real; p:real) return real is begin return sin(x+y+z+t+u+v+w+p); end f; -- problem boundaries (and analytic solution for checking) function Ub(x:real; y:real; z:real; t:real; u:real; v:real; w:Real; p:real) return real is begin return sin(x+y+z+t+u+v+w+p); end Ub; function eval_derivative(xord:integer; yord:integer; zord:integer; tord:integer; uord:integer; vord:integer; word:integer; pord:integer; I :integer; ii:integer; i3:integer; i4:integer; i5:integer; i6:integer; i7:Integer; i8:integer ; UU:real_vector) return real is cx : real_vector(1..nx); -- can be used for any order cy : real_vector(1..ny); cz : real_vector(1..nz); ct : real_vector(1..nt); cu : real_vector(1..nu); cv : real_vector(1..nv); cw : real_vector(1..nw); cp : real_vector(1..np); discrete : real; x, y, z, t, u, v, W, p : real; begin x := xg(i); y := yg(ii); z := zg(i3); t := tg(i4); u := ug(i5); v := vg(i6); w := wg(i7); p := pg(i8); discrete := 0.0; -- eight cases of one partial x, y, z, t, u, v, w, p to any order if xord/=0 and yord=0 and zord=0 and tord=0 and uord=0 and vord=0 and word=0 and pord=0 then rderiv(xord, nx, i, hx, cx); for j in 1..nx loop if j=1 or j=nx then discrete := discrete + cx(j)*Ub(xg(j),y,z,t,u,v,w,p); else discrete := discrete + cx(j)*UU(s(j,ii,i3,i4,i5, i6,I7,i8)); end if; end loop; -- j x elsif xord=0 and yord/=0 and zord=0 and tord=0 and uord=0 and vord=0 and word=0 and pord=0 then rderiv(yord, ny, ii, hy, cy); for jj in 1..ny loop if jj=1 or jj=ny then discrete := discrete + cy(jj)*Ub(x,yg(jj),z,t,u,v,w,p); else discrete := discrete + cy(jj)*UU(s(i,jj,i3,i4,i5, i6,I7,i8)); end if; end loop; -- jj y elsif xord=0 and yord=0 and zord/=0 and tord=0 and uord=0 and vord=0 and word=0 and pord=0 then rderiv(zord, nz, i3, hz, cz); for j3 in 1..nz loop if j3=1 or j3=nz then discrete := discrete + cz(j3)*Ub(x,y,zg(j3),t,u,v,w,p); else discrete := discrete + cz(j3)*UU(s(i,ii,j3,i4,i5, i6,I7,i8)); end if; end loop; -- j3 z elsif xord=0 and yord=0 and zord=0 and tord/=0 and uord=0 and vord=0 and word=0 and pord=0 then rderiv(tord, nt, i4, ht, ct); for j4 in 1..nt loop if j4=1 or j4=nt then discrete := discrete + ct(j4)*Ub(x,y,z,tg(j4),u,v,w,p); else discrete := discrete + ct(j4)*UU(s(i,ii,i3,j4,i5, i6,I7,i8)); end if; end loop; -- j4 t elsif xord=0 and yord=0 and zord=0 and tord=0 and uord/=0 and vord=0 and word=0 and pord=0 then rderiv(uord, nu, i5, hu, cu); for j5 in 1..nu loop if j5=1 or j5=nu then discrete := discrete + cu(j5)*Ub(x,y,z,t,ug(j5),v,w,p); else discrete := discrete + cu(j5)*UU(s(i,ii,i3,i4,j5, i6,i7,i8)); end if; end loop; -- j5 u elsif xord=0 and yord=0 and zord=0 and tord=0 and uord=0 and vord/=0 and word=0 and pord=0 then rderiv(vord, nv, i6, hv, cv); for j6 in 1..nv loop if j6=1 or j6=nv then discrete := discrete + cv(j6)*Ub(x,y,z,t,u,vg(j6),w,p); else discrete := discrete + cv(j6)*UU(s(i,ii,i3,i4,i5, j6,i7,i8)); end if; end loop; -- j6 v elsif xord=0 and yord=0 and zord=0 and tord=0 and uord=0 and vord=0 and word/=0 and pord=0 then rderiv(word, nw, i7, hw, cw); for j7 in 1..nw loop if j7=1 or j7=nw then discrete := discrete + cw(j7)*Ub(x,y,z,t,u,v,wg(j7),p); else discrete := discrete + cw(j7)*UU(s(i,ii,i3,i4,i5, i6,j7,i8)); end if; end loop; -- j7 w elsif xord=0 and yord=0 and zord=0 and tord=0 and uord=0 and vord=0 and word=0 and Pord/=0 then rderiv(pord, np, i8, hp, cp); for j8 in 1..np loop if j8=1 or j8=np then discrete := discrete + cp(j8)*Ub(x,y,z,t,u,v,w,pg(j8)); else discrete := discrete + cp(j8)*UU(s(i,ii,i3,i4,i5, i6,i7,j8)); end if; end loop; -- j8 p end if; -- other cases not needed for this PDE return discrete; end eval_derivative; procedure check_soln(UUU:real_vector) is val, err, smaxerr : real; bi, bii, bi3, bi4, bi5, bi6, bi7, bi8 : integer; begin smaxerr := 0.0; for i in 2..nx-1 loop for ii in 2..ny-1 loop for i3 in 2..nz-1 loop for i4 in 2..nt-1 loop for i5 in 2..nu-1 loop for i6 in 2..nv-1 loop for i7 in 2..nw-1 loop for i8 in 2..np-1 loop val := 0.0; -- + Uxxxx(x,y,z,t,u,v,w,p) val := val + eval_derivative(4,0,0,0,0,0,0,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + Uyyyy(x,y,z,t,u,v,w,p) val := val + eval_derivative(0,4,0,0,0,0,0,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + Uzzzz(x,y,z,t,u,v,w,p) val := val + eval_derivative(0,0,4,0,0,0,0,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + Utttt(x,y,z,t,u,v,w,p) val := val + eval_derivative(0,0,0,4,0,0,0,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + Uuuuu(x,y,z,t,u,v,w,p) val := val + eval_derivative(0,0,0,0,4,0,0,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + Uvvvv(x,y,z,t,u,v,w,p) val := val + eval_derivative(0,0,0,0,0,4,0,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + Uwwww(x,y,z,t,u,v,w,p) val := val + eval_derivative(0,0,0,0,0,0,4,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + Upppp(x,y,z,t,u,v,w,p) val := val + eval_derivative(0,0,0,0,0,0,0,4, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + 2 Uxx(x,y,z,t,u,v,w,p) val := val + 2.0*eval_derivative(2,0,0,0,0,0,0,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + 2 Uyy(x,y,z,t,u,v,w,p) val := val + 2.0*eval_derivative(0,2,0,0,0,0,0,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + 2 Uzz(x,y,z,t,u,v,w,p) val := val + 2.0*eval_derivative(0,0,2,0,0,0,0,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + 2 Utt(x,y,z,t,u,v,w,p) val := val + 2.0*eval_derivative(0,0,0,2,0,0,0,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + 2 Uuu(x,y,z,t,u,v,w,p) val := val + 2.0*eval_derivative(0,0,0,0,2,0,0,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + 2 Uvv(x,y,z,t,u,v,w,p) val := val + 2.0*eval_derivative(0,0,0,0,0,2,0,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + 2 Uww(x,y,z,t,u,v,w,p) val := val + 2.0*eval_derivative(0,0,0,0,0,0,2,0, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + 2 Upp(x,y,z,t,u,v,w,p) val := val + 2.0*eval_derivative(0,0,0,0,0,0,0,2, i,ii,i3,i4,i5,i6,i7,i8,UUU); -- + 8 u(x,y,z,t,u,v,w,p) val := val + 8.0*Ub(xg(i),yg(ii),zg(i3),tg(i4),ug(i5), vg(i6),wg(i7),pg(i8)); -- - RHS val := val - f(xg(i),yg(ii),zg(i3),tg(i4),ug(i5), vg(i6),wg(i7),pg(i8)); if(i=nx-1 and ii=ny-1 and i3=nz-1 and i4=nt-1 and i5=nu-1 and i6=nv-1 and i7=nw-1 and i8=np-1) then put_line("UUU="&real'image(UUU(s(i,ii,i3,i4,i5,i6,i7,i8)))); put_line("p''''="&real'image(eval_derivative(0,0,0,0,0,0,0,4, i,ii,i3,i4,i5,i6,i7,i8,UUU))); put_line("2p''="&real'image(2.0*eval_derivative(0,0,0,0,0,0,0,2, i,ii,i3,i4,i5,i6,i7,i8,UUU))); put_line("8Ub="&real'image(8.0*Ub(xg(i),yg(ii),zg(i3),tg(i4),ug(i5), vg(i6),wg(i7),pg(i8)))); put_line("f="&real'image(f(xg(i),yg(ii),zg(i3),tg(i4),ug(i5), vg(i6),wg(i7),pg(i8)))); end if; -- should be zero err := abs(val); if err>smaxerr then smaxerr := err; bi := i; bii := ii; bi3 := i3; bi4 := i4; bi5 := i5; bi6 := i6; bi7 := i7; bi8 := i8; end if; end loop; --i8 end loop; --i7 end loop; --i6 end loop; --i5 end loop; --i4 end loop; --i3 end loop; -- ii end loop; -- i put_line("check_soln maxerr="&real'image(smaxerr)); put_line("at="&integer'image(bi)&integer'image(bii)& integer'image(bi4)&integer'image(bi4)& integer'image(bi5)&integer'image(bi6)& integer'image(bi7)&integer'image(bi8)); end Check_soln; procedure write_soln(Which:Integer; tag:String) is fout : File_Type; begin create(fout, out_file, "pde48h_eq_ada"&tag&".dat"); put_line("writing pde48h_eq_ada"&tag&".dat"); put_line(fout,"8 "&integer'image(nx)&integer'image(ny)& integer'image(nz)&integer'image(nt)& integer'image(nu)&integer'image(nv)& integer'image(nw)&integer'image(np)); for i in 1..nx loop for ii in 1..ny loop for i3 in 1..nz loop for i4 in 1..nt loop for i5 in 1..nu loop for i6 in 1..nv loop for i7 in 1..nw loop for i8 in 1..np loop put(fout," "); put(fout,xg(i),2,6,0); put(fout," "); put(fout,yg(ii),2,6,0); put(fout," "); put(fout,zg(i3),2,6,0); put(fout," "); put(fout,tg(i4),2,6,0); put(fout," "); put(fout,ug(i5),2,6,0); put(fout," "); put(fout,vg(i6),2,6,0); put(fout," "); put(fout,wg(i7),2,6,0); put(fout," "); put(fout,pg(i8),2,6,0); put(fout," "); if Which=1 then -- expected solution put(fout,Ub(xg(i),yg(ii),zg(i3),tg(i4), ug(i5),vg(i6),wg(i7),pg(i8)),4,6,0); else -- computed solution if i=1 or i=nx or ii=1 or ii=ny or i3=1 or i3=nz or i4=1 or i4=nt or i5=1 or i5=nu or i6=1 or i6=nv or i7=1 or i7=nw or i8=1 or i8=np then -- boundary put(fout,Ub(xg(i),yg(ii),zg(i3),tg(i4), ug(i5),vg(i6),wg(i7),pg(i8)),4,6,0); else put(fout,UU(S(i,ii,i3,i4,i5,i6,i7,i8)),4,6,0); end if; end if; new_line(fout); end loop; -- i8 end loop; -- i7 end loop; -- i6 end loop; -- i5 end loop; -- i4 end loop; -- i3 end loop; -- ii end loop; -- i close(fout); put_line("finished writing pde48h_eq_ada"&tag&".dat"); end write_soln; procedure initializeA is begin put_line("initialize A "); -- initialize A for i in 2..nx-1 loop for ii in 2..ny-1 loop for i3 in 2..nz-1 loop for i4 in 2..nt-1 loop for i5 in 2..nu-1 loop for i6 in 2..nv-1 loop for i7 in 2..nw-1 loop for i8 in 2..np-1 loop for j in 2..nx-1 loop for jj in 2..ny-1 loop for j3 in 2..nz-1 loop for j4 in 2..nt-1 loop for j5 in 2..nu-1 loop for j6 in 2..nv-1 loop for j7 in 2..nw-1 loop for j8 in 2..np-1 loop A(s(i,ii,i3,i4,i5,i6,i7,i8), s(j,jj,j3,j4,j5,j6,j7,j8)):=0.0; end loop; -- j8 end loop; -- j7 end loop; -- j6 end loop; -- j5 end loop; -- j4 end loop; -- j3 end loop; -- jj end loop; -- j R(s(i,ii,i3,i4,i5,i6,i7,i8)) := 0.0; end loop; -- i8 end loop; -- i7 end loop; -- i6 end loop; -- i5 end loop; -- i4 end loop; -- i3 end loop; -- ii end loop; -- i end initializeA; procedure buildA is -- system of equations val : real := 0.0; k: integer; cxx, cxxxx : real_vector(1..nx); cyy, cyyyy : real_vector(1..ny); czz, czzzz : real_vector(1..nz); ctt, ctttt : real_vector(1..nt); cuu, cuuuu : real_vector(1..nu); cvv, cvvvv : real_vector(1..nv); cww, cwwww : real_vector(1..nw); cpp, cpppp : real_vector(1..np); begin -- compute entries in A matrix, inside boundary put_line("compute A matrix and Y RHS "); put_line("check rderiv(4,nx,1,hx,cxxxx"); rderiv(4,nx,1,hx,cxxxx); for i in 1..nx loop put_line("cxxxx("&integer'image(i)& ")="&real'image(cxxxx(i))); end loop; put_line("check rderiv(4,nx,nx,hx,cxxxx"); rderiv(4,nx,nx,hx,cxxxx); for i in 1..nx loop put_line("cxxxx("&integer'image(i)& ")="&real'image(cxxxx(i))); end loop; put_line(" "); val := 0.0; for i in 2..nx-1 loop for ii in 2..ny-1 loop for i3 in 2..nz-1 loop for i4 in 2..nt-1 loop for i5 in 2..nu-1 loop for i6 in 2..nv-1 loop for i7 in 2..nw-1 loop for i8 in 2..np-1 loop k := s(i,ii,i3,i4,i5,i6,i7,i8); -- row index -- for each term in first equation -- for d^4U/dx^4 rderiv(4,nx,i,hx,cxxxx); for j in 1..nx loop val := cxxxx(j); if j=1 or j=nx then R(k) := R(k) - val*Ub(xg(j),yg(ii),zg(i3),tg(i4), ug(i5),vg(i6),wg(i7),pg(i8)); else A(k,s(j,ii,i3,i4,i5,i6,i7,i8)) := A(k,s(j,ii,i3,i4,i5,i6,i7,i8)) + val; end if; end loop; -- j x -- for d^4U/dy^4 rderiv(4,ny,ii,hy,cyyyy); for jj in 1..ny loop val := cyyyy(jj); if jj=1 or jj=ny then R(k) := R(k) - val*Ub(xg(i),yg(jj),zg(i3),tg(i4), ug(i5),vg(i6),wg(i7),pg(i8)); else A(k,s(i,jj,i3,i4,i5,i6,i7,i8)) := A(k,s(i,jj,i3,i4,i5,i6,i7,i8)) + val; end if; end loop; -- jj y -- for d^4U/dz^4 rderiv(4,nz,i3,hz,czzzz); for j3 in 1..nz loop val := czzzz(j3); if j3=1 or j3=nz then R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(j3),tg(i4), ug(i5),vg(i6),wg(i7),pg(i8)); else A(k,s(i,ii,j3,i4,i5,i6,i7,i8)) := A(k,s(i,ii,j3,i4,i5,i6,i7,i8)) + val; end if; end loop; -- j3 z -- for d^4U/dt^4 rderiv(4,nt,i4,ht,ctttt); for j4 in 1..nt loop val := ctttt(j4); if j4=1 or j4=nt then R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(i3),tg(j4), ug(i5),vg(i6),wg(i7),pg(i8)); else A(k,s(i,ii,i3,j4,i5,i6,i7,i8)) := A(k,s(i,ii,i3,j4,i5,i6,i7,i8)) + val; end if; end loop; -- j4 t -- for d^4U/du^4 rderiv(4,nu,i5,hu,cuuuu); for j5 in 1..nu loop val := cuuuu(j5); if j5=1 or j5=nu then R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(i3),tg(i4), ug(j5),vg(i6),wg(i7),pg(i8)); else A(k,s(i,ii,i3,i4,j5,i6,i7,i8)) := A(k,s(i,ii,i3,i4,j5,i6,i7,i8)) + val; end if; end loop; -- j5 u -- for d^4U/dv^4 rderiv(4,nv,i6,hv,cvvvv); for j6 in 1..nv loop val := cvvvv(j6); if j6=1 or j6=nv then R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(i3),tg(i4), ug(i5),vg(j6),wg(i7),pg(i8)); else A(k,s(i,ii,i3,i4,i5,j6,i7,i8)) := A(k,s(i,ii,i3,i4,i5,j6,i7,i8)) + val; end if; end loop; -- j6 v -- for d^4U/dw^4 rderiv(4,nw,i7,hw,cwwww); for j7 in 1..nw loop val := cwwww(j7); if j7=1 or j7=nw then R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(i3),tg(i4), ug(i5),vg(i6),wg(j7),pg(i8)); else A(k,s(i,ii,i3,i4,i5,i6,j7,i8)) := A(k,s(i,ii,i3,i4,i5,i6,j7,i8)) + val; end if; end loop; -- j7 w -- for d^4U/dp^4 rderiv(4,np,i8,hp,cpppp); for j8 in 1..np loop val := cpppp(j8); if j8=1 or j8=np then R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(i3),tg(i4), ug(i5),vg(i6),wg(i7),pg(j8)); else A(k,s(i,ii,i3,i4,i5,i6,i7,j8)) := A(k,s(i,ii,i3,i4,i5,i6,i7,j8)) + val; end if; end loop; -- j8 p -- for d^2U/dx^2 rderiv(2,nx,i,hx,cxx); for j in 1..nx loop val := 2.0*cxx(j); if j=1 or j=nx then R(k) := R(k) - val*Ub(xg(j),yg(ii),zg(i3),tg(i4), ug(i5),vg(i6),wg(i7),pg(i8)); else A(k,s(j,ii,i3,i4,i5,i6,i7,i8)) := A(k,s(j,ii,i3,i4,i5,i6,i7,i8)) + val; end if; end loop; -- j x -- for d^2U/dy^2 rderiv(2,ny,ii,hy,cyy); for jj in 1..ny loop val := 2.0*cyy(jj); if jj=1 or jj=ny then R(k) := R(k) - val*Ub(xg(i),yg(jj),zg(i3),tg(i4), ug(i5),vg(i6),wg(i7),pg(i8)); else A(k,s(i,jj,i3,i4,i5,i6,i7,i8)) := A(k,s(i,jj,i3,i4,i5,i6,i7,i8)) + val; end if; end loop; -- jj y -- for d^2U/dz^2 rderiv(2,nz,i3,hz,czz); for j3 in 1..nz loop val := 2.0*czz(j3); if j3=1 or j3=nz then R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(j3),tg(i4), ug(i5),vg(i6),wg(i7),pg(i8)); else A(k,s(i,ii,j3,i4,i5,i6,i7,i8)) := A(k,s(i,ii,j3,i4,i5,i6,i7,i8)) + val; end if; end loop; -- j3 z -- for d^2U/dt^2 rderiv(2,nt,i4,ht,ctt); for j4 in 1..nt loop val := 2.0*ctt(j4); if j4=1 or j4=nt then R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(i3),tg(j4), ug(i5),vg(i6),wg(i7),pg(i8)); else A(k,s(i,ii,i3,j4,i5,i6,i7,i8)) := A(k,s(i,ii,i3,j4,i5,i6,i7,i8)) + val; end if; end loop; -- j4 t -- for d^2U/du^2 rderiv(2,nu,i5,hu,cuu); for j5 in 1..nu loop val := 2.0*cuu(j5); if j5=1 or j5=nu then R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(i3),tg(i4), ug(j5),vg(i6),wg(i7),pg(i8)); else A(k,s(i,ii,i3,i4,j5,i6,i7,i8)) := A(k,s(i,ii,i3,i4,j5,i6,i7,i8)) + val; end if; end loop; -- j5 u -- for d^2U/dv^2 rderiv(2,nv,i6,hv,cvv); for j6 in 1..nv loop val := 2.0*cvv(j6); if j6=1 or j6=nv then R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(i3),tg(i4), ug(i5),vg(j6),wg(i7),pg(i8)); else A(k,s(i,ii,i3,i4,i5,j6,i7,i8)) := A(k,s(i,ii,i3,i4,i5,j6,i7,i8)) + val; end if; end loop; -- j6 v -- for d^2U/dw^2 rderiv(2,nw,i7,hw,cww); for j7 in 1..nw loop val := 2.0*cww(j7); if j7=1 or j7=nw then R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(i3),tg(i4), ug(i5),vg(i6),wg(j7),pg(i8)); else A(k,s(i,ii,i3,i4,i5,i6,j7,i8)) := A(k,s(i,ii,i3,i4,i5,i6,j7,i8)) + val; end if; end loop; -- j7 w -- for d^2U/dp^2 rderiv(2,np,i8,hp,cpp); for j8 in 1..np loop val := 2.0*cpp(j8); if j8=1 or j8=np then R(k) := R(k) - val*Ub(xg(i),yg(ii),zg(i3),tg(i4), ug(i5),vg(i6),wg(i7),pg(j8)); else A(k,s(i,ii,i3,i4,i5,i6,i7,j8)) := A(k,s(i,ii,i3,i4,i5,i6,i7,j8)) + val; end if; end loop; -- j8 p -- for U(x,y,z,t,u,v,w,p) A(k,s(i,ii,i3,i4,i5,i6,i7,i8)) := A(k,s(i,ii,i3,i4,i5,i6,i7,i8)) + 8.0; -- for f(x,y,z,t,u,v,w,p) RHS constant R(k) := R(k) + f(xg(i),yg(ii),zg(i3),tg(i4), ug(i5),vg(i6),wg(i7),pg(i8)); -- end terms end loop; -- i8 end loop; -- i7 end loop; -- i6 end loop; -- i5 end loop; -- i4 end loop; -- i3 end loop; -- ii end loop; -- i end BuildA; procedure main is x, y, z, t, u, v, w, p : real; err, avgerr, maxerr : real; time_start : duration; time_now : duration; time_last : duration; begin -- main put_line("pde48h_eq.adb running "); put_line("Given Uxxxx(x,y,z,t,u,v,w,p) + Uyyyy(x,y,z,t,u,v,w,p) + Uzzzz(x,y,z,t,u,v,w,p) + "); put_line(" Utttt(x,y,z,t,u,v,w,p) + Uuuuu(x,y,z,t,u,v,w,p) + Uvvvv(x,y,z,t,u,v,w,p) + "); put_line(" Uwwww(x,y,z,t,u,v,w,p) + Upppp(x,y,z,t,u,v,w,p) + "); put_line(" 2 Uxx(x,y,z,t,u,v,w,p) + 2 Uyy(x,y,z,t,u,v,w,p) + 2 Uzz(x,y,z,t,u,v,w,p) + "); put_line(" 2 Utt(x,y,z,t,u,v,w,p) + 2 Uuu(x,y,z,t,u,v,w,p) + 2 Uvv(x,y,z,t,u,v,w,p) + "); put_line(" 2 Uww(x,y,z,t,u,v,w,p) + 2 Upp(x,y,z,t,u,v,w,p) + "); put_line(" 8 U(x,y,z,t,u,v,w,p) = f(x,y,z,t,u,v,w,p) "); put_line("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries "); put_line("tmin<=t<=tmax umin<=u<=umax vmin<=v<=vmax Boundaries "); put_line("wmin<=w<=wmax pmin<=p<=pmax Boundaries "); put_line("Analytic solution U(x,y,z,t,u,v,w,p) = sin(x+y+z+t+u+v+w+p)"); put_line(" "); put_line("xmin="&real'image(xmin)&", xmax="&real'image(xmax)); put_line("ymin="&real'image(ymin)&", ymax="&real'image(ymax)); put_line("zmin="&real'image(zmin)&", zmax="&real'image(zmax)); put_line("tmin="&real'image(tmin)&", tmax="&real'image(tmax)); put_line("umin="&real'image(umin)&", umax="&real'image(umax)); put_line("vmin="&real'image(vmin)&", vmax="&real'image(vmax)); put_line("wmin="&real'image(wmin)&", wmax="&real'image(wmax)); put_line("pmin="&real'image(pmin)&", pmax="&real'image(pmax)); put_line("nx="&integer'image(nx)&", ny="&integer'image(ny)& ", nz="&integer'image(nz)&", nt="&integer'image(nt)& ", nu="&integer'image(nu)&", nv="&integer'image(nv)& ", nw="&integer'image(nw)&", np="&integer'image(np)); put_line(" "); time_start := seconds(clock); time_last := time_start; hx := (xmax-xmin)/real(nx-1); put_line("x divisor, hx="&real'image(hx)); xg(1) := xmin; put_line("xg( 1)="&real'image(xg(1))); for i in 2..nx loop -- decreasing step size xg(i) := xmin+real(i-1)*hx; put_line("xg("&integer'image(i)&")="&real'image(xg(i))); end loop; -- i x put_line(" "); hy := (ymax-ymin)/real(ny-1); put_line("y divisor, hy="&real'image(hy)); yg(1) := ymin; put_line("yg( 1)="&real'image(yg(1))); for ii in 2..ny loop -- decreasing step size yg(ii) := ymin+real(ii-1)*hy; put_line("yg("&integer'image(ii)&")="&real'image(yg(ii))); end loop; -- ii y put_line(" "); hz := (zmax-zmin)/real(nz-1); put_line("z divisor, hz="&real'image(hz)); zg(1) := zmin; put_line("zg( 1)="&real'image(zg(1))); for i3 in 2..nz loop -- decreasing step size zg(i3) := zmin+real(i3-1)*hz; put_line("zg("&integer'image(i3)&")="&real'image(zg(i3))); end loop; -- i3 z put_line(" "); ht := (tmax-tmin)/real(nt-1); put_line("t divisor, hz="&real'image(ht)); tg(1) := tmin; put_line("tg( 1)="&real'image(tg(1))); for i4 in 2..nt loop -- decreasing step size tg(i4) := tmin+real(i4-1)*ht; put_line("tg("&integer'image(i4)&")="&real'image(tg(i4))); end loop; -- i4 t put_line(" "); hu := (umax-umin)/real(nu-1); put_line("u divisor, hu="&real'image(hu)); ug(1) := umin; put_line("ug( 1)="&real'image(ug(1))); for i5 in 2..nu loop -- decreasing step size ug(i5) := umin+real(i5-1)*hu; put_line("ug("&integer'image(i5)&")="&real'image(ug(i5))); end loop; -- i5 u put_line(" "); hv := (vmax-vmin)/real(nv-1); put_line("v divisor, hv="&real'image(hv)); vg(1) := vmin; put_line("vg( 1)="&real'image(vg(1))); for i6 in 2..nv loop -- decreasing step size vg(i6) := vmin+real(i6-1)*hv; put_line("vg("&integer'image(i6)&")="&real'image(vg(i6))); end loop; -- i6 v put_line(" "); hw := (wmax-wmin)/real(nw-1); put_line("w divisor, hw="&real'image(hw)); wg(1) := wmin; put_line("wg( 1)="&real'image(wg(1))); for i7 in 2..nw loop -- decreasing step size wg(i7) := wmin+real(i7-1)*hw; put_line("wg("&integer'image(i7)&")="&real'image(wg(i7))); end loop; -- i7 w put_line(" "); hp := (pmax-pmin)/real(np-1); put_line("p divisor, hp="&real'image(hp)); pg(1) := pmin; put_line("pg( 1)="&real'image(pg(1))); for i8 in 2..np loop -- decreasing step size pg(i8) := pmin+real(i8-1)*hp; put_line("pg("&integer'image(i8)&")="&real'image(pg(i8))); end loop; -- i8 p put_line(" "); put_line("put solution in Ua vector "); -- put solution in Ua vector for i in 2..nx-1 loop x := xg(i); for ii in 2..ny-1 loop y := yg(ii); for i3 in 2..nz-1 loop z := zg(i3); for i4 in 2..nt-1 loop t := tg(i4); for i5 in 2..nu-1 loop u := ug(i5); for i6 in 2..nv-1 loop v := vg(i6); for i7 in 2..nw-1 loop w := wg(i7); for i8 in 2..np-1 loop p := pg(i8); Ua(s(i,ii,i3,i4,i5,i6,i7,i8)) := Ub(x,y,z,t,u,v,w,p); end loop; -- i8 end loop; -- i7 end loop; -- i6 end loop; -- i5 end loop; -- i4 end loop; -- i3 end loop; -- ii end loop; -- i time_now := seconds(clock); put_line("time for previous section = "&Duration'image(time_now-time_last)& " seconds "); time_last := time_now; put_line(" "); initializeA; time_now := seconds(clock); put_line("time for previous section = "&Duration'image(time_now-time_last)& " seconds "); time_last := time_now; put_line(" "); buildA; -- system of equations put_line("A computed stiffness matrix "); time_now := seconds(clock); put_line("time for previous section = "&Duration'image(time_now-time_last)& " seconds "); time_last := time_now; put_line(" "); put_line("solve "&integer'image(nxyztuvwp)&" equations in "& integer'image(nxyztuvwp)&" unknowns "); UU := simeq(A, R); put_line("equations solved "); time_now := seconds(clock); put_line("time for previous section = "&Duration'image(time_now-time_last)& " seconds "); time_last := time_now; put_line(" "); put_line("check pde against known solution"); check_soln(Ua); put_line("check pde against computed solution"); check_soln(UU); put_line(" "); avgerr := 0.0; maxerr := 0.0; put_line(" U computed, Ua analytic, error "); for i in 2..nx-1 loop for ii in 2..ny-1 loop for i3 in 2..nz-1 loop for i4 in 2..nt-1 loop for i5 in 2..nu-1 loop for i6 in 2..nv-1 loop for i7 in 2..nw-1 loop for i8 in 2..np-1 loop err := abs(UU(s(i,ii,i3,i4,i5,i6,i7,i8))- Ua(s(i,ii,i3,i4,i5,i6,i7,i8))); if err>maxerr then maxerr := err; end if; avgerr := avgerr + err; put("ug("&integer'image(i)&","&integer'image(ii)& ","&integer'image(i3)&","&integer'image(i4)& ","&integer'image(i5)&","&integer'image(i6)& ","&integer'image(i7)&","&Integer'Image(i8)&")="); put(UU(s(i,ii,i3,i4,i5,i6,i7,i8)),3,6); put(", Ua="); put(Ua(s(i,ii,i3,i4,i5,i6,i7,i8)),3,6); put(", err="); put(err); new_line; end loop; -- i8 end loop; -- i7 end loop; -- i6 end loop; -- i5 end loop; -- i4 end loop; -- i3 end loop; -- ii end loop; -- i time_now := seconds(clock); put_line("total CPU time = "&Duration'image(time_now-time_start)& " seconds"); put_line(" "); put_line(" maxerr="&real'image(maxerr)& ", avgerr="&real'image(avgerr/real(nxyztuvwp))); put_line(" "); end Main; begin -- pde48h_eq main; -- write_soln(1,"U"); -- analytic solution -- write_soln(2,""); -- computed solution for plot8d put_line("end pde48h_eq.adb"); end pde48h_eq;