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 850659 : 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 1701318 : ln_lambdas_ar, lambdas_ar, dlambdas_dlnT_ar, &
50 1701318 : ln_rlambdas_ar, rlambdas_ar, drlambdas_dlnT_ar
51 : real(dp), dimension(:), pointer :: &
52 : ln_lambdas, lambdas, dlambdas_dlnT, &
53 850659 : ln_rlambdas, rlambdas, drlambdas_dlnT
54 :
55 : include 'formats'
56 :
57 850659 : ierr = 0
58 :
59 850659 : ln_lambdas => ln_lambdas_ar
60 850659 : lambdas => lambdas_ar
61 850659 : dlambdas_dlnT => dlambdas_dlnT_ar
62 850659 : ln_rlambdas => ln_rlambdas_ar
63 850659 : rlambdas => rlambdas_ar
64 850659 : drlambdas_dlnT => drlambdas_dlnT_ar
65 :
66 : call compute_some_lambdas( &
67 850659 : num_lambdas, lo, hi, T9, rates, ln_lambdas, lambdas, dlambdas_dlnT)
68 2680200 : lambda = sum(lambdas(1:num_lambdas))
69 2680200 : dlambda_dlnT = sum(dlambdas_dlnT(1:num_lambdas))
70 :
71 850659 : 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 850659 : rlambdas, drlambdas_dlnT)
79 2680200 : rlambda = sum(rlambdas(1:num_lambdas))
80 2680200 : drlambda_dlnT = sum(drlambdas_dlnT(1:num_lambdas))
81 : end if
82 :
83 850659 : end subroutine reaction_rates
84 :
85 :
86 850659 : subroutine compute_some_lambdas( &
87 850659 : 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 : real(dp) :: T9inv, ln1
96 : real(dp), dimension(7) :: T9fac, dT9fac_dT9, dT9fac_dlnT
97 : integer :: i, j
98 :
99 : include 'formats'
100 :
101 850659 : T9inv = 1d0/T9
102 :
103 850659 : T9fac(1) = 1d0
104 850659 : dT9fac_dT9(1) = 0d0
105 :
106 850659 : T9fac(2) = T9inv
107 850659 : dT9fac_dT9(2) = -T9inv*T9inv
108 :
109 850659 : T9fac(3) = pow(T9inv,one_third)
110 850659 : dT9fac_dT9(3) = -one_third*pow(T9inv,four_thirds)
111 :
112 850659 : T9fac(4) = pow(T9,one_third)
113 850659 : dT9fac_dT9(4) = one_third*pow(T9inv,two_thirds)
114 :
115 850659 : T9fac(5) = T9
116 850659 : dT9fac_dT9(5) = 1d0
117 :
118 850659 : T9fac(6) = pow(T9,five_thirds)
119 850659 : dT9fac_dT9(6) = five_thirds*pow(T9,two_thirds)
120 :
121 850659 : T9fac(7) = log(T9)
122 850659 : dT9fac_dT9(7) = T9inv
123 :
124 6805272 : dT9fac_dlnT = T9*dT9fac_dT9
125 :
126 2680200 : do i = lo, hi
127 1829541 : j = i+1-lo
128 14636328 : ln1 = dot_product(T9fac(:), rates% coefficients(:,i))
129 2680200 : 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 1829541 : ln_lambda(j) = ln1
135 1829541 : lambda(j) = exp(ln_lambda(j))
136 1829541 : dlambda_dlnT(j) = &
137 16465869 : dot_product(dT9fac_dlnT(:), rates% coefficients(:,i))*lambda(j)
138 : end if
139 : end do
140 :
141 850659 : end subroutine compute_some_lambdas
142 :
143 :
144 850659 : subroutine compute_some_inverse_lambdas( &
145 : num_lambdas, lo, hi, T9, rates, &
146 850659 : ln_lambda, lambda, dlambda_dlnT, &
147 850659 : 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 1701318 : real(dp), dimension(num_lambdas) :: A, Qratio, dQratio_dlnT
160 : 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 850659 : indx = get_partition_fcn_indx(T9)
166 850659 : if (indx >= npart) indx = npart-1
167 850659 : if (indx < 1) then
168 844906 : Qratio(1:num_lambdas) = 1d0
169 844906 : dQratio_dlnT(1:num_lambdas) = 0d0
170 : else
171 1152880 : indxp = indx+1
172 576440 : tfac = (T9-Tpart(indx))/(Tpart(indxp)-Tpart(indx))
173 576440 : dtfac_dlnT = T9/(Tpart(indxp)-Tpart(indx))
174 1835294 : do i = lo, hi
175 1258854 : j = i+1-lo
176 1258854 : A(j) = rates% inverse_part(indxp,i)/rates% inverse_part(indx,i)
177 1258854 : blurp = pow(A(j),tfac)
178 1258854 : Qratio(j) = rates% inverse_part(indx,i) * blurp
179 1835294 : dQratio_dlnT(j) = Qratio(j)*log(A(j))*dtfac_dlnT
180 : end do
181 : end if
182 :
183 850659 : lnT9 = log(T9)
184 850659 : T9i = 1d0/T9
185 850659 : dT9i_dlnT = -T9i
186 :
187 2680200 : do i = lo, hi
188 :
189 1829541 : j = i+1-lo
190 :
191 1829541 : ln1 = ln_lambda(j) + &
192 0 : rates% inverse_coefficients(1,i) + &
193 1829541 : rates% inverse_coefficients(2,i)*T9i + &
194 5488623 : 1.5d0*rates% inverse_exp(i)*lnT9
195 :
196 1829541 : if (ln1 < ln1_max) then
197 1829541 : fac1 = exp(ln1)
198 :
199 1829541 : dln1_dlnT = dlambda_dlnT(j)/max(1d-99,lambda(j)) + &
200 : rates% inverse_coefficients(2,i)*dT9i_dlnT + &
201 1829541 : 1.5d0*rates% inverse_exp(i)
202 :
203 1829541 : 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 1829541 : inv_lambda(j) = fac1*Qratio(j)
213 1829541 : if (lambda(j) < 1d-99) then
214 145050 : dinv_lambda_dlnT(j) = 0
215 145050 : cycle
216 : end if
217 :
218 2535150 : dinv_lambda_dlnT(j) = dfac1_dlnT*Qratio(j) + fac1*dQratio_dlnT(j)
219 :
220 : end do
221 :
222 850659 : end subroutine compute_some_inverse_lambdas
223 :
224 :
225 231368 : 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 850659 : 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 231368 : call integer_dict_lookup(rates_dict, handle, indx, ierr)
235 231368 : if (ierr /= 0) indx = 0
236 231368 : end function do_reaclib_lookup
237 :
238 :
239 231358 : 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 231358 : ierr = 0
247 231358 : lo = 0
248 231358 : hi = 0
249 231358 : lo = do_reaclib_lookup(handle, rates% reaction_dict)
250 231358 : if (lo == 0) then
251 189 : ierr = -1
252 189 : return
253 : end if
254 231169 : hi = rates% nreactions
255 231169 : do ir = lo+1, rates% nreactions
256 231169 : if (rates% reaction_handle(ir) /= handle) then
257 231169 : hi = ir-1
258 231169 : exit
259 : end if
260 : end do
261 632349 : do ir = lo-1, 1, -1
262 632349 : if (rates% reaction_handle(ir) /= handle) then
263 231169 : lo = ir+1
264 231169 : exit
265 : end if
266 : end do
267 : end subroutine do_reaclib_indices_for_reaction
268 :
269 :
270 850659 : 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 850659 : ierr)
288 850659 : end subroutine do_reaclib_reaction_rates
289 :
290 :
291 :
292 : end module reaclib_eval
|
