-- pde46h_nu.adb homogeneous Biharmonic 4th order PDE in six dimensions -- solve Uvvvv(v,w,x,y,z,t) + Uwwww(v,w,x,y,z,t) + Uxxxx(v,w,x,y,z,t) + -- Uyyyy(v,w,x,y,z,t) + Uzzzz(v,w,x,y,z,t) + Utttt(v,w,x,y,z,t) + -- 2 Uvv(v,w,x,y,z,t) + 2 Uww(v,w,x,y,z,t) + 2 Uxx(v,w,x,y,z,t) + -- 2 Uyy(v,w,x,y,z,t) + 2 Uzz(v,w,x,y,z,t) + 2 Utt(v,w,x,y,z,t) + -- 6 U(v,w,x,y,z,t) = 0 -- -- boundary conditions computed using U(v,w,x,y,z,t) -- analytic solution is U(v,w,x,y,z,t) = sin(v+w+x+y+z+t) 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 nuderiv; with ada.calendar; use ada.Calendar; procedure pde46h_nu is nv : integer := 6; -- problem parameters nw : integer := 6; nx : integer := 6; ny : integer := 6; nz : integer := 6; nt : integer := 6; nvwxyzt : integer := (nv-2)*(nw-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2); A: real_matrix(1..nvwxyzt,1..nvwxyzt); -- A * U = R R : real_vector(1..nvwxyzt); -- RHS known U : real_vector(1..nvwxyzt); -- solution being computed Ua : real_vector(1..nvwxyzt); -- solution for checking vg : real_vector(1..nv); wg : real_vector(1..nw); xg : real_vector(1..nx); yg : real_vector(1..ny); zg : real_vector(1..nz); tg : real_vector(1..nt); vmin : real := 0.0; -- problem parameters vmax : real := 1.5708; wmin : real := 0.0; wmax : real := 1.5708; xmin : real := 0.0; xmax : real := 1.5708; ymin : real := 0.0; ymax : real := 1.5708; zmin : real := 0.0; zmax : real := 1.5708; tmin : real := 0.0; tmax : real := 1.5708; hv, hw, hx, hy, hz, ht : real; package real_io is new float_io(real); use real_io; -- internal functions function s(i:integer; ii:integer; iii:integer; iiii:Integer; iiiii:Integer; iiiiii:Integer) return integer is -- for internal x,y,z,t begin return (i-2)*(nw-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2) + (ii-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2) + (iii-2)*(ny-2)*(nz-2)*(nt-2) + (iiii-2)*(nz-2)*(nt-2) + (iiiii-2)*(nt-2) + (iiiiii-1); end S; -- problem RHS, 0 -- problem boundaries (and analytic solution for checking) function uana(v:real; w:real; x:real; y:real; z:real; t:real) return real is begin return sin(v+w+x+y+z+t); end uana; function eval_derivative(vord:integer; word:integer; xord:integer; yord:integer; zord:integer; tord:integer; i:integer; ii:integer; iii:integer; iiii:integer; iiiii:integer; iiiiii:integer; u:real_vector) return real is cv : real_vector(1..nv); -- can be used for any order cw : real_vector(1..nw); cx : real_vector(1..nx); cy : real_vector(1..ny); cz : real_vector(1..nz); ct : real_vector(1..nt); discrete : real; v, w, x, y, z, t : real; begin v := vg(i); w := wg(ii); x := xg(iii); y := yg(iiii); z := zg(iiiii); t := tg(iiiiii); discrete := 0.0; -- six cases of one partial v, w, x, y, z, t to any order if vord/=0 and word=0 and xord=0 and yord=0 and zord=0 and tord=0 then nuderiv(vord, nv, i, vg, cv); for j in 1..nv loop if j=1 or j=nv then discrete := discrete + cv(j)*uana(vg(j),w,x,y,z,t); else discrete := discrete + cv(j)*u(s(j,ii,iii,iiii,iiiii,iiiiii)); end if; end loop; -- j elsif vord=0 and word/=0 and xord=0 and yord=0 and zord=0 and tord=0 then nuderiv(word, nw, ii, wg, cw); for jj in 1..nw loop if jj=1 or jj=nw then discrete := discrete + cw(jj)*uana(v,wg(jj),x,y,z,t); else discrete := discrete + cw(jj)*u(s(i,jj,iii,iiii,iiiii,iiiiii)); end if; end loop; -- jj elsif vord=0 and word=0 and xord/=0 and yord=0 and zord=0 and tord=0 then nuderiv(xord, nx, iii, xg, cx); for jjj in 1..nx loop if jjj=1 or jjj=nx then discrete := discrete + cx(jjj)*uana(v,w,xg(jjj),y,z,t); else discrete := discrete + cx(jjj)*u(s(i,ii,jjj,iiii,iiiii,iiiiii)); end if; end loop; -- jjj elsif vord=0 and word=0 and xord=0 and yord/=0 and zord=0 and tord=0 then nuderiv(yord, ny, iiii, yg, cy); for jjjj in 1..ny loop if jjjj=1 or jjjj=ny then discrete := discrete + cy(jjjj)*uana(v,w,x,yg(jjjj),z,t); else discrete := discrete + cy(jjjj)*u(s(i,ii,iii,jjjj,iiiii,iiiiii)); end if; end loop; -- jjjj elsif vord=0 and word=0 and xord=0 and yord=0 and zord/=0 and tord=0 then nuderiv(zord, nz, iiiii, zg, cz); for jjjjj in 1..nz loop if jjjjj=1 or jjjjj=nz then discrete := discrete + cz(jjjjj)*uana(v,w,x,y,zg(jjjjj),t); else discrete := discrete + cz(jjjjj)*u(s(i,ii,iii,iiii,jjjjj,iiiiii)); end if; end loop; -- jjjjj elsif vord=0 and word=0 and xord=0 and yord=0 and zord=0 and tord/=0 then nuderiv(tord, nt, iiiiii, tg, ct); for jjjjjj in 1..nt loop if jjjjjj=1 or jjjjjj=nt then discrete := discrete + ct(jjjjjj)*uana(v,w,x,y,z,tg(jjjjjj)); else discrete := discrete + ct(jjjjjj)*u(s(i,ii,iii,iiii,iiiii,jjjjjj)); end if; end loop; -- jjjjjj end if; -- other cases not needed for this PDE return discrete; end eval_derivative; procedure check_soln(U:real_vector) is val, err, smaxerr : real; bi, bii, biii, biiii, biiiii, biiiiii : integer; begin smaxerr := 0.0; for i in 2..nv-1 loop for ii in 2..nw-1 loop for iii in 2..nx-1 loop for iiii in 2..ny-1 loop for iiiii in 2..nz-1 loop for iiiiii in 2..nt-1 loop val := 0.0; -- + Uvvvv(v,w,x,y,z,t) val := val + eval_derivative(4,0,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u); -- + Uwwww(v,w,x,y,z,t) val := val + eval_derivative(0,4,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u); -- + Uxxxx(v,w,x,y,z,t) val := val + eval_derivative(0,0,4,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u); -- + Uyyyy(v,w,x,y,z,t) val := val + eval_derivative(0,0,0,4,0,0,i,ii,iii,iiii,iiiii,iiiiii,u); -- + Uzzzz(v,w,x,y,z,t) val := val + eval_derivative(0,0,0,0,4,0,i,ii,iii,iiii,iiiii,iiiiii,u); -- + Utttt(v,w,x,y,z,t) val := val + eval_derivative(0,0,0,0,0,4,i,ii,iii,iiii,iiiii,iiiiii,u); -- + 2 Uvv(v,w,x,y,z,t) val := val + 2.0*eval_derivative(2,0,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u); -- + 2 Uww(v,w,x,y,z,t) val := val + 2.0*eval_derivative(0,2,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u); -- + 2 Uxx(v,w,x,y,z,t) val := val + 2.0*eval_derivative(0,0,2,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u); -- + 2 Uyy(v,w,x,y,z,t) val := val + 2.0*eval_derivative(0,0,0,2,0,0,i,ii,iii,iiii,iiiii,iiiiii,u); -- + 2 Uzz(v,w,x,y,z,t) val := val + 2.0*eval_derivative(0,0,0,0,2,0,i,ii,iii,iiii,iiiii,iiiiii,u); -- + 2 Utt(v,w,x,y,z,t) val := val + 2.0*eval_derivative(0,0,0,0,0,2,i,ii,iii,iiii,iiiii,iiiiii,u); -- + 2 u(v,w,x,y,z,t) val := val + 6.0*uana(vg(i),wg(ii),xg(iii),yg(iiii),zg(iiiii),tg(iiiiii)); -- should be zero err := abs(val); if err>smaxerr then smaxerr := err; bi := i; bii := ii; biii := iii; biiii := iiii; biiiii := iiiii; biiiiii := iiiiii; end if; end loop; --iiiiii end loop; --iiiii end loop; --iiii end loop; --iii 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(biiii)&integer'image(biiii)& integer'image(biiiii)&integer'image(biiiiii)); end Check_soln; procedure write_soln(Which:Integer; tag:String) is fout : File_Type; begin create(fout, out_file, "pde46h_nu_ada"&tag&".dat"); put_line("writing pde46h_nu_ada"&tag&".dat"); put_line(fout,"6 "&integer'image(nv)&integer'image(nw)& integer'image(nx)&integer'image(ny)& integer'image(nz)&integer'image(nt)); for i in 1..nv loop for ii in 1..nw loop for iii in 1..nx loop for iiii in 1..ny loop for iiiii in 1..nz loop for iiiiii in 1..nt loop put(fout," "); put(fout,vg(i),2,6,0); put(fout," "); put(fout,wg(ii),2,6,0); put(fout," "); put(fout,xg(iii),2,6,0); put(fout," "); put(fout,yg(iiii),2,6,0); put(fout," "); put(fout,zg(iiiii),2,6,0); put(fout," "); put(fout,tg(iiiiii),2,6,0); put(fout," "); if Which=1 then -- expected solution put(fout,uana(vg(i),wg(ii),xg(iii),yg(iiii), zg(iiiii),tg(iiiiii)),4,6,0); else -- computed solution if i=1 or i=nv or ii=1 or ii=nw or iii=1 or iii=nx or iiii=1 or iiii=ny or iiiii=1 or iiiii=nz or iiiiii=1 or iiiiii=nt then -- boundary put(fout,uana(vg(i),wg(ii),xg(iii),yg(iiii), zg(iiiii),tg(iiiiii)),4,6,0); else put(fout,U(S(i,ii,iii,iiii,iiiii,iiiiii)),4,6,0); end if; end if; new_line(fout); end loop; -- iiiiii end loop; -- iiiii end loop; -- iiii end loop; -- iii end loop; -- ii end loop; -- i close(fout); put_line("finished writing pde46h_nu_ada"&tag&".dat"); end write_soln; procedure initializeA is begin put_line("initialize A "); -- initialize A for i in 2..nv-1 loop for ii in 2..nw-1 loop for iii in 2..nx-1 loop for iiii in 2..ny-1 loop for iiiii in 2..nz-1 loop for iiiiii in 2..nt-1 loop for j in 2..nv-1 loop for jj in 2..nw-1 loop for jjj in 2..nx-1 loop for jjjj in 2..ny-1 loop for jjjjj in 2..nz-1 loop for jjjjjj in 2..nt-1 loop A(s(i,ii,iii,iiii,iiiii,iiiiii),s(j,jj,jjj,jjjj,jjjjj,jjjjjj)) := 0.0; end loop; -- jjjjjj end loop; -- jjjjj end loop; -- jjjj end loop; -- jjj end loop; -- jj end loop; -- j R(s(i,ii,iii,iiii,iiiii,iiiiii)) := 0.0; end loop; -- iiiiii end loop; -- iiiii end loop; -- iiii end loop; -- iii end loop; -- ii end loop; -- i end initializeA; procedure buildA is -- system of equations val : real := 0.0; k: integer; cvv, cvvvv : real_vector(1..nv); cww, cwwww : real_vector(1..nw); cxx, cxxxx : real_vector(1..nx); cyy, cyyyy : real_vector(1..ny); czz, czzzz : real_vector(1..nz); ctt, ctttt : real_vector(1..nt); begin -- compute entries in A matrix, inside boundary put_line("compute A matrix and Y RHS "); put_line("check nuderiv(4,nv,1,vg,cvvvv"); nuderiv(4,nv,1,vg,cvvvv); for i in 1..nv loop put_line("cvvvv("&integer'image(i)& ")="&real'image(cvvvv(i))); end loop; put_line("check nuderiv(4,nv,nv,vg,cvvvv"); nuderiv(4,nv,nv,vg,cvvvv); for i in 1..nv loop put_line("cvvvv("&integer'image(i)& ")="&real'image(cvvvv(i))); end loop; put_line(" "); val := 0.0; for i in 2..nv-1 loop for ii in 2..nw-1 loop for iii in 2..nx-1 loop for iiii in 2..ny-1 loop for iiiii in 2..nz-1 loop for iiiiii in 2..nt-1 loop k := s(i,ii,iii,iiii,iiiii,iiiiii); -- row index -- for each term in first equation -- for d^4u/dv^4 nuderiv(4,nv,i,vg,cvvvv); for j in 1..nv loop val := cvvvv(j); if j=1 or j=nv then R(k) := R(k) - val*uana(vg(j),wg(ii),xg(iii),yg(iiii), zg(iiiii),tg(iiiiii)); else A(k,s(j,ii,iii,iiii,iiiii,iiiiii)) := A(k,s(j,ii,iii,iiii,iiiii,iiiiii)) + val; end if; end loop; -- j -- for d^4u/dw^4 nuderiv(4,nw,ii,wg,cwwww); for jj in 1..nw loop val := cwwww(jj); if jj=1 or jj=nw then R(k) := R(k) - val*uana(vg(i),wg(jj),xg(iii),yg(iiii), zg(iiiii),tg(iiiiii)); else A(k,s(i,jj,iii,iiii,iiiii,iiiiii)) := A(k,s(i,jj,iii,iiii,iiiii,iiiiii)) + val; end if; end loop; -- jj -- for d^4u/dx^4 nuderiv(4,nx,iii,xg,cxxxx); for jjj in 1..nx loop val := cxxxx(jjj); if jjj=1 or jjj=nx then R(k) := R(k) - val*uana(vg(i),wg(ii),xg(jjj),yg(iiii), zg(iiiii),tg(iiiiii)); else A(k,s(i,ii,jjj,iiii,iiiii,iiiiii)) := A(k,s(i,ii,jjj,iiii,iiiii,iiiiii)) + val; end if; end loop; -- jjj -- for d^4u/dy^4 nuderiv(4,ny,iiii,yg,cyyyy); for jjjj in 1..ny loop val := cyyyy(jjjj); if jjjj=1 or jjjj=ny then R(k) := R(k) - val*uana(vg(i),wg(ii),xg(iii),yg(jjjj), zg(iiiii),tg(iiiiii)); else A(k,s(i,ii,iii,jjjj,iiiii,iiiiii)) := A(k,s(i,ii,iii,jjjj,iiiii,iiiiii)) + val; end if; end loop; -- jjjj -- for d^4u/dz^4 nuderiv(4,nz,iiiii,zg,czzzz); for jjjjj in 1..nz loop val := czzzz(jjjjj); if jjjjj=1 or jjjjj=nz then R(k) := R(k) - val*uana(vg(i),wg(ii),xg(iii),yg(iiii), zg(jjjjj),tg(iiiiii)); else A(k,s(i,ii,iii,iiii,jjjjj,iiiiii)) := A(k,s(i,ii,iii,iiii,jjjjj,iiiiii)) + val; end if; end loop; -- jjjjj -- for d^4u/dt^4 nuderiv(4,nt,iiiiii,tg,ctttt); for jjjjjj in 1..nt loop val := ctttt(jjjjjj); if jjjjjj=1 or jjjjjj=nt then R(k) := R(k) - val*uana(vg(i),wg(ii),xg(iii),yg(iiii), zg(iiiii),tg(jjjjjj)); else A(k,s(i,ii,iii,iiii,iiiii,jjjjjj)) := A(k,s(i,ii,iii,iiii,iiiii,jjjjjj)) + val; end if; end loop; -- jjjjjj -- for d^2u/dv^2 nuderiv(2,nv,i,vg,cvv); for j in 1..nv loop val := 2.0*cvv(j); if j=1 or j=nv then R(k) := R(k) - val*uana(vg(j),wg(ii),xg(iii),yg(iiii), zg(iiiii),tg(iiiiii)); else A(k,s(j,ii,iii,iiii,iiiii,iiiiii)) := A(k,s(j,ii,iii,iiii,iiiii,iiiiii)) + val; end if; end loop; -- j -- for d^2u/dw^2 nuderiv(2,nw,ii,wg,cww); for jj in 1..nw loop val := 2.0*cww(jj); if jj=1 or jj=nw then R(k) := R(k) - val*uana(vg(i),wg(jj),xg(iii),yg(iiii), zg(iiiii),tg(iiiiii)); else A(k,s(i,jj,iii,iiii,iiiii,iiiiii)) := A(k,s(i,jj,iii,iiii,iiiii,iiiiii)) + val; end if; end loop; -- jj -- for d^2u/dx^2 nuderiv(2,nx,iii,xg,cxx); for jjj in 1..nx loop val := 2.0*cxx(jjj); if jjj=1 or jjj=nx then R(k) := R(k) - val*uana(vg(i),wg(ii),xg(jjj),yg(iiii), zg(iiiii),tg(iiiiii)); else A(k,s(i,ii,jjj,iiii,iiiii,iiiiii)) := A(k,s(i,ii,jjj,iiii,iiiii,iiiiii)) + val; end if; end loop; -- jjj -- for d^2u/dy^2 nuderiv(2,ny,iiii,yg,cyy); for jjjj in 1..ny loop val := 2.0*cyy(jjjj); if jjjj=1 or jjjj=ny then R(k) := R(k) - val*uana(vg(i),wg(ii),xg(iii),yg(jjjj), zg(iiiii),tg(iiiiii)); else A(k,s(i,ii,iii,jjjj,iiiii,iiiiii)) := A(k,s(i,ii,iii,jjjj,iiiii,iiiiii)) + val; end if; end loop; -- jjjj -- for d^2u/dz^2 nuderiv(2,nz,iiiii,zg,czz); for jjjjj in 1..nz loop val := 2.0*czz(jjjjj); if jjjjj=1 or jjjjj=nz then R(k) := R(k) - val*uana(vg(i),wg(ii),xg(iii),yg(iiii), zg(jjjjj),tg(iiiiii)); else A(k,s(i,ii,iii,iiii,jjjjj,iiiiii)) := A(k,s(i,ii,iii,iiii,jjjjj,iiiiii)) + val; end if; end loop; -- jjjjj -- for d^2u/dt^2 nuderiv(2,nt,iiiiii,tg,ctt); for jjjjjj in 1..nt loop val := 2.0*ctt(jjjjjj); if jjjjjj=1 or jjjjjj=nt then R(k) := R(k) - val*uana(vg(i),wg(ii),xg(iii),yg(iiii), zg(iiiii),tg(jjjjjj)); else A(k,s(i,ii,iii,iiii,iiiii,jjjjjj)) := A(k,s(i,ii,iii,iiii,iiiii,jjjjjj)) + val; end if; end loop; -- jjjjjj -- for u(v,w,x,y,z,t) A(k,s(i,ii,iii,iiii,iiiii,iiiiii)) := A(k,s(i,ii,iii,iiii,iiiii,iiiiii)) + 6.0; -- for f(v,w,x,y,z,t) RHS constant, zero here -- R(k) := R(k) + f(vg(i),wg(ii),xg(iii),yg(iiii), -- zg(iiiii),tg(iiiiii)); -- end terms end loop; -- iiiiii end loop; -- iiiii end loop; -- iiii end loop; -- iii end loop; -- ii end loop; -- i end BuildA; procedure main is v, w, x, y, z, t : real; err, avgerr, maxerr : real; time_start : duration; time_now : duration; time_last : duration; begin -- main put_line("pde46h_nu.adb running "); put_line("Given Uvvvv(v,w,x,y,z,t) + Uwwww(v,w,x,y,z,t) + Uxxxx(v,w,x,y,z,t) + "); put_line(" Uyyyy(v,w,x,y,z,t) + Uzzzz(v,w,x,y,z,t) + Utttt(v,w,x,y,z,t) + "); put_line(" 2 Uvv(v,w,x,y,z,t) + 2 Uww(v,w,x,y,z,t) + 2 Uxx(v,w,x,y,z,t) + "); put_line(" 2 Uyy(v,w,x,y,z,t) + 2 Uzz(v,w,x,y,z,t) + 2 Utt(v,w,x,y,z,t) + "); put_line(" 6 u(v,w,x,y,z,t) = 0 "); put_line("vmin<=v<=vmax wmin<=w<=wmax xmin<=x<=xmax Boundaries "); put_line("ymin<=y<=ymax zmin<=z<=zmax tmin<=t<=tmax Boundaries "); put_line("Analytic solution u(v,w,x,y,z,t) = sin(v+w+x+y+z+t)"); put_line(" "); put_line("vmin="&real'image(vmin)&", vmax="&real'image(vmax)); put_line("wmin="&real'image(wmin)&", wmax="&real'image(wmax)); 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("nv="&integer'image(nv)&", nw="&integer'image(nw)& ", nx="&integer'image(nx)&", ny="&integer'image(ny)& ", nz="&integer'image(nz)&", nt="&integer'image(nt)); put_line(" "); time_start := seconds(clock); time_last := time_start; hv := real(nv*(nv-1))/2.0; put_line("v divisor, hv="&real'image(hv)); vg(1) := vmin; put_line("vg(1))="&real'image(vg(1))); for i in 2..nv loop -- decreasing step size vg(i) := vg(i-1)+(vmax-vmin)*real(nv+1-i)/hv; put_line("vg("&integer'image(i)&")="&real'image(vg(i))); end loop; -- i put_line(" "); hw := real(nw*(nw-1))/2.0; put_line("w divisor, hw="&real'image(hw)); wg(1) := wmin; put_line("wg(1))="&real'image(wg(1))); for ii in 2..nw loop -- decreasing step size wg(ii) := wg(ii-1)+(wmax-wmin)*real(nw+1-ii)/hw; put_line("wg("&integer'image(ii)&")="&real'image(wg(ii))); end loop; -- ii put_line(" "); hx := real(nx*(nx-1))/2.0; put_line("x divisor, hx="&real'image(hx)); xg(1) := xmin; put_line("xg(1))="&real'image(xg(1))); for iii in 2..nx loop -- decreasing step size xg(iii) := xg(iii-1)+(xmax-xmin)*real(nx+1-iii)/hx; put_line("xg("&integer'image(iii)&")="&real'image(xg(iii))); end loop; -- iii put_line(" "); hy := real(ny*(ny-1))/2.0; put_line("y divisor, hy="&real'image(hy)); yg(1) := ymin; put_line("yg(1))="&real'image(yg(1))); for iiii in 2..ny loop -- decreasing step size yg(iiii) := yg(iiii-1)+(ymax-ymin)*real(ny+1-iiii)/hy; put_line("yg("&integer'image(iiii)&")="&real'image(yg(iiii))); end loop; -- iiii put_line(" "); hz := real(nz*(nz-1))/2.0; put_line("z divisor, hz="&real'image(hz)); put_line("zg(1))="&real'image(zg(1))); zg(1) := zmin; for iiiii in 2..nz loop -- decreasing step size zg(iiiii) := zg(iiiii-1)+(zmax-zmin)*real(nz+1-iiiii)/hz; put_line("zg("&integer'image(iiiii)&")="&real'image(zg(iiiii))); end loop; -- iiiii put_line(" "); ht := real(nt*(nt-1))/2.0; put_line("t divisor, hz="&real'image(ht)); tg(1) := tmin; put_line("tg(1))="&real'image(tg(1))); for iiiiii in 2..nt loop -- decreasing step size tg(iiiiii) := tg(iiiiii-1)+(tmax-tmin)*real(nt+1-iiiiii)/ht; put_line("tg("&integer'image(iiiiii)&")="&real'image(tg(iiiiii))); end loop; -- iiiiii put_line(" "); put_line("put solution in Ua vector "); -- put solution in Ua vector for i in 2..nv-1 loop v := vg(i); for ii in 2..nw-1 loop w := wg(ii); for iii in 2..nx-1 loop x := xg(iii); for iiii in 2..ny-1 loop y := yg(iiii); for iiiii in 2..nz-1 loop z := zg(iiiii); for iiiiii in 2..nt-1 loop t := tg(iiiiii); Ua(s(i,ii,iii,iiii,iiiii,iiiiii)) := uana(v,w,x,y,z,t); end loop; -- iiiiii end loop; -- iiiii end loop; -- iiii end loop; -- iii 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(nvwxyzt)&" equations in "& integer'image(nvwxyzt)&" unknowns "); U := 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(U); put_line(" "); avgerr := 0.0; maxerr := 0.0; put_line(" U computed, Ua analytic, error "); for i in 2..nv-1 loop for ii in 2..nw-1 loop for iii in 2..nx-1 loop for iiii in 2..ny-1 loop for iiiii in 2..nz-1 loop for iiiiii in 2..nt-1 loop err := abs(U(s(i,ii,iii,iiii,iiiii,iiiiii))- Ua(s(i,ii,iii,iiii,iiiii,iiiiii))); 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)& ","&integer'image(iiiii)&","&integer'image(iiiiii)&")="); put(U(s(i,ii,iii,iiii,iiiii,iiiiii)),3,6); put(", Ua="); put(Ua(s(i,ii,iii,iiii,iiiii,iiiiii)),3,6); put(", err="); put(err); new_line; end loop; -- iiiiii end loop; -- iiiii end loop; -- iiii end loop; -- iii 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(nx*ny*nz*nt))); put_line(" "); end Main; begin -- pde46h_nu main; -- write_soln(1,"u"); -- analytic solution write_soln(2,""); -- computed solution for plot6d_gl put_line("end pde46h_nu.adb"); end pde46h_nu;