Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 Pablo Marchant & 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 :
21 : module binary_ce
22 :
23 : use const_def, only: dp, avo, secyer, boltzm, standard_cgrav, crad, ev2erg, rsun, msun
24 : use math_lib
25 : use star_lib
26 : use star_def
27 : use binary_def
28 : use interp_1d_def, only: pm_work_size
29 : use interp_1d_lib, only: interp_pm, interp_values, interp_value
30 :
31 : implicit none
32 :
33 : contains
34 :
35 0 : subroutine CE_init(b, restart, ierr)
36 : use interp_1d_def, only: pm_work_size
37 : use interp_1d_lib, only: interp_pm
38 : type (binary_info), pointer :: b
39 : logical, intent(in) :: restart
40 : integer, intent(out) :: ierr
41 : type (star_info), pointer :: s
42 0 : real(dp), pointer :: interp_work(:), adjusted_energy(:)
43 : integer :: k, op_err
44 0 : real(dp) :: rec_energy_HII_to_HI, &
45 0 : rec_energy_HeII_to_HeI, &
46 0 : rec_energy_HeIII_to_HeII, &
47 0 : diss_energy_H2, &
48 0 : frac_HI, frac_HII, &
49 0 : frac_HeI, frac_HeII, frac_HeIII, &
50 0 : avg_charge_He, energy_comp
51 : include 'formats'
52 :
53 : ! TODO: no care is taken in here when model_twins_flag is true
54 : ! not a priority, but needs to be sorted out whenever double core
55 : ! evolution is implemented
56 :
57 0 : ierr = 0
58 :
59 0 : if (b% use_other_CE_init) then
60 0 : call b% other_CE_init(b% binary_id, restart, ierr)
61 0 : return
62 : end if
63 :
64 0 : write(*,*) "TURNING ON CE"
65 :
66 0 : b% s_donor% mix_factor = 0d0
67 0 : b% s_donor% dxdt_nuc_factor = 0d0
68 0 : if (b% point_mass_i == 0) then
69 0 : b% s_accretor% mix_factor = 0d0
70 0 : b% s_accretor% dxdt_nuc_factor = 0d0
71 : end if
72 :
73 0 : if (b% d_i == 1) then
74 0 : b% CE_num1 = b% CE_num1 + 1
75 : else
76 0 : b% CE_num2 = b% CE_num2 + 1
77 : end if
78 :
79 0 : b% keep_donor_fixed = .true.
80 :
81 0 : b% CE_init = .true.
82 :
83 0 : if (.not. restart) then
84 0 : b% CE_years_detached = 0d0
85 0 : s => b% s_donor
86 :
87 0 : write(*,*) "Initiating common envelope phase"
88 :
89 0 : allocate(b% CE_m(s% nz), b% CE_entropy(4*s% nz), stat=ierr)
90 0 : if(ierr /= 0) then
91 0 : write(*,*) "CE_init: Error during allocation"
92 0 : return
93 : end if
94 0 : allocate(b% CE_U_in(4*s% nz), b% CE_U_out(4*s% nz), b% CE_Omega_in(4*s% nz), b% CE_Omega_out(4*s% nz), stat=ierr)
95 0 : if(ierr /= 0) then
96 0 : write(*,*) "CE_init: Error during allocation"
97 0 : return
98 : end if
99 :
100 : ! get energy from the EOS and adjust the different contributions from recombination/dissociation
101 0 : allocate(adjusted_energy(s% nz))
102 0 : do k=1, s% nz
103 : ! the following lines compute the fractions of HI, HII, HeI, HeII and HeIII
104 : ! things like ion_ifneut_H are defined in $MESA_DIR/ionization/public/ionization.def
105 : ! this file can be checked for additional ionization output available
106 0 : frac_HI = 0d0 !get_ion_info(s,ion_ifneut_H,k)
107 0 : frac_HII = 1 - frac_HI
108 :
109 : ! ionization module provides neutral fraction and average charge of He.
110 : ! use these two to compute the mass fractions of HeI and HeII
111 0 : frac_HeI = 0d0 !get_ion_info(s,ion_ifneut_He,k)
112 0 : avg_charge_He = 2d0 !get_ion_info(s,ion_iZ_He,k)
113 : ! the following is the solution to the equations
114 : ! avg_charge_He = 2*fracHeIII + 1*fracHeII
115 : ! 1 = fracHeI + fracHeII + fracHeIII
116 0 : frac_HeII = 2 - 2*frac_HeI - avg_charge_He
117 0 : frac_HeIII = 1 - frac_HeII - frac_HeI
118 :
119 : ! recombination energies from https://physics.nist.gov/PhysRefData/ASD/ionEnergy.html
120 0 : rec_energy_HII_to_HI = avo*13.59843449d0*frac_HII*ev2erg*s% X(k)
121 0 : diss_energy_H2 = avo*4.52d0/2d0*ev2erg*s% X(k)
122 0 : rec_energy_HeII_to_HeI = avo*24.58738880d0*(frac_HeII+frac_HeIII)*ev2erg*s% Y(k)/4d0
123 0 : rec_energy_HeIII_to_HeII = avo*54.4177650d0*frac_HeIII*ev2erg*s% Y(k)/4d0
124 :
125 : adjusted_energy(k) = s% energy(k) &
126 : - (1d0-b% CE_energy_factor_HII_toHI)*rec_energy_HII_to_HI &
127 : - (1d0-b% CE_energy_factor_HeII_toHeI)*rec_energy_HeII_to_HeI &
128 : - (1d0-b% CE_energy_factor_HeIII_toHeII)*rec_energy_HeIII_to_HeII &
129 0 : - (1d0-b% CE_energy_factor_H2)*diss_energy_H2
130 :
131 0 : if (adjusted_energy(k) < 0d0 .or. adjusted_energy(k) > s% energy(k)) then
132 0 : write(*,*) "Error when computing adjusted energy in CE, ", &
133 0 : "s% energy(k):", s% energy(k), " adjusted_energy, ", adjusted_energy(k)
134 0 : stop
135 : end if
136 :
137 0 : if(.false.) then
138 : ! for debug, check the mismatch between the EOS energy and that of a gas+radiation
139 : energy_comp = 3*avo*boltzm*s% T(k)/(2*s% mu(k)) + crad*pow4(s% T(k))/s% rho(k) &
140 : + rec_energy_HII_to_HI &
141 : + rec_energy_HeII_to_HeI &
142 : + rec_energy_HeIII_to_HeII &
143 : + diss_energy_H2
144 :
145 : write(*,*) "compare energies", k, s%m(k)/Msun, s% energy(k), energy_comp, &
146 : (s% energy(k)-energy_comp)/s% energy(k)
147 : end if
148 :
149 : end do
150 :
151 0 : do k=1, s% nz
152 0 : b% CE_m(:) = s% m(:s% nz)
153 : end do
154 :
155 : ! setup values of starting model for the interpolant
156 0 : do k=1, s% nz
157 0 : b% CE_entropy(4*k-3) = exp(s% lnS(k))
158 : end do
159 :
160 : ! Compute internal and potential energies from the inside out, and in the opposite direction.
161 0 : b% CE_U_out(1) = adjusted_energy(1)*s% dm(1)
162 0 : b% CE_Omega_out(1) = - standard_cgrav*s% m(1)*s% dm_bar(1)/s% r(1)
163 0 : do k=2, s% nz
164 0 : b% CE_U_out(4*k-3) = b% CE_U_out(4*(k-1)-3) + adjusted_energy(k)*s% dm(k)
165 0 : b% CE_Omega_out(4*k-3) = b% CE_Omega_out(4*(k-1)-3) - standard_cgrav*s% m(k)*s% dm_bar(k)/s% r(k)
166 : end do
167 0 : b% CE_U_in(4*s% nz-3) = adjusted_energy(s% nz)*s% dm(s% nz)
168 0 : b% CE_Omega_in(4*s% nz-3) = - standard_cgrav*s% m(s% nz)*s% dm_bar(s% nz)/s% r(s% nz)
169 0 : do k=s% nz-1, 1, -1
170 0 : b% CE_U_in(4*k-3) = b% CE_U_in(4*(k+1)-3) + adjusted_energy(k)*s% dm(k)
171 0 : b% CE_Omega_in(4*k-3) = b% CE_Omega_in(4*(k+1)-3) - standard_cgrav*s% m(k)*s% dm_bar(k)/s% r(k)
172 : end do
173 :
174 0 : b% CE_initial_radius = b% r(b% d_i)
175 0 : b% CE_initial_separation = b% separation
176 0 : b% CE_initial_Mdonor = b% m(b% d_i)
177 0 : b% CE_initial_Maccretor = b% m(b% a_i)
178 0 : b% CE_initial_age = s% star_age
179 0 : b% CE_initial_model_number = s% model_number
180 0 : b% CE_b_initial_age = b% binary_age
181 0 : b% CE_b_initial_model_number = b% model_number
182 0 : b% CE_nz = s% nz
183 :
184 0 : allocate(interp_work(s% nz*pm_work_size), stat=ierr)
185 0 : call interp_pm(b% CE_m, s% nz, b% CE_entropy, pm_work_size, interp_work, 'entropy interpolant', op_err)
186 0 : call interp_pm(b% CE_m, s% nz, b% CE_U_in, pm_work_size, interp_work, 'U_in interpolant', op_err)
187 0 : call interp_pm(b% CE_m, s% nz, b% CE_U_out, pm_work_size, interp_work, 'U_out interpolant', op_err)
188 0 : call interp_pm(b% CE_m, s% nz, b% CE_Omega_in, pm_work_size, interp_work, 'Omega_in interpolant', op_err)
189 0 : call interp_pm(b% CE_m, s% nz, b% CE_Omega_out, pm_work_size, interp_work, 'Omega_out interpolant', op_err)
190 0 : if(op_err /= 0) then
191 0 : ierr = -1
192 0 : write(*,*) "CE_init: Error while creating interpolants"
193 0 : return
194 : end if
195 0 : deallocate(adjusted_energy,interp_work)
196 : end if
197 :
198 0 : end subroutine CE_init
199 :
200 0 : subroutine CE_rlo_mdot(binary_id, rlo_mdot, ierr)
201 : integer, intent(in) :: binary_id
202 : real(dp), intent(out) :: rlo_mdot
203 : integer, intent(out) :: ierr
204 : type (binary_info), pointer :: b
205 0 : real(dp) :: exp_factor
206 :
207 : ierr = 0
208 :
209 0 : call binary_ptr(binary_id, b, ierr)
210 :
211 0 : if (b% use_other_CE_rlo_mdot) then
212 0 : call b% other_CE_rlo_mdot(b% binary_id, rlo_mdot, ierr)
213 0 : return
214 : end if
215 :
216 0 : exp_factor = -log(b% CE_mass_loss_rate_low/b% CE_mass_loss_rate_high)
217 :
218 0 : if (b% r(b% d_i)-b% rl(b% d_i) > 0d0) then
219 0 : rlo_mdot = -Msun/secyer*b% CE_mass_loss_rate_high
220 0 : else if (b% r(b% d_i)-b% rl(b% d_i) < -b% CE_rel_rlo_for_detachment*b% rl(b% d_i)) then
221 0 : rlo_mdot = -Msun/secyer*b% CE_mass_loss_rate_low
222 : else
223 : rlo_mdot = -Msun/secyer*b% CE_mass_loss_rate_high * &
224 0 : exp(exp_factor*(b% r(b% d_i)-b% rl(b% d_i))/(b% rl(b% d_i)*b% CE_rel_rlo_for_detachment))
225 : end if
226 :
227 0 : end subroutine CE_rlo_mdot
228 :
229 0 : integer function CE_binary_evolve_step(b)
230 : use binary_utils, only:set_separation_eccentricity
231 : type (binary_info), pointer :: b
232 : type (star_info), pointer :: s
233 0 : real(dp) :: Ebind, Ecore, Ecore_i, lambda, &
234 0 : U_removed ,Omega_removed, U_inold, Omega_inold
235 0 : real(dp) :: separation, initial_Eorb
236 : integer :: ierr, op_err, k
237 :
238 0 : if (b% use_other_CE_binary_evolve_step) then
239 0 : CE_binary_evolve_step = b% other_CE_binary_evolve_step(b% binary_id)
240 0 : return
241 : end if
242 :
243 : ! setup mdot_system_wind to get output right, it is not actually used.
244 : ! setup jdots to zero as well
245 0 : b% mdot_system_wind(b% d_i) = b% s_donor% mstar_dot - b% mtransfer_rate
246 0 : if (b% point_mass_i == 0) then
247 0 : b% mdot_system_wind(b% a_i) = b% s_accretor% mstar_dot
248 : else
249 0 : b% mdot_system_wind(b% a_i) = 0.0d0
250 : end if
251 0 : b% jdot = 0d0
252 0 : b% jdot_mb = 0d0
253 0 : b% jdot_gr = 0d0
254 0 : b% jdot_ml = 0d0
255 0 : b% jdot_missing_wind = 0d0
256 0 : b% extra_jdot = 0d0
257 0 : b% jdot_ls = 0d0
258 :
259 0 : s => b% s_donor
260 :
261 0 : if (b% CE_fixed_lambda < 0d0) then
262 0 : Ecore = 0
263 0 : do k=1, s% nz
264 : Ecore = Ecore + s% energy(k)*s% dm(k) - standard_cgrav*s% m(k)*s% dm_bar(k)/s% r(k)
265 : end do
266 :
267 0 : call interp_value(b% CE_m, b% CE_nz, b% CE_U_out, s% m(1), U_removed, op_err)
268 0 : call interp_value(b% CE_m, b% CE_nz, b% CE_Omega_out, s% m(1), Omega_removed, op_err)
269 0 : call interp_value(b% CE_m, b% CE_nz, b% CE_U_in, s% m(1), U_inold, op_err)
270 0 : call interp_value(b% CE_m, b% CE_nz, b% CE_Omega_in, s% m(1), Omega_inold, op_err)
271 :
272 0 : Ecore_i = U_inold + Omega_inold
273 0 : Ebind = b% CE_alpha_th*U_removed + Omega_removed
274 : lambda = -(standard_cgrav*b% CE_initial_Mdonor*(b% CE_initial_Mdonor - s% m(1))) &
275 0 : /(Ebind*b% CE_initial_radius)
276 : else
277 0 : lambda = b% CE_fixed_lambda
278 : Ebind = -(standard_cgrav*b% CE_initial_Mdonor*(b% CE_initial_Mdonor - s% m(1))) &
279 0 : /(lambda*b% CE_initial_radius)
280 : end if
281 :
282 0 : if (b% d_i == 1) then
283 0 : b% CE_Ebind1 = Ebind
284 0 : b% CE_lambda1 = lambda
285 : else
286 0 : b% CE_Ebind2 = Ebind
287 0 : b% CE_lambda2 = lambda
288 : end if
289 :
290 0 : initial_Eorb = -standard_cgrav*b% CE_initial_Mdonor*b% CE_initial_Maccretor/(2*b% CE_initial_separation)
291 :
292 : separation = -b% CE_alpha*standard_cgrav*s% m(1)*b% CE_initial_Maccretor &
293 0 : /(2*(Ebind+b% CE_alpha*initial_Eorb))
294 :
295 0 : b% m(b% d_i) = b% s_donor% mstar
296 0 : b% time_step = b% s_donor% time_step
297 0 : if (b% point_mass_i == 0) then
298 0 : b% m(b% a_i) = b% s_accretor% mstar
299 : end if
300 :
301 0 : if (b% point_mass_i /= 1) then
302 0 : b% r(1) = Rsun*b% s1% photosphere_r
303 : else
304 0 : b% r(1) = 0
305 : end if
306 0 : if (b% point_mass_i /= 2) then
307 0 : b% r(2) = Rsun*b% s2% photosphere_r
308 : else
309 0 : b% r(2) = 0
310 : end if
311 :
312 : !only change separation if its reduced from the initial value
313 : call set_separation_eccentricity(b% binary_id, &
314 0 : min(b% CE_initial_separation, separation), b% eccentricity, ierr)
315 :
316 0 : b% model_number = b% model_number + 1
317 0 : b% time_step = b% s_donor% time_step
318 0 : b% binary_age = b% binary_age + b% time_step
319 :
320 0 : if (b% r(b% d_i)-b% rl(b% d_i) < 0d0) then
321 0 : b% CE_years_detached = b% CE_years_detached + b% time_step
322 : else
323 0 : b% CE_years_detached = 0d0
324 : end if
325 :
326 0 : CE_binary_evolve_step = keep_going
327 :
328 0 : end function CE_binary_evolve_step
329 :
330 0 : integer function CE_binary_finish_step(b)
331 0 : use binary_utils, only: eval_rlobe
332 : type (binary_info), pointer :: b
333 0 : real(dp) :: h_diff, he_diff, rlobe
334 : logical :: terminate_CE
335 : integer :: k
336 0 : CE_binary_finish_step = keep_going
337 :
338 0 : if (b% use_other_CE_binary_finish_step) then
339 0 : CE_binary_finish_step = b% other_CE_binary_finish_step(b% binary_id)
340 0 : return
341 : end if
342 :
343 0 : terminate_CE = .false.
344 0 : if ((b% r(b% d_i)-b% rl(b% d_i))/(b% rl(b% d_i)*b% CE_rel_rlo_for_detachment) < -1d0) then
345 0 : terminate_CE = .true.
346 0 : write(*,*) "Have reached CE_rel_rlo_for_detachment"
347 0 : else if (b% CE_years_detached > b% CE_years_detached_to_terminate) then
348 0 : terminate_CE = .true.
349 0 : write(*,*) "Have reached CE_years_detached_to_terminate"
350 : end if
351 :
352 : if (terminate_CE) then
353 0 : b% CE_flag = .false.
354 0 : b% mtransfer_rate = 0d0
355 0 : b% s_donor% mix_factor = 1d0
356 0 : b% s_donor% dxdt_nuc_factor = 1d0
357 0 : b% s_donor% timestep_hold = b% s_donor% model_number
358 0 : if (b% point_mass_i == 0) then
359 0 : b% s_accretor% mix_factor = 1d0
360 0 : b% s_accretor% dxdt_nuc_factor = 1d0
361 0 : b% s_accretor% timestep_hold = b% s_accretor% model_number
362 : end if
363 :
364 0 : b% keep_donor_fixed = .false.
365 0 : write(*,*) "TURNING OFF CE"
366 : end if
367 :
368 : ! termination conditions
369 :
370 0 : h_diff = abs(b% s_donor% center_h1 - b% s_donor% surface_h1)
371 0 : he_diff = abs(b% s_donor% center_he4 - b% s_donor% surface_he4)
372 : if (h_diff < b% CE_xa_diff_to_terminate &
373 0 : .and. he_diff < b% CE_xa_diff_to_terminate) then
374 0 : write(*,*) "Central and surface abundances below CE_xa_diff_to_terminate"
375 0 : write(*,*) "Terminating evolution"
376 0 : CE_binary_finish_step = terminate
377 0 : return
378 : end if
379 :
380 : !terminate if, for the current orbital separation, stripping the star down to the point
381 : !where CE_xa_diff_to_terminate would apply, the remaining layers would be Roche lobe
382 : !overflowing.
383 0 : if (b% CE_terminate_when_core_overflows) then
384 0 : do k = 1, b% s_donor% nz
385 0 : h_diff = abs(b% s_donor% center_h1 - b% s_donor% X(k))
386 0 : he_diff = abs(b% s_donor% center_he4 - b% s_donor% Y(k))
387 : if (h_diff < b% CE_xa_diff_to_terminate &
388 0 : .and. he_diff < b% CE_xa_diff_to_terminate) then
389 0 : rlobe = eval_rlobe(b% s_donor% m(k), b% m(b% a_i), b% separation)
390 0 : if (b% s_donor% r(k) > rlobe) then
391 0 : write(*,*) "Terminate due to CE_terminate_when_core_overflows"
392 0 : write(*,*) "Terminating evolution"
393 0 : CE_binary_finish_step = terminate
394 0 : return
395 : end if
396 : exit
397 : end if
398 : end do
399 : end if
400 :
401 0 : if (b% period < b% CE_min_period_in_minutes*60d0) then
402 0 : write(*,*) "Orbital period is below CE_min_period_in_minutes"
403 0 : write(*,*) "Terminating evolution"
404 0 : CE_binary_finish_step = terminate
405 0 : return
406 : end if
407 :
408 0 : end function CE_binary_finish_step
409 :
410 : end module binary_ce
411 :
|
