All COCAL Reader Implemention

CMakeLists.txt

@@ -170,7 +170,7 @@ include(UpdateAddExecutables)
 # dependencies on the PCH if necessary.
 include(SetupClangFormat)
 include(SetupClangTidy)
-include(SetupIwyu)
+#include(SetupIwyu) // This somehow doesn't work for COCAL reader
 if(BUILD_DOCS)
   include(SetupDoxygen)
   include(SetupSphinx)

cmake/FindCOCAL.cmake

@@ -9,39 +9,45 @@
 # - COCAL::Exporter: Functionality to load COCAL volume data and interpolate to
 #   arbitrary points.
 
-if(NOT COCAL_ROOT)
-  set(COCAL_ROOT "")
+if(NOT DEFINED COCAL_ROOT OR COCAL_ROOT STREQUAL "")
   set(COCAL_ROOT $ENV{COCAL_ROOT})
 endif()
 
-find_library(
-  COCAL_LIB
-  NAMES libcocal.a
-  PATHS ${COCAL_ROOT}/lib
-  NO_DEFAULT_PATHS
-)
+if (DEFINED COCAL_ROOT AND NOT COCAL_ROOT STREQUAL "")
+
+  # Find the COCAL library
+  find_library(
+    COCAL_LIB
+    NAMES libcocal.a
+    PATHS ${COCAL_ROOT}/lib
+    NO_DEFAULT_PATHS
+  )
+
+  # Find the COCAL include directory
+  find_path(
+    COCAL_INCLUDE_DIR
+    NAMES coc2cac_bns.f90
+    PATHS ${COCAL_ROOT}/include
+    NO_DEFAULT_PATHS
+  )
 
-find_path(
-  COCAL_INCLUDE_DIR
-  NAMES coc2cac_bns.f90
-  PATHS ${COCAL_ROOT}/include
-  NO_DEFAULT_PATHS
-)
+  # Find MPI
+  find_package(MPI COMPONENTS CXX)
 
-# Optionally link MPI if COCAL is built with MPI dependencies
-find_package(MPI COMPONENTS CXX)
+  if (COCAL_LIB AND MPI_CXX_FOUND)
+    add_library(COCAL::Exporter INTERFACE IMPORTED)
+    set_target_properties(COCAL::Exporter PROPERTIES
+      INTERFACE_INCLUDE_DIRECTORIES "${COCAL_INCLUDE_DIR}"
+      INTERFACE_LINK_LIBRARIES "${COCAL_LIB};MPI::MPI_CXX"
+    )
+  endif()
 
-if (COCAL_LIB AND MPI_CXX_FOUND)
-  add_library(COCAL::Exporter INTERFACE IMPORTED)
-  set_target_properties(COCAL::Exporter PROPERTIES
-    INTERFACE_INCLUDE_DIRECTORIES "${COCAL_INCLUDE_DIR}"
-    INTERFACE_LINK_LIBRARIES "${COCAL_LIB};MPI::MPI_CXX;${OtherDependency_LIBRARIES}"
+  include(FindPackageHandleStandardArgs)
+  find_package_handle_standard_args(
+    COCAL
+    FOUND_VAR COCAL_FOUND
+    REQUIRED_VARS COCAL_LIB COCAL_INCLUDE_DIR
   )
+else()
+  set(COCAL_FOUND FALSE)
 endif()
-
-include(FindPackageHandleStandardArgs)
-find_package_handle_standard_args(
-  COCAL
-  FOUND_VAR COCAL_FOUND
-  REQUIRED_VARS COCAL_LIB COCAL_INCLUDE_DIR
-)

cmake/SetupCocal.cmake

@@ -1,12 +1,15 @@
 # Distributed under the MIT License.
 # See LICENSE.txt for details.
 
-find_package(COCAL REQUIRED)
+
+find_package(COCAL)
 
 if (COCAL_FOUND)
-  message("Found cocal")
+  message("-- Found Cocal at ${COCAL_ROOT}")
   file(APPEND
     "${CMAKE_BINARY_DIR}/BuildInfo.txt"
     "COCAL: ${COCAL_ROOT}\n"
     )
+else()
+  message("-- Could not find Cocal")
 endif()

src/IO/External/CMakeLists.txt

@@ -60,10 +60,8 @@ if (TARGET FUKA::Exporter)
     ${LIBRARY} INTERFACE HAS_FUKA_EXPORTER)
 endif()
 
