model.cu

This section contains model dependent functions defined in the model.cu file related to the iharm model.

Macros

Defines

NUMIN 1.e8

Lowest sampled frequency for superphotons in the simulation in Hz.

NUMAX 1.e24

Highest sampled frequency for superphotons in the simulation in Hz.

Note

Superphotons with frequencies above NUMAX can be generated through scattering events.

LNUMIN log(NUMIN)

Natural logarithm of the lower sampling bound.

Derived from NUMIN as log(NUMIN), this value represents the log-space minimum used for energy/frequency grid construction.

LNUMAX log(NUMAX)

Natural logarithm of the upper sampling bound.

Derived from NUMAX as log(NUMAX), this value represents the log-space maximum used for energy/frequency grid construction.

DLNU ((LNUMAX-LNUMIN)/N_ESAMP)

Log-space step size between samples.

Computed from the log-range [LNUMIN, LNUMAX] divided by N_ESAMP, giving the increment used to traverse the sampling domain uniformly in natural-log space.

IHARM (1)

Macro to identify the iharm GRMHD model.

DO_NOT_USE_TEXTURE_MEMORY 1

Macro not to use texture memory for fluid variables table in device code.

THETAE_MIN 0.3

Minimum dimensionless electron temperature ( \(\Theta_{e, min}\)) for physical validity.

This is used to characterize the lower bound for superphoton generation.

THETAE_MAX 1000.

Maximum dimensionless electron temperature ( \(\Theta_{e, max}\)) for physical validity.

dlE (0.25)

Size of the energy bin in logarithmic scale for the spectral output binning.

lE0 (log(1.e-12))

Minimum energy of the energy bin in logarithmic scale.

N_ESAMP 800

Number of energy samples for the emissivity and weight tables.

N_EBINS 800

Number of energy bins for the spectral output.

WEIGHT_MIN (1.e28)

Minimum superphoton weight to be considered to the spectrum. Superphotons that fall below this weight will be terminated.

RMAX 1000.

Maximum radius to track photons in code units. If photons exceed this radius, their tracking will stop and they will be recorded.

ROULETTE 1.e4

Roulette factor for photon termination based on weight.

NPRIM 10

Number of primitive variables in the iharm model. This should match the value read from the HDF5 file.

MODEL_FUNCTIONS
USE_FIXED_TPTE (0)

Activates the constant temperature ratio model. When enabled (1), the code assumes a globally uniform ratio between proton and electron temperatures ( \( T_p/T_e \)), defined by the params.tp_over_te variable.

  • Sets with_electrons = 0.

  • Calculates Thetae_unit based on a fixed partitioning of the internal energy, often assuming a non-relativistic monoatomic gas ( \( \gamma = 5/3 \)).

USE_MIXED_TPTE (1)

Activates the \( R-\beta \) (Mixed) temperature ratio model. When enabled (1), the code uses the plasma \( \beta \) to interpolate between two temperature ratios: trat_small (for highly magnetized regions like the jet) and trat_large (for the accretion disk).

  • Sets with_electrons = 2.

WITH_ELECTRONS (2)

Selects the sub-grid physics model for electron thermodynamics.

  • Available Models:

  • 0: Fixed Ratio * Assumes a constant ratio between proton and electron temperatures throughout the entire simulation domain. \( \Theta_e \propto \frac{u}{\rho} \).

  • 1: Howes/Kawazura Model * A kinetic-based prescription that tracks electron heating based on local plasma turbulence and entropy ( \( \kappa_{el} \)).

  • 2: R-high / R-low Model (Default) * The standard \(\beta\)-dependent model for Black Hole imaging. It varies the temperature ratio \( R = T_p/T_e \) based on the local plasma \( \beta \) (the ratio of gas to magnetic pressure):

  • In the Disk (high \(\beta\)): Protons are much hotter than electrons ( \( R \to R_{high} \)).

  • In the Jet (low \(\beta\)): Electrons and protons are closer in temperature ( \( R \to R_{low} \)).

ZLOOP for (int i = 0; i < N1; i++) for (int j = 0; j < N2; j++) for (int k = 0; k < N3; k++)

Macro for looping over all grid zones in 3D.

METRIC_eKS 0
METRIC_MKS 1
METRIC_MKS3 2
METRIC_FMKS 3

Functions

Functions

__host__ void init_storage(void)

Allocates memory on the host for fluid primitives and grid geometry.

This function reserves space in the system RAM for the 4D primitive variables array

and the 2D spacetime metric (geometry) array. This must be called before reading any simulation snapshots or calculating the metric.

Returns:

void

__host__ void init_data()

Main initialization routine for physical data and simulation grid.

