fix: repaired bug in thin_disk model that caused it to miss the disk

.gitignore

@@ -19,3 +19,6 @@ test-resources
 *.ppm
 *.mp4
 *.vtk
+
+# vscode
+.vscode/*

LICENSE

@@ -0,0 +1,28 @@
+BSD 3-Clause License
+
+Copyright (c) 2024, AFD Group at UIUC
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+   contributors may be used to endorse or promote products derived from
+   this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

README.md

@@ -22,11 +22,17 @@ $ sudo dnf install hdf5-devel hdf5-static gsl-devel
 ```
 The packages are straightforward to compile from source as well.
 
+To install ipole on macOS using ```gcc``` via homebrew, install 
+```bash
+$ brew install gcc@13 hdf5 gsl
+```
+
 The native macOS version of ```clang``` does not support OpenMP parallelization, so on macOS with ```brew```, install
 ```bash
 $ brew install llvm gsl hdf5
 ```
 these should all be binary packages, no compiling should be required.
+For the clang installation, change the ```HDF5_CC```, ```HDF5_CLINKER``` and ```GSL_DIR``` variables in the makefile. Please raise an issue in the event that neither of these methods work on macOS.
 
 # Building
 
@@ -139,11 +145,11 @@ polarized transfer equation assuming constant coefficients. This
 is intended to produce sensible results even when absorption or Faraday rotation
 and conversion are large in a single step.  A comparison of
 polarized relativistic radiative transport schemes can be found in
-Prather (2021, in prep).
+Prather [2023](https://ui.adsabs.harvard.edu/abs/2023ApJ...950...35P/abstract).
 
 If you are imaging GRMHD simulations, you will find guidance on how a GRMHD
 simulation and analysis pipeline fits together in the PATOKA pipeline paper,
-Wong (2021, in prep).
+Wong [2022](https://ui.adsabs.harvard.edu/abs/2022ApJS..259...64W/abstract).
 
 ```ipole``` treats synchrotron emission and absorption but not bremsstrahlung or Compton
 scattering.  The transfer coefficients (emissivities, absorptivities, and rotativities)

machines/callisto.make

@@ -0,0 +1,5 @@
+GSL_DIR = /opt/sns/gsl/gcc-4.8.5/2.5/
+
+CFLAGS_CUSTOM += -I$(GSL_DIR)/include
+SYSTEM_LIBDIR += -L$(GSL_DIR)/lib
+MATH_LIB += -lgsl -lgslcblas -lm

machines/delta.make

@@ -0,0 +1,34 @@
+# Before building, make sure you have got the required modules,
+# and have set the appropriate paths
+#
+# 1. Load gsl.
+# `module load gsl`
+#
+# 2. Update LD_LIBRARY_PATH so that ipole knows where to look for shared libraries.
+# `export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/sw/spack/deltas11-2023-03/apps/linux-rhel8-zen3/gcc-11.4.0/gsl-2.7.1-ytg74v2/lib`
+#
+# 3. Download and compile your own HDF5 (Delta's HDF5 is perpetually broken).
+# `mkdir -p ~/software/hdf5; cd ~/software/hdf5`
+# `wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.12/hdf5-1.12.2/src/hdf5-1.12.2.tar.gz`
+# `tar -xvzf hdf5-1.12.2.tar.gz`
+# `cd hdf5-1.12.2`
+# `module load gcc`
+# `module load cmake`
+# `mkdir build; cd build`
+# `cmake -DCMAKE_INSTALL_PREFIX=~/software/hdf5/install -DBUILD_SHARED_LIBS=ON ..`
+# `make -j$(nproc)`
+# `make install`
+# 
+# 4. Set environment variables for custom HDF5
+# `export HDF5_DIR=$HOME/software/hdf5/install`
+# `export PATH=$HDF5_DIR/bin:$PATH`
+# `export LD_LIBRARY_PATH=$HDF5_DIR/lib:$LD_LIBRARY_PATH`
+# `export CPATH=$HDF5_DIR/include:$CPATH`
+
+CC=h5cc
+GSL_DIR=$(GSL_HOME)
+#assumes you have loaded cray-hdf5 module which automatically sets HDF5_DIR
+HDF5_DIR=$(HDF5_ROOT)
+
+CFLAGS = -march=native -mtune=native -std=gnu11 -O3 -flto -fopenmp -funroll-loops
+$(info Using Delta settings: remember to prepend GSL lib location $(GSL_HOME)/lib to LD_LIBRARY_PATH for runtime)

makefile

@@ -7,7 +7,8 @@ HDF5_DIR =
 GSL_DIR =
 # System /lib equivalent (can be /usr/lib, /lib64, /usr/lib64)
 # Can leave this blank if it's included automatically by GCC
-SYSTEM_LIBDIR = /lib64
+#SYSTEM_LIBDIR = /lib64
+SYSTEM_LIBDIR =
 
 # Try pointing this to h5pcc or h5cc on your machine, before hunting down libraries
 CC=h5cc
@@ -27,10 +28,10 @@ ECHO=echo -e
 
 # Overrides of the above for macOS
 ifneq (,$(findstring Darwin,$(shell uname)))
-	export HDF5_CC = /usr/local/opt/llvm/bin/clang
-	export HDF5_CLINKER = /usr/local/opt/llvm/bin/clang
+	export HDF5_CC = /opt/homebrew/bin/gcc-13
+	export HDF5_CLINKER = /opt/homebrew/bin/gcc-13
 
-	GSL_DIR=/usr/local
+	GSL_DIR=/opt/homebrew/opt/gsl
 	SYSTEM_LIBDIR=
 
 	MD5=md5
@@ -47,6 +48,12 @@ endif
 ifneq (,$(findstring frontera,$(HOST)))
 	-include $(MAKEFILE_PATH)/machines/frontera.make
 endif
+ifneq (,$(findstring delta,$(HOST)))
+	-include $(MAKEFILE_PATH)/machines/delta.make
+endif
+ifneq (,$(findstring callisto,$(HOST)))
+	-include $(MAKEFILE_PATH)/machines/callisto.make
+endif
 # Hack to check only whether host begins with bh*
 ifneq (,$(findstring beginsbh,begins$(HOST)))
         -include $(MAKEFILE_PATH)/machines/bh-cluster.make

model/iharm/example.par

@@ -77,3 +77,7 @@ trace 0
 trace_stride 1
 # Trace file name
 trace_outf trace.h5
+
+# Histogram of origin of observed emissivity 
+histo 0
+histo_outf histo.h5

model/iharm/model.c

@@ -46,7 +46,7 @@ static double mu_i, mu_e, mu_tot;
 
 // MODEL PARAMETERS: PUBLIC
 double DTd;
-double sigma_cut = 1.0;
+double sigma_cut = 1.;
 double beta_crit = 1.0;
 double sigma_cut_high = -1.0;
 
@@ -71,7 +71,8 @@ static double Ladv_dump;
 
 static int reverse_field = 0;
 
-double tf;
+// Need to share this with main.c!
+extern double tf;
 
 // MAYBES
 //static double t0;
@@ -85,7 +86,7 @@ double tf;
 // TODO the way this is selected is horrid.  Make it a parameter.
 #define ELECTRONS_TFLUID (3)
 static int RADIATION, ELECTRONS;
-static double gam = 1.444444, game = 1.333333, gamp = 1.666667;
+static double gam = 1.444444444444444, game = 1.333333333333333, gamp = 1.666666666666667;
 static double Thetae_unit, Mdotedd;
 
 // Ignore radiation interactions within one degree of polar axis
@@ -96,6 +97,10 @@ static int nloaded = 0;
 
 static hdf5_blob fluid_header = { 0 };
 
+
+// Debug KHARMA reader
+#define DEBUG_READER (0)
+
 struct of_data {
   double t;
   double ****p;
@@ -328,6 +333,59 @@ void init_model(double *tA, double *tB)
       th_beg = 0.0174 * 2;
     }
   }
+
+  #if DEBUG_READER
+    /* Set filename */
+    char debug_fname[256];
+    snprintf(debug_fname, sizeof(debug_fname), "debug_reader_iharm.h5");
+
+    /* Create HDF5 file*/
+    hdf5_create(debug_fname);
+
+    /* Compute gcov, gcon */
+    size_t total_elements = NDIM * NDIM * (N2 + 2) * (N1 + 2);
+    double *gcov_global = malloc(total_elements * sizeof(double));
+    double *gcon_global = malloc(total_elements * sizeof(double));
+    // Use the arrays via indexing. For example, to access element [mu][nu][j][i]:
+    #define IDX(mu, nu, j, i) (((mu) * NDIM * (N2+2) * (N1+2)) + ((nu) * (N2+2) * (N1+2)) + ((j) * (N1+2)) + (i))
+
+#pragma omp parallel for collapse(2)
+  for (int i = 0; i < N1+2; i++) {
+    for (int j = 0; j < N2+2; j++) {
+
+      double X[NDIM] = {0.};
+      ijktoX(i, j, 0, X);
+      double gcov[NDIM][NDIM], gcon[NDIM][NDIM];
+      gcov_func(X, gcov);
+      gcon_func(gcov, gcon);
+
+      for (int mu = 0; mu < NDIM; mu++) {
+        for (int nu = 0; nu < NDIM; nu++) {
+          gcov_global[IDX(mu, nu, j, i)] = gcov[mu][nu];
+          gcon_global[IDX(mu, nu, j, i)] = gcon[mu][nu];
+        }
+      }
+
+    }
+  }
+
+    /* Write gcov, gcon to file */
+    hsize_t dims_grid[4] = { NDIM, NDIM, N2+2, N1+2 };
+    hdf5_write_full_array(gcov_global, "gcov", 4, dims_grid, H5T_NATIVE_DOUBLE);
+    hdf5_write_full_array(gcon_global, "gcon", 4, dims_grid, H5T_NATIVE_DOUBLE);
+
+    /* Write physical quantities */
+    hsize_t dims_phys[3] = { N1+2, N2+2, N3+2 };
+    hdf5_write_full_array(data[0]->ne[0][0], "ne", 3, dims_phys, H5T_NATIVE_DOUBLE);
+    hdf5_write_full_array(data[0]->thetae[0][0], "thetae", 3, dims_phys, H5T_NATIVE_DOUBLE);
+    hdf5_write_full_array(data[0]->b[0][0], "b", 3, dims_phys, H5T_NATIVE_DOUBLE);
+    hdf5_write_full_array(data[0]->sigma[0][0], "sigma", 3, dims_phys, H5T_NATIVE_DOUBLE);
+    hdf5_write_full_array(data[0]->beta[0][0], "beta", 3, dims_phys, H5T_NATIVE_DOUBLE);
+
+    /* Close HDF5 file */
+    hdf5_close();
+
+  #endif // DEBUG_READER
 }
 
 /*
@@ -566,6 +624,8 @@ void init_physical_quantities(int n, double rescale_factor)
 
         data[n]->b[i][j][k] *= rescale_factor;
 
+        // here b is in gauss; this takes it back to code units to 
+        // calculate beta and sigma
         double bsq = data[n]->b[i][j][k] / B_unit;
         bsq = bsq*bsq;
 
@@ -862,7 +922,7 @@ void init_iharm_grid(char *fnam, int dumpidx)
   hdf5_read_single_val(&N3, "n3", H5T_STD_I32LE);
   hdf5_read_single_val(&gam, "gam", H5T_IEEE_F64LE);
 
-  if (hdf5_exists("gam_e")) {
+  if (hdf5_exists("gam_e")) { /* pyharm-converted files save gam_e, gam_p even when electrons are not run. Which overwrites the default values*/
     fprintf(stderr, "custom electron model loaded from dump file...\n");
     hdf5_read_single_val(&game, "gam_e", H5T_IEEE_F64LE);
     hdf5_read_single_val(&gamp, "gam_p", H5T_IEEE_F64LE);

