Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2011-2020 Bill Paxton & 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 test_weak
21 :
22 : use rates_lib
23 : use rates_def
24 : use math_lib
25 : use const_def
26 : use utils_lib, only: mesa_error
27 : use num_lib, only: dfridr
28 :
29 : implicit none
30 :
31 : contains
32 :
33 2 : subroutine do_test_weak
34 : use chem_lib
35 : use chem_def, only: iso_name_length, chem_isos
36 : use rates_def, only: Coulomb_Info
37 : integer :: ierr, i, ir, nr
38 2 : integer, pointer :: ids(:), reaction_ids(:)
39 : type(Coulomb_Info), pointer :: cc
40 : real(dp), dimension(:), pointer :: &
41 2 : lambda, dlambda_dlnT, dlambda_dlnRho, &
42 2 : Q, dQ_dlnT, dQ_dlnRho, &
43 2 : Qneu, dQneu_dlnT, dQneu_dlnRho
44 2 : real(dp) :: logT, T, T9, YeRho, &
45 2 : ye, rho, logRho, eta, d_eta_dlnT, d_eta_dlnRho
46 : character(len=iso_name_length) :: weak_lhs, weak_rhs
47 : character(len=2*iso_name_length + 1) :: key
48 :
49 2 : real(dp) :: dvardx, dvardx_0, dx_0, err, var_0, xdum
50 : logical :: doing_d_dlnd
51 :
52 : include 'formats'
53 :
54 2 : ierr = 0
55 :
56 2 : write (*, *) 'check weak_info_list'
57 2 : weak_lhs = 'o14'
58 2 : weak_rhs = 'n14'
59 2 : i = get_weak_info_list_id(weak_lhs, weak_rhs)
60 2 : if (i <= 0 .or. i > num_weak_info_list_reactions) then
61 0 : write (*, *) 'get_weak_info_list_id failed'
62 0 : call mesa_error(__FILE__, __LINE__)
63 : end if
64 2 : call create_weak_dict_key(weak_lhs, weak_rhs, key)
65 2 : write (*, 1) trim(key)
66 2 : write (*, 1) 'halflife', weak_info_list_halflife(i)
67 2 : write (*, 1) 'Qneu', weak_info_list_Qneu(i)
68 2 : write (*, '(A)')
69 :
70 2 : d_eta_dlnT = 0
71 2 : d_eta_dlnRho = 0
72 :
73 : if (.false.) then ! TESTING
74 : logT = 7.5904236599874348D+00
75 : logRho = 1.0657946486820271D+00
76 : ye = 8.2724691280321605D-01
77 : eta = -5.3262903257381922D+00
78 : d_eta_dlnT = -1.5299344982339016D+00
79 : d_eta_dlnRho = 9.9482489248846617D-01
80 : else if (.true.) then ! TESTING
81 2 : logT = log10(9.0d8)
82 2 : ye = 0.5d0
83 2 : logRho = log10(4.5d5)
84 2 : eta = 10d0
85 :
86 : ! call for cell 565
87 : ! logT = 8.3534130765231005
88 : ! logRho = 2.4507395003327828
89 : ! T = 225638433.79026267
90 : ! Rho = 282.31860559981624
91 : ! abar = 4.0424973056746829
92 : ! zbar = 2.0238390731055702
93 : ! z2bar = 4.2987387813744071
94 : ! ye = 0.50064079703023978
95 : ! eta = -5.3287260155378711
96 : ! d_eta_dlnT= -1.5713400060794886
97 : ! d_eta_dlnRho = 1.0016532307086357
98 :
99 : else
100 : logT = 7.5d0
101 : logRho = 4d0
102 : Ye = 0.5d0
103 : eta = 0d0
104 : end if
105 :
106 2 : T = exp10(logT)
107 2 : T9 = T*1d-9
108 2 : rho = exp10(logRho)
109 2 : YeRho = Ye*rho
110 :
111 2 : write (*, 1) 'logT', logT
112 2 : write (*, 1) 'logRho', logRho
113 2 : write (*, 1) 'ye', ye
114 2 : write (*, 1) 'eta', eta
115 2 : write (*, 1) 'T9', T9
116 2 : write (*, 1) 'lYeRho', log10(YeRho)
117 2 : write (*, '(A)')
118 :
119 2 : nr = num_weak_reactions
120 : allocate ( &
121 : ids(nr), reaction_ids(nr), &
122 : lambda(nr), dlambda_dlnT(nr), dlambda_dlnRho(nr), &
123 : Q(nr), dQ_dlnT(nr), dQ_dlnRho(nr), &
124 : Qneu(nr), dQneu_dlnT(nr), dQneu_dlnRho(nr), &
125 2 : stat=ierr)
126 2 : if (ierr /= 0) return
127 928 : do i = 1, nr
128 926 : ids(i) = i
129 928 : reaction_ids(i) = i
130 : end do
131 :
132 2 : write (*, '(A)')
133 2 : write (*, 2) 'nr', nr
134 2 : write (*, '(A)')
135 :
136 : call eval_weak_reaction_info( &
137 : nr, ids, reaction_ids, cc, T9, YeRho, &
138 : eta, d_eta_dlnT, d_eta_dlnRho, &
139 : lambda, dlambda_dlnT, dlambda_dlnRho, &
140 : Q, dQ_dlnT, dQ_dlnRho, &
141 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
142 2 : ierr)
143 2 : if (ierr /= 0) then
144 0 : write (*, *) 'failed in eval_weak_reaction_info'
145 0 : call mesa_error(__FILE__, __LINE__)
146 : end if
147 :
148 : if (.true.) then
149 : write (*, '(30x,99a16)') &
150 2 : 'halflife', 'Qneu', 'Qtotal'
151 928 : do i = 1, nr
152 926 : ir = ids(i)
153 926 : if (ir <= 0 .or. ir > size(weak_lhs_nuclide_id)) then
154 0 : write (*, *) 'ir', ir
155 0 : call mesa_error(__FILE__, __LINE__)
156 : end if
157 926 : if (Qneu(i) < 1d-20) cycle
158 662 : if (lambda(i) < 1d-20) cycle
159 564 : weak_lhs = chem_isos%name(weak_lhs_nuclide_id(ir))
160 564 : weak_rhs = chem_isos%name(weak_rhs_nuclide_id(ir))
161 564 : write (*, '(a30,99(1pe16.6))') weak_lhs//weak_rhs, &
162 1130 : ln2/lambda(i), Qneu(i), Q(i)
163 : end do
164 2 : write (*, '(A)')
165 : else
166 : write (*, '(30x,5a12,a20)') 'Q', 'Qneu', 'lambda'
167 : do i = 1, nr
168 : ir = ids(i)
169 : if (ir <= 0 .or. ir > size(weak_lhs_nuclide_id)) then
170 : write (*, *) 'ir', ir
171 : call mesa_error(__FILE__, __LINE__)
172 : end if
173 : weak_lhs = chem_isos%name(weak_lhs_nuclide_id(ir))
174 : weak_rhs = chem_isos%name(weak_rhs_nuclide_id(ir))
175 : write (*, '(a30,5f12.6,e20.12)') weak_lhs//weak_rhs, &
176 : Q(i), Qneu(i), lambda(i)
177 : end do
178 : write (*, '(A)')
179 :
180 : if (.false.) then
181 : write (*, '(a30,5a12,a20)') 'd_dT9', 'Q', 'Qneu', 'lambda'
182 : do i = 1, nr
183 : ir = ids(i)
184 : weak_lhs = chem_isos%name(weak_lhs_nuclide_id(ir))
185 : weak_rhs = chem_isos%name(weak_rhs_nuclide_id(ir))
186 : write (*, '(a30,5f12.6,e20.12)') weak_lhs//weak_rhs, &
187 : dQ_dlnT(i), dQneu_dlnT(i), dlambda_dlnT(i)
188 : end do
189 : write (*, '(A)')
190 : end if
191 :
192 : if (.false.) then
193 : write (*, '(a30,5a12,a20)') 'd_d_rho', 'Q', 'Qneu', 'lambda'
194 : do i = 1, nr
195 : ir = ids(i)
196 : weak_lhs = chem_isos%name(weak_lhs_nuclide_id(ir))
197 : weak_rhs = chem_isos%name(weak_rhs_nuclide_id(ir))
198 : write (*, '(a30,5f12.6,e20.12)') weak_lhs//weak_rhs, &
199 : dQ_dlnRho(i), dQneu_dlnRho(i), dlambda_dlnRho(i)
200 : end do
201 : write (*, '(A)')
202 : end if
203 :
204 : end if
205 :
206 2 : write (*, *) 'done'
207 :
208 : if (.false.) then ! dfridr tests for partials
209 :
210 : do i = 1, num_weak_reactions
211 :
212 : ir = ids(i)
213 : weak_lhs = chem_isos%name(weak_lhs_nuclide_id(ir))
214 : weak_rhs = chem_isos%name(weak_rhs_nuclide_id(ir))
215 : write (*, '(a30)') weak_lhs//weak_rhs
216 :
217 : doing_d_dlnd = .true.
218 : doing_d_dlnd = .false.
219 :
220 : var_0 = lambda(i)
221 :
222 : if (doing_d_dlnd) then
223 : dx_0 = max(1d-14, abs(logRho*ln10*1d-6))
224 : dvardx_0 = dlambda_dlnRho(i)
225 : else
226 : dx_0 = max(1d-14, abs(logT*ln10*1d-6))
227 : dvardx_0 = dlambda_dlnT(i)
228 : end if
229 : err = 0d0
230 : dvardx = dfridr(dx_0, dfridr_weak_reaction_info, err)
231 : xdum = (dvardx - dvardx_0)/max(abs(dvardx_0), 1d-50)
232 : write (*, 1) 'analytic, numeric, est err in numeric, rel diff', &
233 : dvardx_0, dvardx, err, xdum
234 : if (doing_d_dlnd) then
235 : write (*, *) 'doing dlnd'
236 : else ! doing d_dlnT
237 : write (*, *) 'doing dlnT'
238 : end if
239 : write (*, *) 'test net'
240 : write (*, '(A)')
241 :
242 : end do
243 :
244 : call mesa_error(__FILE__, __LINE__, 'test rate')
245 :
246 : end if
247 :
248 0 : deallocate ( &
249 0 : ids, reaction_ids, &
250 0 : lambda, dlambda_dlnT, dlambda_dlnRho, &
251 0 : Q, dQ_dlnT, dQ_dlnRho, &
252 4 : Qneu, dQneu_dlnT, dQneu_dlnRho)
253 :
254 : contains
255 :
256 0 : real(dp) function dfridr_weak_reaction_info(delta_x) result(val)
257 : real(dp), intent(in) :: delta_x
258 : integer :: ierr
259 0 : real(dp) :: logYeRho, logT, var, log_var
260 : include 'formats'
261 0 : ierr = 0
262 :
263 0 : if (doing_d_dlnd) then
264 :
265 0 : logYeRho = log10(YeRho)
266 0 : log_var = logYeRho + delta_x/ln10
267 0 : var = exp10(log_var)
268 :
269 : call eval_weak_reaction_info( &
270 : nr, ids, reaction_ids, cc, T9, var, &
271 : eta, d_eta_dlnT, d_eta_dlnRho, &
272 : lambda, dlambda_dlnT, dlambda_dlnRho, &
273 : Q, dQ_dlnT, dQ_dlnRho, &
274 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
275 0 : ierr)
276 :
277 : else
278 :
279 0 : logT = log10(T9) + 9d0
280 0 : log_var = logT + delta_x/ln10
281 0 : var = exp10(log_var - 9d0)
282 :
283 : call eval_weak_reaction_info( &
284 : nr, ids, reaction_ids, cc, var, YeRho, &
285 : eta, d_eta_dlnT, d_eta_dlnRho, &
286 : lambda, dlambda_dlnT, dlambda_dlnRho, &
287 : Q, dQ_dlnT, dQ_dlnRho, &
288 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
289 0 : ierr)
290 :
291 : end if
292 :
293 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'failed in call eval_weak_reaction_info')
294 0 : val = lambda(i)
295 :
296 2 : end function dfridr_weak_reaction_info
297 :
298 : end subroutine do_test_weak
299 :
300 : end module test_weak
|
