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.e16

Highest sampled frequency for superphotons in the simulation in Hz.

Note

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

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.12)

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 8

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

MODEL_FUNCTIONS
GAME (4./3.)

Adiabatic index for electrons.

GAMP (5./3.)

Adiabatic index for the protons.

GAM (1.4444440)

Adiabatic index for the total fluid.

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.

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 gcov[NDIM][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.

  • gcov – [out] Local covariant metric tensor.

  • 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__ 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__ void record_super_photon(struct of_photonSOA ph, struct of_spectrum *d_spect, unsigned long long photon_index)

Records a photon’s final properties into the global spectrum array. This function maps a photon’s energy and polar exit angle to a specific energy and angular bin in the spectrum. Because many threads may try to update the same bin simultaneously, it uses CUDA atomic operations to ensure thread safety.

Note

Energy is binned logarithmically: \( i_E \approx \frac{\ln(E) - \ln(E_0)}{\Delta \ln E} \).

Note

Polar angle \( \theta \) is binned into angular zones, mirrored across the equator.

Parameters:

  • ph – The Structure of Arrays (SOA) containing all photon data.

  • d_spect – Pointer to the global spectrum structure array.

  • photon_index – The specific index of the photon being recorded.

Returns:

void

__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.