cooling_qso_lightup.c File Reference

#include <math.h>
#include <mpi.h>
#include <string.h>
#include <omp.h>
#include <gsl/gsl_interp.h>
#include "physconst.h"
#include "slotsmanager.h"
#include "partmanager.h"
#include "treewalk.h"
#include "hydra.h"
#include "drift.h"
#include "walltime.h"
#include "fof.h"
#include "utils/endrun.h"
#include "utils/paramset.h"
#include "utils/mymalloc.h"
#include "cooling_qso_lightup.h"

Include dependency graph for cooling_qso_lightup.c:

Go to the source code of this file.

Classes
struct	TreeWalkQueryQSOLightup
struct	qso_lightup_params
struct	QSOPriv

Macros
#define	E0_HeII   54.4 /* HeII ionization potential in eV*/
#define	HEMASS   4.002602 /* Helium mass in amu*/
#define	QSO_GET_PRIV(tw)   ((struct QSOPriv *) (tw->priv))

Functions
void	set_qso_lightup_par (struct qso_lightup_params qso)
void	set_qso_lightup_params (ParameterSet *ps)
static double	Q_inst (double Emax, double alpha_q)
static void	load_heii_reion_hist (const char *reion_hist_file)
void	init_qso_lightup (char *reion_hist_file)
double	get_long_mean_free_path_heating (double redshift)
static double	gaussian_rng (double mu, double sigma, const int64_t seed)
static int	build_qso_candidate_list (int **qso_cand, FOFGroups *fof)
static int	count_QSO_halos (int ncand, int64_t *ncand_tot, MPI_Comm Comm)
static int	choose_QSO_halo (int64_t ncand, int64_t *ncand_before, int64_t *ncand_tot, int64_t randseed)
static double	gas_ionization_fraction (void)
static int	ionize_single_particle (int other, double a3inv, double uu_in_cgs)
static void	ionize_ngbiter (TreeWalkQueryQSOLightup *I, TreeWalkResultBase *O, TreeWalkNgbIterBase *iter, LocalTreeWalk *lv)
static void	ionize_copy (int place, TreeWalkQueryQSOLightup *I, TreeWalk *tw)
static int64_t	ionize_all_part (int qso_ind, int *qso_cand, struct QSOPriv priv, ForceTree *tree)
static void	turn_on_quasars (double atime, FOFGroups *fof, DomainDecomp *ddecomp, Cosmology *CP, double uu_in_cgs, FILE *FdHelium)
void	do_heiii_reionization (double atime, FOFGroups *fof, DomainDecomp *ddecomp, Cosmology *CP, double uu_in_cgs, FILE *FdHelium)
int	need_change_helium_ionization_fraction (double atime)
int	during_helium_reionization (double redshift)
int	qso_lightup_on (void)

Variables
static struct qso_lightup_params	QSOLightupParams
static double	qso_inst_heating
static int	Nreionhist
static double *	He_zz
static double *	XHeIII
static double *	LMFP
static gsl_interp *	HeIII_intp
static gsl_interp *	LMFP_intp
static double	last_zz
static double	last_long_mfp_heating

Macro Definition Documentation

E0_HeII

#define E0_HeII   54.4 /* HeII ionization potential in eV*/

This file implements the quasar heating model of Upton-Sanderbeck et al 2019 (in prep). A black hole particle with the right mass is chosen at random and an ionizing bubble created around it. Particles within that bubble are marked as ionized and heated according to emissivity. New bubbles are created until the total HeIII fraction matches the value in the external table. There is also a uniform background heating rate for long mean-free-path photons with very high energies. This heating is added only to not-yet-ionized particles, and is done in cooling_rates.c. Whatever UVB you use should not include these photons.

The text file contains the reionization history and is generated from various physical processes implemented in a python file. The text file fixes the end of helium reionization (which is reasonably well-known). The code has the start time of reionization as a free parameter: if this start of reionization is late then there will be a strong sudden burst of heating.

(Loosely) based on lightup_QSOs.py from https://github.com/uptonsanderbeck/helium_reionization HeII_heating.py contains the details of the reionization history and how it is generated.

This code should run only during a PM timestep, when all particles are active We lose time resolution but I cannot think of another way to ensure cooling is modelled correctly.

Definition at line 47 of file cooling_qso_lightup.c.

HEMASS

#define HEMASS   4.002602 /* Helium mass in amu*/

Definition at line 48 of file cooling_qso_lightup.c.

QSO_GET_PRIV

#define QSO_GET_PRIV

(

tw

)

((struct QSOPriv *) (tw->priv))

Definition at line 413 of file cooling_qso_lightup.c.

