[FLASH-USERS] Fwd: Error in ut_hunt ?
Klaus Weide
klaus at flash.uchicago.edu
Mon Jan 11 16:03:24 EST 2021
On Fri, 8 Jan 2021, RIGON Gabriel wrote:
> I was working on some simulations and I saw something that might be an error
> in the ut_hunt.F90 subroutine.
> Directory: ./source/flashUtilities/interpolation/oneDim/ut_hunt.F90
>
> Take a monotonous array (xx) an a value (x) as argument return the highest
> index (i) such as xx <= x < xx(i)
More precisely (dealing here only with the strictly monotonically
increasing case), ut_hunt attempts to implement:
Take a monotonous array xx(1:n) and a value (x) as argument, return the
highest index (i) such that
x < xx(i+1) if x < xx(1),
xx(i) <= x < xx(i+1) if x < xx(n),
xx(n) otherwise.
Therefore I agree that your suggested change appears more appropriate.
Possibly the case 'x.eq.xx(1)' should be handled explicitly in addition,
as in the three lines commented out with '!!$'.
Ususally this doesn't matter in practice; the subroutine is often used
like this
call ut_hunt(xx,n,x,i)
i = min(max(i,0),n-1)
dx = x - xx(i)
! Apply linear or higher polynomial interpolation:
Approximate slope_at_i based on xx(i) and neighboring point(s)
x_interpolated = x(i) + slope_at_i * dx [ + higher_terms(dx) ]
(Interpolation becomes extrapolation for x outside of the interval
[x(1),x(n)].)
So it doesn't matter what index is returned by ut_hunt when x falls
exactly on the lower or upper boundary, since it is clamped anyway.
Klaus
More information about the flash-users
mailing list