MP-Gadget  5.0.1.dev1-76bc7d4726-dirty
Functions
test_metal_return.c File Reference
#include <stdarg.h>
#include <stddef.h>
#include <setjmp.h>
#include <cmocka.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_interp2d.h>
#include <stdint.h>
#include "stub.h"
#include "libgadget/utils/endrun.h"
#include "libgadget/metal_return.h"
#include "libgadget/slotsmanager.h"
#include "libgadget/metal_tables.h"
Include dependency graph for test_metal_return.c:

Go to the source code of this file.

Functions

void test_yields (void **state)
 
int main (void)
 

Function Documentation

◆ main()

int main ( void  )

Definition at line 66 of file test_metal_return.c.

66  {
67  const struct CMUnitTest tests[] = {
68  cmocka_unit_test(test_yields),
69  };
70  return cmocka_run_group_tests_mpi(tests, NULL, NULL);
71 }
void test_yields(void **state)

References test_yields().

Here is the call graph for this function:

◆ test_yields()

void test_yields ( void **  state)

Definition at line 20 of file test_metal_return.c.

21 {
22  gsl_integration_workspace * gsl_work = gsl_integration_workspace_alloc(GSL_WORKSPACE);
23  set_metal_params(1.3e-3);
24 
25  struct interps interp;
27  /* Compute factor to normalise the total mass in the IMF to unity.*/
28  double imf_norm = compute_imf_norm(gsl_work);
29  assert_true(fabs(imf_norm - 0.624632) < 0.01);
30 
31  double agbyield = compute_agb_yield(interp.agb_mass_interp, agb_total_mass, 0.01, 1, 40, gsl_work);
32  double agbyield2 = compute_agb_yield(interp.agb_mass_interp, agb_total_mass, 0.01, 1, SNAGBSWITCH, gsl_work);
33  assert_true(fabs(agbyield / agbyield2 - 1) < 1e-3);
34  /* Lifetime is about 200 Myr*/
35  double agbyield3 = compute_agb_yield(interp.agb_mass_interp, agb_total_mass, 0.01, 5, 40, gsl_work);
36 
37  /* Integrate the region of the IMF which contains SNII and AGB stars. The yields should never be larger than this*/
38  gsl_function ff = {chabrier_mass, NULL};
39  double agbmax, sniimax, abserr;
40  gsl_integration_qag(&ff, agb_total_mass[0], SNAGBSWITCH, 1e-4, 1e-3, GSL_WORKSPACE, GSL_INTEG_GAUSS61, gsl_work, &agbmax, &abserr);
41  gsl_integration_qag(&ff, SNAGBSWITCH, snii_masses[SNII_NMASS-1], 1e-4, 1e-3, GSL_WORKSPACE, GSL_INTEG_GAUSS61, gsl_work, &sniimax, &abserr);
42 
43  double sniiyield = compute_snii_yield(interp.snii_mass_interp, snii_total_mass, 0.01, 1, 40, gsl_work);
44 
45  double sn1a = sn1a_number(0, 1500, 0.679)*sn1a_total_metals;
46  assert_true(sn1a < 1.3e-3);
47 
48  message(0, "agbyield %g max %g (in 200 Myr: %g)\n", agbyield, agbmax, agbyield3);
49  message(0, "sniiyield %g max %g sn1a %g\n", sniiyield, sniimax, sn1a);
50  message(0, "Total fraction of mass returned %g\n", (sniiyield + sn1a + agbyield)/imf_norm);
51  assert_true(agbyield < agbmax);
52  assert_true(sniiyield < sniimax);
53  assert_true((sniiyield + sn1a + agbyield)/imf_norm < 1.);
54 
55  double masslow1, masshigh1;
56  double masslow2, masshigh2;
57  double masslowsum, masshighsum;
58  find_mass_bin_limits(&masslow1, &masshigh1, 0, 30, 0.02, interp.lifetime_interp);
59  find_mass_bin_limits(&masslow2, &masshigh2, 30, 60, 0.02, interp.lifetime_interp);
60  find_mass_bin_limits(&masslowsum, &masshighsum, 0, 60, 0.02, interp.lifetime_interp);
61  message(0, "0 - 30: %g %g 30 - 60 %g %g 0 - 60 %g %g\n", masslow1, masshigh1, masslow2, masshigh2, masslowsum, masshighsum);
62  assert_true(fabs(masslow1 - masshigh2) < 0.01);
63  assert_true(fabs(masslowsum - masslow2) < 0.01);
64 }
Interp interp
void message(int where, const char *fmt,...)
Definition: endrun.c:175
void set_metal_params(double Sn1aN0)
Definition: metal_return.c:55
void find_mass_bin_limits(double *masslow, double *masshigh, const double dtstart, const double dtend, double stellarmetal, gsl_interp2d *lifetime_tables)
Definition: metal_return.c:230
void setup_metal_table_interp(struct interps *interp)
Definition: metal_return.c:76
double sn1a_number(double dtmyrstart, double dtmyrend, double hub)
Definition: metal_return.c:324
double compute_snii_yield(gsl_interp2d *snii_interp, const double *snii_weights, double stellarmetal, double masslow, double masshigh, gsl_integration_workspace *gsl_work)
Definition: metal_return.c:368
double compute_agb_yield(gsl_interp2d *agb_interp, const double *agb_weights, double stellarmetal, double masslow, double masshigh, gsl_integration_workspace *gsl_work)
Definition: metal_return.c:342
double compute_imf_norm(gsl_integration_workspace *gsl_work)
Definition: metal_return.c:314
double chabrier_mass(double mass, void *params)
Definition: metal_return.c:308
#define SNII_NMASS
Definition: metal_tables.h:309
#define SNAGBSWITCH
Definition: metal_tables.h:13
#define GSL_WORKSPACE
Definition: metal_tables.h:426
static const double snii_total_mass[SNII_NMET *SNII_NMASS]
Definition: metal_tables.h:312
static const double agb_total_mass[AGB_NMET *AGB_NMASS]
Definition: metal_tables.h:72
static const double sn1a_total_metals
Definition: metal_tables.h:60
static const double snii_masses[SNII_NMASS]
Definition: metal_tables.h:310

References agb_total_mass, chabrier_mass(), compute_agb_yield(), compute_imf_norm(), compute_snii_yield(), find_mass_bin_limits(), GSL_WORKSPACE, interp, message(), set_metal_params(), setup_metal_table_interp(), sn1a_number(), sn1a_total_metals, SNAGBSWITCH, snii_masses, SNII_NMASS, and snii_total_mass.

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function: