hombit / mesa2py

Yet another Python binding to MESA
0 stars 1 forks source link

'kap' source file and 'kap' settings choice #6

Closed AndreyTavleev closed 2 years ago

AndreyTavleev commented 2 years ago

Maybe, user wants to manually set the opacity source table. It's controlled by kap_file_prefix for higher temperature, hydrogen-rich conditions, kap_lowT_prefix for low temperatures. Also there is kap_CO_prefix which select the opacity source for higher temperature, hydrogen-poor/metal-rich conditions.

These 3 prefixes are in opacity.f90. Note, that its values are non-default.

Actually there are a lot of prefixes, which define opacity source and blend options, see MESA documentation. Note, that this documentation is for the latest MESA release.

hombit commented 2 years ago

It looks like we can set which prefixes to use once only - when we initialize the first Opac instance. I propose two options: 1) use default options always (updated to modern and self-consistent versions in #8), 2) make it customizable via env variables like MESA_KAP_FILE_PREFIX

AndreyTavleev commented 2 years ago

About #8: see the opacity.f90:kap_init()

The main parameters of kap_init() are (from kap_lib.f90):

subroutine kap_init( &
            kap_file_prefix, CO_prefix, lowT_prefix, &
            blend_logT_upper_bdy, blend_logT_lower_bdy, &
            use_cache, kap_cache_dir, &
            kap_config_file, kap_show_info, ierr)

But blend_logT_upper_bdy=0d0 and blend_logT_lower_bdy=0d0. Default in latest release:

kap_blend_logT_upper_bdy = 3.88d0
kap_blend_logT_lower_bdy = 3.80d0

These parameters define the boundary between tables. Seems that we don't use lowT tables!

AndreyTavleev commented 2 years ago

About opacity.f90:kap_set_choices(). The parameters are (from kap_lib.f90):

subroutine kap_set_choices( &
            handle, cubic_interpolation_in_X, cubic_interpolation_in_Z, &
            include_electron_conduction, &
            use_Zbase_for_Type1, use_Type2_opacities, &
            kap_Type2_full_off_X, kap_Type2_full_on_X, &
            kap_Type2_full_off_dZ, kap_Type2_full_on_dZ, &
            ierr)

By default cubic_interpolation_in_X=.false. and cubic_interpolation_in_Z=.false.. I think our values (.true.) are better. But kap_Type2_full_off_X=0.71_dp, kap_Type2_full_on_X=0.70_dp, so Type2 opacities (with CO) are used for X in [0.7, 0.71], but these tables are for low X! Default in latest release:

kap_Type2_full_off_X = 1d-3
kap_Type2_full_on_X = 1d-6

kap_Type2_full_off_dZ and kap_Type2_full_on_dZ are the same as default (0.001 < (dZ = Z - Zbase=Z+1) < 0.01). Also we can see, that Type2 opacities are used not very often.

AndreyTavleev commented 2 years ago

Important not: The OPAL Type 2 opacity tables are on by default (use_Type2_opacities = .true.) since r15140 release. It's the next release after one for mesa2py (r12778).

AndreyTavleev commented 2 years ago

About opacity.f90:kap_get(). The parameters are (from kap_lib.f90):

subroutine kap_get( &
            handle, zbar, X, Z, Zbase, XC, XN, XO, XNe, logRho, logT, &
            lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
            frac_Type2, kap, dlnkap_dlnRho, dlnkap_dlnT,   ierr)

First, I forget what zbar means.

Second, we have Z=Zbase=Z, but by default Zbase=-1, see documentation:

Zbase

The base metallicity for the opacity tables. 
This provides the reference metallicity necessary to calculate element variations. 
Physically, this usually corresponds to the initial metallicity of the star.

If `use_Type2_opacities=.true.`, then if use_Type2_opacities = .true., 
Type1 opacities will be computed using Zbase instead of Z when Z > Zbase. 
This helps with blending from Type1 to Type2. Ignored if use_Type2_opacities = .false..

I think, we should use the default value Zbase=-1!=Z, because (in theory) we may use Type2 opacity tables.

Third, I don't know, what frac_Type2 means. I think it's an output or internal parameter, so there are no questions about it. UPD: Find the definition of frac_Type2 in the documentation after defining kap_Type2_full_on_dZ parameters.

AndreyTavleev commented 2 years ago

Also please check me, using documentation for the latest MESA release.

AndreyTavleev commented 2 years ago

It looks like we can set which prefixes to use once only - when we initialize the first Opac instance. I propose two options: 1) use default options always (updated to modern and self-consistent versions in #8), 2) make it customizable via env variables like MESA_KAP_FILE_PREFIX

Can we make two classes: Opac_OPAL and Opac_OP with different table prefixes. So user can change tables by changing the opacity module.

AndreyTavleev commented 2 years ago

Don't forget: let's create the variables like kap_Type2_full_off_X=... instead of just writing its values to function call. For convenience.

AndreyTavleev commented 2 years ago

FUCK! In our release of MESA there are very different defaults. I don't understand, what should we use. kap_def.f90 in our release:

  subroutine init_kap_handle_data(i) ! set defaults
    integer, intent(in) :: i
    kap_handles(i)% cubic_interpolation_in_X = .true.
    kap_handles(i)% cubic_interpolation_in_Z = .false. 
    kap_handles(i)% include_electron_conduction = .true.
    kap_handles(i)% use_Zbase_for_Type1 = .true.
    kap_handles(i)% use_Type2_opacities = .false.
    kap_handles(i)% min_logT_for_logR_gt_1 = 3.5d0
    kap_handles(i)% kap_Type2_full_off_X = 0.71d0 ! Type2 full off for X >= this
    kap_handles(i)% kap_Type2_full_on_X = 0.70d0 ! Type2 full on for X <= this
    kap_handles(i)% kap_Type2_full_off_dZ = 0.001d0 ! Type2 is full off for dZ <= this
    kap_handles(i)% kap_Type2_full_on_dZ = 0.01d0 ! Type2 can be full on for dZ >= this
  end subroutine init_kap_handle_data
AndreyTavleev commented 2 years ago

About opacity.f90:kap_get(). The parameters are (from kap_lib.f90):

subroutine kap_get( &
            handle, zbar, X, Z, Zbase, XC, XN, XO, XNe, logRho, logT, &
            lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
            frac_Type2, kap, dlnkap_dlnRho, dlnkap_dlnT,   ierr)

First, I forget what zbar means.

Second, we have Z=Zbase=Z, but by default Zbase=-1, see documentation:

Zbase

The base metallicity for the opacity tables. 
This provides the reference metallicity necessary to calculate element variations. 
Physically, this usually corresponds to the initial metallicity of the star.

If `use_Type2_opacities=.true.`, then if use_Type2_opacities = .true., 
Type1 opacities will be computed using Zbase instead of Z when Z > Zbase. 
This helps with blending from Type1 to Type2. Ignored if use_Type2_opacities = .false..

I think, we should use the default value Zbase=-1!=Z, because (in theory) we may use Type2 opacity tables.

Third, I don't know, what frac_Type2 means. I think it's an output or internal parameter, so there are no questions about it. UPD: Find the definition of frac_Type2 in the documentation after defining kap_Type2_full_on_dZ parameters.

I realized, that default Zbase=-1 is just for variable declaration. We should use the solar value, i.e. Zbase=0.018987109631893317d0 != Z, because if Zbase==Z we never use Type2 opacity tables.