Port metal_return test, last of the libgadget tests

libgadget/Makefile

@@ -21,8 +21,8 @@ TESTED = cooling \
 	neutrinos_lra \
 	omega_nu_single \
 	density \
-	gravity
-#	metal_return \
+	gravity \
+	metal_return
 
 MPI_TESTED = fof exchange
 

libgadget/tests/test_metal_return.c

@@ -1,23 +1,16 @@
 /*Tests for the drift factor module.*/
-#include <stdarg.h>
-#include <stddef.h>
-#include <setjmp.h>
-#include <cmocka.h>
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
+#define BOOST_TEST_MODULE metal_return
+#include "booststub.h"
+
 #include <gsl/gsl_integration.h>
-#include <gsl/gsl_interp2d.h>
-#include <stdint.h>
+#include <boost/math/quadrature/gauss_kronrod.hpp>
 
-#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<double, 61>::integrate(chabrier_mass, agb_total_mass[0], SNAGBSWITCH);
+    const double sniimax = boost::math::quadrature::gauss_kronrod<double, 61>::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));
 }