Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 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 kap_eval_support
21 :
22 : use const_def, only: dp, one_sixth
23 : use math_lib
24 :
25 : implicit none
26 :
27 : private
28 : public :: do_kap_interpolations
29 : public :: locate_logT
30 : public :: locate_logR
31 :
32 : contains
33 :
34 520556 : subroutine locate_log( &
35 : rq, num_logs, log_min_in, log_max_in, ili_logs, logs, log_find, i, log0, log1, ierr)
36 : use kap_def
37 : use utils_lib, only: is_bad
38 : use num_lib, only: binary_search
39 : type (Kap_General_Info), pointer :: rq
40 : integer, intent(in) :: num_logs, ili_logs
41 : real(dp), intent(in) :: log_min_in, log_max_in
42 : real(dp), intent(in), pointer :: logs(:) ! (num_logs)
43 : real(dp), intent(inout) :: log_find ! can change log_find if clipping to table boundaries
44 : integer, intent(out) :: i ! index in logs s.t. logs(i) <= log_find < logs(i+1)
45 : ! one exception: if log_find == log_max then will get i = num_logs-1
46 : real(dp), intent(out) :: log0, log1
47 : integer, intent(out) :: ierr
48 520556 : real(dp) :: dlog, log_min, log_max
49 : integer :: j
50 520556 : ierr = 0
51 520556 : log_min = max(log_min_in, logs(1))
52 520556 : log_max = min(log_max_in, logs(num_logs))
53 520556 : if (num_logs == 1) then
54 0 : i = 1
55 0 : log_find = log_min
56 0 : log0 = log_find
57 0 : log1 = log_find
58 0 : return
59 : end if
60 520556 : if (log_find < log_min .or. log_find > log_max) then
61 0 : if (.not. clip_to_kap_table_boundaries) then
62 0 : ierr = -1
63 0 : return
64 : end if
65 0 : if (log_find < log_min) then
66 0 : i = 1
67 0 : log_find = log_min
68 : else
69 0 : i = num_logs-1
70 0 : log_find = log_max
71 : end if
72 520556 : else if (abs(log_find-log_max) < 1d-7) then
73 0 : i = num_logs-1
74 0 : log_find = log_max
75 520556 : else if (ili_logs == 1) then ! logs equally spaced
76 260278 : dlog = (log_max-log_min)/(num_logs-1)
77 260278 : i = int((log_find-log_min) / dlog) + 1
78 : ! might not be exactly evenly spaced, so minor fixup if necessary
79 260278 : if (logs(i) > log_find .and. i > 1) then
80 0 : i = i-1
81 260278 : else if (log_find >= logs(i+1) .and. i+1 < num_logs) then
82 0 : i = i+1
83 : end if
84 : else
85 260278 : i = binary_search(num_logs, logs, 0, log_find)
86 260278 : if (i >= num_logs) then
87 0 : ierr = -1
88 0 : return
89 : !$OMP critical (kap_eval_crit1)
90 : write(*,*) 'i', i
91 : write(*,*) 'num_logs', num_logs
92 : call mesa_error(__FILE__,__LINE__,'locate_log')
93 : !$OMP end critical (kap_eval_crit1)
94 : end if
95 : end if
96 :
97 520556 : if (i < 1 .or. i >= num_logs) then
98 0 : ierr = -1
99 0 : return
100 : write(*,*) 'i', i
101 : write(*,*) 'num_logs', num_logs
102 : call mesa_error(__FILE__,__LINE__,'locate_log')
103 : end if
104 :
105 520556 : if (logs(i) > log_find .or. log_find > logs(i+1)) then
106 0 : ierr = -1
107 0 : !$OMP critical (kap_eval_crit2)
108 0 : write(*,*) 'dlog', (log_max-log_min)/(num_logs-1)
109 0 : write(*,*) 'log_max', log_max
110 0 : write(*,*) 'log_min', log_min
111 0 : write(*,*) 'log_max_in', log_max_in
112 0 : write(*,*) 'log_min_in', log_min_in
113 0 : write(*,*) 'logs(i)', logs(i)
114 0 : write(*,*) 'log_find', log_find
115 0 : write(*,*) 'logs(i+1)', logs(i+1)
116 0 : write(*,*) 'i', i
117 0 : write(*,*) 'num_logs', num_logs
118 0 : write(*,*) 'ili_logs', ili_logs
119 : !$OMP end critical (kap_eval_crit2)
120 0 : return
121 : call mesa_error(__FILE__,__LINE__,'error in locate_log')
122 : end if
123 520556 : log0 = logs(i)
124 520556 : log1 = logs(i+1)
125 520556 : if (is_bad(log0) .or. is_bad(log1)) then
126 0 : ierr = -1
127 0 : return
128 : !$OMP critical (kap_eval_crit3)
129 : write(*,*) 'logs(i)', logs(i)
130 : write(*,*) 'log_find', log_find
131 : write(*,*) 'logs(i+1)', logs(i+1)
132 : write(*,*) 'i', i
133 : write(*,*) 'num_logs', num_logs
134 : call mesa_error(__FILE__,__LINE__,'error in locate_log')
135 : !$OMP end critical (kap_eval_crit3)
136 : end if
137 520556 : end subroutine locate_log
138 :
139 :
140 260278 : subroutine locate_logT( &
141 : rq, num_logTs, logT_min, logT_max, ili_logTs, logTs, logT, iT, logT0, logT1, ierr)
142 520556 : use kap_def
143 : type (Kap_General_Info), pointer :: rq
144 : integer, intent(in) :: num_logTs, ili_logTs
145 : real(dp), intent(in) :: logT_min, logT_max
146 : real(dp), intent(in), pointer :: logTs(:) ! (num_logTs)
147 : real(dp), intent(inout) :: logT ! can change logT if clipping to table boundaries
148 : integer, intent(out) :: iT ! index in logTs s.t. logTs(i) <= logT < logTs(i+1)
149 : real(dp), intent(out) :: logT0, logT1
150 : integer, intent(out) :: ierr
151 : call locate_log( &
152 260278 : rq, num_logTs, logT_min, logT_max, ili_logTs, logTs, logT, iT, logT0, logT1, ierr)
153 260278 : end subroutine locate_logT
154 :
155 :
156 260278 : subroutine locate_logR( &
157 : rq, num_logRs, logR_min, logR_max, ili_logRs, logRs, logR, iR, logR0, logR1, ierr)
158 : use kap_def
159 : type (Kap_General_Info), pointer :: rq
160 : integer, intent(in) :: num_logRs, ili_logRs
161 : real(dp), intent(in) :: logR_min, logR_max
162 : real(dp), intent(in), pointer :: logRs(:) ! (num_logRs)
163 : real(dp), intent(inout) :: logR ! can change logR if clipping to table boundaries
164 : integer, intent(out) :: iR ! index in logRs s.t. logRs(i) <= logR < logRs(i+1)
165 : real(dp), intent(out) :: logR0, logR1
166 : integer, intent(out) :: ierr
167 : call locate_log( &
168 260278 : rq, num_logRs, logR_min, logR_max, ili_logRs, logRs, logR, iR, logR0, logR1, ierr)
169 260278 : end subroutine locate_logR
170 :
171 :
172 432710 : subroutine do_kap_interpolations( &
173 : fin1, nx, ny, i, j, x0, xget, x1, y0, yget, y1, fval, df_dx, df_dy)
174 : ! derived from routines in the PSPLINE package written by Doug McCune
175 :
176 : real(dp), dimension(:), pointer :: fin1 ! the spline data array, dimensions (4, nx, ny)
177 : integer, intent(in) :: nx, ny, i, j ! target cell in the spline data
178 : real(dp), intent(in) :: x0, xget, x1 ! x0 <= xget <= x1; x0 = xs(i), x1 = xs(i+1)
179 : real(dp), intent(in) :: y0, yget, y1 ! y0 <= yget <= y1; y0 = ys(j), y1 = ys(j+1)
180 : real(dp), intent(out) :: fval, df_dx, df_dy
181 :
182 : real(dp), parameter :: z36th = 1d0 / 36d0
183 :
184 432710 : real(dp), pointer :: fin(:,:,:)
185 :
186 432710 : real(dp) :: xp, xpi, xp2, xpi2, cx, cxi, hx2, cxd, cxdi, hx, hxi
187 432710 : real(dp) :: yp, ypi, yp2, ypi2, cy, cyi, hy2, cyd, cydi, hy, hyi
188 :
189 : logical, parameter :: dbg = .false.
190 :
191 : include 'formats'
192 :
193 432710 : fin(1:4,1:nx,1:ny) => fin1(1:4*nx*ny)
194 :
195 432710 : hx=x1-x0
196 432710 : hxi=1d0/hx
197 432710 : hx2=hx*hx
198 :
199 432710 : xp=(xget-x0)*hxi
200 432710 : xpi=1d0-xp
201 432710 : xp2=xp*xp
202 432710 : xpi2=xpi*xpi
203 :
204 432710 : cx=xp*(xp2-1d0)
205 432710 : cxi=xpi*(xpi2-1d0)
206 432710 : cxd=3d0*xp2-1d0
207 432710 : cxdi=-3d0*xpi2+1d0
208 :
209 432710 : hy=y1-y0
210 432710 : hyi=1d0/hy
211 432710 : hy2=hy*hy
212 :
213 432710 : yp=(yget-y0)*hyi
214 432710 : ypi=1d0-yp
215 432710 : yp2=yp*yp
216 432710 : ypi2=ypi*ypi
217 :
218 432710 : cy=yp*(yp2-1d0)
219 432710 : cyi=ypi*(ypi2-1d0)
220 432710 : cyd=3d0*yp2-1d0
221 432710 : cydi=-3d0*ypi2+1d0
222 :
223 : ! bicubic spline interpolation
224 : fval = &
225 : xpi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1)) &
226 : +xp*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)) &
227 : +one_sixth*hx2*( &
228 : cxi*(ypi*fin(2,i,j) +yp*fin(2,i,j+1))+ &
229 : cx*(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1))) &
230 : +one_sixth*hy2*( &
231 : xpi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+ &
232 : xp*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1))) &
233 : +z36th*hx2*hy2*( &
234 : cxi*(cyi*fin(4,i,j) +cy*fin(4,i,j+1))+ &
235 432710 : cx*(cyi*fin(4,i+1,j)+cy*fin(4,i+1,j+1)))
236 :
237 : if (.false.) then
238 : write(*,3) 'fin(1,i,j)', i, j, fin(1,i,j)
239 : write(*,3) 'fin(1,i,j+1)', i, j+1, fin(1,i,j+1)
240 : write(*,3) 'fin(1,i+1,j)', i+1, j, fin(1,i+1,j)
241 : write(*,3) 'fin(1,i+1,j+1)', i+1, j+1, fin(1,i+1,j+1)
242 : write(*,3) 'fin(2,i,j)', i, j, fin(2,i,j)
243 : write(*,3) 'fin(2,i,j+1)', i, j+1, fin(2,i,j+1)
244 : write(*,3) 'fin(2,i+1,j)', i+1, j, fin(2,i+1,j)
245 : write(*,3) 'fin(2,i+1,j+1)', i+1, j+1, fin(2,i+1,j+1)
246 : write(*,3) 'fin(3,i,j)', i, j, fin(3,i,j)
247 : write(*,3) 'fin(3,i,j+1)', i, j+1, fin(3,i,j+1)
248 : write(*,3) 'fin(3,i+1,j)', i+1, j, fin(3,i+1,j)
249 : write(*,3) 'fin(3,i+1,j+1)', i+1, j+1, fin(3,i+1,j+1)
250 : write(*,3) 'fin(4,i,j)', i, j, fin(4,i,j)
251 : write(*,3) 'fin(4,i,j+1)', i, j+1, fin(4,i,j+1)
252 : write(*,3) 'fin(4,i+1,j)', i+1, j, fin(4,i+1,j)
253 : write(*,3) 'fin(4,i+1,j+1)', i+1, j+1, fin(4,i+1,j+1)
254 : write(*,1) 'ypi', ypi
255 : write(*,1) 'yp', yp
256 : write(*,1) 'xp', xp
257 : write(*,1) 'yp', yp
258 : write(*,1) 'hx2', hx2
259 : write(*,1) 'cxi', cxi
260 : write(*,1) 'cx', cx
261 : write(*,1) 'hy2', hy2
262 : write(*,1) 'cy', cy
263 : write(*,1) 'cyi', cyi
264 : write(*,1) 'one_sixth', one_sixth
265 : write(*,1) 'z36th', z36th
266 : write(*,1) 't1', xpi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))
267 : write(*,1) 't2', xp*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1))
268 : write(*,1) 't3', one_sixth*hx2*( &
269 : cxi*(ypi*fin(2,i,j) +yp*fin(2,i,j+1))+ &
270 : cx*(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
271 : write(*,1) 't4', one_sixth*hy2*( &
272 : xpi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+ &
273 : xp*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
274 : write(*,1) 't5', z36th*hx2*hy2*( &
275 : cxi*(cyi*fin(4,i,j) +cy*fin(4,i,j+1))+ &
276 : cx*(cyi*fin(4,i+1,j)+cy*fin(4,i+1,j+1)))
277 : write(*,1) 'fval', fval
278 : end if
279 :
280 : ! derivatives of bicubic splines
281 : df_dx = &
282 : hxi*( &
283 : -(ypi*fin(1,i,j) +yp*fin(1,i,j+1)) &
284 : +(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1))) &
285 : +one_sixth*hx*( &
286 : cxdi*(ypi*fin(2,i,j) +yp*fin(2,i,j+1))+ &
287 : cxd*(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1))) &
288 : +one_sixth*hxi*hy2*( &
289 : -(cyi*fin(3,i,j) +cy*fin(3,i,j+1)) &
290 : +(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1))) &
291 : +z36th*hx*hy2*( &
292 : cxdi*(cyi*fin(4,i,j) +cy*fin(4,i,j+1))+ &
293 432710 : cxd*(cyi*fin(4,i+1,j)+cy*fin(4,i+1,j+1)))
294 :
295 : df_dy = &
296 : hyi*( &
297 : xpi*(-fin(1,i,j) +fin(1,i,j+1))+ &
298 : xp*(-fin(1,i+1,j)+fin(1,i+1,j+1))) &
299 : +one_sixth*hx2*hyi*( &
300 : cxi*(-fin(2,i,j) +fin(2,i,j+1))+ &
301 : cx*(-fin(2,i+1,j)+fin(2,i+1,j+1))) &
302 : +one_sixth*hy*( &
303 : xpi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+ &
304 : xp*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1))) &
305 : +z36th*hx2*hy*( &
306 : cxi*(cydi*fin(4,i,j) +cyd*fin(4,i,j+1))+ &
307 432710 : cx*(cydi*fin(4,i+1,j)+cyd*fin(4,i+1,j+1)))
308 :
309 432710 : end subroutine do_kap_interpolations
310 :
311 : end module kap_eval_support
|
