Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2018-2020 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 cubic_interpolator
21 :
22 : use const_def, only: dp
23 :
24 : implicit none
25 :
26 : private
27 : public :: interpolator
28 :
29 : !4 points on either dimension for cubic interpolation
30 : integer, parameter :: num_points = 16
31 : integer, parameter :: max_order = 3
32 :
33 : type interpolator
34 : real(dp) :: matrix(num_points, num_points)
35 : real(dp) :: vector(num_points)
36 : real(dp) :: solution(num_points)
37 : logical :: solved
38 : integer :: points_added
39 : contains
40 : private
41 : procedure, public :: initialize
42 : procedure, public :: add_point
43 : procedure, public :: solve_matrix
44 : procedure, public :: evaluate
45 : procedure, public :: evaluate_deriv
46 : end type interpolator
47 :
48 : contains
49 :
50 0 : subroutine initialize(self)
51 : class(interpolator), intent(inout) :: self
52 0 : self%points_added = 0
53 0 : self%solved = .false.
54 0 : end subroutine initialize
55 :
56 0 : subroutine add_point(self, logT, logRho, val)
57 : class(interpolator), intent(inout) :: self
58 : real(dp), intent(in) :: logT
59 : real(dp), intent(in) :: logRho
60 : real(dp), intent(in) :: val
61 : integer :: column
62 : integer :: m, n
63 :
64 0 : column = 1
65 0 : do m = 0, max_order
66 0 : do n = 0, max_order
67 0 : self%matrix(self%points_added + 1, column) = logT**m*logRho**n
68 0 : column = column + 1
69 : end do
70 : end do
71 :
72 0 : self%vector(self%points_added + 1) = val
73 0 : self%points_added = self%points_added + 1
74 0 : end subroutine add_point
75 :
76 0 : subroutine solve_matrix(self)
77 : class(interpolator), intent(inout) :: self
78 : integer :: info, ipiv(num_points)
79 0 : real(dp) :: matrix_copy(num_points, num_points)
80 0 : real(dp) :: vector_copy(num_points)
81 :
82 0 : if (self%points_added /= num_points) then
83 0 : print *, "Need 16 points for cubic interpolation"
84 0 : stop
85 : end if
86 :
87 0 : matrix_copy = self%matrix
88 0 : vector_copy = self%vector
89 :
90 : !write(*,*) "Matrix", matrix_copy
91 : !write(*,*) "Vector", vector_copy
92 0 : call dgesv(num_points, 1, matrix_copy, num_points, ipiv, vector_copy, num_points, info)
93 0 : if (info /= 0) then
94 0 : write (*, *) "Matrix solve failed in interpolator", info
95 0 : write (*, *) "Matrix(:,2)", self%matrix(:, 2)
96 0 : write (*, *) "Matrix(:,5)", self%matrix(:, 5)
97 0 : stop
98 : end if
99 :
100 0 : self%solved = .true.
101 0 : self%solution = vector_copy
102 0 : end subroutine solve_matrix
103 :
104 0 : function evaluate(self, logT, logRho)
105 : class(interpolator), intent(inout) :: self
106 : real(dp), intent(in) :: logT
107 : real(dp), intent(in) :: logRho
108 : real(dp) :: evaluate
109 0 : real(dp) :: sum
110 : integer :: m, n, column
111 :
112 0 : if (.not. self%solved) then
113 0 : call self%solve_matrix()
114 : end if
115 :
116 0 : sum = 0
117 0 : column = 1
118 :
119 0 : do m = 0, max_order
120 0 : do n = 0, max_order
121 0 : sum = sum + self%solution(column)*logT**m*logRho**n
122 0 : column = column + 1
123 : end do
124 : end do
125 :
126 0 : evaluate = sum
127 0 : end function evaluate
128 :
129 0 : function evaluate_deriv(self, logT, logRho, deriv_logT, deriv_logRho)
130 : class(interpolator), intent(inout) :: self
131 : real(dp), intent(in) :: logT
132 : real(dp), intent(in) :: logRho
133 : real(dp) :: evaluate_deriv
134 0 : real(dp) :: result, product
135 : integer :: m, n, column
136 : logical :: deriv_logT, deriv_logRho
137 :
138 0 : if (.not. self%solved) then
139 0 : call self%solve_matrix()
140 : end if
141 :
142 0 : if (deriv_logT .and. deriv_logRho) then
143 0 : print *, "May choose only one derivative"
144 0 : stop
145 : end if
146 :
147 0 : if (.not. deriv_logT .and. .not. deriv_logRho) then
148 0 : print *, "Must choose one derivative"
149 0 : stop
150 : end if
151 :
152 : result = 0
153 : column = 1
154 :
155 0 : do m = 0, max_order
156 0 : do n = 0, max_order
157 0 : product = 1
158 :
159 0 : if (deriv_logT) then
160 0 : if (m == 0) then
161 : product = 0
162 : else
163 0 : product = product*m*logT**(m - 1)*logRho**n
164 : end if
165 : end if
166 :
167 0 : if (deriv_logRho) then
168 0 : if (n == 0) then
169 : product = 0
170 : else
171 0 : product = product*n*logT**m*logRho**(n - 1)
172 : end if
173 : end if
174 :
175 0 : product = product*self%solution(column)
176 0 : result = result + product
177 0 : column = column + 1
178 : end do
179 : end do
180 :
181 0 : evaluate_deriv = result
182 :
183 0 : end function evaluate_deriv
184 :
185 0 : end module cubic_interpolator
|
