test_chameleon.cpp File Reference

#include <catch.hpp>
#include "../test.hpp"
#include "chameleon.cpp"

Include dependency graph for test_chameleon.cpp:

Go to the source code of this file.

Namespaces
	anonymous_namespace{test_chameleon.cpp}

Functions
template<typename T >
T	anonymous_namespace{test_chameleon.cpp}::mean (const std::vector< T > &data)
double	anonymous_namespace{test_chameleon.cpp}::mean (const Mesh &data)
void	anonymous_namespace{test_chameleon.cpp}::get_neighbor_gridindex (std::vector< size_t > &index_list, size_t i, size_t N)
void	anonymous_namespace{test_chameleon.cpp}::init_overdensity (Sim_Param &sim, Mesh &rho, MultiGrid< 3, CHI_PREC_t > &rho_grid)
void	anonymous_namespace{test_chameleon.cpp}::get_grav_pot (Mesh &rho, const FFTW_PLAN_TYPE &p_F, const FFTW_PLAN_TYPE &p_B, double box_size, double phi_prefactor)
template<typename T >
static int	anonymous_namespace{test_chameleon.cpp}::sgn (T val)
void	anonymous_namespace{test_chameleon.cpp}::print_mesh (const std::string &file_name, const Mesh &pot, const double mod=-1-CHI_MIN)
	TEST_CASE ("UNIT TEST: create Multigrid and copy data to/from Mesh","[multigrid]")
	TEST_CASE ("UNIT TEST: create and initialize ChiSolver, check bulk field","[chameleon]")
	TEST_CASE ("UNIT TEST: create and initialize ChiSolver, solve sphere","[chameleon]")

Function Documentation

TEST_CASE	(	"UNIT TEST: create Multigrid and copy data to/from Mesh"	,
		""	[multigrid]
	)

Definition at line 111 of file test_chameleon.cpp.

References CHECK, anonymous_namespace{chameleon.cpp}::min(), anonymous_namespace{chameleon.cpp}::transform_Mesh_to_MultiGrid(), and anonymous_namespace{chameleon.cpp}::transform_MultiGrid_to_Mesh().

                     : create Multigrid and copy data to/from Mesh", "[multigrid]" )
{
    print_unit_msg("create Multigrid and copy data to/from Mesh");

    constexpr size_t N = 32;
    srand(time(0));
    MultiGrid<3, FTYPE_t> grid(N);
    Mesh mesh_from(N);
    Mesh mesh_to(N);

    const std::vector<size_t> some_indices = // Mesh indices = N*N*(N+2)
                                            {4, // (0, 0, 4)
                                            (N+2)*2+5, // (0, 2, 5)
                                            N*(N+2)*4+8*(N+2)+5}; // (4, 8, 5)
    for(size_t i : some_indices) mesh_from[i] = 2.0*rand()/RAND_MAX - 1;

    transform_Mesh_to_MultiGrid(mesh_from, grid);
    CHECK(min(mesh_from) == min(grid));
    transform_MultiGrid_to_Mesh(mesh_to, grid);

    for(size_t i : some_indices) CHECK( mesh_from[i] == mesh_to[i] );
    CHECK(min(mesh_from) == min(grid));
    CHECK(min(mesh_to) == min(grid));
}

TEST_CASE	(	"UNIT TEST: create and initialize	ChiSolver,
		check bulk field"	,
		""	[chameleon]
	)

Definition at line 136 of file test_chameleon.cpp.

References ccl_test_power::argv, Mesh_base< T >::assign(), Chi_Opt::beta, CHECK, CHECK_THROWS_AS, Sim_Param::chi_opt, Sim_Param::cosmo, anonymous_namespace{test_chameleon.cpp}::get_neighbor_gridindex(), ccl_test_distances::h, Chi_Opt::n, REQUIRE, sqrt(), and anonymous_namespace{chameleon.cpp}::transform_Mesh_to_MultiGrid().

                     : create and initialize ChiSolver, check bulk field", "[chameleon]" )
{
    print_unit_msg("create and initialize ChiSolver, check bulk field");

    constexpr size_t N = 32;
    constexpr FTYPE_t a = 0.5;
    constexpr FTYPE_t rho_0 = 0.3;
    const std::vector<size_t> some_indices = // Mesh indices = N*N*(N+2)
                                            {4, // (0, 0, 4)
                                            (N+2)*2+5, // (0, 2, 5)
                                            N*(N+2)*4+8*(N+2)+5}; // (4, 8, 5)

    // initialize Sim_Param
    const int argc = 1;
    const char* const argv[1] = {"test"};
    Sim_Param sim(argc, argv);

    // check thowing of exception with unphysical values of chameleon parameters
    auto constructor = [&sim](){ ChiSolver<FTYPE_t> _sol(N, sim, false); };
    sim.chi_opt.n = -1;
    CHECK_THROWS_AS( constructor(), std::out_of_range );
    sim.chi_opt.n = 0.5;
    sim.chi_opt.beta = -1;
    CHECK_THROWS_AS( constructor(), std::out_of_range );

    // initialize ChiSolver
    sim.chi_opt.beta = 1/sqrt(6.);
    ChiSolver<FTYPE_t> sol(N, sim, false);

    // check thowing of exceptions in uninitialised state
    CHECK_THROWS_AS( sol.set_bulk_field(), std::out_of_range );

    // initialize overdensity -- constant density
    MultiGrid<3, FTYPE_t> rho_grid(N);
    {
        Mesh rho(N);
        rho.assign(rho_0);
        transform_Mesh_to_MultiGrid(rho, rho_grid);
    }

    // check density
    for(size_t i : some_indices) CHECK( rho_grid[0][i] == rho_0 );

    // set ChiSolver -- bulk field
    sol.set_time(a, sim.cosmo);
    sol.add_external_grid(&rho_grid);
    sol.set_bulk_field();

    // check bulk field
    {
        const FTYPE_t chi_bulk = sol.chi_min(rho_0);
        FTYPE_t const* const chi = sol.get_y();
        for(size_t i : some_indices) REQUIRE( chi[i] == Approx(chi_bulk));
    }

    // check that EOM is satisfied
    std::vector<size_t> index_list;
    size_t level = 0;
    const FTYPE_t h = 1.0/FTYPE_t( sol.get_N(level) );
    for(size_t i : some_indices)
    {
        get_neighbor_gridindex(index_list, i, N);
        CHECK( sol.l_operator(level, index_list, true, h) == Approx(0.));
    }
}

