<jupyter_start><jupyter_text>This notebook produces IMF integrated metallicity dependent net yields. As an example we produce the default and alternative tables for the single sun+ best fit (which fix the SSP parameter). This notebook is just an example showing illustratively how to produce these yield tables. If you actually want to use these for your Chemical evolution model you will probably want to increase the time_steps and the mass_steps. You could also consult me to get an idea of sensible Parameters and yield sets that one could use. Mostly those ideas will be based on: http://adsabs.harvard.edu/abs/2018ApJ...861...40P<jupyter_code>%pylab inline
import multiprocessing as mp
from Chempy.parameter import ModelParameters
a = ModelParameters()
# Load solar abundances

from Chempy.solar_abundance import solar_abundances
basic_solar = solar_abundances()
getattr(basic_solar, 'Asplund09')()

# Initialise the SSP class with time-steps

time_steps = np.linspace(0.,13.5,521)

# Load the default yields

from Chempy.yields import SN2_feedback, AGB_feedback, SN1a_feedback
basic_sn2 = SN2_feedback()
getattr(basic_sn2, 'Nomoto2013_net')()
basic_1a = SN1a_feedback()
getattr(basic_1a, "Seitenzahl")()
basic_agb = AGB_feedback()
getattr(basic_agb, "Karakas_net_yield")()

# Load the alternative yields
alt_sn2 = SN2_feedback()
getattr(alt_sn2, 'chieffi04_net')()
alt_1a = SN1a_feedback()
getattr(alt_1a, "Thielemann")()
alt_agb = AGB_feedback()
getattr(alt_agb, "Ventura")()

# Use all elements that are traced

elements_to_trace = list(np.unique(basic_agb.elements+basic_sn2.elements+basic_1a.elements+alt_agb.elements+alt_sn2.elements+alt_1a.elements))
#elements_to_trace = ['H','He','C','N','O','Na','Al','Mg','Si','Ca','Ti','Mn','Fe','Ba','Ne']
print('all the traced elements: ',elements_to_trace)

# Producing the SSP birth elemental fractions (here we use solar)

solar_fractions = []
elements = np.hstack(basic_solar.all_elements)
for item in elements_to_trace:
    solar_fractions.append(float(basic_solar.fractions[np.where(elements==item)]))


# yieldset default
a.yield_table_name_sn2 = 'Nomoto2013_net'
a.yield_table_name_agb = 'Karakas_net_yield'
a.yield_table_name_1a = 'Seitenzahl'
# yieldset alternative
#a.yield_table_name_sn2 = 'chieffi04_net'
#a.yield_table_name_agb = 'Ventura'
#a.yield_table_name_1a = 'Thielemann'

# imf parameters
a.only_net_yields_in_process_tables = True
a.imf_type_name = 'Chabrier_2'
#default
a.high_mass_slope = -2.46
#alternative
#a.high_mass_slope = -2.51
a.imf_parameter = (22.8978, 716.4, 0.25, a.high_mass_slope)
a.mmin = 0.1
a.mmax = 100
# 100,000,000 mass steps are smooth enough for 1000 time steps
a.mass_steps = 200000 #2000 # 200000
a.sn2mmin = 8.
a.sn2mmax = 100.
a.bhmmin = float(a.sn2mmax) ## maximum of hypernova
a.bhmmax = float(a.mmax) ## maximum of the IMF

# sn1a delay parameters for maoz
#default
a.N_0 = np.power(10,-3.07)
a.sn1a_time_delay = np.power(10,-0.8)
#alternative
#a.N_0 = np.power(10,-3.49)
#a.sn1a_time_delay = np.power(10,-0.88)

a.sn1a_exponent = 1.12
a.dummy = 0.0
a.sn1a_parameter = [a.N_0,a.sn1a_time_delay,a.sn1a_exponent,a.dummy]
######################## END OF SETTING CHEMPY PARAMETER FOR SSP YIELD TABLE PRODUCTION
list_of_metallicities = np.logspace(-4,-2,2)
list_of_SSP_tables = []
list_of_SSP_tables.append(list_of_metallicities)
list_of_SSP_tables.append(time_steps)

