Updating the main branch in Daniel's root repo with my changes after completing 2606.30729

.gitignore

@@ -3,6 +3,11 @@
 ##
 ## Get latest from https://github.com/github/gitignore/blob/master/VisualStudio.gitignore
 
+# Additional
+FOORT/cfgs/*
+FOORT/Output/*
+Benchmarking/Seppe/*
+
 # FOORT compilation files
 FOORT/*.o
 FOORT/FOORT

FOORT/src/CMakeLists.txt

@@ -7,6 +7,8 @@ add_library(config_reader_lib STATIC ConfigReader.cpp ConfigReader.h)
 add_library(diagnostics_lib STATIC Diagnostics.cpp Diagnostics.h Geometry.h)
 add_library(diagnostics_emission_lib STATIC DiagnosticsEmission.cpp DiagnosticsEmission.h Geometry.h)
 add_library(geodesic_lib STATIC Geodesic.cpp Geodesic.h Geometry.h)
+add_library(grid_lib STATIC Grid.cpp Grid.h)
+add_library(interpolator_lib STATIC Interpolator.cpp Interpolator.h)
 add_library(input_output_lib STATIC InputOutput.cpp InputOutput.h Geometry.h)
 add_library(integrators_lib STATIC Integrators.cpp Integrators.h Geometry.h)
 add_library(mesh_lib STATIC Mesh.cpp Mesh.h Geometry.h)
@@ -24,7 +26,7 @@ target_link_libraries(geodesic_lib PUBLIC diagnostics_lib input_output_lib integ
 target_link_libraries(input_output_lib stdc++fs)
 target_link_libraries(integrators_lib PUBLIC geodesic_lib metric_lib)
 target_link_libraries(mesh_lib PUBLIC diagnostics_lib input_output_lib utilities_lib)
-target_link_libraries(metric_lib PUBLIC input_output_lib integrators_lib spline_lib)
+target_link_libraries(metric_lib PUBLIC grid_lib input_output_lib integrators_lib interpolator_lib spline_lib)
 target_link_libraries(terminations_lib PUBLIC geodesic_lib input_output_lib)
 target_link_libraries(utilities_lib PUBLIC diagnostics_lib geodesic_lib integrators_lib metric_lib terminations_lib viewscreen_lib)
 target_link_libraries(viewscreen_lib PUBLIC metric_lib mesh_lib)

FOORT/src/Config.cpp

@@ -294,6 +294,28 @@ std::unique_ptr<Metric> Config::GetMetric(const ConfigCollection &theCfg)
 			// All settings complete; create Metric object!
 			TheMetric = std::unique_ptr<Metric>(new BosonStarMetric(Phi_infinity, num_lines, rLogScale, Phi_filename, m_filename));
 		}
+		else if (MetricName == "rotatingbosonstar")
+		{
+			// The rotating boson star with solitonic potential
+
+			// First setting to look up: using a logarithmic r coordinate or not.
+			// Don't need to output message if setting not found
+			bool rLogScale{false};
+			bool flipAngularMomentum{false};
+			std::string MetricFolder{"RotatingBosonStar/data_Will/"};
+			int NumX{500};
+			int NumTh{399};
+			real L{1.};
+			MetricSettings.LookupValue("RLogScale", rLogScale);
+			MetricSettings.LookupValue("FlipAngularMomentum", flipAngularMomentum);
+			MetricSettings.LookupValue("MetricFolder", MetricFolder);
+			MetricSettings.LookupValue("NumX", NumX);
+			MetricSettings.LookupValue("NumTh", NumTh);
+			MetricSettings.LookupValue("L", L);
+
+			// All settings complete; create Metric object!
+			TheMetric = std::unique_ptr<Metric>(new RotatingBosonStarMetric(rLogScale, MetricFolder, NumX, NumTh, L, flipAngularMomentum));
+		}
 		//// METRIC ADD POINT B ////
 		// Add an else if clause to check for your new Metric object!
 		// To look for additional options in the metric configuration, use
@@ -602,26 +624,32 @@ void Config::InitializeDiagnostics(const ConfigCollection &theCfg, DiagBitflag &
 
 			//// Fluid four-velocity model selection and initialization ////
 
-			// Default fluid model and default parameters
-			real defaultxi{1.0};
-			real defaultbetar{1.0};
-			real defaultbetaphi{1.0};
-			std::unique_ptr<FluidVelocityModel> theFluidModel{new GeneralCircularRadialFluid(defaultxi, defaultbetar, defaultbetaphi, theMetric)};
+			// Fluid four-velocity model parameters (with sensible defaults)
+			real subKeplerianparam{1.0};
+			real betaR{1.0};
+			real betaPhi{1.0};
+			real iscolowerbound{0.05};
+			real iscoupperbound{1000.0};
 
 			// Read in fluid velocity model
-			std::string fluidmodelstring{""};
+			std::string fluidmodelstring{"GeneralCircularRadial"};
 			AllDiagSettings["EquatorialEmission"].LookupValue("FluidVelocityModel", fluidmodelstring);
 			if (fluidmodelstring == "GeneralCircularRadial")
 			{
-				real subKeplerianparam{defaultxi};
-				real betaR{defaultbetar};
-				real betaPhi{defaultbetaphi};
 				AllDiagSettings["EquatorialEmission"].LookupValue("xi", subKeplerianparam);
 				AllDiagSettings["EquatorialEmission"].LookupValue("betar", betaR);
 				AllDiagSettings["EquatorialEmission"].LookupValue("betaphi", betaPhi);
-
-				theFluidModel = std::unique_ptr<FluidVelocityModel>{new GeneralCircularRadialFluid(subKeplerianparam, betaR, betaPhi, theMetric)};
+				AllDiagSettings["EquatorialEmission"].LookupValue("ISCOLowerBound", iscolowerbound);
+				AllDiagSettings["EquatorialEmission"].LookupValue("ISCOUpperBound", iscoupperbound);
 			}
+			else
+			{
+				ScreenOutput("Unknown FluidVelocityModel \"" + fluidmodelstring + "\". Using default GeneralCircularRadial model.",
+							 Output_Other_Default);
+			}
+
+			std::unique_ptr<FluidVelocityModel> theFluidModel{
+				new GeneralCircularRadialFluid(subKeplerianparam, betaR, betaPhi, theMetric, iscolowerbound, iscoupperbound)};
 			// Other fluid velocity models can be checked for here...
 
 			// Set EquatorialEmissionDiagnostic options struct

FOORT/src/DiagnosticsEmission.cpp

@@ -180,7 +180,7 @@ std::string GeneralCircularRadialFluid::getFullDescriptionStr() const
 	if (m_ISCOexists && m_theMetric->getrLogScale())
 		trueISCOradius = exp(m_ISCOr);
 
-	return "Circular/radial flow (sub-Keplerian parameter xi = " + std::to_string(m_subKeplerParam) + ", beta_r = " + std::to_string(m_betaR) + ", beta_phi = " + std::to_string(m_betaPhi) + "; " + (m_ISCOexists ? "ISCO = " + std::to_string(trueISCOradius) : "no ISCO found") + ")";
+	return "Circular/radial flow (sub-Keplerian parameter xi = " + std::to_string(m_subKeplerParam) + ", beta_r = " + std::to_string(m_betaR) + ", beta_phi = " + std::to_string(m_betaPhi) + "; " + (m_ISCOexists ? "ISCO = " + std::to_string(trueISCOradius) : "no ISCO found") + "; " + "ISCO lower bound = " + std::to_string(m_ISCOlowerbound) + ", ISCO upper bound = " + std::to_string(m_ISCOupperbound) + ")";
 }
 
 /**
@@ -352,14 +352,20 @@ OneIndex GeneralCircularRadialFluid::GetRadialVelocityd(const Point &p) const
  */
 void GeneralCircularRadialFluid::FindISCO()
 {
-	real lowerbound{0.0};
-	real upperbound{1000.0};
+	real lowerbound{m_ISCOlowerbound};
+	real upperbound{m_ISCOupperbound};
 	const SphericalHorizonMetric *sphermetric = dynamic_cast<const SphericalHorizonMetric *>(m_theMetric);
 	if (sphermetric)
 	{
 		lowerbound = sphermetric->getrLogScale() ? log(sphermetric->getHorizonRadius()) : sphermetric->getHorizonRadius();
 		upperbound = sphermetric->getrLogScale() ? log(10.0 * sphermetric->getHorizonRadius()) : 10.0 * sphermetric->getHorizonRadius();
 	}
+	else if (m_theMetric->getrLogScale())
+	{
+		// If not a spherical horizon metric, but we are using log(r) coordinates, then set the lower bound to 0.0
+		lowerbound = log(m_ISCOlowerbound); // ln(0.05) = -3.912023005428146
+		upperbound = log(m_ISCOupperbound); // ln(1000.0) = 6.907755278982137
+	}
 
 	// Perform binary search for ISCO
 	// Allow only 1000 iterations max

FOORT/src/DiagnosticsEmission.h

@@ -95,8 +95,10 @@ struct FluidVelocityModel
 struct GeneralCircularRadialFluid final : public FluidVelocityModel
 {
 	// Constructor with three parameters and Metric pointer (which is passed to base class constructor)
-	GeneralCircularRadialFluid(real subKeplerParam, real betar, real betaphi, const Metric *const theMetric) : m_subKeplerParam{fmin(fmax(subKeplerParam, 0.0), 1.0)}, m_betaR{fmin(fmax(betar, 0.0), 1.0)},
-																											   m_betaPhi{fmin(fmax(betaphi, 0.0), 1.0)}, FluidVelocityModel(theMetric)
+	GeneralCircularRadialFluid(real subKeplerParam, real betar, real betaphi, const Metric *const theMetric,
+							   real ISCO_lowerbound, real ISCO_upperbound) : m_subKeplerParam{fmin(fmax(subKeplerParam, 0.0), 1.0)}, m_betaR{fmin(fmax(betar, 0.0), 1.0)},
+																			 m_betaPhi{fmin(fmax(betaphi, 0.0), 1.0)}, FluidVelocityModel(theMetric),
+																			 m_ISCOlowerbound{ISCO_lowerbound}, m_ISCOupperbound{ISCO_upperbound}
 	{
 		// Do some checks on three params, which must lie between 0.0 and 1.0 (note that they are adjusted as such in
 		// initializer above)
@@ -150,6 +152,10 @@ struct GeneralCircularRadialFluid final : public FluidVelocityModel
 	bool m_ISCOexists{false};
 	//! ISCO radius
 	real m_ISCOr{-1.0};
+	//! Lower bound for ISCO radius search
+	real m_ISCOlowerbound;
+	//! Upper bound for ISCO radius search
+	real m_ISCOupperbound;
 	//! ISCO t momentum
 	real m_ISCOpt{};
 	//! ISCO phi momentum

FOORT/src/Grid.cpp

@@ -0,0 +1,46 @@
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include "Grid.h"
+
+//! A class for a 2D grid that contains the values of the necessary functions.
+Grid::Grid(int N_row, int N_col)
+{
+    this->N_row = N_row;
+    this->N_col = N_col;
+    this->data = new double[N_row * N_col];
+    size = N_row * N_col;
+}
+
+void Grid::initialize_from_file(std::string file)
+{
+    std::ifstream inputFile(file);
+
+    if (!inputFile)
+    {
+        std::cerr << "Error opening file!" << std::endl;
+    }
+
+    std::string line;
+    int i = 0;
+    int j = 0;
+
+    while (std::getline(inputFile, line))
+    {
+        std::istringstream iss(line);
+
+        double val;
+        while (iss >> val)
+        {
+            data[i * N_col + j] = val;
+            j += 1;
+        }
+        j = 0;
+        i += 1;
+    }
+    std::cout << "Grid initialized from " << file << std::endl
+              << std::endl;
+
+    // close the file
+    inputFile.close();
+}

FOORT/src/Grid.h

@@ -0,0 +1,32 @@
+#include <fstream>
+#include <iostream>
+
+#ifndef GRID_H
+#define GRID_H
+
+//! A class for a 2D grid that contains the values of the necessary functions.
+class Grid
+{
+public:
+    //! Number of rows in the grid
+    int N_row;
+    //! Number of columns in the grid
+    int N_col;
+    //! Pointer to the data
+    double *data;
+    //! Number of rows times the number of columns
+    int size;
+    Grid(int N_row, int N_col);
+
+    //! Destructor
+    ~Grid() { delete[] this->data; }
+
+    //! Overload the () operator to access the data
+    double &operator()(int i, int j) { return this->data[i * N_col + j]; }
+    //! Overload the () operator to access the data (const version)
+    const double &operator()(int i, int j) const { return this->data[i * N_col + j]; }
+
+    void initialize_from_file(std::string file);
+};
+
+#endif // GRID_H