Function Documentation

build_qso_candidate_list()

static int build_qso_candidate_list ( int **  qso_cand, FOFGroups *  fof  )

static

Definition at line 282 of file cooling_qso_lightup.c.

{
    /*Loop over all halos, building the candidate list.*/
    int i, ncand=0;
    *qso_cand = (int *) mymalloc("Quasar_candidates", sizeof(int) * (fof->Ngroups+1));
    for(i = 0; i < fof->Ngroups; i++)
    {
        /* Check that it has the right mass*/
        if(fof->Group[i].Mass < QSOLightupParams.qso_candidate_min_mass)
            continue;
        if(fof->Group[i].Mass > QSOLightupParams.qso_candidate_max_mass)
            continue;
        /*Add to the candidate list*/
        (*qso_cand)[ncand] = i;
        ncand++;
    }
    /*Poison value at the end for safety.*/
    (*qso_cand)[ncand] = -1;
    *qso_cand = (int *) myrealloc(*qso_cand, (ncand+1) * sizeof(int));
    return ncand;
}

References FOFGroups::Group, Group::Mass, mymalloc, myrealloc, FOFGroups::Ngroups, qso_lightup_params::qso_candidate_max_mass, qso_lightup_params::qso_candidate_min_mass, and QSOLightupParams.

Referenced by turn_on_quasars().

Here is the caller graph for this function:

choose_QSO_halo()

static int choose_QSO_halo ( int64_t  ncand, int64_t *  ncand_before, int64_t *  ncand_tot, int64_t  randseed  )

static

Definition at line 340 of file cooling_qso_lightup.c.

{
    double drand = get_random_number(randseed);
    int64_t qso = drand * (*ncand_tot);
    (*ncand_tot)--;
    /* No quasar on this processor*/
    if(qso < *ncand_before)
        (*ncand_before)--;
    if(qso < *ncand_before || qso >= *ncand_before + ncand)
        return -1;
 
    /* If the quasar is on this processor, return the
     * index of the quasar in the current candidate array.*/
    return qso - *ncand_before;
}

References get_random_number().

Referenced by turn_on_quasars().

Here is the call graph for this function:

Here is the caller graph for this function:

count_QSO_halos()

static int count_QSO_halos ( int  ncand, int64_t *  ncand_tot, MPI_Comm  Comm  )

static

Definition at line 309 of file cooling_qso_lightup.c.

{
    int64_t ncand_total = 0, ncand_before = 0;
    int NTask, i, ThisTask;
    MPI_Comm_size(Comm, &NTask);
    MPI_Comm_rank(Comm, &ThisTask);
    int * candcounts = (int*) ta_malloc("qso_cand_counts", int, NTask);
 
    /* Get how many candidates are on each processor. TODO: Make a generic MPIU for this.*/
    MPI_Allgather(&ncand, 1, MPI_INT, candcounts, 1, MPI_INT, Comm);
 
    for(i = 0; i < NTask; i++)
    {
        if(i < ThisTask)
            ncand_before += candcounts[i];
        ncand_total += candcounts[i];
    }
 
    ta_free(candcounts);
    *ncand_tot = ncand_total;
    return ncand_before;
}

References NTask, ta_free, ta_malloc, and ThisTask.

Referenced by turn_on_quasars().

Here is the caller graph for this function:

do_heiii_reionization()

void do_heiii_reionization	(	double	atime,
		FOFGroups *	fof,
		DomainDecomp *	ddecomp,
		Cosmology *	CP,
		double	uu_in_cgs,
		FILE *	FdHelium
	)

Definition at line 640 of file cooling_qso_lightup.c.

{
    if(!QSOLightupParams.QSOLightupOn)
        return;
    if(atime < 1/(1+QSOLightupParams.heIIIreion_start))
        return;
 
    /* Do nothing if we are past the end of the table.*/
    if(atime > He_zz[Nreionhist-1])
        return;
 
    walltime_measure("/Misc");
    //message(0, "HeII: Reionization initiated.\n");
    turn_on_quasars(atime, fof, ddecomp, CP, uu_in_cgs, FdHelium);
}

References CP, QSOPriv::fof, He_zz, qso_lightup_params::heIIIreion_start, Nreionhist, qso_lightup_params::QSOLightupOn, QSOLightupParams, turn_on_quasars(), QSOPriv::uu_in_cgs, and walltime_measure.

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

during_helium_reionization()

int during_helium_reionization

(

double

redshift

)

Definition at line 667 of file cooling_qso_lightup.c.

