From 16743238165295c432beb675bfaa5d81d5f2e6dd Mon Sep 17 00:00:00 2001 From: ignis Date: Tue, 6 Aug 2019 05:00:18 +0900 Subject: [PATCH] convert compact.f90 float literals double precision and separte rhs and tdslv and added python code validataions to pycompact.py --- code/Compact.f90 | 38 ++++--- code/pycompact/pycompact.py | 200 ++++++++++++++++++++++++++++++++---- 2 files changed, 207 insertions(+), 31 deletions(-) diff --git a/code/Compact.f90 b/code/Compact.f90 index 0fd1a9e..9744fcb 100644 --- a/code/Compact.f90 +++ b/code/Compact.f90 @@ -235,19 +235,21 @@ SUBROUTINE stdlu(a,n,l) INTEGER :: n - REAL(KIND=8) :: a(n),l(n) + REAL(KIND=8), INTENT(IN) :: a(n) + REAL(KIND=8), INTENT(OUT) :: l(n) REAL(KIND=8) :: d INTEGER :: i - l(1)=1.0/a(1) + l(1)=1.0d0/a(1) DO i=2,n d=a(i)-l(i-1) - l(i)=1.0/d + l(i)=1.0d0/d ENDDO END SUBROUTINE stdlu SUBROUTINE ptdlu(a,n,l,w) INTEGER :: n - REAL(KIND=8) :: a,l(n),w(n) + REAL(KIND=8), INTENT(OUT) :: a + REAL(KIND=8), INTENT(OUT) :: l(n),w(n) INTEGER :: i REAL(KIND=8) :: aa(n),d @@ -271,21 +273,19 @@ l(n)=1./d END SUBROUTINE ptdlu - SUBROUTINE dfnonp(n,h,x,dx,nd,dir) - INTEGER,INTENT(IN) :: n,nd,dir + + SUBROUTINE rhs1np(n,h,x,dx,nd) + INTEGER,INTENT(IN) :: n,nd REAL(KIND=8),INTENT(IN) :: h REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx INTEGER :: i,j REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4 + h1=1.d0/h - ! print *, "dfnonp received (nd,n)", nd, n - - h1=1./h - - r1=7./3. - r2=1./12. + r1=7.d0/3.d0 + r2=1.d0/12.d0 r3=3. a=-1.25 b=1. @@ -311,6 +311,20 @@ dx(j,i)=h1*(r1*t1+r2*t2) ENDDO ENDDO + + END SUBROUTINE rhs1np + + + SUBROUTINE dfnonp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: n,nd,dir + REAL(KIND=8),INTENT(IN) :: h + REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x + REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx + INTEGER :: i,j + REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4 + + CALL rhs1np (n,h,x,dx,nd) + IF (dir.eq.1) CALL tdslv(dx,n,lxf,nd) ! x-direction IF (dir.eq.2) CALL tdslv(dx,n,lyf,nd) ! y-direction IF (dir.eq.3) CALL tdslv(dx,n,lzf,nd) ! z-direction diff --git a/code/pycompact/pycompact.py b/code/pycompact/pycompact.py index b9d6f3f..9dcc7d4 100644 --- a/code/pycompact/pycompact.py +++ b/code/pycompact/pycompact.py @@ -9,22 +9,6 @@ class CompactScheme: self.shape = (nz, ny, nx) - self.coefx1 = np.zeros((nx)) - self.coefy1 = np.zeros((ny)) - self.coefz1 = np.zeros((nz)) - - self.coefx2 = np.zeros((nx)) - self.coefy2 = np.zeros((ny)) - self.coefz2 = np.zeros((nz)) - - self.coefx3 = np.zeros((nx)) - self.coefy3 = np.zeros((ny)) - self.coefz3 = np.zeros((nz)) - - self.coefx4 = np.zeros((nx)) - self.coefy4 = np.zeros((ny)) - self.coefz4 = np.zeros((nz)) - self.px = px self.py = py self.pz = pz @@ -34,6 +18,7 @@ class CompactScheme: self.hy = h self.hz = h + # Allocate LU compact.lxf = np.zeros(nx, dtype=np.float64) compact.lxs = np.zeros(nx, dtype=np.float64) compact.wxf = np.zeros(nx, dtype=np.float64) @@ -55,6 +40,163 @@ class CompactScheme: compact.ludcmp_calculate(nx, ny, nz, bcx, bcy, bcz) + + def test_ludcmp (self): + import pprint + pp = pprint.PrettyPrinter(indent=4) + + # First Derivative Non-periodic BC + l1 = compact.test_nonp_lud1(self.shape[-1]) + + # Second Derivative Non-periodic BC + l2 = compact.test_nonp_lud2(self.shape[-1]) + + print ("Test Internally Calculated Non-periodic Coefs") + print (np.linalg.norm((l1 - compact.lxf)/compact.lxf)) + print (np.linalg.norm((l2 - compact.lxs)/compact.lxs)) + + + def py_rhs_1_np (self, x): + dx = np.zeros(x.shape) + + h1 = 1./self.hx + + r1 = 7./3. + r2 = 1./12. + r3 = 3. + a = -1.25 + b = 1. + c = 0.25 + + nd, n = x.shape + + dx[:, -2] = x[:, -1] - x[:, -3] + dx[:, -1] = - (a*x[:, -1] + b*x[:, -2] + c*x[:, -3]) + dx[:, 0] = (a*x[:, 0] + b*x[:, 1] + c*x[:, 2]) + dx[:, 1] = x[:, 2] - x[:, 0] + + dx[:,-2] = dx[:,-2]*h1*r3 + dx[:,-1] = dx[:,-1]*h1 + dx[:,0] = dx[:,0]*h1 + dx[:,1] = dx[:,1]*h1*r3 + + for i in range(2,n-2): + t1=x[:,i+1]-x[:,i-1] + t2=x[:,i+2]-x[:,i-2] + dx[:,i]=h1*(r1*t1+r2*t2) + + return dx + + + def py_tdslv(self, r, l): + nd, n = r.shape + + r[:,0] = r[:,0] * l[0] + + for i in range(1,n): + r[:,i] = l[i] * (r[:,i] - r[:,i-1]) + + for i in range(n-1)[::-1]: + r[:,i] = r[:,i] - l[i] * r[:,i+1] + + + def test_dfnonp (self): + + x = np.sin(1.1 * np.arange(512) * self.hx).reshape((1,-1)) + + exact = 1.1 * np.cos(1.1 * np.arange(512) * self.hx).reshape((1,-1)) + + + + print ("First Non-periodic RHS Test") + + dx = self.py_rhs_1_np(x) + + dx_fortran = compact.rhs1np(self.hx, x) + + print (np.linalg.norm((dx - dx_fortran) / dx_fortran)) + print (((dx - dx_fortran) / dx_fortran).min()) + print (((dx - dx_fortran) / dx_fortran).max()) + + + + + print ("First Non-periodic TD SOLVE Test") + + l1 = compact.test_nonp_lud1(512) + + self.py_tdslv(dx, l1) + + compact.tdslv(dx_fortran,l1) + + print ("dx - exact") + print (np.linalg.norm((dx - exact) / exact)) + + print (((dx - exact) / exact).min()) + print (((dx - exact) / exact).max()) + + print ("dx_fortran - exact") + print (np.linalg.norm((dx_fortran - exact) / exact)) + + print (((dx_fortran - exact) / exact).min()) + print (((dx_fortran - exact) / exact).max()) + + ''' + import pprint + pp = pprint.PrettyPrinter(indent=4) + + pp.pprint ((dx - exact) / exact) + + pp.pprint ((zip ((dx - exact).ravel(), dx.ravel(), exact.ravel()))) + ''' + + def verify_nonp_lud1(self): + + print ("Non-periodic coef first derivative") + + nx = 512 + aa = np.ones(nx) * 3. + aa[0] = 0.5 + aa[1] = 4. + aa[-2] = 4. + aa[-1] = 0.5 + + coef = compact.stdlu(aa) + + coef_verify = self.py_stdlu(aa) + + print (np.linalg.norm((coef - coef_verify)/coef_verify)) + + + def verify_nonp_lud2(self): + + print ("Non-periodic coef second derivative") + + nx = 512 + aa = np.ones(nx) * 3. + aa[0] = 2./11. + aa[1] = 10. + aa[-2] = 10. + aa[-1] = 2./11. + + coef = compact.stdlu(aa) + + coef_verify = self.py_stdlu(aa) + + print (np.linalg.norm((coef - coef_verify)/coef_verify)) + + + def py_stdlu(self, aa): + coef = np.ones(aa.shape)/aa[0] + + print ("coef.size = ", coef.size) + + for i in range(1,coef.size): + coef[i]=1.0/(aa[i]-coef[i-1]) + + return coef + + def ddx (self, src): if src.shape != self.shape: @@ -223,8 +365,6 @@ def validate_trigonometric(): pi8 = np.arccos(-1.) - print ("calculated pi value = ", pi8) - l_0 = 2.0 hyp=l_0*pi8/ny hxp=hyp @@ -233,16 +373,38 @@ def validate_trigonometric(): cs = CompactScheme(nx, ny, nz, False, True, True, 4., 2., 2.) + cs.test_ludcmp() + + cs.verify_nonp_lud1() + + cs.verify_nonp_lud2() + + cs.test_dfnonp() + + + + print ("Test ddx") Y1 = np.zeros(shape) true = np.zeros(shape) - XX = np.arange(nx) * hxp YY = np.arange(ny) * hyp ZZ = np.arange(nz) * hzp + + + print ("1-D sine test") + cos_fortran = compact.dfnonp(hxp, np.sin(1.1*XX).reshape((1,-1)), 1) + cos_exact = 1.1 * np.cos(1.1*XX).reshape((1,-1)) + print (np.linalg.norm((cos_fortran - cos_exact)/cos_exact)) + # print (((cos_fortran - cos_exact)/cos_exact)) + + + + print ("3-D trigonometric test") + zz, yy, xx = np.meshgrid(ZZ, YY, XX) Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:]