metal_return.h File Reference

#include "forcetree.h"
#include "timestep.h"
#include "utils/paramset.h"
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_integration.h>
#include "slotsmanager.h"

Include dependency graph for metal_return.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	interps
struct	MetalReturnPriv

Functions
void	setup_metal_table_interp (struct interps *interp)
void	metal_return (const ActiveParticles *act, DomainDecomp *const ddecomp, Cosmology *CP, const double atime, const double AvgGasMass)
void	set_metal_return_params (ParameterSet *ps)
int64_t	metal_return_init (const ActiveParticles *act, Cosmology *CP, struct MetalReturnPriv *priv, const double atime)
void	metal_return_priv_free (struct MetalReturnPriv *priv)
void	stellar_density (const ActiveParticles *act, MyFloat *StarVolumeSPH, MyFloat *MassReturn, const ForceTree *const tree)
int	metals_haswork (int i, MyFloat *MassReturn)

Function Documentation

metal_return()

void metal_return	(	const ActiveParticles *	act,
		DomainDecomp *const	ddecomp,
		Cosmology *	CP,
		const double	atime,
		const double	AvgGasMass
	)

This function is the driver routine for the calculation of metal return.

Definition at line 517 of file metal_return.c.

{
    /* Do nothing if no stars yet*/
    int64_t totstar;
    MPI_Allreduce(&SlotsManager->info[4].size, &totstar, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
    if(totstar == 0)
        return;
 
    struct MetalReturnPriv priv[1];
 
    int64_t nwork = metal_return_init(act, CP, priv, atime);
 
    /* Maximum mass of a gas particle after enrichment: cap it at a few times the initial mass.
     * FIXME: Ideally we should here fork a new particle with a smaller gas mass. We should
     * figure out then how set the gas entropy. A possibly better idea is to add
     * a generic routine to split gas particles into the density code.*/
    priv->MaxGasMass = 4* AvgGasMass;
 
    int64_t totwork;
    MPI_Allreduce(&nwork, &totwork, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
 
    walltime_measure("/SPH/Metals/Init");
 
    if(totwork == 0) {
        metal_return_priv_free(priv);
        return;
    }
 
    ForceTree gasTree = {0};
    /* Just gas, no moments*/
    force_tree_rebuild_mask(&gasTree, ddecomp, GASMASK, 0, NULL);
 
    /* Compute total number of weights around each star for actively returning stars*/
    stellar_density(act, priv->StarVolumeSPH, priv->MassReturn, &gasTree);
 
    /* Do the metal return*/
    TreeWalk tw[1] = {{0}};
 
    tw->ev_label = "METALS";
    tw->visit = (TreeWalkVisitFunction) treewalk_visit_ngbiter;
    tw->ngbiter = (TreeWalkNgbIterFunction) metal_return_ngbiter;
    tw->ngbiter_type_elsize = sizeof(TreeWalkNgbIterMetals);
    tw->haswork = metal_return_haswork;
    tw->fill = (TreeWalkFillQueryFunction) metal_return_copy;
    tw->reduce = NULL;
    tw->postprocess = (TreeWalkProcessFunction) metal_return_postprocess;
    tw->query_type_elsize = sizeof(TreeWalkQueryMetals);
    tw->result_type_elsize = sizeof(TreeWalkResultMetals);
    tw->repeatdisallowed = 1;
    tw->tree = &gasTree;
    tw->priv = priv;
 
    priv->spin = init_spinlocks(SlotsManager->info[0].size);
    treewalk_run(tw, act->ActiveParticle, act->NumActiveParticle);
    free_spinlocks(priv->spin);
 
    force_tree_free(&gasTree);
    metal_return_priv_free(priv);
 
    /* collect some timing information */
    walltime_measure("/SPH/Metals/Yield");
}

References ActiveParticles::ActiveParticle, CP, TreeWalk::ev_label, TreeWalk::fill, force_tree_free(), force_tree_rebuild_mask(), free_spinlocks(), GASMASK, TreeWalk::haswork, slots_manager_type::info, init_spinlocks(), MetalReturnPriv::MassReturn, MetalReturnPriv::MaxGasMass, metal_return_copy(), metal_return_haswork(), metal_return_init(), metal_return_ngbiter(), metal_return_postprocess(), metal_return_priv_free(), MPI_INT64, TreeWalk::ngbiter, TreeWalk::ngbiter_type_elsize, ActiveParticles::NumActiveParticle, TreeWalk::postprocess, TreeWalk::priv, TreeWalk::query_type_elsize, TreeWalk::reduce, TreeWalk::repeatdisallowed, TreeWalk::result_type_elsize, slot_info::size, SlotsManager, MetalReturnPriv::spin, MetalReturnPriv::StarVolumeSPH, stellar_density(), TreeWalk::tree, treewalk_run(), treewalk_visit_ngbiter(), TreeWalk::visit, and walltime_measure.

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

metal_return_init()

int64_t metal_return_init	(	const ActiveParticles *	act,
		Cosmology *	CP,
		struct MetalReturnPriv *	priv,
		const double	atime
	)

Definition at line 437 of file metal_return.c.

{
    int nthread = omp_get_max_threads();
    priv->gsl_work = ta_malloc("gsl_work", gsl_integration_workspace *, nthread);
    int i;
    /* Allocate a workspace for each thread*/
    for(i=0; i < nthread; i++)
        priv->gsl_work[i] = gsl_integration_workspace_alloc(GSL_WORKSPACE);
    priv->hub = CP->HubbleParam;
 
    /* Initialize*/
    setup_metal_table_interp(&priv->interp);
    priv->StellarAges = (MyFloat *) mymalloc("StellarAges", SlotsManager->info[4].size * sizeof(MyFloat));
    priv->MassReturn = (MyFloat *) mymalloc("MassReturn", SlotsManager->info[4].size * sizeof(MyFloat));
    priv->LowDyingMass = (MyFloat *) mymalloc("LowDyingMass", SlotsManager->info[4].size * sizeof(MyFloat));
    priv->HighDyingMass = (MyFloat *) mymalloc("HighDyingMass", SlotsManager->info[4].size * sizeof(MyFloat));
    priv->StarVolumeSPH = (MyFloat *) mymalloc("StarVolumeSPH", SlotsManager->info[4].size * sizeof(MyFloat));
 
    priv->imf_norm = compute_imf_norm(priv->gsl_work[0]);
    /* Maximum possible mass return for below*/
    double maxmassfrac = mass_yield(0, 1/(CP->HubbleParam*HUBBLE * SEC_PER_MEGAYEAR), snii_metallicities[SNII_NMET-1], CP->HubbleParam, &priv->interp, priv->imf_norm, priv->gsl_work[0],agb_masses[0], MAXMASS);
 
    int64_t haswork = 0;
    /* First find the mass return as a fraction of the total mass and the age of the star.
     * This is done first so we can skip density computation for not active stars*/
    #pragma omp parallel for reduction(+: haswork)
    for(i=0; i < act->NumActiveParticle;i++)
    {
        int p_i = act->ActiveParticle ? act->ActiveParticle[i] : i;
        if(P[p_i].Type != 4)
            continue;
        int tid = omp_get_thread_num();
        const int slot = P[p_i].PI;
        priv->StellarAges[slot] = atime_to_myr(CP, STARP(p_i).FormationTime, atime, priv->gsl_work[tid]);
        /* Note this takes care of units*/
        double initialmass = P[p_i].Mass + STARP(p_i).TotalMassReturned;
        find_mass_bin_limits(&priv->LowDyingMass[slot], &priv->HighDyingMass[slot], STARP(p_i).LastEnrichmentMyr, priv->StellarAges[P[p_i].PI], STARP(p_i).Metallicity, priv->interp.lifetime_interp);
 
        priv->MassReturn[slot] = initialmass * mass_yield(STARP(p_i).LastEnrichmentMyr, priv->StellarAges[P[p_i].PI], STARP(p_i).Metallicity, CP->HubbleParam, &priv->interp, priv->imf_norm, priv->gsl_work[tid],priv->LowDyingMass[slot], priv->HighDyingMass[slot]);
        //message(3, "Particle %d PI %d massgen %g mass %g initmass %g\n", p_i, P[p_i].PI, priv->MassReturn[P[p_i].PI], P[p_i].Mass, initialmass);
        /* Guard against making a zero mass particle and warn since this should not happen.*/
        if(STARP(p_i).TotalMassReturned + priv->MassReturn[slot] > initialmass * maxmassfrac) {
            if(priv->MassReturn[slot] / STARP(p_i).TotalMassReturned > 0.01)
                message(1, "Large mass return id %ld %g from %d mass %g initial %g (maxfrac %g) age %g lastenrich %g metal %g dymass %g %g\n",
                    P[p_i].ID, priv->MassReturn[slot], p_i, STARP(p_i).TotalMassReturned, initialmass, maxmassfrac, priv->StellarAges[P[p_i].PI], STARP(p_i).LastEnrichmentMyr, STARP(p_i).Metallicity, priv->LowDyingMass[slot], priv->HighDyingMass[slot]);
            priv->MassReturn[slot] = initialmass * maxmassfrac - STARP(p_i).TotalMassReturned;
            if(priv->MassReturn[slot] < 0) {
                priv->MassReturn[slot] = 0;
            }
            /* Ensure that we skip this step*/
            if(!metals_haswork(p_i, priv->MassReturn))
                STARP(p_i).LastEnrichmentMyr = priv->StellarAges[P[p_i].PI];
 
        }
        /* Keep count of how much work we need to do*/
        if(metals_haswork(p_i, priv->MassReturn))
            haswork++;
    }
    return haswork;
}

References ActiveParticles::ActiveParticle, agb_masses, atime_to_myr(), compute_imf_norm(), CP, find_mass_bin_limits(), MetalReturnPriv::gsl_work, GSL_WORKSPACE, MetalReturnPriv::HighDyingMass, MetalReturnPriv::hub, HUBBLE, Cosmology::HubbleParam, MetalReturnPriv::imf_norm, slots_manager_type::info, MetalReturnPriv::interp, interps::lifetime_interp, MetalReturnPriv::LowDyingMass, mass_yield(), MetalReturnPriv::MassReturn, MAXMASS, message(), metals_haswork(), mymalloc, ActiveParticles::NumActiveParticle, P, SEC_PER_MEGAYEAR, setup_metal_table_interp(), slot_info::size, SlotsManager, snii_metallicities, SNII_NMET, STARP, MetalReturnPriv::StarVolumeSPH, MetalReturnPriv::StellarAges, and ta_malloc.

Referenced by metal_return().

Here is the call graph for this function:

Here is the caller graph for this function:

metal_return_priv_free()

void metal_return_priv_free

(

struct MetalReturnPriv *

priv

)

Definition at line 500 of file metal_return.c.

{
    myfree(priv->StarVolumeSPH);
    myfree(priv->HighDyingMass);
    myfree(priv->LowDyingMass);
    myfree(priv->MassReturn);
    myfree(priv->StellarAges);
 
    int i;
    for(i=0; i < omp_get_max_threads(); i++)
        gsl_integration_workspace_free(priv->gsl_work[i]);
 
    ta_free(priv->gsl_work);
}

References MetalReturnPriv::gsl_work, MetalReturnPriv::HighDyingMass, MetalReturnPriv::LowDyingMass, MetalReturnPriv::MassReturn, myfree, MetalReturnPriv::StarVolumeSPH, MetalReturnPriv::StellarAges, and ta_free.

Referenced by metal_return().

Here is the caller graph for this function:

metals_haswork()

int metals_haswork	(	int	i,
		MyFloat *	MassReturn
	)

Definition at line 708 of file metal_return.c.

{
    if(P[i].Type != 4)
        return 0;
    int pi = P[i].PI;
    /* Don't do enrichment from all stars, just those with significant enrichment*/
    if(MassReturn[pi] < 1e-3 * (P[i].Mass + STARP(i).TotalMassReturned))
        return 0;
    return 1;
}

References MetalReturnPriv::MassReturn, P, and STARP.

Referenced by metal_return_haswork(), metal_return_init(), and stellar_density_haswork().

Here is the caller graph for this function:

set_metal_return_params()

void set_metal_return_params

(

ParameterSet *

ps

)

Definition at line 62 of file metal_return.c.

{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0) {
        MetalParams.Sn1aN0 = param_get_double(ps, "MetalsSn1aN0");
        MetalParams.SPHWeighting = param_get_int(ps, "MetalsSPHWeighting");
        MetalParams.MaxNgbDeviation = param_get_double(ps, "MetalsMaxNgbDeviation");
    }
    MPI_Bcast(&MetalParams, sizeof(struct metal_return_params), MPI_BYTE, 0, MPI_COMM_WORLD);
}

References metal_return_params::MaxNgbDeviation, MetalParams, param_get_double(), param_get_int(), metal_return_params::Sn1aN0, metal_return_params::SPHWeighting, and ThisTask.

Referenced by read_parameter_file().

Here is the call graph for this function:

Here is the caller graph for this function:

setup_metal_table_interp()

void setup_metal_table_interp

(

struct interps *

interp

)

Definition at line 76 of file metal_return.c.

{
    interp->lifetime_interp = gsl_interp2d_alloc(gsl_interp2d_bilinear, LIFE_NMET, LIFE_NMASS);
    gsl_interp2d_init(interp->lifetime_interp, lifetime_metallicity, lifetime_masses, lifetime, LIFE_NMET, LIFE_NMASS);
    interp->agb_mass_interp = gsl_interp2d_alloc(gsl_interp2d_bilinear, AGB_NMET, AGB_NMASS);
    gsl_interp2d_init(interp->agb_mass_interp, agb_metallicities, agb_masses, agb_total_mass, AGB_NMET, AGB_NMASS);
    interp->agb_metallicity_interp = gsl_interp2d_alloc(gsl_interp2d_bilinear, AGB_NMET, AGB_NMASS);
    gsl_interp2d_init(interp->agb_metallicity_interp, agb_metallicities, agb_masses, agb_total_metals, AGB_NMET, AGB_NMASS);
    int i;
    for(i=0; i<NMETALS; i++) {
        interp->agb_metals_interp[i] = gsl_interp2d_alloc(gsl_interp2d_bilinear, AGB_NMET, AGB_NMASS);
        gsl_interp2d_init(interp->agb_metals_interp[i], agb_metallicities, agb_masses, agb_yield[i], AGB_NMET, AGB_NMASS);
    }
    interp->snii_mass_interp = gsl_interp2d_alloc(gsl_interp2d_bilinear, SNII_NMET, SNII_NMASS);
    gsl_interp2d_init(interp->snii_mass_interp, snii_metallicities, snii_masses, snii_total_mass, SNII_NMET, SNII_NMASS);
    interp->snii_metallicity_interp = gsl_interp2d_alloc(gsl_interp2d_bilinear, SNII_NMET, SNII_NMASS);
    gsl_interp2d_init(interp->snii_metallicity_interp, snii_metallicities, snii_masses, snii_total_metals, SNII_NMET, SNII_NMASS);
    for(i=0; i<NMETALS; i++) {
        interp->snii_metals_interp[i] = gsl_interp2d_alloc(gsl_interp2d_bilinear, SNII_NMET, SNII_NMASS);
        gsl_interp2d_init(interp->snii_metals_interp[i], snii_metallicities, snii_masses, snii_yield[i], SNII_NMET, SNII_NMASS);
    }
}

References agb_masses, agb_metallicities, AGB_NMASS, AGB_NMET, agb_total_mass, agb_total_metals, agb_yield, interp, LIFE_NMASS, LIFE_NMET, lifetime, lifetime_masses, lifetime_metallicity, NMETALS, snii_masses, snii_metallicities, SNII_NMASS, SNII_NMET, snii_total_mass, snii_total_metals, and snii_yield.

Referenced by metal_return_init(), and test_yields().

Here is the caller graph for this function:

stellar_density()

void stellar_density	(	const ActiveParticles *	act,
		MyFloat *	StarVolumeSPH,
		MyFloat *	MassReturn,
		const ForceTree *const	tree
	)

Definition at line 923 of file metal_return.c.

{
    TreeWalk tw[1] = {{0}};
    struct StellarDensityPriv priv[1];
 
    tw->ev_label = "STELLAR_DENSITY";
    tw->visit = treewalk_visit_nolist_ngbiter;
    tw->NoNgblist = 1;
    tw->ngbiter_type_elsize = sizeof(TreeWalkNgbIterStellarDensity);
    tw->ngbiter = (TreeWalkNgbIterFunction) stellar_density_ngbiter;
    tw->haswork = stellar_density_haswork;
    tw->fill = (TreeWalkFillQueryFunction) stellar_density_copy;
    tw->reduce = (TreeWalkReduceResultFunction) stellar_density_reduce;
    tw->postprocess = (TreeWalkProcessFunction) stellar_density_check_neighbours;
    tw->query_type_elsize = sizeof(TreeWalkQueryStellarDensity);
    tw->result_type_elsize = sizeof(TreeWalkResultStellarDensity);
    tw->priv = priv;
    tw->tree = tree;
 
    int i;
 
    priv->MassReturn = MassReturn;
 
    priv->Left = (MyFloat *) mymalloc("DENS_PRIV->Left", SlotsManager->info[4].size * sizeof(MyFloat));
    priv->Right = (MyFloat *) mymalloc("DENS_PRIV->Right", SlotsManager->info[4].size * sizeof(MyFloat));
    priv->NumNgb = (MyFloat (*) [NHSML]) mymalloc("DENS_PRIV->NumNgb", SlotsManager->info[4].size * sizeof(priv->NumNgb[0]));
    priv->VolumeSPH = (MyFloat (*) [NHSML]) mymalloc("DENS_PRIV->VolumeSPH", SlotsManager->info[4].size * sizeof(priv->VolumeSPH[0]));
    priv->maxcmpte = (int *) mymalloc("maxcmpte", SlotsManager->info[4].size * sizeof(int));
 
    priv->DesNumNgb = GetNumNgb(GetDensityKernelType());
 
    #pragma omp parallel for
    for(i = 0; i < act->NumActiveParticle; i++) {
        int a = act->ActiveParticle ? act->ActiveParticle[i] : i;
        /* Skip the garbage particles */
        if(P[a].IsGarbage)
            continue;
        if(!stellar_density_haswork(a, tw))
            continue;
        int pi = P[a].PI;
        priv->Left[pi] = 0;
        priv->Right[pi] = tree->BoxSize;
    }
 
    /* allocate buffers to arrange communication */
 
    treewalk_do_hsml_loop(tw, act->ActiveParticle, act->NumActiveParticle, 1);
    #pragma omp parallel for
    for(i = 0; i < act->NumActiveParticle; i++) {
        int a = act->ActiveParticle ? act->ActiveParticle[i] : i;
        /* Skip the garbage particles */
        if(P[a].IsGarbage)
            continue;
        if(!stellar_density_haswork(a, tw))
            continue;
        /* Copy the Star Volume SPH*/
        StarVolumeSPH[P[a].PI] = priv->VolumeSPH[P[a].PI][0];
        if(priv->VolumeSPH[P[a].PI][0] == 0)
            endrun(3, "i = %d pi = %d StarVolumeSPH %g hsml %g\n", a, P[a].PI, priv->VolumeSPH[P[a].PI][0], P[a].Hsml);
    }
 
    myfree(priv->maxcmpte);
    myfree(priv->VolumeSPH);
    myfree(priv->NumNgb);
    myfree(priv->Right);
    myfree(priv->Left);
 
    double timeall = walltime_measure(WALLTIME_IGNORE);
 
    double timecomp = tw->timecomp3 + tw->timecomp1 + tw->timecomp2;
    double timewait = tw->timewait1 + tw->timewait2;
    double timecomm = tw->timecommsumm1 + tw->timecommsumm2;
    walltime_add("/SPH/Metals/Density/Compute", timecomp);
    walltime_add("/SPH/Metals/Density/Wait", timewait);
    walltime_add("/SPH/Metals/Density/Comm", timecomm);
    walltime_add("/SPH/Metals/Density/Misc", timeall - (timecomp + timewait + timecomm));
 
    return;
}

References ActiveParticles::ActiveParticle, ForceTree::BoxSize, StellarDensityPriv::DesNumNgb, endrun(), TreeWalk::ev_label, TreeWalk::fill, GetDensityKernelType(), GetNumNgb(), TreeWalk::haswork, slots_manager_type::info, StellarDensityPriv::Left, StellarDensityPriv::MassReturn, StellarDensityPriv::maxcmpte, myfree, mymalloc, TreeWalk::ngbiter, TreeWalk::ngbiter_type_elsize, NHSML, TreeWalk::NoNgblist, ActiveParticles::NumActiveParticle, StellarDensityPriv::NumNgb, P, TreeWalk::postprocess, TreeWalk::priv, TreeWalk::query_type_elsize, TreeWalk::reduce, TreeWalk::result_type_elsize, StellarDensityPriv::Right, slot_info::size, SlotsManager, stellar_density_check_neighbours(), stellar_density_copy(), stellar_density_haswork(), stellar_density_ngbiter(), stellar_density_reduce(), TreeWalk::timecommsumm1, TreeWalk::timecommsumm2, TreeWalk::timecomp1, TreeWalk::timecomp2, TreeWalk::timecomp3, TreeWalk::timewait1, TreeWalk::timewait2, TreeWalk::tree, treewalk_do_hsml_loop(), treewalk_visit_nolist_ngbiter(), TreeWalk::visit, StellarDensityPriv::VolumeSPH, walltime_add, WALLTIME_IGNORE, and walltime_measure.

Referenced by metal_return().

Here is the call graph for this function:

Here is the caller graph for this function: