Gizmod

Introduction

GIZMOD is a fork of the popular physics simulation engine GIZMO, written by Phil Hopkins. The modifications it makes to the core of GIZMO are primarily focused on a new method for analyzing dark matter self-interactions between a host halo and its subhaloes.

The GIZMO codebase is written in C. It also has good documentation, which still holds mostly true for the changes made in GIZMOD and is a good place to start when troubleshooting.

Modifications

Many portions of the GIZMO source code has been modified to implement these changes. For those unfamiliar, I have added the changes in the form of diff code blocks. These represents regions in the codebase that have been changed. Each subheadings under the Modifications section represents a file which has been modified from the original source code. The line numbers on the left side represent where in the file the modifications can be found. Lines which have been added to the code are marked with a +, while any deleted lines are marked with a -.

All these modifications can also be found in the GIZMOD source code in various commits. However, since the commits took place over many months, and not all of them are strictly useful for a new reader (some are simple fixes which don't need to be documented) I have curated the relevant core changes to the code below.

sidm/sidm_core.c

This file contains all of the required functions for the simulating dark matter self-interactions. All of the below modifications to the file are only enabled when the DM_SIDM_ANALYTIC_HOST_TOGGLE flag is set in Config.sh.

Additional Header Files

#include "dispersion_profile.h"
#include <gsl/gsl_interp2d.h>

Ensure we include both our new sidm/dispersion_profile.h header file, as well as the GSL 2D interpolation library. This allows for the use of a 2D interpolation grid and gives access to the variables defined in the new header file.

Integrand Function

double integrand(double x, void *p) {
    struct int_params *params = (struct int_params *)p;
    return g_func(params->A*sqrt(x))/g_func(params->A) * 1.0/sqrt((1 + params->B - x)*(x + params->B - 1));
}

This function will later be set as a gsl function. It takes two inputs: x and a void pointer p. Here x is the input variable of the function and p is a pointer which we cast as a struct int_params pointer, containing the A and B parameters needed to evaluate the function at the point x.

G Function

double g_func(double x) {
    if(All.DM_InteractionVelocityScale>0) {double V = x/All.DM_InteractionVelocityScale;}
    return x*All.DM_InteractionCrossSection/(1 + V*V*V*V);
}

The g_func takes the input: x and returns the result of the following mathematical operation:

\[ g(x):=\frac{x\sigma_0}{1+(\frac{x}{2})^4} \]

Note that we only account for velocity dependent cross sections in the event that All.DM_InteractionVelocityScale is larger than zero.

2D Interpolation Grid Initialization

void init_interp_grid()
{
    A_pts = malloc(A_res*sizeof(double));
    B_pts = malloc(B_res*sizeof(double));
    I_pts = malloc(A_res*B_res*sizeof(double));
    for(int i = 0; i < A_res; i++) {*(A_pts + i) = i*((double)A_max/A_res);}
    for(int j = 0; j < B_res; j++) {*(B_pts + j) = j*((double)B_max/B_res);}
    sh_interp2d_integral = gsl_interp2d_alloc(gsl_interp2d_bilinear, A_res, B_res);

    A_interp_acc = gsl_interp_accel_alloc();
    B_interp_acc = gsl_interp_accel_alloc();

    gsl_function F;
    F.function = &integrand;
    struct int_params params;

    for(size_t i = 1; i < A_res; i++) {
        for(size_t j = 1; j < B_res; j++) {
        /*Calculate Integral*/

        params.A = *(A_pts + i);
        params.B = *(B_pts + j);

        F.params = (void *)&params;
        double result, error;
        size_t neval;
        double lbound = 1 - *(B_pts + j);
        double ubound = 1 + *(B_pts + j);

        int code = gsl_integration_qng(&F, lbound, ubound, epsabs, epsrel, &result, &error, &neval);
        /*Add the integral times g(A)/Pi to the interpolation table*/
        gsl_interp2d_set(sh_interp2d_integral, I_pts, i, j, g_func(params.A)*result/M_PI);
        }
    }
    gsl_interp2d_init(sh_interp2d_integral, A_pts, B_pts, I_pts, A_res, B_res);
}