Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module interp_1d_mp ! high accuracy monotonicity preserving algorithms
21 :
22 : use const_def, only: dp
23 :
24 : implicit none
25 : private
26 : public :: m3, m3_on_uniform_grid
27 :
28 : contains
29 :
30 :
31 1315 : subroutine m3(x, nx, f1, which, slope_only, nwork, work1, str, ierr)
32 : use interp_1d_def
33 : use interp_1d_misc
34 : use interp_1d_pm, only: mk_pmlinear, mk_pmquad
35 : integer, intent(in) :: nx ! length of x vector
36 : real(dp), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
37 : real(dp), intent(inout), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
38 : integer, intent(in) :: which
39 : logical, intent(in) :: slope_only
40 : integer, intent(in) :: nwork
41 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
42 : character (len=*) :: str
43 : integer, intent(out) :: ierr
44 :
45 1315 : real(dp), dimension(:), pointer :: h, s_mid, s, d, d_mid, e_mid, hd_mid, &
46 1315 : spL, spR, t, tmax, tmp, tmp1, tmp2
47 : real(dp), parameter :: tiny = 1d-20
48 : integer :: i
49 1315 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
50 1315 : f(1:4,1:nx) => f1(1:4*nx)
51 :
52 1315 : ierr = 0
53 :
54 1315 : if (nx < 2) then
55 : return
56 : end if
57 :
58 1315 : if (nx == 2) then
59 0 : call mk_pmlinear(x, f1, slope_only, nwork, work1, str, ierr)
60 0 : return
61 : end if
62 :
63 1315 : if (nx == 3) then
64 0 : call mk_pmquad(x, f1, slope_only, nwork, work1, str, ierr)
65 0 : return
66 : end if
67 :
68 1315 : if (nwork < mp_work_size) then
69 0 : ierr = -1
70 0 : return
71 : end if
72 :
73 : i = 0
74 1315 : h(1:nx) => work1(i+1:i+nx); i = i+nx
75 1315 : s_mid(1:nx) => work1(i+1:i+nx); i = i+nx
76 1315 : s(1:nx) => work1(i+1:i+nx); i = i+nx
77 1315 : d(1:nx) => work1(i+1:i+nx); i = i+nx
78 1315 : d_mid(1:nx) => work1(i+1:i+nx); i = i+nx
79 1315 : e_mid(1:nx) => work1(i+1:i+nx); i = i+nx
80 1315 : hd_mid(1:nx) => work1(i+1:i+nx); i = i+nx
81 1315 : spL(1:nx) => work1(i+1:i+nx); i = i+nx
82 1315 : spR(1:nx) => work1(i+1:i+nx); i = i+nx
83 1315 : t(1:nx) => work1(i+1:i+nx); i = i+nx
84 1315 : tmax(1:nx) => work1(i+1:i+nx); i = i+nx
85 1315 : tmp(1:nx) => work1(i+1:i+nx); i = i+nx
86 1315 : tmp1(1:nx) => work1(i+1:i+nx); i = i+nx
87 1315 : tmp2(1:nx) => work1(i+1:i+nx); i = i+nx
88 :
89 : ! grid spacing
90 13111597 : do i=1,nx-1
91 13111597 : h(i) = x(i+1) - x(i)
92 : end do
93 13111597 : do i = 1, nx-1
94 13111597 : if (h(i) == 0) then
95 : write(*, '(a,1x,2i5,1x,a)') &
96 0 : 'same interpolation x values at', i, i+1, 'for ' // trim(str)
97 0 : ierr = -1
98 0 : return
99 : end if
100 : end do
101 :
102 : ! divided differences
103 13111597 : do i=1,nx-1
104 13111597 : s_mid(i) = (f(1,i+1) - f(1,i)) / h(i) ! eqn 2.1
105 : end do
106 13110282 : do i=2,nx-1
107 13110282 : d(i) = (s_mid(i) - s_mid(i-1)) / (x(i+1) - x(i-1)) ! eqn 3.1
108 : end do
109 : ! need to extend d to full range. simplest way is just to copy from neighbor
110 1315 : d(1) = d(2)
111 1315 : d(nx) = d(nx-1)
112 :
113 : ! d_mid eqn(3.4) -- modified according to eqn (2.27) of Suresh & Huynh, 1997
114 13111597 : do i=1,nx-1
115 13110282 : tmp1(i) = 4*d(i+1) - d(i)
116 13111597 : tmp2(i) = 4*d(i) - d(i+1)
117 : end do
118 1315 : call minmod4(d_mid(1:nx-1), nx-1, d(2:nx), d(1:nx-1), tmp1(1:nx-1), tmp2(1:nx-1))
119 :
120 13111597 : do i=1,nx-1
121 13111597 : hd_mid(i) = h(i)*d_mid(i)
122 : end do
123 :
124 : ! spL(i) = p'(i-1/2)(xi) = smid(i-1) + h(i-1)*d_mid(i-1) from Theorem 1
125 13111597 : do i=1,nx-1
126 13111597 : spL(i+1) = s_mid(i) + hd_mid(i)
127 : end do
128 :
129 : ! spR(i) = p'(i+1/2)(xi) = smid(i) - h(i)*d_mid(i) from Theorem 1
130 13111597 : do i=1,nx-1
131 13111597 : spR(i) = s_mid(i) - hd_mid(i)
132 : end do
133 :
134 1315 : call minmod(s(2:nx-1), nx-2, s_mid(1:nx-2), s_mid(2:nx-1)) ! eqn (2.8)
135 1315 : call minmod(t(2:nx-1), nx-2, spL(2:nx-1), spR(2:nx-1))
136 :
137 1315 : if (which == average) then
138 :
139 0 : do i=2,nx-1
140 : f(2,i) = sign(1d0, t(i))* &
141 : min((dabs(spL(i)+spR(i)))/2d0, &
142 0 : max(3*dabs(s(i)), 1.5d0*dabs(t(i))))
143 : end do
144 :
145 1315 : else if (which == quartic) then
146 :
147 13108967 : do i=2,nx-2
148 13108967 : e_mid(i) = (d(i+1) - d(i)) / (x(i+2) - x(i-1)) ! eqn 4.1
149 : end do
150 : ! eqn 3.5b for p'(i); eqn 4.20 for quadratic f'(i)
151 13107652 : do i=3,nx-2
152 : f(2,i) = s_mid(i) - &
153 : h(i) * (d(i) + h(i-1)* &
154 : (e_mid(i-1)*(x(i+2)-x(i)) &
155 13107652 : + e_mid(i)*(x(i)-x(i-2)))/(x(i+2)-x(i-2)))
156 : end do
157 : ! finish off ends with average
158 : f(2,2) = sign(1d0, t(2))* &
159 1315 : min((dabs(spL(2)+spR(2)))/2d0, max(3*dabs(s(2)), 1.5d0*dabs(t(2))))
160 : f(2,nx-1) = sign(1d0, t(nx-1))* &
161 : min((dabs(spL(nx-1)+spR(nx-1)))/2d0, &
162 1315 : max(3*dabs(s(nx-1)), 1.5d0*dabs(t(nx-1))))
163 13110282 : do i=2,nx-1
164 13108967 : tmp1(i) = f(2,i)
165 13110282 : tmp2(i) = tmp1(i)
166 : end do
167 1315 : call median(tmp1(2:nx-1), nx-2, tmp2(2:nx-1), spL(2:nx-1), spR(2:nx-1))
168 13110282 : do i=2,nx-1
169 13110282 : f(2,i) = tmp1(i)
170 : end do
171 13110282 : do i=2,nx-1
172 : tmax(i) = sign(1d0, t(i))* &
173 13110282 : max(3*dabs(s(i)), 1.5d0*dabs(t(i)))
174 : end do
175 13110282 : do i=2,nx-1
176 13110282 : tmp1(i) = f(2,i)
177 : end do
178 1315 : call minmod(tmp2(2:nx-1), nx-2, tmp1(2:nx-1), tmax(2:nx-1))
179 13110282 : do i=2,nx-1
180 13110282 : f(2,i) = tmp2(i)
181 : end do
182 :
183 : else !if (which == super_bee) then
184 :
185 0 : do i=2,nx-1
186 : f(2,i) = sign(1d0, t(i))* &
187 : min(max(dabs(spL(i)), dabs(spR(i))), &
188 0 : max(3*dabs(s(i)), 1.5d0*dabs(t(i))))
189 : end do
190 :
191 : end if
192 :
193 : ! slope at i=1
194 : !f(2, 1) = minmod1(spR(1), 3*s_mid(1)) ! eqn (5.2)
195 1315 : f(2,1) = minmod1(s_mid(1), s_mid(2)) ! stabilize the ends
196 :
197 : ! slope at i=nx
198 : !f(2, nx) = minmod1(spL(nx), 3*s_mid(nx-1)) ! eqn (5.2)
199 1315 : f(2,nx) = minmod1(s_mid(nx-2), s_mid(nx-1)) ! stabilize the ends
200 :
201 1315 : if (slope_only) return
202 :
203 : ! 2nd and 3rd derivatives
204 13111597 : do i=1,nx-1
205 13110282 : f(3,i) = (3*s_mid(i) - 2*f(2,i) - f(2,i+1)) / h(i)
206 13111597 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s_mid(i)) / (h(i)*h(i))
207 : end do
208 1315 : f(3,nx) = (3*f(1, nx-1) - 3*f(1, nx) + (f(2, nx-1) + 2*f(2, nx)) * h(nx-1)) / (h(nx-1)*h(nx-1))
209 1315 : f(4,nx) = (-2*f(1, nx-1) + 2*f(1, nx) - (f(2, nx-1) + f(2, nx))*h(nx-1)) / (h(nx-1)*h(nx-1)*h(nx-1))
210 :
211 2630 : end subroutine m3
212 :
213 :
214 0 : subroutine m3_on_uniform_grid(dx, nx, f1, which, slope_only, nwork, work1, str, ierr)
215 1315 : use interp_1d_def
216 : use interp_1d_misc
217 : use interp_1d_pm, only: mk_pmlinear, mk_pmquad
218 : ! make piecewise monotonic cubic interpolant
219 : real(dp), intent(in) :: dx ! the grid spacing
220 : integer, intent(in) :: nx ! length of x vector
221 : real(dp), intent(inout), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
222 : integer, intent(in) :: which
223 : logical, intent(in) :: slope_only
224 : integer, intent(in) :: nwork
225 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
226 : character (len=*) :: str
227 : integer, intent(out) :: ierr
228 :
229 0 : real(dp), dimension(:), pointer :: s_mid, s, d, d_mid, e_mid, spL, spR, t, tmax, &
230 0 : tmp, tmp1, tmp2
231 : real(dp), parameter :: tiny = 1d-20
232 0 : real(dp) :: x(3)
233 : integer :: i
234 0 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
235 0 : f(1:4,1:nx) => f1(1:4*nx)
236 :
237 0 : ierr = 0
238 :
239 0 : if (nx < 2) then
240 : return
241 : end if
242 :
243 0 : if (nx == 2) then
244 0 : x(1) = 0
245 0 : x(2) = dx
246 0 : call mk_pmlinear(x, f1, slope_only, nwork, work1, str, ierr)
247 0 : return
248 : end if
249 :
250 0 : if (nx == 3) then
251 0 : x(1) = 0
252 0 : x(2) = dx
253 0 : x(3) = 2*dx
254 0 : call mk_pmquad(x, f1, slope_only, nwork, work1, str, ierr)
255 0 : return
256 : end if
257 :
258 0 : if (dx == 0) then
259 0 : ierr = -1
260 0 : return
261 : end if
262 :
263 0 : if (nwork < mp_work_size) then
264 0 : ierr = -1
265 0 : return
266 : end if
267 :
268 0 : i = 0
269 0 : s_mid(1:nx) => work1(i+1:i+nx); i = i+nx
270 0 : s(1:nx) => work1(i+1:i+nx); i = i+nx
271 0 : d(1:nx) => work1(i+1:i+nx); i = i+nx
272 0 : d_mid(1:nx) => work1(i+1:i+nx); i = i+nx
273 0 : e_mid(1:nx) => work1(i+1:i+nx); i = i+nx
274 0 : spL(1:nx) => work1(i+1:i+nx); i = i+nx
275 0 : spR(1:nx) => work1(i+1:i+nx); i = i+nx
276 0 : t(1:nx) => work1(i+1:i+nx); i = i+nx
277 0 : tmax(1:nx) => work1(i+1:i+nx); i = i+nx
278 0 : tmp(1:nx) => work1(i+1:i+nx); i = i+nx
279 0 : tmp1(1:nx) => work1(i+1:i+nx); i = i+nx
280 0 : tmp2(1:nx) => work1(i+1:i+nx); i = i+nx
281 :
282 : ! divided differences
283 0 : do i=1,nx-1
284 0 : s_mid(i) = (f(1,i+1) - f(1,i)) / dx ! eqn 2.1
285 : end do
286 0 : do i=2,nx-1
287 0 : d(i) = (s_mid(i) - s_mid(i-1)) / (2*dx) ! eqn 3.1
288 : end do
289 : ! need to extend d to full range. simplest way is just to copy from neighbor
290 0 : d(1) = d(2)
291 0 : d(nx) = d(nx-1)
292 :
293 : ! d_mid eqn(3.4) -- modified according to eqn (2.27) of Suresh & Huynh, 1997
294 0 : do i=1,nx-1
295 0 : tmp1(i) = 4*d(i+1) - d(i)
296 0 : tmp2(i) = 4*d(i) - d(i+1)
297 : end do
298 0 : call minmod4(d_mid(1:nx-1), nx-1, d(2:nx), d(1:nx-1), tmp1(1:nx-1), tmp2(1:nx-1))
299 :
300 : ! spL(i) = p'(i-1/2)(xi) = smid(i-1) + h(i-1)*d_mid(i-1) from Theorem 1
301 0 : do i=1,nx-1
302 0 : spL(i+1) = s_mid(i) + dx*d_mid(i)
303 : end do
304 :
305 : ! spR(i) = p'(i+1/2)(xi) = smid(i) - h(i)*d_mid(i) from Theorem 1
306 0 : do i=1,nx-1
307 0 : spR(i) = s_mid(i) - dx*d_mid(i)
308 : end do
309 :
310 0 : call minmod(s(2:nx-1), nx-2, s_mid(1:nx-2), s_mid(2:nx-1)) ! eqn (2.8)
311 0 : call minmod(t(2:nx-1), nx-2, spL(2:nx-1), spR(2:nx-1))
312 :
313 0 : if (which == average) then
314 :
315 0 : do i=2,nx-1
316 : f(2,i) = sign(1d0, t(i))* &
317 : min((dabs(spL(i)+spR(i)))/2d0, &
318 0 : max(3*dabs(s(i)), 1.5d0*dabs(t(i))))
319 : end do
320 :
321 0 : else if (which == quartic) then
322 :
323 0 : do i=2,nx-2
324 0 : e_mid(i) = (d(i+1) - d(i)) / (3*dx) ! eqn 4.1
325 : end do
326 : ! eqn 3.5b for p'(i); eqn 4.20 for quadratic f'(i)
327 0 : do i=3,nx-2
328 : f(2,i) = s_mid(i) - &
329 0 : dx * (d(i) + dx*(e_mid(i-1)*(2*dx) + e_mid(i)*(2*dx))/(4*dx))
330 : end do
331 : ! finish off ends with average
332 : f(2,2) = sign(1d0, t(2))* &
333 0 : min((dabs(spL(2)+spR(2)))/2d0, max(3*dabs(s(2)), 1.5d0*dabs(t(2))))
334 : f(2,nx-1) = sign(1d0, t(nx-1))* &
335 : min((dabs(spL(nx-1)+spR(nx-1)))/2d0, &
336 0 : max(3*dabs(s(nx-1)), 1.5d0*dabs(t(nx-1))))
337 0 : do i=2,nx-1
338 0 : tmp1(i) = f(2,i)
339 0 : tmp2(i) = tmp1(i)
340 : end do
341 0 : call median(tmp1(2:nx-1), nx-2, tmp2(2:nx-1), spL(2:nx-1), spR(2:nx-1))
342 0 : do i=2,nx-1
343 0 : f(2,i) = tmp1(i)
344 : end do
345 0 : do i=2,nx-1
346 : tmax(i) = sign(1d0, t(i))* &
347 0 : max(3*dabs(s(i)), 1.5d0*dabs(t(i)))
348 : end do
349 0 : do i=2,nx-1
350 0 : tmp1(i) = f(2,i)
351 : end do
352 0 : call minmod(tmp2(2:nx-1), nx-2, tmp1(2:nx-1), tmax(2:nx-1))
353 0 : do i=2,nx-1
354 0 : f(2,i) = tmp2(i)
355 : end do
356 :
357 : else !if (which == super_bee) then
358 :
359 0 : do i=2,nx-1
360 : f(2,i) = sign(1d0, t(i))* &
361 : min(max(dabs(spL(i)), dabs(spR(i))), &
362 0 : max(3*dabs(s(i)), 1.5d0*dabs(t(i))))
363 : end do
364 :
365 : end if
366 :
367 : ! slope at i=1
368 : !f(2, 1) = minmod1(spR(1), 3*s_mid(1)) ! eqn (5.2)
369 0 : f(2,1) = minmod1(s_mid(1), s_mid(2)) ! stabilize the ends
370 :
371 : ! slope at i=nx
372 : !f(2, nx) = minmod1(spL(nx), 3*s_mid(nx-1)) ! eqn (5.2)
373 0 : f(2, nx) = minmod1(s_mid(nx-2), s_mid(nx-1)) ! stabilize the ends
374 :
375 0 : if (slope_only) return
376 :
377 : ! 2nd and 3rd derivatives
378 0 : do i=1,nx-1
379 0 : f(3,i) = (3*s_mid(i) - 2*f(2,i) - f(2,i+1)) / dx
380 0 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s_mid(i)) / (dx*dx)
381 : end do
382 0 : f(3, nx) = (3*f(1, nx-1) - 3*f(1, nx) + (f(2, nx-1) + 2*f(2, nx)) * dx) / (dx*dx)
383 0 : f(4, nx) = (-2*f(1, nx-1) + 2*f(1, nx) - (f(2, nx-1) + f(2, nx))*dx) / (dx*dx*dx)
384 :
385 0 : end subroutine m3_on_uniform_grid
386 :
387 :
388 : end module interp_1d_mp
|
