: class to solve chameleon equations of motion More...

Inheritance diagram for anonymous_namespace{chameleon.cpp}::ChiSolver< T >:

Collaboration diagram for anonymous_namespace{chameleon.cpp}::ChiSolver< T >:

Public Member Functions
	ChiSolver (size_t N, unsigned int Nmin, const Sim_Param &sim, bool verbose=true)

	ChiSolver (size_t N, const Sim_Param &sim, bool verbose=true)

T	chi_a (T a) const

T	chi_force_units (T a) const

void	get_chi_k (Mesh &rho_k)

bool	check_surr_dens (T const *const rho_grid, std::vector< size_t > index_list, size_t i, size_t N)

void	set_time (T a, const Cosmo_Param &cosmo)

T	get_chi_prefactor () const

T	get_phi_prefactor () const

T	l_operator (const T chi_i, const size_t level, const std::vector< size_t > &index_list, const bool addsource, const T h) const

T	l_operator (const size_t level, const std::vector< size_t > &index_list, const bool addsource, const T h) const override

T	dl_operator (const size_t level, const std::vector< size_t > &index_list, const T h) const override

bool	find_opposite_l_sign (const T f1, const T l1, T df, T &f2, T &l2, const size_t level, const std::vector< size_t > &index_list, const T h) const

bool	check_bisection_convergence (const T df_new, const T l_new) const

T	bisection_step (T &f1, T &l1, T &f2, T &l2, const size_t level, const std::vector< size_t > &index_list, const T h) const

T	bisection (T f1, T l1, const T df, const size_t level, const std::vector< size_t > &index_list, const T h) const

T	upd_operator (const T f, const size_t level, const std::vector< size_t > &index_list, const T h) const override

void	correct_sol (Grid< 3, T > &f, const Grid< 3, T > &corr, const size_t level) override

ES	check_convergence () override
	check if solution already converged More...

void	check_solution (size_t level, Grid< 3, T > &chi) override

void	set_convergence (double eps, double err_stop, double err_stop_min, double rms_stop_min, size_t num_fail)

void	set_bisection_convergence (size_t max_bi_step, T dchi_stop, T l_stop)

void	set_def_convergence ()

void	set_bulk_field ()

void	set_linear_sol_at_level (Mesh &rho, const FFTW_PLAN_TYPE &p_F, const FFTW_PLAN_TYPE &p_B, size_t level)

void	set_linear_recursively (size_t level)

void	set_linear (Mesh &rho, const FFTW_PLAN_TYPE &p_F, const FFTW_PLAN_TYPE &p_B)

void	set_screened (size_t level=0)

T	chi_min (T delta) const

Public Member Functions inherited from MultiGridSolver< 3, T >
	MultiGridSolver ()

	MultiGridSolver (size_t N)

	MultiGridSolver (size_t N, bool verbose)

	MultiGridSolver (size_t N, size_t Nmin, bool verbose)

size_t	get_istep () const

T *	get_y (size_t level=0)

T const *	get_y (size_t level=0) const

Grid< NDIM, T > &	get_grid (size_t level=0)

const Grid< NDIM, T > &	get_grid (size_t level=0) const

MultiGrid< NDIM, T > &	get_mlt_grid (size_t level=0)

const MultiGrid< NDIM, T > &	get_mlt_grid (size_t level=0) const

T *	get_external_field (size_t level, size_t field)

T const *	get_external_field (size_t level, size_t field) const

Grid< NDIM, T > &	get_external_grid (size_t level, size_t field)

const Grid< NDIM, T > &	get_external_grid (size_t level, size_t field) const

size_t	get_external_field_size () const

T	get_multigrid_source (size_t level, size_t i) const

void	set_epsilon (double eps_converge)

void	set_maxsteps (size_t maxsteps)

void	set_ngs_sweeps (size_t ngs_fine, size_t ngs_coarse)

void	set_convergence_criterion_residual (bool use_residual)

void	set_Nlevel (size_t N)

size_t	get_N (size_t level=0) const

size_t	get_Ntot (size_t level=0) const

size_t	get_Nlevel () const

void	add_external_grid (MultiGrid< NDIM, T > *field)

void	set_initial_guess (T guess)

void	set_initial_guess (T *guess)

void	set_initial_guess (Grid< NDIM, T > &guess)

Exit_Status	solve ()

void	clear ()

virtual void	correct_sol (Grid< NDIM, T > &f, const Grid< NDIM, T > &corr, const size_t level)

