Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2021 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 phase_separation
21 :
22 : use star_private_def
23 : use const_def, only: dp
24 :
25 : implicit none
26 :
27 : private
28 : public :: do_phase_separation
29 :
30 : logical, parameter :: dbg = .false.
31 :
32 : ! offset to higher phase than 0.5 to avoid interference
33 : ! between phase separation mixing and latent heat for Skye.
34 : real(dp), parameter :: eos_phase_boundary = 0.9d0
35 :
36 : contains
37 :
38 0 : subroutine do_phase_separation(s, dt, ierr)
39 : type (star_info), pointer :: s
40 : real(dp), intent(in) :: dt
41 : integer, intent(out) :: ierr
42 :
43 0 : if(s% phase_separation_option == 'CO') then
44 0 : call do_2component_phase_separation(s, dt, 'CO', ierr)
45 0 : else if(s% phase_separation_option == 'ONe') then
46 0 : call do_2component_phase_separation(s, dt, 'ONe', ierr)
47 : else
48 0 : write(*,*) 'invalid phase_separation_option'
49 0 : stop
50 : end if
51 0 : end subroutine do_phase_separation
52 :
53 0 : subroutine do_2component_phase_separation(s, dt, components, ierr)
54 : use chem_def, only: ic12, io16, ine20
55 : use chem_lib, only: chem_get_iso_id
56 : type (star_info), pointer :: s
57 : real(dp), intent(in) :: dt
58 : character (len=*), intent(in) :: components
59 : integer, intent(out) :: ierr
60 :
61 0 : real(dp) :: XNe, XO, XC, pad
62 : integer :: k, k_bound, kstart, net_ic12, net_io16, net_ine20
63 : logical :: save_Skye_use_ion_offsets
64 :
65 : ! Set phase separation mixing mass negative at beginning of phase separation
66 0 : s% phase_sep_mixing_mass = -1d0
67 0 : s% eps_phase_separation(1:s%nz) = 0d0
68 :
69 0 : if(s% phase(s% nz) < eos_phase_boundary) then
70 0 : s% crystal_core_boundary_mass = 0d0
71 0 : return
72 : end if
73 :
74 0 : net_ic12 = s% net_iso(ic12)
75 0 : net_io16 = s% net_iso(io16)
76 0 : net_ine20 = s% net_iso(ine20)
77 :
78 : ! Find zone of phase transition from liquid to solid
79 0 : k_bound = -1
80 0 : do k = s%nz,1,-1
81 0 : if(s% phase(k-1) <= eos_phase_boundary .and. s% phase(k) > eos_phase_boundary) then
82 : k_bound = k
83 : exit
84 : end if
85 : end do
86 :
87 0 : XC = s% xa(net_ic12,k_bound)
88 0 : XO = s% xa(net_io16,k_bound)
89 0 : XNe = s% xa(net_ine20,k_bound)
90 : ! Check that we're still in C/O or O/Ne dominated material as appropriate,
91 : ! otherwise skip phase separation
92 0 : if(components == 'CO'.and. XO + XC < 0.9d0) return
93 0 : if(components == 'ONe'.and. XNe + XO < 0.8d0) return ! O/Ne mixtures tend to have more byproducts of burning mixed in
94 :
95 : ! If there is a phase transition, reset the composition at the boundary
96 0 : if(k_bound > 0) then
97 :
98 : ! core boundary needs to be padded by a minimal amount (less than a zone worth of mass)
99 : ! to account for loss of precision during remeshing.
100 0 : pad = s% min_dq * s% m(1) * 0.5d0
101 0 : do k = s%nz,1,-1
102 0 : if(s% m(k) > s% crystal_core_boundary_mass + pad) then
103 : kstart = k
104 : exit
105 : end if
106 : end do
107 :
108 : ! calculate energy associated with phase separation, ignoring the ionization
109 : ! energy term that Skye sometimes calculates
110 0 : save_Skye_use_ion_offsets = s% eos_rq% Skye_use_ion_offsets
111 0 : s% eos_rq% Skye_use_ion_offsets = .false.
112 0 : call update_model_(s,1,s%nz,.false.)
113 0 : do k=1,s% nz
114 0 : s% eps_phase_separation(k) = s% energy(k)
115 : end do
116 :
117 : ! loop runs outward starting at previous crystallization boundary
118 0 : do k = kstart,1,-1
119 : ! Start by checking if this material should be crystallizing
120 0 : if(s% phase(k) <= eos_phase_boundary) then
121 0 : s% crystal_core_boundary_mass = s% m(k+1)
122 0 : exit
123 : end if
124 :
125 0 : call move_one_zone(s,k,components)
126 : ! crystallized out to k now, liquid starts at k-1.
127 : ! now mix the liquid material outward until stably stratified
128 0 : call mix_outward(s, k-1, 0)
129 :
130 : end do
131 :
132 0 : call update_model_(s,1,s%nz,.false.)
133 :
134 : ! phase separation heating term for use by energy equation
135 0 : do k=1,s% nz
136 0 : s% eps_phase_separation(k) = (s% eps_phase_separation(k) - s% energy(k)) / dt
137 : end do
138 0 : s% eos_rq% Skye_use_ion_offsets = save_Skye_use_ion_offsets
139 0 : s% need_to_setvars = .true.
140 : end if
141 :
142 0 : ierr = 0
143 0 : end subroutine do_2component_phase_separation
144 :
145 0 : subroutine move_one_zone(s,k,components)
146 0 : use chem_def, only: ic12, io16, ine20
147 : use chem_lib, only: chem_get_iso_id
148 : type(star_info), pointer :: s
149 : integer, intent(in) :: k
150 : character (len=*), intent(in) :: components
151 :
152 0 : real(dp) :: XC, XO, XNe, XC1, XO1, XNe1, dXO, dXNe, Xfac
153 : integer :: net_ic12, net_io16, net_ine20
154 :
155 0 : net_ic12 = s% net_iso(ic12)
156 0 : net_io16 = s% net_iso(io16)
157 0 : net_ine20 = s% net_iso(ine20)
158 :
159 0 : if(components == 'CO') then
160 0 : XO = s% xa(net_io16,k)
161 0 : XC = s% xa(net_ic12,k)
162 :
163 : ! Call Blouin phase diagram.
164 : ! Need to rescale temporarily because phase diagram assumes XO + XC = 1
165 0 : Xfac = XO + XC
166 0 : XO = XO/Xfac
167 0 : XC = XC/Xfac
168 :
169 0 : dXO = blouin_delta_xo(XO)
170 :
171 0 : s% xa(net_io16,k) = Xfac*(XO + dXO)
172 0 : s% xa(net_ic12,k) = Xfac*(XC - dXO)
173 :
174 : ! Redistribute change in C,O into zone k-1,
175 : ! conserving total mass of C,O
176 0 : XC1 = s% xa(net_ic12,k-1)
177 0 : XO1 = s% xa(net_io16,k-1)
178 0 : s% xa(net_ic12,k-1) = XC1 + Xfac*dXO * s% dq(k) / s% dq(k-1)
179 0 : s% xa(net_io16,k-1) = XO1 - Xfac*dXO * s% dq(k) / s% dq(k-1)
180 0 : else if(components == 'ONe') then
181 0 : XNe = s% xa(net_ine20,k)
182 0 : XO = s% xa(net_io16,k)
183 :
184 : ! Call Blouin phase diagram.
185 : ! Need to rescale temporarily because phase diagram assumes XO + XNe = 1
186 0 : Xfac = XO + XNe
187 0 : XO = XO/Xfac
188 0 : XNe = XNe/Xfac
189 :
190 0 : dXNe = blouin_delta_xne(XNe)
191 :
192 0 : s% xa(net_ine20,k) = Xfac*(XNe + dXNe)
193 0 : s% xa(net_io16,k) = Xfac*(XO - dXNe)
194 :
195 : ! Redistribute change in Ne,O into zone k-1,
196 : ! conserving total mass of Ne,O
197 0 : XO1 = s% xa(net_io16,k-1)
198 0 : XNe1 = s% xa(net_ine20,k-1)
199 0 : s% xa(net_io16,k-1) = XO1 + Xfac*dXNe * s% dq(k) / s% dq(k-1)
200 0 : s% xa(net_ine20,k-1) = XNe1 - Xfac*dXNe * s% dq(k) / s% dq(k-1)
201 : else
202 0 : write(*,*) 'invalid components option in phase separation'
203 0 : stop
204 : end if
205 :
206 0 : call update_model_(s,k-1,s%nz,.true.)
207 :
208 0 : end subroutine move_one_zone
209 :
210 : ! mix composition outward until reaching stable composition profile
211 0 : subroutine mix_outward(s,kbot,min_mix_zones)
212 : type(star_info), pointer :: s
213 : integer, intent(in) :: kbot, min_mix_zones
214 :
215 0 : real(dp) :: avg_xa(s%species)
216 0 : real(dp) :: mass, B_term, grada, gradr
217 : integer :: k, l, ktop
218 : logical :: use_brunt
219 :
220 0 : use_brunt = s% phase_separation_mixing_use_brunt
221 :
222 0 : do k=kbot-min_mix_zones,1,-1
223 0 : ktop = k
224 :
225 0 : if (s% m(ktop) > s% phase_sep_mixing_mass) then
226 0 : s% phase_sep_mixing_mass = s% m(ktop)
227 : end if
228 :
229 0 : mass = SUM(s%dm(ktop:kbot))
230 0 : do l = 1, s%species
231 0 : avg_xa(l) = SUM(s%dm(ktop:kbot)*s%xa(l,ktop:kbot))/mass
232 : end do
233 :
234 : ! some potential safeguards from conv_premix
235 : ! avg_xa = MAX(MIN(avg_xa, 1._dp), 0._dp)
236 : ! avg_xa = avg_xa/SUM(avg_xa)
237 :
238 0 : do l = 1, s%species
239 0 : s%xa(l,ktop:kbot) = avg_xa(l)
240 : end do
241 :
242 : ! updates, eos, opacities, mu, etc now that abundances have changed,
243 : ! but only in the cells near the boundary where we need to check here.
244 : ! Will call full update over mixed region after exiting loop.
245 0 : call update_model_(s, ktop-1, ktop+1, use_brunt)
246 :
247 0 : if(use_brunt) then
248 0 : B_term = s% unsmoothed_brunt_B(ktop)
249 0 : grada = s% grada_face(ktop)
250 0 : gradr = s% gradr(ktop)
251 0 : if(B_term + grada - gradr > 0d0) then
252 : ! stable against further mixing, so exit loop
253 : exit
254 : end if
255 : else ! simpler calculation based on mu gradient
256 0 : if(s% mu(ktop) >= s% mu(ktop-1)) then
257 : ! stable against further mixing, so exit loop
258 : exit
259 : end if
260 : end if
261 :
262 : end do
263 :
264 : ! Call a final update over all mixed cells now.
265 0 : call update_model_(s, ktop, kbot+1, .true.)
266 :
267 0 : end subroutine mix_outward
268 :
269 0 : real(dp) function blouin_delta_xo(Xin)
270 : real(dp), intent(in) :: Xin ! mass fraction
271 0 : real(dp) :: Xnew ! mass fraction
272 0 : real(dp) :: xo, dxo ! number fractions
273 0 : real(dp) :: a0, a1, a2, a3, a4, a5
274 :
275 : ! Convert input mass fraction to number fraction, assuming C/O mixture
276 0 : xo = (Xin/16d0)/(Xin/16d0 + (1d0 - Xin)/12d0)
277 :
278 0 : a0 = 0d0
279 0 : a1 = -0.311540d0
280 0 : a2 = 2.114743d0
281 0 : a3 = -1.661095d0
282 0 : a4 = -1.406005d0
283 0 : a5 = 1.263897d0
284 :
285 : dxo = &
286 : a0 + &
287 : a1*xo + &
288 : a2*xo*xo + &
289 : a3*xo*xo*xo + &
290 : a4*xo*xo*xo*xo + &
291 0 : a5*xo*xo*xo*xo*xo
292 :
293 0 : xo = xo + dxo
294 :
295 : ! Convert back to mass fraction
296 0 : Xnew = 16d0*xo/(16d0*xo + 12d0*(1d0-xo))
297 :
298 0 : blouin_delta_xo = Xnew - Xin
299 0 : end function blouin_delta_xo
300 :
301 0 : real(dp) function blouin_delta_xne(Xin)
302 : real(dp), intent(in) :: Xin ! mass fraction
303 0 : real(dp) :: Xnew ! mass fraction
304 0 : real(dp) :: xne, dxne ! number fractions
305 0 : real(dp) :: a0, a1, a2, a3, a4, a5
306 :
307 : ! Convert input mass fraction to number fraction, assuming O/Ne mixture
308 0 : xne = (Xin/20d0)/(Xin/20d0 + (1d0 - Xin)/16d0)
309 :
310 0 : a0 = 0d0
311 0 : a1 = -0.120299d0
312 0 : a2 = 1.304399d0
313 0 : a3 = -1.722625d0
314 0 : a4 = 0.393996d0
315 0 : a5 = 0.144529d0
316 :
317 : dxne = &
318 : a0 + &
319 : a1*xne + &
320 : a2*xne*xne + &
321 : a3*xne*xne*xne + &
322 : a4*xne*xne*xne*xne + &
323 0 : a5*xne*xne*xne*xne*xne
324 :
325 0 : xne = xne + dxne
326 :
327 : ! Convert back to mass fraction
328 0 : Xnew = 20d0*xne/(20d0*xne + 16d0*(1d0-xne))
329 :
330 0 : blouin_delta_xne = Xnew - Xin
331 0 : end function blouin_delta_xne
332 :
333 0 : subroutine update_model_ (s, kc_t, kc_b, do_brunt)
334 :
335 : use turb_info, only: set_mlt_vars
336 : use brunt, only: do_brunt_B
337 : use micro
338 :
339 : type(star_info), pointer :: s
340 : integer, intent(in) :: kc_t
341 : integer, intent(in) :: kc_b
342 : logical, intent(in) :: do_brunt
343 :
344 : integer :: ierr
345 : integer :: kf_t
346 : integer :: kf_b
347 :
348 0 : logical :: mask(s%nz)
349 :
350 0 : mask(:) = .true.
351 :
352 : ! Update the model to reflect changes in the abundances across
353 : ! cells kc_t:kc_b (the mask part of this call is unused, mask=true for all zones).
354 : ! Do updates at constant (P,T) rather than constant (rho,T).
355 0 : s%fix_Pgas = .true.
356 0 : call set_eos_with_mask(s, kc_t, kc_b, mask, ierr)
357 0 : if (ierr /= 0) then
358 0 : write(*,*) 'phase_separation: error from call to set_eos_with_mask'
359 0 : stop
360 : end if
361 0 : s%fix_Pgas = .false.
362 :
363 : ! Update opacities across cells kc_t:kc_b (this also sets rho_face
364 : ! and related quantities on faces kc_t:kc_b)
365 : call set_micro_vars(s, kc_t, kc_b, &
366 0 : skip_eos=.TRUE., skip_net=.TRUE., skip_neu=.TRUE., skip_kap=.FALSE., ierr=ierr)
367 0 : if (ierr /= 0) then
368 0 : write(*,*) 'phase_separation: error from call to set_micro_vars'
369 0 : stop
370 : end if
371 :
372 : ! This is expensive, so only do it if we really need to.
373 0 : if(do_brunt) then
374 : ! Need to make sure we can set brunt for mix_outward calculation.
375 0 : if(.not. s% calculate_Brunt_B) then
376 0 : stop "phase separation requires s% calculate_Brunt_B = .true."
377 : end if
378 0 : call do_brunt_B(s, kc_t, kc_b, ierr) ! for unsmoothed_brunt_B
379 0 : if (ierr /= 0) then
380 0 : write(*,*) 'phase_separation: error from call to do_brunt_B'
381 0 : stop
382 : end if
383 : end if
384 :
385 : ! Finally update MLT for interior faces
386 :
387 0 : kf_t = kc_t
388 0 : kf_b = kc_b + 1
389 :
390 0 : call set_mlt_vars(s, kf_t+1, kf_b-1, ierr)
391 0 : if (ierr /= 0) then
392 0 : write(*,*) 'phase_separation: failed in call to set_mlt_vars during update_model_'
393 0 : stop
394 : end if
395 :
396 0 : return
397 :
398 0 : end subroutine update_model_
399 :
400 : end module phase_separation
|
