All cocal implementation

CMakeLists.txt

@@ -49,7 +49,6 @@ endif()
 if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.30)
   cmake_policy(SET CMP0167 NEW)
 endif ()
-
 # FetchContent can be used to download and build some dependencies. We can
 # consider turning this on by default later.
 option(SPECTRE_FETCH_MISSING_DEPS "Download missing dependencies" OFF)
@@ -76,6 +75,7 @@ include(SpectreSetSiteName)
 # We need Python in the build system and throughout the code. Setting it up
 # early so we can rely on a consistent Python version in the build system.
 find_package(Python 3.8 REQUIRED)
+# find_package(COCAL REQUIRED)
 
 # Define the location of Python code in the build directory. Unit tests etc. can
 # add this location to their `PYTHONPATH`.
@@ -112,6 +112,7 @@ include(ProhibitInSourceBuild)
 include(SpectreSetupFlagsTarget)
 include(SetupFortran)
 include(SetupNinjaColors)
+include(SetupCocal)  # Add this line
 include(SetOutputDirectory)
 include(SpectreAddInterfaceLibraryHeaders)
 include(SpectreTargetHeaders)
@@ -170,8 +171,6 @@ include(SetupParaView)
 
 add_subdirectory(external)
 
-include(SetupKokkos)
-
 include(SetupLIBCXXCharm)
 # The precompiled header must be setup after all libraries have been found
 if (USE_PCH)
@@ -190,6 +189,7 @@ include(UpdateAddExecutables)
 # dependencies on the PCH if necessary.
 include(SetupClangFormat)
 include(SetupClangTidy)
+#include(SetupIwyu) // This somehow doesn't work for COCAL reader
 if(BUILD_DOCS)
   include(SetupDoxygen)
   include(SetupSphinx)
@@ -216,6 +216,8 @@ include_directories(${CMAKE_SOURCE_DIR}/src)
 include_directories(SYSTEM ${CMAKE_BINARY_DIR}/src)
 
 add_subdirectory(src)
+# add_subdirectory(src/PointwiseFunctions/AnalyticData/GrMhd/COCAL_reader)
+# target_link_libraries(EvolveGhValenciaDivCleanBns PRIVATE COCALReader)
 add_subdirectory(support)
 if(BUILD_TESTING)
   add_subdirectory(tests)

cmake/FindCOCAL.cmake

@@ -0,0 +1,53 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+# Find the COCAL initial data code (assumed repository or documentation link)
+#
+# Pass `COCAL_ROOT` to the CMake build configuration to set up the following
+# targets:
+#
+# - COCAL::Exporter: Functionality to load COCAL volume data and interpolate to
+#   arbitrary points.
+
+if(NOT DEFINED COCAL_ROOT OR COCAL_ROOT STREQUAL "")
+  set(COCAL_ROOT $ENV{COCAL_ROOT})
+endif()
+
+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 MPI
+  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()
+
+  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()

cmake/SetupCocal.cmake

@@ -0,0 +1,15 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+
+find_package(COCAL)
+
+if (COCAL_FOUND)
+  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/Evolution/Systems/GrMhd/GhValenciaDivClean/AllSolutions.hpp

@@ -24,6 +24,14 @@ using FukaInitialDataList = tmpl::list<grmhd::AnalyticData::FukaInitialData>;
 using FukaInitialDataList = NoSuchType;
 #endif
 
+// Check if COCAL is linked and therefore we can load COCAL initial data
+#ifdef HAS_COCAL_EXPORTER
+#include "PointwiseFunctions/AnalyticData/GrMhd/CocalInitialData.hpp"
+using CocalInitialDataList = tmpl::list<grmhd::AnalyticData::CocalInitialData>;
+#else
+using CocalInitialDataList = NoSuchType;
+#endif
+
 namespace ghmhd::GhValenciaDivClean::InitialData {
 // These are solutions that can be used for analytic prescriptions
 using analytic_solutions_and_data_list = gh::ghmhd_solutions;
@@ -34,5 +42,7 @@ using initial_data_list = tmpl::flatten<tmpl::list<
         tmpl::conditional_t<std::is_same_v<SpecInitialDataList, NoSuchType>,
                             tmpl::list<>, SpecInitialDataList>,
         tmpl::conditional_t<std::is_same_v<FukaInitialDataList, NoSuchType>,
-                            tmpl::list<>, FukaInitialDataList>>>>>;
+                            tmpl::list<>, FukaInitialDataList>,
+        tmpl::conditional_t<std::is_same_v<CocalInitialDataList, NoSuchType>,
+                            tmpl::list<>, CocalInitialDataList>>>>>;
 }  // namespace ghmhd::GhValenciaDivClean::InitialData

src/IO/External/CMakeLists.txt

@@ -8,6 +8,10 @@ if (TARGET FUKA::Exporter)
   set(_LIB_TYPE "")
 endif()
 
+if (TARGET COCAL::Exporter)
+  set(_LIB_TYPE "")
+endif()
+
 add_spectre_library(${LIBRARY} ${_LIB_TYPE})
 
 if (TARGET SpEC::Exporter)
@@ -55,3 +59,32 @@ if (TARGET FUKA::Exporter)
   target_compile_definitions(
     ${LIBRARY} INTERFACE HAS_FUKA_EXPORTER)
 endif()