virtual void	check_solution (size_t level, Grid< NDIM, T > &sol)

void	check_solution (size_t level)

Public Attributes
const T	n
	Hu-Sawicki paramater. More...

const T	beta
	Chameleon coupling constant. More...

const T	chi_0
	2betaMpl*phi_scr More...

const T	phi_prefactor
	prefactor for Poisson equation for gravitational potential More...

const T	chi_prefactor_0
	dimensionless prefactor to poisson equation at a = 1 More...

T	chi_prefactor
	time-dependent prefactor More...

size_t	m_conv_stop = 0

double	m_rms_stop_min

double	m_err_stop

double	m_err_stop_min

size_t	m_num_fail

size_t	m_max_bisection_steps

T	m_dchi_stop

T	m_l_stop

std::vector< std::map< size_t, T > >	fix_vals
	<index, value> More...

Public Attributes inherited from MultiGridSolver< 3, T >
bool	_store_all_residual

std::vector< double >	_res_domain_array

Additional Inherited Members
Public Types inherited from MultiGridSolver< 3, T >
enum	Exit_Status

Protected Member Functions inherited from MultiGridSolver< 3, T >
void	get_neighbor_gridindex (std::vector< size_t > &index_list, size_t i, size_t ngrid)

Protected Attributes inherited from MultiGridSolver< 3, T >
bool	_verbose

bool	_conv_criterion_residual

double	_eps_converge

size_t	_maxsteps

size_t	_istep_vcycle

double	_rms_res

double	_rms_res_i

double	_rms_res_old

Detailed Description

template<typename T>
class anonymous_namespace{chameleon.cpp}::ChiSolver< T >

: class to solve chameleon equations of motion

Definition at line 156 of file chameleon.cpp.

Constructor & Destructor Documentation

template<typename T >

anonymous_namespace{chameleon.cpp}::ChiSolver< T >::ChiSolver	(	size_t	N,
		unsigned int	Nmin,
		const Sim_Param &	sim,
		bool	verbose = `true`
	)

inline

Definition at line 183 of file chameleon.cpp.

                                                                                       :
         MultiGridSolver<3, T>(N, Nmin, verbose), n(sim.chi_opt.n), beta(sim.chi_opt.beta), chi_0(2*beta*MPL*sim.chi_opt.phi),
         phi_prefactor( // prefactor for Poisson equation for gravitational potential, [] = (h/Mpc)^2
             // 4*pi*G*rho_m,0 + computing units [Mpc/h]
             3/2.*sim.cosmo.Omega_m*pow(sim.cosmo.H0 * sim.cosmo.h / c_kms ,2)
         ),
         chi_prefactor_0( // dimensionless prefactor to poisson equation
             // beta*rho_m,0 / (Mpl*chi_0) + computing units [Mpc/h]; additional a^(-3 -3/(1-n)) at each timestep
             phi_prefactor*pow(sim.box_opt.box_size ,2) / sim.chi_opt.phi
         ),
         fix_vals(this->get_Nlevel())
         {
             if ((n <= 0) || (n >= 1) || (chi_0 <= 0)) throw std::out_of_range("invalid values of chameleon power-law potential parameters");
         }

template<typename T >

anonymous_namespace{chameleon.cpp}::ChiSolver< T >::ChiSolver	(	size_t	N,
		const Sim_Param &	sim,
		bool	verbose = `true`
	)

inline

Definition at line 198 of file chameleon.cpp.

198 : ChiSolver(N, 2, sim, verbose) {}

anonymous_namespace{chameleon.cpp}::ChiSolver::ChiSolver

ChiSolver(size_t N, unsigned int Nmin, const Sim_Param &sim, bool verbose=true)

Definition: chameleon.cpp:183

Member Function Documentation

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::bisection	(	T	f1,
		T	l1,
		const T	df,
		const size_t	level,
		const std::vector< size_t > &	index_list,
		const T	h
	)		const

inline

Definition at line 356 of file chameleon.cpp.

References growth_allz::T.

     {/* initialize bisection solver -- find two values with opposite value of l_operator -- and start iterating */
         T f_new, f2, l2;
 
         // get initial guess -- return when failed
         if (!find_opposite_l_sign(f1, l1, df, f2, l2, level, index_list, h)) return f2;
 
         // iterate
         for (size_t i = 0; i < m_max_bisection_steps; ++i){
             f_new = bisection_step(f1, l1, f2, l2, level, index_list, h);
             if (f_new != CHI_MIN)
             {
                 return f_new;
             }
         }
 
         // failed iteration
         return f1;
     }

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::bisection_step	(	T &	f1,
		T &	l1,
		T &	f2,
		T &	l2,
		const size_t	level,
		const std::vector< size_t > &	index_list,
		const T	h
	)		const

