From a6b3bac2f305c13e9b00424f00e9ad3dbbeeaaab Mon Sep 17 00:00:00 2001 From: Simeon Bird Date: Thu, 26 Dec 2024 11:55:14 -0800 Subject: [PATCH] Port metal_return test, last of the libgadget tests --- libgadget/Makefile | 4 +- libgadget/tests/test_metal_return.c | 59 +++++++++++++---------------- 2 files changed, 29 insertions(+), 34 deletions(-) diff --git a/libgadget/Makefile b/libgadget/Makefile index cdac24f8..4faf55fd 100644 --- a/libgadget/Makefile +++ b/libgadget/Makefile @@ -21,8 +21,8 @@ TESTED = cooling \ neutrinos_lra \ omega_nu_single \ density \ - gravity -# metal_return \ + gravity \ + metal_return MPI_TESTED = fof exchange diff --git a/libgadget/tests/test_metal_return.c b/libgadget/tests/test_metal_return.c index 055dfc2e..e2f2bd48 100644 --- a/libgadget/tests/test_metal_return.c +++ b/libgadget/tests/test_metal_return.c @@ -1,23 +1,16 @@ /*Tests for the drift factor module.*/ -#include -#include -#include -#include -#include -#include -#include -#include +#define BOOST_TEST_MODULE metal_return +#include "booststub.h" + #include -#include -#include +#include -#include "stub.h" #include "libgadget/utils/endrun.h" #include "libgadget/metal_return.h" #include "libgadget/slotsmanager.h" #include "libgadget/metal_tables.h" -void test_yields(void ** state) +BOOST_AUTO_TEST_CASE(test_yields) { gsl_integration_workspace * gsl_work = gsl_integration_workspace_alloc(GSL_WORKSPACE); set_metal_params(1.3e-3); @@ -26,31 +19,40 @@ void test_yields(void ** state) setup_metal_table_interp(&interp); /* Compute factor to normalise the total mass in the IMF to unity.*/ double imf_norm = compute_imf_norm(gsl_work); - assert_true(fabs(imf_norm - 0.624632) < 0.01); + BOOST_TEST(imf_norm == 0.624632, tt::tolerance(0.01)); double agbyield = compute_agb_yield(interp.agb_mass_interp, agb_total_mass, 0.01, 1, 40, gsl_work); double agbyield2 = compute_agb_yield(interp.agb_mass_interp, agb_total_mass, 0.01, 1, SNAGBSWITCH, gsl_work); - assert_true(fabs(agbyield / agbyield2 - 1) < 1e-3); + BOOST_TEST(agbyield == agbyield2, tt::tolerance(1e-3)); /* Lifetime is about 200 Myr*/ double agbyield3 = compute_agb_yield(interp.agb_mass_interp, agb_total_mass, 0.01, 5, 40, gsl_work); - /* Integrate the region of the IMF which contains SNII and AGB stars. The yields should never be larger than this*/ - gsl_function ff = {chabrier_mass, NULL}; - double agbmax, sniimax, abserr; - gsl_integration_qag(&ff, agb_total_mass[0], SNAGBSWITCH, 1e-4, 1e-3, GSL_WORKSPACE, GSL_INTEG_GAUSS61, gsl_work, &agbmax, &abserr); - gsl_integration_qag(&ff, SNAGBSWITCH, snii_masses[SNII_NMASS-1], 1e-4, 1e-3, GSL_WORKSPACE, GSL_INTEG_GAUSS61, gsl_work, &sniimax, &abserr); + /* Integrate the region of the IMF which contains SNII and AGB stars. The yields should never be larger than this + * The Chabrier IMF used for computing SnII and AGB yields. + * See 1305.2913 eq 3*/ + auto chabrier_mass = [](const double mass) { + double imf; + if(mass <= 1) + imf = 0.852464 / mass * exp(- pow(log(mass / 0.079)/ 0.69, 2)/2); + else + imf = 0.237912 * pow(mass, -2.3); + return mass * imf; + }; + // Gauss-Kronrod integration for smooth functions. Boost uses by default the machine precision for accuracy and a max depth of 15. + const double agbmax = boost::math::quadrature::gauss_kronrod::integrate(chabrier_mass, agb_total_mass[0], SNAGBSWITCH); + const double sniimax = boost::math::quadrature::gauss_kronrod::integrate(chabrier_mass, SNAGBSWITCH, snii_masses[SNII_NMASS-1]); double sniiyield = compute_snii_yield(interp.snii_mass_interp, snii_total_mass, 0.01, 1, 40, gsl_work); double sn1a = sn1a_number(0, 1500, 0.679)*sn1a_total_metals; - assert_true(sn1a < 1.3e-3); + BOOST_TEST(sn1a < 1.3e-3); message(0, "agbyield %g max %g (in 200 Myr: %g)\n", agbyield, agbmax, agbyield3); message(0, "sniiyield %g max %g sn1a %g\n", sniiyield, sniimax, sn1a); message(0, "Total fraction of mass returned %g\n", (sniiyield + sn1a + agbyield)/imf_norm); - assert_true(agbyield < agbmax); - assert_true(sniiyield < sniimax); - assert_true((sniiyield + sn1a + agbyield)/imf_norm < 1.); + BOOST_TEST(agbyield < agbmax); + BOOST_TEST(sniiyield < sniimax); + BOOST_TEST((sniiyield + sn1a + agbyield)/imf_norm < 1.); double masslow1, masshigh1; double masslow2, masshigh2; @@ -59,13 +61,6 @@ void test_yields(void ** state) find_mass_bin_limits(&masslow2, &masshigh2, 30, 60, 0.02, interp.lifetime_interp); find_mass_bin_limits(&masslowsum, &masshighsum, 0, 60, 0.02, interp.lifetime_interp); message(0, "0 - 30: %g %g 30 - 60 %g %g 0 - 60 %g %g\n", masslow1, masshigh1, masslow2, masshigh2, masslowsum, masshighsum); - assert_true(fabs(masslow1 - masshigh2) < 0.01); - assert_true(fabs(masslowsum - masslow2) < 0.01); -} - -int main(void) { - const struct CMUnitTest tests[] = { - cmocka_unit_test(test_yields), - }; - return cmocka_run_group_tests_mpi(tests, NULL, NULL); + BOOST_TEST(masslow1 == masshigh2, tt::tolerance(0.01)); + BOOST_TEST(masslowsum == masslow2, tt::tolerance(0.01)); }