WD CO Phase Separation Updates (#480)

* simplify phase separation EOS update interface

* do phase separation EOS update at constant pressure

* update Skye interface to make computing ion offsets optional

* ignore ion offset energy when computing phase separation heating

* simplify mixing procedure now that energetics are better

Author: evbauer

eos/defaults/eos.defaults

@@ -180,6 +180,14 @@
 
    mass_fraction_limit_for_Skye = 1d-4
 
+   ! Skye_use_ion_offsets
+   ! ~~~~~~~~~~~~~~~~~~~~
+   ! If true, set neutral ground state as zero point for each ion species energy.
+   ! If false, zero point assumes fully ionized.
+   ! ::
+
+   Skye_use_ion_offsets = .true.
+
    ! Skye Gamma Limits
    ! ~~~~~~~~~~~~~~~~~
    ! The one-component plasma fits in Skye are only valid in certain ranges of the ion interaction parameter Gamma.

eos/private/eos_ctrls_io.f90

@@ -106,6 +106,7 @@ module eos_ctrls_io
 
    ! limits for Skye
    logical :: use_Skye
+   logical :: Skye_use_ion_offsets
    real(dp) :: mass_fraction_limit_for_Skye
    real(dp) :: Skye_min_gamma_for_solid ! The minimum Gamma_i at which to use the solid free energy fit (below this, extrapolate).
    real(dp) :: Skye_max_gamma_for_liquid ! The maximum Gamma_i at which to use the liquid free energy fit (above this, extrapolate).
@@ -231,6 +232,7 @@ module eos_ctrls_io
 
       ! controls for Skye
       use_Skye, &
+      Skye_use_ion_offsets, &
       mass_fraction_limit_for_Skye, &
       Skye_min_gamma_for_solid, & ! The minimum Gamma_i at which to use the solid free energy fit (below this, extrapolate).
       Skye_max_gamma_for_liquid, & ! The maximum Gamma_i at which to use the liquid free energy fit (above this, extrapolate).
@@ -486,6 +488,7 @@ subroutine store_controls(rq)
       rq% PC_Gamma_full_crystal = PC_Gamma_full_crystal
       ! controls for Skye
       rq% use_Skye = use_Skye
+      rq% Skye_use_ion_offsets = Skye_use_ion_offsets
       rq% mass_fraction_limit_for_Skye = mass_fraction_limit_for_Skye
       rq%Skye_min_gamma_for_solid = Skye_min_gamma_for_solid
       rq%Skye_max_gamma_for_liquid = Skye_max_gamma_for_liquid
@@ -622,6 +625,7 @@ subroutine set_controls_for_writing(rq)
       PC_Gamma_full_crystal = rq% PC_Gamma_full_crystal
       ! controls for Skye
       use_Skye = rq% use_Skye
+      Skye_use_ion_offsets = rq% Skye_use_ion_offsets
       mass_fraction_limit_for_Skye = rq% mass_fraction_limit_for_Skye   
       Skye_min_gamma_for_solid = rq% Skye_min_gamma_for_solid
       Skye_max_gamma_for_liquid = rq% Skye_max_gamma_for_liquid  

eos/private/skye.f90

@@ -211,6 +211,7 @@ subroutine Get_Skye_EOS_Results( &
             T, Rho, X, abar, zbar, &
             rq%Skye_min_gamma_for_solid, rq%Skye_max_gamma_for_liquid, &
             rq%Skye_solid_mixing_rule, rq%mass_fraction_limit_for_Skye, &
+            rq%Skye_use_ion_offsets, &
             species, chem_id, xa, &
             res, d_dlnd, d_dlnT, d_dxa, ierr)
 
@@ -259,7 +260,8 @@ subroutine skye_eos( &
             temp_in, den_in, Xfrac, abar, zbar,  &
             Skye_min_gamma_for_solid, Skye_max_gamma_for_liquid, &
             Skye_solid_mixing_rule, &
-            mass_fraction_limit, species, chem_id, xa, &
+            mass_fraction_limit, use_ion_offsets, &
+            species, chem_id, xa, &
             res, d_dlnd, d_dlnT, d_dxa, ierr)
 
          use eos_def
@@ -279,6 +281,7 @@ subroutine skye_eos( &
          real(dp), intent(in) :: xa(:)
          real(dp), intent(in) :: temp_in, den_in, mass_fraction_limit, Skye_min_gamma_for_solid, Skye_max_gamma_for_liquid
          real(dp), intent(in) :: Xfrac, abar, zbar
+         logical, intent(in) :: use_ion_offsets
          character(len=128), intent(in) :: Skye_solid_mixing_rule
          integer, intent(out) :: ierr
          real(dp), intent(out), dimension(nv) :: res, d_dlnd, d_dlnT
@@ -351,7 +354,9 @@ subroutine skye_eos( &
          ! Ideal ion free energy, only depends on abar
          F_ideal_ion = compute_F_ideal_ion(temp, den, abar, relevant_species, ACMI, ya)
 
-         F_ideal_ion = F_ideal_ion + compute_ion_offset(species, xa, chem_id) ! Offset so ion ground state energy is zero.
+         if (use_ion_offsets) then
+            F_ideal_ion = F_ideal_ion + compute_ion_offset(species, xa, chem_id) ! Offset so ion ground state energy is zero.
+         end if
 
          ! Ideal electron-positron thermodynamics (s, e, p)
          ! Derivatives are handled by HELM code, so we don't pass *in* any auto_diff types (just get them as return values).

eos/public/eos_def.f90

@@ -293,6 +293,7 @@ end subroutine other_eos_interface
 
          ! limits for Skye
          logical :: use_Skye
+         logical :: Skye_use_ion_offsets
          real(dp) :: mass_fraction_limit_for_Skye
          real(dp) :: Skye_min_gamma_for_solid ! The minimum Gamma_i at which to use the solid free energy fit (below this, extrapolate).
          real(dp) :: Skye_max_gamma_for_liquid ! The maximum Gamma_i at which to use the liquid free energy fit (above this, extrapolate).

star/private/evolve.f90

@@ -624,16 +624,11 @@ integer function do_step_part2(id, first_try)
                end if
 
                if(s% do_phase_separation) then
-                  do k=1,s% nz
-                     s% eps_phase_separation(k) = s% energy(k)
-                  end do
-                  call do_phase_separation(s, ierr)
+                  call do_phase_separation(s, dt, ierr)
                   if (failed('do_phase_separation')) return
+                  
                   call set_vars_if_needed(s, dt, 'after phase separation', ierr)
                   if (failed('set_vars_if_needed after phase separation')) return
-                  do k=1,s% nz ! for use by energy equation
-                     s% eps_phase_separation(k) = (s% eps_phase_separation(k) - s% energy(k)) / dt
-                  end do
                else
                   s% crystal_core_boundary_mass = -1d0
                end if

star/private/phase_separation.f90

@@ -36,24 +36,25 @@ module phase_separation
 
       logical, parameter :: dbg = .false.
 
-      integer, parameter :: FIXED_PT_MODE = 5
-      integer, parameter :: FIXED_DT_MODE = 6      
-
       ! offset to higher phase than 0.5 to avoid interference
       ! between phase separation mixing and latent heat for Skye.
       real(dp), parameter :: eos_phase_boundary = 0.9d0
 
       contains
 
 
-      subroutine do_phase_separation(s, ierr)
+      subroutine do_phase_separation(s, dt, ierr)
          use chem_def, only: chem_isos, ic12, io16
          use chem_lib, only: chem_get_iso_id
          type (star_info), pointer :: s
+         real(dp), intent(in) :: dt
          integer, intent(out) :: ierr
 
          real(dp) :: dq_crystal, XO, XC, pad
          integer :: k, k_bound, kstart, net_ic12, net_io16
+         logical :: save_Skye_use_ion_offsets
+
+         s% eps_phase_separation(:) = 0d0
 
          if(s% phase(s% nz) < eos_phase_boundary) then
             s% crystal_core_boundary_mass = 0d0
@@ -93,6 +94,15 @@ subroutine do_phase_separation(s, ierr)
                   exit
                end if
             end do
+
+            ! calculate energy associated with phase separation, ignoring the ionization
+            ! energy term that Skye sometimes calculates
+            save_Skye_use_ion_offsets = s% eos_rq% Skye_use_ion_offsets
+            s% eos_rq% Skye_use_ion_offsets = .false.
+            call update_model_(s,1,s%nz,.false.)
+            do k=1,s% nz
+               s% eps_phase_separation(k) = s% energy(k)
+            end do
 
             ! loop runs outward starting at previous crystallization boundary
             do k = kstart,1,-1
@@ -111,6 +121,13 @@ subroutine do_phase_separation(s, ierr)
 
             end do
 
+            call update_model_(s,1,s%nz,.false.)
+            
+            ! phase separation heating term for use by energy equation
+            do k=1,s% nz
+               s% eps_phase_separation(k) = (s% eps_phase_separation(k) - s% energy(k)) / dt
+            end do
+            s% eos_rq% Skye_use_ion_offsets = save_Skye_use_ion_offsets
             s% need_to_setvars = .true.
          end if
 
@@ -126,9 +143,6 @@ subroutine move_one_zone(s,k,dq_crystal)
 
         real(dp) :: XC, XO, XC1, XO1, dXO, Xfac, dqsum
         integer :: net_ic12, net_io16
-        integer :: update_mode(s%nz)
-        
-        update_mode(:) = FIXED_DT_MODE
 
         dq_crystal = dq_crystal + s% dq(k)
 
@@ -156,7 +170,7 @@ subroutine move_one_zone(s,k,dq_crystal)
         s% xa(net_ic12,k-1) = XC1 + Xfac*dXO * s% dq(k) / s% dq(k-1)
         s% xa(net_io16,k-1) = XO1 - Xfac*dXO * s% dq(k) / s% dq(k-1)
 
-        call update_model_(s,update_mode,k-1,s%nz,.true.)
+        call update_model_(s,k-1,s%nz,.true.)
 
       end subroutine move_one_zone
 
@@ -168,12 +182,10 @@ subroutine mix_outward(s,kbot)
         integer, intent(in)      :: kbot
 
         real(dp) :: avg_xa(s%species)
-        real(dp) :: mass, XHe_out, dXC_top, dXC_bot, dXO_top, dXO_bot, B_term, grada, gradr
+        real(dp) :: mass, dXC_top, dXC_bot, dXO_top, dXO_bot, B_term, grada, gradr
         integer :: k, l, ktop, net_ihe4, net_ic12, net_io16
-        integer :: update_mode(s%nz)
         logical :: use_brunt
 
-        update_mode(:) = FIXED_DT_MODE
         use_brunt = s% phase_separation_mixing_use_brunt
         net_ihe4 = s% net_iso(ihe4)
         net_ic12 = s% net_iso(ic12)
@@ -195,32 +207,14 @@ subroutine mix_outward(s,kbot)
            ! avg_xa = MAX(MIN(avg_xa, 1._dp), 0._dp)
            ! avg_xa = avg_xa/SUM(avg_xa)
 
-           XHe_out = max(s%xa(net_ihe4,ktop),s%xa(net_ihe4,ktop-1))
-           if(XHe_out < s% eos_rq% mass_fraction_limit_for_Skye) then
-              ! ok to mix all species
-              do l = 1, s%species
-                 s%xa(l,ktop:kbot) = avg_xa(l)
-              end do
-           else
-              ! Mixing He can cause energy problems for eps_phase_separation
-              ! when using Skye, so once we encounter enough He that it is
-              ! included in Skye energy calculation, stop mixing all species.
-              ! Instead, just flatten out the O16 profile, and mix in exchange
-              ! for C12.
-              dXO_top = avg_xa(net_io16) - s%xa(net_io16,ktop)
-              dXO_bot = avg_xa(net_io16) - s%xa(net_io16,kbot)
-              s%xa(net_io16,ktop:kbot) = avg_xa(net_io16)
-              dXC_top = -dXO_top
-              dXC_bot = -dXO_bot
-              s%xa(net_ic12,ktop) = s%xa(net_ic12,ktop) + dXC_top
-              s%xa(net_ic12,ktop+1:kbot) = s%xa(net_ic12,kbot) + dXC_bot
-           end if
-
+           do l = 1, s%species
+              s%xa(l,ktop:kbot) = avg_xa(l)
+           end do
 
            ! updates, eos, opacities, mu, etc now that abundances have changed,
            ! but only in the cells near the boundary where we need to check here.
            ! Will call full update over mixed region after exiting loop.
-           call update_model_(s, update_mode, ktop-1, ktop+1, use_brunt)
+           call update_model_(s, ktop-1, ktop+1, use_brunt)
 
            if(use_brunt) then
               B_term = s% unsmoothed_brunt_B(ktop)
@@ -240,7 +234,7 @@ subroutine mix_outward(s,kbot)
         end do
 
         ! Call a final update over all mixed cells now.
-        call update_model_(s, update_mode, ktop, kbot, .true.)
+        call update_model_(s, ktop, kbot, .true.)
 
       end subroutine mix_outward
 
@@ -276,30 +270,35 @@ real(dp) function blouin_delta_xo(Xin)
         blouin_delta_xo = Xnew - Xin
       end function blouin_delta_xo
 
-      subroutine update_model_ (s, update_mode, kc_t, kc_b, do_brunt)
+      subroutine update_model_ (s, kc_t, kc_b, do_brunt)
 
         use turb_info, only: set_mlt_vars
         use brunt, only: do_brunt_B
         use micro
 
         type(star_info), pointer :: s
-        integer, intent(in)      :: update_mode(:)
         integer, intent(in)      :: kc_t
         integer, intent(in)      :: kc_b
         logical, intent(in)      :: do_brunt
 
         integer  :: ierr
         integer  :: kf_t
         integer  :: kf_b
-        
+
+        logical :: mask(s%nz)
+
+        mask(:) = .true.
+
         ! Update the model to reflect changes in the abundances across
-        ! cells kc_t:kc_b
-        
-        call set_eos_with_mask(s, kc_t, kc_b, update_mode==FIXED_DT_MODE, ierr)
+        ! cells kc_t:kc_b (the mask part of this call is unused, mask=true for all zones).
+        ! Do updates at constant (P,T) rather than constant (rho,T).
+        s%fix_Pgas = .true.
+        call set_eos_with_mask(s, kc_t, kc_b, mask, ierr)
         if (ierr /= 0) then
            write(*,*) 'phase_separation: error from call to set_eos_with_mask'
            stop
         end if
+        s%fix_Pgas = .false.
 
         ! Update opacities across cells kc_t:kc_b (this also sets rho_face
         ! and related quantities on faces kc_t:kc_b)        

star/test_suite/wd_cool_0.6M/inlist_wd_cool_0.6M

@@ -69,6 +69,10 @@
       atm_option = 'table'
       atm_table = 'WD_tau_25'
 
+      ! C/O phase separation
+      ! do_phase_separation = .true.
+      ! do_phase_separation_heating = .true.
+
       ! Diffusion controls
       do_element_diffusion = .true.
       diffusion_use_full_net = .true.