This function performs the following critical tasks:

  1. Opens and parses the HDF5 dump file.

  2. Sets up the computational grid dimensions (N1, N2, N3).

  3. Configures the electron temperature model (Fixed, R-beta, or Custom).

  4. Loads the spacetime geometry (MKS/KS) and calculates the event horizon.

  5. Allocates memory and populates the primitive variable arrays.

  6. Calculates initial diagnostics (accretion rate, bias normalization).

Returns:

void

__device__ int record_criterion(double X1)

Determines if a photon has exited the simulation domain and should be recorded.

This function checks if the photon has crossed a specific outer boundary.

Parameters:

X1 – The current internal radial coordinate ( \(X^1\)) of the photon.

Returns:

Returns 1 if the photon is outside the recording radius, 0 otherwise.

__device__ int stop_criterion(double X1, double *w, curandState *localState)

Evaluates whether to terminate the integration of a photon’s path. This function checks for two types of termination:

  1. Physical: Falling into the Black Hole or leaving the domain.

  2. Statistical: Using a “Roulette” to terminate photons with negligible weights.

Parameters:

  • X1 – Current radial coordinate in code units.

  • w – Pointer to the photon’s weight (energy/intensity).

  • localState – Pointer to the GPU’s random number generator state.

Returns:

Returns 1 if the photon should be stopped/deleted, 0 if tracking should continue.

__device__ void Xtoijk(const double X[NDIM], int *i, int *j, int *k, double del[NDIM])

Maps continuous coordinates to discrete grid indices and interpolation weights.

Bridge between the geodesic integrator and the fluid grid. It identifies the nearest cell and calculates the fractional offset (del) required for trilinear interpolation of fluid properties.

Parameters:

  • X – Input continuous internal coordinates \( X^\mu \).

  • i, j, k – Output integer indices for the grid cell.

  • del – Output array [1..3] containing the fractional distance within the cell [0, 1].

Returns:

void

__host__ __device__ void coord(const int i, const int j, const int k, double *X)

Maps discrete grid indices to continuous internal coordinates.

Calculates the coordinate of a cell center.

Parameters:

  • i, j, k – Integer grid indices.

  • X – Output array to store the resulting coordinates \( X^\mu \).

Returns:

void

__host__ __device__ void gcov_func(const double *X, double gcov[][NDIM])

Computes the covariant metric tensor \( g_μν \) in Modified Kerr-Schild (MKS) coordinates.

It transforms the standard Kerr-Schild metric into the simulation’s specific Modified Kerr-Schild (MKS) coordinate system using Jacobian factors (tfac, rfac, etc.).

Parameters:

  • X – Input internal coordinates \( X^\mu \).

  • gcov – Output 4x4 array where the metric components will be stored.

__host__ double dOmega_func(double x2i, double x2f)

Calculates the integrated solid angle between two polar grid boundaries. This function computes the solid angle by evaluating the integral:

\[ \Delta\Omega = \int_{0}^{2\pi} d\phi \int_{\theta_i}^{\theta_f} \sin\theta \, d\theta = 2\pi (\cos\theta_i - \cos\theta_f) \]
where \( \theta \) is the physical polar angle derived from the internal coordinate \( X^2 \).

Parameters:

  • x2i – Starting internal polar coordinate (X2).

  • x2f – Ending internal polar coordinate (X2).

Returns:

The solid angle in steradians.

__host__ __device__ void bl_coord(const double *X, double *r, double *th)

Transforms internal simulation coordinates to physical Boyer-Lindquist coordinates. This function maps the logarithmic radial coordinate and the stretched polar coordinate to their physical counterparts:

  • \[ r = e^{X^1} \]
    \[ \theta = \pi X^2 + \frac{1 - h}{2} \sin(2\pi X^2) \]
    where \( h \) is the hslope parameter controlling equatorial refinement.

Note

\( X^3 \) (azimuthal angle) remains unchanged in this transformation.

Parameters:

  • X – Input array of internal coordinates \( X^\mu \).

  • r – Output pointer for the Boyer-Lindquist radius.

  • th – Output pointer for the Boyer-Lindquist polar angle (theta).

__host__ __device__ void get_fluid_zone(const int i, const int j, const int k, double *Ne, double *Thetae, double *B, double Ucon[NDIM], double Bcon[NDIM], const struct of_geom *d_geom, const double *d_p)

Recovers physical fluid and magnetic field properties from primitive variables.

This function performs the transformation from the GRMHD code’s primitive variables (relative to a normal observer) to the physical quantities in the fluid’s local frame.

It calculates the 4-velocity \( u^\mu \) ensuring the normalization condition:

\[ u_\mu u^\mu = -1 \]

And the magnetic 4-vector \( b^\mu \) satisfying the ideal MHD condition:

\[ b^\mu u_\mu = 0 \]

