-
Notifications
You must be signed in to change notification settings - Fork 37
Expand file tree
/
Copy pathtestinterp3.f90
More file actions
116 lines (90 loc) · 3.05 KB
/
Copy pathtestinterp3.f90
File metadata and controls
116 lines (90 loc) · 3.05 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
program testinterp3
!! Need program statement for FORD
use phys_consts, only: wp,pi
use interpolation
use h5fortran, only: hdf5_file
implicit none
type(hdf5_file) :: hout
integer, parameter :: lx1=80, lx2=90, lx3=100
integer, parameter :: lx1i=256, lx2i=256, lx3i=256
real(wp), parameter :: stride=0.5_wp
real(wp) :: x1(lx1),x2(lx2),x3(lx3),f(lx1,lx2,lx3)
real(wp) :: x1i(lx1i),x2i(lx2i),x3i(lx3i),fi(lx1i,lx2i,lx3i)
real(wp), dimension(1:lx1i*lx2i*lx3i) :: x1ilist,x2ilist,x3ilist,filist
integer :: ix1,ix2,ix3,ik, ierr
!grid for original data
x1=[ ((real(ix1,wp)-1._wp)*stride, ix1=1,lx1) ]
x2=[ ((real(ix2,wp)-1._wp)*stride, ix2=1,lx2) ]
x3=[ ((real(ix3,wp)-1._wp)*stride, ix3=1,lx3) ]
!center grid points at zero
x1=x1-sum(x1)/size(x1,1)
x2=x2-sum(x2)/size(x2,1)
x3=x3-sum(x3)/size(x3,1)
!test function
do ix3=1,lx3
do ix2=1,lx2
do ix1=1,lx1
f(ix1,ix2,ix3)=sin(2._wp*pi/5._wp*x1(ix1))*cos(2._wp*pi/20._wp*x2(ix2))*sin(2._wp*pi/15._wp*x3(ix3))
end do
end do
end do
!> grid for interpolated data
x1i=[ ((real(ix1,wp)-1)*stride/(lx1i/lx1), ix1=1,lx1i) ]
x2i=[ ((real(ix2,wp)-1)*stride/(lx2i/lx2), ix2=1,lx2i) ]
x3i=[ ((real(ix3,wp)-1)*stride/(lx3i/lx3), ix3=1,lx3i) ]
!> center grid points at zero
x1i=x1i-sum(x1i)/size(x1i,1)
x2i=x2i-sum(x2i)/size(x2i,1)
x3i=x3i-sum(x3i)/size(x3i,1)
!> try a 333d interpolation
do ix3=1,lx3i
do ix2=1,lx2i
do ix1=1,lx1i
ik=(ix3-1)*lx1i*lx2i+(ix2-1)*lx1i+ix1
x1ilist(ik)=x1i(ix1)
x2ilist(ik)=x2i(ix2)
x3ilist(ik)=x3i(ix3)
end do
end do
end do
filist=interp3(x1,x2,x3,f,x1ilist,x2ilist,x3ilist)
fi=reshape(filist,[lx1i,lx2i,lx3i])
!> dump results to a file so we can check things
call hout%initialize("input3d.h5", ierr, status="replace", action="write")
if(ierr/=0) error stop 'interp3: could not open input file'
call hout%write("/lx1", lx1, ierr)
call hout%write("/lx2", lx2, ierr)
call hout%write("/lx3", lx3, ierr)
call hout%write("/x1", x1, ierr)
call hout%write("/x2", x2, ierr)
call hout%write("/x3", x3, ierr)
call hout%write("/f", f, ierr)
call hout%finalize(ierr)
if(ierr/=0) error stop 'interp3: could not close input file'
call hout%initialize("output3d.h5", ierr, status="replace", action="write")
if(ierr/=0) error stop 'interp3: could not open output file'
call hout%write("/lx1", lx1i, ierr)
call hout%write("/lx2", lx2i, ierr)
call hout%write("/lx3", lx3i, ierr)
call hout%write("/x1", x1i, ierr)
call hout%write("/x2", x2i, ierr)
call hout%write("/x3", x3i, ierr)
call hout%write("/f", fi, ierr)
call hout%finalize(ierr)
if(ierr/=0) error stop 'interp3: could not close output file'
end program
!> has no problem with > 2GB output files
! block
! integer :: u
! open(newunit=u,file='input3D.dat',status='replace',form='unformatted',access='stream', action='write')
! write(u) lx1,lx2,lx3
! write(u) x1,x2,x3,f
! close(u)
! end block
! block
! integer :: u
! open(newunit=u,file='output3D.dat',status='replace',form='unformatted',access='stream', action='write')
! write(u) lx1i,lx2i,lx3i
! write(u) x1i,x2i,x3i,fi !< since only interpolating in x1
! close(u)
! end block