From 3f2034f0dbb0d1637502f06d4368a61eae8706c4 Mon Sep 17 00:00:00 2001 From: ignis Date: Sat, 13 Mar 2021 07:41:20 +0900 Subject: [PATCH] pycompact - added dns-field test --- code/makefile | 3 + code/pycompact/pycompact.py | 16 ++- code/test_field.f90 | 225 ++++++++++++++++++++++++++++++++++++ 3 files changed, 240 insertions(+), 4 deletions(-) create mode 100644 code/test_field.f90 diff --git a/code/makefile b/code/makefile index e3edb3c..8eaaad9 100644 --- a/code/makefile +++ b/code/makefile @@ -61,6 +61,9 @@ clean: test : test_calculate ./test_calculate +test_field : test_field.o m_openmpi.o m_parameters.o Compact.o m_calculate.o m_arrays.o + ${compiler} -o test_field ${flags} m_openmpi.o m_parameters.o Compact.o m_calculate.o m_arrays.o test_field.o + test_calculate : test_calculate.o ${compiler} -o test_calculate ${flags} m_openmpi.o m_parameters.o Compact.o m_calculate.o test_calculate.o diff --git a/code/pycompact/pycompact.py b/code/pycompact/pycompact.py index dae1419..0b46dab 100644 --- a/code/pycompact/pycompact.py +++ b/code/pycompact/pycompact.py @@ -497,24 +497,32 @@ def test_dns_data(): file_name = "./fort.1000" - answer = "export-ddxc.dat" + answer = "./fort.2000" t, nx, ny, nz, u, v, w, Y0, Y1 = read_data(file_name) - ddx_answer = np.fromfile(answer, dtype=np.float64, count=Y1.size).reshape(Y1.shape) + with open(answer, 'rb') as ans_file: + ans_file.seek(4) + ddx_answer = np.fromfile(ans_file, dtype=np.float64, count=Y1.size, ).reshape(Y1.shape) cs = CompactScheme(nx, ny, nz, False, True, True, 4, 2, 2) - ddx = cs.ddx(1.-Y1) + ddx = cs.ddx(Y1) print (ddx.min(), ddx.max()) print (ddx_answer.min(), ddx_answer.max()) + relerr = (ddx - ddx_answer) / (ddx_answer) + err = (ddx - ddx_answer) + print ("Absolute Error", np.nanmin(err), np.nanmax(err)) + print ("Relative Error", np.nanmin(relerr), np.nanmax(relerr)) if __name__ == "__main__": - validate_trigonometric() + # validate_trigonometric() + + print("DNS Field Test") test_dns_data() diff --git a/code/test_field.f90 b/code/test_field.f90 new file mode 100644 index 0000000..4425ade --- /dev/null +++ b/code/test_field.f90 @@ -0,0 +1,225 @@ +program test_calculate + +use m_parameters + +use m_calculate +use m_arrays + +implicit none + +real*8, allocatable, dimension(:,:,:) :: aaa,bbb,ccc,ddd + +integer :: ierr + +real*8 :: xx,yy,zz,fxyz + +integer :: i,j,k + +real*4, dimension(2) :: startt +real*4, dimension(2) :: endt +real*4 :: result + +nxp = 512 +nyp = 256 +nzp = 256 + +l_0 = 2.0 +hyp=l_0*pi/REAL(nyp) +hxp=hyp +hzp=hyp + +allocate(aaa(nxp, nyp, nzp), stat=ierr); aaa=0.0 +allocate(bbb(nxp, nyp, nzp), stat=ierr); bbb=0.0 +allocate(ccc(nxp, nyp, nzp), stat=ierr); ccc=0.0 +allocate(ddd(nxp, nyp, nzp), stat=ierr); ddd=0.0 + +call m_arrays_init +call m_calculate_init + +call READ_FILE(1000) + + +write(*,*) "dimension - ", nxp, nyp, nzp +write(*,*) "spacing - ", hxp, hyp, hzp + +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 + +write(*,*) "test ddx" + +call ETIME(startt, result) + + +call ddx(bbb, y) +call ETIME(endt, result) +print *, "usert", endt(1) - startt(1) +print *, "systt", endt(2) - startt(2) + +do k = 1,nzp +do j = 1,nyp +do i = 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) + +OPEN(2000,FORM='unformatted',STATUS='unknown') +write(2000) bbb +CLOSE(2000) + +call m_calculate_finalize + + +deallocate(aaa) +deallocate(bbb) +deallocate(ccc) +deallocate(ddd) + +end program + +!=============================================================================== + +SUBROUTINE READ_FILE(num) + +use m_parameters +use m_arrays + +REAL*8 :: old_scalar(2,nxp,nyp,nzp) +REAL*8 :: new_scalar(nxp,nyp,nzp,2) + + INTEGER, INTENT(IN) :: num + REAL, DIMENSION(2) :: tmpr + INTEGER :: nx, ny, nz + REAL :: tmp1,tmp2 + REAL :: dt,dummyu + INTEGER :: ncyc + + REAL*4 :: e, t(2) + REAL*4 :: e0, t0(2) + REAL*4 :: e1, t1(2) + + e0 = etime(t0) + + OPEN(num,FORM='unformatted',STATUS='unknown') + + READ (num) tnow,nx,ny,nz,tmp1,tmp2 + + IF ((nx .ne. nxp) .or. (ny .ne. nyp) .or. (nz .ne. nzp)) THEN + WRITE(0,*) "Array dimension mismatch", nx, ny, nz, " != ", nxp, nyp, nzp + STOP -1 + ENDIF + + READ (num) ncyc,dt,dummyu + READ (num) tmpr(1:2) + READ (num) tmpr(1:2) + READ (num) tmpr(1:2) + + num_=num + + IF(num.le.shiftnum) THEN + ! WRITE(*,*) ' with an old fort data from Nueman-0X' + READ (num) u,v,w,old_scalar + y(:,:,:) = old_scalar(2,:,:,:) + ELSE + ! WRITE(*,*) ' with a new fort data from Comb-Cluster' + READ (num) u,v,w,new_scalar + u = u + dummyu + y = new_scalar(:,:,:,2) + ENDIF + + CLOSE (num) + + e = etime(t) + + WRITE(*,'(a,i2,a,i1,a,i4,a,i4,a,i4,a,f8.2,a,f8.2,a,f8.2)') & + "PID", myid, "|PASS", ipass, "|fort.", num, & + "|", nprogress, "/", endnum-startnum+1, & + "|E:", e, "|U:", t(1), "|S:", t(2) + + + + +END SUBROUTINE READ_FILE + +SUBROUTINE WRITE_FILE(num) + +use m_parameters +use m_arrays + +REAL*8 :: old_scalar(2,nxp,nyp,nzp) +REAL*8 :: new_scalar(nxp,nyp,nzp,2) + + INTEGER, INTENT(IN) :: num + REAL, DIMENSION(2) :: tmpr + INTEGER :: nx, ny, nz + REAL :: tmp1,tmp2 + REAL :: dt,dummyu + INTEGER :: ncyc + + REAL*4 :: e, t(2) + REAL*4 :: e0, t0(2) + REAL*4 :: e1, t1(2) + + e0 = etime(t0) + + OPEN(num,FORM='unformatted',STATUS='unknown') + + READ (num) tnow,nx,ny,nz,tmp1,tmp2 + + IF ((nx .ne. nxp) .or. (ny .ne. nyp) .or. (nz .ne. nzp)) THEN + WRITE(0,*) "Array dimension mismatch", nx, ny, nz, " != ", nxp, nyp, nzp + STOP -1 + ENDIF + + READ (num) ncyc,dt,dummyu + READ (num) tmpr(1:2) + READ (num) tmpr(1:2) + READ (num) tmpr(1:2) + + num_=num + + IF(num.le.shiftnum) THEN + ! WRITE(*,*) ' with an old fort data from Nueman-0X' + READ (num) u,v,w,old_scalar + y(:,:,:) = old_scalar(2,:,:,:) + ELSE + ! WRITE(*,*) ' with a new fort data from Comb-Cluster' + READ (num) u,v,w,new_scalar + u = u + dummyu + y = new_scalar(:,:,:,2) + ENDIF + + CLOSE (num) + + e = etime(t) + + WRITE(*,'(a,i2,a,i1,a,i4,a,i4,a,i4,a,f8.2,a,f8.2,a,f8.2)') & + "PID", myid, "|PASS", ipass, "|fort.", num, & + "|", nprogress, "/", endnum-startnum+1, & + "|E:", e, "|U:", t(1), "|S:", t(2) + + + + +END SUBROUTINE WRITE_FILE