inline

Definition at line 333 of file chameleon.cpp.

References anonymous_namespace{chameleon.cpp}::CHI_MIN, and growth_allz::T.

     {/* given 'f1' and 'f2' with different signs of l_operator(f_i) perform one step of bisection:
         f_new = (f1 + f2) / 2
         change whichever l_operator(f_i) has the same sign as l_operator(f_new)
         return value indicates convergence -- 0 (unphysical for chameleon) not, otherwise yes*/
 
         const T f_new = (f1 + f2) / 2;
         const T l_new = l_operator(f_new, level, index_list, true, h);
 
         if (check_bisection_convergence(f2 - f_new, l_new)) return f_new;
 
         if (l1*l_new > 0) {
             f1 = f_new;
             l1 = l_new;
             
         } else {
             f2 = f_new;
             l2 = l_new;
         }
 
         return CHI_MIN;
     }

template<typename T >

bool anonymous_namespace{chameleon.cpp}::ChiSolver< T >::check_bisection_convergence	(	const T	df_new,
		const T	l_new
	)		const

inline

Definition at line 328 of file chameleon.cpp.

     {
         return ((std::abs(l_new) < m_l_stop) || (std::abs(df_new) < m_dchi_stop));
     }

template<typename T >

ES anonymous_namespace{chameleon.cpp}::ChiSolver< T >::check_convergence ( )

inlineoverridevirtual

check if solution already converged

Returns: ES exit status (SUCCESS, FAILURE, SLOW, ITERATE, MAX_STEPS)

bring variables from this-> namespace here
Compute ratio of residual to previous residual
convergence criteria
print some information
Print out some information
check individual conditions

Reimplemented from MultiGridSolver< 3, T >.

Definition at line 422 of file chameleon.cpp.

     {
         const double _rms_res_i = this->_rms_res_i;
         const double _rms_res = this->_rms_res;
         const double _rms_res_old = this->_rms_res_old;
         const double _verbose = this->_verbose;
         const double _istep_vcycle = this->_istep_vcycle;
         const double _eps_converge = this->_eps_converge;
         const double _maxsteps = this->_maxsteps;
 
         const double err = _rms_res_old != 0.0 ? _rms_res/_rms_res_old : 0.0;
         ES converged = ES::ITERATE;
 
         bool res_below_optimal_tresh = (_rms_res < _eps_converge) || (_eps_converge == 0.0);
         bool res_below_minimal_tresh = (_rms_res < m_rms_stop_min) || (m_rms_stop_min == 0.0);
         bool fast_convergence = err < m_err_stop_min;
         bool slow_convergence = err > m_err_stop;
         bool worse_solution = err > 1;
         bool over_max_steps = _istep_vcycle >= _maxsteps;
 
         auto print_success = [=](size_t m_conv_stop){
             BOOST_LOG_TRIVIAL(debug) << "\n\tSUCCESS: res = " << _rms_res << ", err = " << err << ", num_err = " << m_conv_stop << " (istep = " << _istep_vcycle << ")";
             };
         auto print_failure = [=](size_t m_conv_stop){
             BOOST_LOG_TRIVIAL(debug) << "\tFAILURE: res = " << _rms_res << ", err = " << err << ", num_err = " << m_conv_stop << " (istep = " << _istep_vcycle << ")";
             };
 
         auto print_iterate = [=](size_t m_conv_stop){
             BOOST_LOG_TRIVIAL(debug) << "\tITERATE: res = " << _rms_res << ", err = " << err << ", num_err = " << m_conv_stop << " (istep = " << _istep_vcycle << ")";
             };
 
         if (_verbose){
             BOOST_LOG_TRIVIAL(debug) << "    Checking for convergence at step = " << _istep_vcycle;
             BOOST_LOG_TRIVIAL(debug) << "        Residual = " << _rms_res << "  Residual_old = " <<  _rms_res_old;
             BOOST_LOG_TRIVIAL(debug) << "        Residual_i = " << _rms_res_i << "  Err = " << err;
         }
 
         if (res_below_optimal_tresh && !fast_convergence){
             print_success(m_conv_stop);
             converged = ES::SUCCESS;
         }
         else if (worse_solution)
         {
             if(_verbose) print_failure(m_conv_stop);
             m_conv_stop++;
             converged = ES::ITERATE;
         }
         else if (slow_convergence)
         {
             if (res_below_minimal_tresh)
             {
                 print_success(m_conv_stop);
                 converged = ES::SLOW;
             }
             else
             {
                 m_conv_stop++;
                 if (m_conv_stop >= m_num_fail)
                 {
                     print_failure(m_conv_stop);
                     converged = ES::FAILURE;
                 }
                 else
                 {
                      if(_verbose) print_iterate(m_conv_stop);
                     converged = ES::ITERATE;
                 }
             }
         }
         else
         {
             if (m_conv_stop >= m_num_fail)
             {
                 print_failure(m_conv_stop);
                 converged = ES::FAILURE;
             }
             else
             {
                 if(_verbose) print_iterate(m_conv_stop);
                converged = ES::ITERATE;
             }
         }
 
         if(over_max_steps)
         {
             print_failure(m_conv_stop);
             converged = ES::MAX_STEPS;
         }
         
         if (converged != ES::ITERATE)
         {
             m_conv_stop = 0;
         }
         return converged;
     }

