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