FML::COSMOLOGY::SPHERICALCOLLAPSE Namespace Reference

Data Structures
class	SphericalCollapse
	Class for doing general spherical collapse calculations, i.e. More...
class	SphericalCollapseModel
	The class contains all the info about the model needed to run the spherical collapse and halomodel Cosmology/gravity functions provided using lambda functions. More...

Typedefs
using	Timer = FML::UTILS::Timings
using	DVector = FML::SOLVERS::ODESOLVER::DVector
using	DVector2D = FML::SOLVERS::ODESOLVER::DVector2D
using	ODEFunction = FML::SOLVERS::ODESOLVER::ODEFunction
using	ODEFunctionJacobian = FML::SOLVERS::ODESOLVER::ODEFunctionJacobian
using	ODESolver = FML::SOLVERS::ODESOLVER::ODESolver
using	Spline = FML::INTERPOLATION::SPLINE::Spline
using	Func1 = std::function< double(double)>

Functions
void	compute_growthfactor (const SphericalCollapseModel &model, Spline &growthfactor_spline, Spline &growthrate_spline, const double xini=std::log(1e-5), const double xend=0.0, const int npts=1000)
	Computes the growth factor (not normalized) and growth rate for a given model and returns splines with \(D(x)\) and \(f(x)\) where \(x=\log(a)\). More...
double	roverR_of_delta_and_x (double delta, double x, double RH0=1.0, double delta_ini=0.0)
	This is the dimensionless quantity \(r(x)/R\) (the radius of the tophat) for spherical collapse evolution,. More...
