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 3 : subroutine do_test_ecapture
45 :
46 1 : call Init_Composition
47 1 : call Setup_eos(handle)
48 :
49 1 : log10Rho = 9.5d0
50 1 : Rho = exp10(log10Rho)
51 1 : log10T = 8.5d0
52 1 : T = exp10(log10T)
53 :
54 : ! check that the coulomb corrections are behaving
55 :
56 1 : write (*, '(A)')
57 1 : write (*, *) 'do_test_coulomb'
58 :
59 1 : call do_test_coulomb
60 :
61 1 : write (*, *) 'done'
62 1 : write (*, '(A)')
63 :
64 : ! check that the special weak reactions are working
65 :
66 1 : write (*, '(A)')
67 1 : write (*, *) 'do_test_special_weak'
68 :
69 1 : 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 1 : 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 1 : call do_test_special_weak(.true.)
79 :
80 1 : write (*, '(A)')
81 1 : write (*, *) 'do_test_suzuki'
82 1 : 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 1 : call do_test_special_weak(.false.)
87 :
88 1 : write (*, *) 'done'
89 1 : write (*, '(A)')
90 :
91 : ! deallocate the eos tables
92 1 : call free_eos_handle(handle)
93 1 : call eos_shutdown
94 :
95 1 : end subroutine do_test_ecapture
96 :
97 1 : 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 : real(dp) :: z1, z2
105 : type(Coulomb_Info), pointer :: cc
106 :
107 : include 'formats'
108 :
109 1 : which_mui_coulomb = PCR2009
110 1 : which_vs_coulomb = Itoh2002
111 :
112 1 : z1 = 12
113 1 : z2 = 11
114 :
115 1 : allocate (cc)
116 :
117 : call coulomb_set_context(cc, T, Rho, log10T, log10Rho, &
118 1 : zbar, abar, z2bar)
119 :
120 1 : mu1 = do_mui_coulomb(cc, z1)
121 1 : mu2 = do_mui_coulomb(cc, z2)
122 :
123 1 : vs = do_vs_coulomb(cc, z1)
124 :
125 1 : deallocate (cc)
126 :
127 1 : write (*, '(6X, A4, 3F26.16)') 'mu', mu1%val, mu2%val, abs(mu1%val - mu2%val)*kev*T
128 1 : write (*, '(6X, A4, F26.16)') 'vs', vs%val
129 :
130 1 : end subroutine do_test_coulomb
131 :
132 3 : subroutine do_test_special_weak(use_special)
133 :
134 1 : 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 : real(dp) :: Rho, T, log10Rho, log10T
146 : real(dp), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT
147 : real(dp), dimension(num_eos_d_dxa_results, species) :: d_dxa
148 : integer :: ierr
149 :
150 : integer :: i, nr
151 3 : integer, pointer :: ids(:), reaction_ids(:)
152 : type(Coulomb_Info), pointer :: cc
153 :
154 : real(dp), dimension(:), pointer :: &
155 3 : lambda, dlambda_dlnT, dlambda_dlnRho, &
156 3 : Q, dQ_dlnT, dQ_dlnRho, &
157 3 : Qneu, dQneu_dlnT, dQneu_dlnRho
158 : real(dp) :: logT, T9, YeRho, &
159 : ye, logRho, eta, d_eta_dlnT, d_eta_dlnRho
160 : character(len=iso_name_length), dimension(2) :: weak_lhs, weak_rhs
161 :
162 0 : allocate (cc)
163 :
164 3 : 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 3 : stat=ierr)
172 :
173 : ! pick reactions
174 3 : weak_lhs(1) = 'mg24'
175 3 : weak_rhs(1) = 'na24'
176 :
177 3 : weak_lhs(2) = 'na24'
178 3 : weak_rhs(2) = 'mg24'
179 :
180 9 : do i = 1, nr
181 6 : ids(i) = get_weak_rate_id(weak_lhs(i), weak_rhs(i))
182 9 : reaction_ids(i) = 0
183 : end do
184 :
185 3 : logT = 8.3d0
186 3 : logRho = 9.8d0
187 3 : Ye = 0.5d0
188 :
189 3 : T = exp10(logT)
190 3 : T9 = T*1d-9
191 3 : rho = exp10(logRho)
192 3 : 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 3 : res, d_dlnd, d_dlnT, d_dxa, ierr)
200 :
201 : call coulomb_set_context(cc, T, Rho, log10T, log10Rho, &
202 3 : zbar, abar, z2bar)
203 :
204 3 : eta = res(i_eta)
205 3 : d_eta_dlnT = d_dlnT(i_eta)
206 3 : d_eta_dlnRho = d_dlnd(i_eta)
207 :
208 3 : if (use_special) then
209 :
210 1 : do_ecapture = .true.
211 1 : ecapture_states_file = 'test_special.states'
212 1 : ecapture_transitions_file = 'test_special.transitions'
213 1 : which_mui_coulomb = PCR2009
214 1 : 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 1 : ierr)
223 :
224 1 : write (*, *) "special weak rates"
225 : else
226 2 : 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 2 : ierr)
234 :
235 2 : if (use_suzuki_tables) then
236 1 : write (*, *) "suzuki weak rates"
237 : else
238 1 : write (*, *) "weaklib weak rates"
239 : end if
240 :
241 : end if
242 :
243 9 : do i = 1, nr
244 9 : write (*, '(6X, 2A6, ES26.16)') weak_lhs(i), weak_rhs(i), lambda(i)
245 : end do
246 3 : 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 3 : Qneu, dQneu_dlnT, dQneu_dlnRho, cc)
253 :
254 6 : end subroutine do_test_special_weak
255 :
256 1 : subroutine Init_Composition
257 3 : use chem_def
258 : use chem_lib
259 :
260 : real(dp) :: dabar_dx(species), dzbar_dx(species), &
261 : sumx, xh, xhe, xz, mass_correction, dmc_dx(species)
262 :
263 1 : allocate (net_iso(num_chem_isos), chem_id(species), stat=ierr)
264 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'allocate failed')
265 1 : X = 0.70d0
266 1 : Z = 0.02d0
267 :
268 7857 : net_iso(:) = 0
269 :
270 1 : chem_id(h1) = ih1; net_iso(ih1) = h1
271 1 : chem_id(he4) = ihe4; net_iso(ihe4) = he4
272 1 : chem_id(c12) = ic12; net_iso(ic12) = c12
273 1 : chem_id(n14) = in14; net_iso(in14) = n14
274 1 : chem_id(o16) = io16; net_iso(io16) = o16
275 1 : chem_id(ne20) = ine20; net_iso(ine20) = ne20
276 1 : chem_id(mg24) = img24; net_iso(img24) = mg24
277 :
278 1 : xa(h1) = 0.0d0
279 1 : xa(he4) = 0.0d0
280 1 : xa(c12) = 0.0d0
281 1 : xa(n14) = 0.0d0
282 1 : xa(o16) = 0.0d0
283 1 : xa(ne20) = 0.5d0
284 1 : xa(mg24) = 0.5d0
285 :
286 8 : do i = 1, species
287 7 : za(i) = chem_isos%Z(chem_id(i))
288 7 : aa(i) = chem_isos%W(chem_id(i))
289 7 : ya(i) = xa(i)/aa(i)
290 8 : 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 1 : mass_correction, sumx, dabar_dx, dzbar_dx, dmc_dx)
296 :
297 1 : end subroutine Init_Composition
298 :
299 1 : subroutine Setup_eos(handle)
300 : ! allocate and load the eos tables
301 1 : use eos_def
302 : use eos_lib
303 : integer, intent(out) :: handle
304 :
305 : integer :: ierr
306 : logical, parameter :: use_cache = .true.
307 :
308 1 : call eos_init('', use_cache, ierr)
309 1 : 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 1 : handle = alloc_eos_handle(ierr)
315 1 : 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 1 : end subroutine Setup_eos
321 :
322 : end module test_ecapture
|
