Three component phase separation for ultramassive white dwarfs, part 2#863
Three component phase separation for ultramassive white dwarfs, part 2#863
Conversation
Updated code in phase_separation.f90 and a new wd_o_ne_3_phase test suite.
|
@mcastrotapia I've pushed a few "housekeeping" fixes to get the static analysis script to pass, so you might want to merge this branch into your fork to stay up to date. Anyway, this PR will eventually be the one that merges the code into Thanks for submitting this code! |
|
Thanks for the fixes and for reviewing the PR! I will merge the branch into my fork, and I'll be attentive to comments/questions and updates. |
| if(s% phase(k) <= eos_phase_boundary) then | ||
| s% crystal_core_boundary_mass = s% m(k+1) | ||
| exit | ||
| if (s% crystal_core_boundary_mass>s% m(k+1)) then |
There was a problem hiding this comment.
What is this meant to catch? It seems there must be some subtle bug here I was missing.
There was a problem hiding this comment.
I think I was trying to prevent the core from moving inwards in case the extra heating suddenly changes the energy in both phases and the phase parameter. However, I think it's not necessary since the lines where I added the following:
if (s% phase(s% nz) < eos_phase_boundary) then !!! prevent to move the core size inwards if the core is suddently "melted" leaving everything liquid under phi<0.9
`if (s% crystal_core_boundary_mass>0d0)then
s% crystal_core_boundary_mass=s% crystal_core_boundary_mass
return
else
s% crystal_core_boundary_mass = 0d0
return
end if
end if`
must be enough for that
| end if | ||
| end subroutine do_phase_separation | ||
|
|
||
| subroutine do_2component_phase_separation(s, dt, components, ierr) |
There was a problem hiding this comment.
If we keep the current code structure, we'll need to rename this subroutine to something other than "do_2component...", since it is now explicitly being used to do some 3 component cases as well. I had another branch a while ago that was attempting a different branch for 3-component distillation, but that is now somewhat stale and never got to a point of being ready to merge in. I'll need to reconcile the naming scheme somehow...
There was a problem hiding this comment.
Okay, I agree, if you want me to help with that, please let me know :)
| if(s% phase(s% nz) < eos_phase_boundary) then !!! prevent to move the core size inwards if the core is suddently "melted" leaving everything liquid under phi<0.9 | ||
| if (s% crystal_core_boundary_mass>0d0)then | ||
| s% crystal_core_boundary_mass=s% crystal_core_boundary_mass | ||
| return |
There was a problem hiding this comment.
I see why you might need something like this, but I'm not sure this is the best approach in general. One alternative I've experimented with is just redistributing the phase separation heating over a significant fraction of the interior. That way you avoid transient local heating spikes that might melt the core, but still get the net right amount of energy deposited in the interior, which seems to work for capturing the effect on cooling rate. @mcastrotapia do you think that could be a viable alternative here?
There was a problem hiding this comment.
I think it is a good alternative. At some point, I tried something similar, but I didn't have success, and I ended up doing this because it seemed easier and faster to implement.
| ! Check that we're still in C/O or O/Ne dominated material as appropriate, | ||
| ! otherwise skip phase separation | ||
| if(components == 'CO'.and. XO + XC < 0.9d0) return | ||
| if(components == 'ONe'.and. XNe + XO < 0.8d0) return ! O/Ne mixtures tend to have more byproducts of burning mixed in |
There was a problem hiding this comment.
I believe we need to bring these lines back. Is there an argument for deleting them?
They shouldn't affect the new 3-componenent functionality that we're trying to add here, and I'd like to keep the behavior of the 2-component options unchanged unless there is a compelling reason to change them.
There was a problem hiding this comment.
I deleted them when trying to let the code decide for itself if doing C/O or O/Ne in the case of 2-component phase separation, this is because in ultramassive WD models the core is usually O/Ne, but towards the outer layers there is a C/O dominated zone before reaching the H/He atmosphere (for example, Figure 2 of Camisassa et al., 2019). In the 3-component, I was also letting the code check which were the 3 most abundant elements in the liquid zone to decide which phase diagram must be used.
| call set_micro_vars(s, kc_t, kc_b, & | ||
| skip_eos=.TRUE., skip_net=.TRUE., skip_neu=.TRUE., skip_kap=.FALSE., ierr=ierr) | ||
| skip_eos=.TRUE., skip_net=.TRUE., skip_neu=.TRUE., skip_kap=.TRUE., ierr=ierr) |
There was a problem hiding this comment.
I think we need to revert this change. This is now skipping the opacity update that is the point of this call. In fact, as written, I think this call to set_micro_vars skips everything and does nothing at all.
There was a problem hiding this comment.
I agree, I was mostly having convergence problems with this when using a small timestep and (I think) with the heating in the solid phase for such timesteps. That's why it would be good to implement the heating term as you suggested in one of the previous comments
| do j=1,num_x1 | ||
| do i=1,num_x2 | ||
| read(iounit,*) x1l, x2l, deltax1 | ||
| x1_l(j)=x1l | ||
| if (j == 1) then | ||
| x2_l(i) =x2l | ||
| end if | ||
| deltax1_sob_f(1,j,i) = deltax1 | ||
| end do | ||
| end do | ||
| close(iounit) | ||
| ! just use "not a knot" bc's at edges of tables | ||
| ibcxmin = 0; bcxmin(1:num_x1) = 0 | ||
| ibcxmax = 0; bcxmax(1:num_x1) = 0 | ||
| ibcymin = 0; bcymin(1:num_x2) = 0 | ||
| ibcymax = 0; bcymax(1:num_x2) = 0 | ||
| call interp_mkbicub_db( & | ||
| x1_l, num_x1, x2_l, num_x2, deltax1_sob_f1, num_x1, & | ||
| ibcxmin,bcxmin,ibcxmax,bcxmax, & | ||
| ibcymin,bcymin,ibcymax,bcymax, & | ||
| ilinx,iliny,ierr) | ||
| if (ierr /= 0) then | ||
| write(*,*) 'interp_mkbicub_db error' | ||
| ierr = -1 | ||
| call mesa_error(__FILE__,__LINE__) | ||
| end if | ||
| do j=1,num_x1 | ||
| do i=1,num_x2 | ||
| do k=1,4 | ||
| if (is_bad(deltax1_sob_f(k,j,i))) then | ||
| write(*,*) 'deltax1_sob_f', i, j, k, deltax1_sob_f(k,j,i) | ||
| end if | ||
| end do | ||
| end do | ||
| end do | ||
| call interp_evbicub_db( & | ||
| x1_, x2_, x1_l, num_x1, x2_l, num_x2, & | ||
| ilinx, iliny, deltax1_sob_f1, num_x1, ict, fval, ier) | ||
| dx1_=fval(1) ! delta_x1 from 2d interpolation |
There was a problem hiding this comment.
Indentation of this block should be fixed. (Just a minor aesthetic note.)
There was a problem hiding this comment.
There's another block later on with the same indentation problem as well.
There was a problem hiding this comment.
Note to self for later: update run_star_extras, rn script, etc to make this a fully-fledged test case and integrate it into do1_test_source for the test suite.
| ierr = 0 | ||
| end subroutine do_2component_phase_separation | ||
|
|
||
| subroutine move_one_zone(s,k,components) |
There was a problem hiding this comment.
The move_one_zone subroutine has developed a lot of different branches that I think make it rather confusing. I can see why it developed into what it is, but I'd probably like to refactor this into two separate routines: move_one_zone_for_2comp that basically just retains the old functionality (plus the nice addition of accounting for Ne22), and move_one_zone_for_3comp that implements all the new 3-component cases.
|
We'll also need to add some documentation of the new options in controls.defaults. |
3 participants
Updated code in phase_separation.f90 and a new wd_o_ne_3_phase test suite.