template<typename T >

void anonymous_namespace{chameleon.cpp}::ChiSolver< T >::check_solution	(	size_t	level,
		Grid< 3, T > &	chi
	)

inlineoverride

Definition at line 524 of file chameleon.cpp.

     {
         for (auto fv : fix_vals[level]) chi[fv.first] = fv.second;
     }

template<typename T >

bool anonymous_namespace{chameleon.cpp}::ChiSolver< T >::check_surr_dens	(	T const *const	rho_grid,
		std::vector< size_t >	index_list,
		size_t	i,
		size_t	N
	)

inline

Definition at line 237 of file chameleon.cpp.

References anonymous_namespace{test_chameleon.cpp}::get_neighbor_gridindex().

     {/* internal method for finding highest density in nearby points */
         // never fix bulk field in under-dense region
         if (rho_grid[i] <= 0) return false;
 
         // check surrounding points if theres is higher density
         this->get_neighbor_gridindex(index_list, i, N);
         for(size_t i_s : index_list) if (rho_grid[i_s] > rho_grid[i]) return false;
 
         // if current point has the highest density, fix chameleon value
         return true;
     }

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::chi_a ( T a ) const

inline

Definition at line 200 of file chameleon.cpp.

References pow().

     {
         return chi_0*pow(a, 3/(1-n));
     }

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::chi_force_units ( T a ) const

inline

Definition at line 205 of file chameleon.cpp.

     {
         return beta/MPL*chi_a(a)/phi_prefactor;
     }

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::chi_min ( T delta ) const

inline

Definition at line 664 of file chameleon.cpp.

References anonymous_namespace{chameleon.cpp}::CHI_MIN, and pow().

     {/* get chi_bulk for given overdensity */
         if (delta > -1) return std::pow(1+delta, 1/(n-1)) + CHI_MIN;
         else if (delta == -1) return 1 + CHI_MIN;
         else throw std::range_error("invalid values of density: " + std::to_string(delta) + " < -1)");        
     }

template<typename T >

void anonymous_namespace{chameleon.cpp}::ChiSolver< T >::correct_sol	(	Grid< 3, T > &	f,
		const Grid< 3, T > &	corr,
		const size_t	level
	)

inlineoverride

Definition at line 391 of file chameleon.cpp.

References anonymous_namespace{test_chameleon.cpp}::get_neighbor_gridindex(), Grid< NDIM, T >::get_Ntot(), ccl_test_distances::h, cl_cmbl_bm::l, and growth_allz::T.

     {/* Method for correcting solution when going up,
         check for unphysical values */
 
         const size_t Ntot  = f.get_Ntot();
         const T * const rho = this->get_external_field(level, 0);
 
         // bisection variables
         const size_t N = this->get_N(level);
         const T h = 1.0/T( N );
         std::vector<size_t> index_list;
         
         #pragma omp parallel for private(index_list)
         for(size_t i = 0; i < Ntot; i++)
         {
             // do not change values in screened regions
             if (rho[i] == MARK_CHI_BOUND_COND) continue;
             else if (std::abs(corr[i]/(f[i] - CHI_MIN)) < SWITCH_BIS_NEW) f[i] += corr[i];
             else{
                 this->get_neighbor_gridindex(index_list, i, N);
                 T l =  l_operator(level, index_list, true, h);
                 f[i] = bisection(f[i], l, corr[i] / 2, level, index_list, h);
             }   
         }
     }

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::dl_operator	(	const size_t	level,
		const std::vector< size_t > &	index_list,
		const T	h
	)		const