+
+
+if (TARGET COCAL::Exporter)
+  spectre_target_sources(
+    ${LIBRARY}
+    PRIVATE
+    InterpolateFromCocal.cpp
+    )
+  spectre_target_headers(
+    ${LIBRARY}
+    INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
+    HEADERS
+    InterpolateFromCocal.hpp
+    )
+  target_link_libraries(
+    ${LIBRARY}
+    PUBLIC
+    COCAL::Exporter
+    DataStructures
+    GeneralRelativity
+    Hydro
+    Utilities
+    PRIVATE
+    ErrorHandling
+    )
+  target_compile_definitions(
+    ${LIBRARY} INTERFACE HAS_COCAL_EXPORTER)
+endif()
+

src/IO/External/InterpolateFromCocal.cpp

@@ -0,0 +1,204 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "IO/External/InterpolateFromCocal.hpp"
+
+#include <array>
+#include <string>
+#include <vector>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+#include "Utilities/TMPL.hpp"
+#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, 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,
+                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);
+}
+
+namespace io {
+
+namespace {
+DataVector to_datavector(std::vector<double> vec) {
+  DataVector result(vec.size());
+  std::copy(vec.begin(), vec.end(), result.begin());
+  return result;
+}
+}  // namespace
+
+tuples::tagged_tuple_from_typelist<cocal_tags> interpolate_from_cocal(
+    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);
+  std::vector<double> z_coords(num_points);
+  for (size_t i = 0; i < num_points; i++) {
+    x_coords[i] = x.get(0)[i];
+    y_coords[i] = x.get(1)[i];
+    z_coords[i] = x.get(2)[i];
+  }
+
+  std::vector<double> lapse(num_points);
+  std::vector<double> shift_x(num_points);
+  std::vector<double> shift_y(num_points);
+  std::vector<double> shift_z(num_points);
+  std::vector<double> spatial_metric_xx(num_points);
+  std::vector<double> spatial_metric_xy(num_points);
+  std::vector<double> spatial_metric_xz(num_points);
+  std::vector<double> spatial_metric_yy(num_points);
+  std::vector<double> spatial_metric_yz(num_points);
+  std::vector<double> spatial_metric_zz(num_points);
+  std::vector<double> extrinsic_curvature_xx(num_points);
+  std::vector<double> extrinsic_curvature_xy(num_points);
+  std::vector<double> extrinsic_curvature_xz(num_points);
+  std::vector<double> extrinsic_curvature_yy(num_points);
+  std::vector<double> extrinsic_curvature_yz(num_points);
+  std::vector<double> extrinsic_curvature_zz(num_points);
+  std::vector<double> rest_mass_density(num_points);
+  std::vector<double> spatial_velocity_x(num_points);
+  std::vector<double> spatial_velocity_y(num_points);
+  std::vector<double> spatial_velocity_z(num_points);
+  std::vector<double> pressure(num_points);
+  std::vector<double> specific_internal_energy(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);
+  get<gr::Tags::Shift<DataVector, 3>>(result).get(1) = to_datavector(shift_y);
+  get<gr::Tags::Shift<DataVector, 3>>(result).get(2) = to_datavector(shift_z);
+  get<gr::Tags::SpatialMetric<DataVector, 3>>(result).get(0, 0) =
+      to_datavector(spatial_metric_xx);
+
+  get<gr::Tags::SpatialMetric<DataVector, 3>>(result).get(0, 1) =
+      to_datavector(spatial_metric_xy);
+  get<gr::Tags::SpatialMetric<DataVector, 3>>(result).get(0, 2) =
+      to_datavector(spatial_metric_xz);
+  get<gr::Tags::SpatialMetric<DataVector, 3>>(result).get(1, 1) =
+      to_datavector(spatial_metric_yy);
+  get<gr::Tags::SpatialMetric<DataVector, 3>>(result).get(1, 2) =
+      to_datavector(spatial_metric_yz);
+  get<gr::Tags::SpatialMetric<DataVector, 3>>(result).get(2, 2) =
+      to_datavector(spatial_metric_zz);
+  get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(result).get(0, 0) =
+      to_datavector(extrinsic_curvature_xx);
+  get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(result).get(0, 1) =
+      to_datavector(extrinsic_curvature_xy);
+  get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(result).get(0, 2) =
+      to_datavector(extrinsic_curvature_xz);
+  get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(result).get(1, 1) =
+      to_datavector(extrinsic_curvature_yy);
+  get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(result).get(1, 2) =
+      to_datavector(extrinsic_curvature_yz);
+  get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(result).get(2, 2) =
+      to_datavector(extrinsic_curvature_zz);
+  get(get<hydro::Tags::RestMassDensity<DataVector>>(result)) =
+      to_datavector(rest_mass_density);
+  get<hydro::Tags::SpatialVelocity<DataVector, 3>>(result).get(0) =
+      to_datavector(spatial_velocity_x);
+  get<hydro::Tags::SpatialVelocity<DataVector, 3>>(result).get(1) =
+      to_datavector(spatial_velocity_y);
+  get<hydro::Tags::SpatialVelocity<DataVector, 3>>(result).get(2) =
+      to_datavector(spatial_velocity_z);
+  get(get<hydro::Tags::Pressure<DataVector>>(result)) = to_datavector(pressure);
+  get(get<hydro::Tags::SpecificInternalEnergy<DataVector>>(result)) =
+      to_datavector(specific_internal_energy);
+
+
+
+  return result;
+}
+
+}  // namespace io