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