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_ecapture
21 :
22 : use const_def
23 : use eos_def
24 : use eos_lib
25 : use rates_def
26 : use rates_lib
27 : use num_lib
28 : use math_lib
29 : use utils_lib, only: mesa_error
30 :
31 : implicit none
32 :
33 : real(dp) :: X, Z, Y, abar, zbar, z2bar, z53bar, ye
34 : integer, parameter :: species = 7
35 : integer, parameter :: h1 = 1, he4 = 2, c12 = 3, n14 = 4, o16 = 5, ne20 = 6, mg24 = 7
36 : integer, pointer, dimension(:) :: net_iso, chem_id
37 : real(dp), dimension(species) :: xa, ya, za, aa, za52
38 : integer :: handle, i, ierr
39 :
40 : real(dp) :: log10T, T, log10Rho, Rho
41 :
42 : contains
43 :
44 6 : subroutine do_test_ecapture
45 :
46 2 : call Init_Composition
47 2 : call Setup_eos(handle)
48 :
49 2 : log10Rho = 9.5d0
50 2 : Rho = exp10(log10Rho)
51 2 : log10T = 8.5d0
52 2 : T = exp10(log10T)
53 :
54 : ! check that the coulomb corrections are behaving
55 :
56 2 : write (*, '(A)')
57 2 : write (*, *) 'do_test_coulomb'
58 :
59 2 : call do_test_coulomb
60 :
61 2 : write (*, *) 'done'
62 2 : write (*, '(A)')
63 :
64 : ! check that the special weak reactions are working
65 :
66 2 : write (*, '(A)')
67 2 : write (*, *) 'do_test_special_weak'
68 :
69 2 : use_suzuki_tables = .false.
70 : ! this checks the weaklib tables
71 : ! they are extremely sparsely sampled
72 : ! so these are bad estimates of the rates
73 : ! they are shown for comparison purposes
74 2 : call do_test_special_weak(.false.)
75 :
76 : ! this shows the results using the
77 : ! on-the-fly weak rates discussed in MESA III (Section 8)
78 2 : call do_test_special_weak(.true.)
79 :
80 2 : write (*, '(A)')
81 2 : write (*, *) 'do_test_suzuki'
82 2 : use_suzuki_tables = .true.
83 : ! this shows results from a set of denser tables
84 : ! compiled by Suzuki et al. (2016)
85 : ! and mentioned in MESA V (Appendix A.2.1)
86 2 : call do_test_special_weak(.false.)
87 :
88 2 : write (*, *) 'done'
89 2 : write (*, '(A)')
90 :
91 : ! deallocate the eos tables
92 2 : call free_eos_handle(handle)
93 2 : call eos_shutdown
94 :
95 2 : end subroutine do_test_ecapture
96 :
97 2 : subroutine do_test_coulomb
98 :
99 : use auto_diff
100 : use eval_coulomb
101 : use rates_def, only: Coulomb_Info, which_mui_coulomb, which_vs_coulomb
102 :
103 : type(auto_diff_real_2var_order1) :: mu1, mu2, vs
104 2 : real(dp) :: z1, z2
105 : type(Coulomb_Info), pointer :: cc
106 :
107 : include 'formats'
108 :
109 2 : which_mui_coulomb = PCR2009
110 2 : which_vs_coulomb = Itoh2002
111 :
112 2 : z1 = 12
113 2 : z2 = 11
114 :
115 2 : allocate (cc)
116 :
117 : call coulomb_set_context(cc, T, Rho, log10T, log10Rho, &
118 2 : zbar, abar, z2bar)
119 :
120 2 : mu1 = do_mui_coulomb(cc, z1)
121 2 : mu2 = do_mui_coulomb(cc, z2)
122 :
123 2 : vs = do_vs_coulomb(cc, z1)
124 :
125 2 : deallocate (cc)
126 :
127 2 : write (*, '(6X, A4, 3F26.16)') 'mu', mu1%val, mu2%val, abs(mu1%val - mu2%val)*kev*T
128 2 : write (*, '(6X, A4, F26.16)') 'vs', vs%val
129 :
130 2 : end subroutine do_test_coulomb
131 :
132 6 : subroutine do_test_special_weak(use_special)
133 :
134 2 : use const_lib, only: const_init
135 : use utils_lib
136 : use eos_def
137 : use eos_lib
138 : use chem_def
139 : use chem_lib
140 : use rates_def, only: Coulomb_Info
141 : use eval_coulomb
142 :
143 : logical, intent(in) :: use_special
144 :
145 12 : real(dp) :: Rho, T, log10Rho, log10T
146 480 : real(dp), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT
147 132 : real(dp), dimension(num_eos_d_dxa_results, species) :: d_dxa
148 : integer :: ierr
149 :
150 : integer :: i, nr
151 6 : integer, pointer :: ids(:), reaction_ids(:)
152 : type(Coulomb_Info), pointer :: cc
153 :
154 : real(dp), dimension(:), pointer :: &
155 6 : lambda, dlambda_dlnT, dlambda_dlnRho, &
156 6 : Q, dQ_dlnT, dQ_dlnRho, &
157 6 : Qneu, dQneu_dlnT, dQneu_dlnRho
158 12 : real(dp) :: logT, T9, YeRho, &
159 6 : ye, logRho, eta, d_eta_dlnT, d_eta_dlnRho
160 : character(len=iso_name_length), dimension(2) :: weak_lhs, weak_rhs
161 :
162 6 : allocate (cc)
163 :
164 6 : nr = 2
165 :
166 : allocate ( &
167 : ids(nr), reaction_ids(nr), &
168 : lambda(nr), dlambda_dlnT(nr), dlambda_dlnRho(nr), &
169 : Q(nr), dQ_dlnT(nr), dQ_dlnRho(nr), &
170 : Qneu(nr), dQneu_dlnT(nr), dQneu_dlnRho(nr), &
171 6 : stat=ierr)
172 :
173 : ! pick reactions
174 6 : weak_lhs(1) = 'mg24'
175 6 : weak_rhs(1) = 'na24'
176 :
177 6 : weak_lhs(2) = 'na24'
178 6 : weak_rhs(2) = 'mg24'
179 :
180 18 : do i = 1, nr
181 12 : ids(i) = get_weak_rate_id(weak_lhs(i), weak_rhs(i))
182 18 : reaction_ids(i) = 0
183 : end do
184 :
185 6 : logT = 8.3d0
186 6 : logRho = 9.8d0
187 6 : Ye = 0.5d0
188 :
189 6 : T = exp10(logT)
190 6 : T9 = T*1d-9
191 6 : rho = exp10(logRho)
192 6 : YeRho = Ye*rho
193 :
194 : ! get a set of results for given temperature and density
195 : call eosDT_get( &
196 : handle, &
197 : species, chem_id, net_iso, xa, &
198 : Rho, logRho, T, logT, &
199 6 : res, d_dlnd, d_dlnT, d_dxa, ierr)
200 :
201 : call coulomb_set_context(cc, T, Rho, log10T, log10Rho, &
202 6 : zbar, abar, z2bar)
203 :
204 6 : eta = res(i_eta)
205 6 : d_eta_dlnT = d_dlnT(i_eta)
206 6 : d_eta_dlnRho = d_dlnd(i_eta)
207 :
208 6 : if (use_special) then
209 :
210 2 : do_ecapture = .true.
211 2 : ecapture_states_file = 'test_special.states'
212 2 : ecapture_transitions_file = 'test_special.transitions'
213 2 : which_mui_coulomb = PCR2009
214 2 : which_vs_coulomb = Itoh2002
215 :
216 : call eval_ecapture_reaction_info( &
217 : nr, ids, cc, T9, YeRho, &
218 : eta, d_eta_dlnT, d_eta_dlnRho, &
219 : lambda, dlambda_dlnT, dlambda_dlnRho, &
220 : Q, dQ_dlnT, dQ_dlnRho, &
221 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
222 2 : ierr)
223 :
224 2 : write (*, *) "special weak rates"
225 : else
226 4 : do_ecapture = .false.
227 : call eval_weak_reaction_info( &
228 : nr, ids, reaction_ids, cc, T9, YeRho, &
229 : eta, d_eta_dlnT, d_eta_dlnRho, &
230 : lambda, dlambda_dlnT, dlambda_dlnRho, &
231 : Q, dQ_dlnT, dQ_dlnRho, &
232 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
233 4 : ierr)
234 :
235 4 : if (use_suzuki_tables) then
236 2 : write (*, *) "suzuki weak rates"
237 : else
238 2 : write (*, *) "weaklib weak rates"
239 : end if
240 :
241 : end if
242 :
243 18 : do i = 1, nr
244 18 : write (*, '(6X, 2A6, ES26.16)') weak_lhs(i), weak_rhs(i), lambda(i)
245 : end do
246 6 : write (*, '(A)')
247 :
248 0 : deallocate ( &
249 0 : ids, reaction_ids, &
250 0 : lambda, dlambda_dlnT, dlambda_dlnRho, &
251 0 : Q, dQ_dlnT, dQ_dlnRho, &
252 6 : Qneu, dQneu_dlnT, dQneu_dlnRho, cc)
253 :
254 12 : end subroutine do_test_special_weak
255 :
256 2 : subroutine Init_Composition
257 6 : use chem_def
258 : use chem_lib
259 :
260 32 : real(dp) :: dabar_dx(species), dzbar_dx(species), &
261 18 : sumx, xh, xhe, xz, mass_correction, dmc_dx(species)
262 :
263 2 : allocate (net_iso(num_chem_isos), chem_id(species), stat=ierr)
264 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'allocate failed')
265 2 : X = 0.70d0
266 2 : Z = 0.02d0
267 :
268 15714 : net_iso(:) = 0
269 :
270 2 : chem_id(h1) = ih1; net_iso(ih1) = h1
271 2 : chem_id(he4) = ihe4; net_iso(ihe4) = he4
272 2 : chem_id(c12) = ic12; net_iso(ic12) = c12
273 2 : chem_id(n14) = in14; net_iso(in14) = n14
274 2 : chem_id(o16) = io16; net_iso(io16) = o16
275 2 : chem_id(ne20) = ine20; net_iso(ine20) = ne20
276 2 : chem_id(mg24) = img24; net_iso(img24) = mg24
277 :
278 2 : xa(h1) = 0.0d0
279 2 : xa(he4) = 0.0d0
280 2 : xa(c12) = 0.0d0
281 2 : xa(n14) = 0.0d0
282 2 : xa(o16) = 0.0d0
283 2 : xa(ne20) = 0.5d0
284 2 : xa(mg24) = 0.5d0
285 :
286 16 : do i = 1, species
287 14 : za(i) = chem_isos%Z(chem_id(i))
288 14 : aa(i) = chem_isos%W(chem_id(i))
289 14 : ya(i) = xa(i)/aa(i)
290 16 : za52(i) = pow(real(chem_isos%Z(chem_id(i)), kind=dp), 5.0d0/2.d0)
291 : end do
292 :
293 : call composition_info( &
294 : species, chem_id, xa, xh, xhe, xz, abar, zbar, z2bar, z53bar, ye, &
295 2 : mass_correction, sumx, dabar_dx, dzbar_dx, dmc_dx)
296 :
297 2 : end subroutine Init_Composition
298 :
299 2 : subroutine Setup_eos(handle)
300 : ! allocate and load the eos tables
301 2 : use eos_def
302 : use eos_lib
303 : integer, intent(out) :: handle
304 :
305 : integer :: ierr
306 : logical, parameter :: use_cache = .true.
307 :
308 2 : call eos_init('', use_cache, ierr)
309 2 : if (ierr /= 0) then
310 0 : write (*, *) 'eos_init failed in Setup_eos'
311 0 : call mesa_error(__FILE__, __LINE__)
312 : end if
313 :
314 2 : handle = alloc_eos_handle(ierr)
315 2 : if (ierr /= 0) then
316 0 : write (*, *) 'failed trying to allocate eos handle'
317 0 : call mesa_error(__FILE__, __LINE__)
318 : end if
319 :
320 2 : end subroutine Setup_eos
321 :
322 : end module test_ecapture
|
