DISPATCH
interpolate_mod.f90
1 !===============================================================================
2 !===============================================================================
3 MODULE interpolate_mod
4  USE io_mod
5  USE trace_mod
6  USE mpi_mod
7  USE kinds_mod
8  USE index_mod
9  implicit none
10  private
11 
12  ! a generic interface for interpolation functions.
13  abstract interface
14  function interpolator_interface (q,i,j,k,w) result(qtilde)
15  import kindscalarvar
16  real(kind=KindScalarVar), intent(in):: q(:,:,:)
17  integer, intent(in):: i, j, k
18  real, intent(in):: w(3)
19  real(kind=KindScalarVar):: qtilde
20  end function
21  end interface
22 
23  type:: interpolator_t
24  integer:: order_interpolator
25  procedure(interpolator_interface), pointer, nopass:: interpolator => null()
26  contains
27  procedure, nopass:: trilinear_log
28  procedure, nopass:: trilinear_pv
29  procedure, nopass:: four_d
30  procedure, nopass:: four_d_log
31  procedure, nopass:: four_d_pv
32  end type
33  type(interpolator_t):: interpolator
34 
35 PUBLIC selectinterpolator, interpolator, interpolator_interface, interpolator_unsigned
36 CONTAINS
37 
38 !===============================================================================
39 !> Create and return a procedure pointer for an interpolation function.
40 !===============================================================================
41 FUNCTION selectinterpolator (iorder) result(interp)
42  procedure(interpolator_interface), pointer :: interp
43  integer, intent(in) :: iorder
44  !-----------------------------------------------------------------------------
45  call trace_begin('interpolate_t%SelectInterpolator')
46 
47  select case (iorder)
48  case (0) ! donor cell
49  interp => interpolator_donor_cell
50  interpolator%interpolator => interpolator_donor_cell
51  case (1) ! linear van Leer
52  interp => interpolator_van_leer_3d
53  interpolator%interpolator => interpolator_van_leer_3d
54  case (-1) ! straight, un-limited trilinear interpolation
55  interp => interpolator_trilinear
56  interpolator%interpolator => interpolator_trilinear
57  case default
58  call mpi%abort("The order you have selected has not been implemented. Abort!")
59  end select
60  interpolator%order_interpolator = iorder
61 
62  call trace_end
63 END FUNCTION selectinterpolator
64 
65 !===============================================================================
66 !> Donor cell (direct injection) interpolation.
67 !===============================================================================
68 FUNCTION interpolator_donor_cell (q,i,j,k,w) result(qtilde)
69  real(kind=KindScalarVar), intent(in):: q(:,:,:)
70  integer, intent(in):: i, j, k
71  real, intent(in):: w(3)
72  real(kind=KindScalarVar):: qtilde
73  !-----------------------------------------------------------------------------
74  call trace_begin('interpolate_mod::donor_cell', 2)
75 
76  qtilde = q(i,j,k)
77 
78  call trace_end
79 END FUNCTION interpolator_donor_cell
80 
81 !===============================================================================
82 !> Piece-wise monotonic, conservative linear interpolation (PLI; van Leer 1977).
83 !> aa, bb, cc are the linear interpolation weights in each direction.
84 !> p is the patch we are interpolating *FROM*.
85 !>
86 !> Note that, for an interpolation in 3-D, three 1-D, linear van Leer
87 !> interpolations are *not* guaranteed to be monotonic!
88 !===============================================================================
89 FUNCTION interpolator_van_leer_3d (q,i,j,k,w) result(qtilde)
90  real(kind=KindScalarVar), intent(in):: q(:,:,:)
91  integer, intent(in):: i, j, k
92  real, intent(in):: w(3)
93  real(kind=KindScalarVar):: qtilde
94  !-----------------------------------------------------------------------------
95  logical :: ldim(3), lzc=.false., lfc(3)=.false.
96  real(8) :: qijk, dq1, dq2, dq3, qp, qm, dp, dm, aych1, aych2, aych3
97  integer:: itimer=0
98  !-----------------------------------------------------------------------------
99  call trace_begin('interpolate_mod::van_Leer', 1, itimer=itimer)
100  ldim(:) = .false.
101  if (lzc) ldim(:) = .true.
102  if (lfc(2) .or. lfc(3)) ldim(1) = .true.
103  if (lfc(3) .or. lfc(1)) ldim(2) = .true.
104  if (lfc(1) .or. lfc(2)) ldim(3) = .true.
105  if (size(q,1) <= 1) ldim(1) = .false.
106  if (size(q,2) <= 1) ldim(2) = .false.
107  if (size(q,3) <= 1) ldim(3) = .false.
108 
109  qijk = q(i,j,k)
110 
111  ! --- First, determine the van Leer slopes in each direction (as needed) ---
112  dq1 = 0.0d0
113  dq2 = 0.0d0
114  dq3 = 0.0d0
115 
116  ! 1st dimension
117  if (ldim(1)) then
118  qp = q(i+1,j,k)
119  qm = q(i-1,j,k)
120  dp = qp - qijk
121  dm = qijk - qm
122  if ( dp * dm > 0.0d0) then
123  dq1 = (dp * dm) / (qp - qm)
124  else
125  dq1 = 0.0d0
126  end if
127  end if
128 
129  ! 2nd dimension
130  if (ldim(2)) then
131  qp = q(i,j+1,k)
132  qm = q(i,j-1,k)
133  dp = qp - qijk
134  dm = qijk - qm
135  if ( dp * dm > 0.0d0) then
136  dq2 = (dp * dm) / (qp - qm)
137  else
138  dq2 = 0.0d0
139  end if
140  end if
141 
142  ! 3rd dimension
143  if (ldim(3)) then
144  qp = q(i,j,k+1)
145  qm = q(i,j,k-1)
146  dp = qp - qijk
147  dm = qijk - qm
148  if ( dp * dm .gt. 0.0d0) then
149  dq3 = (dp * dm) / (qp - qm)
150  else
151  dq3 = 0.0d0
152  end if
153  end if
154 
155  ! --- Second, apply a correction (if any) for curvilinear coords. ---
156  ! **Currently disabled**
157  aych1 = 0.0d0
158  aych2 = 0.0d0
159  aych3 = 0.0d0
160 
161  ! Finally, determine the interpolated value.
162  qtilde = qijk + dq1 * ( w(1) - aych1 ) &
163  + dq2 * ( w(2) - aych2 ) &
164  + dq3 * ( w(3) - aych3 )
165 
166  call trace_end (itimer)
167 END FUNCTION interpolator_van_leer_3d
168 
169 !===============================================================================
170 !> Plain, un-limited, trilinear interpolation.
171 !===============================================================================
172 FUNCTION interpolator_trilinear (q,i,j,k,w) result(qtilde)
173  real(kind=KindScalarVar), intent(in):: q(:,:,:)
174  integer, intent(in):: i, j, k
175  real, intent(in):: w(3)
176  real(kind=KindScalarVar):: qtilde, m(3)
177  integer:: i1, j1, k1
178  integer:: itimer=0
179  !-----------------------------------------------------------------------------
180  !call trace_begin('interpolate_mod::trilinear', 2, itimer=itimer)
181  i1 = merge(1,0,size(q,1) > 1)
182  j1 = merge(1,0,size(q,2) > 1)
183  k1 = merge(1,0,size(q,3) > 1)
184  m = 1.0 - w
185  qtilde = m(3) * (m(2) * (m(1) * q(i ,j ,k ) + w(1) * q(i+i1,j ,k )) + &
186  w(2) * (m(1) * q(i ,j+j1,k ) + w(1) * q(i+i1,j+j1,k ))) + &
187  w(3) * (m(2) * (m(1) * q(i ,j ,k+k1) + w(1) * q(i+i1,j ,k+k1)) + &
188  w(2) * (m(1) * q(i ,j+j1,k+k1) + w(1) * q(i+i1,j+j1,k+k1)))
189  !call trace_end (itimer)
190 END FUNCTION interpolator_trilinear
191 
192 !===============================================================================
193 !> Plain, un-limited, trilinear interpolation.
194 !===============================================================================
195 FUNCTION trilinear_log (q,i,j,k,w) result(qtilde)
196  real(kind=KindScalarVar), intent(in):: q(:,:,:)
197  integer, intent(in):: i, j, k
198  real, intent(in):: w(3)
199  real(kind=KindScalarVar):: qtilde, m(3)
200  integer:: i1, j1, k1
201  integer:: itimer=0
202  !-----------------------------------------------------------------------------
203  !call trace_begin('interpolate_mod::trilinear', 2, itimer=itimer)
204  i1 = merge(1,0,size(q,1) > 1)
205  j1 = merge(1,0,size(q,2) > 1)
206  k1 = merge(1,0,size(q,3) > 1)
207  m = 1.0 - w
208  qtilde = m(3)*(m(2)*(m(1)*log(q(i ,j ,k )) + w(1)*log(q(i+i1,j ,k ))) + &
209  w(2)*(m(1)*log(q(i ,j+j1,k )) + w(1)*log(q(i+i1,j+j1,k )))) + &
210  w(3)*(m(2)*(m(1)*log(q(i ,j ,k+k1)) + w(1)*log(q(i+i1,j ,k+k1))) + &
211  w(2)*(m(1)*log(q(i ,j+j1,k+k1)) + w(1)*log(q(i+i1,j+j1,k+k1))))
212  !call trace_end (itimer)
213 END FUNCTION trilinear_log
214 
215 !===============================================================================
216 !> Plain, un-limited, trilinear interpolation.
217 !===============================================================================
218 FUNCTION trilinear_pv (q,d,i,j,k,w) result(qtilde)
219  real(kind=KindScalarVar), intent(in):: d(:,:,:), q(:,:,:)
220  integer, intent(in):: i, j, k
221  real, intent(in):: w(3)
222  real(kind=KindScalarVar):: qtilde, m(3)
223  integer:: i1, j1, k1
224  integer:: itimer=0
225  !-----------------------------------------------------------------------------
226  !call trace_begin('interpolate_mod::trilinear', 2, itimer=itimer)
227  i1 = merge(1,0,size(q,1) > 1)
228  j1 = merge(1,0,size(q,2) > 1)
229  k1 = merge(1,0,size(q,3) > 1)
230  m = 1.0 - w
231  qtilde = m(3)*(m(2)*(m(1)*q(i ,j ,k )/d(i ,j ,k ) + &
232  w(1)*q(i+i1,j ,k )/d(i+i1,j ,k )) + &
233  w(2)*(m(1)*q(i ,j+j1,k )/d(i ,j+j1,k ) + &
234  w(1)*q(i+i1,j+j1,k )/d(i+i1,j+j1,k ))) + &
235  w(3)*(m(2)*(m(1)*q(i ,j ,k+k1)/d(i ,j ,k+k1) + &
236  w(1)*q(i+i1,j ,k+k1)/d(i+i1,j ,k+k1)) + &
237  w(2)*(m(1)*q(i ,j+j1,k+k1)/d(i ,j+j1,k+k1) + &
238  w(1)*q(i+i1,j+j1,k+k1)/d(i+i1,j+j1,k+k1)))
239  !call trace_end (itimer)
240 END FUNCTION trilinear_pv
241 
242 !===============================================================================
243 !> Plain, un-limited, trilinear interpolation.
244 !===============================================================================
245 FUNCTION interpolator_trilinear_3d (q,i,j,k,w) result(qtilde)
246  real(kind=KindScalarVar), intent(in):: q(:,:,:)
247  integer, intent(in):: i, j, k
248  real, intent(in):: w(3)
249  real(kind=KindScalarVar):: qtilde, m(3)
250  integer:: i1, j1, k1
251  integer:: itimer=0
252  !-----------------------------------------------------------------------------
253  !call trace_begin('interpolate_mod::trilinear', 2, itimer=itimer)
254  m = 1.0 - w
255  qtilde = m(3) * (m(2) * (m(1) * q(i ,j ,k ) + w(1) * q(i+1,j ,k )) + &
256  w(2) * (m(1) * q(i ,j+1,k ) + w(1) * q(i+1,j+1,k ))) + &
257  w(3) * (m(2) * (m(1) * q(i ,j ,k+1) + w(1) * q(i+1,j ,k+1)) + &
258  w(2) * (m(1) * q(i ,j+1,k+1) + w(1) * q(i+1,j+1,k+1)))
259  !call trace_end (itimer)
260 END FUNCTION interpolator_trilinear_3d
261 
262 !===============================================================================
263 !> Do the interpolation under the assumption that the variable is
264 !> positive-definite (i.e. `patch%unsigned(iv) = .true.`).
265 !===============================================================================
266 FUNCTION interpolator_unsigned (q,i,j,k,w) result(qtilde)
267  real(kind=KindScalarVar), intent(in):: q(:,:,:)
268  integer, intent(in):: i, j, k
269  real, intent(in):: w(3)
270  real(kind=KindScalarVar):: qtilde
271  real(kind=KindScalarVar), allocatable:: logq(:,:,:)
272  integer:: itimer=0
273  !-----------------------------------------------------------------------------
274  call trace_begin('interpolate_mod::unsigned', 1, itimer=itimer)
275  !
276  allocate(logq(size(q,1),size(q,2),size(q,3)))
277  logq(:,:,:) = log(q(:,:,:))
278  qtilde = exp(interpolator%interpolator(logq,i,j,k,w))
279  !
280  deallocate (logq)
281  call trace_end (itimer)
282 END FUNCTION interpolator_unsigned
283 
284 !===============================================================================
285 !> Plain, un-limited, trilinear interpolation.
286 !===============================================================================
287 FUNCTION four_d (q, i, j, k, l, w) RESULT (qtilde)
288  real(kind=KindScalarVar), intent(in) :: q(:,:,:,:)
289  real(kind=KindScalarVar) :: qtilde
290  integer, intent(in) :: i, j, k, l(2)
291  real, intent(in) :: w(4)
292  !.............................................................................
293  integer:: l1, l2
294  real(kind=KindScalarVar):: m(4)
295  integer:: itimer=0
296  !-----------------------------------------------------------------------------
297  !call trace_begin('interpolate_mod::trilinear', 2, itimer=itimer)
298  l1 = l(1)
299  l2 = l(2)
300  m = 1.0 - w
301  qtilde = &
302  m(4)*(m(3)*(m(2)*(m(1)*q(i ,j ,k ,l1) + w(1)*q(i+1,j ,k ,l1)) + &
303  w(2)*(m(1)*q(i ,j+1,k ,l1) + w(1)*q(i+1,j+1,k ,l1))) + &
304  w(3)*(m(2)*(m(1)*q(i ,j ,k+1,l1) + w(1)*q(i+1,j ,k+1,l1)) + &
305  w(2)*(m(1)*q(i ,j+1,k+1,l1) + w(1)*q(i+1,j+1,k+1,l1)))) + &
306  w(4)*(m(3)*(m(2)*(m(1)*q(i ,j ,k ,l2) + w(1)*q(i+1,j ,k ,l2)) + &
307  w(2)*(m(1)*q(i ,j+1,k ,l2) + w(1)*q(i+1,j+1,k ,l2))) + &
308  w(3)*(m(2)*(m(1)*q(i ,j ,k+1,l2) + w(1)*q(i+1,j ,k+1,l2)) + &
309  w(2)*(m(1)*q(i ,j+1,k+1,l2) + w(1)*q(i+1,j+1,k+1,l2))))
310  !call trace_end (itimer)
311 END FUNCTION four_d
312 
313 !===============================================================================
314 !> Plain, un-limited, trilinear interpolation.
315 !===============================================================================
316 FUNCTION four_d_pv (q, d1, d2, i, j, k, l, w) RESULT (qtilde)
317  real(kind=KindScalarVar), intent(in) :: q(:,:,:,:), d1(:,:,:), d2(:,:,:)
318  real(kind=KindScalarVar) :: qtilde
319  integer, intent(in) :: i, j, k, l(2)
320  real, intent(in) :: w(4)
321  !.............................................................................
322  integer:: l1, l2
323  real(kind=KindScalarVar):: m(4)
324  integer:: itimer=0
325  !-----------------------------------------------------------------------------
326  !call trace_begin('interpolate_mod::trilinear', 2, itimer=itimer)
327  l1 = l(1)
328  l2 = l(2)
329  m = 1.0 - w
330  qtilde = &
331  m(4)*(m(3)*(m(2)*(m(1)*q(i ,j ,k ,l1)/d1(i ,j ,k ) + &
332  w(1)*q(i+1,j ,k ,l1)/d1(i+1,j ,k )) + &
333  w(2)*(m(1)*q(i ,j+1,k ,l1)/d1(i ,j+1,k ) + &
334  w(1)*q(i+1,j+1,k ,l1)/d1(i+1,j+1,k ))) + &
335  w(3)*(m(2)*(m(1)*q(i ,j ,k+1,l1)/d1(i ,j ,k+1) + &
336  w(1)*q(i+1,j ,k+1,l1)/d1(i+1,j ,k+1)) + &
337  w(2)*(m(1)*q(i ,j+1,k+1,l1)/d1(i ,j+1,k+1) + &
338  w(1)*q(i+1,j+1,k+1,l1)/d1(i+1,j+1,k+1)))) + &
339  w(4)*(m(3)*(m(2)*(m(1)*q(i ,j ,k ,l2)/d2(i ,j ,k ) + &
340  w(1)*q(i+1,j ,k ,l2)/d2(i+1,j ,k )) + &
341  w(2)*(m(1)*q(i ,j+1,k ,l2)/d2(i ,j+1,k ) + &
342  w(1)*q(i+1,j+1,k ,l2)/d2(i+1,j+1,k ))) + &
343  w(3)*(m(2)*(m(1)*q(i ,j ,k+1,l2)/d2(i ,j ,k+1) + &
344  w(1)*q(i+1,j ,k+1,l2)/d2(i+1,j ,k+1)) + &
345  w(2)*(m(1)*q(i ,j+1,k+1,l2)/d2(i ,j+1,k+1) + &
346  w(1)*q(i+1,j+1,k+1,l2)/d2(i+1,j+1,k+1))))
347  !call trace_end (itimer)
348 END FUNCTION four_d_pv
349 
350 !===============================================================================
351 !> Plain, un-limited, trilinear interpolation.
352 !===============================================================================
353 FUNCTION four_d_log (q, i, j, k, l, w) RESULT (qtilde)
354  real(kind=KindScalarVar), intent(in) :: q(:,:,:,:)
355  real(kind=KindScalarVar) :: qtilde
356  integer, intent(in) :: i, j, k, l(2)
357  real, intent(in) :: w(4)
358  !.............................................................................
359  integer:: l1, l2
360  real(kind=KindScalarVar):: m(4)
361  integer:: itimer=0
362  !-----------------------------------------------------------------------------
363  !call trace_begin('interpolate_mod::trilinear', 2, itimer=itimer)
364  l1 = l(1)
365  l2 = l(2)
366  m = 1.0 - w
367  qtilde = &
368  m(4)*(m(3)*(m(2)*(m(1)*log(q(i ,j ,k ,l1)) + w(1)*log(q(i+1,j ,k ,l1))) + &
369  w(2)*(m(1)*log(q(i ,j+1,k ,l1)) + w(1)*log(q(i+1,j+1,k ,l1)))) + &
370  w(3)*(m(2)*(m(1)*log(q(i ,j ,k+1,l1)) + w(1)*log(q(i+1,j ,k+1,l1))) + &
371  w(2)*(m(1)*log(q(i ,j+1,k+1,l1)) + w(1)*log(q(i+1,j+1,k+1,l1))))) + &
372  w(4)*(m(3)*(m(2)*(m(1)*log(q(i ,j ,k ,l2)) + w(1)*log(q(i+1,j ,k ,l2))) + &
373  w(2)*(m(1)*log(q(i ,j+1,k ,l2)) + w(1)*log(q(i+1,j+1,k ,l2)))) + &
374  w(3)*(m(2)*(m(1)*log(q(i ,j ,k+1,l2)) + w(1)*log(q(i+1,j ,k+1,l2))) + &
375  w(2)*(m(1)*log(q(i ,j+1,k+1,l2)) + w(1)*log(q(i+1,j+1,k+1,l2)))))
376  !call trace_end (itimer)
377 END FUNCTION four_d_log
378 
379 END MODULE interpolate_mod
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...
Definition: index_mod.f90:7
Definition: io_mod.f90:4