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 reaclib_eval
21 : use rates_def
22 : use math_lib
23 :
24 : implicit none
25 :
26 : real(dp), parameter :: &
27 : lam_max = 1d99, &
28 : ln1_max = 227.955924206411d0 ! log(lam_max)
29 :
30 :
31 : contains
32 :
33 :
34 852157 : subroutine reaction_rates( &
35 : num_lambdas, lo, hi, T9, rates, nuclides, forward_only, &
36 : lambda, dlambda_dlnT, &
37 : rlambda, drlambda_dlnT, &
38 : ierr)
39 : integer, intent(in) :: num_lambdas, lo, hi
40 : real(dp), intent(in) :: T9
41 : type(reaction_data), intent(in) :: rates
42 : type(nuclide_data), intent(in) :: nuclides
43 : logical, intent(in) :: forward_only
44 : real(dp), intent(out) :: lambda, dlambda_dlnT
45 : real(dp), intent(out) :: rlambda, drlambda_dlnT
46 : integer, intent(out) :: ierr
47 :
48 : real(dp), dimension(num_lambdas), target :: &
49 8054502 : ln_lambdas_ar, lambdas_ar, dlambdas_dlnT_ar, &
50 8054502 : ln_rlambdas_ar, rlambdas_ar, drlambdas_dlnT_ar
51 : real(dp), dimension(:), pointer :: &
52 : ln_lambdas, lambdas, dlambdas_dlnT, &
53 852157 : ln_rlambdas, rlambdas, drlambdas_dlnT
54 :
55 : include 'formats'
56 :
57 852157 : ierr = 0
58 :
59 852157 : ln_lambdas => ln_lambdas_ar
60 852157 : lambdas => lambdas_ar
61 852157 : dlambdas_dlnT => dlambdas_dlnT_ar
62 852157 : ln_rlambdas => ln_rlambdas_ar
63 852157 : rlambdas => rlambdas_ar
64 852157 : drlambdas_dlnT => drlambdas_dlnT_ar
65 :
66 : call compute_some_lambdas( &
67 852157 : num_lambdas, lo, hi, T9, rates, ln_lambdas, lambdas, dlambdas_dlnT)
68 2684834 : lambda = sum(lambdas(1:num_lambdas))
69 2684834 : dlambda_dlnT = sum(dlambdas_dlnT(1:num_lambdas))
70 :
71 852157 : if (forward_only) then
72 0 : rlambda = 0
73 0 : drlambda_dlnT = 0
74 : else
75 : call compute_some_inverse_lambdas( &
76 : num_lambdas, lo, hi, T9, rates, &
77 : ln_lambdas, lambdas, dlambdas_dlnT, &
78 852157 : rlambdas, drlambdas_dlnT)
79 2684834 : rlambda = sum(rlambdas(1:num_lambdas))
80 2684834 : drlambda_dlnT = sum(drlambdas_dlnT(1:num_lambdas))
81 : end if
82 :
83 852157 : end subroutine reaction_rates
84 :
85 :
86 852157 : subroutine compute_some_lambdas( &
87 852157 : num_lambdas, lo, hi, T9, rates, ln_lambda, lambda, dlambda_dlnT)
88 : use utils_lib, only: is_bad
89 : use const_def, only: dp, one_third, two_thirds, four_thirds, five_thirds
90 : integer, intent(in) :: num_lambdas, lo, hi ! range of rates to do
91 : real(dp), intent(in) :: T9
92 : type(reaction_data), intent(in) :: rates
93 : real(dp), dimension(:), intent(out) :: ln_lambda, lambda, dlambda_dlnT
94 :
95 6817256 : real(dp) :: T9inv, ln1
96 18747454 : real(dp), dimension(7) :: T9fac, dT9fac_dT9, dT9fac_dlnT
97 : integer :: i, j
98 :
99 : include 'formats'
100 :
101 852157 : T9inv = 1d0/T9
102 :
103 852157 : T9fac(1) = 1d0
104 852157 : dT9fac_dT9(1) = 0d0
105 :
106 852157 : T9fac(2) = T9inv
107 852157 : dT9fac_dT9(2) = -T9inv*T9inv
108 :
109 852157 : T9fac(3) = pow(T9inv,one_third)
110 852157 : dT9fac_dT9(3) = -one_third*pow(T9inv,four_thirds)
111 :
112 852157 : T9fac(4) = pow(T9,one_third)
113 852157 : dT9fac_dT9(4) = one_third*pow(T9inv,two_thirds)
114 :
115 852157 : T9fac(5) = T9
116 852157 : dT9fac_dT9(5) = 1d0
117 :
118 852157 : T9fac(6) = pow(T9,five_thirds)
119 852157 : dT9fac_dT9(6) = five_thirds*pow(T9,two_thirds)
120 :
121 852157 : T9fac(7) = log(T9)
122 852157 : dT9fac_dT9(7) = T9inv
123 :
124 6817256 : dT9fac_dlnT = T9*dT9fac_dT9
125 :
126 2684834 : do i = lo, hi
127 1832677 : j = i+1-lo
128 14661416 : ln1 = dot_product(T9fac(:), rates% coefficients(:,i))
129 2684834 : if (ln1 > ln1_max) then
130 0 : ln_lambda(j) = ln1_max
131 0 : lambda(j) = lam_max ! == exp(ln1_max)
132 0 : dlambda_dlnT(j) = 0
133 : else
134 1832677 : ln_lambda(j) = ln1
135 1832677 : lambda(j) = exp(ln_lambda(j))
136 : dlambda_dlnT(j) = &
137 14661416 : dot_product(dT9fac_dlnT(:), rates% coefficients(:,i))*lambda(j)
138 : end if
139 : end do
140 :
141 852157 : end subroutine compute_some_lambdas
142 :
143 :
144 852157 : subroutine compute_some_inverse_lambdas( &
145 : num_lambdas, lo, hi, T9, rates, &
146 852157 : ln_lambda, lambda, dlambda_dlnT, &
147 852157 : inv_lambda, dinv_lambda_dlnT)
148 : use utils_lib, only: is_bad
149 : use chem_def, only: Tpart, npart
150 : use chem_lib, only: get_partition_fcn_indx
151 : integer, intent(in) :: num_lambdas, lo, hi ! range of rates to do
152 : real(dp), intent(in) :: T9
153 : type(reaction_data), intent(in) :: rates
154 : real(dp), dimension(:), intent(in) :: ln_lambda, lambda, dlambda_dlnT
155 : real(dp), dimension(:), intent(out) :: inv_lambda, dinv_lambda_dlnT
156 :
157 : integer :: indx,indxp
158 : integer :: i, j
159 7202345 : real(dp), dimension(num_lambdas) :: A, Qratio, dQratio_dlnT
160 4389148 : real(dp) :: tfac, dtfac_dlnT, lnT9, T9i, dT9i_dlnT, ln1, fac1, dfac1_dlnT, dln1_dlnT,blurp
161 :
162 : include 'formats'
163 :
164 : ! find index of partition function and logarithmically interpolate
165 852157 : indx = get_partition_fcn_indx(T9)
166 852157 : if (indx >= npart) indx = npart-1
167 852157 : if (indx < 1) then
168 846240 : Qratio(1:num_lambdas) = 1d0
169 846240 : dQratio_dlnT(1:num_lambdas) = 0d0
170 : else
171 577505 : indxp = indx+1
172 577505 : tfac = (T9-Tpart(indx))/(Tpart(indxp)-Tpart(indx))
173 577505 : dtfac_dlnT = T9/(Tpart(indxp)-Tpart(indx))
174 1838594 : do i = lo, hi
175 1261089 : j = i+1-lo
176 1261089 : A(j) = rates% inverse_part(indxp,i)/rates% inverse_part(indx,i)
177 1261089 : blurp = pow(A(j),tfac)
178 1261089 : Qratio(j) = rates% inverse_part(indx,i) * blurp
179 1838594 : dQratio_dlnT(j) = Qratio(j)*log(A(j))*dtfac_dlnT
180 : end do
181 : end if
182 :
183 852157 : lnT9 = log(T9)
184 852157 : T9i = 1d0/T9
185 852157 : dT9i_dlnT = -T9i
186 :
187 2684834 : do i = lo, hi
188 :
189 1832677 : j = i+1-lo
190 :
191 : ln1 = ln_lambda(j) + &
192 : rates% inverse_coefficients(1,i) + &
193 : rates% inverse_coefficients(2,i)*T9i + &
194 1832677 : 1.5d0*rates% inverse_exp(i)*lnT9
195 :
196 1832677 : if (ln1 < ln1_max) then
197 1832677 : fac1 = exp(ln1)
198 :
199 : dln1_dlnT = dlambda_dlnT(j)/max(1d-99,lambda(j)) + &
200 : rates% inverse_coefficients(2,i)*dT9i_dlnT + &
201 1832677 : 1.5d0*rates% inverse_exp(i)
202 :
203 1832677 : dfac1_dlnT = dln1_dlnT*fac1
204 :
205 : else
206 0 : ln1 = ln1_max
207 0 : fac1 = lam_max ! == exp(ln1_max)
208 0 : dln1_dlnT = 0
209 0 : dfac1_dlnT = 0
210 : end if
211 :
212 1832677 : inv_lambda(j) = fac1*Qratio(j)
213 1832677 : if (lambda(j) < 1d-99) then
214 145358 : dinv_lambda_dlnT(j) = 0
215 145358 : cycle
216 : end if
217 :
218 2539476 : dinv_lambda_dlnT(j) = dfac1_dlnT*Qratio(j) + fac1*dQratio_dlnT(j)
219 :
220 : end do
221 :
222 852157 : end subroutine compute_some_inverse_lambdas
223 :
224 :
225 232616 : integer function do_reaclib_lookup(handle, rates_dict) result(indx)
226 : ! returns first reaction index that matches handle.
227 : ! there may be several following that one having the same handle.
228 : ! returns 0 if handle doesn't match any of the reactions
229 852157 : use utils_lib, only: integer_dict_lookup
230 : character(len=*), intent(in) :: handle ! as in rates% reaction_handle
231 : type (integer_dict), pointer :: rates_dict ! from create_reaclib_rates_dict
232 : integer :: ierr
233 : ierr = 0
234 232616 : call integer_dict_lookup(rates_dict, handle, indx, ierr)
235 232616 : if (ierr /= 0) indx = 0
236 232616 : end function do_reaclib_lookup
237 :
238 :
239 232596 : subroutine do_reaclib_indices_for_reaction(handle, rates, lo, hi, ierr)
240 : character(len=*), intent(in) :: handle ! as in rates% reaction_handle
241 : type(reaction_data), intent(in) :: rates
242 : integer, intent(out) :: lo, hi
243 : integer, intent(out) :: ierr
244 : integer :: ir
245 : include 'formats'
246 232596 : ierr = 0
247 232596 : lo = 0
248 232596 : hi = 0
249 232596 : lo = do_reaclib_lookup(handle, rates% reaction_dict)
250 232596 : if (lo == 0) then
251 339 : ierr = -1
252 339 : return
253 : end if
254 232257 : hi = rates% nreactions
255 232257 : do ir = lo+1, rates% nreactions
256 232257 : if (rates% reaction_handle(ir) /= handle) then
257 232257 : hi = ir-1
258 232257 : exit
259 : end if
260 : end do
261 634517 : do ir = lo-1, 1, -1
262 634517 : if (rates% reaction_handle(ir) /= handle) then
263 232257 : lo = ir+1
264 232257 : exit
265 : end if
266 : end do
267 : end subroutine do_reaclib_indices_for_reaction
268 :
269 :
270 852157 : subroutine do_reaclib_reaction_rates( &
271 : lo, hi, T9, rates, nuclides, forward_only, &
272 : lambda, dlambda_dlnT, &
273 : rlambda, drlambda_dlnT, &
274 : ierr)
275 : integer, intent(in) :: lo, hi ! from reaclib_indices_for_reaction
276 : real(dp), intent(in) :: T9
277 : type(reaction_data), intent(in) :: rates
278 : type(nuclide_data), intent(in) :: nuclides
279 : logical, intent(in) :: forward_only
280 : real(dp), intent(out) :: lambda, dlambda_dlnT
281 : real(dp), intent(out) :: rlambda, drlambda_dlnT
282 : integer, intent(out) :: ierr
283 : call reaction_rates( &
284 : hi-lo+1, lo, hi, T9, rates, nuclides, forward_only, &
285 : lambda, dlambda_dlnT, &
286 : rlambda, drlambda_dlnT, &
287 852157 : ierr)
288 852157 : end subroutine do_reaclib_reaction_rates
289 :
290 :
291 :
292 : end module reaclib_eval
|
