-- check_cylinder_deriv.adb cylindrical coordinates -- x = r cos theta 0 <= theta < 2 Pi -- y = r sin theta r < 0 restriction for derivative -- z = z z real -- r = sqrt(x^2+y^2) -- theta = arctan(y,x) if theta<0 theta = 2 Pi + arctan(y,x) -- z = z -- -- compute del U(radius,theta,z) r=radius, t=theta angle, z=height -- del U = (d/dr U, 1/r d/dt U, d/dz U) -- -- compute del^2 U(radius,theta,z) r=radius, t=theta angle, z=height -- del^2 U = d/dr(r d/dr U) + 1/r^2 d^2/dt^2 U + d^2/dz^2 U -- = 2/r d/dr U + d^2/dr^2 U + 1/r^2 d^2/dt^2 U + d^2/dz^2 U with Ada.Text_IO; use Ada.Text_IO; with Ada.Numerics.Long_Elementary_Functions; use Ada.Numerics.Long_Elementary_Functions; with Real_Arrays; use Real_Arrays; with Nuderiv; -- procedure procedure check_cylinder_deriv is package Real_IO is new Ada.Text_IO.Float_IO(Real); use Real_IO; Pi : Real := 3.1415926535897932384626433832795028841971; Nr : Integer := 5; -- number radius values Rmin : Real := 0.2; Rmax : Real := 1.0; Hr : Real := (Rmax-Rmin)/Real(Nr-1); Nt : Integer := 7; -- number angle values Tmin : Real := 0.0; Tmax : Real := 300.0*Pi/180.0; -- Ht not applicable Nz : Integer := 5; -- number Z coordinates Zmin : Real := -1.0; Zmax : Real := 1.0; Hz : Real := (Zmax-Zmin)/Real(Nz-1); Rg : Real_Vector(1..Nr); -- r grid values Tg : Real_Vector(1..Nt); -- theta grid values Zg : Real_Vector(1..Nz); -- z grid values Nrtz : Integer := Nr*Nt*Nz; Ua : Real_Vector(1..Nrtz); Cr : Real_Vector(1..Nr); -- first deriv wrt r Crr : Real_Vector(1..Nr); -- second deriv wrt r Ctt : Real_Vector(1..Nt); -- second deriv wrt t Czz : Real_Vector(1..Nz); -- second deriv wrt z R, T, Z, Ana, Cmp, Err : Real; Urv, Urrv, Uttv, Uzzv : Real; -- computed derivatives Maxerrur, Maxerrurr, Maxerrutt, Maxerruzz, Maxerr : Real := 0.0; X, Y, Tmp : Real; function U(R:Real; T:Real; Z:Real) return Real is begin return R**3 + T**3 + Z**3; end U; function Ur(R:Real; T:Real; Z:Real) return Real is begin return 3.0*R**2; end Ur; function Urr(R:Real; T:Real; Z:Real) return Real is begin return 6.0*R; end Urr; function Ut(R:Real; T:Real; Z:Real) return Real is begin return 3.0*T**2; end Ut; function Utt(R:Real; T:Real; Z:Real) return Real is begin return 6.0*T; end Utt; function Uz(R:Real; T:Real; Z:Real) return Real is begin return 3.0*Z**2; end Uz; function Uzz(R:Real; T:Real; Z:Real) return Real is begin return 6.0*Z; end Uzz; function Uuu(R:Real; T:Real; Z:Real) return Real is begin return 2.0*Ur(R,T,Z)/R + Urr(R,T,Z) + Utt(R,T,Z)/(R*R) + Uzz(R,T,Z); end Uuu; function D2(Urv:Real; Urrv:Real; Uttv:Real; Uzzv:Real; R:Real) return Real is begin return 2.0*Urv/R + Urrv + Uttv/(R*R)+ Uzzv; end D2; function S(I:Integer; J:Integer; K:Integer) return Integer is begin -- index of each grid point return (I-1)*Nt*Nz + (J-1)*Nz + K; end S; begin Put_Line("check_cylinder_deriv.adb starting"); Put_Line("r gird values, not necessarily equally spaced"); for I in 1..Nr loop Rg(I) := Rmin + Hr*Real(I-1); Put_Line("Rg("&Integer'Image(I)&")="&Real'Image(Rg(I))); end loop; Put_Line("theta gird values, not equally spaced"); Tg(1) := 0.0; Tg(2) := 40.0*Pi/180.0; Tg(3) := 85.0*Pi/180.0; Tg(4) := 130.0*Pi/180.0; Tg(5) := 180.0*Pi/180.0; Tg(6) := 240.0*Pi/180.0; Tg(7) := 300.0*Pi/180.0; for J in 1..Nt loop Put_Line("Tg("&Integer'Image(J)&")="&Real'Image(Tg(J))); end loop; Put_Line("z gird values, not necessarily equally spaced"); for K in 1..Nz loop Zg(K) := Zmin+ Hz*Real(K-1); Put_Line("Tg("&Integer'Image(K)&")="&Real'Image(Zg(K))); end loop; Put_Line("analytic partial second derivative"); for I in 1..Nr loop R := Rg(I); for J in 1..Nt loop T := Tg(J); for K in 1..Nz loop Z := Zg(K); Ua(S(I,J,K)) := U(Rg(I),Tg(J),Zg(K)); Put_Line("Ua("&Integer'Image(I)& ","&Integer'Image(J)& ","&Integer'Image(K)& ")="&Real'Image(Ua(S(I,J,K)))& ", D2="&Real'Image(Uuu(Rg(I),Tg(J),Zg(K)))); X := R*Cos(T); Y := R*Sin(T); if abs(Sqrt(X**2+Y**2)-R)>0.0001 then Put_Line("R error"); end if; Tmp := Arctan(Y,X); if Tmp<0.0 then Tmp := 2.0*Pi + Tmp; end if; if abs(Tmp-T)>0.0001 then Put_Line("Theta error"); end if; -- Z = Z end loop; end loop; end loop; for I in 1..Nr loop R := Rg(I); for J in 1..Nt loop T := Tg(J); for K in 1..Nz loop Z := Zg(K); -- compute Ur, Urr, Utt, Uzz and call D2 Nuderiv(1, Nr, I, Rg, Cr); Nuderiv(2, Nr, I, Rg, Crr); Urv := 0.0; Urrv := 0.0; for Ii in 1..Nr loop Urv := Urv + Cr(Ii)*U(Rg(Ii),T,Z); Urrv := Urrv + Crr(Ii)*U(Rg(Ii),T,Z); end loop; Uttv := 0.0; Nuderiv(2, Nt, J, Tg, Ctt); for Jj in 1..Nt loop Uttv := Uttv + Ctt(Jj)*U(R,Tg(Jj),Z); end loop; Uzzv := 0.0; Nuderiv(2, Nz, K, Zg, Czz); for Kk in 1..Nz loop Uzzv := Uzzv + Czz(Kk)*U(R,T,Zg(Kk)); end loop; Cmp := D2(Urv, Urrv, Uttv, Uzzv, R); Ana := Uuu(R,T,Z); Err := Ana - Cmp; Maxerr := Real'Max(Maxerr, Err); Put_Line("i="&Integer'Image(I)& ", j="&Integer'Image(J)& ", k="&Integer'Image(K)); Put_Line("Ana="&Real'Image(Ana)& ", Cmp="&Real'Image(Cmp)& ", err="&Real'Image(Err)); Err := Ur(R,T,Z)-Urv; Maxerrur := Real'Max(Maxerrur, Err); Put_Line("Ur="&Real'Image(Ur(R,T,Z))& ", Urv="&Real'Image(Urv)& ", err="&Real'Image(Err)); Err := Urr(R,T,Z)-Urrv; Maxerrurr := Real'Max(Maxerrurr, Err); Put_Line("Urr="&Real'Image(Urr(R,T,Z))& ", Urrv="&Real'Image(Urrv)& ", err="&Real'Image(Err)); Err := Utt(R,T,Z)-Uttv; Maxerrutt := Real'Max(Maxerrutt, Err); Put_Line("Utt="&Real'Image(Utt(R,T,Z))& ", Uttv="&Real'Image(Uttv)& ", err="&Real'Image(Err)); Err := Uzz(R,T,Z)-Uzzv; Maxerruzz := Real'Max(Maxerruzz, Err); Put_Line("Uzz="&Real'Image(Uzz(R,T,Z))& ", Uzzv="&Real'Image(Uzzv)& ", err="&Real'Image(Err)); Put_Line(" "); end loop; end loop; end loop; Put_Line("Maxerr Ur="&Real'Image(Maxerrur)); Put_Line("Maxerr Urr="&Real'Image(Maxerrurr)); Put_Line("Maxerr Utt="&Real'Image(Maxerrutt)); Put_Line("Maxerr Uzz="&Real'Image(Maxerruzz)); Put_Line("Maxerr deriv="&Real'Image(Maxerr)); Put_line("check_cylinder_deriv.adb finishing"); end check_cylinder_deriv;