inlineoverridevirtual

Reimplemented from MultiGridSolver< 3, T >.

Definition at line 293 of file chameleon.cpp.

References ccl_test_distances::h, pow(), and growth_allz::T.

     {
         // solution
         const T chi = this->get_y(level)[ index_list[0] ];
 
         // Derivative of source
         const T dsource = chi_prefactor*(1-n)*pow(chi - CHI_MIN, n-2);
 
         // Derivative of kinetic term
         const T dkinetic = -2.0*3;
 
         return dkinetic/(h*h) - dsource;
     }

template<typename T >

bool anonymous_namespace{chameleon.cpp}::ChiSolver< T >::find_opposite_l_sign	(	const T	f1,
		const T	l1,
		T	df,
		T &	f2,
		T &	l2,
		const size_t	level,
		const std::vector< size_t > &	index_list,
		const T	h
	)		const

inline

Definition at line 307 of file chameleon.cpp.

References anonymous_namespace{chameleon.cpp}::CONVERGENCE_BI_STEPS_INIT, and growth_allz::T.

     {/* find such 'f2' that 'l_operator(f2)' has opposite sign than l1
         use df as a guess in which direction to be looking */
         f2 = f1;
         for (size_t j = 0; j < CONVERGENCE_BI_STEPS_INIT; ++j)
         {
             f2 += df;
             if (f2 <= CHI_MIN)
             {
                 const size_t i = index_list[0];
                 const T rho = this->get_external_field(level, 0)[i];
 
                 f2 = chi_min(rho);
                 j = CONVERGENCE_BI_STEPS_INIT;//< end of loop but first check for an improvement
             }
             l2 = l_operator(f2, level, index_list, true, h);
             if (l1*l2 <= 0) return true;
         }
         return false;
     }

template<typename T >

void anonymous_namespace{chameleon.cpp}::ChiSolver< T >::get_chi_k ( Mesh & rho_k )

inline

Definition at line 210 of file chameleon.cpp.

References get_k_sq(), Mesh_base< T >::length, Mesh::N, PI, pow(), and growth_allz::T.

     {/* transform input density in k-space into linear prediction for chameleon field,
         includes w(k) corrections for interpolation of particles */
         const size_t N = rho_k.N;
         const size_t l_half = rho_k.length/2;
         const T mass_sq = (1-n)*chi_prefactor/pow(2*PI, 2); // dimensionless square mass, with derivative factor k* = 2*PI / L
         const T chi_a_n = -1/(1-n); // prefactor for chi(k), in chi_a units
         
         T k2, g_k;
 
         #pragma omp parallel for private(k2, g_k)
         for(size_t i=0; i < l_half;i++){
             k2 = get_k_sq(N, i);
             if (k2 == 0)
             {
                 rho_k[2*i] = 0;
                 rho_k[2*i+1] = 0;
             }
             else
             {
                 g_k = chi_a_n/(k2+mass_sq)*mass_sq; // Green function
                 rho_k[2*i] *= g_k;
                 rho_k[2*i+1] *= g_k;
             }
         }
     }

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::get_chi_prefactor ( ) const

inline

Definition at line 255 of file chameleon.cpp.

255 { return chi_prefactor; }

anonymous_namespace{chameleon.cpp}::ChiSolver::chi_prefactor

T chi_prefactor

time-dependent prefactor

Definition: chameleon.cpp:165

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::get_phi_prefactor ( ) const

inline

Definition at line 256 of file chameleon.cpp.

256 { return phi_prefactor; }

anonymous_namespace{chameleon.cpp}::ChiSolver::phi_prefactor

const T phi_prefactor

prefactor for Poisson equation for gravitational potential

Definition: chameleon.cpp:163

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::l_operator	(	const T	chi_i,
		const size_t	level,
		const std::vector< size_t > &	index_list,
		const bool	addsource,
		const T	h
	)		const

inline

Definition at line 258 of file chameleon.cpp.