TEST_CASE	(	"UNIT TEST: create and initialize	ChiSolver,
		solve sphere"	,
		""	[chameleon]
	)

Definition at line 202 of file test_chameleon.cpp.

References ccl_test_power::argv, Sim_Param::box_opt, Box_Opt::box_size, Test_Opt::coarse_sweeps, Mesh::complex(), Sim_Param::cosmo, create_dir(), FFTW_DEST_PLAN, FFTW_PLAN_C2R, FFTW_PLAN_OMP, FFTW_PLAN_OMP_CLEAN, FFTW_PLAN_OMP_INIT, FFTW_PLAN_R2C, FFTW_PLAN_TYPE, Test_Opt::fine_sweeps, anonymous_namespace{test_chameleon.cpp}::get_grav_pot(), anonymous_namespace{test_chameleon.cpp}::init_overdensity(), Test_Opt::max_steps, Test_Opt::N_grid, Test_Opt::N_min, Run_Opt::nt, Out_Opt::out_dir, Sim_Param::out_opt, Sim_Param::print_info(), anonymous_namespace{test_chameleon.cpp}::print_mesh(), Mesh_base< T >::real(), remove_all_files(), Sim_Param::run_opt, Test_Opt::step_per_iter, Sim_Param::test_opt, anonymous_namespace{chameleon.cpp}::transform_MultiGridSolver_to_Mesh(), and Test_Opt::verbose.

                     : create and initialize ChiSolver, solve sphere", "[chameleon]" )
{
    print_unit_msg("create and initialize ChiSolver, solve sphere");

    // initialize Sim_Param
    const char* const argv[1] = {"test"};
    Sim_Param sim(1, argv);
    const size_t N = sim.test_opt.N_grid;
    const size_t N_min = sim.test_opt.N_min;
    const bool verbose = sim.test_opt.verbose;

    // initialize ChiSolver
    ChiSolver<CHI_PREC_t> sol(N, N_min, sim, verbose);

    // initialize overdensity -- constant density in sphere of radius R, center at x0, y0, z0
    MultiGrid<3, CHI_PREC_t> rho_grid(N);
    Mesh rho(N);
    init_overdensity(sim, rho, rho_grid);

    // FFTW preparation
    if (!FFTW_PLAN_OMP_INIT()){
        throw std::runtime_error("Errors during multi-thread initialization");
    }
    FFTW_PLAN_OMP(sim.run_opt.nt);
    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 ChiSolver
    sol.set_time(1, sim.cosmo);
    sol.set_def_convergence();
    sol.add_external_grid(&rho_grid);
    sol.set_ngs_sweeps(sim.test_opt.fine_sweeps, sim.test_opt.coarse_sweeps); //< fine, coarse

    // compute gravitational potential
    Mesh phi_pot(rho); //< copy density
    get_grav_pot(phi_pot, p_F, p_B, sim.box_opt.box_size, sol.get_phi_prefactor());

    // get linear prediction
    sol.set_linear(rho, p_F, p_B);
    sol.set_screened();

    // full solution on Mesh
    Mesh chi_full(N);

    // create directory structure
    std::string out_dir = sim.out_opt.out_dir + "test_ChiSolver/";
    remove_all_files(out_dir);
    create_dir(out_dir);
    sim.print_info(out_dir, "test");

    // print gravitational potential
    print_mesh(out_dir + "grav_pot.dat", phi_pot, 0);

    // Solve the equation -- full V-cycles
    BOOST_LOG_TRIVIAL(info) << "Starting V-cycles with intermediate output...\n";
    
    int istep = 0;
    const size_t step_per_iter = sim.test_opt.step_per_iter;
    while(1)
    {
        // istep, max_step (use int for the last iteration -- negative max_step)
        istep += sol.get_istep();
        int max_step = istep + step_per_iter > sim.test_opt.max_steps ? sim.test_opt.max_steps - istep : step_per_iter;

        // print chi_full
        transform_MultiGridSolver_to_Mesh(chi_full, sol);
        print_mesh(out_dir + "chi_istep_" + std::to_string(istep) + ".dat", chi_full);

        // check max_step and convergence
        if ((max_step <= 0) || ((istep != 0) && (sol.get_istep() < step_per_iter))) break;

        // set & solve
        sol.set_maxsteps(max_step);
        sol.solve();
    }

    // FFTW CLEANUP
    FFTW_DEST_PLAN(p_F);
    FFTW_DEST_PLAN(p_B);
    FFTW_PLAN_OMP_CLEAN();
}