-# Add COCAL interpolation support
-message("\n\n\n\n\nHERE\n\n\n\n\n")
+
 if (TARGET COCAL::Exporter)
-  message("\n\n\n\n\nHERE\n\n\n\n\n")
   spectre_target_sources(
     ${LIBRARY}
     PRIVATE

src/IO/External/InterpolateFromCocal.cpp

@@ -15,9 +15,39 @@
 #include "Utilities/TaggedTuple.hpp"
 
 extern "C" {
+// Define the COCAL Fortran reader functions here for different id_types.
+void coc2cac_co(const int& npoints, const double* x, const double* y,
+                const double* z, const char* data_directory, const int& len_dir,
+                double* lapse, double* shift_x, double* shift_y,
+                double* shift_z, double* spatial_metric_xx,
+                double* spatial_metric_xy, double* spatial_metric_xz,
+                double* spatial_metric_yy, double* spatial_metric_yz,
+                double* spatial_metric_zz, double* extrinsic_curvature_xx,
+                double* extrinsic_curvature_xy, double* extrinsic_curvature_xz,
+                double* extrinsic_curvature_yy, double* extrinsic_curvature_yz,
+                double* extrinsic_curvature_zz, double* rest_mass_density,
+                double* spatial_velocity_x, double* spatial_velocity_y,
+                double* spatial_velocity_z, double* pressure,
+                double* specific_internal_energy);
+
 void coc2cac_ir(const int& npoints, const double* x, const double* y,
-                const double* z, double* lapse, double* shift_x,
-                double* shift_y, double* shift_z, double* spatial_metric_xx,
+                const double* z, const char* data_directory, const int& len_dir,
+                double* lapse, double* shift_x, double* shift_y,
+                double* shift_z, double* spatial_metric_xx,
+                double* spatial_metric_xy, double* spatial_metric_xz,
+                double* spatial_metric_yy, double* spatial_metric_yz,
+                double* spatial_metric_zz, double* extrinsic_curvature_xx,
+                double* extrinsic_curvature_xy, double* extrinsic_curvature_xz,
+                double* extrinsic_curvature_yy, double* extrinsic_curvature_yz,
+                double* extrinsic_curvature_zz, double* rest_mass_density,
+                double* spatial_velocity_x, double* spatial_velocity_y,
+                double* spatial_velocity_z, double* pressure,
+                double* specific_internal_energy);
+
+void coc2cac_sp(const int& npoints, const double* x, const double* y,
+                const double* z, const char* data_directory, const int& len_dir,
+                double* lapse, double* shift_x, double* shift_y,
+                double* shift_z, double* spatial_metric_xx,
                 double* spatial_metric_xy, double* spatial_metric_xz,
                 double* spatial_metric_yy, double* spatial_metric_yz,
                 double* spatial_metric_zz, double* extrinsic_curvature_xx,
@@ -48,11 +78,12 @@ DataVector to_datavector(std::vector<double> vec) {
 }  // namespace
 
 tuples::tagged_tuple_from_typelist<cocal_tags> interpolate_from_cocal(
-    const gsl::not_null<std::mutex*> cocal_lock,
-    /*const std::string& data_directory,*/
+    const gsl::not_null<std::mutex*> cocal_lock, const CocalIdType id_type,
+    const std::string& data_directory,
     const tnsr::I<DataVector, 3, Frame::Inertial>& x) {
   tuples::tagged_tuple_from_typelist<cocal_tags> result{};
   const std::lock_guard lock{*cocal_lock};
+  const ScopedFpeState disable_fpes(false);  // test for fpe issues
   const size_t num_points = get<0>(x).size();
   std::vector<double> x_coords(num_points);
   std::vector<double> y_coords(num_points);
@@ -88,60 +119,52 @@ tuples::tagged_tuple_from_typelist<cocal_tags> interpolate_from_cocal(
   std::vector<double> pressure(num_points);
   std::vector<double> specific_internal_energy(num_points);
 
-  // std::cout << "Calling coc2cac_ir " /* << num_points << " points." */ <<
-  // std::endl;
-
-  coc2cac_ir(num_points, x_coords.data(), y_coords.data(), z_coords.data(),
-             lapse.data(), shift_x.data(), shift_y.data(), shift_z.data(),
-             spatial_metric_xx.data(), spatial_metric_xy.data(),
-             spatial_metric_xz.data(), spatial_metric_yy.data(),
-             spatial_metric_yz.data(), spatial_metric_zz.data(),
-             extrinsic_curvature_xx.data(), extrinsic_curvature_xy.data(),
-             extrinsic_curvature_xz.data(), extrinsic_curvature_yy.data(),
-             extrinsic_curvature_yz.data(), extrinsic_curvature_zz.data(),
-             rest_mass_density.data(), spatial_velocity_x.data(),
-             spatial_velocity_y.data(), spatial_velocity_z.data(),
-             pressure.data(), specific_internal_energy.data());
-
-  get(get<gr::Tags::Lapse<DataVector>>(result)) = DataVector(num_points);
-  get<gr::Tags::Shift<DataVector, 3>>(result).get(0) = DataVector(num_points);
-  get<gr::Tags::Shift<DataVector, 3>>(result).get(1) = DataVector(num_points);
-  get<gr::Tags::Shift<DataVector, 3>>(result).get(2) = DataVector(num_points);
-  get<gr::Tags::SpatialMetric<DataVector, 3>>(result).get(0, 0) =
-      DataVector(num_points);
-  get<gr::Tags::SpatialMetric<DataVector, 3>>(result).get(0, 1) =
-      DataVector(num_points);
-  get<gr::Tags::SpatialMetric<DataVector, 3>>(result).get(0, 2) =
-      DataVector(num_points);
-  get<gr::Tags::SpatialMetric<DataVector, 3>>(result).get(1, 1) =
-      DataVector(num_points);
-  get<gr::Tags::SpatialMetric<DataVector, 3>>(result).get(1, 2) =
-      DataVector(num_points);
-  get<gr::Tags::SpatialMetric<DataVector, 3>>(result).get(2, 2) =
-      DataVector(num_points);
-  get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(result).get(0, 0) =
-      DataVector(num_points);
-  get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(result).get(0, 1) =
-      DataVector(num_points);
-  get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(result).get(0, 2) =
-      DataVector(num_points);
-  get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(result).get(1, 1) =
-      DataVector(num_points);
-  get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(result).get(1, 2) =
-      DataVector(num_points);
-  get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(result).get(2, 2) =
-      DataVector(num_points);
-  get(get<hydro::Tags::RestMassDensity<DataVector>>(result)) =
-      DataVector(num_points);
-  get<hydro::Tags::SpatialVelocity<DataVector, 3>>(result).get(0) =
-      DataVector(num_points);
-  get<hydro::Tags::SpatialVelocity<DataVector, 3>>(result).get(1) =
-      DataVector(num_points);
-  get<hydro::Tags::SpatialVelocity<DataVector, 3>>(result).get(2) =
-      DataVector(num_points);
-  get(get<hydro::Tags::SpecificInternalEnergy<DataVector>>(result)) =
-      DataVector(num_points);
-  get(get<hydro::Tags::Pressure<DataVector>>(result)) = DataVector(num_points);
+  // Pass the string and its length to Fortran
+  const char* dir_cstr =
+      data_directory.c_str();  // Convert std::string to C string
+  int len_dir = static_cast<int>(
+      data_directory.length());  // Get the length of the string
+
+  if (id_type == CocalIdType::Co) {
+    coc2cac_co(num_points, x_coords.data(), y_coords.data(), z_coords.data(),
+               dir_cstr, len_dir, lapse.data(), shift_x.data(), shift_y.data(),
+               shift_z.data(), spatial_metric_xx.data(),
+               spatial_metric_xy.data(), spatial_metric_xz.data(),
+               spatial_metric_yy.data(), spatial_metric_yz.data(),
+               spatial_metric_zz.data(), extrinsic_curvature_xx.data(),
+               extrinsic_curvature_xy.data(), extrinsic_curvature_xz.data(),
+               extrinsic_curvature_yy.data(), extrinsic_curvature_yz.data(),
+               extrinsic_curvature_zz.data(), rest_mass_density.data(),
+               spatial_velocity_x.data(), spatial_velocity_y.data(),
+               spatial_velocity_z.data(), pressure.data(),
+               specific_internal_energy.data());
+  } else if (id_type == CocalIdType::Ir) {
+    coc2cac_ir(num_points, x_coords.data(), y_coords.data(), z_coords.data(),
+               dir_cstr, len_dir, lapse.data(), shift_x.data(), shift_y.data(),
+               shift_z.data(), spatial_metric_xx.data(),
+               spatial_metric_xy.data(), spatial_metric_xz.data(),
+               spatial_metric_yy.data(), spatial_metric_yz.data(),
+               spatial_metric_zz.data(), extrinsic_curvature_xx.data(),
+               extrinsic_curvature_xy.data(), extrinsic_curvature_xz.data(),
+               extrinsic_curvature_yy.data(), extrinsic_curvature_yz.data(),
+               extrinsic_curvature_zz.data(), rest_mass_density.data(),
+               spatial_velocity_x.data(), spatial_velocity_y.data(),
+               spatial_velocity_z.data(), pressure.data(),
+               specific_internal_energy.data());
+  } else if (id_type == CocalIdType::Sp) {
+    coc2cac_sp(num_points, x_coords.data(), y_coords.data(), z_coords.data(),
+               dir_cstr, len_dir, lapse.data(), shift_x.data(), shift_y.data(),
+               shift_z.data(), spatial_metric_xx.data(),
+               spatial_metric_xy.data(), spatial_metric_xz.data(),
+               spatial_metric_yy.data(), spatial_metric_yz.data(),
+               spatial_metric_zz.data(), extrinsic_curvature_xx.data(),
+               extrinsic_curvature_xy.data(), extrinsic_curvature_xz.data(),
+               extrinsic_curvature_yy.data(), extrinsic_curvature_yz.data(),
+               extrinsic_curvature_zz.data(), rest_mass_density.data(),
+               spatial_velocity_x.data(), spatial_velocity_y.data(),
+               spatial_velocity_z.data(), pressure.data(),
+               specific_internal_energy.data());
+  }
 
   get(get<gr::Tags::Lapse<DataVector>>(result)) = to_datavector(lapse);
   get<gr::Tags::Shift<DataVector, 3>>(result).get(0) = to_datavector(shift_x);

src/IO/External/InterpolateFromCocal.hpp

@@ -16,6 +16,9 @@
 
 namespace io {
 
+/// Type of COCAL initial data
+enum class CocalIdType { Co, Ir, Sp };
+
 // List of tags supplied by COCAL initial data
 using cocal_tags =
     tmpl::list<gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, 3>,
@@ -29,13 +32,15 @@ using cocal_tags =
 
 /*!
  * \brief Interpolate numerical COCAL initial data to arbitrary points
- * \param cocal_lock
+ * \param cocal_lock A simple lock to lock the thread
+ * \param id_type Type of COCAL initial data (e.g., Co, Ir, Sp)
+ * \param data_directory Directory containing COCAL data
  * \param x Coordinates of points to interpolate to
  * \return Data interpolated to the given points
  */
 tuples::tagged_tuple_from_typelist<cocal_tags> interpolate_from_cocal(
-    gsl::not_null<std::mutex*> cocal_lock,
-    /*const std::string& data_directory,*/
+    gsl::not_null<std::mutex*> cocal_lock, const CocalIdType id_type,
+    const std::string& data_directory,
     const tnsr::I<DataVector, 3, Frame::Inertial>& x);
 
 }  // namespace io

src/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt

@@ -94,6 +94,20 @@ if (TARGET FUKA::Exporter)
 endif()
 
 if (TARGET COCAL::Exporter)
+add_custom_command(
+  OUTPUT ${CMAKE_SOURCE_DIR}/src/PointwiseFunctions/AnalyticData/GrMhd/COCAL_reader/libcocal.a
+  COMMAND make -C ${CMAKE_SOURCE_DIR}/src/PointwiseFunctions/AnalyticData/GrMhd/COCAL_reader clean
+  COMMAND make -C ${CMAKE_SOURCE_DIR}/src/PointwiseFunctions/AnalyticData/GrMhd/COCAL_reader
+  COMMENT "Building COCAL library..."
+  VERBATIM
+)
+
+add_custom_target(BuildCocalLibrary
+  DEPENDS ${CMAKE_SOURCE_DIR}/src/PointwiseFunctions/AnalyticData/GrMhd/COCAL_reader/libcocal.a
+)
+
+add_dependencies(${LIBRARY} BuildCocalLibrary)
+
   spectre_target_sources(
     ${LIBRARY}
     PRIVATE

src/PointwiseFunctions/AnalyticData/GrMhd/COCAL_reader/IO_input_CF_flco_export.f90

@@ -0,0 +1,22 @@
+subroutine IO_input_CF_flco_export(filenm,emd,ome,ber,radi)
+  use phys_constant, only : long, nnrg, nntg, nnpg
+  implicit none
+  integer :: ir, it, ip, nrtmp, nttmp, nptmp
+  real(8), pointer :: emd(:,:,:)
+  real(8) :: ome, ber, radi
+  character(len=*) :: filenm
+!
+! --- Matter
+    open(12,file=trim(filenm),status='old')
+    read(12,'(5i5)')  nrtmp, nttmp, nptmp
+    do ip = 0, nptmp
+      do it = 0, nttmp
+        do ir = 0, nrtmp
+          read(12,'(1p,6e20.12)') emd(ir,it,ip)
+        end do
+      end do
+    end do
+    read(12,'(1p,6e20.12)') ome, ber, radi
+    close(12)
+!
+end subroutine IO_input_CF_flco_export

src/PointwiseFunctions/AnalyticData/GrMhd/COCAL_reader/IO_input_CF_flir_export.f90

@@ -0,0 +1,25 @@
+subroutine IO_input_CF_flir_export(filenm,emd,vep,ome,ber,radi)
+  use phys_constant, only : long, nnrg, nntg, nnpg
+  implicit none
+  integer :: ir, it, ip, nrtmp, nttmp, nptmp
+  real(8), pointer :: emd(:,:,:), vep(:,:,:), rs(:,:)
+  real(8) :: ome, ber, radi
+  character(len=*) :: filenm
+  character*400 :: tmp_file
+!
+! --- Matter
+!  tmp_file = "/home/fs01/ml2847/SPECTRE/BNS/ID/IRE3.0_SLy_SLy_010_M1.96/work_area_BNS/bnsflu_3D_mpt1.las"
+  open(12,file=trim(filenm),status='old')
+!  open(12,file=trim(tmp_file),status='old')
+  read(12,'(5i5)')  nrtmp, nttmp, nptmp
+  do ip = 0, nptmp
+    do it = 0, nttmp
+      do ir = 0, nrtmp
+        read(12,'(1p,6e20.12)') emd(ir,it,ip), vep(ir,it,ip)
+      end do
+    end do
+  end do
+  read(12,'(1p,6e20.12)') ome, ber, radi
+  close(12)
+!
+end subroutine IO_input_CF_flir_export