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_sg ! high accuracy monotonicity preserving algorithms
21 :
22 : implicit none
23 : private
24 : public :: m3_sg, m3_on_uniform_grid_sg
25 :
26 : contains
27 :
28 :
29 4 : subroutine m3_sg(x,nx, f1, which, slope_only, nwork, work1, str, ierr)
30 : ! make piecewise monotonic cubic interpolant
31 : use interp_1d_def
32 : use interp_1d_misc_sg
33 : use interp_1d_pm_sg, only: mk_pmlinear_sg, mk_pmquad_sg
34 : integer, intent(in) :: nx ! length of x vector
35 : real, intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
36 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
37 : integer, intent(in) :: which
38 : logical, intent(in) :: slope_only
39 : integer, intent(in) :: nwork
40 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
41 : character (len=*) :: str
42 : integer, intent(out) :: ierr
43 :
44 4 : real, dimension(:), pointer :: h, s_mid, s, d, d_mid, e_mid, hd_mid, &
45 4 : spL, spR, t, tmax, tmp, tmp1, tmp2
46 : real, parameter :: tiny = 1e-20
47 : integer :: i
48 4 : real, pointer :: f(:,:) ! (4,nx) ! data & interpolation coefficients
49 4 : f(1:4,1:nx) => f1(1:4*nx)
50 :
51 4 : ierr = 0
52 :
53 4 : if (nx < 2) then
54 0 : return
55 : end if
56 :
57 4 : if (nx == 2) then
58 0 : call mk_pmlinear_sg(x, f1, slope_only, nwork, work1, str, ierr)
59 0 : return
60 : end if
61 :
62 4 : if (nx == 3) then
63 0 : call mk_pmquad_sg(x, f1, slope_only, nwork, work1, str, ierr)
64 0 : return
65 : end if
66 :
67 4 : if (nwork < mp_work_size) then
68 0 : ierr = -1
69 0 : return
70 : end if
71 :
72 : i = 0
73 4 : h(1:nx) => work1(i+1:i+nx); i = i+nx
74 4 : s_mid(1:nx) => work1(i+1:i+nx); i = i+nx
75 4 : s(1:nx) => work1(i+1:i+nx); i = i+nx
76 4 : d(1:nx) => work1(i+1:i+nx); i = i+nx
77 4 : d_mid(1:nx) => work1(i+1:i+nx); i = i+nx
78 4 : e_mid(1:nx) => work1(i+1:i+nx); i = i+nx
79 4 : hd_mid(1:nx) => work1(i+1:i+nx); i = i+nx
80 4 : spL(1:nx) => work1(i+1:i+nx); i = i+nx
81 4 : spR(1:nx) => work1(i+1:i+nx); i = i+nx
82 4 : t(1:nx) => work1(i+1:i+nx); i = i+nx
83 4 : tmax(1:nx) => work1(i+1:i+nx); i = i+nx
84 4 : tmp(1:nx) => work1(i+1:i+nx); i = i+nx
85 4 : tmp1(1:nx) => work1(i+1:i+nx); i = i+nx
86 4 : tmp2(1:nx) => work1(i+1:i+nx); i = i+nx
87 :
88 : ! grid spacing
89 286 : do i=1,nx-1
90 286 : h(i) = x(i+1) - x(i)
91 : end do
92 286 : do i = 1,nx-1
93 286 : if (h(i) == 0) then
94 : write(*, '(a,1x,2i5,1x,a)') &
95 0 : 'same interpolation x values at', i, i+1, 'for ' // trim(str)
96 0 : ierr = -1
97 0 : return
98 : end if
99 : end do
100 :
101 : ! divided differences
102 286 : do i=1,nx-1
103 286 : s_mid(i) = (f(1,i+1) - f(1,i)) / h(i) ! eqn 2.1
104 : end do
105 282 : do i=2,nx-1
106 282 : d(i) = (s_mid(i) - s_mid(i-1)) / (x(i+1) - x(i-1)) ! eqn 3.1
107 : end do
108 : ! need to extend d to full range. simplest way is just to copy from neighbor
109 4 : d(1) = d(2)
110 4 : d(nx) = d(nx-1)
111 :
112 : ! d_mid eqn(3.4) -- modified according to eqn (2.27) of Suresh & Huynh, 1997
113 286 : do i=1,nx-1
114 282 : tmp1(i) = 4*d(i+1) - d(i)
115 286 : tmp2(i) = 4*d(i) - d(i+1)
116 : end do
117 4 : call minmod4_sg(d_mid(1:nx-1),nx-1, d(2:nx), d(1:nx-1), tmp1(1:nx-1), tmp2(1:nx-1))
118 :
119 286 : do i=1,nx-1
120 286 : hd_mid(i) = h(i)*d_mid(i)
121 : end do
122 : ! spL(i) = p'(i-1/2)(xi) = smid(i-1) + h(i-1)*d_mid(i-1) from Theorem 1
123 286 : do i=1,nx-1
124 286 : spL(i+1) = s_mid(i) + hd_mid(i)
125 : end do
126 : ! spR(i) = p'(i+1/2)(xi) = smid(i) - h(i)*d_mid(i) from Theorem 1
127 286 : do i=1,nx-1
128 286 : spR(i) = s_mid(i) - hd_mid(i)
129 : end do
130 :
131 4 : call minmod_sg(s(2:nx-1),nx-2, s_mid(1:nx-2), s_mid(2:nx-1)) ! eqn (2.8)
132 4 : call minmod_sg(t(2:nx-1),nx-2, spL(2:nx-1), spR(2:nx-1))
133 :
134 4 : if (which == average) then
135 :
136 0 : do i=2,nx-1
137 : f(2,i) = sign(1.0, t(i))* &
138 : min((abs(spL(i)+spR(i)))/2, &
139 0 : max(3*abs(s(i)), 1.5*abs(t(i))))
140 : end do
141 :
142 4 : else if (which == quartic) then
143 :
144 278 : do i=2,nx-2
145 278 : e_mid(i) = (d(i+1) - d(i)) / (x(i+2) - x(i-1)) ! eqn 4.1
146 : end do
147 : ! eqn 3.5b for p'(i); eqn 4.20 for quadratic f'(i)
148 274 : do i=3,nx-2
149 : f(2,i) = s_mid(i) - &
150 : h(i) * (d(i) + h(i-1)* &
151 : (e_mid(i-1)*(x(i+2)-x(i)) &
152 274 : + e_mid(i)*(x(i)-x(i-2)))/(x(i+2)-x(i-2)))
153 : end do
154 : ! finish off ends with average
155 : f(2,2) = sign(1.0, t(2))* &
156 4 : min((abs(spL(2)+spR(2)))/2, max(3*abs(s(2)), 1.5*abs(t(2))))
157 : f(2,nx-1) = sign(1.0, t(nx-1))* &
158 : min((abs(spL(nx-1)+spR(nx-1)))/2, &
159 4 : max(3*abs(s(nx-1)), 1.5*abs(t(nx-1))))
160 282 : do i=2,nx-1
161 278 : tmp1(i) = f(2,i)
162 282 : tmp2(i) = tmp1(i)
163 : end do
164 4 : call median_sg(tmp1(2:nx-1),nx-2, tmp2(2:nx-1), spL(2:nx-1), spR(2:nx-1))
165 282 : do i=2,nx-1
166 282 : f(2,i) = tmp1(i)
167 : end do
168 282 : do i=2,nx-1
169 : tmax(i) = sign(1.0, t(i))* &
170 282 : max(3*abs(s(i)), 1.5*abs(t(i)))
171 : end do
172 282 : do i=2,nx-1
173 282 : tmp1(i) = f(2,i)
174 : end do
175 4 : call minmod_sg(tmp2(2:nx-1),nx-2, tmp1(2:nx-1), tmax(2:nx-1))
176 282 : do i=2,nx-1
177 282 : f(2,i) = tmp2(i)
178 : end do
179 :
180 : else !if (which == super_bee) then
181 :
182 0 : do i=2,nx-1
183 : f(2,i) = sign(1.0, t(i))* &
184 : min(max(abs(spL(i)), abs(spR(i))), &
185 0 : max(3*abs(s(i)), 1.5*abs(t(i))))
186 : end do
187 :
188 : end if
189 :
190 : ! slope at i=1
191 : !f(2, 1) = minmod1_sg(spR(1), 3*s_mid(1)) ! eqn (5.2)
192 4 : f(2,1) = minmod1_sg(s_mid(1), s_mid(2)) ! stabilize the ends
193 :
194 : ! slope at i=nx
195 : !f(2,nx) = minmod1_sg(spL(nx), 3*s_mid(nx-1)) ! eqn (5.2)
196 4 : f(2,nx) = minmod1_sg(s_mid(nx-2), s_mid(nx-1)) ! stabilize the ends
197 :
198 4 : if (slope_only) return
199 :
200 : ! 2nd and 3rd derivatives
201 286 : do i=1,nx-1
202 282 : f(3,i) = (3*s_mid(i) - 2*f(2,i) - f(2,i+1)) / h(i)
203 286 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s_mid(i)) / (h(i)*h(i))
204 : end do
205 4 : 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))
206 4 : 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))
207 :
208 4 : end subroutine m3_sg
209 :
210 :
211 0 : subroutine m3_on_uniform_grid_sg(dx,nx, f1, which, slope_only, nwork, work1, str, ierr)
212 : use interp_1d_def
213 : use interp_1d_misc_sg
214 : use interp_1d_pm_sg, only: mk_pmlinear_sg, mk_pmquad_sg
215 : ! make piecewise monotonic cubic interpolant
216 : real, intent(in) :: dx ! the grid spacing
217 : integer, intent(in) :: nx ! length of x vector
218 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
219 : integer, intent(in) :: which
220 : logical, intent(in) :: slope_only
221 : integer, intent(in) :: nwork
222 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
223 : character (len=*) :: str
224 : integer, intent(out) :: ierr
225 :
226 0 : real, dimension(:), pointer :: s_mid, s, d, d_mid, e_mid, spL, spR, t, tmax, &
227 0 : tmp, tmp1, tmp2
228 : real, parameter :: tiny = 1e-20
229 0 : real :: x(3)
230 : integer :: i
231 0 : real, pointer :: f(:,:) ! (4,nx) ! data & interpolation coefficients
232 0 : f(1:4,1:nx) => f1(1:4*nx)
233 :
234 0 : ierr = 0
235 :
236 0 : if (nx < 2) then
237 0 : return
238 : end if
239 :
240 0 : if (nx == 2) then
241 0 : x(1) = 0
242 0 : x(2) = dx
243 0 : call mk_pmlinear_sg(x, f1, slope_only, nwork, work1, str, ierr)
244 0 : return
245 : end if
246 :
247 0 : if (nx == 3) then
248 0 : x(1) = 0
249 0 : x(2) = dx
250 0 : x(3) = 2*dx
251 0 : call mk_pmquad_sg(x, f1, slope_only, nwork, work1, str, ierr)
252 0 : return
253 : end if
254 :
255 0 : if (dx == 0) then
256 0 : ierr = -1
257 0 : return
258 : end if
259 :
260 0 : if (nwork < mp_work_size) then
261 0 : ierr = -1
262 0 : return
263 : end if
264 :
265 0 : i = 0
266 0 : s_mid(1:nx) => work1(i+1:i+nx); i = i+nx
267 0 : s(1:nx) => work1(i+1:i+nx); i = i+nx
268 0 : d(1:nx) => work1(i+1:i+nx); i = i+nx
269 0 : d_mid(1:nx) => work1(i+1:i+nx); i = i+nx
270 0 : e_mid(1:nx) => work1(i+1:i+nx); i = i+nx
271 0 : spL(1:nx) => work1(i+1:i+nx); i = i+nx
272 0 : spR(1:nx) => work1(i+1:i+nx); i = i+nx
273 0 : t(1:nx) => work1(i+1:i+nx); i = i+nx
274 0 : tmax(1:nx) => work1(i+1:i+nx); i = i+nx
275 0 : tmp(1:nx) => work1(i+1:i+nx); i = i+nx
276 0 : tmp1(1:nx) => work1(i+1:i+nx); i = i+nx
277 0 : tmp2(1:nx) => work1(i+1:i+nx); i = i+nx
278 :
279 : ! divided differences
280 0 : do i=1,nx-1
281 0 : s_mid(i) = (f(1,i+1) - f(1,i)) / dx ! eqn 2.1
282 : end do
283 0 : do i=2,nx-1
284 0 : d(i) = (s_mid(i) - s_mid(i-1)) / (2*dx) ! eqn 3.1
285 : end do
286 : ! need to extend d to full range. simplest way is just to copy from neighbor
287 0 : d(1) = d(2)
288 0 : d(nx) = d(nx-1)
289 :
290 : ! d_mid eqn(3.4) -- modified according to eqn (2.27) of Suresh & Huynh, 1997
291 0 : do i=1,nx-1
292 0 : tmp1(i) = 4*d(i+1) - d(i)
293 0 : tmp2(i) = 4*d(i) - d(i+1)
294 : end do
295 0 : call minmod4_sg(d_mid(1:nx-1),nx-1, d(2:nx), d(1:nx-1), tmp1(1:nx-1), tmp2(1:nx-1))
296 :
297 : ! spL(i) = p'(i-1/2)(xi) = smid(i-1) + h(i-1)*d_mid(i-1) from Theorem 1
298 0 : do i=1,nx-1
299 0 : spL(i+1) = s_mid(i) + dx*d_mid(i)
300 : end do
301 : ! spR(i) = p'(i+1/2)(xi) = smid(i) - h(i)*d_mid(i) from Theorem 1
302 0 : do i=1,nx-1
303 0 : spR(i) = s_mid(i) - dx*d_mid(i)
304 : end do
305 :
306 0 : call minmod_sg(s(2:nx-1),nx-2, s_mid(1:nx-2), s_mid(2:nx-1)) ! eqn (2.8)
307 0 : call minmod_sg(t(2:nx-1),nx-2, spL(2:nx-1), spR(2:nx-1))
308 :
309 0 : if (which == average) then
310 :
311 0 : do i=2,nx-1
312 : f(2,i) = sign(1.0, t(i))* &
313 : min((abs(spL(i)+spR(i)))/2, &
314 0 : max(3*abs(s(i)), 1.5*abs(t(i))))
315 : end do
316 :
317 0 : else if (which == quartic) then
318 :
319 0 : do i=2,nx-2
320 0 : e_mid(i) = (d(i+1) - d(i)) / (3*dx) ! eqn 4.1
321 : end do
322 : ! eqn 3.5b for p'(i); eqn 4.20 for quadratic f'(i)
323 0 : do i=3,nx-2
324 : f(2,i) = s_mid(i) - &
325 0 : dx * (d(i) + dx*(e_mid(i-1)*(2*dx) + e_mid(i)*(2*dx))/(4*dx))
326 : end do
327 : ! finish off ends with average
328 : f(2,2) = sign(1.0, t(2))* &
329 0 : min((abs(spL(2)+spR(2)))/2, max(3*abs(s(2)), 1.5*abs(t(2))))
330 : f(2,nx-1) = sign(1.0, t(nx-1))* &
331 : min((abs(spL(nx-1)+spR(nx-1)))/2, &
332 0 : max(3*abs(s(nx-1)), 1.5*abs(t(nx-1))))
333 0 : do i=2,nx-1
334 0 : tmp1(i) = f(2,i)
335 0 : tmp2(i) = tmp1(i)
336 : end do
337 0 : call median_sg(tmp1(2:nx-1),nx-2, tmp2(2:nx-1), spL(2:nx-1), spR(2:nx-1))
338 0 : do i=2,nx-1
339 0 : f(2,i) = tmp1(i)
340 : end do
341 0 : do i=2,nx-1
342 : tmax(i) = sign(1.0, t(i))* &
343 0 : max(3*abs(s(i)), 1.5*abs(t(i)))
344 : end do
345 0 : do i=2,nx-1
346 0 : tmp1(i) = f(2,i)
347 : end do
348 0 : call minmod_sg(tmp2(2:nx-1),nx-2, tmp1(2:nx-1), tmax(2:nx-1))
349 0 : do i=2,nx-1
350 0 : f(2,i) = tmp2(i)
351 : end do
352 :
353 : else !if (which == super_bee) then
354 :
355 0 : do i=2,nx-1
356 : f(2,i) = sign(1.0, t(i))* &
357 : min(max(abs(spL(i)), abs(spR(i))), &
358 0 : max(3*abs(s(i)), 1.5*abs(t(i))))
359 : end do
360 :
361 : end if
362 :
363 : ! slope at i=1
364 : !f(2, 1) = minmod1_sg(spR(1), 3*s_mid(1)) ! eqn (5.2)
365 0 : f(2,1) = minmod1_sg(s_mid(1), s_mid(2)) ! stabilize the ends
366 :
367 : ! slope at i=nx
368 : !f(2,nx) = minmod1_sg(spL(nx), 3*s_mid(nx-1)) ! eqn (5.2)
369 0 : f(2,nx) = minmod1_sg(s_mid(nx-2), s_mid(nx-1)) ! stabilize the ends
370 :
371 0 : if (slope_only) return
372 :
373 : ! 2nd and 3rd derivatives
374 0 : do i=1,nx-1
375 0 : f(3,i) = (3*s_mid(i) - 2*f(2,i) - f(2,i+1)) / dx
376 0 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s_mid(i)) / (dx*dx)
377 : end do
378 0 : f(3,nx) = (3*f(1,nx-1) - 3*f(1,nx) + (f(2,nx-1) + 2*f(2,nx)) * dx) / (dx*dx)
379 0 : f(4,nx) = (-2*f(1,nx-1) + 2*f(1,nx) - (f(2,nx-1) + f(2,nx))*dx) / (dx*dx*dx)
380 :
381 0 : end subroutine m3_on_uniform_grid_sg
382 :
383 :
384 : end module interp_1d_mp_sg
|