from Chempy.wrapper import SSP_wrap


def create_one_SSP_table(parameters):
    differential_table = True # True is the default Chempy behaviour
    metallicity = parameters
    print(metallicity,a.yield_table_name_sn2)
    basic_ssp = SSP_wrap(a)
    basic_ssp.calculate_feedback(metallicity,list(elements_to_trace),list(solar_fractions),np.copy(time_steps))

    x = basic_ssp.agb_table
    y = basic_ssp.sn1a_table
    z = basic_ssp.sn2_table
    s = basic_ssp.bh_table
    d = basic_ssp.table

    u = np.zeros_like(x)
    names = list(u.dtype.names)

    # here we still add all processes, but we can also make individual process contribution to element
    for j,jtem in enumerate(names):
        u[jtem] = x[jtem] + y[jtem] + z[jtem] + s[jtem]
    if differential_table:
        for el in elements_to_trace:
            d[el] = u[el]
    else:
        for el in elements_to_trace:
            d[el] = np.cumsum(u[el])
        for name in ['mass_of_ms_stars_dying', 'mass_in_remnants', 'sn2', 'sn1a', 'pn', 'bh', 'hydrogen_mass_accreted_onto_white_dwarfs', 'unprocessed_ejecta']:
            d[name] = np.cumsum(d[name])
    return(d)
print("There are %d CPUs on this machine" % mp.cpu_count())
number_processes = max(1,mp.cpu_count() - 1)
pool = mp.Pool(number_processes)
results = pool.map(create_one_SSP_table, list_of_metallicities)
pool.close()
pool.join()

list_of_SSP_tables.append(results)
np.save('data/paper_1_ssp_tables/SSP_tables_default', list_of_SSP_tables)
print('the metallicities for which SSP yield tables were calculated: ',list_of_SSP_tables[0])
print('metallicities: ',len(list_of_SSP_tables[2]))
print('timesteps: ', len(list_of_SSP_tables[2][0]), '= ', len(list_of_SSP_tables[1]))
print('the data type of the SSP yield table: ',list_of_SSP_tables[2][0].dtype)
x = list_of_SSP_tables[2][0]
z = list_of_SSP_tables[2][1]
alpha = 0.5
factor = 1.05

## Actual plotting

fig = plt.figure(figsize=(8.69,6.69), dpi=100)
ax = fig.add_subplot(111)
ax.plot(time_steps,np.cumsum(x["Fe"]),'b', label = 'Z = 0.0001')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(x["Fe"])*0.9) ,s = 'Fe',color = 'b')
ax.plot(time_steps,np.cumsum(z["Fe"]), 'r', label = 'Z = 0.01')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(z["Fe"])) ,s = 'Fe',color = 'r')

ax.plot(time_steps,np.cumsum(x["Mg"]),'b')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(x["Mg"])*0.9) ,s = 'Mg',color = 'b')
ax.plot(time_steps,np.cumsum(z["Mg"]), 'r')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(z["Mg"])) ,s = 'Mg',color = 'r')

ax.plot(time_steps,np.cumsum(x["Al"]),'b')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(x["Al"])*0.9) ,s = 'Al',color = 'b')
ax.plot(time_steps,np.cumsum(z["Al"]), 'r')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(z["Al"])) ,s = 'Al',color = 'r')

ax.plot(time_steps,np.cumsum(x["C"]),'b')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(x["C"])*0.9) ,s = 'C',color = 'b')
ax.plot(time_steps,np.cumsum(z["C"]), 'r')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(z["C"])) ,s = 'C',color = 'r')

ax.plot(time_steps,np.ones_like(time_steps)*5e-3,marker = '|', markersize = 10, linestyle = '', color = 'k', alpha = 2*alpha)#, label = 'time-steps')
ax.annotate(xy = (time_steps[1],2.7e-3),s = 'model time-steps', color = 'k', alpha = 2*alpha)
ax.legend(loc = 'best')

ax.set_ylim(2e-5,6e-3)
ax.set_xlim(7e-3,25)
ax.set_title(r'default yield of SSP with mass = 1M$_\odot$ for different metallicies')
ax.set_ylabel(r"net yield in M$_\odot$")
ax.set_xlabel("time in Gyr")

