From b6084a48b0d3943d6b3920e31cc8ae16c27c595f Mon Sep 17 00:00:00 2001
From: Yeongdo Park <ignis@postech.ac.kr>
Date: Thu, 8 Aug 2019 12:32:23 +0900
Subject: [PATCH] test_compact and pycompact results agree

---
 code/pycompact/pycompact.py |  5 +++
 code/test_compact.f90       | 71 ++++++++++++++++++++++++++++++++-----
 2 files changed, 68 insertions(+), 8 deletions(-)

diff --git a/code/pycompact/pycompact.py b/code/pycompact/pycompact.py
index 9dcc7d4..73e704a 100644
--- a/code/pycompact/pycompact.py
+++ b/code/pycompact/pycompact.py
@@ -398,6 +398,11 @@ def validate_trigonometric():
     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 (cos_fortran.min(), cos_fortran.max())
+
+    print (cos_exact.min(), cos_exact.max())
+
     print (np.linalg.norm((cos_fortran - cos_exact)/cos_exact))
     # print (((cos_fortran - cos_exact)/cos_exact))
 
diff --git a/code/test_compact.f90 b/code/test_compact.f90
index f16ade2..71b53eb 100644
--- a/code/test_compact.f90
+++ b/code/test_compact.f90
@@ -4,8 +4,8 @@ use Compact
 
 implicit none
 
-real*8, parameter :: pi=3.14159265358979323846
-real*8, parameter :: me=1.00e-20
+real*8, parameter :: pi=3.14159265358979323846d0
+real*8, parameter :: me=1.00d-20
 
 integer :: nxp,nyp,nzp
 
@@ -15,7 +15,7 @@ real*8 :: l_0
 
 real*8, allocatable, dimension(:,:,:) :: aaa,bbb,ccc,ddd
 
-real*8, allocatable, dimension(:,:) :: aa,bb,cc,dd
+real*8, allocatable, dimension(:,:) :: aa,bb,cc,dd,ee
 
 integer :: ierr
 
@@ -31,7 +31,7 @@ nxp = 512
 nyp = 256
 nzp = 256
 
-l_0 = 2.0
+l_0 = 2.0d0
 hyp=l_0*pi/REAL(nyp)
 hxp=hyp
 hzp=hyp
@@ -40,6 +40,7 @@ allocate(aa(1, nxp), stat=ierr); aa=0.0
 allocate(bb(1, nxp), stat=ierr); bb=0.0
 allocate(cc(1, nxp), stat=ierr); cc=0.0
 allocate(dd(1, nxp), stat=ierr); dd=0.0
+allocate(ee(1, nxp), stat=ierr); ee=0.0
 
 
 allocate(aaa(nxp, nyp, nzp), stat=ierr); aaa=0.0
@@ -51,9 +52,9 @@ allocate(ddd(nxp, nyp, nzp), stat=ierr); ddd=0.0
 call ludcmp(nxp, nyp, nzp, 1, 0, 0)
 
 do i = 1,nxp
-xx = i * hxp
-aa(1,i) = sin(1.1 * xx)
-dd(1,i) = 1.1 * cos(1.1 * xx)
+xx = dfloat(i-1) * hxp
+aa(1,i) = sin(1.1d0 * xx)
+dd(1,i) = 1.1d0 * cos(1.1d0 * xx)
 end do
 
 call rhs1np(nxp, hxp, aa, bb, 1)
@@ -73,8 +74,61 @@ print *, cc
 
 print *, dd
 
+ee = ((cc - dd)/dd)
+
+
+write(*,*) "relerr_min_max", minval(ee(1,nxp/4:3*nxp/4)), maxval(ee(1,nxp/4:3*nxp/4))
+write(*,*) "relerr_min_max", minval(ee), maxval(ee)
+write(*,*) "  exact_min_max", minval(dd), maxval(dd)
+write(*,*) "compact_min_max", minval(cc), maxval(cc)
+
+
+write(*,*) "test ddx"
+
+do k = 1,nzp
+do j = 1,nyp
+do i = 1,nxp
+
+xx = i * hxp
+yy = j * hyp
+zz = k * hzp
+
+aaa(i,j,k) = sin(1.1 * xx) * sin(3.0 * yy) * sin(2.0 * zz)
+
+end do
+end do
+end do
+
+
+
+
+
+do k = 1,nzp
+do j = 1,nyp
+
+call dfnonp(nxp, hxp, aaa(:,j,k), bbb(:,j,k), 1, 1)
+
+end do
+end do
+
+do k = 1,nzp
+do j = 1,nyp
+do i = 1,nxp ! nxp/4,3*nxp/4
+
+xx = i * hxp
+yy = j * hyp
+zz = k * hzp
+fxyz = (1.1 * cos(1.1 * xx) * sin(3.0 * yy) * sin(2.0 * zz))
+ccc(i,j,k) = (bbb(i,j,k) - fxyz) / (fxyz)
+
+if ((abs(ccc(i,j,k)) > 0.01 ) .and. (abs(bbb(i,j,k)) > 1.0e-14 )) write(*,*) bbb(i,j,k), fxyz
+
+end do
+end do
+end do
+
+write(*,*) "relerr_min_max", minval(ccc), maxval(ccc)
 
-write(*,*) "relerr_min_max", minval((cc - dd)/dd), maxval(((cc - dd)/dd))
 
 
 
@@ -85,6 +139,7 @@ deallocate(aa)
 deallocate(bb)
 deallocate(cc)
 deallocate(dd)
+deallocate(ee)
 
 deallocate(aaa)
 deallocate(bbb)