{
    if(!QSOLightupParams.QSOLightupOn)
        return 0;
    if(redshift > QSOLightupParams.heIIIreion_start)
        return 0;
 
    /* Past the end of the table, it has finished.*/
    if(redshift < 1./He_zz[Nreionhist-1] - 1)
        return 0;
 
    return 1;
}

References He_zz, qso_lightup_params::heIIIreion_start, Nreionhist, qso_lightup_params::QSOLightupOn, and QSOLightupParams.

Referenced by run().

Here is the caller graph for this function:

gas_ionization_fraction()

static double gas_ionization_fraction ( void  )

static

Definition at line 359 of file cooling_qso_lightup.c.

{
    int64_t n_ionized_tot = 0, n_gas_tot = 0;
    int i, n_ionized = 0;
    #pragma omp parallel for reduction(+:n_ionized)
    for (i = 0; i < PartManager->NumPart; i++){
        if (P[i].Type == 0 && P[i].HeIIIionized == 1){
            n_ionized ++;
        }
    }
    /* Get total ionization fraction: note this is only the current gas particles.
     * Particles that become stars are not counted.*/
    sumup_large_ints(1, &n_ionized, &n_ionized_tot);
    MPI_Allreduce(&SlotsManager->info[0].size, &n_gas_tot, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
    return (double) n_ionized_tot / (double) n_gas_tot;
}

References slots_manager_type::info, MPI_INT64, part_manager_type::NumPart, P, PartManager, slot_info::size, SlotsManager, and sumup_large_ints().

Referenced by need_change_helium_ionization_fraction(), and turn_on_quasars().

Here is the call graph for this function:

Here is the caller graph for this function:

gaussian_rng()

static double gaussian_rng ( double  mu, double  sigma, const int64_t  seed  )

static

Definition at line 271 of file cooling_qso_lightup.c.

{
    double u1 = get_random_number(seed);
    double u2 = get_random_number(seed + 1);
    double z1 = sqrt(-2 * log(u1) ) * cos(2 * M_PI * u2);
    return mu + sigma * z1;
}

References get_random_number(), and sigma.

Referenced by ionize_ngbiter().

Here is the call graph for this function:

Here is the caller graph for this function:

get_long_mean_free_path_heating()

double get_long_mean_free_path_heating

(

double

redshift

)

Definition at line 249 of file cooling_qso_lightup.c.

{
    if(!QSOLightupParams.QSOLightupOn)
        return 0;
    if(redshift > QSOLightupParams.heIIIreion_start)
        return 0;
    if(redshift == last_zz)
        return last_long_mfp_heating;
    double atime = 1/(1+redshift);
 
    /* Guard against the end of the table*/
    if(atime > He_zz[Nreionhist-1])
        return 0;
 
    double long_mfp_heating = gsl_interp_eval(LMFP_intp, He_zz, LMFP, atime, NULL);
 
    last_zz = redshift;
    last_long_mfp_heating = long_mfp_heating;
    return long_mfp_heating;
}

References He_zz, qso_lightup_params::heIIIreion_start, last_long_mfp_heating, last_zz, LMFP, LMFP_intp, Nreionhist, qso_lightup_params::QSOLightupOn, and QSOLightupParams.

Referenced by get_lambdanet().

Here is the caller graph for this function:

init_qso_lightup()

void init_qso_lightup

(

char *

reion_hist_file

)

Definition at line 238 of file cooling_qso_lightup.c.

{
    if(QSOLightupParams.QSOLightupOn)
        load_heii_reion_hist(reion_hist_file);
}

References load_heii_reion_hist(), qso_lightup_params::QSOLightupOn, and QSOLightupParams.

Referenced by init_cooling().

Here is the call graph for this function:

Here is the caller graph for this function:

ionize_all_part()

static int64_t ionize_all_part ( int  qso_ind, int *  qso_cand, struct QSOPriv  priv, ForceTree *  tree  )

static

Definition at line 473 of file cooling_qso_lightup.c.

{
    /* This treewalk finds not yet ionized particles within the radius of the black hole, ionizes them and
     * adds an instantaneous heating to them. */
    TreeWalk tw[1] = {{0}};
 
    tw->ev_label = "HELIUM";
    /* This could select the black holes to be made quasars, but we do it below.*/
    tw->haswork = NULL;
    tw->tree = tree;
 
    /* We set Hsml to a constant in ngbiter, so this
     * searches a constant distance from the halo.*/
    tw->visit = (TreeWalkVisitFunction) treewalk_visit_ngbiter;
    tw->ngbiter_type_elsize = sizeof(TreeWalkNgbIterBase);
    tw->ngbiter = (TreeWalkNgbIterFunction) ionize_ngbiter;
 
    tw->fill = (TreeWalkFillQueryFunction) ionize_copy;
    tw->reduce = NULL;
    tw->postprocess = NULL;
    tw->query_type_elsize = sizeof(TreeWalkQueryQSOLightup);
    tw->result_type_elsize = sizeof(TreeWalkResultBase);
 
    priv.N_ionized = ta_malloc("n_ionized", int64_t, omp_get_max_threads());
    memset(priv.N_ionized, 0, sizeof(int64_t) * omp_get_max_threads());
    tw->priv = &priv;
 
    /* This runs only on one BH*/
    if(qso_ind > 0)
        treewalk_run(tw, &qso_cand[qso_ind], 1);
    else
        treewalk_run(tw, NULL, 0);
 
    int64_t N_ionized = 0;
    int i;
    for(i = 0; i < omp_get_max_threads(); i++)
        N_ionized += priv.N_ionized[i];
 
    ta_free(priv.N_ionized);
 
    return N_ionized;
}

References TreeWalk::ev_label, TreeWalk::fill, TreeWalk::haswork, ionize_copy(), ionize_ngbiter(), QSOPriv::N_ionized, TreeWalk::ngbiter, TreeWalk::ngbiter_type_elsize, TreeWalk::postprocess, TreeWalk::priv, TreeWalk::query_type_elsize, TreeWalk::reduce, TreeWalk::result_type_elsize, ta_free, ta_malloc, TreeWalk::tree, treewalk_run(), treewalk_visit_ngbiter(), and TreeWalk::visit.

Referenced by turn_on_quasars().

Here is the call graph for this function:

Here is the caller graph for this function:

ionize_copy()

static void ionize_copy ( int  place, TreeWalkQueryQSOLightup *  I, TreeWalk *  tw  )

static

Definition at line 453 of file cooling_qso_lightup.c.

{
    int k;
    FOFGroups * fof = QSO_GET_PRIV(tw)->fof;
    /* Strictly speaking this is inefficient:
     * we are also copying the properties of the *particle*
     * in place in treewalk.c. However, this does not matter unless
     * there are more local groups than particles!*/
    for(k = 0; k < 3; k++)
    {
        I->base.Pos[k] = fof->Group[place].CM[k];
    }
    I->ID = fof->Group[place].base.MinID;
}

References TreeWalkQueryQSOLightup::base, Group::base, Group::CM, FOFGroups::Group, TreeWalkQueryQSOLightup::ID, BaseGroup::MinID, TreeWalkQueryBase::Pos, and QSO_GET_PRIV.

Referenced by ionize_all_part().

Here is the caller graph for this function:

ionize_ngbiter()

static void ionize_ngbiter ( TreeWalkQueryQSOLightup *  I, TreeWalkResultBase *  O, TreeWalkNgbIterBase *  iter, LocalTreeWalk *  lv  )

static

Ionize and heat the particles

Definition at line 419 of file cooling_qso_lightup.c.

{
 
    if(iter->other == -1) {
        /* Gas only ( 1 == 1 << 0, the bit for type 0)*/
        iter->mask = 1;
        /* Bubble size*/
        double bubble = gaussian_rng(QSOLightupParams.mean_bubble, sqrt(QSOLightupParams.var_bubble), I->ID);
        iter->Hsml = bubble;
        /* Don't care about gas HSML */
        iter->symmetric = NGB_TREEFIND_ASYMMETRIC;
        return;
    }
 
    int other = iter->other;
 
    /* Only ionize gas*/
    if(P[other].Type != 0)
        return;
 
    int ionized = ionize_single_particle(other, QSO_GET_PRIV(lv->tw)->a3inv, QSO_GET_PRIV(lv->tw)->uu_in_cgs);
 
    if(!ionized)
        return;
 
    int tid = omp_get_thread_num();
    /* Add to the ionization counter for this thread*/
    QSO_GET_PRIV(lv->tw)->N_ionized[tid] ++;
}

References gaussian_rng(), TreeWalkNgbIterBase::Hsml, TreeWalkQueryQSOLightup::ID, ionize_single_particle(), TreeWalkNgbIterBase::mask, qso_lightup_params::mean_bubble, NGB_TREEFIND_ASYMMETRIC, TreeWalkNgbIterBase::other, P, QSO_GET_PRIV, QSOLightupParams, TreeWalkNgbIterBase::symmetric, LocalTreeWalk::tw, and qso_lightup_params::var_bubble.

Referenced by ionize_all_part().

Here is the call graph for this function:

Here is the caller graph for this function:

ionize_single_particle()

static int ionize_single_particle ( int  other, double  a3inv, double  uu_in_cgs  )

static

Definition at line 380 of file cooling_qso_lightup.c.

{
    /* Mark it ionized if not done so already.*/
    int done;
    #pragma omp atomic capture
    {
        done = P[other].HeIIIionized;
        P[other].HeIIIionized = 1;
    }
    if(done)
        return 0;
 
    /* Heat the particle*/
 
    /* Number of helium atoms per g in the particle*/
    double nheperg = (1 - HYDROGEN_MASSFRAC) / (PROTONMASS * HEMASS);
    /* Total heating per unit mass ergs/g for the particle*/
    double deltau = qso_inst_heating * nheperg;
 
    /* Conversion factor between internal energy and entropy.*/
    double entropytou = pow(SPHP(other).Density * a3inv, GAMMA_MINUS1) / GAMMA_MINUS1;
    /* Convert to entropy in internal units*/
    /* Only one thread may get here*/
    SPHP(other).Entropy += deltau / uu_in_cgs / entropytou;
    return 1;
}

References GAMMA_MINUS1, HEMASS, HYDROGEN_MASSFRAC, P, PROTONMASS, qso_inst_heating, and SPHP.

Referenced by ionize_ngbiter(), and turn_on_quasars().

Here is the caller graph for this function:

load_heii_reion_hist()

static void load_heii_reion_hist ( const char *  reion_hist_file )

static

Definition at line 137 of file cooling_qso_lightup.c.

{
    int ThisTask;
    FILE * fd = NULL;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
 
    message(0, "HeII: Loading HeII reionization history from file: %s\n",reion_hist_file);
    if(ThisTask == 0) {
        fd = fopen(reion_hist_file, "r");
        if(!fd)
            endrun(456, "HeII: Could not open reionization history file at: '%s'\n", reion_hist_file);
 
        /*Find size of file*/
        Nreionhist = 0;
        while(1)
        {
            char buffer[1024];
            char * retval = fgets(buffer, 1024, fd);
            /*Happens on end of file*/
            if(!retval)
                break;
            retval = strtok(buffer, " \t\n");
            /*Discard comments*/
            if(!retval || retval[0] == '#')
                continue;
            Nreionhist++;
        }
        rewind(fd);
        /* Discard first two lines*/
        Nreionhist -=2;
    }
 
    MPI_Bcast(&Nreionhist, 1, MPI_INT, 0, MPI_COMM_WORLD);
 
    if(Nreionhist<= 2)
        endrun(1, "HeII: Reionization history contains: %d entries, not enough.\n", Nreionhist);
 
    /*Allocate memory for the reionization history table.*/
    He_zz = (double *) mymalloc("ReionizationTable", 3 * Nreionhist * sizeof(double));
    XHeIII = He_zz + Nreionhist;
    LMFP = He_zz + 2 * Nreionhist;
 
    if(ThisTask == 0)
    {
        double qso_spectral_index = 0, photon_threshold_energy = 0;
        int prei = 0;
        int i = 0;
        while(i < Nreionhist)
        {
            char buffer[1024];
            char * line = fgets(buffer, 1024, fd);
            /*Happens on end of file*/
            if(!line)
                break;
            char * retval = strtok(line, " \t\n");
            if(!retval || retval[0] == '#')
                continue;
            if(prei == 0)
            {
                qso_spectral_index = atof(retval);
                prei++;
                continue;
            }
            else if(prei == 1)
            {
                photon_threshold_energy = atof(retval);
                prei++;
                continue;
            }
            /* First column: redshift. Convert to scale factor so it is increasing.*/
            He_zz[i] = 1./(1+atof(retval));
            /* Second column: HeIII fraction.*/
            retval = strtok(NULL, " \t");
            if(!retval)
                endrun(12, "HeII: Line %s of reionization table was incomplete!\n", line);
            XHeIII[i] = atof(retval);
            /* Third column: long mean free path photons.*/
            retval = strtok(NULL, " \t");
            if(!retval)
                endrun(12, "HeII: Line %s of reionization table was incomplete!\n", line);
            LMFP[i] = atof(retval);
            i++;
        }
        fclose(fd);
        qso_inst_heating = Q_inst(photon_threshold_energy, qso_spectral_index);
    }
    /*Broadcast data to other processors*/
    MPI_Bcast(He_zz, 3 * Nreionhist, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(&qso_inst_heating, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    /* Initialize the interpolators*/
    HeIII_intp = gsl_interp_alloc(gsl_interp_linear,Nreionhist);
    LMFP_intp = gsl_interp_alloc(gsl_interp_linear,Nreionhist);
    gsl_interp_init(HeIII_intp, He_zz, XHeIII, Nreionhist);
    gsl_interp_init(LMFP_intp, He_zz, LMFP, Nreionhist);
 
    QSOLightupParams.heIIIreion_start = 1/He_zz[0]-1;
 
    message(0, "HeII: Read %d lines z_reion = %g - %g from file %s\n", Nreionhist, 1/He_zz[0] -1, 1/He_zz[Nreionhist-1]-1, reion_hist_file);
}

References endrun(), He_zz, HeIII_intp, qso_lightup_params::heIIIreion_start, LMFP, LMFP_intp, message(), mymalloc, Nreionhist, Q_inst(), qso_inst_heating, QSOLightupParams, ThisTask, and XHeIII.

Referenced by init_qso_lightup().

Here is the call graph for this function:

Here is the caller graph for this function:

need_change_helium_ionization_fraction()

int need_change_helium_ionization_fraction

(

double

atime

)

Definition at line 657 of file cooling_qso_lightup.c.

{
    double desired_ion_frac = gsl_interp_eval(HeIII_intp, He_zz, XHeIII, atime, NULL);
    double curionfrac = gas_ionization_fraction();
    if(curionfrac < desired_ion_frac)
        return 1;
    return 0;
}

References gas_ionization_fraction(), He_zz, HeIII_intp, and XHeIII.

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

Q_inst()

static double Q_inst ( double  Emax, double  alpha_q  )

static

Definition at line 113 of file cooling_qso_lightup.c.

{
    /* Total ionizing flux for the short mean free path photons*/
    double intflux = (pow(Emax,-alpha_q+1.)-pow(E0_HeII,-alpha_q+1.))/(pow(Emax,-alpha_q)- pow(E0_HeII,-alpha_q));
    /* Heating is input per unit mass, so quasars in denser areas will provide more heating.*/
    double Q_inst = (alpha_q/(alpha_q - 1.0))*intflux -E0_HeII;
    return eVinergs * Q_inst;
}

References E0_HeII, and eVinergs.

Referenced by load_heii_reion_hist().

Here is the caller graph for this function:

qso_lightup_on()

int qso_lightup_on

(

void

)

Definition at line 682 of file cooling_qso_lightup.c.

{
    return QSOLightupParams.QSOLightupOn;
}

References qso_lightup_params::QSOLightupOn, and QSOLightupParams.

Referenced by open_outputfiles().

Here is the caller graph for this function:

set_qso_lightup_par()

void set_qso_lightup_par

(

struct qso_lightup_params

qso

)

Definition at line 88 of file cooling_qso_lightup.c.

{
    QSOLightupParams = qso;
}

References QSOLightupParams.

set_qso_lightup_params()

void set_qso_lightup_params

(

ParameterSet *

ps

)

Definition at line 94 of file cooling_qso_lightup.c.

{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0) {
        QSOLightupParams.QSOLightupOn = param_get_int(ps, "QSOLightupOn");
        QSOLightupParams.qso_candidate_max_mass = param_get_double(ps, "QSOMaxMass");
        QSOLightupParams.qso_candidate_min_mass = param_get_double(ps, "QSOMinMass");
        QSOLightupParams.mean_bubble = param_get_double(ps, "QSOMeanBubble");
        QSOLightupParams.var_bubble = param_get_double(ps, "QSOVarBubble");
        QSOLightupParams.heIIIreion_finish_frac = param_get_double(ps, "QSOHeIIIReionFinishFrac");
    }
    MPI_Bcast(&QSOLightupParams, sizeof(struct qso_lightup_params), MPI_BYTE, 0, MPI_COMM_WORLD);
}

References qso_lightup_params::heIIIreion_finish_frac, qso_lightup_params::mean_bubble, param_get_double(), param_get_int(), qso_lightup_params::qso_candidate_max_mass, qso_lightup_params::qso_candidate_min_mass, qso_lightup_params::QSOLightupOn, QSOLightupParams, ThisTask, and qso_lightup_params::var_bubble.

Referenced by read_parameter_file().

Here is the call graph for this function:

Here is the caller graph for this function:

turn_on_quasars()

static void turn_on_quasars ( double  atime, FOFGroups *  fof, DomainDecomp *  ddecomp, Cosmology *  CP, double  uu_in_cgs, FILE *  FdHelium  )

static

Definition at line 520 of file cooling_qso_lightup.c.

{
    int ncand = 0;
    int * qso_cand = NULL;
    int64_t n_gas_tot=0, tot_n_ionized=0, ncand_tot=0;
    MPI_Allreduce(&SlotsManager->info[0].size, &n_gas_tot, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
    double desired_ion_frac = gsl_interp_eval(HeIII_intp, He_zz, XHeIII, atime, NULL);
    struct QSOPriv priv;
    priv.fof = fof;
    priv.uu_in_cgs = uu_in_cgs;
    priv.a3inv = 1/pow(atime, 3);
 
    /* If the desired ionization fraction is above a threshold (by default 0.95)
     * ionize all particles*/
    if(desired_ion_frac > QSOLightupParams.heIIIreion_finish_frac) {
        int i, nionized=0;
        int64_t nion_tot=0;
        #pragma omp parallel for reduction(+: nionized)
        for (i = 0; i < PartManager->NumPart; i++){
            if (P[i].Type == 0)
                nionized += ionize_single_particle(i, priv.a3inv, priv.uu_in_cgs);
        }
        sumup_large_ints(1, &nionized, &nion_tot);
        message(0, "HeII: Helium ionization finished, flash-ionizing %ld particles (%g of total)\n", nion_tot, (double) nion_tot /(double) n_gas_tot);
    }
 
    double rhobar = CP->OmegaBaryon * (3 * HUBBLE * CP->HubbleParam * HUBBLE * CP->HubbleParam)/ (8 * M_PI * GRAVITY) * priv.a3inv;
    double totbubblegasmass = 4 * M_PI / 3. * pow(QSOLightupParams.mean_bubble, 3) * rhobar;
    /* Total expected ionizations if the bubbles do not overlap at all
     * and the bubble is at mean density.*/
    int64_t non_overlapping_bubble_number = n_gas_tot * totbubblegasmass / CP->OmegaBaryon;
    double initionfrac = gas_ionization_fraction();
    double curionfrac = initionfrac;
 
    message(0, "HeII: Started helium reionization model with ionization fraction %g\n", initionfrac);
    if(curionfrac < desired_ion_frac) {
        ncand = build_qso_candidate_list(&qso_cand, fof);
        walltime_measure("/HeIII/Find");
    }
 
    int64_t ncand_before = count_QSO_halos(ncand, &ncand_tot, MPI_COMM_WORLD);
    int iteration;
 
    /* If there are no quasars this will be tough*/
    if(ncand_tot == 0) {
        if(qso_cand)
            myfree(qso_cand);
        return;
    }
    message(0, "HeII: Built quasar candidate list from %d quasars\n", ncand_tot);
    ForceTree gasTree = {0};
    force_tree_rebuild_mask(&gasTree, ddecomp, GASMASK, 0, NULL);
    for(iteration = 0; curionfrac < desired_ion_frac; iteration++){
        /* Get a new quasar*/
        int new_qso = choose_QSO_halo(ncand, &ncand_before, &ncand_tot, fof->TotNgroups+iteration);
        if(new_qso >= ncand)
            endrun(12, "HeII: QSO %d > no. candidates %d! Cannot happen\n", new_qso, ncand);
        /* Make sure someone has a quasar*/
        if(ncand_tot <= 0) {
            if(desired_ion_frac - curionfrac > 0.1)
                message(0, "HeII: Ionization fraction %g less than desired ionization fraction of %g because not enough quasars\n", curionfrac, desired_ion_frac);
            break;
        }
 
        int64_t n_ionized = ionize_all_part(new_qso, qso_cand, priv, &gasTree);
        int64_t tot_qso_ionized = 0;
        /* Check that the ionization fraction changed*/
        MPI_Allreduce(&n_ionized, &tot_qso_ionized, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
        curionfrac += (double) tot_qso_ionized / (double) n_gas_tot;
        tot_n_ionized += tot_qso_ionized;
        /* Get the quasar position on rank 0*/
        double qso_pos[3] = {0};
        if(new_qso > 0) {
            int qplace = qso_cand[new_qso];
            int k;
            for(k = 0; k < 3; k++) {
                qso_pos[k] = fof->Group[qplace].CM[k]- PartManager->CurrentParticleOffset[k];
                while(qso_pos[k] > PartManager->BoxSize) qso_pos[k] -= PartManager->BoxSize;
                while(qso_pos[k] <= 0) qso_pos[k] += PartManager->BoxSize;
            }
            message(1, "HeII: Quasar %d changed the HeIII ionization fraction to %g, ionizing %ld\n", qplace, curionfrac, tot_qso_ionized);
        }
        /* Format: Time = current scale factor,
         * ID of the quasar (the index of the FOF halo)
         * FOF halo position, x,y,z,
         * Current ionized fraction
         * total number of particles ionized by this quasar*/
        MPI_Allreduce(MPI_IN_PLACE, qso_pos, 3, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
        if(FdHelium) {
            fprintf(FdHelium, "%g %g %g %g %g %ld\n", atime, qso_pos[0], qso_pos[1], qso_pos[2], curionfrac, tot_qso_ionized);
            fflush(FdHelium);
        }
 
        /* Break the loop if we do not ionize enough particles this round.
         * Try again next timestep when we will hopefully have new BHs.*/
        if(tot_qso_ionized < 0.01 * non_overlapping_bubble_number && iteration > 10) {
            message(0, "HeII: Stopping ionization at iteration %d because insufficient ionization happened.\n", iteration);
            break;
        }
 
        /* Remove this candidate from the list by moving the list down.*/
        if( new_qso >= 0) {
            memmove(qso_cand+new_qso, qso_cand+new_qso+1, ncand - new_qso+1);
            ncand--;
        }
    }
    force_tree_free(&gasTree);
    if(qso_cand) {
        myfree(qso_cand);
    }
 
    if(tot_n_ionized > 0)
        message(0, "HeII: HeIII fraction from %g -> %g, ionizing %ld. Wanted %g.\n", initionfrac, curionfrac, tot_n_ionized, desired_ion_frac);
    else
        message(0, "HeII: HeIII fraction unchanged at %g. Wanted %g\n", curionfrac, desired_ion_frac);
    walltime_measure("/HeIII/Ionize");
}

References QSOPriv::a3inv, part_manager_type::BoxSize, build_qso_candidate_list(), choose_QSO_halo(), Group::CM, count_QSO_halos(), CP, part_manager_type::CurrentParticleOffset, endrun(), QSOPriv::fof, force_tree_free(), force_tree_rebuild_mask(), gas_ionization_fraction(), GASMASK, GRAVITY, FOFGroups::Group, He_zz, HeIII_intp, qso_lightup_params::heIIIreion_finish_frac, HUBBLE, Cosmology::HubbleParam, slots_manager_type::info, ionize_all_part(), ionize_single_particle(), qso_lightup_params::mean_bubble, message(), MPI_INT64, myfree, part_manager_type::NumPart, Cosmology::OmegaBaryon, P, PartManager, QSOLightupParams, slot_info::size, SlotsManager, sumup_large_ints(), FOFGroups::TotNgroups, QSOPriv::uu_in_cgs, walltime_measure, and XHeIII.

Referenced by do_heiii_reionization().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

He_zz

double* He_zz

static

Definition at line 81 of file cooling_qso_lightup.c.

Referenced by do_heiii_reionization(), during_helium_reionization(), get_long_mean_free_path_heating(), load_heii_reion_hist(), need_change_helium_ionization_fraction(), and turn_on_quasars().

HeIII_intp

gsl_interp* HeIII_intp

static

Definition at line 84 of file cooling_qso_lightup.c.

Referenced by load_heii_reion_hist(), need_change_helium_ionization_fraction(), and turn_on_quasars().

last_long_mfp_heating

double last_long_mfp_heating

static

Definition at line 245 of file cooling_qso_lightup.c.

Referenced by get_long_mean_free_path_heating().

last_zz

double last_zz

static

Definition at line 244 of file cooling_qso_lightup.c.

Referenced by get_long_mean_free_path_heating().

LMFP

double* LMFP

static

Definition at line 83 of file cooling_qso_lightup.c.

Referenced by get_long_mean_free_path_heating(), and load_heii_reion_hist().

LMFP_intp

gsl_interp* LMFP_intp

static

Definition at line 85 of file cooling_qso_lightup.c.

Referenced by get_long_mean_free_path_heating(), and load_heii_reion_hist().

Nreionhist

int Nreionhist

static

Definition at line 80 of file cooling_qso_lightup.c.

Referenced by do_heiii_reionization(), during_helium_reionization(), get_long_mean_free_path_heating(), and load_heii_reion_hist().

qso_inst_heating

double qso_inst_heating

static

Definition at line 79 of file cooling_qso_lightup.c.

Referenced by ionize_single_particle(), and load_heii_reion_hist().

QSOLightupParams

struct qso_lightup_params QSOLightupParams

static

Definition at line 1 of file cooling_qso_lightup.c.

Referenced by build_qso_candidate_list(), do_heiii_reionization(), during_helium_reionization(), get_long_mean_free_path_heating(), init_qso_lightup(), ionize_ngbiter(), load_heii_reion_hist(), qso_lightup_on(), set_qso_lightup_par(), set_qso_lightup_params(), and turn_on_quasars().

XHeIII

double* XHeIII

static

Definition at line 82 of file cooling_qso_lightup.c.

Referenced by load_heii_reion_hist(), need_change_helium_ionization_fraction(), and turn_on_quasars().