double	evolve_spherical_collapse_equations (const SphericalCollapseModel &model, const double delta_ini, DVector &x_array, DVector &delta_array, DVector &delta_prime_array, DVector &delta_lin_array, DVector &delta_lin_prime_array, const bool get_arrays=true, const double xini=std::log(1e-3), const double xcollapse=std::log(1.0), const int npts=1000, const double delta_infinity=1e16)
	Evolve the spherical collapse equation + linear for a given set of IC delta_ini over a given range [xini,xcollapse] and return the arrays + the x-value where delta(x) > delta_infinity (i.e. More...
double	find_delta_ini (const SphericalCollapseModel &model, const double xini=std::log(1e-3), const double xcollapse=std::log(1.0), const int npts=1000, const double epsilon_delta=1e-8)
	Do root finding to figure out the delta_ini that gives collapse at x=xcollapse. More...
auto	compute_times (const double xini, const double xcollapse, const SphericalCollapseModel &model, const Spline &delta_spline, const Spline &delta_prime_spline) -> std::tuple< double, double, double, double >
	Extract useful times from spherical collapse evolution like turnaround, virialization and non-linear onset. More...
void	compute_sphericalcollapse_splines (const SphericalCollapseModel &model, Spline &deltac_of_x_spline, Spline &DeltaVir_of_x_spline, Spline &xnl_of_x_spline, Spline &xta_of_x_spline, Spline &xvir_of_x_spline, Spline &delta_ini_of_x_spline, const double xini=std::log(1e-3), const int npts=1000, bool verbose=false)
	Do spherical collapse for all times. More...

Typedef Documentation

DVector

using FML::COSMOLOGY::SPHERICALCOLLAPSE::DVector = typedef FML::SOLVERS::ODESOLVER::DVector

Definition at line 21 of file SphericalCollapse.h.

DVector2D

using FML::COSMOLOGY::SPHERICALCOLLAPSE::DVector2D = typedef FML::SOLVERS::ODESOLVER::DVector2D

Definition at line 22 of file SphericalCollapse.h.

Func1

using FML::COSMOLOGY::SPHERICALCOLLAPSE::Func1 = typedef std::function<double(double)>

Definition at line 14 of file SphericalCollapseModel.h.

ODEFunction

typedef FML::SOLVERS::ODESOLVER::ODEFunction FML::COSMOLOGY::SPHERICALCOLLAPSE::ODEFunction

Definition at line 23 of file SphericalCollapse.h.

ODEFunctionJacobian

using FML::COSMOLOGY::SPHERICALCOLLAPSE::ODEFunctionJacobian = typedef FML::SOLVERS::ODESOLVER::ODEFunctionJacobian

Definition at line 24 of file SphericalCollapse.h.

ODESolver

using FML::COSMOLOGY::SPHERICALCOLLAPSE::ODESolver = typedef FML::SOLVERS::ODESOLVER::ODESolver

Definition at line 25 of file SphericalCollapse.h.

Spline

using FML::COSMOLOGY::SPHERICALCOLLAPSE::Spline = typedef FML::INTERPOLATION::SPLINE::Spline

Definition at line 26 of file SphericalCollapse.h.

Timer

using FML::COSMOLOGY::SPHERICALCOLLAPSE::Timer = typedef FML::UTILS::Timings

Definition at line 20 of file SphericalCollapse.h.

Function Documentation

compute_growthfactor()

void FML::COSMOLOGY::SPHERICALCOLLAPSE::compute_growthfactor	(	const SphericalCollapseModel &	model,
		Spline &	growthfactor_spline,
		Spline &	growthrate_spline,
		const double	xini = std::log(1e-5),
		const double	xend = 0.0,
		const int	npts = 1000
	)

Computes the growth factor (not normalized) and growth rate for a given model and returns splines with \(D(x)\) and \(f(x)\) where \(x=\log(a)\).

Parameters

[in]	model	The spherical-collapse model
[in,out]	growthfactor_spline	Spline of the growthfactor \(D(\log a)\)
[in,out]	growthrate_spline	Spline of the growth-rate \(d\log D(\log a)/d\log a\)
[in]	xini	The time \(x=\log(a)\) we start the integration. Fiducial \(a=10^{-5}\).
[in]	xend	The time we end the integration.
[in]	npts	Number of points (linear between xini and xend) to compute it on.

Definition at line 12 of file SphericalCollapse.cpp.

                         {
 
        auto x_array = FML::MATH::linspace(xini, xend, npts);
        ODEFunction deriv = [&model](double x, const double * y, double * dydx) {
          return model.ode_system_linear_growth(x, y, dydx);
        };
 
        // EdS D ~ a so D'=D
        DVector y_initial{std::exp(xini), std::exp(xini)};
        ODESolver ode;
        ode.solve(deriv, x_array, y_initial);
        const auto growthfactor_array = ode.get_data_by_component(0);
        const auto delta_prime_lin_array = ode.get_data_by_component(1);
        auto growthrate_array = delta_prime_lin_array;
        for(size_t i = 0; i < growthfactor_array.size(); i++){
          growthrate_array[i] /= growthfactor_array[i];
        }
        growthfactor_of_x_spline = Spline(x_array, growthfactor_array, "Growthfactor D(loga)"); 
        growthrate_of_x_spline = Spline(x_array, growthrate_array, "Growthrate f(loga)"); 
      }

compute_sphericalcollapse_splines()

void FML::COSMOLOGY::SPHERICALCOLLAPSE::compute_sphericalcollapse_splines	(	const SphericalCollapseModel &	model,
		Spline &	deltac_of_x_spline,
		Spline &	DeltaVir_of_x_spline,
		Spline &	xnl_of_x_spline,
		Spline &	xta_of_x_spline,
		Spline &	xvir_of_x_spline,
		Spline &	delta_ini_of_x_spline,
		const double	xini,
		const int	npts,
		bool	verbose
	)

Do spherical collapse for all times.

Definition at line 206 of file SphericalCollapse.cpp.

                       {
 
        // Run spherical collapse and extract splines
        auto run_and_generate_delta_splines = [](
            const double xini,
            const double xcollapse,
            const int npts,
            const SphericalCollapseModel & model,
            Spline & delta_spline, 
            Spline & delta_prime_spline, 
            Spline & delta_lin_spline,
            double delta_ini = -1.0){
          if(delta_ini < 0.0)
            delta_ini = find_delta_ini(model, xini, xcollapse, npts);
          DVector x_array, delta_array, delta_prime_array, delta_lin_array, delta_lin_prime_array;
          evolve_spherical_collapse_equations(
              model, 
              delta_ini,
              x_array,
              delta_array, 
              delta_prime_array, 
              delta_lin_array, 
              delta_lin_prime_array, 
              true,
              xini, 
              xcollapse, 
              npts);
          delta_spline = Spline(x_array, delta_array, "delta(loga)");
          delta_prime_spline = Spline(x_array, delta_prime_array, "ddelta/dloga(loga)");
          delta_lin_spline = Spline(x_array, delta_lin_array, "deltaLin(loga)");
          return delta_ini;
        };
 
        // Range we want deltac(z) etc on
        const double zini = 10.0;
        const double zend = 0.0;
        const int npts_z = 100;
        const auto xcollapse_array = FML::MATH::linspace(std::log(1.0/(1.0+zini)),std::log(1.0/(1.0+zend)),npts_z);
 
        DVector deltac_array(xcollapse_array.size());
        DVector DeltaVir_array(xcollapse_array.size());
        DVector xta_array(xcollapse_array.size());
        DVector xvir_array(xcollapse_array.size());
        DVector xnl_array(xcollapse_array.size());
        DVector delta_ini_array(xcollapse_array.size());
 
#ifdef USE_OMP
#pragma omp parallel for schedule(dynamic,1)
#endif
        for(size_t i = 0; i < xcollapse_array.size(); i++){
          auto xcollapse = xcollapse_array[i];
 
          // Generate solution
          Spline delta_spline;
          Spline delta_prime_spline;
          Spline delta_lin_spline;
          const double delta_ini = run_and_generate_delta_splines(
              xini, 
              xcollapse, 
              npts, 
              model, 
              delta_spline, 
              delta_prime_spline, 
              delta_lin_spline);
 
          // Compute quantities
          const auto times = compute_times(xini, xcollapse, model, delta_spline, delta_prime_spline);
          const double xta  = std::get<0>(times);
          const double xvir = std::get<1>(times);
          const double xnl  = std::get<2>(times);
          const double xta2  = std::get<3>(times);
          const double DeltaVir = (1.0+delta_spline(xvir))*std::exp(3.0*xcollapse-3.0*xvir);
          const double deltac = delta_lin_spline(xcollapse);
 
          if(xcollapse == 0.0 and verbose){
            std::cout << "\n#=====================================================\n";
            std::cout << "# For collapse at z=0 we get:\n";
            std::cout << "#=====================================================\n";
            std::cout << " deltac       = " << deltac << "\n";
            std::cout << " DeltaVir     = " << DeltaVir << "\n";
            std::cout << " Turnaround r = " << roverR_of_delta_and_x(delta_spline(xta),xta,1.0,0.0) << "\n";
            std::cout << " Virial     r = " << roverR_of_delta_and_x(delta_spline(xvir),xvir,1.0,0.0) << "\n";
            std::cout << " rta/2        = " << roverR_of_delta_and_x(delta_spline(xta2),xta2,1.0,0.0) << "\n";
            std::cout << "#=====================================================\n";
          }
 
          xta_array[i] = xta;
          xnl_array[i] = xnl;
          xvir_array[i] = xvir;
          deltac_array[i] = deltac;
          DeltaVir_array[i] = DeltaVir;
          delta_ini_array[i] = delta_ini;
 
          // In this case the problem is not well posed so return nan
          if(model.muofx(xcollapse) <= 0.0){
            throw std::runtime_error("Spherical collapse. Negative Geff/G encountere. Problem not well posed!");
          }
        }
        xta_of_x_spline = Spline(xcollapse_array, xta_array, "loga_turnaround(loga_collapse)");
        xnl_of_x_spline = Spline(xcollapse_array, xnl_array, "loga_nonlinear(loga_collapse)");
        xvir_of_x_spline = Spline(xcollapse_array, xvir_array, "loga_virialization(loga_collapse)");
        deltac_of_x_spline = Spline(xcollapse_array, deltac_array, "deltac(loga_collapse)");
        DeltaVir_of_x_spline = Spline(xcollapse_array, DeltaVir_array, "DeltaVir(loga_collapse)");
        delta_ini_of_x_spline = Spline(xcollapse_array, delta_ini_array, "delta_ini(loga_collapse)");
      }

compute_times()

auto FML::COSMOLOGY::SPHERICALCOLLAPSE::compute_times	(	const double	xini,
		const double	xcollapse,
		const SphericalCollapseModel &	model,
		const Spline &	delta_spline,
		const Spline &	delta_prime_spline
	)		-> std::tuple< double, double, double, double >

Extract useful times from spherical collapse evolution like turnaround, virialization and non-linear onset.

Definition at line 147 of file SphericalCollapse.cpp.

                                                                                         {
 
        // Accuracy in x for finding times
        const double epsilon_x = 1e-8;
 
        // Compute turnaround time
        auto turnaround_condition = [&delta_spline, &delta_prime_spline](double x){
          const double delta = delta_spline(x), delta_prime = delta_prime_spline(x);
          return 1.0 - delta_prime/3.0/(1.0+delta);
        };
        const double xturnaround = FML::MATH::find_root_bisection(turnaround_condition, {xini, xcollapse}, epsilon_x);
 
        // Compute virialization time
        auto virial_condition = [&model, &delta_spline, &delta_prime_spline](double x) {
          const double delta = delta_spline(x), delta_prime = delta_prime_spline(x);
          const double roverR = roverR_of_delta_and_x(delta, x);
          return model.virial_condition(delta, delta_prime, roverR, x);
        };
        // If we have several times where 2T+W=0 pick the one with largest redshift
        double x_search = xturnaround;
        const double deltax = 0.01;
        const int sign = virial_condition(x_search) > 0 ? 1 : -1;
        while(true){
          const int sign1 = virial_condition(x_search) > 0 ? 1 : -1;
          if(sign * sign1 == -1) break;
          x_search += deltax;
        }
        const double xvir = FML::MATH::find_root_bisection(virial_condition, {xturnaround, x_search}, epsilon_x);
 
        // Compute non-linear time
        auto nonlinear_condition = [&model, &delta_spline, &delta_prime_spline](double x){
          const double delta = delta_spline(x), delta_prime = delta_prime_spline(x);
          const double roverR = roverR_of_delta_and_x(delta, x);
          return model.nonlinear_condition(delta, delta_prime, roverR, x);;
        };
        const double xnl = FML::MATH::find_root_bisection(nonlinear_condition, {xini, xcollapse}, epsilon_x);
 
        // Compute time for which r=rturnaround/2
        const double roverR_turnaround = roverR_of_delta_and_x(delta_spline(xturnaround), xturnaround);
        auto halfturnaround_condition = [&roverR_turnaround, &delta_spline](double x) {
          const double delta = delta_spline(x);
          const double roverR = roverR_of_delta_and_x(delta, x);
          return roverR - roverR_turnaround/2.0;
        };
        const double xhalfturnaround = FML::MATH::find_root_bisection(halfturnaround_condition, {xturnaround, xcollapse}, epsilon_x);
 
        return std::tuple(xturnaround, xvir, xnl, xhalfturnaround);
      }

evolve_spherical_collapse_equations()

double FML::COSMOLOGY::SPHERICALCOLLAPSE::evolve_spherical_collapse_equations	(	const SphericalCollapseModel &	model,
		const double	delta_ini,
		DVector &	x_array,
		DVector &	delta_array,
		DVector &	delta_prime_array,
		DVector &	delta_lin_array,
		DVector &	delta_lin_prime_array,
		const bool	get_arrays = true,
		const double	xini = std::log(1e-3),
		const double	xcollapse = std::log(1.0),
		const int	npts = 1000,
		const double	delta_infinity = 1e16
	)

Evolve the spherical collapse equation + linear for a given set of IC delta_ini over a given range [xini,xcollapse] and return the arrays + the x-value where delta(x) > delta_infinity (i.e.

it has collapsed completely). If it haven't collapsed return xcollapse This value can be used to more easily figure out delta_ini that corresponds to collapse at exactly xcollapse via root finding

Definition at line 61 of file SphericalCollapse.cpp.

                                       {
 
        // ODE solver parameters
        const double hstart = 1e-3;
        const double epsilon_abs = 1e-6;
        const double epsilon_rel = 1e-6;
 
        // The first x for which delta > delta_infinity
        double xminover = xcollapse;
 
        // Set x array (linear in 'a' seems to work well)
        x_array = FML::MATH::linspace(std::exp(xini),std::exp(xcollapse),npts);
        for(auto & x : x_array) x = std::log(x); 
 
        ODEFunction deriv = [&model,&xminover,&delta_infinity](double x, const double * y, double * dydx) -> int {
          // Check if we have collapsed
          if(y[0] > delta_infinity){
            if(x < xminover) xminover = x;
            dydx[0] = dydx[1] = 0.0;
            return GSL_SUCCESS;
          }
          return model.ode_system_non_linear_growth(x, y, dydx);
        };
 
        // The growing mode has D ~ a = exp(x) so D' = D
        DVector y_initial{delta_ini, delta_ini, delta_ini, delta_ini};
        ODESolver ode;
        ode.set_accuracy(hstart, epsilon_abs, epsilon_rel);
        ode.solve(deriv, x_array, y_initial, gsl_odeiv2_step_rkck);
 
        // Fetch solution arrays
        if(get_arrays){
          delta_array = ode.get_data_by_component(0);
          delta_prime_array = ode.get_data_by_component(1);
          delta_lin_array = ode.get_data_by_component(2);
          delta_lin_prime_array = ode.get_data_by_component(3);
        }
        return xminover;
      }

find_delta_ini()

double FML::COSMOLOGY::SPHERICALCOLLAPSE::find_delta_ini	(	const SphericalCollapseModel &	model,
		const double	xini,
		const double	xcollapse,
		const int	npts,
		const double	epsilon_delta
	)

Do root finding to figure out the delta_ini that gives collapse at x=xcollapse.

Definition at line 117 of file SphericalCollapse.cpp.

                                     {
 
        auto condition = [&](double delta_ini){
          DVector x_array, delta_array, delta_prime_array, delta_lin_array, delta_lin_prime_array;
          const double xminover = evolve_spherical_collapse_equations(
              model, 
              delta_ini, 
              x_array,
              delta_array, 
              delta_prime_array, 
              delta_lin_array, 
              delta_lin_prime_array, 
              false,
              xini, 
              xcollapse, 
              npts);
          return xminover - xcollapse;
        };
        return FML::MATH::find_root_bisection(condition, {0.0, 1.0}, epsilon_delta);
      }

roverR_of_delta_and_x()

double FML::COSMOLOGY::SPHERICALCOLLAPSE::roverR_of_delta_and_x	(	double	delta,
		double	x,
		double	RH0,
		double	delta_ini
	)

This is the dimensionless quantity \(r(x)/R\) (the radius of the tophat) for spherical collapse evolution,.

Definition at line 43 of file SphericalCollapse.cpp.

                            { 
        return (RH0) * std::pow((1.0+delta_ini)/(1.0+delta),1.0/3.0)*std::exp(x); 
      };