diff --git a/code/m_calculate.f90 b/code/m_calculate.f90 index cc8fa43..b1cf20e 100644 --- a/code/m_calculate.f90 +++ b/code/m_calculate.f90 @@ -10,6 +10,9 @@ module m_calculate real*8, allocatable, dimension(:,:) :: xsrc real*8, allocatable, dimension(:,:) :: xdst + real*8, allocatable, dimension(:,:) :: rsrc + real*8, allocatable, dimension(:,:) :: rdst + real*8, allocatable, dimension(:,:,:) :: zsrc real*8, allocatable, dimension(:,:,:) :: zdst @@ -31,6 +34,9 @@ contains allocate(xsrc(nyp, nxp), stat=ierr) allocate(xdst(nyp, nxp), stat=ierr) + allocate(rsrc(nxp, nzp), stat=ierr) + allocate(rdst(nxp, nzp), stat=ierr) + allocate(zsrc(nyp, nzp, nxp), stat=ierr) allocate(zdst(nyp, nzp, nxp), stat=ierr) @@ -44,6 +50,9 @@ contains deallocate(xsrc) deallocate(xdst) + deallocate(rsrc) + deallocate(rdst) + deallocate(zsrc) deallocate(zdst) @@ -52,7 +61,8 @@ contains subroutine ddx(dst, src) - real*8, dimension(nxp,nyp,nzp) :: src, dst + real*8, dimension(nxp,nyp,nzp), intent(in) :: src + real*8, dimension(nxp,nyp,nzp), intent(out) :: dst integer :: i, j ,k @@ -72,7 +82,8 @@ contains subroutine ddy(dst, src) - real*8, dimension(nxp,nyp,nzp) :: src, dst + real*8, dimension(nxp,nyp,nzp), intent(in) :: src + real*8, dimension(nxp,nyp,nzp), intent(out) :: dst integer :: i, j ,k @@ -88,7 +99,8 @@ contains subroutine ddz(dst, src) - real*8, dimension(nxp,nyp,nzp) :: src, dst + real*8, dimension(nxp,nyp,nzp), intent(in) :: src + real*8, dimension(nxp,nyp,nzp), intent(out) :: dst integer :: i, j ,k @@ -106,9 +118,38 @@ contains end subroutine ddz + subroutine ddz_ref(dst, src) + + real*8, dimension(nxp,nyp,nzp), intent(in) :: src + real*8, dimension(nxp,nyp,nzp), intent(out) :: dst + + integer :: i, j ,k + + do j = 1,nyp + + do k = 1,nzp + do i = 1,nxp + rsrc(i,k) = src(i,j,k) + end do + end do + + call dfp(nzp, hzp, rsrc, rdst, nxp, 3) + + do k = 1,nzp + do i = 1,nxp + dst(i,j,k) = rdst(i,k) + end do + end do + + end do + + end subroutine ddz_ref + + subroutine d2dx(dst, src) - real*8, dimension(nxp,nyp,nzp) :: src, dst + real*8, dimension(nxp,nyp,nzp), intent(in) :: src + real*8, dimension(nxp,nyp,nzp), intent(out) :: dst integer :: i, j ,k @@ -128,7 +169,8 @@ contains subroutine d2dy(dst, src) - real*8, dimension(nxp,nyp,nzp) :: src, dst + real*8, dimension(nxp,nyp,nzp), intent(in) :: src + real*8, dimension(nxp,nyp,nzp), intent(out) :: dst integer :: i, j ,k @@ -144,7 +186,8 @@ contains subroutine d2dz(dst, src) - real*8, dimension(nxp,nyp,nzp) :: src, dst + real*8, dimension(nxp,nyp,nzp), intent(in) :: src + real*8, dimension(nxp,nyp,nzp), intent(out) :: dst integer :: i, j ,k @@ -161,6 +204,35 @@ contains end subroutine d2dz + + subroutine d2dz_ref(dst, src) + + real*8, dimension(nxp,nyp,nzp), intent(in) :: src + real*8, dimension(nxp,nyp,nzp), intent(out) :: dst + + integer :: i, j ,k + + do j = 1,nyp + + do k = 1,nzp + do i = 1,nxp + rsrc(i,k) = src(i,j,k) + end do + end do + + call d2fp(nzp, hzp, rsrc, rdst, nxp, 3) + + do k = 1,nzp + do i = 1,nxp + dst(i,j,k) = rdst(i,k) + end do + end do + + end do + + end subroutine d2dz_ref + + subroutine tp(a, b, nx) ! a(nb,nx) = transpose(b(nx,nb)) diff --git a/code/makefile b/code/makefile index 2d4c554..596aab9 100644 --- a/code/makefile +++ b/code/makefile @@ -42,7 +42,7 @@ m_terms.f90 : code_gen/code_gen.py code_gen/terms.input code_gen/resources/m_tem ${compiler} -c ${flags} $< clean: - rm -f *.o *.mod x-edge-cold-bc-uPrime-hybrid + rm -f *.o *.mod x-edge-cold-bc-uPrime-hybrid test_calculate testc : test_calculate diff --git a/code/test_calculate.f90 b/code/test_calculate.f90 index e774341..35fec3b 100644 --- a/code/test_calculate.f90 +++ b/code/test_calculate.f90 @@ -10,10 +10,14 @@ real*8, allocatable, dimension(:,:,:) :: aaa,bbb,ccc,ddd integer :: ierr -real*8 :: xx,yy,zz +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 @@ -23,15 +27,18 @@ hyp=l_0*pi/REAL(nyp) hxp=hyp hzp=hyp -allocate(aaa(nxp, nyp, nzp), stat=ierr) -allocate(bbb(nxp, nyp, nzp), stat=ierr) -allocate(ccc(nxp, nyp, nzp), stat=ierr) -allocate(ddd(nxp, nyp, nzp), stat=ierr) +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_calculate_init +write(*,*) "dimension - ", nxp, nyp, nzp +write(*,*) "spacing - ", hxp, hyp, hzp + do k = 1,nzp do j = 1,nyp do i = 1,nxp @@ -40,35 +47,248 @@ xx = i * hxp yy = j * hyp zz = k * hzp -aaa(i,j,k) = sin(2.0 * xx) * sin(3.0 * yy) * sin(4.0 * zz) +aaa(i,j,k) = sin(1.1 * xx) * sin(3.0 * yy) * sin(2.0 * zz) end do end do end do -write(*,*) minval(aaa), maxval(aaa) +write(*,*) "test ddx" + +call ETIME(startt, result) call ddx(bbb, aaa) -write(*,*) minval(bbb), maxval(bbb) +call ETIME(endt, result) +print *, "usert", endt(1) - startt(1) +print *, "systt", endt(2) - startt(2) -call ddy(ccc, aaa) -write(*,*) minval(ccc), maxval(ccc) +do k = 1,nzp +do j = 1,nyp +do i = nxp/4,3*nxp/4 -call ddz(ddd, aaa) -write(*,*) minval(ddd), maxval(ddd) +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(*,*) "test ddy" + +call ETIME(startt, result) + + +call ddy(bbb, aaa) +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 = (3.0 * sin(1.1 * xx) * cos(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(*,*) "test ddz" + +call ETIME(startt, result) + + +call ddz(bbb, aaa) +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 = (2.0 * sin(1.1 * xx) * sin(3.0 * yy) * cos(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(*,*) "test ddz_ref" + +call ETIME(startt, result) + + +call ddz_ref(bbb, aaa) +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 = (2.0 * sin(1.1 * xx) * sin(3.0 * yy) * cos(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(*,*) "test d2dx" + +call ETIME(startt, result) + call d2dx(bbb, aaa) -write(*,*) minval(bbb), maxval(bbb) +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.21 * sin(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(fxyz) > 1.0e-14 )) write(*,*) bbb(i,j,k), fxyz + +end do +end do +end do + +write(*,*) "relerr_min_max", minval(ccc), maxval(ccc) + + + +write(*,*) "test d2dy" + +call ETIME(startt, result) + + +call d2dy(bbb, aaa) +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 = -(9.0 * sin(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(fxyz) > 1.0e-14 )) write(*,*) bbb(i,j,k), fxyz +end do +end do +end do + +write(*,*) "relerr_min_max", minval(ccc), maxval(ccc) + + +write(*,*) "test d2dz" + +call ETIME(startt, result) + + +call d2dz(bbb, aaa) +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 = -(4.0 * sin(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(fxyz) > 1.0e-14 )) write(*,*) bbb(i,j,k), fxyz +end do +end do +end do + +write(*,*) "relerr_min_max", minval(ccc), maxval(ccc) + + +write(*,*) "test d2dz_ref" + +call ETIME(startt, result) + + +call d2dz_ref(bbb, aaa) +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 = -(4.0 * sin(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(fxyz) > 1.0e-14 )) write(*,*) bbb(i,j,k), fxyz +end do +end do +end do + +write(*,*) "relerr_min_max", minval(ccc), maxval(ccc) -call d2dy(ccc, aaa) -write(*,*) minval(ccc), maxval(ccc) -call d2dz(ddd, aaa) -write(*,*) minval(ddd), maxval(ddd)