model/kharma/example.par

@@ -0,0 +1,83 @@
+# example.par -- ipole example parameter file for iharm problem
+
+# Common
+dump tests/test-resources/sample_dump_SANE_a+0.94_MKS_0900.h5
+outfile image.h5
+nx 160
+ny 160
+
+# Emission 1=Pandya, 4=Dexter, others are debug/custom
+emission_type 4
+
+# M87 parameters
+MBH 6.2e9
+dsource 16.9e6
+freqcgs 230.e9
+# M_unit scales the accretion rate to match a known object flux
+# These are example values
+# SANE:
+M_unit 3e26
+# MAD:
+#M_unit 6e24
+
+# e- Temperature, via the Rhigh model described 
+trat_small 1
+trat_large 20
+# Constant e- temperature ratio
+#tp_over_te 3
+
+# Adaptive res
+# enable by setting a minimum "base" image size
+#nx_min 40
+#ny_min 40
+refine_abs = 2e-2
+refine_rel = 1e-2
+refine_cut = 0
+
+# Camera
+rcam 1000
+# Values in degrees
+thetacam 17
+phicam 0
+rotcam 0
+
+# FOV from Earth
+fovx_dsource 160
+fovy_dsource 160
+
+# Geodesic accuracy parameters
+#eps 0.01
+#maxnstep 50000
+
+# Options
+# Convention for EVPA defining stokes Q,U:
+# 0 is measured East of North, 1 is North of West
+qu_conv 0
+# Don't produce an output file
+quench_output 0
+# Add a .ppm image of the unpolarized flux
+add_ppm 0
+# Only calculate the unpolarized image
+only_unpolarized 0
+# 1 to emit only for th>PI/2, 2 for other hemisphere
+# 0 emits everywhere
+counterjet 0
+# Offset for each geodesic in pixels, used to prevent
+# overtraining ML models to pixel locations
+xoff 0
+yoff 0
+
+# Path trace -- save emissivities & local state
+# for every step along a geodesic, or several
+trace 0
+# Pixel to trace i,j (rightward from left, upward from bottom)
+#trace_i 40
+#trace_j 40
+# Or trace every N pixels in each direction
+trace_stride 1
+# Trace file name
+trace_outf trace.h5
+
+# Histogram of origin of observed emissivity 
+histo 0
+histo_outf histo.h5