References ccl_test_distances::h, pow(), and growth_allz::T.

     {/* The dicretized equation L(phi) */
         // Solution and pde-source grid at current level
         const size_t i = index_list[0];
         T const* const chi = this->get_y(level); // solution
         const T rho = this->get_external_field(level, 0)[i];
         
         // do not change values in screened regions
         if (rho == MARK_CHI_BOUND_COND) return 0;
 
         // The right hand side of the PDE 
         T source = (1 + rho - pow(chi_i - CHI_MIN, n - 1)) * chi_prefactor;
 
         // Compute the standard kinetic term [D^2 phi] (in 1D this is phi''_i =  phi_{i+1} + phi_{i-1} - 2 phi_{i} )
         T kinetic = -2*chi_i*3; // chi, '-2*3' is factor in 3D discrete laplacian
         // go through all surrounding points
         for(auto it = index_list.begin() + 1; it < index_list.end(); ++it) kinetic += chi[*it];
 
         // source term arising rom restricting the equation down to the lower level
         if( level > 0 && addsource ) source += this->get_multigrid_source(level, i);
 
         // The discretized equation of motion L_{ijk...}(phi) = 0
         return kinetic/(h*h) - source;
     }

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::l_operator	(	const size_t	level,
		const std::vector< size_t > &	index_list,
		const bool	addsource,
		const T	h
	)		const

inlineoverridevirtual

Reimplemented from MultiGridSolver< 3, T >.

Definition at line 283 of file chameleon.cpp.

References growth_allz::T.

     {/* The dicretized equation L(phi) */
         // Solution and pde-source grid at current level
         const size_t i = index_list[0];
         T const* const chi = this->get_y(level); // solution
         const T chi_i = chi[i];
         return l_operator(chi_i, level, index_list, addsource, h);
     }

template<typename T >

void anonymous_namespace{chameleon.cpp}::ChiSolver< T >::set_bisection_convergence	(	size_t	max_bi_step,
		T	dchi_stop,
		T	l_stop
	)

inline

Definition at line 538 of file chameleon.cpp.

     {
         m_max_bisection_steps = max_bi_step;
         m_dchi_stop = dchi_stop;
         m_l_stop = l_stop;
     }

template<typename T >

void anonymous_namespace{chameleon.cpp}::ChiSolver< T >::set_bulk_field ( )

inline

Definition at line 553 of file chameleon.cpp.

References growth_allz::T.

     {
         if (!chi_prefactor) throw std::out_of_range("invalid value of scale factor");
         if (!this->get_external_field_size()) throw std::out_of_range("overdensity not set"); 
         BOOST_LOG_TRIVIAL(debug) << "Setting initial guess for chameleon field...";
 
         T* const f = this->get_y(); // initial guess
         T const* const rho = this->get_external_field(0, 0); // overdensity
         const size_t N_tot = this->get_Ntot();
 
         #pragma omp parallel for
         for (size_t i = 0; i < N_tot; ++i)
         {
             f[i] = chi_min(rho[i]);
         }
     }

template<typename T >

void anonymous_namespace{chameleon.cpp}::ChiSolver< T >::set_convergence	(	double	eps,
		double	err_stop,
		double	err_stop_min,
		double	rms_stop_min,
		size_t	num_fail
	)

inline

Definition at line 529 of file chameleon.cpp.

     {
         this->set_epsilon(eps);
         m_err_stop = err_stop;
         m_err_stop_min = err_stop_min;
         m_rms_stop_min = rms_stop_min;
         m_num_fail = num_fail;
     }

template<typename T >

void anonymous_namespace{chameleon.cpp}::ChiSolver< T >::set_def_convergence ( )

inline

Definition at line 545 of file chameleon.cpp.

     {// set convergence -- compensate for units in which we compute laplace operator
         const double err_mod = get_chi_prefactor();
         set_convergence(err_mod*CONVERGENCE_RES, CONVERGENCE_ERR, CONVERGENCE_ERR_MIN, err_mod*CONVERGENCE_RES_MIN, CONVERGENCE_NUM_FAIL);
         set_bisection_convergence(CONVERGENCE_BI_STEPS_INIT, CONVERGENCE_BI_DCHI, CONVERGENCE_BI_L);
     }

template<typename T >

void anonymous_namespace{chameleon.cpp}::ChiSolver< T >::set_linear	(	Mesh &	rho,
		const FFTW_PLAN_TYPE &	p_F,
		const FFTW_PLAN_TYPE &	p_B
	)

inline