ax.set_yscale('log')
ax.set_xscale('log')
plt.show()
# yieldset alternative
a.yield_table_name_sn2 = 'chieffi04_net'
a.yield_table_name_agb = 'Ventura'
a.yield_table_name_1a = 'Thielemann'

a.high_mass_slope = -2.51
a.imf_parameter = (22.8978, 716.4, 0.25, a.high_mass_slope)
a.N_0 = np.power(10,-3.49)
a.sn1a_time_delay = np.power(10,-0.88)
a.sn1a_parameter = [a.N_0,a.sn1a_time_delay,a.sn1a_exponent,a.dummy]
######################## END OF SETTING CHEMPY PARAMETER FOR SSP YIELD TABLE PRODUCTION

list_of_SSP_tables_alt = []
list_of_SSP_tables_alt.append(list_of_metallicities)
list_of_SSP_tables_alt.append(time_steps)

print("There are %d CPUs on this machine" % mp.cpu_count())
number_processes = max(1,mp.cpu_count() - 1)
pool = mp.Pool(number_processes)
results = pool.map(create_one_SSP_table, list_of_metallicities)
pool.close()
pool.join()

list_of_SSP_tables_alt.append(results)
np.save('data/paper_1_ssp_tables/SSP_tables_alternative', list_of_SSP_tables_alt)
print('the data type of the SSP yield table: ',list_of_SSP_tables[2][0].dtype)
y = list_of_SSP_tables_alt[2][1]
## Actual plotting

fig = plt.figure(figsize=(8.69,6.69), dpi=100)
ax = fig.add_subplot(111)
ax.plot(time_steps,np.cumsum(y["Fe"]),'b', label = 'alternative')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(y["Fe"])*0.9) ,s = 'Fe',color = 'b')
ax.plot(time_steps,np.cumsum(z["Fe"]), 'r', label = 'default')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(z["Fe"])) ,s = 'Fe',color = 'r')

ax.plot(time_steps,np.cumsum(y["Mg"]),'b')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(y["Mg"])*0.9) ,s = 'Mg',color = 'b')
ax.plot(time_steps,np.cumsum(z["Mg"]), 'r')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(z["Mg"])) ,s = 'Mg',color = 'r')

ax.plot(time_steps,np.cumsum(y["Al"]),'b')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(y["Al"])*0.9) ,s = 'Al',color = 'b')
ax.plot(time_steps,np.cumsum(z["Al"]), 'r')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(z["Al"])) ,s = 'Al',color = 'r')

ax.plot(time_steps,np.cumsum(y["C"]),'b')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(y["C"])*0.9) ,s = 'C',color = 'b')
ax.plot(time_steps,np.cumsum(z["C"]), 'r')
ax.annotate(xy = (time_steps[-1]*factor,np.sum(z["C"])) ,s = 'C',color = 'r')

ax.plot(time_steps,np.ones_like(time_steps)*5e-3,marker = '|', markersize = 10, linestyle = '', color = 'k', alpha = 2*alpha)#, label = 'time-steps')
ax.annotate(xy = (time_steps[1],2.7e-3),s = 'model time-steps', color = 'k', alpha = 2*alpha)
ax.legend(loc = 'best')

ax.set_ylim(2e-5,6e-3)
ax.set_xlim(7e-3,25)
ax.set_title(r'yield of SSP with mass = 1M$_\odot$ for the different yieldsets for Z=0.01')
ax.set_ylabel(r"net yield in M$_\odot$")
ax.set_xlabel("time in Gyr")

ax.set_yscale('log')
ax.set_xscale('log')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17763704318936877, 'fraction_numerical': 0.03343023255813953, 'mean_word_length': 5.885632594710508, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '139854', 'license': 'MIT', 'n_tokens_mistral': 4209, 'n_tokens_neox': 3999, 'n_words': 718, 'path': 'tutorials/8-Yield tables for SPH simulations and comparison to other tables.ipynb', 'repo_name': 'oplevne/Chempy', 'size': '9632'}