Skip to content

Commit

Permalink
Add COCAL initial data and related modifications
Browse files Browse the repository at this point in the history
  • Loading branch information
MaiLi-UIUC committed Oct 19, 2024
1 parent 13bc1ea commit 5215be7
Show file tree
Hide file tree
Showing 10 changed files with 635 additions and 18 deletions.
22 changes: 5 additions & 17 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -40,21 +40,6 @@ endif ()
if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.27)
cmake_policy(SET CMP0144 NEW)
endif()
# - The `NEW` behavior of `DOWNLOAD_EXTRACT_TIMESTAMP` is recommended.
if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.24)
cmake_policy(SET CMP0135 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)
# - Mark FetchContent as 'SYSTEM' to suppress warnings.
# Supported since CMake version 3.25.
if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.25)
set(SPECTRE_FETCHCONTENT_BASE_ARGS SYSTEM)
else()
set(SPECTRE_FETCHCONTENT_BASE_ARGS "")
endif()

# Define standard installation directories
include(GNUInstallDirs)
Expand All @@ -71,6 +56,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`.
Expand Down Expand Up @@ -107,6 +93,7 @@ include(ProhibitInSourceBuild)
include(SpectreSetupFlagsTarget)
include(SetupFortran)
include(SetupNinjaColors)
include(SetupCocal) # Add this line
include(SetOutputDirectory)
include(SpectreAddInterfaceLibraryHeaders)
include(SpectreTargetHeaders)
Expand Down Expand Up @@ -165,8 +152,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)
Expand All @@ -185,6 +170,7 @@ include(UpdateAddExecutables)
# dependencies on the PCH if necessary.
include(SetupClangFormat)
include(SetupClangTidy)
include(SetupIwyu)
if(BUILD_DOCS)
include(SetupDoxygen)
include(SetupSphinx)
Expand All @@ -211,6 +197,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)
Expand Down
47 changes: 47 additions & 0 deletions cmake/FindCOCAL.cmake
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
# 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 COCAL_ROOT)
set(COCAL_ROOT "")
set(COCAL_ROOT $ENV{COCAL_ROOT})
endif()

find_library(
COCAL_LIB
NAMES libcocal.a
PATHS ${COCAL_ROOT}/lib
NO_DEFAULT_PATHS
)

find_path(
COCAL_INCLUDE_DIR
NAMES coc2cac_bns.f90
PATHS ${COCAL_ROOT}/include
NO_DEFAULT_PATHS
)

# 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;${OtherDependency_LIBRARIES}"
)
endif()

include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(
COCAL
FOUND_VAR COCAL_FOUND
REQUIRED_VARS COCAL_LIB COCAL_INCLUDE_DIR
)
12 changes: 12 additions & 0 deletions cmake/SetupCocal.cmake
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
# Distributed under the MIT License.
# See LICENSE.txt for details.

find_package(COCAL REQUIRED)

if (COCAL_FOUND)
message("Found cocal")
file(APPEND
"${CMAKE_BINARY_DIR}/BuildInfo.txt"
"COCAL: ${COCAL_ROOT}\n"
)
endif()
12 changes: 11 additions & 1 deletion src/Evolution/Systems/GrMhd/GhValenciaDivClean/AllSolutions.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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
35 changes: 35 additions & 0 deletions src/IO/External/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -55,3 +59,34 @@ if (TARGET FUKA::Exporter)
target_compile_definitions(
${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
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()

192 changes: 192 additions & 0 deletions src/IO/External/InterpolateFromCocal.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,192 @@
// Distributed under the MIT License.
// See LICENSE.txt for details.

#include "IO/External/InterpolateFromCocal.hpp"

#include <array>
#include <iostream>
#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" {
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,
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(const double* data, size_t size) {
// DataVector result(size);
// std::copy(data, data + size, result.begin());
// return result;
// }

// } // namespace
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 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 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::cout << "Coordinates prepared " << std::endl;

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);

// 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);

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
Loading

0 comments on commit 5215be7

Please sign in to comment.