Definition at line 614 of file chameleon.cpp.

     {/* set chameleon guess to linear prediction */
 
         // solve level = 0, use already allocated space and created plans
         set_linear_sol_at_level(rho, p_F, p_B, 0);
 
         // recursively solve at level > 0, create new grids and plans (small)
         set_linear_recursively(1);
     }

template<typename T >

void anonymous_namespace{chameleon.cpp}::ChiSolver< T >::set_linear_recursively ( size_t level )

inline

Definition at line 591 of file chameleon.cpp.

References Mesh::complex(), FFTW_PLAN_C2R, FFTW_PLAN_R2C, FFTW_PLAN_TYPE, Mesh_base< T >::real(), and anonymous_namespace{chameleon.cpp}::transform_Grid_to_Mesh().

     {
         // we are at the bottom level
         if (level >= this->get_Nlevel()) return;
 
         // extract parameters
         const size_t N = this->get_N(level);
 
         // create temporary Mesh to compute linear potential, copy density
         Mesh rho(N);
         transform_Grid_to_Mesh(rho, this->get_external_grid(level, 0));
 
         // initialize FFTW
         const FFTW_PLAN_TYPE p_F = FFTW_PLAN_R2C(N, N, N, rho.real(), rho.complex(), FFTW_ESTIMATE);
         const FFTW_PLAN_TYPE p_B = FFTW_PLAN_C2R(N, N, N, rho.complex(), rho.real(), FFTW_ESTIMATE);
 
         // set linear prediction
         set_linear_sol_at_level(rho, p_F, p_B, level);
 
         // solve linear prediction for the next level
         set_linear_recursively(level + 1);
     }

template<typename T >

void anonymous_namespace{chameleon.cpp}::ChiSolver< T >::set_linear_sol_at_level	(	Mesh &	rho,
		const FFTW_PLAN_TYPE &	p_F,
		const FFTW_PLAN_TYPE &	p_B,
		size_t	level
	)

inline

Definition at line 570 of file chameleon.cpp.

References anonymous_namespace{chameleon.cpp}::CHI_MIN, fftw_execute_dft_c2r(), fftw_execute_dft_r2c(), Mesh::N, and anonymous_namespace{chameleon.cpp}::transform_Mesh_to_Grid().

     {/* set chameleon guess to linear prediction at specific level */
         // check we have the right grids
         if (this->get_N(level) != rho.N) throw std::range_error("Mesh of a different size than Grid at current level!");
 
         // get delta(k)
         fftw_execute_dft_r2c(p_F, rho);
 
         // get dchi(k)
         get_chi_k(rho);
 
         // get dchi(x)
         fftw_execute_dft_c2r(p_B, rho);
 
         // transform dchi into chi, in chi_a units this means only 'chi = 1 + dchi' */
         rho += 1 + CHI_MIN;
 
         // copy dchi(x) onto Grid at level
         transform_Mesh_to_Grid(rho, this->get_grid(level));
     }

template<typename T >

void anonymous_namespace{chameleon.cpp}::ChiSolver< T >::set_screened ( size_t level = 0 )

inline

< we are at the bottom level

< CHI_MIN to indicate high-density region

< recursive call to fix all levels

Definition at line 624 of file chameleon.cpp.

References anonymous_namespace{chameleon.cpp}::CHI_MIN, anonymous_namespace{chameleon.cpp}::MARK_CHI_BOUND_COND, and growth_allz::T.

     {/* check solution for invalid values (non-linear regime), fix values in high density regions, try to improve guess in others */
         if (level >= this->get_Nlevel()) return; 
 
         Grid<3, T>& chi = this->get_grid(level); // guess
         const size_t N_tot = this->get_Ntot(level);
         const size_t N = this->get_N(level);
         T* const rho_grid = this->get_external_field(level, 0); // overdensity
         std::vector<size_t> index_list;
 
         #pragma omp parallel for private(index_list)
         for (size_t i = 0; i < N_tot; ++i)
         {
             if (chi[i] <= CHI_MIN) // non-linear regime
             {
                 if (check_surr_dens(rho_grid, index_list, i, N)) chi[i] = CHI_MIN; 
                 else chi[i] = rho_grid[i] > 0 ? chi_min(rho_grid[i]) : 1/2. + CHI_MIN;//< phi_s / 2 in underdense region as starting point
             }
         }
 
         fix_vals[level].clear();
 
         // writing into map, do not use omp
         for (size_t i = 0; i < N_tot; ++i)
         {
             if (chi[i] == CHI_MIN) // fix chameleon to bulk value, set unphysical density to indicate screened regime
             {
 
                 chi[i] = chi_min(rho_grid[i]);
                 rho_grid[i] = MARK_CHI_BOUND_COND;
                 fix_vals[level].emplace(i, chi[i]);
             }
         }
 
         size_t num_high_density = fix_vals[level].size();
         BOOST_LOG_TRIVIAL(debug) << "Identified and fixed " << num_high_density << "(" << std::setprecision(2) << num_high_density*100.0/N_tot <<  "%) points at level " << level;
 
         set_screened(level + 1); 
     }