Parameters:

  • i, j, k – Grid indices.

  • Ne – [out] Physical electron number density \( n_e \).

  • Thetae – [out] Dimensionless electron temperature \( \Theta_e \).

  • B – [out] Fluid-frame magnetic field strength.

  • Ucon – [out] Contravariant 4-velocity \( u^\mu \).

  • Bcon – [out] Contravariant magnetic 4-vector \( b^\mu \).

  • d_geom – Pointer to the geometry metric structure.

  • d_p – Pointer to the primitive variable array.

Returns:

void

__device__ void get_fluid_params(double X[NDIM], double *Ne, double *Thetae, double *B, double Ucon[NDIM], double Ucov[NDIM], double Bcon[NDIM], double Bcov[NDIM], double *d_p)

Samples fluid and magnetic properties at an arbitrary spatial point using trilinear interpolation.

  • This function is called during geodesic integration to determine local plasma conditions. It performs the following sequence:

  1. Validates the photon’s position against the grid boundaries.

  2. Computes weights for trilinear interpolation between the 8 neighboring cells.

  3. Interpolates primitive variables (density, energy, velocity, magnetic field).

  4. Reconstructs the 4-velocity and magnetic 4-vector using the local metric.

  5. Applies safety cuts.

Parameters:

  • X – Continuous internal coordinates.

  • Ne – [out] Interpolated electron number density.

  • Thetae – [out] Interpolated dimensionless electron temperature.

  • B – [out] Physical magnetic field magnitude.

  • Ucon, Ucov – [out] Reconstructed 4-velocity vectors.

  • Bcon, Bcov – [out] Reconstructed magnetic 4-vectors.

  • d_p – Pointer to the global primitive variable array in GPU memory.

Returns:

void

__device__ void gcov_func_row0(const double *X, double gcov_row0[NDIM])

Computes the first row of the covariant metric tensor \( g_{μν} \).

Parameters:

  • X – Input internal coordinates \( X^\mu \).

  • gcov_row0 – Output array to store the first row of the covariant metric tensor \( g_{0ν} \).

Returns:

void

__device__ double bias_func(double Te, double w, int round_scatt)

Calculates the scattering bias parameter based on local temperature. This function determines the “bias factor” used to adjust scattering probabilities. It prioritizes high-temperature regions (where \( \Theta_e \) is large) to improve the signal-to-noise ratio of the final spectrum.

Parameters:

  • Te – Dimensionless electron temperature ( \( \Theta_e \)).

  • w – Current statistical weight of the photon.

  • round_scatt – The current scattering iteration count.

Returns:

The calculated bias factor used to rescale Monte Carlo probabilities.

__device__ double stepsize(const double X[NDIM], const double K[NDIM])

Performs the calculation of the stepsize \( d\lambda \) for the geodesic integration. The stepsize follows this formula:

\( d\lambda = \left[ \frac{|K^1| + \epsilon_{s}}{\epsilon_{e}} + \frac{|K^2| + \epsilon_{s}}{\epsilon_{e} \cdot \min(X^2, 1-X^2)} + \frac{|K^3| + \epsilon_{s}}{\epsilon_{e}} \right]^{-1} \)

Where:

  • \( \epsilon_{e} \) is the tolerance parameter (EPS).

  • \( \epsilon_{s} \) is the safety floor (SMALL) to prevent division by zero.

  • \( \min(X^2, 1-X^2) \) is the proximity factor that reduces the step size as the photon approaches the polar coordinate singularities at 0 or 1.

__host__ __device__ double thetae_func(double uu, double rho, double B, double kel)

Calculates the dimensionless electron temperature \( \Theta_e \). This function provides several models for the electron-to-proton temperature ratio. It is primarily used to map the single-fluid internal energy from GRMHD onto the electron populations that actually produce the radiation.

Parameters:

  • uu – Local internal energy density.

  • rho – Local mass density.

  • B – Local magnetic field strength.

  • kel – Local entropy/electron constant (used for mode 1).

Returns:

The effective dimensionless electron temperature, soft-clamped by Thetae_max.

__host__ void report_spectrum_h5(unsigned long long N_superph_made, struct of_spectrum ***spect, const char *filename)

Processes energy and angled binned simulation data to generate and save the final spectrum in hdf5 format.

This function converts the raw energy and photon counts accumulated in the spectral grid into physical units. It calculates the SED across different inclination angles, determines the average optical depths, and computes global simulation diagnostics such as total luminosity and accretion efficiency. Contains the fluid_header, params and output groups.

Parameters:

  • N_superph_made – Total number of superphotons generated during the run.

  • spect – 3D array containing the accumulated spectral data (Photon physical origin,Energy,Theta bins).

  • filename – Name of the output file (saved in the ./output/ directory).