PageRenderTime 206ms CodeModel.GetById 104ms RepoModel.GetById 1ms app.codeStats 0ms

/tools/flat_table_2D.f90

https://bitbucket.org/NeilMiller/h5tbuilder
FORTRAN Modern | 401 lines | 304 code | 71 blank | 26 comment | 36 complexity | 5e568400da773267598c783436bd8bea MD5 | raw file
  1. module table_2D_mod
  2. use flat_interpolator_mod
  3. implicit none
  4. integer, private, parameter :: NDIM = 2
  5. type :: table_2D
  6. private
  7. character(LEN=80) :: name
  8. logical :: initialized = .false.
  9. integer :: nderiv
  10. integer :: nmax1
  11. integer :: nmax2
  12. integer, dimension(NDIM) :: nx
  13. logical, dimension(NDIM) :: log_spc
  14. real(8), dimension(NDIM) :: x0
  15. real(8), dimension(NDIM) :: xmax
  16. real(8), dimension(NDIM) :: dx
  17. logical, dimension(NDIM) :: zero = .false. ! Allows for zero point on logarithmic grids
  18. ! If true, the first position in the table
  19. ! along a dimension corresponds to zero and
  20. ! the second point has the value x0 (instead
  21. ! of the first)
  22. real(8), dimension(:,:), allocatable :: tab
  23. real(8), dimension(:,:), allocatable :: xx
  24. type(flat_interpolator) :: interpol
  25. contains
  26. procedure :: init => init_table_2D
  27. procedure :: get_xmin => get_xmin_table_2D
  28. procedure :: get_xmax => get_xmax_table_2D
  29. procedure :: load_slab_der => load_slab_der_table_2D
  30. procedure :: load_slab_pos => load_slab_pos_table_2D
  31. procedure :: interp => interp_table_2D
  32. end type table_2D
  33. private :: init_table_2D
  34. private :: get_xmin_table_2D
  35. private :: get_xmax_table_2D
  36. private :: load_slab_der_table_2D
  37. private :: load_slab_pos_table_2D
  38. private :: interp_table_2D
  39. private :: table_position
  40. contains
  41. function get_xmin_table_2D(this,id) result(xmin)
  42. real(8) :: xmin
  43. class(table_2D), intent(in) :: this
  44. integer, intent(in) :: id
  45. if (id<=NDIM.and.id>0) then
  46. xmin = this%x0(id)
  47. else
  48. xmin = this%x0(NDIM)
  49. endif
  50. return
  51. end function get_xmin_table_2D
  52. function get_xmax_table_2D(this,id) result(xmax)
  53. real(8) :: xmax
  54. class(table_2D), intent(in) :: this
  55. integer, intent(in) :: id
  56. if (id<=NDIM.and.id>0) then
  57. xmax = this%xmax(id)
  58. else
  59. xmax = this%xmax(NDIM)
  60. endif
  61. return
  62. end function get_xmax_table_2D
  63. subroutine init_table_2D(this,nx,nderiv,x0,dx,log_spc,tab,zero)
  64. class(table_2D), intent(out) :: this
  65. integer, dimension(NDIM), intent(in) :: nx
  66. integer, dimension(NDIM), intent(in) :: nderiv
  67. real(8), dimension(NDIM), intent(in) :: x0,dx
  68. logical, dimension(NDIM), intent(in), optional :: log_spc
  69. logical, dimension(NDIM), intent(in), optional :: zero
  70. real(8), dimension(:,:,0:,0:), intent(in), optional :: tab
  71. integer :: ierr,i,j,i1,i2,i3,i4
  72. ! Set the array properties
  73. this%nx = nx
  74. this%x0 = x0
  75. this%dx = dx
  76. this%nderiv = MINVAL(nderiv(1:NDIM))
  77. if (present(zero).and.present(log_spc)) then
  78. do i=1,NDIM
  79. this%zero(i) = zero(i) .and. log_spc(i)
  80. this%log_spc(i) = log_spc(i)
  81. enddo
  82. elseif (present(log_spc)) then
  83. this%log_spc = log_spc
  84. this%zero = .false.
  85. else
  86. this%zero = .false.
  87. this%log_spc = .false.
  88. this%xmax = x0 + dx*(nx-1.d0)
  89. endif
  90. ! Build tables of point values
  91. ! This is done only to prevent numerous
  92. ! recalculations of 10.d0**x for
  93. ! logarithmic tables (i.e. a hash is used
  94. ! to find the table positions, not a table search).
  95. if (allocated(this%xx)) deallocate(this%xx)
  96. print *,1
  97. allocate(this%xx(MAXVAL(this%nx),NDIM))
  98. print *,2
  99. do i=1,NDIM
  100. if (this%log_spc(i).and.this%zero(i)) then
  101. this%xx(1,i) = 0.d0
  102. do j=2,this%nx(i)
  103. this%xx(j,i) = this%x0(i)*10.d0**(REAL(j-2)*this%dx(i))
  104. enddo
  105. elseif (this%log_spc(i)) then
  106. do j=1,this%nx(i)
  107. this%xx(j,i) = this%x0(i)*10.d0**(REAL(j-1)*this%dx(i))
  108. enddo
  109. else
  110. do j=1,this%nx(i)
  111. this%xx(j,i) = this%x0(i) + REAL(j-1)*this%dx(i)
  112. enddo
  113. endif
  114. enddo
  115. ! Set the flat array sizes
  116. this%nmax1 = (this%nderiv+1)**NDIM-1
  117. this%nmax2 = 1
  118. do i=1,NDIM
  119. this%nmax2 = this%nx(i)*this%nmax2
  120. enddo
  121. this%nmax2 = this%nmax2-1
  122. ! Allocate the table
  123. if (allocated(this%tab)) deallocate(this%tab)
  124. allocate(this%tab(0:this%nmax1,0:this%nmax2))
  125. this%initialized = .true.
  126. ! Load table if present
  127. if (present(tab)) then
  128. do i1=0,this%nderiv
  129. do i2=0,this%nderiv
  130. ierr = this%load_slab_der(tab(:,:,i1,i2),(/i1,i2/))
  131. if (ierr .ne. 0) print *, 'Table load error ',ierr
  132. enddo
  133. enddo
  134. endif
  135. return
  136. end subroutine
  137. function load_slab_der_table_2D(this,dat,ideriv) result(ierr)
  138. integer :: ierr
  139. class(table_2D), intent(inout) :: this
  140. real(8), dimension(:,:), intent(in) :: dat
  141. integer, dimension(:), intent(in) :: ideriv
  142. integer :: i,id,ii,i1,i2,i3,i4
  143. ierr = 0
  144. ! Can't load an uninitialized table
  145. if (.not.this%initialized) then
  146. ierr = -1
  147. return
  148. endif
  149. ! Check that the array is large enough
  150. if (SIZE(dat)<this%nmax2+1) then
  151. ierr = i
  152. return
  153. endif
  154. ! Check derivative position specification
  155. if (SIZE(ideriv)<NDIM) then
  156. ierr = -2
  157. return
  158. endif
  159. do i=1,NDIM
  160. if (ideriv(i)<0 .or. ideriv(i)>this%nderiv) then
  161. ierr = -i-2
  162. return
  163. endif
  164. enddo
  165. ! Put it in the table
  166. id = 0
  167. do i=NDIM,1,-1
  168. id = (this%nderiv+1)*id + ideriv(i)
  169. enddo
  170. do i1=0,this%nx(1)-1
  171. do i2=0,this%nx(2)-1
  172. ii = this%nx(1)*i2 + i1
  173. this%tab(id,ii) = dat(i1+1,i2+1)
  174. enddo
  175. enddo
  176. return
  177. end function load_slab_der_table_2D
  178. function load_slab_pos_table_2D(this,dat,ipos) result(ierr)
  179. integer :: ierr
  180. class(table_2D), intent(inout) :: this
  181. real(8), dimension(:,:), intent(in) :: dat
  182. integer, dimension(:), intent(in) :: ipos
  183. integer :: i,id,ii,i1,i2,i3,i4
  184. ierr = 0
  185. ! Can't load an uninitialized table
  186. if (.not.this%initialized) then
  187. ierr = -1
  188. return
  189. endif
  190. ! Check that the array is large enough
  191. if (SIZE(dat)<this%nmax1+1) then
  192. ierr = i
  193. return
  194. endif
  195. ! Check derivative position specification
  196. if (SIZE(ipos)<NDIM) then
  197. ierr = -2
  198. return
  199. endif
  200. do i=1,NDIM
  201. if (ipos(i)<1 .or. ipos(i)>this%nx(i)) then
  202. ierr = -i-2
  203. return
  204. endif
  205. enddo
  206. ! Put it in the table
  207. ii = 0
  208. do i=NDIM,1,-1
  209. ii = this%nx(i)*id + ipos(i) - 1
  210. enddo
  211. do i1=0,this%nderiv
  212. do i2=0,this%nderiv
  213. id = (this%nderiv+1)*i2 + i1
  214. this%tab(id,ii) = dat(i1+1,i2+1)
  215. enddo
  216. enddo
  217. return
  218. end function load_slab_pos_table_2D
  219. function interp_table_2D(this,x,der,order) result(f)
  220. real(8) :: f
  221. class(table_2D) :: this
  222. real(8) :: x(NDIM),dd(NDIM),delta(NDIM)
  223. integer :: idx(NDIM)
  224. integer, optional :: der(NDIM)
  225. integer, optional :: order
  226. integer :: ierr,ipos(0:2**NDIM-1)
  227. integer :: ii,i1,i2,i3,i4
  228. if (present(order)) then
  229. print *,'Not allowing you to choose the order for now'
  230. endif
  231. if (this%nderiv==0) then
  232. ierr = this%interpol%set_order('linear')
  233. elseif (this%nderiv==1) then
  234. ierr = this%interpol%set_order('cubic')
  235. else
  236. ierr = this%interpol%set_order('quintic')
  237. endif
  238. call table_position(this,x,dd,idx,delta)
  239. ! write(*,*) "luketbl: dd = ", dd
  240. ! write(*,*) "luketbl: delta = ", delta
  241. do i1=0,1
  242. do i2=0,1
  243. ii = 2*i2 + i1
  244. ipos(ii) = this%nx(1)*(idx(2)-1 + i2) &
  245. + (idx(1)-1 + i1)
  246. enddo
  247. enddo
  248. if (present(der)) then
  249. f = this%interpol%f2D(dd,delta,this%tab(:,ipos),der,this%nderiv+1)
  250. else
  251. f = this%interpol%f2D(dd,delta,this%tab(:,ipos),(/0,0/),this%nderiv+1)
  252. endif
  253. return
  254. end function interp_table_2D
  255. subroutine table_position(this,x,dd,idx,delta)
  256. class(table_2D), intent(in) :: this
  257. real(8), dimension(NDIM), intent(in) :: x
  258. real(8), dimension(NDIM), intent(out) :: dd
  259. real(8), dimension(NDIM), intent(out), optional :: delta
  260. integer, dimension(NDIM), intent(out) :: idx
  261. real(8), dimension(NDIM) :: xt
  262. integer :: i
  263. real(8) :: xl,xh,x0
  264. integer, save :: nprint = 0
  265. ! Transform logarithmic dimensions
  266. do i=1,NDIM
  267. if (this%log_spc(i)) then
  268. xt(i) = log10(x(i))
  269. else
  270. xt(i) = x(i)
  271. endif
  272. enddo
  273. ! Hash into table
  274. do i=1,NDIM
  275. if (this%log_spc(i)) then
  276. x0 = log10(this%x0(i))
  277. else
  278. x0 = this%x0(i)
  279. endif
  280. if (xt(i)<x0) then
  281. idx(i) = 1
  282. if (this%log_spc(i) .and. this%zero(i)) then
  283. continue
  284. elseif (this%log_spc(i)) then
  285. if (nprint<20) write(6,'(a,i2,a,2es12.3)')'Value ',i,' off bottom of table.',10.d0**xt(i),10.d0**x0
  286. else
  287. if (nprint<20) write(6,'(a,i2,a,2es12.3)')'Value ',i,' off bottom of table.',xt(i),x0
  288. endif
  289. nprint = nprint + 1
  290. elseif (xt(i)==(x0 + (this%nx(i)-1)*this%dx(i))) then
  291. idx(i) = this%nx(i)-1
  292. elseif (xt(i)>(x0 + (this%nx(i)-1)*this%dx(i))) then
  293. idx(i) = this%nx(i)-1
  294. if (nprint<20) write(6,'(a,i2,a,2f5.1,i3)')'Value ',i,' off top of table.',xt(i),this%dx(i),this%nx(i)
  295. nprint = nprint + 1
  296. elseif (xt(i).ne.xt(i)) then
  297. idx(i) = 1
  298. if (nprint<20) write(6,'(a,i2,a,es12.3)')'Value ',i,' is NaN.',xt(i)
  299. nprint = nprint + 1
  300. else
  301. idx(i) = INT(FLOOR((xt(i)-x0)/this%dx(i))) + 1
  302. endif
  303. enddo
  304. ! Calculate fractional distances
  305. do i=1,NDIM
  306. if (idx(i)==1 .and. this%log_spc(i) .and. this%zero(i)) then
  307. dd(i) = x(i)/this%x0(i)
  308. if (present(delta)) delta(i) = this%x0(i)
  309. elseif (this%log_spc(i) .and. this%zero(i)) then
  310. idx(i) = idx(i) + 1 ! Set to the correct grid position
  311. xl = this%xx(idx(i),i)
  312. xh = this%xx(idx(i)+1,i)
  313. dd(i) = (x(i)-xl)/(xh-xl)
  314. if (present(delta)) delta(i) = xh-xl
  315. elseif (this%log_spc(i)) then
  316. xl = this%xx(idx(i),i)
  317. xh = this%xx(idx(i)+1,i)
  318. dd(i) = (x(i)-xl)/(xh-xl)
  319. if (present(delta)) delta(i) = xh-xl
  320. else
  321. dd(i) = (x(i)-(idx(i)-1)*this%dx(i) - this%x0(i))/this%dx(i)
  322. if (present(delta)) delta(i) = this%dx(i)
  323. endif
  324. if (dd(i).ne.dd(i)) dd(i) = 0.d0
  325. dd(i) = min(dd(i),1.0)
  326. dd(i) = max(dd(i),0.0)
  327. enddo
  328. return
  329. end subroutine table_position
  330. end module table_2D_mod