template<typename T >

void anonymous_namespace{chameleon.cpp}::ChiSolver< T >::set_time	(	T	a,
		const Cosmo_Param &	cosmo
	)

inline

Definition at line 250 of file chameleon.cpp.

References pow().

     {
         chi_prefactor = chi_prefactor_0*pow(a, -3.*(2.-n)/(1.-n));
     }

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::upd_operator	(	const T	f,
		const size_t	level,
		const std::vector< size_t > &	index_list,
		const T	h
	)		const

inlineoverridevirtual

Reimplemented from MultiGridSolver< 3, T >.

Definition at line 377 of file chameleon.cpp.

References cl_cmbl_bm::l, and growth_allz::T.

     {/* Method for updating solution:
         if df is large, try bisection, otherwise Newton`s method
         try Newton`s method and check for unphysical values */
         T l  =  l_operator(level, index_list, true, h);
         T dl = dl_operator(level, index_list, h);
         T df = -l/dl;
         T df_rel = std::abs(df/(f - CHI_MIN));
 
         static_assert((SWITCH_BIS_NEW < 1), "Newton`s method cannot be allowed with negative values. Adjust 'SWITCH_BIS_NEW < 1'.");
 
         return df_rel < SWITCH_BIS_NEW ? f + df : bisection(f, l, df / 2, level, index_list, h);
     }

Member Data Documentation

template<typename T >

const T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::beta

Chameleon coupling constant.

Definition at line 161 of file chameleon.cpp.

template<typename T >

const T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::chi_0

2*beta*Mpl*phi_scr

Definition at line 162 of file chameleon.cpp.

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::chi_prefactor

time-dependent prefactor

Definition at line 165 of file chameleon.cpp.

template<typename T >

const T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::chi_prefactor_0

dimensionless prefactor to poisson equation at a = 1

Definition at line 164 of file chameleon.cpp.

template<typename T >

std::vector<std::map<size_t, T> > anonymous_namespace{chameleon.cpp}::ChiSolver< T >::fix_vals

<index, value>

Definition at line 180 of file chameleon.cpp.

template<typename T >

size_t anonymous_namespace{chameleon.cpp}::ChiSolver< T >::m_conv_stop = 0

Definition at line 168 of file chameleon.cpp.

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::m_dchi_stop

Definition at line 176 of file chameleon.cpp.

template<typename T >

double anonymous_namespace{chameleon.cpp}::ChiSolver< T >::m_err_stop

Definition at line 170 of file chameleon.cpp.

template<typename T >

double anonymous_namespace{chameleon.cpp}::ChiSolver< T >::m_err_stop_min

Definition at line 171 of file chameleon.cpp.

template<typename T >

T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::m_l_stop

Definition at line 177 of file chameleon.cpp.

template<typename T >

size_t anonymous_namespace{chameleon.cpp}::ChiSolver< T >::m_max_bisection_steps

Definition at line 175 of file chameleon.cpp.

template<typename T >

size_t anonymous_namespace{chameleon.cpp}::ChiSolver< T >::m_num_fail

Definition at line 172 of file chameleon.cpp.

template<typename T >

double anonymous_namespace{chameleon.cpp}::ChiSolver< T >::m_rms_stop_min

Definition at line 169 of file chameleon.cpp.

template<typename T >

const T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::n

Hu-Sawicki paramater.

Definition at line 160 of file chameleon.cpp.

template<typename T >

const T anonymous_namespace{chameleon.cpp}::ChiSolver< T >::phi_prefactor

prefactor for Poisson equation for gravitational potential

Definition at line 163 of file chameleon.cpp.

The documentation for this class was generated from the following file:

src/ApproximationsSchemes/chameleon.cpp

Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

template<typename T> class anonymous_namespace{chameleon.cpp}::ChiSolver< T >

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation

template<typename T>
class anonymous_namespace{chameleon.cpp}::ChiSolver< T >