-- pde_abc_eq.adb see WEB page for development CMSC 455 supplemental -- The PDE to be solved is: -- a1(x,y)*uxx(x,y) + b1(x,y)*uxy(x,y) + c1(x,y)*uyy(x,y) + -- d1(x,y)*ux(x,y) + e1(x,y)*uy(x,y) + f1(x,y)*u(x,y) = c(x,y) -- -- The functions and other user information is in the files: -- abc.ads abc.adb -- also compile with simeq.adb deriv.adb -- real_arrays.ads real_arrays.adb with Ada.Text_IO; use Ada.Text_IO; with Real_Arrays; use Real_Arrays; with Ada.Numerics.Long_Elementary_Functions; use Ada.Numerics.Long_Elementary_Functions; with abc; use abc; with simeq; with deriv; with rderiv; procedure pde_abc_eq is -- nx and ny are defined in abc.ads us:Real_Vector(1..(nx-2)*(ny-2)); -- solution being computed, X uu:Real_Vector(1..(nx-2)*(ny-2)); -- solution if checking ut:Real_Matrix(1..(nx-2)*(ny-2),1..(nx-2)*(ny-2)); -- solution matrix, A uc:Real_Vector(1..(nx-2)*(ny-2)); -- equation constant, Y -- nx-2 by ny-2 internal, non boundary cells -- all elements 1..nx by 1..ny include boundary elements -- xmin, xmax, ymin, ymax are defined in abc.ads hx:Real := (xmax-xmin)/Real(nx-1); -- x step hy:Real := (ymax-ymin)/Real(ny-1); -- y step hx2:Real := hx*hx; hy2:Real := hy*hy; hxy:Real := hx*hy; -- cs; -- index of constant vector now in uc(ij) pdxx:Real_Matrix(1..nd,1..nd); pdxy:Real_Matrix(1..nd,1..nd); pdyy:Real_Matrix(1..nd,1..nd); pdx :Real_Matrix(1..nd,1..nd); pdy :Real_Matrix(1..nd,1..nd); pd :Real_Matrix(1..nd,1..nd); md:integer := (nd+1)/2 ; -- mid point of 1..nd, nd is odd error:Real; -- sum of absolute exact-computed avgerror:Real; -- average error maxerror:Real; -- maximum cell error package Real_IO is new Float_IO(Real); use Real_IO; function S(i:Integer; j:Integer) return Integer is begin return i+(nx-2)*(j-1); end S; function max(i:Integer; j:Integer) return Integer is begin if i>j then return i; end if; return j; end max; function max(x:real; y:real) return real is begin if x>y then return x; end if; return y; end max; procedure clear is -- just zero arrays ii : Integer; begin for i in 1..nx-2 loop -- solution, ut, coordinates -- solution coordinates OK for xmin+i*hx, ymin+j*hy for j in 1..ny-2 loop ii := S(I,J); us(ii) := 0.0; uc(ii) := 0.0; if ifcheck>0 then -- all element coordinates uu(ii) := u(xmin+Real(i)*hx,ymin+Real(j)*hy); end if; end loop; end loop; end clear; procedure printus is -- print the computed solution ii:Integer; begin Put_Line("computed solution"); Put_Line("us(x,y) "); for i in 1..nx-2 loop for j in 1..ny-2 loop ii := S(i,J); Put("i="&Integer'Image(i)& ", j="&Integer'Image(j)& ", at x="); Put(xmin+Real(i)*hx, 5,3,0); Put(", y="); Put(ymin+Real(j)*hy, 5,3,0); Put(", us(i,j)="); Put(us(ii), 6,3,0); New_Line; end loop; end loop; New_Line; end printus; procedure check is -- typically not known, instrumentation ii:Integer; uexact, err:Real; begin error := 0.0; maxerror := 0.0; for i in 1..nx-2 loop -- solution coordinates for j in 1..ny-2 loop ii := i+(nx-2)*(j-1); uexact := u(xmin+Real(i)*hx, ymin+Real(j)*hy); err := abs(us(ii)-uexact); error := error + err; if err>maxerror then maxerror:=err; end if; end loop; end loop; end check; procedure printexact is -- typically not known begin Put_Line("exact solution"); Put_Line("u(x,y) "); for i in 1..nx-2 loop for j in 1..ny-2 loop Put("i="&Integer'Image(i)& ", j="&Integer'Image(j)& ", at x="); Put(xmin+Real(i)*hx, 5,3,0); Put(", y="); Put(ymin+Real(j)*hy, 5,3,0); Put(", u(i,j)="); Put(u(xmin+Real(i)*hx,ymin+Real(j)*hy), 6,3,0); New_Line; end loop; end loop; New_Line; end printexact; procedure printdiff is -- only when checking ii:Integer; begin Put_Line("difference exact-computed solution"); Put_Line("u(x,y)-us(x,y) "); for i in 1..nx-2 loop for j in 1..ny-2 loop ii := i+(nx-2)*(j-1); Put("i="&Integer'Image(i)& ", j="&Integer'Image(j)& ", at x="); Put(xmin+Real(i)*hx, 5,3,0); Put(", y="); Put(ymin+Real(j)*hy, 5,3,0); Put(", diff(i,j)="); Put(u(xmin+Real(i)*hx,ymin+Real(j)*hy)-us(ii), 6,3,0); New_Line; end loop; end loop; New_Line; end printdiff; procedure print_coef is -- for debugging begin Put_Line("pdxx j=1 j=2 j=3 j=4 j=5 y"); for i in 1..nd loop Put("i="&Integer'Image(i)); for j in 1..nd loop Put(pdxx(i,j), 6,2,0); end loop; New_Line; end loop; Put_Line("x"); New_Line; Put_Line("pdxy j=1 j=2 j=3 j=4 j=5 y"); for i in 1..nd loop Put("i="&Integer'Image(i)); for j in 1..nd loop Put(pdxy(i,j), 6,2,0); end loop; New_Line; end loop; New_Line; Put_Line("pdyy j=1 j=2 j=3 j=4 j=5 y"); for i in 1..nd loop Put("i="&Integer'Image(i)); for j in 1..nd loop Put(pdyy(i,j), 6,2,0); end loop; New_Line; end loop; New_Line; Put_Line("pdx j=1 j=2 j=3 j=4 j=5 y"); for i in 1..nd loop Put("i="&Integer'Image(i)); for j in 1..nd loop Put(pdx(i,j), 6,2,0); end loop; New_Line; end loop; New_Line; Put_Line("pdy j=1 j=2 j=3 j=4 j=5 y"); for i in 1..nd loop Put("i="&Integer'Image(i)); for j in 1..nd loop Put(pdy(i,j), 6,2,0); end loop; New_Line; end loop; New_Line; end print_coef; procedure build_coef(ix:Integer; iy:Integer) is -- ix and iy in 1..nd are 'point' to "deriv" and index a:Real_Vector(1..nd*nd); begin -- build discretization matrices for i in 1..nd loop for j in 1..nd loop pdxx(i,j):=0.0; -- coefficient of a1(x,y) pdxy(i,j):=0.0; -- coefficient of b1(x,y) pdyy(i,j):=0.0; -- coefficient of c1(x,y) pdx (i,j):=0.0; -- coefficient of d1(x,y) pdy (i,j):=0.0; -- coefficient of e1(x,y) pd (i,j):=0.0; -- coefficient of f1(x,y) end loop; end loop; pd(ix,iy) := 1.0; -- u(i,j) term deriv(2, nd, ix, a); -- for pdxx for i in 1..nd loop pdxx(i,iy) := a(i)/hx2; end loop; deriv(2, nd, iy, a); -- for pdyy for j in 1..nd loop pdyy(ix,j) := a(j)/hy2; end loop; deriv(1, nd, ix, a); -- for pdx for i in 1..nd loop pdx(i,iy) := a(i)/hx; end loop; deriv(1, nd, iy, a); -- for pdy for j in 1..nd loop pdy(ix,j) := a(j)/hy; end loop; for i in 1..nd loop for j in 1..nd loop pdxy(i,j) := pdx(i,iy)*pdy(ix,j); end loop; end loop; if ifcheck>2 then print_coef; end if; end build_coef; procedure check_row(row : Integer) is -- only for checking -- row is in ut solution coordinates ii : Integer; sum : Real; begin sum := 0.0; for i in 1..nx-2 loop for j in 1..ny-2 loop ii := i+(nx-2)*(j-1); sum := sum + ut(row,ii)*uu(ii); -- uu is check solution end loop; end loop; sum := sum - uc(row); if abs(sum)>0.001 then Put_Line("row="&Integer'Image(row)&" is bad err="& Real'Image(sum)); end if; end check_row; procedure build_row(i:Integer; j:Integer; bi:Integer; bj:integer) is -- i and j are solution coordinates, ir and jr are all element ii, ij, ir, jr, ia, ja: integer; p, a1c, b1c, c1c, d1c, e1c, f1c, cc : real; begin ir := i + 1; -- all element coordinates jr := j + 1; ii := i+(nx-2)*(j-1); -- row of solution, ut, matrix if ifcheck>2 then Put_Line("working row="&Integer'Image(ii)& ", i="&Integer'Image(i)& ", j="&Integer'Image(j)& ", bi="&Integer'Image(bi)& ", bj="&Integer'Image(bj)& ", ir="&Integer'Image(ir)& ", jr="&Integer'Image(jr)); end if; -- check a1c := a1(xmin+Real(i)*hx,ymin+Real(j)*hy); -- yes, ir-1, jr-1 b1c := b1(xmin+Real(i)*hx,ymin+Real(j)*hy); c1c := c1(xmin+Real(i)*hx,ymin+Real(j)*hy); d1c := d1(xmin+Real(i)*hx,ymin+Real(j)*hy); e1c := e1(xmin+Real(i)*hx,ymin+Real(j)*hy); f1c := f1(xmin+Real(i)*hx,ymin+Real(j)*hy); cc := c(xmin+Real(i)*hx,ymin+Real(j)*hy); if ifcheck>1 then Put_Line("a1="&Real'Image(a1c)); Put_Line("b1="&Real'Image(b1c)); Put_Line("c1="&Real'Image(c1c)); Put_Line("d1="&Real'Image(d1c)); Put_Line("e1="&Real'Image(e1c)); Put_Line("f1="&Real'Image(f1c)); Put_Line("c ="&Real'Image(cc)); end if; -- check for ki in 1..nd loop -- discretization for kj in 1..nd loop ia := ir-bi+ki-md; -- all element coordinates 1 and nx are boundary ja := jr-bj+kj-md; -- may be outside solution, thus boundary ij := (ia-1)+(nx-2)*(ja-2); -- solution, ut, coordinates p := a1c*pdxx(ki,kj) + b1c*pdxy(ki,kj) + c1c*pdyy(ki,kj)+ d1c* pdx(ki,kj) + e1c* pdy(ki,kj) + f1c* pd(ki,kj); if ifcheck>1 then Put_Line("ia="&Integer'Image(ia)&", ja="&Integer'Image(ja)& ", ij="&Integer'Image(ij)&", p="&Real'Image(p)); end if; if ia=1 then -- boundary value tests (correct for starting with 1) uc(ii) := uc(ii) - p * u(xmin+real(ia-1)*hx, ymin+real(ja-1)*hy); elsif ia=nx then uc(ii) := uc(ii) - p * u(xmin+real(ia-1)*hx, ymin+real(ja-1)*hy); elsif ja=1 then uc(ii) := uc(ii) - p * u(xmin+real(ia-1)*hx, ymin+real(ja-1)*hy); elsif ja=ny then uc(ii) := uc(ii) - p * u(xmin+real(ia-1)*hx, ymin+real(ja-1)*hy); else -- non boundary case ut(ii,ij) := ut(ii,ij) + p; Put_Line("ut(ii,"&Integer'Image(ij)&")="&Real'Image(ut(ii,ij))); end if; end loop; end loop; uc(ii) := uc(ii) + cc; Put_Line("row ii="&Integer'Image(ii)&", uc(ii)="&Real'Image(uc(ii))); check_row(ii); end build_row; procedure init_matrix is -- build the ut matrix i, j : integer; begin -- zero internal cells, the matrix does not include boundary cells for i in 1..nx-2 loop for j in 1..ny-2 loop for ii in 1..nx-2 loop for jj in 1..ny-2 loop ut(i+(nx-2)*(j-1),ii+(nx-2)*(jj-1)) := 0.0; end loop; end loop; end loop; end loop; Put_Line("internal cells zeroed, do four corners"); -- do four corners (special case, solution, ut, coordinates) i:=1; j:=1; build_coef(md-1, md-1); build_row(i, j, -1, -1); i:=1; j:=ny-2; build_coef(md-1, md+1); build_row(i, j, -1, 1); i:=nx-2; j:=1; build_coef(md+1, md-1); build_row(i, j, 1, -1); i:=nx-2; j:=ny-2; build_coef(md+1, md+1); build_row(i, j, 1, 1); Put_Line("four corners generated. do central cells"); -- central cells in solution, ut, coordinates build_coef(md, md); for i in 2..nx-3 loop -- compute standard central cells for j in 2..ny-3 loop build_row(i, j, 0, 0); end loop; end loop; Put_Line("standard central cells generated, do sides"); -- do four sides, one element in. special cases build_coef(md-1, md); for j in 2..ny-3 loop -- do i=1 i := 1; build_row(i, j, -1, 0); end loop; -- on edge build_coef(md+1, md); for j in 2..ny-3 loop -- do i=nx-2 i := nx-2; build_row(i, j, 1, 0); end loop; -- on edge build_coef(md, md-1); for i in 2..nx-3 loop -- do j=1 j := 1; build_row(i, j, 0 , -1); end loop; -- on edge build_coef(md, Md+1); for i in 2..nx-3 loop -- do j=ny-2 j := ny-2; build_row(i, j, 0 , 1); end loop; -- on edge Put_Line("four sides generated"); end init_matrix; procedure print_mat is sum : Real; begin Put_Line("initial matrix"); for i in 1..(nx-2)*(ny-2) loop for J in 1..(nx-2)*(ny-2) loop Put(ut(i,j),5,2,0); if (nx-2)*(j/(nx-2))=J and j/=(nx-2)*(ny-2) then New_Line; end if; end loop; Put(uc(i),5,2,0); sum := 0.0; for ii in 1..(nx-2)*(ny-2) loop sum := sum + ut(i,ii)*uu(ii); -- uu is check solution end loop; Put("="); Put(sum,5,2,0); New_Line; end loop; New_Line; end print_mat; procedure check_soln is -- against PDE, added later, used rderiv and xg,yg -- for all interior points x, y, rmserr, sumerr, sumsqerr, maxerr, Avgerr, err : real; cx : Real_Matrix(1..nx,1..nx); cxx : Real_Matrix(1..nx,1..nx); cy : Real_Matrix(1..ny,1..ny); cyy : Real_Matrix(1..ny,1..ny); tcxy : Real_Vector(1..max(nx,ny)); dxx, dxy, dyy, dx, dy : Real; -- a1, b1, c1, d1, e1 coefficients ug : Real_Matrix(1..nx,1..ny); -- full, including boundaries, -- copied from us() xg : Real_Vector(1..nx); yg : Real_Vector(1..ny); cnt : Integer := 0; ii : Integer; begin If ifcheck>5 then Put_Line("check_soln setup"); end if; for i in 1..nx loop xg(i) := xmin + real(i-1)*Hx; end loop; for j in 1..ny loop yg(j) := ymin + real(j-1)*Hy; end loop; -- transform us() into ug, add boundaries, make code simpler for i in 2..nx-1 loop -- get from us() for j in 2..ny-1 loop ii := S(i-1,j-1); -- ii for us(ii) computed solution at i,j ug(i,j) := us(ii); if ifcheck>5 then Put_Line("ug("&Integer'Image(i)&","&Integer'Image(j)& ")="&Real'Image(ug(i,j))); end if; end loop; end loop; for i in 1..nx loop -- fill boundaries, loops below easier ug(i,1) := u(xg(i),yg(1)); ug(i,ny) := u(xg(i),yg(ny)); end loop; for j in 1..ny loop -- fill boundaries ug(1,j) := u(xg(1), yg(j)); ug(nx,j) := u(xg(nx),yg(j)); end loop; for i in 1..nx loop rderiv(1,nx,i,hx,tcxy); -- fills row of derivative matrix for k in 1..nx loop cx(i,k) := tcxy(k); end loop; rderiv(2,nx,i,hx,tcxy); for k in 1..nx loop cxx(i,k) := tcxy(k); end loop; end loop; for j in 1..ny loop rderiv(1,ny,j,hy,tcxy); for k in 1..ny loop cy(j,k) := tcxy(k); end loop; rderiv(2,ny,j,hy,tcxy); for k in 1..ny loop cyy(j,k) := tcxy(k); end loop; end loop; if ifcheck>5 then for i in 1..nx loop for k in 1..nx loop Put_Line("cx(i,k)="&Real'Image(cx(i,k))); end loop; New_Line; end loop; New_Line; for i in 1..nx loop for k in 1..nx loop Put_Line("cxx(i,k)="&Real'Image(cxx(i,k))); end loop; New_Line; end loop; New_Line; for j in 1..ny loop for k in 1..ny loop Put_Line("cy(j,k)="&Real'Image(cy(j,k))); end loop; New_Line; end loop; New_Line; for j in 1..ny loop for k in 1..ny loop Put_Line("cyy(j,k)="&Real'Image(cyy(j,k))); end loop; New_Line; end loop; New_Line; end if; -- ifcheck sumerr := 0.0; sumsqerr := 0.0; maxerr := 0.0; for i in 2..nx-1 loop -- check full equation, at every point x := xg(i); for j in 2..ny-1 loop y := yg(j); -- compute derivatives, plug into PDE, subtract RHS, should be zero dxx := 0.0; dx := 0.0; for k in 1..nx loop -- PDE a1(x,y)*uxx(x,y) term dxx := dxx + cxx(i,k)*ug(k,j); -- PDE d1(x,y)*ux(x,y) term dx := dx + cx(i,k)*ug(k,j); end loop; dyy := 0.0; dy := 0.0; for k in 1..ny loop -- PDE c1(x,y)*uyy(x,y) term dyy := dyy + cyy(j,k)*ug(i,k); -- PDE e1(x,y)*uy(x,y) term dy := dy + cy(j,k)*ug(i,k); end loop; -- PDE b1(x,y)*uxy(x,y) term dxy := 0.0; for jj in 1..ny loop -- walk jj in y, computing uxx tcxy(jj) := 0.0; -- numeric ux part of uxy for K in 1..nx loop tcxy(jj) := tcxy(jj) + cx(i,k)*ug(k,jj); end loop; end loop; for k in 1..ny loop -- uy part of uxy dxy := dxy + cy(j,k)*tcxy(k); -- PDE uxy term end loop; -- check if ifcheck>5 then Put("numerical dxx="); Put(dxx,4,3,0); Put(", dxy="); Put(dxy,4,3,0); Put(", dyy="); Put(dyy,4,3,0); New_Line; Put("analytic dxx="); Put(6.0*x+6.0*Y,4,3,0); Put(", dxy="); Put(6.0*x+8.0*y+5.0,4,3,0); Put(", dyy="); Put(12.0*y+8.0*X,4,3,0); New_Line; Put("numerical dx= "); Put(dx,4,3,0); Put(", dy= "); Put(dy,4,3,0); Put(", u= "); Put(ug(i,j),4,3,0); New_Line; Put("analytic dx= "); Put(3.0*x*x+6.0*x*y+4.0*y*y+5.0*y+6.0,4,3,0); Put(", dy= "); Put(6.0*y*y+3.0*x*x+8.0*x*y+5.0*x+7.0,4,3,0); Put(", u= "); Put(u(x,y),4,3,0); New_Line; end if; err := a1(xg(i),yg(j))*dxx + b1(xg(i),yg(j))*dxy + c1(xg(i),yg(j))*dyy + d1(xg(i),yg(j))*dx + e1(xg(i),yg(j))*dy + f1(xg(i),yg(j))*ug(i,j) - c(xg(i),yg(j)); -- RHS err := abs(err); sumerr := sumerr + Err; sumsqerr := sumsqerr + err*Err; maxerr := max(err,maxerr); cnt := cnt + 1; end loop; -- j end loop; -- i rmserr := sqrt(sumsqerr/Real(cnt)); avgerr := sumerr/Real(cnt); Put_Line("check_soln values against PDE "); Put("maxerr="); Put(maxerr,3,6); Put(", rmserr="); Put(rmserr,3,6); Put(", avgerr="); Put(avgerr,3,6); New_Line; end check_soln; begin -- pde_abc_eq Put_Line("pde_abc_eq.adb running"); Put_Line("nx="&Integer'Image(nx)& ", ny="&Integer'Image(ny)); Put_Line("xmin="&Real'Image(xmin)& ", xmax="&Real'Image(xmax)); Put_Line("ymin="&Real'Image(ymin)& ", ymax="&Real'Image(ymax)); Put_Line("hx="&Real'Image(hx)& ", hy="&Real'Image(hy)); Put_Line("hx2="&Real'Image(hx2)& ", hy2="&Real'Image(hy2)); clear; if Ifcheck>0 then Put_Line("md="&Integer'Image(md)); end if; init_matrix; print_mat; us := simeq(ut, uc); printus; printexact; printdiff; if ifcheck>0 then check; avgerror := error/Real((nx-2)*(ny-2)); Put_Line("average error = "&Real'Image(avgerror)); Put_Line("maximum error = "&Real'Image(maxerror)); end if; check_soln; end pde_abc_eq;