diff --git a/CMakeLists.txt b/CMakeLists.txt index 445955dce..439514f15 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,7 @@ # ------------------------------------------------------------------------------ # CMake version requirement # ------------------------------------------------------------------------------ -cmake_minimum_required(VERSION 3.24) +cmake_minimum_required(VERSION 3.26) # ------------------------------------------------------------------------------ # Policies - use the latest of everything @@ -59,6 +59,10 @@ option(IPPL_MARK_FAILING_TESTS OFF) option(IPPL_ENABLE_SCRIPTS "Generate job script templates for some benchmarks/tests" OFF) + +option(IPPL_ENABLE_CATALYST "Enable ParaView Catalyst" OFF) +set(Catalyst_DIR "" CACHE PATH "Catalyst cmake package directory or install prefix (catalyst_DIR and CATALYST_DIR also accepted)") +set(Catalyst_VERSION "" CACHE STRING "Catalyst version or git. (default 2.1.0)") # "Build IPPL as a shared library (ON) or static library (OFF)" OFF) if(IPPL_DYL) # set(BUILD_SHARED_LIBS ON CACHE BOOL "" FORCE) message(WARNING "IPPL_DYL is deprecated; use # -DBUILD_SHARED_LIBS=ON instead.") endif() diff --git a/cmake/Dependencies.cmake b/cmake/Dependencies.cmake index fa7c770a3..e6d340d79 100644 --- a/cmake/Dependencies.cmake +++ b/cmake/Dependencies.cmake @@ -2,11 +2,12 @@ # Dependencies.cmake # ~~~ # -# Resolves third-party libraries: Kokkos and Heffte. +# Resolves third-party libraries: Kokkos, Heffte, and Catalyst. # # Responsibilities: # - Fetch or find Kokkos, using version and backends from Platforms.cmake # - Fetch Heffte if IPPL_ENABLE_FFT is ON, using CUDA or AVX2 based on platform +# - Fetch or find Catalyst when in-situ support or the FEL demo is enabled # # Not responsible for: # - Selecting platform backends → Platforms.cmake @@ -489,6 +490,114 @@ if(IPPL_ENABLE_FINUFFT) add_compile_definitions(ENABLE_FINUFFT) endif() +# ------------------------------------------------------------------------------ +# Catalyst (libcatalyst SDK and bundled Conduit parser) +# ------------------------------------------------------------------------------ +if(IPPL_ENABLE_CATALYST OR IPPL_ENABLE_FEL) + enable_language(C) + + if(NOT Catalyst_VERSION) + set(Catalyst_VERSION 2.1.0) + endif() + + extract_git_label(Catalyst_VERSION CATALYST_VERSION_GIT) + + # Support Catalyst_DIR, catalyst_DIR, CATALYST_DIR (package dir or install prefix) + foreach(_cat_dir_var IN ITEMS Catalyst_DIR catalyst_DIR CATALYST_DIR) + if(DEFINED ${_cat_dir_var} AND ${_cat_dir_var}) + if(EXISTS "${${_cat_dir_var}}/catalyst-config.cmake") + set(catalyst_DIR "${${_cat_dir_var}}") + else() + list(APPEND CMAKE_PREFIX_PATH "${${_cat_dir_var}}") + endif() + endif() + endforeach() + unset(_cat_dir_var) + + # Do not mistake the package generated by an earlier FetchContent configure + # for a user-provided Catalyst installation on a subsequent configure. + set(_ippl_reuse_fetched_catalyst OFF) + set(_ippl_catalyst_dir_is_in_build_tree OFF) + if(catalyst_DIR) + cmake_path(IS_PREFIX CMAKE_BINARY_DIR "${catalyst_DIR}" NORMALIZE + _ippl_catalyst_dir_is_in_build_tree) + endif() + if(NOT Catalyst_DIR AND NOT CATALYST_DIR + AND (_ippl_catalyst_dir_is_in_build_tree + OR (IPPL_CATALYST_FETCHED AND NOT catalyst_DIR))) + set(_ippl_reuse_fetched_catalyst ON) + unset(catalyst_DIR CACHE) + unset(catalyst_FOUND) + endif() + + if(NOT CATALYST_VERSION_GIT AND NOT _ippl_reuse_fetched_catalyst) + find_package(catalyst ${Catalyst_VERSION} CONFIG QUIET) + set(CATALYST_VERSION_GIT "v${Catalyst_VERSION}") + endif() + + if(NOT CATALYST_VERSION_GIT) + set(CATALYST_VERSION_GIT "v${Catalyst_VERSION}") + endif() + + if(catalyst_FOUND AND NOT _ippl_reuse_fetched_catalyst) + set(IPPL_CATALYST_FETCHED OFF CACHE INTERNAL "Catalyst was fetched by IPPL" FORCE) + set(IPPL_CATALYST_VERSION "${catalyst_VERSION}") + colour_message(STATUS ${Green} "✅ Catalyst ${catalyst_VERSION} found externally") + else() + set(IPPL_CATALYST_FETCHED ON CACHE INTERNAL "Catalyst was fetched by IPPL" FORCE) + colour_message(STATUS ${Green} "✅ Catalyst ${CATALYST_VERSION_GIT} building from source") + + set(CATALYST_BUILD_SHARED_LIBS ON CACHE BOOL "" FORCE) + set(CATALYST_BUILD_STUB_IMPLEMENTATION ON CACHE BOOL "" FORCE) + set(CATALYST_BUILD_TESTING OFF CACHE BOOL "" FORCE) + set(CATALYST_BUILD_TOOLS OFF CACHE BOOL "" FORCE) + set(CATALYST_WRAP_PYTHON OFF CACHE BOOL "" FORCE) + set(CATALYST_WRAP_FORTRAN OFF CACHE BOOL "" FORCE) + set(CATALYST_WITH_EXTERNAL_CONDUIT OFF CACHE BOOL "" FORCE) + set(CATALYST_USE_MPI ON CACHE BOOL "" FORCE) + + FetchContent_Declare( + catalyst + GIT_REPOSITORY "https://gitlab.kitware.com/paraview/catalyst.git" + GIT_TAG "${CATALYST_VERSION_GIT}" + DOWNLOAD_EXTRACT_TIMESTAMP ON) + FetchContent_MakeAvailable(catalyst) + + if(NOT TARGET catalyst::catalyst) + message(FATAL_ERROR "Catalyst FetchContent did not provide catalyst::catalyst") + endif() + + file( + STRINGS "${catalyst_SOURCE_DIR}/CMakeLists.txt" _catalyst_project_line + REGEX "^[ \t]*project\\(CATALYST VERSION [0-9]+\\.[0-9]+") + string(REGEX MATCH "VERSION[ \t]+([0-9]+\\.[0-9]+(\\.[0-9]+)?)" + _catalyst_version_match "${_catalyst_project_line}") + if(NOT CMAKE_MATCH_1) + message(FATAL_ERROR "Could not determine the fetched Catalyst project version") + endif() + set(IPPL_CATALYST_VERSION "${CMAKE_MATCH_1}") + + unset(_catalyst_project_line) + unset(_catalyst_version_match) + endif() + + string(REGEX MATCH "^([0-9]+)\\.([0-9]+)" _catalyst_major_minor + "${IPPL_CATALYST_VERSION}") + if(NOT _catalyst_major_minor) + message(FATAL_ERROR "Could not determine Catalyst's major/minor version") + endif() + set(IPPL_CATALYST_PACKAGE_DIRNAME "catalyst-${CMAKE_MATCH_1}.${CMAKE_MATCH_2}") + if(IPPL_CATALYST_FETCHED) + set(IPPL_CATALYST_BUILD_PACKAGE_DIR + "${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_LIBDIR}/cmake/${IPPL_CATALYST_PACKAGE_DIRNAME}") + endif() + + set(IPPL_CATALYST_USE_MPI "${CATALYST_USE_MPI}") + unset(_catalyst_major_minor) + unset(_ippl_catalyst_dir_is_in_build_tree) + unset(_ippl_reuse_fetched_catalyst) +endif() + # ------------------------------------------------------------------------------ # GoogleTest # ------------------------------------------------------------------------------ @@ -510,22 +619,3 @@ if(IPPL_ENABLE_UNIT_TESTS) message(STATUS "✅ GoogleTest built from source (${GTest_VERSION})") endif() endif() - -# ------------------------------------------------------------------------------ -# FEL module header-only dependencies (nlohmann/json for config parsing, stb_image_write for the -# Poynting-flux visualization). -# ------------------------------------------------------------------------------ -if(IPPL_ENABLE_FEL) - # Fetch the CMake package instead of downloading the release header directly. CMake's - # file(DOWNLOAD) does not fail by default and can leave a zero-byte json.hpp behind when a - # release-asset host is unavailable, which only surfaces later as a confusing compile error. - set(JSON_BuildTests OFF CACHE BOOL "Disable nlohmann/json tests" FORCE) - FetchContent_Declare( - nlohmann_json - GIT_REPOSITORY https://github.com/nlohmann/json.git - GIT_TAG v3.11.3 - GIT_SHALLOW ON) - FetchContent_MakeAvailable(nlohmann_json) - - message(STATUS "✅ nlohmann/json loaded for the FEL module.") -endif() diff --git a/cmake/IPPLConfig.cmake.in b/cmake/IPPLConfig.cmake.in index ab10e60e1..40f3f9fa8 100644 --- a/cmake/IPPLConfig.cmake.in +++ b/cmake/IPPLConfig.cmake.in @@ -44,6 +44,45 @@ if(@IPPL_ENABLE_FFT@ AND "@Heffte_FOUND@" AND NOT TARGET Heffte::Heffte) find_dependency(Heffte CONFIG HINTS "${_ipplPrefix}" "@Heffte_DIR@") endif() +if(@IPPL_ENABLE_CATALYST@) + if("@IPPL_CATALYST_USE_MPI@" AND NOT CMAKE_C_COMPILER_LOADED) + set(${CMAKE_FIND_PACKAGE_NAME}_FOUND FALSE) + set(${CMAKE_FIND_PACKAGE_NAME}_NOT_FOUND_MESSAGE + "IPPL was built with MPI-enabled Catalyst, which requires the C language. Enable it in the consuming project with project(... LANGUAGES C CXX).") + return() + endif() + + if(NOT TARGET catalyst::catalyst) + message(STATUS "Finding IPPL dependency: Catalyst") + + if(_ipplBuildTree) + if("@IPPL_CATALYST_FETCHED@") + find_dependency(catalyst @IPPL_CATALYST_VERSION@ CONFIG + HINTS "@IPPL_CATALYST_BUILD_PACKAGE_DIR@") + else() + find_dependency(catalyst @IPPL_CATALYST_VERSION@ CONFIG HINTS "@catalyst_DIR@") + endif() + elseif("@IPPL_CATALYST_FETCHED@") + find_dependency( + catalyst @IPPL_CATALYST_VERSION@ CONFIG + HINTS + "${_ipplPrefix}/@CMAKE_INSTALL_LIBDIR@/cmake/@IPPL_CATALYST_PACKAGE_DIRNAME@") + else() + find_dependency(catalyst @IPPL_CATALYST_VERSION@ CONFIG + HINTS "${_ipplPrefix}" "@catalyst_DIR@") + endif() + + message(STATUS "Catalyst version ${catalyst_VERSION}") + endif() +endif() + include("${CMAKE_CURRENT_LIST_DIR}/IPPLTargets.cmake") + +if("@IPPL_ENABLE_CATALYST@" AND NOT _ipplBuildTree) + set_property( + TARGET IPPL::ippl APPEND PROPERTY INTERFACE_COMPILE_DEFINITIONS + "IPPL_CATALYST_SCRIPTS_DIR=\"${_ipplPrefix}/@CMAKE_INSTALL_DATADIR@/ippl/catalyst_scripts\"") +endif() + unset(_ipplPrefix) unset(_ipplBuildTree) diff --git a/cmake/InstallIppl.cmake b/cmake/InstallIppl.cmake index b3ba241e2..7e95d22be 100644 --- a/cmake/InstallIppl.cmake +++ b/cmake/InstallIppl.cmake @@ -29,6 +29,7 @@ set(_ippl_install_excludes PATTERN "*.cc" EXCLUDE PATTERN "*.cpp" EXCLUDE PATTERN "*.cu" EXCLUDE + PATTERN "catalyst_scripts" EXCLUDE # match your filenames if needed ) diff --git a/cmake/SetupCatalyst.cmake b/cmake/SetupCatalyst.cmake new file mode 100644 index 000000000..ebd88172e --- /dev/null +++ b/cmake/SetupCatalyst.cmake @@ -0,0 +1,43 @@ +# ----------------------------------------------------------------------------- +# SetupCatalyst.cmake +# +# Links Catalyst into the ippl target. Resolution (find vs FetchContent) happens +# in Dependencies.cmake. +# ----------------------------------------------------------------------------- + +if(NOT TARGET ippl) + message(FATAL_ERROR "SetupCatalyst.cmake must be included after the ippl target is created.") +endif() + +if(NOT TARGET catalyst::catalyst) + message( + FATAL_ERROR + "Catalyst enabled but catalyst::catalyst target is missing. " + "Set -DCatalyst_DIR=/path/to/catalyst/cmake or allow FetchContent.") +endif() + +message(STATUS "Catalyst enabled") + +target_compile_definitions( + ippl + PUBLIC IPPL_ENABLE_CATALYST + $) + +target_link_libraries(ippl PUBLIC catalyst::catalyst) + +message(STATUS "Catalyst Summary:") +message(STATUS " External install: ${catalyst_FOUND}") +if(catalyst_VERSION) + message(STATUS " Version: ${catalyst_VERSION}") +elseif(Catalyst_VERSION) + message(STATUS " Version: ${Catalyst_VERSION}") +endif() +if(catalyst_DIR) + message(STATUS " CMake config dir: ${catalyst_DIR}") +elseif(Catalyst_PACKAGE_DIR) + message(STATUS " CMake config dir: ${Catalyst_PACKAGE_DIR}") +endif() + +get_target_property(_cat_type catalyst::catalyst TYPE) +message(STATUS " Target: catalyst::catalyst (${_cat_type})") +unset(_cat_type) diff --git a/demos/alpine/ExamplesWithoutPicManager/ChargedParticles.hpp b/demos/alpine/ExamplesWithoutPicManager/ChargedParticles.hpp index 47d0a4124..af444d3a0 100644 --- a/demos/alpine/ExamplesWithoutPicManager/ChargedParticles.hpp +++ b/demos/alpine/ExamplesWithoutPicManager/ChargedParticles.hpp @@ -262,10 +262,13 @@ class ChargedParticles : public ippl::ParticleBase { } void registerAttributes() { + P.set_name("velocity"); + E.set_name("electric_field"); + q.set_name("charge"); // register the particle attributes - this->addAttribute(q); this->addAttribute(P); this->addAttribute(E); + this->addAttribute(q); } ~ChargedParticles() {} diff --git a/demos/alpine/ParticleContainer.hpp b/demos/alpine/ParticleContainer.hpp index 8af0a0acf..8b183307a 100644 --- a/demos/alpine/ParticleContainer.hpp +++ b/demos/alpine/ParticleContainer.hpp @@ -38,9 +38,9 @@ class ParticleContainer : public ippl::ParticleBase< q.set_name("charge"); E.set_name("electric_field"); // register the particle attributes - this->addAttribute(q); this->addAttribute(P); this->addAttribute(E); + this->addAttribute(q); } void setupBCs() { setBCAllPeriodic(); } diff --git a/demos/alpine/PenningTrap.cpp b/demos/alpine/PenningTrap.cpp index 06cf43153..da95fd45c 100644 --- a/demos/alpine/PenningTrap.cpp +++ b/demos/alpine/PenningTrap.cpp @@ -92,6 +92,10 @@ int main(int argc, char* argv[]) { manager.run(manager.getNt()); + #ifdef IPPL_ENABLE_CATALYST + manager.cat_viz.Finalize(); + #endif + msg << "End." << endl; IpplTimings::stopTimer(mainTimer); diff --git a/demos/alpine/PenningTrapManager.h b/demos/alpine/PenningTrapManager.h index c7e81de43..1c879ed27 100644 --- a/demos/alpine/PenningTrapManager.h +++ b/demos/alpine/PenningTrapManager.h @@ -15,6 +15,10 @@ #include "Random/NormalDistribution.h" #include "Random/Randn.h" +#ifdef IPPL_ENABLE_CATALYST +#include "Stream/InSitu/CatalystAdaptor.h" +#endif + using view_type = typename ippl::detail::ViewType, 1>::view_type; template @@ -45,6 +49,13 @@ class PenningTrapManager : public AlpineManager { double alpha_m; double DrInv_m; + + + #ifdef IPPL_ENABLE_CATALYST + public: + ippl::CatalystAdaptor cat_viz{std::string{TestName}}; + #endif + public: void pre_run() override { Inform m("Pre Run"); @@ -124,6 +135,24 @@ class PenningTrapManager : public AlpineManager { this->grid2par(); + #ifdef IPPL_ENABLE_CATALYST + m << "Catalyst is enabled" << endl; + + std::shared_ptr runtime_steer_registry = ippl::MakeVisRegistryRuntimePtr(); + std::shared_ptr runtime_vis_registry = ippl::MakeVisRegistryRuntimePtr( + "density", this->fcontainer_m->getRho(), + "ions", this->pcontainer_m + ); + // runtime_vis_registry->add("potential", this->fcontainer_m->getRho() ); + runtime_vis_registry->add("electrostatic", this->fcontainer_m->getE() ); + + static IpplTimings::TimerRef CAinit = IpplTimings::getTimer("CAinit"); + IpplTimings::startTimer(CAinit); + cat_viz.Initialize(runtime_vis_registry, runtime_steer_registry); + IpplTimings::stopTimer(CAinit); + + #endif + this->dump(); m << "Done"; @@ -239,6 +268,11 @@ class PenningTrapManager : public AlpineManager { static IpplTimings::TimerRef domainDecomposition = IpplTimings::getTimer("loadBalance"); static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("solve"); + #ifdef IPPL_ENABLE_CATALYST + static IpplTimings::TimerRef TMR_CAremember = IpplTimings::getTimer("CAremember"); + static IpplTimings::TimerRef TMR_CAexecute = IpplTimings::getTimer("CAexecute"); + #endif + double alpha = this->alpha_m; double Bext = this->Bext_m; double DrInv = this->DrInv_m; @@ -298,6 +332,13 @@ class PenningTrapManager : public AlpineManager { // scatter the charge onto the underlying grid this->par2grid(); + // Save deep copy of the density field (for later). + #ifdef IPPL_ENABLE_CATALYST + IpplTimings::startTimer(TMR_CAremember); + cat_viz.rememberNow("density"); + IpplTimings::stopTimer(TMR_CAremember); + #endif + // Field solve IpplTimings::startTimer(SolveTimer); this->fsolver_m->runSolver(); @@ -306,6 +347,13 @@ class PenningTrapManager : public AlpineManager { // gather E field this->grid2par(); + //trigger In Situ pipeline + #ifdef IPPL_ENABLE_CATALYST + IpplTimings::startTimer(TMR_CAexecute); + cat_viz.Execute(it, this->time_m); + IpplTimings::stopTimer(TMR_CAexecute); + #endif + IpplTimings::startTimer(PTimer); auto R2view = pc->R.getView(); auto P2view = pc->P.getView(); diff --git a/demos/collisions/P3MParticleContainer.hpp b/demos/collisions/P3MParticleContainer.hpp index 9d964e3e9..07b3cac8f 100644 --- a/demos/collisions/P3MParticleContainer.hpp +++ b/demos/collisions/P3MParticleContainer.hpp @@ -51,10 +51,13 @@ class P3MParticleContainer : public ippl::ParticleBasesetParticleBC(ippl::BC::PERIODIC); } void registerAttributes() { + P.set_name("velocity"); + E.set_name("electric_field"); + Q.set_name("charge"); // register the particle attributes - this->addAttribute(Q); this->addAttribute(P); this->addAttribute(E); + this->addAttribute(Q); } }; diff --git a/demos/cosmology/GravityParticleContainer.hpp b/demos/cosmology/GravityParticleContainer.hpp index e670da577..59be4e567 100644 --- a/demos/cosmology/GravityParticleContainer.hpp +++ b/demos/cosmology/GravityParticleContainer.hpp @@ -69,9 +69,12 @@ class ParticleContainer : public ippl::ParticleBase< * @brief Register the particle attributes. */ void registerAttributes() { - this->addAttribute(m); + V.set_name("velocity"); + F.set_name("gravitational_field"); + m.set_name("mass"); this->addAttribute(V); this->addAttribute(F); + this->addAttribute(m); } /** @@ -91,4 +94,4 @@ class ParticleContainer : public ippl::ParticleBase< void setBCAllPeriodic() { this->setParticleBC(ippl::BC::PERIODIC); } }; -#endif // IPPL_PARTICLE_CONTAINER_H \ No newline at end of file +#endif // IPPL_PARTICLE_CONTAINER_H diff --git a/demos/electrostaticPIF/ChargedParticlesPIF.hpp b/demos/electrostaticPIF/ChargedParticlesPIF.hpp index caacaf27f..9d8014286 100644 --- a/demos/electrostaticPIF/ChargedParticlesPIF.hpp +++ b/demos/electrostaticPIF/ChargedParticlesPIF.hpp @@ -157,10 +157,13 @@ class ChargedParticlesPIF : public ippl::ParticleBase { */ ChargedParticlesPIF(PLayout& pl) : ippl::ParticleBase(pl) { + P.set_name("velocity"); + E.set_name("electric_field"); + q.set_name("charge"); // register the particle attributes - this->addAttribute(q); this->addAttribute(P); this->addAttribute(E); + this->addAttribute(q); } ChargedParticlesPIF(PLayout& pl, Vector_t hr, Vector_t rmin, Vector_t rmax, @@ -174,10 +177,13 @@ class ChargedParticlesPIF : public ippl::ParticleBase { , Np_m(Np) , useUpsampledInputs_m(useUpsampledInputs) , useFinufft_m(useFinufft) { + P.set_name("velocity"); + E.set_name("electric_field"); + q.set_name("charge"); // register the particle attributes - this->addAttribute(q); this->addAttribute(P); this->addAttribute(E); + this->addAttribute(q); setupBCs(); for (unsigned int i = 0; i < Dim; i++) decomp_m[i] = decomp[i]; diff --git a/demos/fel/CMakeLists.txt b/demos/fel/CMakeLists.txt index 7e4bc3fb1..39a8b87d7 100644 --- a/demos/fel/CMakeLists.txt +++ b/demos/fel/CMakeLists.txt @@ -7,9 +7,23 @@ message(STATUS "Configuring demos/fel/") add_executable(FreeElectronLaser FreeElectronLaser.cpp) +set(IPPL_FEL_CONFIG_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/config.json") +set(IPPL_FEL_CONFIG_BUILD "${CMAKE_CURRENT_BINARY_DIR}/config.json") + +add_custom_command( + OUTPUT "${IPPL_FEL_CONFIG_BUILD}" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different "${IPPL_FEL_CONFIG_SOURCE}" + "${IPPL_FEL_CONFIG_BUILD}" + DEPENDS "${IPPL_FEL_CONFIG_SOURCE}" + COMMENT "Staging the FEL configuration" + VERBATIM) + +add_custom_target(ippl_fel_config DEPENDS "${IPPL_FEL_CONFIG_BUILD}") +add_dependencies(FreeElectronLaser ippl_fel_config) + target_link_libraries( FreeElectronLaser - PRIVATE IPPL::ippl nlohmann_json::nlohmann_json) + PRIVATE IPPL::ippl catalyst::catalyst) target_include_directories(FreeElectronLaser PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) target_include_directories(FreeElectronLaser SYSTEM PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/stb) diff --git a/demos/fel/Config.h b/demos/fel/Config.h index 5efa08614..87136a1f4 100644 --- a/demos/fel/Config.h +++ b/demos/fel/Config.h @@ -8,19 +8,21 @@ // FreeElectronLaser.cpp. #include +#include +#include +#include #include -#include #include +#include +#include #include +#include #include #include "Types/Vector.h" #include "units.h" -#define JSON_HAS_RANGES 0 -#include - struct config { using scalar = double; @@ -40,36 +42,133 @@ struct config { bool space_charge; // Flag for considering space charge effects // BUNCH PARAMETERS - ippl::Vector mean_position; // Mean initial position of the particle bunch - ippl::Vector sigma_position; // Standard deviation of the initial position distribution + ippl::Vector mean_position; // Mean initial position of the particle bunch + ippl::Vector + sigma_position; // Standard deviation of the initial position distribution ippl::Vector position_truncations; // Truncations of the position distribution - ippl::Vector sigma_momentum; // Standard deviation of the initial momentum distribution - scalar bunch_gamma; // Relativistic gamma factor of the bunch + ippl::Vector + sigma_momentum; // Standard deviation of the initial momentum distribution + scalar bunch_gamma; // Relativistic gamma factor of the bunch // UNDULATOR PARAMETERS scalar undulator_K; // Undulator parameter K scalar undulator_period; // Period of the undulator scalar undulator_length; // Length of the undulator - uint32_t output_rhythm; // Frequency of output in timesteps std::string output_path; // Path to output files std::unordered_map experiment_options; // Additional experimental options }; -template -ippl::Vector getVector(const nlohmann::json& j) { - if (j.is_array()) { - assert(j.size() == Dim); - ippl::Vector ret; - for (unsigned i = 0; i < Dim; i++) - ret[i] = (scalar)j[i]; - return ret; - } else { - std::cerr << "Warning: Obtaining Vector from scalar json\n"; - ippl::Vector ret = (scalar)j; - return ret; +namespace fel_config_detail { + + inline conduit_cpp::Node requiredNode(const conduit_cpp::Node& root, const std::string& path) { + if (!root.has_path(path)) { + throw std::runtime_error("Missing required configuration value '" + path + "'"); + } + return root[path]; } -} + + inline long double numericElement(const conduit_cpp::Node& node, conduit_index_t index, + const std::string& path) { + if (!node.dtype().is_number() || index < 0 || index >= node.number_of_elements()) { + throw std::runtime_error("Configuration value '" + path + "' must be numeric"); + } + + using Id = conduit_cpp::DataType::Id; + switch (node.dtype().id()) { + case Id::int8: + return node.as_int8_ptr()[index]; + case Id::int16: + return node.as_int16_ptr()[index]; + case Id::int32: + return node.as_int32_ptr()[index]; + case Id::int64: + return static_cast(node.as_int64_ptr()[index]); + case Id::uint8: + return node.as_uint8_ptr()[index]; + case Id::uint16: + return node.as_uint16_ptr()[index]; + case Id::uint32: + return node.as_uint32_ptr()[index]; + case Id::uint64: + return static_cast(node.as_uint64_ptr()[index]); + case Id::float32: + return node.as_float32_ptr()[index]; + case Id::float64: + return node.as_float64_ptr()[index]; + case Id::unknown: + break; + } + + throw std::runtime_error("Unsupported numeric type for configuration value '" + path + "'"); + } + + inline double requiredNumber(const conduit_cpp::Node& root, const std::string& path) { + const auto node = requiredNode(root, path); + if (node.number_of_elements() != 1) { + throw std::runtime_error("Configuration value '" + path + "' must be a scalar"); + } + return static_cast(numericElement(node, 0, path)); + } + + inline std::string requiredString(const conduit_cpp::Node& root, const std::string& path) { + const auto node = requiredNode(root, path); + if (!node.dtype().is_string()) { + throw std::runtime_error("Configuration value '" + path + "' must be a string"); + } + return node.as_string(); + } + + template + Scalar checkedNumericCast(long double value, const std::string& path) { + if constexpr (std::is_integral_v) { + if (!std::isfinite(value) || std::trunc(value) != value + || value < static_cast(std::numeric_limits::lowest()) + || value > static_cast(std::numeric_limits::max())) { + throw std::runtime_error("Configuration value '" + path + + "' must be an in-range integer"); + } + } + return static_cast(value); + } + + template + Scalar requiredInteger(const conduit_cpp::Node& root, const std::string& path) { + static_assert(std::is_integral_v); + const auto node = requiredNode(root, path); + if (node.number_of_elements() != 1) { + throw std::runtime_error("Configuration value '" + path + "' must be a scalar"); + } + return checkedNumericCast(numericElement(node, 0, path), path); + } + + template + ippl::Vector getVector(const conduit_cpp::Node& root, const std::string& path) { + const auto node = requiredNode(root, path); + if (!node.dtype().is_number()) { + throw std::runtime_error("Configuration value '" + path + "' must be numeric"); + } + + const auto size = node.number_of_elements(); + ippl::Vector result; + if (size == 1) { + std::cerr << "Warning: Obtaining vector from scalar configuration value '" << path + << "'\n"; + result = checkedNumericCast(numericElement(node, 0, path), path); + return result; + } + if (size != static_cast(Dim)) { + throw std::runtime_error("Configuration value '" + path + "' must contain " + + std::to_string(Dim) + " elements"); + } + + for (unsigned i = 0; i < Dim; ++i) { + result[i] = checkedNumericCast(numericElement(node, i, path), path); + } + return result; + } + +} // namespace fel_config_detail // Compile-time / run-time string hashing used to switch on unit-scale names. template @@ -127,7 +226,9 @@ inline size_t chash(const std::string& _val) { } inline std::string lowercase_singular(std::string str) { // Convert string to lowercase - std::transform(str.begin(), str.end(), str.begin(), ::tolower); + std::transform(str.begin(), str.end(), str.begin(), [](unsigned char c) { + return static_cast(std::tolower(c)); + }); // Check if the string ends with "s" and remove it if it does if (!str.empty() && str.back() == 's') { @@ -136,10 +237,11 @@ inline std::string lowercase_singular(std::string str) { return str; } -inline double get_time_multiplier(const nlohmann::json& j) { - std::string length_scale_string = lowercase_singular((std::string)j["mesh"]["time-scale"]); - double time_factor = 1.0; - switch (chash(length_scale_string)) { +inline double get_time_multiplier(const conduit_cpp::Node& root) { + const std::string time_scale = fel_config_detail::requiredString(root, "mesh/time-scale"); + std::string time_scale_string = lowercase_singular(time_scale); + double time_factor = 1.0; + switch (chash(time_scale_string)) { case chash<"planck-time">(): case chash<"plancktime">(): case chash<"pt">(): @@ -162,16 +264,15 @@ inline double get_time_multiplier(const nlohmann::json& j) { time_factor = 1.0; break; default: - std::cerr << "Unrecognized time scale: " << (std::string)j["mesh"]["time-scale"] - << "\n"; + std::cerr << "Unrecognized time scale: " << time_scale << "\n"; break; } return time_factor; } -inline double get_length_multiplier(const nlohmann::json& options) { - std::string length_scale_string = - lowercase_singular((std::string)options["mesh"]["length-scale"]); - double length_factor = 1.0; +inline double get_length_multiplier(const conduit_cpp::Node& root) { + const std::string length_scale = fel_config_detail::requiredString(root, "mesh/length-scale"); + std::string length_scale_string = lowercase_singular(length_scale); + double length_factor = 1.0; switch (chash(length_scale_string)) { case chash<"planck-length">(): case chash<"plancklength">(): @@ -195,78 +296,98 @@ inline double get_length_multiplier(const nlohmann::json& options) { length_factor = 1.0; break; default: - std::cerr << "Unrecognized length scale: " - << (std::string)options["mesh"]["length-scale"] << "\n"; + std::cerr << "Unrecognized length scale: " << length_scale << "\n"; break; } return length_factor; } inline config read_config(const char* filepath) { - std::ifstream cfile(filepath); - nlohmann::json j; - cfile >> j; - config::scalar lmult = get_length_multiplier(j); - config::scalar tmult = get_time_multiplier(j); - config ret; - - ret.extents[0] = ((config::scalar)j["mesh"]["extents"][0] * lmult) / unit_length_in_meters; - ret.extents[1] = ((config::scalar)j["mesh"]["extents"][1] * lmult) / unit_length_in_meters; - ret.extents[2] = ((config::scalar)j["mesh"]["extents"][2] * lmult) / unit_length_in_meters; - ret.resolution = getVector(j["mesh"]["resolution"]); - - if (j.contains("timestep-ratio")) { - ret.timestep_ratio = (config::scalar)j["timestep-ratio"]; - } else { - ret.timestep_ratio = 1; - } - ret.total_time = ((config::scalar)j["mesh"]["total-time"] * tmult) / unit_time_in_seconds; - ret.space_charge = (bool)(j["mesh"]["space-charge"]); - ret.bunch_gamma = (config::scalar)(j["bunch"]["gamma"]); - if (ret.bunch_gamma < config::scalar(1)) { - std::cerr << "Gamma must be >= 1\n"; - exit(1); - } - assert(j.contains("undulator")); - assert(j["undulator"].contains("static-undulator")); - - ret.undulator_K = j["undulator"]["static-undulator"]["undulator-parameter"]; - ret.undulator_period = ((config::scalar)j["undulator"]["static-undulator"]["period"] * lmult) - / unit_length_in_meters; - ret.undulator_length = ((config::scalar)j["undulator"]["static-undulator"]["length"] * lmult) - / unit_length_in_meters; - assert(!std::isnan(ret.undulator_length)); - assert(!std::isnan(ret.undulator_period)); - assert(!std::isnan(ret.extents[0])); - assert(!std::isnan(ret.extents[1])); - assert(!std::isnan(ret.extents[2])); - assert(!std::isnan(ret.total_time)); - ret.length_scale_in_jobfile = get_length_multiplier(j); - ret.temporal_scale_in_jobfile = get_time_multiplier(j); - ret.charge = (config::scalar)j["bunch"]["charge"] * electron_charge_in_unit_charges; - ret.mass = (config::scalar)j["bunch"]["mass"] * electron_mass_in_unit_masses; - ret.num_particles = (uint64_t)j["bunch"]["number-of-particles"]; - ret.mean_position = - getVector(j["bunch"]["position"]) * lmult / unit_length_in_meters; - ret.sigma_position = - getVector(j["bunch"]["sigma-position"]) * lmult / unit_length_in_meters; - ret.position_truncations = getVector(j["bunch"]["distribution-truncations"]) - * lmult / unit_length_in_meters; - ret.sigma_momentum = getVector(j["bunch"]["sigma-momentum"]); - ret.output_rhythm = j["output"].contains("rhythm") ? uint32_t(j["output"]["rhythm"]) : 0; - ret.output_path = "../data/"; - if (j["output"].contains("path")) { - ret.output_path = j["output"]["path"]; - if (!ret.output_path.ends_with('/')) { - ret.output_path.push_back('/'); + try { + conduit_cpp::Node root; + conduit_node_load(conduit_cpp::c_node(&root), filepath, "json"); + + const config::scalar lmult = get_length_multiplier(root); + const config::scalar tmult = get_time_multiplier(root); + config ret{}; + + ret.extents = fel_config_detail::getVector(root, "mesh/extents") * lmult + / unit_length_in_meters; + ret.resolution = fel_config_detail::getVector(root, "mesh/resolution"); + + ret.timestep_ratio = root.has_path("timestep-ratio") + ? fel_config_detail::requiredNumber(root, "timestep-ratio") + : config::scalar(1); + ret.total_time = fel_config_detail::requiredNumber(root, "mesh/total-time") * tmult + / unit_time_in_seconds; + ret.space_charge = + fel_config_detail::requiredNumber(root, "mesh/space-charge") != config::scalar(0); + ret.bunch_gamma = fel_config_detail::requiredNumber(root, "bunch/gamma"); + if (ret.bunch_gamma < config::scalar(1)) { + throw std::runtime_error("Configuration value 'bunch/gamma' must be >= 1"); } - } - if (j.contains("experimentation")) { - nlohmann::json je = j["experimentation"]; - for (auto it = je.begin(); it != je.end(); it++) { - ret.experiment_options[it.key()] = double(it.value()); + + ret.undulator_K = fel_config_detail::requiredNumber( + root, "undulator/static-undulator/undulator-parameter"); + ret.undulator_period = + fel_config_detail::requiredNumber(root, "undulator/static-undulator/period") * lmult + / unit_length_in_meters; + ret.undulator_length = + fel_config_detail::requiredNumber(root, "undulator/static-undulator/length") * lmult + / unit_length_in_meters; + + if (!std::isfinite(ret.undulator_length) || !std::isfinite(ret.undulator_period) + || !std::isfinite(ret.extents[0]) || !std::isfinite(ret.extents[1]) + || !std::isfinite(ret.extents[2]) || !std::isfinite(ret.total_time)) { + throw std::runtime_error("FEL configuration contains a non-finite physical value"); + } + + ret.length_scale_in_jobfile = lmult; + ret.temporal_scale_in_jobfile = tmult; + ret.charge = fel_config_detail::requiredNumber(root, "bunch/charge") + * electron_charge_in_unit_charges; + ret.mass = + fel_config_detail::requiredNumber(root, "bunch/mass") * electron_mass_in_unit_masses; + ret.num_particles = fel_config_detail::requiredInteger( + root, "bunch/number-of-particles"); + ret.mean_position = fel_config_detail::getVector(root, "bunch/position") + * lmult / unit_length_in_meters; + ret.sigma_position = + fel_config_detail::getVector(root, "bunch/sigma-position") * lmult + / unit_length_in_meters; + ret.position_truncations = + fel_config_detail::getVector(root, "bunch/distribution-truncations") + * lmult / unit_length_in_meters; + ret.sigma_momentum = + fel_config_detail::getVector(root, "bunch/sigma-momentum"); + + ret.output_path = "../data/"; + if (root.has_path("output/path")) { + ret.output_path = fel_config_detail::requiredString(root, "output/path"); + if (!ret.output_path.ends_with('/')) { + ret.output_path.push_back('/'); + } + } + + if (root.has_path("experimentation")) { + const auto experimentation = root["experimentation"]; + if (!experimentation.dtype().is_object()) { + throw std::runtime_error("Configuration value 'experimentation' must be an object"); + } + for (conduit_index_t i = 0; i < experimentation.number_of_children(); ++i) { + const auto option = experimentation.child(i); + if (!option.dtype().is_number() || option.number_of_elements() != 1) { + throw std::runtime_error("Experimentation option '" + option.name() + + "' must be a numeric scalar"); + } + ret.experiment_options[option.name()] = option.to_double(); + } } + + return ret; + } catch (const std::exception& error) { + throw std::runtime_error("Failed to read FEL configuration '" + std::string(filepath) + + "': " + error.what()); } - return ret; } #endif diff --git a/demos/fel/FELParticleContainer.hpp b/demos/fel/FELParticleContainer.hpp index 802bc9c4b..3799fea88 100644 --- a/demos/fel/FELParticleContainer.hpp +++ b/demos/fel/FELParticleContainer.hpp @@ -50,13 +50,21 @@ class FELParticleContainer : public ippl::ParticleBase> { void setPL(std::shared_ptr>& pl) { pl_m = pl; } void registerAttributes() { + gamma_beta.set_name("gamma_beta"); + E_gather.set_name("electric_field"); + B_gather.set_name("magnetic_field"); + Q.set_name("charge"); + mass.set_name("mass"); + R_nm1.set_name("previous_position"); + R_np1.set_name("next_position"); + + this->addAttribute(gamma_beta); + this->addAttribute(E_gather); + this->addAttribute(B_gather); this->addAttribute(Q); this->addAttribute(mass); - this->addAttribute(gamma_beta); this->addAttribute(R_nm1); this->addAttribute(R_np1); - this->addAttribute(E_gather); - this->addAttribute(B_gather); } void setupBCs() { setBCAllOpen(); } diff --git a/demos/fel/FreeElectronLaser.cpp b/demos/fel/FreeElectronLaser.cpp index 22e3d4406..2997ee698 100644 --- a/demos/fel/FreeElectronLaser.cpp +++ b/demos/fel/FreeElectronLaser.cpp @@ -5,14 +5,20 @@ // or // mpirun -np [N] ./FreeElectronLaser [] --info [0-5] // -// Reads a MITHRA-style JSON job file (default: ../fel/config.json) describing the -// grid, the relativistic electron bunch, and the undulator. The simulation +// Reads a MITHRA-style JSON job file (by default the config.json staged next to +// the build-tree executable) describing the grid, the relativistic electron +// bunch, and the undulator. The simulation // runs in a Lorentz frame co-moving with the bunch: a charge-conserving // current is deposited onto the grid, Maxwell's equations are advanced with a // standard FDTD solver (absorbing boundaries), and the particles are pushed // with a relativistic Boris pusher that also feels the (frame-transformed) // undulator field. Radiated power is written to a CSV and a -// Poynting-flux video is produced via ffmpeg. +// narrow-band radiation diagnostic is produced alongside it. + + + + // "resolution": [96, 96, 3000], + // "resolution": [48, 48, 1500], constexpr unsigned Dim = 3; using T = double; @@ -24,11 +30,12 @@ using T = double; #include "Utility/IpplTimings.h" -// stb_image_write's implementation must be emitted in exactly one translation -// unit; define the macro before the (transitive) include of the header. -#define STB_IMAGE_WRITE_IMPLEMENTATION #include "FreeElectronLaserManager.h" +#ifndef IPPL_FEL_DEFAULT_CONFIG +#define IPPL_FEL_DEFAULT_CONFIG "config.json" +#endif + int main(int argc, char* argv[]) { ippl::initialize(argc, argv); { @@ -37,16 +44,14 @@ int main(int argc, char* argv[]) { static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total"); IpplTimings::startTimer(mainTimer); - // First positional argument (if any, and not an --option) is the config - // file path; otherwise fall back to the shipped example config. The - // path is relative to the working directory, so this default assumes - // the program is launched from the build directory. - const char* config_path = "../fel/config.json"; + // First positional argument (if any, and not an --option) overrides the + // example configuration staged next to the build-tree executable. + std::string config_path = IPPL_FEL_DEFAULT_CONFIG; if (argc > 1 && argv[1][0] != '-') { config_path = argv[1]; } msg << "Reading configuration from " << config_path << endl; - config cfg = read_config(config_path); + config cfg = read_config(config_path.c_str()); // Create the manager for the FEL application. FreeElectronLaserManager manager(cfg); @@ -58,6 +63,10 @@ int main(int argc, char* argv[]) { manager.run(manager.getNt()); +#ifdef IPPL_ENABLE_CATALYST + manager.cat_viz.Finalize(); +#endif + msg << "End." << endl; IpplTimings::stopTimer(mainTimer); diff --git a/demos/fel/FreeElectronLaserManager.h b/demos/fel/FreeElectronLaserManager.h index 01b443385..7dbf296ef 100644 --- a/demos/fel/FreeElectronLaserManager.h +++ b/demos/fel/FreeElectronLaserManager.h @@ -2,6 +2,7 @@ #define IPPL_FREE_ELECTRON_LASER_MANAGER_H #include +#include #include #include #include @@ -16,10 +17,13 @@ #include "LorentzTransform.h" #include "MithraBunch.h" #include "Undulator.h" -#include "VideoWriter.h" #include "datatypes.h" #include "units.h" +#ifdef IPPL_ENABLE_CATALYST +#include "Stream/InSitu/CatalystAdaptor.h" +#endif + // FEL simulation manager. // // An electromagnetic (FDTD) PIC manager for the Free Electron Laser: it owns the @@ -53,7 +57,7 @@ class FreeElectronLaserManager : public ippl::BaseManager { , frame_m(ippl::UniaxialLorentzframe::from_gamma(frame_gamma_m)) , undulator_m(uparams_m, 2.0 * cfg.sigma_position[2] * frame_gamma_m * frame_gamma_m) {} - ~FreeElectronLaserManager() { video_m.close(); } + ~FreeElectronLaserManager() = default; protected: config m_config; @@ -84,7 +88,6 @@ class FreeElectronLaserManager : public ippl::BaseManager { ippl::undulator_parameters uparams_m; ///< Undulator parameters. ippl::UniaxialLorentzframe frame_m; ///< Boost into the co-moving frame (z-axis). ippl::Undulator undulator_m; ///< Static undulator field model. - FELVideoWriter video_m; ///< Optional ffmpeg Poynting-flux video. // --- narrow-band (resonant) radiation power diagnostic state --- // MITHRA reports the FEL output power as a sliding-window single-frequency @@ -98,6 +101,10 @@ class FreeElectronLaserManager : public ippl::BaseManager { Kokkos::View rp_fdt_m; ///< ring buffer [Nf][nx][ny][4] = (Ex,Ey,Bx,By)_lab public: +#ifdef IPPL_ENABLE_CATALYST + ippl::CatalystAdaptor cat_viz{std::string{"FreeElectronLaser"}}; +#endif + size_type getTotalP() const { return totalP_m; } void setTotalP(size_type totalP_) { totalP_m = totalP_; } @@ -247,6 +254,22 @@ class FreeElectronLaserManager : public ippl::BaseManager { void pre_run() override { Inform m("Pre Run"); + int outputDirectoryReady = 1; + if (ippl::Comm->rank() == 0) { + std::error_code error; + std::filesystem::create_directories(this->m_config.output_path, error); + if (error) { + outputDirectoryReady = 0; + m << "Unable to create output directory '" << this->m_config.output_path + << "': " << error.message() << endl; + } + } + MPI_Bcast(&outputDirectoryReady, 1, MPI_INT, 0, ippl::Comm->getCommunicator()); + if (!outputDirectoryReady) { + throw IpplException("FreeElectronLaserManager::pre_run", + "Unable to create the configured output directory"); + } + // The longitudinal box and the simulated time are measured in the // co-moving frame: stretch z and shorten the time accordingly. this->m_config.extents[2] *= frame_gamma_m; @@ -304,8 +327,15 @@ class FreeElectronLaserManager : public ippl::BaseManager { initializeParticles(); - // Open the ffmpeg pipe (rank 0, only if periodic output was requested). - video_m.open(this->m_config); +#ifdef IPPL_ENABLE_CATALYST + auto runtime_vis_registry = ippl::MakeVisRegistryRuntimePtr( + "Particles", this->pcontainer_m, + "EField", this->fcontainer_m->getE(), + "Bfield", this->fcontainer_m->getB() + ); + auto runtime_steer_registry = ippl::MakeVisRegistryRuntimePtr(); + cat_viz.Initialize(runtime_vis_registry, runtime_steer_registry); +#endif this->dump(); @@ -354,6 +384,12 @@ class FreeElectronLaserManager : public ippl::BaseManager { // 2. Advance the electromagnetic field one FDTD step. this->solver_m->solve(); +#ifdef IPPL_ENABLE_CATALYST + // Execute immediately after the solver has recomputed E/B. The adaptor + // makes its host snapshots synchronously during this call. + cat_viz.Execute(this->it_m, this->time_m); +#endif + // 3. Push particles with the self-consistent field plus the undulator // field transformed into the co-moving frame. auto und = undulator_m; @@ -369,14 +405,6 @@ class FreeElectronLaserManager : public ippl::BaseManager { dumpRadiation(); dumpRadiationBanded(); dumpFELDiagnostics(); - - // Emit a video frame on the configured rhythm. writeFrame is collective, - // so every rank must reach it under the same condition. - if (this->m_config.output_rhythm != 0 - && (this->it_m % (int)this->m_config.output_rhythm) == 0) { - video_m.writeFrame(this->it_m, *this->fcontainer_m, *this->pcontainer_m, frame_m, - this->m_config, this->nr_m); - } } // Radiated power leaving the downstream end of the domain, transformed back diff --git a/demos/fel/README.md b/demos/fel/README.md index 8db95d543..8629f6ba3 100644 --- a/demos/fel/README.md +++ b/demos/fel/README.md @@ -3,7 +3,7 @@ An electromagnetic PIC simulation of a free-electron laser: a relativistic electron bunch is tracked through an undulator in a co-moving Lorentz frame, with the self-consistent field advanced by an FDTD Maxwell solver. The radiated -power is written to a CSV (and, optionally, a Poynting-flux video). +power and FEL diagnostics are written to CSV files. ## Build @@ -12,23 +12,29 @@ cmake -S . -B build -DIPPL_ENABLE_FEL=ON -DCMAKE_CXX_STANDARD=20 cmake --build build --target FreeElectronLaser ``` +The FEL configuration is parsed with the Conduit API shipped by Catalyst. CMake +therefore finds or fetches Catalyst when the FEL demo is enabled, even when +`IPPL_ENABLE_CATALYST` itself is off. + The executable is built at -`build/fel/FreeElectronLaser`. +`build/demos/fel/FreeElectronLaser`. The example configuration is staged beside +it as `build/demos/fel/config.json` whenever the target is built. ## Run ```sh cd build -./fel/FreeElectronLaser ../fel/config.json --info 5 +./demos/fel/FreeElectronLaser --info 5 ``` -The argument is a MITHRA-style JSON job file (defaults to `../fel/config.json`); -see [config.json](config.json) for the available keys. Run on multiple ranks -with `mpirun -np ...`. +An optional first argument can select another MITHRA-style JSON job file. By +default, the executable uses the staged `build/demos/fel/config.json`; see +[config.json](config.json) for the available keys. Run on multiple ranks with +`mpirun -np ...`. Output is written to the directory given by `output.path` in the config: -`radiation_.csv` holds the radiated power. If `output.rhythm > 0`, a -Poynting-flux video is produced and requires **ffmpeg** on the `PATH`. +`radiation_.csv` holds the radiated power. The directory is created +automatically; relative paths are resolved from the process working directory. ## Acknowledgements diff --git a/demos/fel/VideoWriter.h b/demos/fel/VideoWriter.h deleted file mode 100644 index ba9664141..000000000 --- a/demos/fel/VideoWriter.h +++ /dev/null @@ -1,254 +0,0 @@ -#ifndef IPPL_FEL_VIDEO_WRITER_H -#define IPPL_FEL_VIDEO_WRITER_H - -// Poynting-flux visualization for the FEL simulation. -// -// Renders a slice of the (lab-frame) Poynting field plus the projected particle -// positions to a BMP frame and pipes the stream to ffmpeg. Ported from the -// original FreeElectronLaser.cpp main loop and encapsulated here so the manager -// and driver stay readable. stb_image_write's implementation is emitted by the -// translation unit that defines STB_IMAGE_WRITE_IMPLEMENTATION before including -// this header (FreeElectronLaser.cpp). - -#include -#include -#include // popen / pclose / FILE -#include -#include -#include - -#include -#include - -#include "Config.h" -#include "LorentzTransform.h" -#include "datatypes.h" -#include "units.h" - -// Google "Turbo" colormap (Anton Mikhailov), 256 RGB entries in [0,1]. -inline constexpr float turbo_cm[256][3] = { - {0.18995, 0.07176, 0.23217}, {0.19483, 0.08339, 0.26149}, {0.19956, 0.09498, 0.29024}, - {0.20415, 0.10652, 0.31844}, {0.20860, 0.11802, 0.34607}, {0.21291, 0.12947, 0.37314}, - {0.21708, 0.14087, 0.39964}, {0.22111, 0.15223, 0.42558}, {0.22500, 0.16354, 0.45096}, - {0.22875, 0.17481, 0.47578}, {0.23236, 0.18603, 0.50004}, {0.23582, 0.19720, 0.52373}, - {0.23915, 0.20833, 0.54686}, {0.24234, 0.21941, 0.56942}, {0.24539, 0.23044, 0.59142}, - {0.24830, 0.24143, 0.61286}, {0.25107, 0.25237, 0.63374}, {0.25369, 0.26327, 0.65406}, - {0.25618, 0.27412, 0.67381}, {0.25853, 0.28492, 0.69300}, {0.26074, 0.29568, 0.71162}, - {0.26280, 0.30639, 0.72968}, {0.26473, 0.31706, 0.74718}, {0.26652, 0.32768, 0.76412}, - {0.26816, 0.33825, 0.78050}, {0.26967, 0.34878, 0.79631}, {0.27103, 0.35926, 0.81156}, - {0.27226, 0.36970, 0.82624}, {0.27334, 0.38008, 0.84037}, {0.27429, 0.39043, 0.85393}, - {0.27509, 0.40072, 0.86692}, {0.27576, 0.41097, 0.87936}, {0.27628, 0.42118, 0.89123}, - {0.27667, 0.43134, 0.90254}, {0.27691, 0.44145, 0.91328}, {0.27701, 0.45152, 0.92347}, - {0.27698, 0.46153, 0.93309}, {0.27680, 0.47151, 0.94214}, {0.27648, 0.48144, 0.95064}, - {0.27603, 0.49132, 0.95857}, {0.27543, 0.50115, 0.96594}, {0.27469, 0.51094, 0.97275}, - {0.27381, 0.52069, 0.97899}, {0.27273, 0.53040, 0.98461}, {0.27106, 0.54015, 0.98930}, - {0.26878, 0.54995, 0.99303}, {0.26592, 0.55979, 0.99583}, {0.26252, 0.56967, 0.99773}, - {0.25862, 0.57958, 0.99876}, {0.25425, 0.58950, 0.99896}, {0.24946, 0.59943, 0.99835}, - {0.24427, 0.60937, 0.99697}, {0.23874, 0.61931, 0.99485}, {0.23288, 0.62923, 0.99202}, - {0.22676, 0.63913, 0.98851}, {0.22039, 0.64901, 0.98436}, {0.21382, 0.65886, 0.97959}, - {0.20708, 0.66866, 0.97423}, {0.20021, 0.67842, 0.96833}, {0.19326, 0.68812, 0.96190}, - {0.18625, 0.69775, 0.95498}, {0.17923, 0.70732, 0.94761}, {0.17223, 0.71680, 0.93981}, - {0.16529, 0.72620, 0.93161}, {0.15844, 0.73551, 0.92305}, {0.15173, 0.74472, 0.91416}, - {0.14519, 0.75381, 0.90496}, {0.13886, 0.76279, 0.89550}, {0.13278, 0.77165, 0.88580}, - {0.12698, 0.78037, 0.87590}, {0.12151, 0.78896, 0.86581}, {0.11639, 0.79740, 0.85559}, - {0.11167, 0.80569, 0.84525}, {0.10738, 0.81381, 0.83484}, {0.10357, 0.82177, 0.82437}, - {0.10026, 0.82955, 0.81389}, {0.09750, 0.83714, 0.80342}, {0.09532, 0.84455, 0.79299}, - {0.09377, 0.85175, 0.78264}, {0.09287, 0.85875, 0.77240}, {0.09267, 0.86554, 0.76230}, - {0.09320, 0.87211, 0.75237}, {0.09451, 0.87844, 0.74265}, {0.09662, 0.88454, 0.73316}, - {0.09958, 0.89040, 0.72393}, {0.10342, 0.89600, 0.71500}, {0.10815, 0.90142, 0.70599}, - {0.11374, 0.90673, 0.69651}, {0.12014, 0.91193, 0.68660}, {0.12733, 0.91701, 0.67627}, - {0.13526, 0.92197, 0.66556}, {0.14391, 0.92680, 0.65448}, {0.15323, 0.93151, 0.64308}, - {0.16319, 0.93609, 0.63137}, {0.17377, 0.94053, 0.61938}, {0.18491, 0.94484, 0.60713}, - {0.19659, 0.94901, 0.59466}, {0.20877, 0.95304, 0.58199}, {0.22142, 0.95692, 0.56914}, - {0.23449, 0.96065, 0.55614}, {0.24797, 0.96423, 0.54303}, {0.26180, 0.96765, 0.52981}, - {0.27597, 0.97092, 0.51653}, {0.29042, 0.97403, 0.50321}, {0.30513, 0.97697, 0.48987}, - {0.32006, 0.97974, 0.47654}, {0.33517, 0.98234, 0.46325}, {0.35043, 0.98477, 0.45002}, - {0.36581, 0.98702, 0.43688}, {0.38127, 0.98909, 0.42386}, {0.39678, 0.99098, 0.41098}, - {0.41229, 0.99268, 0.39826}, {0.42778, 0.99419, 0.38575}, {0.44321, 0.99551, 0.37345}, - {0.45854, 0.99663, 0.36140}, {0.47375, 0.99755, 0.34963}, {0.48879, 0.99828, 0.33816}, - {0.50362, 0.99879, 0.32701}, {0.51822, 0.99910, 0.31622}, {0.53255, 0.99919, 0.30581}, - {0.54658, 0.99907, 0.29581}, {0.56026, 0.99873, 0.28623}, {0.57357, 0.99817, 0.27712}, - {0.58646, 0.99739, 0.26849}, {0.59891, 0.99638, 0.26038}, {0.61088, 0.99514, 0.25280}, - {0.62233, 0.99366, 0.24579}, {0.63323, 0.99195, 0.23937}, {0.64362, 0.98999, 0.23356}, - {0.65394, 0.98775, 0.22835}, {0.66428, 0.98524, 0.22370}, {0.67462, 0.98246, 0.21960}, - {0.68494, 0.97941, 0.21602}, {0.69525, 0.97610, 0.21294}, {0.70553, 0.97255, 0.21032}, - {0.71577, 0.96875, 0.20815}, {0.72596, 0.96470, 0.20640}, {0.73610, 0.96043, 0.20504}, - {0.74617, 0.95593, 0.20406}, {0.75617, 0.95121, 0.20343}, {0.76608, 0.94627, 0.20311}, - {0.77591, 0.94113, 0.20310}, {0.78563, 0.93579, 0.20336}, {0.79524, 0.93025, 0.20386}, - {0.80473, 0.92452, 0.20459}, {0.81410, 0.91861, 0.20552}, {0.82333, 0.91253, 0.20663}, - {0.83241, 0.90627, 0.20788}, {0.84133, 0.89986, 0.20926}, {0.85010, 0.89328, 0.21074}, - {0.85868, 0.88655, 0.21230}, {0.86709, 0.87968, 0.21391}, {0.87530, 0.87267, 0.21555}, - {0.88331, 0.86553, 0.21719}, {0.89112, 0.85826, 0.21880}, {0.89870, 0.85087, 0.22038}, - {0.90605, 0.84337, 0.22188}, {0.91317, 0.83576, 0.22328}, {0.92004, 0.82806, 0.22456}, - {0.92666, 0.82025, 0.22570}, {0.93301, 0.81236, 0.22667}, {0.93909, 0.80439, 0.22744}, - {0.94489, 0.79634, 0.22800}, {0.95039, 0.78823, 0.22831}, {0.95560, 0.78005, 0.22836}, - {0.96049, 0.77181, 0.22811}, {0.96507, 0.76352, 0.22754}, {0.96931, 0.75519, 0.22663}, - {0.97323, 0.74682, 0.22536}, {0.97679, 0.73842, 0.22369}, {0.98000, 0.73000, 0.22161}, - {0.98289, 0.72140, 0.21918}, {0.98549, 0.71250, 0.21650}, {0.98781, 0.70330, 0.21358}, - {0.98986, 0.69382, 0.21043}, {0.99163, 0.68408, 0.20706}, {0.99314, 0.67408, 0.20348}, - {0.99438, 0.66386, 0.19971}, {0.99535, 0.65341, 0.19577}, {0.99607, 0.64277, 0.19165}, - {0.99654, 0.63193, 0.18738}, {0.99675, 0.62093, 0.18297}, {0.99672, 0.60977, 0.17842}, - {0.99644, 0.59846, 0.17376}, {0.99593, 0.58703, 0.16899}, {0.99517, 0.57549, 0.16412}, - {0.99419, 0.56386, 0.15918}, {0.99297, 0.55214, 0.15417}, {0.99153, 0.54036, 0.14910}, - {0.98987, 0.52854, 0.14398}, {0.98799, 0.51667, 0.13883}, {0.98590, 0.50479, 0.13367}, - {0.98360, 0.49291, 0.12849}, {0.98108, 0.48104, 0.12332}, {0.97837, 0.46920, 0.11817}, - {0.97545, 0.45740, 0.11305}, {0.97234, 0.44565, 0.10797}, {0.96904, 0.43399, 0.10294}, - {0.96555, 0.42241, 0.09798}, {0.96187, 0.41093, 0.09310}, {0.95801, 0.39958, 0.08831}, - {0.95398, 0.38836, 0.08362}, {0.94977, 0.37729, 0.07905}, {0.94538, 0.36638, 0.07461}, - {0.94084, 0.35566, 0.07031}, {0.93612, 0.34513, 0.06616}, {0.93125, 0.33482, 0.06218}, - {0.92623, 0.32473, 0.05837}, {0.92105, 0.31489, 0.05475}, {0.91572, 0.30530, 0.05134}, - {0.91024, 0.29599, 0.04814}, {0.90463, 0.28696, 0.04516}, {0.89888, 0.27824, 0.04243}, - {0.89298, 0.26981, 0.03993}, {0.88691, 0.26152, 0.03753}, {0.88066, 0.25334, 0.03521}, - {0.87422, 0.24526, 0.03297}, {0.86760, 0.23730, 0.03082}, {0.86079, 0.22945, 0.02875}, - {0.85380, 0.22170, 0.02677}, {0.84662, 0.21407, 0.02487}, {0.83926, 0.20654, 0.02305}, - {0.83172, 0.19912, 0.02131}, {0.82399, 0.19182, 0.01966}, {0.81608, 0.18462, 0.01809}, - {0.80799, 0.17753, 0.01660}, {0.79971, 0.17055, 0.01520}, {0.79125, 0.16368, 0.01387}, - {0.78260, 0.15693, 0.01264}, {0.77377, 0.15028, 0.01148}, {0.76476, 0.14374, 0.01041}, - {0.75556, 0.13731, 0.00942}, {0.74617, 0.13098, 0.00851}, {0.73661, 0.12477, 0.00769}, - {0.72686, 0.11867, 0.00695}, {0.71692, 0.11268, 0.00629}, {0.70680, 0.10680, 0.00571}, - {0.69650, 0.10102, 0.00522}, {0.68602, 0.09536, 0.00481}, {0.67535, 0.08980, 0.00449}, - {0.66449, 0.08436, 0.00424}, {0.65345, 0.07902, 0.00408}, {0.64223, 0.07380, 0.00401}, - {0.63082, 0.06868, 0.00401}, {0.61923, 0.06367, 0.00410}, {0.60746, 0.05878, 0.00427}, - {0.59550, 0.05399, 0.00453}, {0.58336, 0.04931, 0.00486}, {0.57103, 0.04474, 0.00529}, - {0.55852, 0.04028, 0.00579}, {0.54583, 0.03593, 0.00638}, {0.53295, 0.03169, 0.00705}, - {0.51989, 0.02756, 0.00780}, {0.50664, 0.02354, 0.00863}, {0.49321, 0.01963, 0.00955}, - {0.47960, 0.01583, 0.01055}}; - -// Write an RGB buffer as a BMP to an open FILE* (the ffmpeg pipe). -inline bool writeBMPToFD(FILE* fd, int width, int height, const unsigned char* data) { - const int channels = 3; // RGB - const int stride = width * channels; - std::vector flippedData(data, data + stride * height); - - if (!stbi_write_bmp_to_func( - [](void* context, void* data, int size) { - FILE* f = reinterpret_cast(context); - fwrite(data, 1, size, f); - }, - fd, width, height, channels, flippedData.data())) { - return false; - } - return true; -} - -// Encapsulates the ffmpeg pipe and per-frame rendering. -template -class FELVideoWriter { -public: - // Open the ffmpeg pipe on rank 0 when periodic output is requested. - void open(const config& cfg) { - if (ippl::Comm->rank() == 0 && cfg.output_rhythm != 0) { - const char* ffmpegCmd = - "ffmpeg -y -f image2pipe -vcodec bmp -r 30 -i - -vf " - "scale=force_original_aspect_ratio=decrease:force_divisible_by=2,format=yuv420p " - "-c:v libx264 -movflags +faststart ffmpeg_popen.mkv"; - ffmpeg_file_ = popen(ffmpegCmd, "w"); - } - } - - void close() { - if (ffmpeg_file_ != nullptr) { - pclose(ffmpeg_file_); - ffmpeg_file_ = nullptr; - } - } - - // Render and emit one frame for time step `it`. Collective over all ranks - // (the per-rank partial images are reduced onto rank 0 before encoding). - template - void writeFrame(int it, FieldContainer_t& fc, ParticleContainer_t& pc, const Frame_t& frame, - const config& cfg, const Vector_t& nr) { - const int rank = ippl::Comm->rank(); - const int size = ippl::Comm->size(); - - const int img_height = 400; - const int img_width = int(400.0 * cfg.extents[2] / cfg.extents[0]); - const int floatcount = img_width * img_height * 3; - - std::vector imagedata(floatcount, 0.0f); - std::vector recvbuffer(floatcount, 0.0f); - - auto& mesh = fc.getMesh(); - const auto origin = mesh.getOrigin(); - const auto ldom = fc.getFL().getLocalNDIndex(); - - // Project local particles onto the (z, x) plane as green dots. - auto phmirror = pc.R.getHostMirror(); - Kokkos::deep_copy(phmirror, pc.R.getView()); - for (size_t hi = 0; hi < pc.getLocalNum(); hi++) { - ippl::Vector ppos = phmirror(hi); - ppos -= origin; - ppos /= vector_cast(cfg.extents); - int x_imgcoord = ppos[2] * img_width; - int y_imgcoord = ppos[0] * img_height; - if (y_imgcoord >= 0 && x_imgcoord >= 0 && x_imgcoord < img_width - && y_imgcoord < img_height) { - const float intensity = - std::min(255.f, (img_width * img_height * 15.f) / cfg.num_particles); - float& g = imagedata[(y_imgcoord * img_width + x_imgcoord) * 3 + 1]; - g = std::min(255.f, g + intensity); - } - } - - // Color the lab-frame Poynting magnitude on the mid-y plane. - auto eh = fc.getE().getHostMirror(); - auto bh = fc.getB().getHostMirror(); - Kokkos::deep_copy(eh, fc.getE().getView()); - Kokkos::deep_copy(bh, fc.getB().getView()); - - for (int i = 1; i < img_width; i++) { - for (int j = 1; j < img_height; j++) { - int i_remap = (double(i) / (img_width - 1)) * (nr[2] - 4) + 2; - int j_remap = (double(j) / (img_height - 1)) * (nr[0] - 4) + 2; - if (i_remap >= ldom.first()[2] && i_remap <= ldom.last()[2] - && j_remap >= ldom.first()[0] && j_remap <= ldom.last()[0]) { - ippl::Vector E = eh(j_remap + 1 - ldom.first()[0], nr[1] / 2, - i_remap + 1 - ldom.first()[2]); - ippl::Vector B = bh(j_remap + 1 - ldom.first()[0], nr[1] / 2, - i_remap + 1 - ldom.first()[2]); - - auto eblab = frame.inverse_transform_EB( - Kokkos::make_pair, ippl::Vector>( - ippl::Vector(E), ippl::Vector(B))); - ippl::Vector poynting = ippl::cross(eblab.first, eblab.second); - - float normalized = std::sqrt(poynting.Pnorm()) * 0.00001f; - int index = (int)std::max(0.0f, std::min(normalized * 255.0f, 255.0f)); - for (int c = 0; c < 3; c++) { - imagedata[(j * img_width + i) * 3 + c] += turbo_cm[index][c] * 255.0f; - } - } - } - } - - // Butterfly reduction of the partial images onto rank 0. - int mask = 1; - while (mask < size) { - int partner = rank ^ mask; - if ((rank & mask) == 0) { - MPI_Recv(recvbuffer.data(), floatcount, MPI_FLOAT, partner, 0, - ippl::Comm->getCommunicator(), MPI_STATUS_IGNORE); - for (int f = 0; f < floatcount; f++) { - imagedata[f] += recvbuffer[f]; - } - } else { - MPI_Send(imagedata.data(), floatcount, MPI_FLOAT, partner, 0, - ippl::Comm->getCommunicator()); - } - mask <<= 1; - } - - if (rank == 0 && ffmpeg_file_ != nullptr) { - std::vector final_img(floatcount); - std::transform(imagedata.begin(), imagedata.end(), final_img.begin(), - [](float x) { return (uint8_t)std::min(255.0f, std::max(0.0f, x)); }); - writeBMPToFD(ffmpeg_file_, img_width, img_height, final_img.data()); - } - (void)it; - } - -private: - FILE* ffmpeg_file_ = nullptr; -}; - -#endif diff --git a/demos/fel/config.json b/demos/fel/config.json index d633a7514..bd1525872 100644 --- a/demos/fel/config.json +++ b/demos/fel/config.json @@ -26,7 +26,6 @@ } }, "output": { - "rhythm": 30, "path": "../renderdata" } -} +} \ No newline at end of file diff --git a/demos/fel/stb/README.md b/demos/fel/stb/README.md deleted file mode 100644 index 078dbc617..000000000 --- a/demos/fel/stb/README.md +++ /dev/null @@ -1,10 +0,0 @@ -# stb_image_write - -`stb_image_write.h` is vendored specifically for the FEL mini-app from the -upstream [`nothings/stb`](https://github.com/nothings/stb) repository at commit -`2c980bb59875b0d32144a71867fbdebb2f77cd20`. - -- Upstream file: `stb_image_write.h` v1.16 -- SHA-256: `cbd5f0ad7a9cf4468affb36354a1d2338034f2c12473cf1a8e32053cb6914a05` -- License: dual-licensed under the MIT License or public domain; the complete - license text is included in the header. diff --git a/demos/fel/stb/stb_image_write.h b/demos/fel/stb/stb_image_write.h deleted file mode 100644 index e4b32ed1b..000000000 --- a/demos/fel/stb/stb_image_write.h +++ /dev/null @@ -1,1724 +0,0 @@ -/* stb_image_write - v1.16 - public domain - http://nothings.org/stb - writes out PNG/BMP/TGA/JPEG/HDR images to C stdio - Sean Barrett 2010-2015 - no warranty implied; use at your own risk - - Before #including, - - #define STB_IMAGE_WRITE_IMPLEMENTATION - - in the file that you want to have the implementation. - - Will probably not work correctly with strict-aliasing optimizations. - -ABOUT: - - This header file is a library for writing images to C stdio or a callback. - - The PNG output is not optimal; it is 20-50% larger than the file - written by a decent optimizing implementation; though providing a custom - zlib compress function (see STBIW_ZLIB_COMPRESS) can mitigate that. - This library is designed for source code compactness and simplicity, - not optimal image file size or run-time performance. - -BUILDING: - - You can #define STBIW_ASSERT(x) before the #include to avoid using assert.h. - You can #define STBIW_MALLOC(), STBIW_REALLOC(), and STBIW_FREE() to replace - malloc,realloc,free. - You can #define STBIW_MEMMOVE() to replace memmove() - You can #define STBIW_ZLIB_COMPRESS to use a custom zlib-style compress function - for PNG compression (instead of the builtin one), it must have the following signature: - unsigned char * my_compress(unsigned char *data, int data_len, int *out_len, int quality); - The returned data will be freed with STBIW_FREE() (free() by default), - so it must be heap allocated with STBIW_MALLOC() (malloc() by default), - -UNICODE: - - If compiling for Windows and you wish to use Unicode filenames, compile - with - #define STBIW_WINDOWS_UTF8 - and pass utf8-encoded filenames. Call stbiw_convert_wchar_to_utf8 to convert - Windows wchar_t filenames to utf8. - -USAGE: - - There are five functions, one for each image file format: - - int stbi_write_png(char const *filename, int w, int h, int comp, const void *data, int stride_in_bytes); - int stbi_write_bmp(char const *filename, int w, int h, int comp, const void *data); - int stbi_write_tga(char const *filename, int w, int h, int comp, const void *data); - int stbi_write_jpg(char const *filename, int w, int h, int comp, const void *data, int quality); - int stbi_write_hdr(char const *filename, int w, int h, int comp, const float *data); - - void stbi_flip_vertically_on_write(int flag); // flag is non-zero to flip data vertically - - There are also five equivalent functions that use an arbitrary write function. You are - expected to open/close your file-equivalent before and after calling these: - - int stbi_write_png_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data, int stride_in_bytes); - int stbi_write_bmp_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data); - int stbi_write_tga_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data); - int stbi_write_hdr_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const float *data); - int stbi_write_jpg_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int quality); - - where the callback is: - void stbi_write_func(void *context, void *data, int size); - - You can configure it with these global variables: - int stbi_write_tga_with_rle; // defaults to true; set to 0 to disable RLE - int stbi_write_png_compression_level; // defaults to 8; set to higher for more compression - int stbi_write_force_png_filter; // defaults to -1; set to 0..5 to force a filter mode - - - You can define STBI_WRITE_NO_STDIO to disable the file variant of these - functions, so the library will not use stdio.h at all. However, this will - also disable HDR writing, because it requires stdio for formatted output. - - Each function returns 0 on failure and non-0 on success. - - The functions create an image file defined by the parameters. The image - is a rectangle of pixels stored from left-to-right, top-to-bottom. - Each pixel contains 'comp' channels of data stored interleaved with 8-bits - per channel, in the following order: 1=Y, 2=YA, 3=RGB, 4=RGBA. (Y is - monochrome color.) The rectangle is 'w' pixels wide and 'h' pixels tall. - The *data pointer points to the first byte of the top-left-most pixel. - For PNG, "stride_in_bytes" is the distance in bytes from the first byte of - a row of pixels to the first byte of the next row of pixels. - - PNG creates output files with the same number of components as the input. - The BMP format expands Y to RGB in the file format and does not - output alpha. - - PNG supports writing rectangles of data even when the bytes storing rows of - data are not consecutive in memory (e.g. sub-rectangles of a larger image), - by supplying the stride between the beginning of adjacent rows. The other - formats do not. (Thus you cannot write a native-format BMP through the BMP - writer, both because it is in BGR order and because it may have padding - at the end of the line.) - - PNG allows you to set the deflate compression level by setting the global - variable 'stbi_write_png_compression_level' (it defaults to 8). - - HDR expects linear float data. Since the format is always 32-bit rgb(e) - data, alpha (if provided) is discarded, and for monochrome data it is - replicated across all three channels. - - TGA supports RLE or non-RLE compressed data. To use non-RLE-compressed - data, set the global variable 'stbi_write_tga_with_rle' to 0. - - JPEG does ignore alpha channels in input data; quality is between 1 and 100. - Higher quality looks better but results in a bigger image. - JPEG baseline (no JPEG progressive). - -CREDITS: - - - Sean Barrett - PNG/BMP/TGA - Baldur Karlsson - HDR - Jean-Sebastien Guay - TGA monochrome - Tim Kelsey - misc enhancements - Alan Hickman - TGA RLE - Emmanuel Julien - initial file IO callback implementation - Jon Olick - original jo_jpeg.cpp code - Daniel Gibson - integrate JPEG, allow external zlib - Aarni Koskela - allow choosing PNG filter - - bugfixes: - github:Chribba - Guillaume Chereau - github:jry2 - github:romigrou - Sergio Gonzalez - Jonas Karlsson - Filip Wasil - Thatcher Ulrich - github:poppolopoppo - Patrick Boettcher - github:xeekworx - Cap Petschulat - Simon Rodriguez - Ivan Tikhonov - github:ignotion - Adam Schackart - Andrew Kensler - -LICENSE - - See end of file for license information. - -*/ - -#ifndef INCLUDE_STB_IMAGE_WRITE_H -#define INCLUDE_STB_IMAGE_WRITE_H - -#include - -// if STB_IMAGE_WRITE_STATIC causes problems, try defining STBIWDEF to 'inline' or 'static inline' -#ifndef STBIWDEF -#ifdef STB_IMAGE_WRITE_STATIC -#define STBIWDEF static -#else -#ifdef __cplusplus -#define STBIWDEF extern "C" -#else -#define STBIWDEF extern -#endif -#endif -#endif - -#ifndef STB_IMAGE_WRITE_STATIC // C++ forbids static forward declarations -STBIWDEF int stbi_write_tga_with_rle; -STBIWDEF int stbi_write_png_compression_level; -STBIWDEF int stbi_write_force_png_filter; -#endif - -#ifndef STBI_WRITE_NO_STDIO -STBIWDEF int stbi_write_png(char const *filename, int w, int h, int comp, const void *data, int stride_in_bytes); -STBIWDEF int stbi_write_bmp(char const *filename, int w, int h, int comp, const void *data); -STBIWDEF int stbi_write_tga(char const *filename, int w, int h, int comp, const void *data); -STBIWDEF int stbi_write_hdr(char const *filename, int w, int h, int comp, const float *data); -STBIWDEF int stbi_write_jpg(char const *filename, int x, int y, int comp, const void *data, int quality); - -#ifdef STBIW_WINDOWS_UTF8 -STBIWDEF int stbiw_convert_wchar_to_utf8(char *buffer, size_t bufferlen, const wchar_t* input); -#endif -#endif - -typedef void stbi_write_func(void *context, void *data, int size); - -STBIWDEF int stbi_write_png_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data, int stride_in_bytes); -STBIWDEF int stbi_write_bmp_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data); -STBIWDEF int stbi_write_tga_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data); -STBIWDEF int stbi_write_hdr_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const float *data); -STBIWDEF int stbi_write_jpg_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int quality); - -STBIWDEF void stbi_flip_vertically_on_write(int flip_boolean); - -#endif//INCLUDE_STB_IMAGE_WRITE_H - -#ifdef STB_IMAGE_WRITE_IMPLEMENTATION - -#ifdef _WIN32 - #ifndef _CRT_SECURE_NO_WARNINGS - #define _CRT_SECURE_NO_WARNINGS - #endif - #ifndef _CRT_NONSTDC_NO_DEPRECATE - #define _CRT_NONSTDC_NO_DEPRECATE - #endif -#endif - -#ifndef STBI_WRITE_NO_STDIO -#include -#endif // STBI_WRITE_NO_STDIO - -#include -#include -#include -#include - -#if defined(STBIW_MALLOC) && defined(STBIW_FREE) && (defined(STBIW_REALLOC) || defined(STBIW_REALLOC_SIZED)) -// ok -#elif !defined(STBIW_MALLOC) && !defined(STBIW_FREE) && !defined(STBIW_REALLOC) && !defined(STBIW_REALLOC_SIZED) -// ok -#else -#error "Must define all or none of STBIW_MALLOC, STBIW_FREE, and STBIW_REALLOC (or STBIW_REALLOC_SIZED)." -#endif - -#ifndef STBIW_MALLOC -#define STBIW_MALLOC(sz) malloc(sz) -#define STBIW_REALLOC(p,newsz) realloc(p,newsz) -#define STBIW_FREE(p) free(p) -#endif - -#ifndef STBIW_REALLOC_SIZED -#define STBIW_REALLOC_SIZED(p,oldsz,newsz) STBIW_REALLOC(p,newsz) -#endif - - -#ifndef STBIW_MEMMOVE -#define STBIW_MEMMOVE(a,b,sz) memmove(a,b,sz) -#endif - - -#ifndef STBIW_ASSERT -#include -#define STBIW_ASSERT(x) assert(x) -#endif - -#define STBIW_UCHAR(x) (unsigned char) ((x) & 0xff) - -#ifdef STB_IMAGE_WRITE_STATIC -static int stbi_write_png_compression_level = 8; -static int stbi_write_tga_with_rle = 1; -static int stbi_write_force_png_filter = -1; -#else -int stbi_write_png_compression_level = 8; -int stbi_write_tga_with_rle = 1; -int stbi_write_force_png_filter = -1; -#endif - -static int stbi__flip_vertically_on_write = 0; - -STBIWDEF void stbi_flip_vertically_on_write(int flag) -{ - stbi__flip_vertically_on_write = flag; -} - -typedef struct -{ - stbi_write_func *func; - void *context; - unsigned char buffer[64]; - int buf_used; -} stbi__write_context; - -// initialize a callback-based context -static void stbi__start_write_callbacks(stbi__write_context *s, stbi_write_func *c, void *context) -{ - s->func = c; - s->context = context; -} - -#ifndef STBI_WRITE_NO_STDIO - -static void stbi__stdio_write(void *context, void *data, int size) -{ - fwrite(data,1,size,(FILE*) context); -} - -#if defined(_WIN32) && defined(STBIW_WINDOWS_UTF8) -#ifdef __cplusplus -#define STBIW_EXTERN extern "C" -#else -#define STBIW_EXTERN extern -#endif -STBIW_EXTERN __declspec(dllimport) int __stdcall MultiByteToWideChar(unsigned int cp, unsigned long flags, const char *str, int cbmb, wchar_t *widestr, int cchwide); -STBIW_EXTERN __declspec(dllimport) int __stdcall WideCharToMultiByte(unsigned int cp, unsigned long flags, const wchar_t *widestr, int cchwide, char *str, int cbmb, const char *defchar, int *used_default); - -STBIWDEF int stbiw_convert_wchar_to_utf8(char *buffer, size_t bufferlen, const wchar_t* input) -{ - return WideCharToMultiByte(65001 /* UTF8 */, 0, input, -1, buffer, (int) bufferlen, NULL, NULL); -} -#endif - -static FILE *stbiw__fopen(char const *filename, char const *mode) -{ - FILE *f; -#if defined(_WIN32) && defined(STBIW_WINDOWS_UTF8) - wchar_t wMode[64]; - wchar_t wFilename[1024]; - if (0 == MultiByteToWideChar(65001 /* UTF8 */, 0, filename, -1, wFilename, sizeof(wFilename)/sizeof(*wFilename))) - return 0; - - if (0 == MultiByteToWideChar(65001 /* UTF8 */, 0, mode, -1, wMode, sizeof(wMode)/sizeof(*wMode))) - return 0; - -#if defined(_MSC_VER) && _MSC_VER >= 1400 - if (0 != _wfopen_s(&f, wFilename, wMode)) - f = 0; -#else - f = _wfopen(wFilename, wMode); -#endif - -#elif defined(_MSC_VER) && _MSC_VER >= 1400 - if (0 != fopen_s(&f, filename, mode)) - f=0; -#else - f = fopen(filename, mode); -#endif - return f; -} - -static int stbi__start_write_file(stbi__write_context *s, const char *filename) -{ - FILE *f = stbiw__fopen(filename, "wb"); - stbi__start_write_callbacks(s, stbi__stdio_write, (void *) f); - return f != NULL; -} - -static void stbi__end_write_file(stbi__write_context *s) -{ - fclose((FILE *)s->context); -} - -#endif // !STBI_WRITE_NO_STDIO - -typedef unsigned int stbiw_uint32; -typedef int stb_image_write_test[sizeof(stbiw_uint32)==4 ? 1 : -1]; - -static void stbiw__writefv(stbi__write_context *s, const char *fmt, va_list v) -{ - while (*fmt) { - switch (*fmt++) { - case ' ': break; - case '1': { unsigned char x = STBIW_UCHAR(va_arg(v, int)); - s->func(s->context,&x,1); - break; } - case '2': { int x = va_arg(v,int); - unsigned char b[2]; - b[0] = STBIW_UCHAR(x); - b[1] = STBIW_UCHAR(x>>8); - s->func(s->context,b,2); - break; } - case '4': { stbiw_uint32 x = va_arg(v,int); - unsigned char b[4]; - b[0]=STBIW_UCHAR(x); - b[1]=STBIW_UCHAR(x>>8); - b[2]=STBIW_UCHAR(x>>16); - b[3]=STBIW_UCHAR(x>>24); - s->func(s->context,b,4); - break; } - default: - STBIW_ASSERT(0); - return; - } - } -} - -static void stbiw__writef(stbi__write_context *s, const char *fmt, ...) -{ - va_list v; - va_start(v, fmt); - stbiw__writefv(s, fmt, v); - va_end(v); -} - -static void stbiw__write_flush(stbi__write_context *s) -{ - if (s->buf_used) { - s->func(s->context, &s->buffer, s->buf_used); - s->buf_used = 0; - } -} - -static void stbiw__putc(stbi__write_context *s, unsigned char c) -{ - s->func(s->context, &c, 1); -} - -static void stbiw__write1(stbi__write_context *s, unsigned char a) -{ - if ((size_t)s->buf_used + 1 > sizeof(s->buffer)) - stbiw__write_flush(s); - s->buffer[s->buf_used++] = a; -} - -static void stbiw__write3(stbi__write_context *s, unsigned char a, unsigned char b, unsigned char c) -{ - int n; - if ((size_t)s->buf_used + 3 > sizeof(s->buffer)) - stbiw__write_flush(s); - n = s->buf_used; - s->buf_used = n+3; - s->buffer[n+0] = a; - s->buffer[n+1] = b; - s->buffer[n+2] = c; -} - -static void stbiw__write_pixel(stbi__write_context *s, int rgb_dir, int comp, int write_alpha, int expand_mono, unsigned char *d) -{ - unsigned char bg[3] = { 255, 0, 255}, px[3]; - int k; - - if (write_alpha < 0) - stbiw__write1(s, d[comp - 1]); - - switch (comp) { - case 2: // 2 pixels = mono + alpha, alpha is written separately, so same as 1-channel case - case 1: - if (expand_mono) - stbiw__write3(s, d[0], d[0], d[0]); // monochrome bmp - else - stbiw__write1(s, d[0]); // monochrome TGA - break; - case 4: - if (!write_alpha) { - // composite against pink background - for (k = 0; k < 3; ++k) - px[k] = bg[k] + ((d[k] - bg[k]) * d[3]) / 255; - stbiw__write3(s, px[1 - rgb_dir], px[1], px[1 + rgb_dir]); - break; - } - /* FALLTHROUGH */ - case 3: - stbiw__write3(s, d[1 - rgb_dir], d[1], d[1 + rgb_dir]); - break; - } - if (write_alpha > 0) - stbiw__write1(s, d[comp - 1]); -} - -static void stbiw__write_pixels(stbi__write_context *s, int rgb_dir, int vdir, int x, int y, int comp, void *data, int write_alpha, int scanline_pad, int expand_mono) -{ - stbiw_uint32 zero = 0; - int i,j, j_end; - - if (y <= 0) - return; - - if (stbi__flip_vertically_on_write) - vdir *= -1; - - if (vdir < 0) { - j_end = -1; j = y-1; - } else { - j_end = y; j = 0; - } - - for (; j != j_end; j += vdir) { - for (i=0; i < x; ++i) { - unsigned char *d = (unsigned char *) data + (j*x+i)*comp; - stbiw__write_pixel(s, rgb_dir, comp, write_alpha, expand_mono, d); - } - stbiw__write_flush(s); - s->func(s->context, &zero, scanline_pad); - } -} - -static int stbiw__outfile(stbi__write_context *s, int rgb_dir, int vdir, int x, int y, int comp, int expand_mono, void *data, int alpha, int pad, const char *fmt, ...) -{ - if (y < 0 || x < 0) { - return 0; - } else { - va_list v; - va_start(v, fmt); - stbiw__writefv(s, fmt, v); - va_end(v); - stbiw__write_pixels(s,rgb_dir,vdir,x,y,comp,data,alpha,pad, expand_mono); - return 1; - } -} - -static int stbi_write_bmp_core(stbi__write_context *s, int x, int y, int comp, const void *data) -{ - if (comp != 4) { - // write RGB bitmap - int pad = (-x*3) & 3; - return stbiw__outfile(s,-1,-1,x,y,comp,1,(void *) data,0,pad, - "11 4 22 4" "4 44 22 444444", - 'B', 'M', 14+40+(x*3+pad)*y, 0,0, 14+40, // file header - 40, x,y, 1,24, 0,0,0,0,0,0); // bitmap header - } else { - // RGBA bitmaps need a v4 header - // use BI_BITFIELDS mode with 32bpp and alpha mask - // (straight BI_RGB with alpha mask doesn't work in most readers) - return stbiw__outfile(s,-1,-1,x,y,comp,1,(void *)data,1,0, - "11 4 22 4" "4 44 22 444444 4444 4 444 444 444 444", - 'B', 'M', 14+108+x*y*4, 0, 0, 14+108, // file header - 108, x,y, 1,32, 3,0,0,0,0,0, 0xff0000,0xff00,0xff,0xff000000u, 0, 0,0,0, 0,0,0, 0,0,0, 0,0,0); // bitmap V4 header - } -} - -STBIWDEF int stbi_write_bmp_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data) -{ - stbi__write_context s = { 0 }; - stbi__start_write_callbacks(&s, func, context); - return stbi_write_bmp_core(&s, x, y, comp, data); -} - -#ifndef STBI_WRITE_NO_STDIO -STBIWDEF int stbi_write_bmp(char const *filename, int x, int y, int comp, const void *data) -{ - stbi__write_context s = { 0 }; - if (stbi__start_write_file(&s,filename)) { - int r = stbi_write_bmp_core(&s, x, y, comp, data); - stbi__end_write_file(&s); - return r; - } else - return 0; -} -#endif //!STBI_WRITE_NO_STDIO - -static int stbi_write_tga_core(stbi__write_context *s, int x, int y, int comp, void *data) -{ - int has_alpha = (comp == 2 || comp == 4); - int colorbytes = has_alpha ? comp-1 : comp; - int format = colorbytes < 2 ? 3 : 2; // 3 color channels (RGB/RGBA) = 2, 1 color channel (Y/YA) = 3 - - if (y < 0 || x < 0) - return 0; - - if (!stbi_write_tga_with_rle) { - return stbiw__outfile(s, -1, -1, x, y, comp, 0, (void *) data, has_alpha, 0, - "111 221 2222 11", 0, 0, format, 0, 0, 0, 0, 0, x, y, (colorbytes + has_alpha) * 8, has_alpha * 8); - } else { - int i,j,k; - int jend, jdir; - - stbiw__writef(s, "111 221 2222 11", 0,0,format+8, 0,0,0, 0,0,x,y, (colorbytes + has_alpha) * 8, has_alpha * 8); - - if (stbi__flip_vertically_on_write) { - j = 0; - jend = y; - jdir = 1; - } else { - j = y-1; - jend = -1; - jdir = -1; - } - for (; j != jend; j += jdir) { - unsigned char *row = (unsigned char *) data + j * x * comp; - int len; - - for (i = 0; i < x; i += len) { - unsigned char *begin = row + i * comp; - int diff = 1; - len = 1; - - if (i < x - 1) { - ++len; - diff = memcmp(begin, row + (i + 1) * comp, comp); - if (diff) { - const unsigned char *prev = begin; - for (k = i + 2; k < x && len < 128; ++k) { - if (memcmp(prev, row + k * comp, comp)) { - prev += comp; - ++len; - } else { - --len; - break; - } - } - } else { - for (k = i + 2; k < x && len < 128; ++k) { - if (!memcmp(begin, row + k * comp, comp)) { - ++len; - } else { - break; - } - } - } - } - - if (diff) { - unsigned char header = STBIW_UCHAR(len - 1); - stbiw__write1(s, header); - for (k = 0; k < len; ++k) { - stbiw__write_pixel(s, -1, comp, has_alpha, 0, begin + k * comp); - } - } else { - unsigned char header = STBIW_UCHAR(len - 129); - stbiw__write1(s, header); - stbiw__write_pixel(s, -1, comp, has_alpha, 0, begin); - } - } - } - stbiw__write_flush(s); - } - return 1; -} - -STBIWDEF int stbi_write_tga_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data) -{ - stbi__write_context s = { 0 }; - stbi__start_write_callbacks(&s, func, context); - return stbi_write_tga_core(&s, x, y, comp, (void *) data); -} - -#ifndef STBI_WRITE_NO_STDIO -STBIWDEF int stbi_write_tga(char const *filename, int x, int y, int comp, const void *data) -{ - stbi__write_context s = { 0 }; - if (stbi__start_write_file(&s,filename)) { - int r = stbi_write_tga_core(&s, x, y, comp, (void *) data); - stbi__end_write_file(&s); - return r; - } else - return 0; -} -#endif - -// ************************************************************************************************* -// Radiance RGBE HDR writer -// by Baldur Karlsson - -#define stbiw__max(a, b) ((a) > (b) ? (a) : (b)) - -#ifndef STBI_WRITE_NO_STDIO - -static void stbiw__linear_to_rgbe(unsigned char *rgbe, float *linear) -{ - int exponent; - float maxcomp = stbiw__max(linear[0], stbiw__max(linear[1], linear[2])); - - if (maxcomp < 1e-32f) { - rgbe[0] = rgbe[1] = rgbe[2] = rgbe[3] = 0; - } else { - float normalize = (float) frexp(maxcomp, &exponent) * 256.0f/maxcomp; - - rgbe[0] = (unsigned char)(linear[0] * normalize); - rgbe[1] = (unsigned char)(linear[1] * normalize); - rgbe[2] = (unsigned char)(linear[2] * normalize); - rgbe[3] = (unsigned char)(exponent + 128); - } -} - -static void stbiw__write_run_data(stbi__write_context *s, int length, unsigned char databyte) -{ - unsigned char lengthbyte = STBIW_UCHAR(length+128); - STBIW_ASSERT(length+128 <= 255); - s->func(s->context, &lengthbyte, 1); - s->func(s->context, &databyte, 1); -} - -static void stbiw__write_dump_data(stbi__write_context *s, int length, unsigned char *data) -{ - unsigned char lengthbyte = STBIW_UCHAR(length); - STBIW_ASSERT(length <= 128); // inconsistent with spec but consistent with official code - s->func(s->context, &lengthbyte, 1); - s->func(s->context, data, length); -} - -static void stbiw__write_hdr_scanline(stbi__write_context *s, int width, int ncomp, unsigned char *scratch, float *scanline) -{ - unsigned char scanlineheader[4] = { 2, 2, 0, 0 }; - unsigned char rgbe[4]; - float linear[3]; - int x; - - scanlineheader[2] = (width&0xff00)>>8; - scanlineheader[3] = (width&0x00ff); - - /* skip RLE for images too small or large */ - if (width < 8 || width >= 32768) { - for (x=0; x < width; x++) { - switch (ncomp) { - case 4: /* fallthrough */ - case 3: linear[2] = scanline[x*ncomp + 2]; - linear[1] = scanline[x*ncomp + 1]; - linear[0] = scanline[x*ncomp + 0]; - break; - default: - linear[0] = linear[1] = linear[2] = scanline[x*ncomp + 0]; - break; - } - stbiw__linear_to_rgbe(rgbe, linear); - s->func(s->context, rgbe, 4); - } - } else { - int c,r; - /* encode into scratch buffer */ - for (x=0; x < width; x++) { - switch(ncomp) { - case 4: /* fallthrough */ - case 3: linear[2] = scanline[x*ncomp + 2]; - linear[1] = scanline[x*ncomp + 1]; - linear[0] = scanline[x*ncomp + 0]; - break; - default: - linear[0] = linear[1] = linear[2] = scanline[x*ncomp + 0]; - break; - } - stbiw__linear_to_rgbe(rgbe, linear); - scratch[x + width*0] = rgbe[0]; - scratch[x + width*1] = rgbe[1]; - scratch[x + width*2] = rgbe[2]; - scratch[x + width*3] = rgbe[3]; - } - - s->func(s->context, scanlineheader, 4); - - /* RLE each component separately */ - for (c=0; c < 4; c++) { - unsigned char *comp = &scratch[width*c]; - - x = 0; - while (x < width) { - // find first run - r = x; - while (r+2 < width) { - if (comp[r] == comp[r+1] && comp[r] == comp[r+2]) - break; - ++r; - } - if (r+2 >= width) - r = width; - // dump up to first run - while (x < r) { - int len = r-x; - if (len > 128) len = 128; - stbiw__write_dump_data(s, len, &comp[x]); - x += len; - } - // if there's a run, output it - if (r+2 < width) { // same test as what we break out of in search loop, so only true if we break'd - // find next byte after run - while (r < width && comp[r] == comp[x]) - ++r; - // output run up to r - while (x < r) { - int len = r-x; - if (len > 127) len = 127; - stbiw__write_run_data(s, len, comp[x]); - x += len; - } - } - } - } - } -} - -static int stbi_write_hdr_core(stbi__write_context *s, int x, int y, int comp, float *data) -{ - if (y <= 0 || x <= 0 || data == NULL) - return 0; - else { - // Each component is stored separately. Allocate scratch space for full output scanline. - unsigned char *scratch = (unsigned char *) STBIW_MALLOC(x*4); - int i, len; - char buffer[128]; - char header[] = "#?RADIANCE\n# Written by stb_image_write.h\nFORMAT=32-bit_rle_rgbe\n"; - s->func(s->context, header, sizeof(header)-1); - -#ifdef __STDC_LIB_EXT1__ - len = sprintf_s(buffer, sizeof(buffer), "EXPOSURE= 1.0000000000000\n\n-Y %d +X %d\n", y, x); -#else - len = sprintf(buffer, "EXPOSURE= 1.0000000000000\n\n-Y %d +X %d\n", y, x); -#endif - s->func(s->context, buffer, len); - - for(i=0; i < y; i++) - stbiw__write_hdr_scanline(s, x, comp, scratch, data + comp*x*(stbi__flip_vertically_on_write ? y-1-i : i)); - STBIW_FREE(scratch); - return 1; - } -} - -STBIWDEF int stbi_write_hdr_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const float *data) -{ - stbi__write_context s = { 0 }; - stbi__start_write_callbacks(&s, func, context); - return stbi_write_hdr_core(&s, x, y, comp, (float *) data); -} - -STBIWDEF int stbi_write_hdr(char const *filename, int x, int y, int comp, const float *data) -{ - stbi__write_context s = { 0 }; - if (stbi__start_write_file(&s,filename)) { - int r = stbi_write_hdr_core(&s, x, y, comp, (float *) data); - stbi__end_write_file(&s); - return r; - } else - return 0; -} -#endif // STBI_WRITE_NO_STDIO - - -////////////////////////////////////////////////////////////////////////////// -// -// PNG writer -// - -#ifndef STBIW_ZLIB_COMPRESS -// stretchy buffer; stbiw__sbpush() == vector<>::push_back() -- stbiw__sbcount() == vector<>::size() -#define stbiw__sbraw(a) ((int *) (void *) (a) - 2) -#define stbiw__sbm(a) stbiw__sbraw(a)[0] -#define stbiw__sbn(a) stbiw__sbraw(a)[1] - -#define stbiw__sbneedgrow(a,n) ((a)==0 || stbiw__sbn(a)+n >= stbiw__sbm(a)) -#define stbiw__sbmaybegrow(a,n) (stbiw__sbneedgrow(a,(n)) ? stbiw__sbgrow(a,n) : 0) -#define stbiw__sbgrow(a,n) stbiw__sbgrowf((void **) &(a), (n), sizeof(*(a))) - -#define stbiw__sbpush(a, v) (stbiw__sbmaybegrow(a,1), (a)[stbiw__sbn(a)++] = (v)) -#define stbiw__sbcount(a) ((a) ? stbiw__sbn(a) : 0) -#define stbiw__sbfree(a) ((a) ? STBIW_FREE(stbiw__sbraw(a)),0 : 0) - -static void *stbiw__sbgrowf(void **arr, int increment, int itemsize) -{ - int m = *arr ? 2*stbiw__sbm(*arr)+increment : increment+1; - void *p = STBIW_REALLOC_SIZED(*arr ? stbiw__sbraw(*arr) : 0, *arr ? (stbiw__sbm(*arr)*itemsize + sizeof(int)*2) : 0, itemsize * m + sizeof(int)*2); - STBIW_ASSERT(p); - if (p) { - if (!*arr) ((int *) p)[1] = 0; - *arr = (void *) ((int *) p + 2); - stbiw__sbm(*arr) = m; - } - return *arr; -} - -static unsigned char *stbiw__zlib_flushf(unsigned char *data, unsigned int *bitbuffer, int *bitcount) -{ - while (*bitcount >= 8) { - stbiw__sbpush(data, STBIW_UCHAR(*bitbuffer)); - *bitbuffer >>= 8; - *bitcount -= 8; - } - return data; -} - -static int stbiw__zlib_bitrev(int code, int codebits) -{ - int res=0; - while (codebits--) { - res = (res << 1) | (code & 1); - code >>= 1; - } - return res; -} - -static unsigned int stbiw__zlib_countm(unsigned char *a, unsigned char *b, int limit) -{ - int i; - for (i=0; i < limit && i < 258; ++i) - if (a[i] != b[i]) break; - return i; -} - -static unsigned int stbiw__zhash(unsigned char *data) -{ - stbiw_uint32 hash = data[0] + (data[1] << 8) + (data[2] << 16); - hash ^= hash << 3; - hash += hash >> 5; - hash ^= hash << 4; - hash += hash >> 17; - hash ^= hash << 25; - hash += hash >> 6; - return hash; -} - -#define stbiw__zlib_flush() (out = stbiw__zlib_flushf(out, &bitbuf, &bitcount)) -#define stbiw__zlib_add(code,codebits) \ - (bitbuf |= (code) << bitcount, bitcount += (codebits), stbiw__zlib_flush()) -#define stbiw__zlib_huffa(b,c) stbiw__zlib_add(stbiw__zlib_bitrev(b,c),c) -// default huffman tables -#define stbiw__zlib_huff1(n) stbiw__zlib_huffa(0x30 + (n), 8) -#define stbiw__zlib_huff2(n) stbiw__zlib_huffa(0x190 + (n)-144, 9) -#define stbiw__zlib_huff3(n) stbiw__zlib_huffa(0 + (n)-256,7) -#define stbiw__zlib_huff4(n) stbiw__zlib_huffa(0xc0 + (n)-280,8) -#define stbiw__zlib_huff(n) ((n) <= 143 ? stbiw__zlib_huff1(n) : (n) <= 255 ? stbiw__zlib_huff2(n) : (n) <= 279 ? stbiw__zlib_huff3(n) : stbiw__zlib_huff4(n)) -#define stbiw__zlib_huffb(n) ((n) <= 143 ? stbiw__zlib_huff1(n) : stbiw__zlib_huff2(n)) - -#define stbiw__ZHASH 16384 - -#endif // STBIW_ZLIB_COMPRESS - -STBIWDEF unsigned char * stbi_zlib_compress(unsigned char *data, int data_len, int *out_len, int quality) -{ -#ifdef STBIW_ZLIB_COMPRESS - // user provided a zlib compress implementation, use that - return STBIW_ZLIB_COMPRESS(data, data_len, out_len, quality); -#else // use builtin - static unsigned short lengthc[] = { 3,4,5,6,7,8,9,10,11,13,15,17,19,23,27,31,35,43,51,59,67,83,99,115,131,163,195,227,258, 259 }; - static unsigned char lengtheb[]= { 0,0,0,0,0,0,0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 0 }; - static unsigned short distc[] = { 1,2,3,4,5,7,9,13,17,25,33,49,65,97,129,193,257,385,513,769,1025,1537,2049,3073,4097,6145,8193,12289,16385,24577, 32768 }; - static unsigned char disteb[] = { 0,0,0,0,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,11,11,12,12,13,13 }; - unsigned int bitbuf=0; - int i,j, bitcount=0; - unsigned char *out = NULL; - unsigned char ***hash_table = (unsigned char***) STBIW_MALLOC(stbiw__ZHASH * sizeof(unsigned char**)); - if (hash_table == NULL) - return NULL; - if (quality < 5) quality = 5; - - stbiw__sbpush(out, 0x78); // DEFLATE 32K window - stbiw__sbpush(out, 0x5e); // FLEVEL = 1 - stbiw__zlib_add(1,1); // BFINAL = 1 - stbiw__zlib_add(1,2); // BTYPE = 1 -- fixed huffman - - for (i=0; i < stbiw__ZHASH; ++i) - hash_table[i] = NULL; - - i=0; - while (i < data_len-3) { - // hash next 3 bytes of data to be compressed - int h = stbiw__zhash(data+i)&(stbiw__ZHASH-1), best=3; - unsigned char *bestloc = 0; - unsigned char **hlist = hash_table[h]; - int n = stbiw__sbcount(hlist); - for (j=0; j < n; ++j) { - if (hlist[j]-data > i-32768) { // if entry lies within window - int d = stbiw__zlib_countm(hlist[j], data+i, data_len-i); - if (d >= best) { best=d; bestloc=hlist[j]; } - } - } - // when hash table entry is too long, delete half the entries - if (hash_table[h] && stbiw__sbn(hash_table[h]) == 2*quality) { - STBIW_MEMMOVE(hash_table[h], hash_table[h]+quality, sizeof(hash_table[h][0])*quality); - stbiw__sbn(hash_table[h]) = quality; - } - stbiw__sbpush(hash_table[h],data+i); - - if (bestloc) { - // "lazy matching" - check match at *next* byte, and if it's better, do cur byte as literal - h = stbiw__zhash(data+i+1)&(stbiw__ZHASH-1); - hlist = hash_table[h]; - n = stbiw__sbcount(hlist); - for (j=0; j < n; ++j) { - if (hlist[j]-data > i-32767) { - int e = stbiw__zlib_countm(hlist[j], data+i+1, data_len-i-1); - if (e > best) { // if next match is better, bail on current match - bestloc = NULL; - break; - } - } - } - } - - if (bestloc) { - int d = (int) (data+i - bestloc); // distance back - STBIW_ASSERT(d <= 32767 && best <= 258); - for (j=0; best > lengthc[j+1]-1; ++j); - stbiw__zlib_huff(j+257); - if (lengtheb[j]) stbiw__zlib_add(best - lengthc[j], lengtheb[j]); - for (j=0; d > distc[j+1]-1; ++j); - stbiw__zlib_add(stbiw__zlib_bitrev(j,5),5); - if (disteb[j]) stbiw__zlib_add(d - distc[j], disteb[j]); - i += best; - } else { - stbiw__zlib_huffb(data[i]); - ++i; - } - } - // write out final bytes - for (;i < data_len; ++i) - stbiw__zlib_huffb(data[i]); - stbiw__zlib_huff(256); // end of block - // pad with 0 bits to byte boundary - while (bitcount) - stbiw__zlib_add(0,1); - - for (i=0; i < stbiw__ZHASH; ++i) - (void) stbiw__sbfree(hash_table[i]); - STBIW_FREE(hash_table); - - // store uncompressed instead if compression was worse - if (stbiw__sbn(out) > data_len + 2 + ((data_len+32766)/32767)*5) { - stbiw__sbn(out) = 2; // truncate to DEFLATE 32K window and FLEVEL = 1 - for (j = 0; j < data_len;) { - int blocklen = data_len - j; - if (blocklen > 32767) blocklen = 32767; - stbiw__sbpush(out, data_len - j == blocklen); // BFINAL = ?, BTYPE = 0 -- no compression - stbiw__sbpush(out, STBIW_UCHAR(blocklen)); // LEN - stbiw__sbpush(out, STBIW_UCHAR(blocklen >> 8)); - stbiw__sbpush(out, STBIW_UCHAR(~blocklen)); // NLEN - stbiw__sbpush(out, STBIW_UCHAR(~blocklen >> 8)); - memcpy(out+stbiw__sbn(out), data+j, blocklen); - stbiw__sbn(out) += blocklen; - j += blocklen; - } - } - - { - // compute adler32 on input - unsigned int s1=1, s2=0; - int blocklen = (int) (data_len % 5552); - j=0; - while (j < data_len) { - for (i=0; i < blocklen; ++i) { s1 += data[j+i]; s2 += s1; } - s1 %= 65521; s2 %= 65521; - j += blocklen; - blocklen = 5552; - } - stbiw__sbpush(out, STBIW_UCHAR(s2 >> 8)); - stbiw__sbpush(out, STBIW_UCHAR(s2)); - stbiw__sbpush(out, STBIW_UCHAR(s1 >> 8)); - stbiw__sbpush(out, STBIW_UCHAR(s1)); - } - *out_len = stbiw__sbn(out); - // make returned pointer freeable - STBIW_MEMMOVE(stbiw__sbraw(out), out, *out_len); - return (unsigned char *) stbiw__sbraw(out); -#endif // STBIW_ZLIB_COMPRESS -} - -static unsigned int stbiw__crc32(unsigned char *buffer, int len) -{ -#ifdef STBIW_CRC32 - return STBIW_CRC32(buffer, len); -#else - static unsigned int crc_table[256] = - { - 0x00000000, 0x77073096, 0xEE0E612C, 0x990951BA, 0x076DC419, 0x706AF48F, 0xE963A535, 0x9E6495A3, - 0x0eDB8832, 0x79DCB8A4, 0xE0D5E91E, 0x97D2D988, 0x09B64C2B, 0x7EB17CBD, 0xE7B82D07, 0x90BF1D91, - 0x1DB71064, 0x6AB020F2, 0xF3B97148, 0x84BE41DE, 0x1ADAD47D, 0x6DDDE4EB, 0xF4D4B551, 0x83D385C7, - 0x136C9856, 0x646BA8C0, 0xFD62F97A, 0x8A65C9EC, 0x14015C4F, 0x63066CD9, 0xFA0F3D63, 0x8D080DF5, - 0x3B6E20C8, 0x4C69105E, 0xD56041E4, 0xA2677172, 0x3C03E4D1, 0x4B04D447, 0xD20D85FD, 0xA50AB56B, - 0x35B5A8FA, 0x42B2986C, 0xDBBBC9D6, 0xACBCF940, 0x32D86CE3, 0x45DF5C75, 0xDCD60DCF, 0xABD13D59, - 0x26D930AC, 0x51DE003A, 0xC8D75180, 0xBFD06116, 0x21B4F4B5, 0x56B3C423, 0xCFBA9599, 0xB8BDA50F, - 0x2802B89E, 0x5F058808, 0xC60CD9B2, 0xB10BE924, 0x2F6F7C87, 0x58684C11, 0xC1611DAB, 0xB6662D3D, - 0x76DC4190, 0x01DB7106, 0x98D220BC, 0xEFD5102A, 0x71B18589, 0x06B6B51F, 0x9FBFE4A5, 0xE8B8D433, - 0x7807C9A2, 0x0F00F934, 0x9609A88E, 0xE10E9818, 0x7F6A0DBB, 0x086D3D2D, 0x91646C97, 0xE6635C01, - 0x6B6B51F4, 0x1C6C6162, 0x856530D8, 0xF262004E, 0x6C0695ED, 0x1B01A57B, 0x8208F4C1, 0xF50FC457, - 0x65B0D9C6, 0x12B7E950, 0x8BBEB8EA, 0xFCB9887C, 0x62DD1DDF, 0x15DA2D49, 0x8CD37CF3, 0xFBD44C65, - 0x4DB26158, 0x3AB551CE, 0xA3BC0074, 0xD4BB30E2, 0x4ADFA541, 0x3DD895D7, 0xA4D1C46D, 0xD3D6F4FB, - 0x4369E96A, 0x346ED9FC, 0xAD678846, 0xDA60B8D0, 0x44042D73, 0x33031DE5, 0xAA0A4C5F, 0xDD0D7CC9, - 0x5005713C, 0x270241AA, 0xBE0B1010, 0xC90C2086, 0x5768B525, 0x206F85B3, 0xB966D409, 0xCE61E49F, - 0x5EDEF90E, 0x29D9C998, 0xB0D09822, 0xC7D7A8B4, 0x59B33D17, 0x2EB40D81, 0xB7BD5C3B, 0xC0BA6CAD, - 0xEDB88320, 0x9ABFB3B6, 0x03B6E20C, 0x74B1D29A, 0xEAD54739, 0x9DD277AF, 0x04DB2615, 0x73DC1683, - 0xE3630B12, 0x94643B84, 0x0D6D6A3E, 0x7A6A5AA8, 0xE40ECF0B, 0x9309FF9D, 0x0A00AE27, 0x7D079EB1, - 0xF00F9344, 0x8708A3D2, 0x1E01F268, 0x6906C2FE, 0xF762575D, 0x806567CB, 0x196C3671, 0x6E6B06E7, - 0xFED41B76, 0x89D32BE0, 0x10DA7A5A, 0x67DD4ACC, 0xF9B9DF6F, 0x8EBEEFF9, 0x17B7BE43, 0x60B08ED5, - 0xD6D6A3E8, 0xA1D1937E, 0x38D8C2C4, 0x4FDFF252, 0xD1BB67F1, 0xA6BC5767, 0x3FB506DD, 0x48B2364B, - 0xD80D2BDA, 0xAF0A1B4C, 0x36034AF6, 0x41047A60, 0xDF60EFC3, 0xA867DF55, 0x316E8EEF, 0x4669BE79, - 0xCB61B38C, 0xBC66831A, 0x256FD2A0, 0x5268E236, 0xCC0C7795, 0xBB0B4703, 0x220216B9, 0x5505262F, - 0xC5BA3BBE, 0xB2BD0B28, 0x2BB45A92, 0x5CB36A04, 0xC2D7FFA7, 0xB5D0CF31, 0x2CD99E8B, 0x5BDEAE1D, - 0x9B64C2B0, 0xEC63F226, 0x756AA39C, 0x026D930A, 0x9C0906A9, 0xEB0E363F, 0x72076785, 0x05005713, - 0x95BF4A82, 0xE2B87A14, 0x7BB12BAE, 0x0CB61B38, 0x92D28E9B, 0xE5D5BE0D, 0x7CDCEFB7, 0x0BDBDF21, - 0x86D3D2D4, 0xF1D4E242, 0x68DDB3F8, 0x1FDA836E, 0x81BE16CD, 0xF6B9265B, 0x6FB077E1, 0x18B74777, - 0x88085AE6, 0xFF0F6A70, 0x66063BCA, 0x11010B5C, 0x8F659EFF, 0xF862AE69, 0x616BFFD3, 0x166CCF45, - 0xA00AE278, 0xD70DD2EE, 0x4E048354, 0x3903B3C2, 0xA7672661, 0xD06016F7, 0x4969474D, 0x3E6E77DB, - 0xAED16A4A, 0xD9D65ADC, 0x40DF0B66, 0x37D83BF0, 0xA9BCAE53, 0xDEBB9EC5, 0x47B2CF7F, 0x30B5FFE9, - 0xBDBDF21C, 0xCABAC28A, 0x53B39330, 0x24B4A3A6, 0xBAD03605, 0xCDD70693, 0x54DE5729, 0x23D967BF, - 0xB3667A2E, 0xC4614AB8, 0x5D681B02, 0x2A6F2B94, 0xB40BBE37, 0xC30C8EA1, 0x5A05DF1B, 0x2D02EF8D - }; - - unsigned int crc = ~0u; - int i; - for (i=0; i < len; ++i) - crc = (crc >> 8) ^ crc_table[buffer[i] ^ (crc & 0xff)]; - return ~crc; -#endif -} - -#define stbiw__wpng4(o,a,b,c,d) ((o)[0]=STBIW_UCHAR(a),(o)[1]=STBIW_UCHAR(b),(o)[2]=STBIW_UCHAR(c),(o)[3]=STBIW_UCHAR(d),(o)+=4) -#define stbiw__wp32(data,v) stbiw__wpng4(data, (v)>>24,(v)>>16,(v)>>8,(v)); -#define stbiw__wptag(data,s) stbiw__wpng4(data, s[0],s[1],s[2],s[3]) - -static void stbiw__wpcrc(unsigned char **data, int len) -{ - unsigned int crc = stbiw__crc32(*data - len - 4, len+4); - stbiw__wp32(*data, crc); -} - -static unsigned char stbiw__paeth(int a, int b, int c) -{ - int p = a + b - c, pa = abs(p-a), pb = abs(p-b), pc = abs(p-c); - if (pa <= pb && pa <= pc) return STBIW_UCHAR(a); - if (pb <= pc) return STBIW_UCHAR(b); - return STBIW_UCHAR(c); -} - -// @OPTIMIZE: provide an option that always forces left-predict or paeth predict -static void stbiw__encode_png_line(unsigned char *pixels, int stride_bytes, int width, int height, int y, int n, int filter_type, signed char *line_buffer) -{ - static int mapping[] = { 0,1,2,3,4 }; - static int firstmap[] = { 0,1,0,5,6 }; - int *mymap = (y != 0) ? mapping : firstmap; - int i; - int type = mymap[filter_type]; - unsigned char *z = pixels + stride_bytes * (stbi__flip_vertically_on_write ? height-1-y : y); - int signed_stride = stbi__flip_vertically_on_write ? -stride_bytes : stride_bytes; - - if (type==0) { - memcpy(line_buffer, z, width*n); - return; - } - - // first loop isn't optimized since it's just one pixel - for (i = 0; i < n; ++i) { - switch (type) { - case 1: line_buffer[i] = z[i]; break; - case 2: line_buffer[i] = z[i] - z[i-signed_stride]; break; - case 3: line_buffer[i] = z[i] - (z[i-signed_stride]>>1); break; - case 4: line_buffer[i] = (signed char) (z[i] - stbiw__paeth(0,z[i-signed_stride],0)); break; - case 5: line_buffer[i] = z[i]; break; - case 6: line_buffer[i] = z[i]; break; - } - } - switch (type) { - case 1: for (i=n; i < width*n; ++i) line_buffer[i] = z[i] - z[i-n]; break; - case 2: for (i=n; i < width*n; ++i) line_buffer[i] = z[i] - z[i-signed_stride]; break; - case 3: for (i=n; i < width*n; ++i) line_buffer[i] = z[i] - ((z[i-n] + z[i-signed_stride])>>1); break; - case 4: for (i=n; i < width*n; ++i) line_buffer[i] = z[i] - stbiw__paeth(z[i-n], z[i-signed_stride], z[i-signed_stride-n]); break; - case 5: for (i=n; i < width*n; ++i) line_buffer[i] = z[i] - (z[i-n]>>1); break; - case 6: for (i=n; i < width*n; ++i) line_buffer[i] = z[i] - stbiw__paeth(z[i-n], 0,0); break; - } -} - -STBIWDEF unsigned char *stbi_write_png_to_mem(const unsigned char *pixels, int stride_bytes, int x, int y, int n, int *out_len) -{ - int force_filter = stbi_write_force_png_filter; - int ctype[5] = { -1, 0, 4, 2, 6 }; - unsigned char sig[8] = { 137,80,78,71,13,10,26,10 }; - unsigned char *out,*o, *filt, *zlib; - signed char *line_buffer; - int j,zlen; - - if (stride_bytes == 0) - stride_bytes = x * n; - - if (force_filter >= 5) { - force_filter = -1; - } - - filt = (unsigned char *) STBIW_MALLOC((x*n+1) * y); if (!filt) return 0; - line_buffer = (signed char *) STBIW_MALLOC(x * n); if (!line_buffer) { STBIW_FREE(filt); return 0; } - for (j=0; j < y; ++j) { - int filter_type; - if (force_filter > -1) { - filter_type = force_filter; - stbiw__encode_png_line((unsigned char*)(pixels), stride_bytes, x, y, j, n, force_filter, line_buffer); - } else { // Estimate the best filter by running through all of them: - int best_filter = 0, best_filter_val = 0x7fffffff, est, i; - for (filter_type = 0; filter_type < 5; filter_type++) { - stbiw__encode_png_line((unsigned char*)(pixels), stride_bytes, x, y, j, n, filter_type, line_buffer); - - // Estimate the entropy of the line using this filter; the less, the better. - est = 0; - for (i = 0; i < x*n; ++i) { - est += abs((signed char) line_buffer[i]); - } - if (est < best_filter_val) { - best_filter_val = est; - best_filter = filter_type; - } - } - if (filter_type != best_filter) { // If the last iteration already got us the best filter, don't redo it - stbiw__encode_png_line((unsigned char*)(pixels), stride_bytes, x, y, j, n, best_filter, line_buffer); - filter_type = best_filter; - } - } - // when we get here, filter_type contains the filter type, and line_buffer contains the data - filt[j*(x*n+1)] = (unsigned char) filter_type; - STBIW_MEMMOVE(filt+j*(x*n+1)+1, line_buffer, x*n); - } - STBIW_FREE(line_buffer); - zlib = stbi_zlib_compress(filt, y*( x*n+1), &zlen, stbi_write_png_compression_level); - STBIW_FREE(filt); - if (!zlib) return 0; - - // each tag requires 12 bytes of overhead - out = (unsigned char *) STBIW_MALLOC(8 + 12+13 + 12+zlen + 12); - if (!out) return 0; - *out_len = 8 + 12+13 + 12+zlen + 12; - - o=out; - STBIW_MEMMOVE(o,sig,8); o+= 8; - stbiw__wp32(o, 13); // header length - stbiw__wptag(o, "IHDR"); - stbiw__wp32(o, x); - stbiw__wp32(o, y); - *o++ = 8; - *o++ = STBIW_UCHAR(ctype[n]); - *o++ = 0; - *o++ = 0; - *o++ = 0; - stbiw__wpcrc(&o,13); - - stbiw__wp32(o, zlen); - stbiw__wptag(o, "IDAT"); - STBIW_MEMMOVE(o, zlib, zlen); - o += zlen; - STBIW_FREE(zlib); - stbiw__wpcrc(&o, zlen); - - stbiw__wp32(o,0); - stbiw__wptag(o, "IEND"); - stbiw__wpcrc(&o,0); - - STBIW_ASSERT(o == out + *out_len); - - return out; -} - -#ifndef STBI_WRITE_NO_STDIO -STBIWDEF int stbi_write_png(char const *filename, int x, int y, int comp, const void *data, int stride_bytes) -{ - FILE *f; - int len; - unsigned char *png = stbi_write_png_to_mem((const unsigned char *) data, stride_bytes, x, y, comp, &len); - if (png == NULL) return 0; - - f = stbiw__fopen(filename, "wb"); - if (!f) { STBIW_FREE(png); return 0; } - fwrite(png, 1, len, f); - fclose(f); - STBIW_FREE(png); - return 1; -} -#endif - -STBIWDEF int stbi_write_png_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int stride_bytes) -{ - int len; - unsigned char *png = stbi_write_png_to_mem((const unsigned char *) data, stride_bytes, x, y, comp, &len); - if (png == NULL) return 0; - func(context, png, len); - STBIW_FREE(png); - return 1; -} - - -/* *************************************************************************** - * - * JPEG writer - * - * This is based on Jon Olick's jo_jpeg.cpp: - * public domain Simple, Minimalistic JPEG writer - http://www.jonolick.com/code.html - */ - -static const unsigned char stbiw__jpg_ZigZag[] = { 0,1,5,6,14,15,27,28,2,4,7,13,16,26,29,42,3,8,12,17,25,30,41,43,9,11,18, - 24,31,40,44,53,10,19,23,32,39,45,52,54,20,22,33,38,46,51,55,60,21,34,37,47,50,56,59,61,35,36,48,49,57,58,62,63 }; - -static void stbiw__jpg_writeBits(stbi__write_context *s, int *bitBufP, int *bitCntP, const unsigned short *bs) { - int bitBuf = *bitBufP, bitCnt = *bitCntP; - bitCnt += bs[1]; - bitBuf |= bs[0] << (24 - bitCnt); - while(bitCnt >= 8) { - unsigned char c = (bitBuf >> 16) & 255; - stbiw__putc(s, c); - if(c == 255) { - stbiw__putc(s, 0); - } - bitBuf <<= 8; - bitCnt -= 8; - } - *bitBufP = bitBuf; - *bitCntP = bitCnt; -} - -static void stbiw__jpg_DCT(float *d0p, float *d1p, float *d2p, float *d3p, float *d4p, float *d5p, float *d6p, float *d7p) { - float d0 = *d0p, d1 = *d1p, d2 = *d2p, d3 = *d3p, d4 = *d4p, d5 = *d5p, d6 = *d6p, d7 = *d7p; - float z1, z2, z3, z4, z5, z11, z13; - - float tmp0 = d0 + d7; - float tmp7 = d0 - d7; - float tmp1 = d1 + d6; - float tmp6 = d1 - d6; - float tmp2 = d2 + d5; - float tmp5 = d2 - d5; - float tmp3 = d3 + d4; - float tmp4 = d3 - d4; - - // Even part - float tmp10 = tmp0 + tmp3; // phase 2 - float tmp13 = tmp0 - tmp3; - float tmp11 = tmp1 + tmp2; - float tmp12 = tmp1 - tmp2; - - d0 = tmp10 + tmp11; // phase 3 - d4 = tmp10 - tmp11; - - z1 = (tmp12 + tmp13) * 0.707106781f; // c4 - d2 = tmp13 + z1; // phase 5 - d6 = tmp13 - z1; - - // Odd part - tmp10 = tmp4 + tmp5; // phase 2 - tmp11 = tmp5 + tmp6; - tmp12 = tmp6 + tmp7; - - // The rotator is modified from fig 4-8 to avoid extra negations. - z5 = (tmp10 - tmp12) * 0.382683433f; // c6 - z2 = tmp10 * 0.541196100f + z5; // c2-c6 - z4 = tmp12 * 1.306562965f + z5; // c2+c6 - z3 = tmp11 * 0.707106781f; // c4 - - z11 = tmp7 + z3; // phase 5 - z13 = tmp7 - z3; - - *d5p = z13 + z2; // phase 6 - *d3p = z13 - z2; - *d1p = z11 + z4; - *d7p = z11 - z4; - - *d0p = d0; *d2p = d2; *d4p = d4; *d6p = d6; -} - -static void stbiw__jpg_calcBits(int val, unsigned short bits[2]) { - int tmp1 = val < 0 ? -val : val; - val = val < 0 ? val-1 : val; - bits[1] = 1; - while(tmp1 >>= 1) { - ++bits[1]; - } - bits[0] = val & ((1<0)&&(DU[end0pos]==0); --end0pos) { - } - // end0pos = first element in reverse order !=0 - if(end0pos == 0) { - stbiw__jpg_writeBits(s, bitBuf, bitCnt, EOB); - return DU[0]; - } - for(i = 1; i <= end0pos; ++i) { - int startpos = i; - int nrzeroes; - unsigned short bits[2]; - for (; DU[i]==0 && i<=end0pos; ++i) { - } - nrzeroes = i-startpos; - if ( nrzeroes >= 16 ) { - int lng = nrzeroes>>4; - int nrmarker; - for (nrmarker=1; nrmarker <= lng; ++nrmarker) - stbiw__jpg_writeBits(s, bitBuf, bitCnt, M16zeroes); - nrzeroes &= 15; - } - stbiw__jpg_calcBits(DU[i], bits); - stbiw__jpg_writeBits(s, bitBuf, bitCnt, HTAC[(nrzeroes<<4)+bits[1]]); - stbiw__jpg_writeBits(s, bitBuf, bitCnt, bits); - } - if(end0pos != 63) { - stbiw__jpg_writeBits(s, bitBuf, bitCnt, EOB); - } - return DU[0]; -} - -static int stbi_write_jpg_core(stbi__write_context *s, int width, int height, int comp, const void* data, int quality) { - // Constants that don't pollute global namespace - static const unsigned char std_dc_luminance_nrcodes[] = {0,0,1,5,1,1,1,1,1,1,0,0,0,0,0,0,0}; - static const unsigned char std_dc_luminance_values[] = {0,1,2,3,4,5,6,7,8,9,10,11}; - static const unsigned char std_ac_luminance_nrcodes[] = {0,0,2,1,3,3,2,4,3,5,5,4,4,0,0,1,0x7d}; - static const unsigned char std_ac_luminance_values[] = { - 0x01,0x02,0x03,0x00,0x04,0x11,0x05,0x12,0x21,0x31,0x41,0x06,0x13,0x51,0x61,0x07,0x22,0x71,0x14,0x32,0x81,0x91,0xa1,0x08, - 0x23,0x42,0xb1,0xc1,0x15,0x52,0xd1,0xf0,0x24,0x33,0x62,0x72,0x82,0x09,0x0a,0x16,0x17,0x18,0x19,0x1a,0x25,0x26,0x27,0x28, - 0x29,0x2a,0x34,0x35,0x36,0x37,0x38,0x39,0x3a,0x43,0x44,0x45,0x46,0x47,0x48,0x49,0x4a,0x53,0x54,0x55,0x56,0x57,0x58,0x59, - 0x5a,0x63,0x64,0x65,0x66,0x67,0x68,0x69,0x6a,0x73,0x74,0x75,0x76,0x77,0x78,0x79,0x7a,0x83,0x84,0x85,0x86,0x87,0x88,0x89, - 0x8a,0x92,0x93,0x94,0x95,0x96,0x97,0x98,0x99,0x9a,0xa2,0xa3,0xa4,0xa5,0xa6,0xa7,0xa8,0xa9,0xaa,0xb2,0xb3,0xb4,0xb5,0xb6, - 0xb7,0xb8,0xb9,0xba,0xc2,0xc3,0xc4,0xc5,0xc6,0xc7,0xc8,0xc9,0xca,0xd2,0xd3,0xd4,0xd5,0xd6,0xd7,0xd8,0xd9,0xda,0xe1,0xe2, - 0xe3,0xe4,0xe5,0xe6,0xe7,0xe8,0xe9,0xea,0xf1,0xf2,0xf3,0xf4,0xf5,0xf6,0xf7,0xf8,0xf9,0xfa - }; - static const unsigned char std_dc_chrominance_nrcodes[] = {0,0,3,1,1,1,1,1,1,1,1,1,0,0,0,0,0}; - static const unsigned char std_dc_chrominance_values[] = {0,1,2,3,4,5,6,7,8,9,10,11}; - static const unsigned char std_ac_chrominance_nrcodes[] = {0,0,2,1,2,4,4,3,4,7,5,4,4,0,1,2,0x77}; - static const unsigned char std_ac_chrominance_values[] = { - 0x00,0x01,0x02,0x03,0x11,0x04,0x05,0x21,0x31,0x06,0x12,0x41,0x51,0x07,0x61,0x71,0x13,0x22,0x32,0x81,0x08,0x14,0x42,0x91, - 0xa1,0xb1,0xc1,0x09,0x23,0x33,0x52,0xf0,0x15,0x62,0x72,0xd1,0x0a,0x16,0x24,0x34,0xe1,0x25,0xf1,0x17,0x18,0x19,0x1a,0x26, - 0x27,0x28,0x29,0x2a,0x35,0x36,0x37,0x38,0x39,0x3a,0x43,0x44,0x45,0x46,0x47,0x48,0x49,0x4a,0x53,0x54,0x55,0x56,0x57,0x58, - 0x59,0x5a,0x63,0x64,0x65,0x66,0x67,0x68,0x69,0x6a,0x73,0x74,0x75,0x76,0x77,0x78,0x79,0x7a,0x82,0x83,0x84,0x85,0x86,0x87, - 0x88,0x89,0x8a,0x92,0x93,0x94,0x95,0x96,0x97,0x98,0x99,0x9a,0xa2,0xa3,0xa4,0xa5,0xa6,0xa7,0xa8,0xa9,0xaa,0xb2,0xb3,0xb4, - 0xb5,0xb6,0xb7,0xb8,0xb9,0xba,0xc2,0xc3,0xc4,0xc5,0xc6,0xc7,0xc8,0xc9,0xca,0xd2,0xd3,0xd4,0xd5,0xd6,0xd7,0xd8,0xd9,0xda, - 0xe2,0xe3,0xe4,0xe5,0xe6,0xe7,0xe8,0xe9,0xea,0xf2,0xf3,0xf4,0xf5,0xf6,0xf7,0xf8,0xf9,0xfa - }; - // Huffman tables - static const unsigned short YDC_HT[256][2] = { {0,2},{2,3},{3,3},{4,3},{5,3},{6,3},{14,4},{30,5},{62,6},{126,7},{254,8},{510,9}}; - static const unsigned short UVDC_HT[256][2] = { {0,2},{1,2},{2,2},{6,3},{14,4},{30,5},{62,6},{126,7},{254,8},{510,9},{1022,10},{2046,11}}; - static const unsigned short YAC_HT[256][2] = { - {10,4},{0,2},{1,2},{4,3},{11,4},{26,5},{120,7},{248,8},{1014,10},{65410,16},{65411,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {12,4},{27,5},{121,7},{502,9},{2038,11},{65412,16},{65413,16},{65414,16},{65415,16},{65416,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {28,5},{249,8},{1015,10},{4084,12},{65417,16},{65418,16},{65419,16},{65420,16},{65421,16},{65422,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {58,6},{503,9},{4085,12},{65423,16},{65424,16},{65425,16},{65426,16},{65427,16},{65428,16},{65429,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {59,6},{1016,10},{65430,16},{65431,16},{65432,16},{65433,16},{65434,16},{65435,16},{65436,16},{65437,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {122,7},{2039,11},{65438,16},{65439,16},{65440,16},{65441,16},{65442,16},{65443,16},{65444,16},{65445,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {123,7},{4086,12},{65446,16},{65447,16},{65448,16},{65449,16},{65450,16},{65451,16},{65452,16},{65453,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {250,8},{4087,12},{65454,16},{65455,16},{65456,16},{65457,16},{65458,16},{65459,16},{65460,16},{65461,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {504,9},{32704,15},{65462,16},{65463,16},{65464,16},{65465,16},{65466,16},{65467,16},{65468,16},{65469,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {505,9},{65470,16},{65471,16},{65472,16},{65473,16},{65474,16},{65475,16},{65476,16},{65477,16},{65478,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {506,9},{65479,16},{65480,16},{65481,16},{65482,16},{65483,16},{65484,16},{65485,16},{65486,16},{65487,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {1017,10},{65488,16},{65489,16},{65490,16},{65491,16},{65492,16},{65493,16},{65494,16},{65495,16},{65496,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {1018,10},{65497,16},{65498,16},{65499,16},{65500,16},{65501,16},{65502,16},{65503,16},{65504,16},{65505,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {2040,11},{65506,16},{65507,16},{65508,16},{65509,16},{65510,16},{65511,16},{65512,16},{65513,16},{65514,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {65515,16},{65516,16},{65517,16},{65518,16},{65519,16},{65520,16},{65521,16},{65522,16},{65523,16},{65524,16},{0,0},{0,0},{0,0},{0,0},{0,0}, - {2041,11},{65525,16},{65526,16},{65527,16},{65528,16},{65529,16},{65530,16},{65531,16},{65532,16},{65533,16},{65534,16},{0,0},{0,0},{0,0},{0,0},{0,0} - }; - static const unsigned short UVAC_HT[256][2] = { - {0,2},{1,2},{4,3},{10,4},{24,5},{25,5},{56,6},{120,7},{500,9},{1014,10},{4084,12},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {11,4},{57,6},{246,8},{501,9},{2038,11},{4085,12},{65416,16},{65417,16},{65418,16},{65419,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {26,5},{247,8},{1015,10},{4086,12},{32706,15},{65420,16},{65421,16},{65422,16},{65423,16},{65424,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {27,5},{248,8},{1016,10},{4087,12},{65425,16},{65426,16},{65427,16},{65428,16},{65429,16},{65430,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {58,6},{502,9},{65431,16},{65432,16},{65433,16},{65434,16},{65435,16},{65436,16},{65437,16},{65438,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {59,6},{1017,10},{65439,16},{65440,16},{65441,16},{65442,16},{65443,16},{65444,16},{65445,16},{65446,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {121,7},{2039,11},{65447,16},{65448,16},{65449,16},{65450,16},{65451,16},{65452,16},{65453,16},{65454,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {122,7},{2040,11},{65455,16},{65456,16},{65457,16},{65458,16},{65459,16},{65460,16},{65461,16},{65462,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {249,8},{65463,16},{65464,16},{65465,16},{65466,16},{65467,16},{65468,16},{65469,16},{65470,16},{65471,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {503,9},{65472,16},{65473,16},{65474,16},{65475,16},{65476,16},{65477,16},{65478,16},{65479,16},{65480,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {504,9},{65481,16},{65482,16},{65483,16},{65484,16},{65485,16},{65486,16},{65487,16},{65488,16},{65489,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {505,9},{65490,16},{65491,16},{65492,16},{65493,16},{65494,16},{65495,16},{65496,16},{65497,16},{65498,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {506,9},{65499,16},{65500,16},{65501,16},{65502,16},{65503,16},{65504,16},{65505,16},{65506,16},{65507,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {2041,11},{65508,16},{65509,16},{65510,16},{65511,16},{65512,16},{65513,16},{65514,16},{65515,16},{65516,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {16352,14},{65517,16},{65518,16},{65519,16},{65520,16},{65521,16},{65522,16},{65523,16},{65524,16},{65525,16},{0,0},{0,0},{0,0},{0,0},{0,0}, - {1018,10},{32707,15},{65526,16},{65527,16},{65528,16},{65529,16},{65530,16},{65531,16},{65532,16},{65533,16},{65534,16},{0,0},{0,0},{0,0},{0,0},{0,0} - }; - static const int YQT[] = {16,11,10,16,24,40,51,61,12,12,14,19,26,58,60,55,14,13,16,24,40,57,69,56,14,17,22,29,51,87,80,62,18,22, - 37,56,68,109,103,77,24,35,55,64,81,104,113,92,49,64,78,87,103,121,120,101,72,92,95,98,112,100,103,99}; - static const int UVQT[] = {17,18,24,47,99,99,99,99,18,21,26,66,99,99,99,99,24,26,56,99,99,99,99,99,47,66,99,99,99,99,99,99, - 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99}; - static const float aasf[] = { 1.0f * 2.828427125f, 1.387039845f * 2.828427125f, 1.306562965f * 2.828427125f, 1.175875602f * 2.828427125f, - 1.0f * 2.828427125f, 0.785694958f * 2.828427125f, 0.541196100f * 2.828427125f, 0.275899379f * 2.828427125f }; - - int row, col, i, k, subsample; - float fdtbl_Y[64], fdtbl_UV[64]; - unsigned char YTable[64], UVTable[64]; - - if(!data || !width || !height || comp > 4 || comp < 1) { - return 0; - } - - quality = quality ? quality : 90; - subsample = quality <= 90 ? 1 : 0; - quality = quality < 1 ? 1 : quality > 100 ? 100 : quality; - quality = quality < 50 ? 5000 / quality : 200 - quality * 2; - - for(i = 0; i < 64; ++i) { - int uvti, yti = (YQT[i]*quality+50)/100; - YTable[stbiw__jpg_ZigZag[i]] = (unsigned char) (yti < 1 ? 1 : yti > 255 ? 255 : yti); - uvti = (UVQT[i]*quality+50)/100; - UVTable[stbiw__jpg_ZigZag[i]] = (unsigned char) (uvti < 1 ? 1 : uvti > 255 ? 255 : uvti); - } - - for(row = 0, k = 0; row < 8; ++row) { - for(col = 0; col < 8; ++col, ++k) { - fdtbl_Y[k] = 1 / (YTable [stbiw__jpg_ZigZag[k]] * aasf[row] * aasf[col]); - fdtbl_UV[k] = 1 / (UVTable[stbiw__jpg_ZigZag[k]] * aasf[row] * aasf[col]); - } - } - - // Write Headers - { - static const unsigned char head0[] = { 0xFF,0xD8,0xFF,0xE0,0,0x10,'J','F','I','F',0,1,1,0,0,1,0,1,0,0,0xFF,0xDB,0,0x84,0 }; - static const unsigned char head2[] = { 0xFF,0xDA,0,0xC,3,1,0,2,0x11,3,0x11,0,0x3F,0 }; - const unsigned char head1[] = { 0xFF,0xC0,0,0x11,8,(unsigned char)(height>>8),STBIW_UCHAR(height),(unsigned char)(width>>8),STBIW_UCHAR(width), - 3,1,(unsigned char)(subsample?0x22:0x11),0,2,0x11,1,3,0x11,1,0xFF,0xC4,0x01,0xA2,0 }; - s->func(s->context, (void*)head0, sizeof(head0)); - s->func(s->context, (void*)YTable, sizeof(YTable)); - stbiw__putc(s, 1); - s->func(s->context, UVTable, sizeof(UVTable)); - s->func(s->context, (void*)head1, sizeof(head1)); - s->func(s->context, (void*)(std_dc_luminance_nrcodes+1), sizeof(std_dc_luminance_nrcodes)-1); - s->func(s->context, (void*)std_dc_luminance_values, sizeof(std_dc_luminance_values)); - stbiw__putc(s, 0x10); // HTYACinfo - s->func(s->context, (void*)(std_ac_luminance_nrcodes+1), sizeof(std_ac_luminance_nrcodes)-1); - s->func(s->context, (void*)std_ac_luminance_values, sizeof(std_ac_luminance_values)); - stbiw__putc(s, 1); // HTUDCinfo - s->func(s->context, (void*)(std_dc_chrominance_nrcodes+1), sizeof(std_dc_chrominance_nrcodes)-1); - s->func(s->context, (void*)std_dc_chrominance_values, sizeof(std_dc_chrominance_values)); - stbiw__putc(s, 0x11); // HTUACinfo - s->func(s->context, (void*)(std_ac_chrominance_nrcodes+1), sizeof(std_ac_chrominance_nrcodes)-1); - s->func(s->context, (void*)std_ac_chrominance_values, sizeof(std_ac_chrominance_values)); - s->func(s->context, (void*)head2, sizeof(head2)); - } - - // Encode 8x8 macroblocks - { - static const unsigned short fillBits[] = {0x7F, 7}; - int DCY=0, DCU=0, DCV=0; - int bitBuf=0, bitCnt=0; - // comp == 2 is grey+alpha (alpha is ignored) - int ofsG = comp > 2 ? 1 : 0, ofsB = comp > 2 ? 2 : 0; - const unsigned char *dataR = (const unsigned char *)data; - const unsigned char *dataG = dataR + ofsG; - const unsigned char *dataB = dataR + ofsB; - int x, y, pos; - if(subsample) { - for(y = 0; y < height; y += 16) { - for(x = 0; x < width; x += 16) { - float Y[256], U[256], V[256]; - for(row = y, pos = 0; row < y+16; ++row) { - // row >= height => use last input row - int clamped_row = (row < height) ? row : height - 1; - int base_p = (stbi__flip_vertically_on_write ? (height-1-clamped_row) : clamped_row)*width*comp; - for(col = x; col < x+16; ++col, ++pos) { - // if col >= width => use pixel from last input column - int p = base_p + ((col < width) ? col : (width-1))*comp; - float r = dataR[p], g = dataG[p], b = dataB[p]; - Y[pos]= +0.29900f*r + 0.58700f*g + 0.11400f*b - 128; - U[pos]= -0.16874f*r - 0.33126f*g + 0.50000f*b; - V[pos]= +0.50000f*r - 0.41869f*g - 0.08131f*b; - } - } - DCY = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, Y+0, 16, fdtbl_Y, DCY, YDC_HT, YAC_HT); - DCY = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, Y+8, 16, fdtbl_Y, DCY, YDC_HT, YAC_HT); - DCY = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, Y+128, 16, fdtbl_Y, DCY, YDC_HT, YAC_HT); - DCY = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, Y+136, 16, fdtbl_Y, DCY, YDC_HT, YAC_HT); - - // subsample U,V - { - float subU[64], subV[64]; - int yy, xx; - for(yy = 0, pos = 0; yy < 8; ++yy) { - for(xx = 0; xx < 8; ++xx, ++pos) { - int j = yy*32+xx*2; - subU[pos] = (U[j+0] + U[j+1] + U[j+16] + U[j+17]) * 0.25f; - subV[pos] = (V[j+0] + V[j+1] + V[j+16] + V[j+17]) * 0.25f; - } - } - DCU = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, subU, 8, fdtbl_UV, DCU, UVDC_HT, UVAC_HT); - DCV = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, subV, 8, fdtbl_UV, DCV, UVDC_HT, UVAC_HT); - } - } - } - } else { - for(y = 0; y < height; y += 8) { - for(x = 0; x < width; x += 8) { - float Y[64], U[64], V[64]; - for(row = y, pos = 0; row < y+8; ++row) { - // row >= height => use last input row - int clamped_row = (row < height) ? row : height - 1; - int base_p = (stbi__flip_vertically_on_write ? (height-1-clamped_row) : clamped_row)*width*comp; - for(col = x; col < x+8; ++col, ++pos) { - // if col >= width => use pixel from last input column - int p = base_p + ((col < width) ? col : (width-1))*comp; - float r = dataR[p], g = dataG[p], b = dataB[p]; - Y[pos]= +0.29900f*r + 0.58700f*g + 0.11400f*b - 128; - U[pos]= -0.16874f*r - 0.33126f*g + 0.50000f*b; - V[pos]= +0.50000f*r - 0.41869f*g - 0.08131f*b; - } - } - - DCY = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, Y, 8, fdtbl_Y, DCY, YDC_HT, YAC_HT); - DCU = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, U, 8, fdtbl_UV, DCU, UVDC_HT, UVAC_HT); - DCV = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, V, 8, fdtbl_UV, DCV, UVDC_HT, UVAC_HT); - } - } - } - - // Do the bit alignment of the EOI marker - stbiw__jpg_writeBits(s, &bitBuf, &bitCnt, fillBits); - } - - // EOI - stbiw__putc(s, 0xFF); - stbiw__putc(s, 0xD9); - - return 1; -} - -STBIWDEF int stbi_write_jpg_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int quality) -{ - stbi__write_context s = { 0 }; - stbi__start_write_callbacks(&s, func, context); - return stbi_write_jpg_core(&s, x, y, comp, (void *) data, quality); -} - - -#ifndef STBI_WRITE_NO_STDIO -STBIWDEF int stbi_write_jpg(char const *filename, int x, int y, int comp, const void *data, int quality) -{ - stbi__write_context s = { 0 }; - if (stbi__start_write_file(&s,filename)) { - int r = stbi_write_jpg_core(&s, x, y, comp, data, quality); - stbi__end_write_file(&s); - return r; - } else - return 0; -} -#endif - -#endif // STB_IMAGE_WRITE_IMPLEMENTATION - -/* Revision history - 1.16 (2021-07-11) - make Deflate code emit uncompressed blocks when it would otherwise expand - support writing BMPs with alpha channel - 1.15 (2020-07-13) unknown - 1.14 (2020-02-02) updated JPEG writer to downsample chroma channels - 1.13 - 1.12 - 1.11 (2019-08-11) - - 1.10 (2019-02-07) - support utf8 filenames in Windows; fix warnings and platform ifdefs - 1.09 (2018-02-11) - fix typo in zlib quality API, improve STB_I_W_STATIC in C++ - 1.08 (2018-01-29) - add stbi__flip_vertically_on_write, external zlib, zlib quality, choose PNG filter - 1.07 (2017-07-24) - doc fix - 1.06 (2017-07-23) - writing JPEG (using Jon Olick's code) - 1.05 ??? - 1.04 (2017-03-03) - monochrome BMP expansion - 1.03 ??? - 1.02 (2016-04-02) - avoid allocating large structures on the stack - 1.01 (2016-01-16) - STBIW_REALLOC_SIZED: support allocators with no realloc support - avoid race-condition in crc initialization - minor compile issues - 1.00 (2015-09-14) - installable file IO function - 0.99 (2015-09-13) - warning fixes; TGA rle support - 0.98 (2015-04-08) - added STBIW_MALLOC, STBIW_ASSERT etc - 0.97 (2015-01-18) - fixed HDR asserts, rewrote HDR rle logic - 0.96 (2015-01-17) - add HDR output - fix monochrome BMP - 0.95 (2014-08-17) - add monochrome TGA output - 0.94 (2014-05-31) - rename private functions to avoid conflicts with stb_image.h - 0.93 (2014-05-27) - warning fixes - 0.92 (2010-08-01) - casts to unsigned char to fix warnings - 0.91 (2010-07-17) - first public release - 0.90 first internal release -*/ - -/* ------------------------------------------------------------------------------- -This software is available under 2 licenses -- choose whichever you prefer. ------------------------------------------------------------------------------- -ALTERNATIVE A - MIT License -Copyright (c) 2017 Sean Barrett -Permission is hereby granted, free of charge, to any person obtaining a copy of -this software and associated documentation files (the "Software"), to deal in -the Software without restriction, including without limitation the rights to -use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies -of the Software, and to permit persons to whom the Software is furnished to do -so, subject to the following conditions: -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. ------------------------------------------------------------------------------- -ALTERNATIVE B - Public Domain (www.unlicense.org) -This is free and unencumbered software released into the public domain. -Anyone is free to copy, modify, publish, use, compile, sell, or distribute this -software, either in source code form or as a compiled binary, for any purpose, -commercial or non-commercial, and by any means. -In jurisdictions that recognize copyright laws, the author or authors of this -software dedicate any and all copyright interest in the software to the public -domain. We make this dedication for the benefit of the public at large and to -the detriment of our heirs and successors. We intend this dedication to be an -overt act of relinquishment in perpetuity of all present and future rights to -this software under copyright law. -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN -ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION -WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ------------------------------------------------------------------------------- -*/ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8727680d1..c557a1500 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -101,6 +101,15 @@ if(IPPL_ENABLE_SOLVERS) add_subdirectory(LinearSolvers) endif() + +message (STATUS "Adding Stream subdirectory") +add_subdirectory (Stream) + +if(IPPL_ENABLE_CATALYST) + include(${PROJECT_SOURCE_DIR}/cmake/SetupCatalyst.cmake) +endif() + + include(${PROJECT_SOURCE_DIR}/cmake/PlatformOptions.cmake) target_link_libraries(ippl PUBLIC Kokkos::kokkos MPI::MPI_CXX) diff --git a/src/Particle/ParticleAttrib.h b/src/Particle/ParticleAttrib.h index 80b9dad20..6d8999008 100644 --- a/src/Particle/ParticleAttrib.h +++ b/src/Particle/ParticleAttrib.h @@ -271,6 +271,17 @@ namespace ippl { */ void internalCopy(const hash_type& indices) override; + #ifdef IPPL_ENABLE_CATALYST + void signConduitBlueprintNode( + const size_type Np_local + , conduit_cpp::Node& node_fields + , ViewRegistry& viewRegistry + , Inform& ca_m + , Inform& ca_warn + , const bool forceHostCopy + ) const override ; + #endif + private: view_type dview_m{"ParticleAttrib::dview", 0}; view_type buf_m{"ParticleAttrib::buf", 0}; diff --git a/src/Particle/ParticleAttrib.hpp b/src/Particle/ParticleAttrib.hpp index ff163818f..a0967f282 100644 --- a/src/Particle/ParticleAttrib.hpp +++ b/src/Particle/ParticleAttrib.hpp @@ -531,4 +531,85 @@ namespace ippl { DefineParticleReduction(Min, min, if (myVal < valL) valL = myVal, std::less) DefineParticleReduction(Prod, prod, valL *= myVal, std::multiplies) + + #ifdef IPPL_ENABLE_CATALYST + ////////////////////////////////////////////////////////////////////////////////////// + // Note: + // In general, for runtime performance, neither function overloading nor if + // constexpr has an inherent advantage when used correctly for compile-time + // dispatch. . + // Function overloading with template parameter extraction or constraints or sfinae + // are all hard to implement in this case, so we switched to if const expr. + ////////////////////////////////////////////////////////////////////////////////////// + template + void ParticleAttrib::signConduitBlueprintNode( + const size_type Np_local, + conduit_cpp::Node& node_fields, + ViewRegistry& viewRegistry, + Inform& ca_m, + Inform& ca_warn, + const bool forceHostCopy + ) const + { + host_mirror_type hostMirror; + if(forceHostCopy){ + hostMirror = this->getHostMirror(); + Kokkos::deep_copy(hostMirror , this->getView()); + } else{ + hostMirror = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), this->getView()); + } + auto field = node_fields[this->name_m]; + field["association"].set_string("vertex"); + field["topology"].set_string("p_unstructured_topo"); + field["volume_dependent"].set_string("false"); + + + if constexpr (std::is_scalar_v) { + // --- SCALAR CASE --- + ca_m << level4 <<"::Execute()excute_entry() for attribute: "<name_m << '\n' + << " call to:\n" + << " ParticleAttribute<" << typeid(T).name() << ">::signConduitBlueprintNode()" << endl; + + field["values"].set_external(hostMirror.data(), Np_local); + + + } else if constexpr (is_vector_v) { + // --- VECTOR CASE --- + ca_m << level4 <<"::Execute()excute_entry() for attribute: "<name_m << '\n' + << " call to:\n" + << " ParticleAttribute>::signConduitBlueprintNode()" << endl; + + + using elem_t = std::remove_pointer_t; + const size_t stride_bytes = sizeof(elem_t); + // static constexpr size_t stride_bytes = sizeof(elem_t); + + if(Np_local>0){ + field["values/x"].set_external(&hostMirror.data()[0][0], Np_local, 0 , stride_bytes ); + if constexpr (T::dim>=2){ + field["values/y"].set_external(&hostMirror.data()[0][1], Np_local, 0 , stride_bytes ); + } + if constexpr (T::dim>=3) { + field["values/z"].set_external(&hostMirror.data()[0][2], Np_local, 0 , stride_bytes ); + } + }else /* (Np_local=0) */ { + // If Np_local is 0. We MUST provide valid, empty arrays for the gather to work. + using component_type = typename T::value_type; + field["values/x"].set_external(static_cast(nullptr), 0); + if constexpr (T::dim>=2) field["values/y"].set_external(static_cast(nullptr), 0); + if constexpr (T::dim>=3) field["values/z"].set_external(static_cast(nullptr), 0); + } + } else { + // --- INVALID CASE --- + ca_warn << "::Execute()excute_entry() for attribute:"<name_m << endl + << " call to:" << endl + << " ParticleAttribute<" << typeid(T).name() << ">::signConduitBlueprintNode()" << endl + << " For this type of Attribute the Conduit Blueprint description wasnt \n" + << " implemented in ippl. Therefore this type of attribute is not \n" + << " supported for visualisation." << endl; + } + viewRegistry.set(hostMirror); + } + #endif + } // namespace ippl diff --git a/src/Particle/ParticleAttribBase.h b/src/Particle/ParticleAttribBase.h index 3e9baa197..761fe69a2 100644 --- a/src/Particle/ParticleAttribBase.h +++ b/src/Particle/ParticleAttribBase.h @@ -21,6 +21,12 @@ #include "Communicate/Archive.h" + + #ifdef IPPL_ENABLE_CATALYST + #include + #include "Stream/Registry/ViewRegistry.h" + #endif + namespace ippl { namespace detail { // Maximum length for attribute names (including null terminator) @@ -94,6 +100,16 @@ namespace ippl { virtual void applyPermutation(const hash_type&) = 0; virtual void internalCopy(const hash_type&) = 0; + #ifdef IPPL_ENABLE_CATALYST + virtual void signConduitBlueprintNode( + const size_type Np_local + , conduit_cpp::Node& node_fields + , ViewRegistry& viewRegistry + , Inform& ca_m + , Inform& ca_warn + , const bool forceHostCopy + ) const = 0; + #endif protected: const size_type* localNum_mp; char name_m[ATTRIB_NAME_MAX_LEN]; diff --git a/src/Particle/ParticleBase.h b/src/Particle/ParticleBase.h index 9793da2a9..8dd3395a5 100644 --- a/src/Particle/ParticleBase.h +++ b/src/Particle/ParticleBase.h @@ -85,9 +85,9 @@ namespace ippl { */ template class ParticleBase : public ParticleBaseBase { + public: constexpr static bool EnableIDs = sizeof...(IDProperties) > 0; - public: using vector_type = typename PLayout::vector_type; using index_type = typename PLayout::index_type; using particle_position_type = typename PLayout::particle_position_type; @@ -200,6 +200,12 @@ namespace ippl { return attributes_m.template get()[i]; } + /*! Const overload — needed when called on a const ParticleBase reference. */ + template + const attribute_type* getAttribute(size_t i) const { + return attributes_m.template get()[i]; + } + /*! * Calls a given function for all attributes in the bunch * @tparam MemorySpace the memory space of the attributes to visit (void to visit all of @@ -245,6 +251,28 @@ namespace ippl { return total; } + /** + * @brief Return whether an attribute is one of ParticleBase's built-in attributes. + * + * Attribute lists are partitioned by memory space, so their iteration order cannot be + * used to identify the built-in position and ID attributes reliably. + */ + template + bool isBuiltinAttribute(const detail::ParticleAttribBase* attribute) const { + if (attribute == nullptr) { + return false; + } + + const void* candidate = dynamic_cast(attribute); + if (candidate == static_cast(&R)) { + return true; + } + if constexpr (EnableIDs) { + return candidate == static_cast(&ID); + } + return false; + } + /*! * Create nLocal rank local particles. This is a collective call, * i.e. all MPI ranks must call this. diff --git a/src/Stream/CMakeLists.txt b/src/Stream/CMakeLists.txt new file mode 100644 index 000000000..7d5cc1c07 --- /dev/null +++ b/src/Stream/CMakeLists.txt @@ -0,0 +1,5 @@ +# ----------------------------------------------------------------------------- +# src/Stream/CMakeLists.txt +# ----------------------------------------------------------------------------- +add_subdirectory(InSitu) +add_subdirectory(Registry) \ No newline at end of file diff --git a/src/Stream/InSitu/CMakeLists.txt b/src/Stream/InSitu/CMakeLists.txt new file mode 100644 index 000000000..8e62bae37 --- /dev/null +++ b/src/Stream/InSitu/CMakeLists.txt @@ -0,0 +1,99 @@ +# ----------------------------------------------------------------------------- +# src/Stream/InSitu/CMakeLists.txt +# ----------------------------------------------------------------------------- +target_include_directories(ippl + PUBLIC + $ + $ +) + +if(IPPL_ENABLE_CATALYST) + # ProxyWriter uses Catalyst's Conduit API and is only needed by the Catalyst adaptor. + target_sources(ippl + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/ProxyWriter.cpp + ) + + set(_ippl_catalyst_resource_source_dir + "${CMAKE_CURRENT_SOURCE_DIR}/catalyst_scripts") + set(_ippl_catalyst_resource_build_dir + "${PROJECT_BINARY_DIR}/${CMAKE_INSTALL_DATADIR}/ippl/catalyst_scripts") + + # CONFIGURE_DEPENDS notices added or removed resources. Each generated rule + # also depends on its source file, so editing a script resynchronizes it on + # the next normal build without copying Python cache files. + file( + GLOB_RECURSE _ippl_catalyst_resources + CONFIGURE_DEPENDS + LIST_DIRECTORIES FALSE + RELATIVE "${_ippl_catalyst_resource_source_dir}" + "${_ippl_catalyst_resource_source_dir}/*.py" + "${_ippl_catalyst_resource_source_dir}/*.yaml" + "${_ippl_catalyst_resource_source_dir}/__init__py" + ) + + set(_ippl_catalyst_staged_resources) + foreach(_ippl_catalyst_resource IN LISTS _ippl_catalyst_resources) + set(_ippl_catalyst_resource_source + "${_ippl_catalyst_resource_source_dir}/${_ippl_catalyst_resource}") + set(_ippl_catalyst_resource_output + "${_ippl_catalyst_resource_build_dir}/${_ippl_catalyst_resource}") + get_filename_component( + _ippl_catalyst_resource_output_dir + "${_ippl_catalyst_resource_output}" + DIRECTORY + ) + + add_custom_command( + OUTPUT "${_ippl_catalyst_resource_output}" + COMMAND "${CMAKE_COMMAND}" -E make_directory + "${_ippl_catalyst_resource_output_dir}" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "${_ippl_catalyst_resource_source}" + "${_ippl_catalyst_resource_output}" + DEPENDS "${_ippl_catalyst_resource_source}" + VERBATIM + ) + list(APPEND _ippl_catalyst_staged_resources + "${_ippl_catalyst_resource_output}") + endforeach() + + add_custom_target( + ippl_catalyst_resources ALL + DEPENDS ${_ippl_catalyst_staged_resources} + COMMENT "Staging Catalyst Python and YAML resources" + ) + add_dependencies(ippl ippl_catalyst_resources) + + unset(_ippl_catalyst_resource_source_dir) + unset(_ippl_catalyst_resource_build_dir) + unset(_ippl_catalyst_resources) + unset(_ippl_catalyst_staged_resources) + unset(_ippl_catalyst_resource) + unset(_ippl_catalyst_resource_source) + unset(_ippl_catalyst_resource_output) + unset(_ippl_catalyst_resource_output_dir) +endif() + + +install(FILES + ProxyWriter.h + ProxyWriter.hpp + CatalystAdaptor.h + CatalystAdaptorSteering.hpp + CatalystAdaptor.hpp + CatalystVisitors.h + DESTINATION include/Stream/InSitu +) + +# Catalyst Python pipelines and their YAML configuration are runtime resources, +# not Python packages. Keep them under the platform-independent data directory. +install( + DIRECTORY catalyst_scripts/ + DESTINATION ${CMAKE_INSTALL_DATADIR}/ippl/catalyst_scripts + FILES_MATCHING + PATTERN "__pycache__" EXCLUDE + PATTERN "*.py" + PATTERN "*.yaml" + PATTERN "__init__py" +) diff --git a/src/Stream/InSitu/CatalystAdaptor.h b/src/Stream/InSitu/CatalystAdaptor.h new file mode 100644 index 000000000..ee9e93bc0 --- /dev/null +++ b/src/Stream/InSitu/CatalystAdaptor.h @@ -0,0 +1,727 @@ +/** + * @file CatalystAdaptor.h + * @brief Declarations and lightweight types for ParaView Catalyst in-situ integration. + * + * This header provides helper types, forward declarations, and includes used by the + * Catalyst adaptor. Heavy implementation details live in `CatalystAdaptor.hpp` and + * `CatalystAdaptorSteering.hpp`. + */ +#ifndef CatalystAdaptor_h +#define CatalystAdaptor_h + +#include "Ippl.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#if defined(MPI_VERSION) +#include +#endif + +#include "Utility/IpplException.h" + +#include "Stream/Registry/ViewRegistry.h" + +#include "Stream/InSitu/ProxyWriter.h" +#include "Stream/Registry/RegistryHelper.h" + + +namespace ippl{ + +/* FORWARD DECLARATION */ +class VisRegistryRuntime; + + + +/** + * @struct Button + * @brief Momentary-action control for steering. + * + * A Button acts like an edge-trigger: it reports a single transition when + * pressed and then resets. Useful for one-shot actions triggered from the + * Catalyst GUI without latching state. + */ +struct Button { + Button() = default; + + // explicit not explicit allows: if(my_btn) to work + operator bool() const { return value_m; } + + // explicit Button(bool v) : value_m(v) {} + explicit Button(bool v){ + + if(v){ // Button is initialized pressd + value_m = true; + priorState_m=false; + } + else /* if(!v) */{ // Button is initialized unpressd + value_m=false; + priorState_m = false; + } + + } + + + // Assignment operator: Button = bool + Button& operator=(bool v) { + if(v) { // Button is being pressed + if(!value_m && !priorState_m) { // True unpressed state + value_m = true; + priorState_m = false; // Fixed typo: was "priot_state" + } + else if(value_m && !priorState_m) { // Button was pressed last iteration and is still "being pushed down" + // Internal button needs to snap back! And mark it as "freshly snapped" + value_m = false; + priorState_m = true; + } + else if(!value_m && priorState_m) { // Button was pressed at some previous iteration and is still "being pushed down" + // Button needs to stay in freshly snapped back state + value_m = false; + priorState_m = true; + } + else if(value_m && priorState_m) { // Impossible state + throw IpplException("CatalystAdaptor::Button Assignment", "Impossible State: Button is malfunctioning!!!"); + } + } + else { // Button is unpressed + value_m = false; + priorState_m = false; + } + return *this; // CRITICAL: Return reference to this object for chaining + } + + // Friend function to overload << operator for output streams + friend std::ostream& operator<<(std::ostream& os, const Button& btn) { + os << (btn.value_m ? "PUSHED" : "not PUSHED"); + return os; + } + + private: + bool value_m = false; + bool priorState_m = false; +};// Button + + +/** + * @brief Host-side 1D mask view (0/1) used for ghost/halo flags. + */ +using HostMaskView1D_t = Kokkos::View; + + +/** + * @brief Cache key for ghosted data. + * + * Tuple components: + * - pointer identifying the mesh/topology + * - pointer identifying the owning view/container + * - size/extent of the ghost region + */ +using GhostKey_t = std::tuple; + + +/** + * @brief Hash functor for GhostKey_t used in unordered caches. + */ +struct GhostKeyHash { + std::size_t operator()(const GhostKey_t& k) const { + // Get the hash for each element in the tuple + auto h1 = std::hash{}(std::get<0>(k)); + auto h2 = std::hash{}(std::get<1>(k)); + auto h3 = std::hash{}(std::get<2>(k)); + + // Combine the hashes. This is a common pattern (based on boost::hash_combine) + // It xors and bit-shifts to mix the bits well. + h1 ^= h2 + 0x9e3779b9 + (h1 << 6) + (h1 >> 2); + h1 ^= h3 + 0x9e3779b9 + (h1 << 6) + (h1 >> 2); + return h1; + } +}; + + +/** + * @class CatalystAdaptor + * @brief High-level orchestrator for Catalyst initialization, execution, and steering. + * + * The adaptor wires up registered fields/particles and steerables to Conduit `node_m` instances, + * invokes Catalyst scripts, and forwards/fetches steering values between the + * simulation and the GUI. Use the runtime registry API (Initialize/Execute) + * for flexible, non-templated integration. + */ +class CatalystAdaptor { +public: + struct InitVisitor; + struct ExecVisitor; + struct SteerInitVisitor; + struct SteerForwardVisitor; + struct SteerFetchVisitor; + + using VisVisitorVariant_t = std::variant; + using SteerVisitorVariant_t = std::variant; + +private: + static std::string environmentValue(const char* name, std::string fallback) { + const char* value = std::getenv(name); + return value && *value ? std::string(value) : std::move(fallback); + } + + static std::filesystem::path catalystOutputDirectory() { + const char* value = std::getenv("IPPL_CATALYST_OUTPUT_DIR"); + return value && *value ? std::filesystem::path(value) + : std::filesystem::current_path() / "catalyst"; + } + + static std::optional normalizeExperimentName( + std::optional experimentName) { + if (experimentName && experimentName->empty()) { + experimentName.reset(); + } + return experimentName; + } + + std::shared_ptr visRegistry_m; + std::shared_ptr steerRegistry_m; + + ViewRegistry viewRegistry_m; + conduit_cpp::Node node_m; + conduit_cpp::Node results_m; + public: + Inform catalystInfo_m; + Inform catalystWarn_m; + private: + + ProxyWriter proxyWriter_m; + std::unordered_map>> enumChoicesByType_m; + std::unordered_map>> enumChoices_m; + + + const std::optional experimentName_m; + const std::string catalystVis_m; + const std::string catalystLive_m; + const std::string catalystSteer_m; + const std::string catalystPng_m; + const std::string catalystVtk_m; + const std::string catalystVerbosity_m; + const std::string catalystGhostMask_m; + const std::string proxyOption_m; + + const bool visEnabled_m; + const bool liveEnabled_m; + const bool steerEnabled_m; + const bool pngExtracts_m ; + const bool vtkExtracts_m ; + const int outputLevel_m ; + const bool useGhostMasks_m; + + const std::filesystem::path resourceDir_m; + const std::filesystem::path outputDir_m; + + std::unordered_map forceHostCopy_m; + + std::unordered_map ghostMaskCache_m; + +public: + + + + + /** + * @brief Construct a Catalyst adaptor. + * + * @param experimentName Optional label forwarded to Catalyst extractor scripts. + * Runtime output is written below IPPL_CATALYST_OUTPUT_DIR, or ./catalyst when unset. + */ + explicit CatalystAdaptor(std::optional experimentName = std::nullopt) : + catalystInfo_m("CatalystAdaptor::", 0), // Only print on rank 0 + catalystWarn_m("CatalystAdaptor_WARNING", std::cerr, INFORM_ALL_NODES), + experimentName_m(normalizeExperimentName(std::move(experimentName))), + catalystVis_m(environmentValue("IPPL_CATALYST_VIS", "ON")), + catalystLive_m(environmentValue("IPPL_CATALYST_LIVE", "OFF")), + catalystSteer_m(environmentValue("IPPL_CATALYST_STEER", "OFF")), + catalystPng_m(environmentValue("IPPL_CATALYST_PNG", "OFF")), + catalystVtk_m(environmentValue("IPPL_CATALYST_VTK", "OFF")), + catalystVerbosity_m(environmentValue( + "IPPL_CATALYST_VERBOSITY", std::to_string(ippl::Info->getOutputLevel()))), + catalystGhostMask_m(environmentValue("IPPL_CATALYST_GHOST_MASKS", "OFF")), + proxyOption_m(environmentValue("IPPL_CATALYST_PROXY_OPTION", "ON")), + visEnabled_m(catalystVis_m != "OFF"), + liveEnabled_m(catalystLive_m == "ON"), + steerEnabled_m(catalystSteer_m == "ON"), + pngExtracts_m(catalystPng_m == "ON"), + vtkExtracts_m(catalystVtk_m == "ON"), + outputLevel_m(std::stoi(catalystVerbosity_m)), + useGhostMasks_m(catalystGhostMask_m == "ON"), + resourceDir_m(IPPL_CATALYST_SCRIPTS_DIR), + outputDir_m(catalystOutputDirectory()) + { + catalystInfo_m.setOutputLevel(outputLevel_m); + + // #if defined(MPI_VERSION) + // MPI_Barrier(MPI_COMM_WORLD); + // if(ippl::Comm->rank()==0) catalystWarn_m << "[rank = 0 size=" << ippl::Comm->size() << "]" << endl; + // MPI_Barrier(MPI_COMM_WORLD); + // if(ippl::Comm->rank()==1) catalystWarn_m << "[rank= 1 size=" << ippl::Comm->size() << "]" << endl; + // MPI_Barrier(MPI_COMM_WORLD); + // #endif + + catalystInfo_m << level4 << "::CatalystAdaptor() Global Output Level setting: " << ippl::Info->getOutputLevel() << endl; + catalystInfo_m << level4 << "::CatalystAdaptor() Catalyst Info Output Level setting: " << catalystInfo_m.getOutputLevel() << endl; + catalystInfo_m << level4 << "::CatalystAdaptor() Catalyst Warn Output Level setting: " << catalystWarn_m.getOutputLevel() << endl; + catalystInfo_m << level4 << "::CatalystAdaptor() using resourceDir_m = " << resourceDir_m.string() << endl; + catalystInfo_m << level4 << "::CatalystAdaptor() using outputDir_m = " << outputDir_m.string() << endl; + if (pngExtracts_m) + { catalystInfo_m << level4 << "::CatalystAdaptor() PNG extraction ACTIVATED" << endl;} + else{ catalystInfo_m << level4 << "::CatalystAdaptor() PNG extraction DEACTIVATED" << endl;} + if (vtkExtracts_m) + { catalystInfo_m << level4 << "::CatalystAdaptor() VTK extraction ACTIVATED" << endl;} + else{ catalystInfo_m << level4 << "::CatalystAdaptor() VTK extraction DEACTIVATED" << endl;} + if (steerEnabled_m) + { catalystInfo_m << level4 << "::CatalystAdaptor() Steering ACTIVATED" << endl;} + else{ catalystInfo_m << level4 << "::CatalystAdaptor() Steering DEACTIVATED" << endl;} + + } + + private: + + // ============================================================================================== + // HELPERS ===================================================================================== + // ============================================================================================== + /** + * @brief Sets a file path to a node, using an environment variable if available, otherwise a default path. + * + * @param nodePath The Conduit `node_m` instance to set the file path in. + * @param envVar The name of the environment variable to check. + * @param defaultFilePath The default file path to use if the environment variable is not set or invalid. + */ + void setNodeScript( + conduit_cpp::Node nodePath, + // const char* envVar, + const std::string envVar, + const std::filesystem::path defaultFilePath + ); + + + // ========================================================== + // VISUALIZATION CHANNEL INITIALIZERS ======================= + // ========================================================== + + + + /* SCALAR FIELDS - handles both reference and shared_ptr */ + // == ippl::Field, Cell>* + /** + * @brief Initializes a Conduit `node_m` instance entry for a scalar field. + * + * @tparam T Field value type. + * @tparam Dim Field dimension. + * @tparam ViewArgs Additional template arguments for the field. + * @param entry The scalar field to initialize. + * @param label The label for the field/channel. + */ + template + void InitVizChannel( + [[maybe_unused]] + const Field& entry + , const std::string label + ); + + + /* VECTOR FIELDS - handles both reference and shared_ptr */ + // == ippl::Field, 3, ippl::UniformCartesian, Cell>* + /** + * @brief Initializes a Conduit `node_m` instance entry for a vector field. + * + * @tparam T Vector value type. + * @tparam Dim Field dimension. + * @tparam Dim_v Vector dimension. + * @tparam ViewArgs Additional template arguments for the field. + * @param entry The vector field to initialize. + * @param label The label for the field/channel. + */ + template + void InitVizChannel( + [[maybe_unused]] + const Field, Dim, ViewArgs...>& entry + , const std::string label + ); + + + + // PARTICLECONTAINERS DERIVED FROM PARTICLEBASE: + // == ippl::ParticleBaseBase -> ParticleBase, ... , ... > + /** + * @brief Initializes a Conduit `node_m` instance entry for a particle container derived from ParticleBaseBase. + * + * @tparam T Particle container type (must derive from ippl::ParticleBaseBase). + * @param entry The particle container to initialize. + * @param label The label for the container/channel. + */ + template + requires std::derived_from, ParticleBaseBase> + void InitVizChannel( + [[maybe_unused]] + const T& entry + , const std::string label + ); + + + /* SHARED_PTR DISPATCHER - automatically unwraps and dispatches to appropriate overload */ + /** + * @brief Dispatcher for InitVizChannel: unwraps shared_ptr and dispatches to the appropriate overload. + * + * @tparam T Entry type. + * @param entry Shared pointer to the entry. + * @param label The label for the entry/channel. + */ + template + void InitVizChannel( + const std::shared_ptr& entry + , const std::string label + ); + + + + // BASE CASE: + // only enabled if EntryT is NOT derived from ippl::ParticleBaseBase + /** + * @brief Fallback for InitVizChannel: handles types not derived from ParticleBaseBase + * and not having specific overloads. Should never be called since these types are already + * filtered out inside a Visitor struct with AllowedVisType_v. + * + * @tparam T Entry type. + * @param entry The entry to initialize (not a particle container). + * @param label The label for the entry/channel. + */ + template + requires (!std::derived_from, ParticleBaseBase>) + void InitVizChannel( + [[maybe_unused]] + // T&& entry + const T& entry + , const std::string label + ); + + + + // ========================================================== + // VISUALIZATION CHANNEL EXECUTIONERS ======================= + // ========================================================== + + + /* SCALAR FIELDS - handles both reference and shared_ptr */ + /* VECTOR FIELDS - handles both reference and shared_ptr */ + // == ippl::Field, Cell> + // == ippl::Field, 3, ippl::UniformCartesian, Cell> + /** + * @brief Executes a scalar/vector field entry, populating the Conduit `node_m` instance and updating the view registry. + * + * @tparam T Field value type. + * @tparam Dim Field dimension. + * @tparam ViewArgs Additional template arguments for the field. + * @param entry The scalar/vector field to execute. + * @param label The label for the field/channel. + * + */ + template + void ExecVizChannel( + const Field& entry + , const std::string label + ); + + + // PARTICLECONTAINERS DERIVED FROM PARTICLEBASE: + /** + * @brief Executes a particle container entry (derived from ParticleBaseBase), populating the Conduit `node_m` instance and updating the view registry. + * + * @tparam T Particle container type (must derive from ippl::ParticleBaseBase). + * @param entry The particle container to execute. + * @param label The label for the container/channel. + */ + template + requires (std::derived_from, ParticleBaseBase>) + void ExecVizChannel( + const T& entry + , const std::string label + ); + + + // BASE CASE: only enabled if EntryT is NOT derived from ippl::ParticleBaseBase + /** + * @brief Fallback for ExecVizChannel: handles types not derived from ParticleBaseBase. + * + * @tparam T Entry type. + * @param label The label for the entry/channel. + * @param entry The entry to execute (not a particle container). + */ + template + requires (!std::derived_from, ParticleBaseBase>) + void ExecVizChannel( + [[maybe_unused]] T&& entry + , const std::string label + ); + + + /* SHARED_PTR DISPATCHER - automatically unwraps and dispatches to appropriate overload */ + /** + * @brief Dispatcher for ExecVizChannel: unwraps shared_ptr and dispatches to the appropriate overload. + * + * @tparam T Entry type. + * @param entry Shared pointer to the entry. + * @param label The label for the entry/channel. + */ + template + void ExecVizChannel( + const std::shared_ptr& entry + , const std::string label + ); + + + + // ========================================================== + // STEERING CHANNEL INITIALIZERS===== ======================= + // ========================================================== + + /** + * @brief Initializes a steerable channel in the Conduit `node_m` instance for runtime parameter adjustment. + * + * @tparam T Type of the steerable parameter. + * @param steerableScalarForwardpass The initial value to set. + * @param label The label for the steerable channel. + */ + template + requires (!std::is_enum_v>) + void InitSteerChannel( const T& steerableScalarForwardpass, const std::string& label ); + + // Enum overloads (arbitrary enum types) + template + requires (std::is_enum_v>) + void InitSteerChannel( const E& e, const std::string& label ); + + // Bool-like switch overload + void InitSteerChannel( const bool& sw, const std::string& label ); + + // Vector overloads for steerable channels + template + void InitSteerChannel( const ippl::Vector& steerableVecForwardpass, const std::string& label ); + + + void InitSteerChannel( const ippl::Button& btn, const std::string& label ); + + // Generic std::vector elements (arithmetic/bool/Button) for array steerables + template + requires (std::is_arithmetic_v> || std::is_enum_v> || std::is_same_v, bool> || std::is_same_v, ippl::Button>) + void InitSteerChannel( const std::vector& arr, const std::string& label ); + + // std::vector of ippl::Vector steerables + template + void InitSteerChannel( const std::vector>& arr, const std::string& label ); + + + // ========================================================== + // STEERING CHANNEL EXECUTIONERS============================= + // ========================================================== + + /** + * @brief Adds a steerable channel to the Conduit `node_m` instance for runtime parameter adjustment. + * + * @tparam T Type of the steerable parameter. + * @param steerableScalarForwardpass The value to pass forward. + * @param steerableSuffix Suffix for the steerable channel name. + */ + template + requires (!std::is_enum_v>) + void ForwardSteerChannel(const T& steerableScalarForwardpass, const std::string& steerableSuffix); + + // Enum overloads (arbitrary enum types) + template + requires (std::is_enum_v>) + void ForwardSteerChannel(const E& e, const std::string& steerableSuffix); + + // Bool-like switch overload + void ForwardSteerChannel(const bool& sw, const std::string& steerableSuffix); + + // Button-like push overloads + void ForwardSteerChannel(const ippl::Button& btn, const std::string& steerableSuffix); + + // Vector overloads for steerable channels + template + void ForwardSteerChannel(const ippl::Vector& steerableVecForwardpass, const std::string& steerableSuffix); + + // Generic std::vector elements (arithmetic/bool/Button) for array steerables + template + requires (std::is_arithmetic_v> || std::is_enum_v> || std::is_same_v, bool> || std::is_same_v, ippl::Button>) + void ForwardSteerChannel(const std::vector& arr, const std::string& label); + + // std::vector of ippl::Vector arrays forward + template + void ForwardSteerChannel(const std::vector>& arr, const std::string& label); + + + + // ========================================================== + // STEERING CHANNEL FETCHERS ============================= + // ========================================================== + + /** + * @brief Fetches the value of a steerable channel from Catalyst results. + * + * @tparam T Type of the steerable parameter. + * @param steerableScalarBackwardpass Reference to store the fetched value. + * @param steerableSuffix Suffix for the steerable channel name. + */ + template + requires (!std::is_enum_v>) + void FetchSteerChannel( T& steerableScalarBackwardpass, const std::string& steerableSuffix); + + // Enum overloads (arbitrary enum types) + template + requires (std::is_enum_v>) + void FetchSteerChannel( E& e, const std::string& steerableSuffix); + + // ippl Vector overload + template + void FetchSteerChannel( ippl::Vector& steerableVecBackwardpass, const std::string& steerableSuffix); + + // standard vector overload for basic types + template + requires (std::is_arithmetic_v> || std::is_enum_v> || std::is_same_v, bool> || std::is_same_v, ippl::Button>) + void FetchSteerChannel( std::vector& out, const std::string& label); + + // std::vector of ippl::Vector arrays backward + template + void FetchSteerChannel( std::vector>& out, const std::string& label); + + + /** + * @brief Retrieves results from Catalyst and populates the given Conduit `node_m` instance. + * + */ + void fetchResults(); + + + // ===================================================================================== + // CatalystAdaptor Public Methods + // ===================================================================================== + public: + + /** + * @brief Struct steering registration. + * + * Expose any user struct composed of already supported steerable member + * types (arithmetic, bool, ippl::Button, ippl::Vector<>, enums). Nested structs are not + * supported yet. Must be called before adding the struct instance to a + * runtime registry. Args must be an even-sized pack of: name(string-like), pointer-to-member. + * Validates each member type at registration (throws IpplException if invalid). + * Stores three lambdas which expand the pack and delegate to visitor overloads. + * + * @tparam T The struct type to register. + * @tparam Args Types of the member pointers. + * @param args Pointers to the struct members. + */ + template + static void RegisterStructMembers(Args&&... args); + + + + /* DEPRECATED atm ... */ + // Optional enum metadata: label -> list of (text,value) choices + // std::unordered_map>> enumChoices_m; + // /** + // * @brief Provides enum choices metadata so the GUI shows a dropdown. + // * + // * Use this before InitializeRuntime. + // * Example: RegisterEnumChoices("mode", {{"Off",0},{"Basic",1},{"Advanced",2}}); + // * + // * @param label The label for the enum choice. + // * @param entries A vector of pairs containing the display name and integer value. + // */ + // void RegisterEnumChoices(const std::string& label, + // const std::vector>& entries) { + // enumChoices_m[label] = entries; + // } + + // /** + // * @brief Provides typed enum choices metadata mapped to a specific label for the GUI dropdown. + // * + // * @tparam E The enumeration type. + // * @param label The specific label of the steerable channel to associate the choices with. + // * @param entries A vector of pairs containing the display name and the enumeration value. + // */ + // template + // requires (std::is_enum_v>) + // void RegisterEnumChoicesTyped(const std::string& label, const std::vector>& entries); + + /** + * @brief Provides typed enum choices metadata registered by the enum type, to be reused across labels. + * + * @tparam E The enumeration type. + * @param entries A vector of pairs containing the display name and the enumeration value. + */ + template + requires (std::is_enum_v>) + void RegisterEnumChoicesTyped(const std::vector>& entries); + + + + + /** + * @brief Initializes Catalyst using runtime registries (visualization and steering). + * + * @param visReg Shared pointer to the visualization runtime registry. + * @param steerReg Shared pointer to the steering runtime registry. + */ + void Initialize( + const std::shared_ptr& visReg, + const std::shared_ptr& steerReg + ); + + /** + * @brief Explicitly forces a host copy for a specifically labelled channel right now. + * + * @param label The label specifying which field or particle to remember currently. + */ + void rememberNow(const std::string label); + + /** + * @brief Executes Catalyst for a given timestep using the runtime registry. + * + * Populates forward steerable values and fetches back updated ones. + * + * @param cycle The current simulation cycle or timestep index. + * @param time The current simulation time. + * @param rank The MPI rank of the executing process (defaults to ippl::Comm->rank()). + */ + void Execute( + int cycle, double time, + int rank = ippl::Comm->rank() + ); + + + /** + * @brief Finalizes Catalyst and releases resources. + */ + void Finalize(); + + + +};//class CatalystAdaptor +} //namespace ippl + +#include "Stream/InSitu/CatalystVisitors.h" +#include "Stream/Registry/VisRegistryRuntime.h" // visitor structs +#include "CatalystAdaptor.hpp" + + + +#endif diff --git a/src/Stream/InSitu/CatalystAdaptor.hpp b/src/Stream/InSitu/CatalystAdaptor.hpp new file mode 100644 index 000000000..b195628d6 --- /dev/null +++ b/src/Stream/InSitu/CatalystAdaptor.hpp @@ -0,0 +1,1227 @@ +#pragma once +#include "Stream/InSitu/CatalystAdaptor.h" + +#include "Ippl.h" +#include +#include +#include +#include +namespace ippl{ + +// ============================================================================================== +// HELPERS ===================================================================================== +// ============================================================================================== +void CatalystAdaptor::setNodeScript( + conduit_cpp::Node nodePath, + const std::string envVar, + const std::filesystem::path defaultFilePath +){ + const char* filePathEnv = std::getenv(envVar.c_str()); + std::filesystem::path filePath; + if (filePathEnv && std::filesystem::exists(filePathEnv)) { + catalystInfo_m << level4 <<"::Initialize()::setNodeScripts(...):\n" + << " Using " << envVar << " from environment:\n" + << " "<< filePathEnv << endl; + filePath = filePathEnv; + } else { + catalystInfo_m << level4 <<"::Initialize()::setNodeScripts(...): No valid " << envVar <<" set.\n" + << " Using default:\n" + << " " << defaultFilePath << endl; + filePath = defaultFilePath; + } + nodePath.set(filePath.string()); +} + + +// ============================================================================================== +// VIS CHANNEL INITIALIZER:====================================================================== +// ============================================================================================== + +// == ippl::Field, Cell> +template +void CatalystAdaptor::InitVizChannel( [[maybe_unused]] const Field& entry , const std::string label) +{ + catalystInfo_m << level4 <<"::Initialize()::InitVizChannel(ippl::Field<" << typeid(T).name() << "," << Dim << ">) called" << endl; + forceHostCopy_m[label] = false; + + const std::string channelName = "ippl_sField_" + label; + if(pngExtracts_m){ + const std::string script = "catalyst/scripts/" + label; + setNodeScript( node_m[script + "/filename"], + "CATALYST_EXTRACTOR_SCRIPT_" +label, + resourceDir_m / "catalyst_extractors" / "png_ext_sfield.py" + ); + conduit_cpp::Node args = node_m[script + "/args"]; + args.append().set_string("--channel_name"); + args.append().set_string(channelName); + args.append().set_string("--label"); + args.append().set_string(label); + if(experimentName_m){ + args.append().set_string("--experiment_name"); + args.append().set_string(*experimentName_m); + } + + args.append().set_string("--verbosity"); + args.append().set_string(std::to_string(catalystInfo_m.getOutputLevel())); + } + + conduit_cpp::Node scriptArgs = node_m["catalyst/scripts/script/args"]; + scriptArgs.append().set_string(channelName); +} + + +// == ippl::Field, 3, ippl::UniformCartesian, Cell> +template +void CatalystAdaptor::InitVizChannel( [[maybe_unused]] const Field, Dim, ViewArgs...>& entry , const std::string label) +{ + catalystInfo_m << level4 << "::Initialize()::InitVizChannel(ippl::Field," << Dim + << ">) called" << endl; + + forceHostCopy_m[label] = false; + + + const std::string channelName = "ippl_vField_" + label; + if(pngExtracts_m){ + const std::string script = "catalyst/scripts/" + label; + + setNodeScript( node_m[script + "/filename"], + "CATALYST_EXTRACTOR_SCRIPT_" + label, + resourceDir_m / "catalyst_extractors" / "png_ext_vfield.py" + + ); + conduit_cpp::Node args = node_m[script + "/args"]; + args.append().set_string("--channel_name"); + args.append().set_string(channelName); + args.append().set_string("--label"); + args.append().set_string(label); + if(experimentName_m){ + args.append().set_string("--experiment_name"); + args.append().set_string(*experimentName_m); + } + args.append().set_string("--verbosity"); + args.append().set_string(std::to_string(catalystInfo_m.getOutputLevel())); + } + + + conduit_cpp::Node scriptArgs = node_m["catalyst/scripts/script/args"]; + scriptArgs.append().set_string(channelName); + +} + +// PARTICLECONTAINERS derived from ParticleBaseBase: +// == ippl::ParticleBase,...>,...> +template +requires (std::derived_from, ParticleBaseBase>) +void CatalystAdaptor::InitVizChannel( [[maybe_unused]] const T& entry, const std::string label) +{ + catalystInfo_m << level4 << "::Initialize()::InitVizChannel(ParticleBase).name() << ","<< particle_dim_v + << ",...>...> [or subclass]) called" << endl; + + forceHostCopy_m[label] = false; + + const std::string channelName = "ippl_particles_" + label; + if(pngExtracts_m){ + const std::string script = "catalyst/scripts/"+ label; + + setNodeScript( + node_m[script + "/filename"], + "CATALYST_EXTRACTOR_SCRIPT_" +label, + resourceDir_m / "catalyst_extractors" / "png_ext_particle.py" + ); + + conduit_cpp::Node args = node_m[script + "/args"]; + args.append().set_string("--channel_name"); + args.append().set_string(channelName); + args.append().set_string("--label"); + args.append().set_string(label); + + if(experimentName_m){ + args.append().set_string("--experiment_name"); + args.append().set_string(*experimentName_m); + } + + args.append().set_string("--verbosity"); + args.append().set_string(std::to_string(catalystInfo_m.getOutputLevel())); + } + + conduit_cpp::Node scriptArgs = node_m["catalyst/scripts/script/args"]; + scriptArgs.append().set_string(channelName); +} + +/* SHARED_PTR DISPATCHER - automatically unwraps and dispatches to appropriate overload */ +template +void CatalystAdaptor::InitVizChannel( const std::shared_ptr& entry, const std::string label) +{ + if (entry) { + InitVizChannel( *entry + , label + ); + } + else { + catalystWarn_m << "::Initialize()InitVizChannel(nullptr): nullptr passed as entry." << endl + << " ID: "<< label << endl + << " ==> Channel will not be registered in Conduit Node." << endl; + } +} + + +// BASE CASE: +template +requires (!std::derived_from, ParticleBaseBase>) +void CatalystAdaptor::InitVizChannel([[maybe_unused]] const T& entry, const std::string label) +{ + catalystWarn_m << "::Initialize()InitVizChannel(nullptr): Entry type can't be processed." << endl + << " ID: "<< label << endl + << " Type: "<< typeid(std::decay_t).name() << endl + << " ==>Channel will not be registered in Conduit Node." << endl + << " If you see this something is wrong with: CatalystAdaptor::InitVisitor!!" << endl; +} + + + + +// == ippl::Field, Cell>* +// == ippl::Field, 3, ippl::UniformCartesian, Cell>* +template +void CatalystAdaptor::ExecVizChannel(const Field& entry, const std::string label) +{ + const bool refreshOnly = viewRegistry_m.contains(label); + if (refreshOnly && !forceHostCopy_m[label]) { + return; + } + + using Field_type = Field; + const Field_type* field = &entry; + + std::string channelName; + if constexpr (std::is_scalar_v) { + channelName = "ippl_sField_" + label; + catalystInfo_m << level4 <<"::Execute()::ExecVizChannel(" << label << ") | Type:ippl::Field<" << typeid(T).name() << "," << Dim << ">) called" << endl; + + } else if constexpr (is_vector_v) { + channelName = "ippl_vField_" + label; + catalystInfo_m << level4 <<"::Execute()::ExecVizChannel(" << label << ") | Type: ippl::Field::dim << ">," << Dim << ">)" << endl; + }else{ + channelName = "ippl_errorField_" + label; + + catalystInfo_m << level4 << "::Execute()::ExecVizChannel(Field<"<)\n" + << " For this type of Field the Conduit Blueprint description wasnt \n" + << " implemented in ippl. Therefore this type of field is not \n" + << " supported for visualisation." << endl; + } + + conduit_cpp::Node field_node; + typename Field_type::Layout_t& Layout_ = field->getLayout(); + typename Field_type::Mesh_t& Mesh_ = field->get_mesh(); + + const auto LocalNDIndex_ = Layout_.getLocalNDIndex(); + const auto Origin_ = Mesh_.getOrigin(); + const auto Spacing_ = Mesh_.getMeshSpacing(); + + const size_t nGhost = field->getNghost(); // returns int + + const size_t extra = (useGhostMasks_m) ? size_t(2*nGhost) : 0 ; + const size_t index_offset = (useGhostMasks_m) ? size_t(nGhost) : 0 ; + + int dims_n=1; + double extra_origin=0; + + static_assert(Dim >= 1 && Dim <= 3, + "Catalyst visualization supports fields with one to three dimensions"); + + const auto Ox = Origin_[0] + + (double(int(LocalNDIndex_[0].first()) - int(index_offset)) + + extra_origin) + * Spacing_[0]; + double Oy = 0.0; + double Oz = 0.0; + if constexpr (Dim >= 2) { + Oy = Origin_[1] + + (double(int(LocalNDIndex_[1].first()) - int(index_offset)) + extra_origin) + * Spacing_[1]; + } + if constexpr (Dim >= 3) { + Oz = Origin_[2] + + (double(int(LocalNDIndex_[2].first()) - int(index_offset)) + extra_origin) + * Spacing_[2]; + } + + const size_t nx = LocalNDIndex_[0].length() + extra; + size_t ny = 1; + size_t nz = 1; + if constexpr (Dim >= 2) { + ny = LocalNDIndex_[1].length() + extra; + } + if constexpr (Dim >= 3) { + nz = LocalNDIndex_[2].length() + extra; + } + + const auto& fullDeviceView = field->getView(); // original view + using DeviceView_t = typename Field::view_type; + using HostView_t = Kokkos::View< + typename DeviceView_t::data_type, + Kokkos::LayoutLeft, + Kokkos::HostSpace + >; + + auto makeHostView = [&](const char* viewLabel) -> HostView_t { + if constexpr (Dim == 1) { + return HostView_t(viewLabel, nx); + } else if constexpr (Dim == 2) { + return HostView_t(viewLabel, nx, ny); + } else { + return HostView_t(viewLabel, nx, ny, nz); + } + }; + + auto copyFieldToHost = [&](HostView_t& hostView, bool includeGhosts) { + if (includeGhosts) { + Kokkos::deep_copy(hostView, fullDeviceView); + } else if constexpr (Dim == 1) { + const auto r0 = Kokkos::make_pair(nGhost, nGhost + nx); + Kokkos::deep_copy(hostView, Kokkos::subview(fullDeviceView, r0)); + } else if constexpr (Dim == 2) { + const auto r0 = Kokkos::make_pair(nGhost, nGhost + nx); + const auto r1 = Kokkos::make_pair(nGhost, nGhost + ny); + Kokkos::deep_copy(hostView, Kokkos::subview(fullDeviceView, r0, r1)); + } else { + const auto r0 = Kokkos::make_pair(nGhost, nGhost + nx); + const auto r1 = Kokkos::make_pair(nGhost, nGhost + ny); + const auto r2 = Kokkos::make_pair(nGhost, nGhost + nz); + Kokkos::deep_copy(hostView, Kokkos::subview(fullDeviceView, r0, r1, r2)); + } + }; + + if (refreshOnly) { + HostView_t* hostMirrorFinal = viewRegistry_m.find(label); + if (!hostMirrorFinal) { + throw IpplException("Stream::InSitu::CatalystAdaptor::ExecVizChannel", + "Missing host mirror for refresh: " + label); + } + copyFieldToHost(*hostMirrorFinal, useGhostMasks_m); + return; + } + + // channel for this field of type mesh adheres to conduits mesh blueprint + auto channel = node_m["catalyst/channels/"+ channelName]; + auto channel_state = channel["state"]; + + channel["type"].set_string("mesh"); + auto data = channel["data"]; + + + auto fields = data["fields"]; + field_node = fields[label]; + data["topologies/fmesh_topo/type"].set_string("uniform"); + data["topologies/fmesh_topo/coordset"].set_string("cart_uniform_coords"); + data["coordsets/cart_uniform_coords/type"].set_string("uniform"); + + const void* meshKey = static_cast(&Mesh_); + const void* layoutKey = static_cast(&Layout_); + auto ghostKey = GhostKey_t{meshKey, layoutKey, nGhost}; + + const size_t localNumCells = nx * ny * nz; + const int rank = ippl::Comm->rank(); + + // using RankViewCells_t = Kokkos::View; + // RankViewCells_t rank_id_view_cells("rank_id_view_cells", localNumCells); + // auto host_policy = getRangePolicy(rank_id_view_cells); + + using RankViewCells_t = Kokkos::View; + using HostExecSpace = Kokkos::DefaultHostExecutionSpace; + RankViewCells_t rank_id_view_cells("rank_id_view_cells_3D", nx, ny, nz); + + if (localNumCells > 0) { + Kokkos::MDRangePolicy> host_policy( + {0, 0, 0}, // Start indices {i, j, k} + {nx, ny, nz} // End indices {i, j, k} + ); + Kokkos::parallel_for("fill_rank_ids_3D", host_policy, + KOKKOS_LAMBDA(const int i, const int j, const int k) { + rank_id_view_cells(i, j, k) = rank; + }); + } + + auto rank_field = fields["RankID"]; + rank_field["association"].set_string("element"); // associate_m); + rank_field["topology"].set_string("fmesh_topo"); + rank_field["volume_dependent"].set_string("false"); + if (localNumCells > 0) { + rank_field["values"].set_external(rank_id_view_cells.data(), localNumCells); + } else { + rank_field["values"].set_external(static_cast(nullptr), 0); + } + data["metadata/vtk_fields/RankID/attribute_type"].set_string("ProcessIds"); + viewRegistry_m.set(label + "_rank_id_cells", rank_id_view_cells); + + + // auto print_ranked_mesh_info = [&](){ + // catalystWarn_m << "[ rank=" << ippl::Comm->rank() << "]" + // << " | dims(points)=" << nx << "x" << ny << "x" << nz + // << " | ghost: " << nGhost + // << " | origin=(" << Ox << "," << Oy << "," << Oz << ")" + // << " | spacing=(" << Spacing_[0] << "," << (Dim>=2?Spacing_[1]:0)<< "," << (Dim>=3?Spacing_[2]:0) << ")" << endl; + // }; + + // #if defined(MPI_VERSION) + // MPI_Barrier(MPI_COMM_WORLD); + // if(ippl::Comm->rank()==0) print_ranked_mesh_info(); + // MPI_Barrier(MPI_COMM_WORLD); + // if(ippl::Comm->rank()==1) print_ranked_mesh_info(); + // MPI_Barrier(MPI_COMM_WORLD); + // #endif + + { + data["coordsets/cart_uniform_coords/dims/i"].set(nx+ dims_n ); + data["coordsets/cart_uniform_coords/spacing/dx"].set(Spacing_[0]); + data["coordsets/cart_uniform_coords/origin/x"].set( Ox ); + data["topologies/fmesh_topo/origin/x"].set( Ox ); + } + if constexpr(Dim >= 2){ + data["coordsets/cart_uniform_coords/dims/j"].set(ny+ dims_n); + data["coordsets/cart_uniform_coords/spacing/dy"].set(Spacing_[1]); + data["coordsets/cart_uniform_coords/origin/y"].set( Oy ); + data["topologies/fmesh_topo/origin/y"].set( Oy ); + } + if constexpr(Dim >= 3){ + data["coordsets/cart_uniform_coords/dims/k"].set(nz+ dims_n); + data["coordsets/cart_uniform_coords/spacing/dz"].set(Spacing_[2]); + data["coordsets/cart_uniform_coords/origin/z"].set( Oz ); + data["topologies/fmesh_topo/origin/z"].set( Oz ); + } + + // ================================== + // Prepare Field Data in a HostMirror + // ================================== + + // Version 1: Cut out Ghost Cells from data during a deep copy into a new Kokkos View. + auto getHostMirrorView_noGhosts = [&]() -> HostView_t { + HostView_t hostMirrorFinal = makeHostView("hostMirrorNoGhosts_LayoutLeft"); + // This single deep_copy now performs: + // - Device-to-Host transfer + // - LayoutRight-to-LayoutLeft transpose + // - Cutting Ghost Cells from data. + copyFieldToHost(hostMirrorFinal, false); + return hostMirrorFinal ; + + ///////////////////////////////////////////////////////////////////////////////////////////////////////// + // NOTE: + // for both lambdas: + // Technically, we need to use deep copy: 1. when explicitly forced 2. In case of Different memory spaces. + // 3.(?) when layouts don't match up + // Even if we can't use a subview directly since data of subview isn't meaningfull accessible in raw format + // with data() (so the subview can't be used by Conduit), but Conduit has it's own methods to access a + // substructures of arrays. See: Conduit Strided Structured Field descriptions. + // more efficient versions without default deep copies should be possible partially relying on shallow copies. + // + // if (!forceHostCopy_m[label] && std::is_same::value) + // HostView_t hostMirrorFinal = HostView_t("hostMirrorNoGhosts",nx,ny,nz); + // -> use data directy if possible eg + // HostView_t hostMirrorFinal = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), fullDeviceView); + // return hostMirrorFinal; + // + // + // (?) create_mirror_and_copy, adapts to spaces! but can convert from LayoutRight to LayoutLeft? + ///////////////////////////////////////////////////////////////////////////////////////////////////////// + }; + + + // Version 2: Mark Ghost Cells in data with VTK meta data. + auto getHostMirrorView_withGhosts = [&]() -> HostView_t + { + using m_t = unsigned char; // Element Type for Field Mask + + /* Define the N-D Host View types we needed for the upcoming section/copying. */ + using DeviceMaskView_t = typename Field::view_type; //Fetch original Field View Type + using HostMaskView1D_t = Kokkos::View; // Host STORAGE + using HostMaskView_t = Kokkos::View< + typename DeviceMaskView_t::data_type, // e.g., unsigned char*** + Kokkos::LayoutLeft, + Kokkos::HostSpace + >; // Host WRAPPER + HostMaskView1D_t hostMaskView1D; // Declare storage view + + + // --- START OF CACHING LOGIC --- + auto it = ghostMaskCache_m.find(ghostKey); + if (it != ghostMaskCache_m.end()) + { + // --- CACHE HIT ---// Re-use the existing ghost view from the cache + catalystInfo_m << level4 <<"::Execute()::ExecVizChannel(" << label << ") | GhostCache HIT" << endl; + hostMaskView1D = it->second; + } + else + { + // --- CACHE MISS --- + catalystInfo_m << level4 <<"::Execute()::ExecVizChannel(" << label << ") | GhostCache MISS" << endl; + + /* Allocate a mask field matching the source field (same mesh/layout/nghost) */ + Field ghostMaskField(Mesh_, Layout_, static_cast(nGhost)); + + /* Fill entire allocation (owned + ghosts) with 1 */ + ghostMaskField = static_cast(1); + DeviceMaskView_t deviceMaskView = ghostMaskField.getView(); + auto interior = ghostMaskField.template getFieldRangePolicy<>(); + + /* Fill inner cells with 0 */ + // TODO(?): use ippl dimension independent iterators + if constexpr (Dim == 1) { + Kokkos::parallel_for("ZeroOwnedMask1D", interior, KOKKOS_LAMBDA(const int i) { + deviceMaskView(i) = static_cast(0); + }); + } else if constexpr (Dim == 2) { + Kokkos::parallel_for("ZeroOwnedMask2D", interior, KOKKOS_LAMBDA(const int i, const int j) { + deviceMaskView(i,j) = static_cast(0); + }); + } else if constexpr (Dim == 3) { + Kokkos::parallel_for("ZeroOwnedMask3D", interior, KOKKOS_LAMBDA(const int i, const int j, const int k) { + deviceMaskView(i,j,k) = static_cast(0); + }); + } + Kokkos::fence(); + + /* Allocate the 1D host view that will own the memory */ + hostMaskView1D = HostMaskView1D_t("hostGhostMask_1D", deviceMaskView.size()); + + // Wrap the flat host allocation in a view with the field's actual rank. This + // gives deep_copy matching source and destination ranks while retaining a single + // flat allocation for Conduit's external array. + HostMaskView_t hostMaskView_N_Rank; + if constexpr (Dim == 1) { + hostMaskView_N_Rank = + HostMaskView_t(hostMaskView1D.data(), deviceMaskView.extent(0)); + } else if constexpr (Dim == 2) { + hostMaskView_N_Rank = HostMaskView_t(hostMaskView1D.data(), + deviceMaskView.extent(0), + deviceMaskView.extent(1)); + } else { + hostMaskView_N_Rank = HostMaskView_t(hostMaskView1D.data(), + deviceMaskView.extent(0), + deviceMaskView.extent(1), + deviceMaskView.extent(2)); + } + + + // The deep_copy now performs: + // - Device-to-Host transfer + // - LayoutRight-to-LayoutLeft + // - Since ND wraps 1D; Copied data is in the hostMaskView1D owned data + Kokkos::deep_copy(hostMaskView_N_Rank, deviceMaskView); + + // Store the 1D view (which owns the memory) in the cache, enough to keep in memory + ghostMaskCache_m[ghostKey] = hostMaskView1D; + // --- END OF CACHING LOGIC --- + } + + + + + // auto ghostMask_field_meta = data["metadata/vtk_fields/GhostMask_field"]; // can't chooses arbitrary name!!!! + auto ghostMask_field_meta = data["metadata/vtk_fields/vtkGhostType"]; + ghostMask_field_meta["attribute_type"] = "Ghosts"; // same as set string??... + // auto ghostMask_field_node = fields["ghostMask_field"]; // can't chooses arbitrary name!!!! must be vtkGhostType + auto ghostMask_field_node = fields["vtkGhostType"]; + ghostMask_field_node["association"].set_string("element"); //associate_m); // vs vertex ... + ghostMask_field_node["topology"].set_string("fmesh_topo"); + ghostMask_field_node["volume_dependent"].set_string("false"); + // ghostMask_field_node["values"].set_external(hostMaskView.data(), hostMaskView.size()); + ghostMask_field_node["values"].set_external(hostMaskView1D.data(), hostMaskView1D.size()); + + + //////////////////////////////////////////////////////////////////////////////////////////// + // Note: + // Field name in the conduit nodes for data/fields and metadata/vtk_fields has to coincide + // so the meta data can be properly associated with the data. + //////////////////////////////////////////////////////////////////////////////////////////// + + HostView_t hostMirrorFinal = makeHostView("hostMirrorWithGhosts_LayoutLeft"); + copyFieldToHost(hostMirrorFinal, true); + + return hostMirrorFinal; + // --- END FIX FOR MAIN FIELD --- + }; + + + HostView_t hostMirrorFinal = (useGhostMasks_m) ? getHostMirrorView_withGhosts() : getHostMirrorView_noGhosts(); + /* FOR BOTH CASES FINAL NODE SETTINGS ARE DONE AND WE HAVE THE DATA INSIDE hostMirrorFinal */ + using elem_t = std::remove_pointer_t; + // will return size of vector and amounts of vectors (not size of multiple doubles) + const auto n_elems = hostMirrorFinal.size(); + // Use true element size as stride (handles padding) + static constexpr size_t stride_bytes = sizeof(elem_t); + // offset is zero?? guaranteed? + const size_t offset = 0; + + + field_node["association"].set_string("element"); //associate_m); + field_node["topology"].set_string("fmesh_topo"); + field_node["volume_dependent"].set_string("false"); + if constexpr (std::is_scalar_v) { + // --- SCALAR FIELD CASE --- + field_node["values"].set_external(hostMirrorFinal.data(), n_elems); + } else if constexpr (is_vector_v) { + // --- VECTOR FIELD CASE --- + if (n_elems > 0) { + field_node["values/x"].set_external(&hostMirrorFinal.data()[0][0], n_elems, + offset, stride_bytes); + if constexpr (T::dim >= 2) { + field_node["values/y"].set_external(&hostMirrorFinal.data()[0][1], n_elems, + offset, stride_bytes); + } + if constexpr (T::dim >= 3) { + field_node["values/z"].set_external(&hostMirrorFinal.data()[0][2], n_elems, + offset, stride_bytes); + } + } else { + using component_type = typename T::value_type; + field_node["values/x"].set_external(static_cast(nullptr), 0); + if constexpr (T::dim >= 2) { + field_node["values/y"].set_external(static_cast(nullptr), 0); + } + if constexpr (T::dim >= 3) { + field_node["values/z"].set_external(static_cast(nullptr), 0); + } + } + } + // else { + // --- INVALID CASE --- + // } + + /* save view so data isn't discarded */ + viewRegistry_m.set(label, hostMirrorFinal); +} + + + +// == PARTICLECONTAINERS derived from ippl::ParticleBase,...>,...> +template +requires (std::derived_from, ParticleBaseBase>) +void CatalystAdaptor::ExecVizChannel(const T& entry, const std::string label) +{ + const bool refreshOnly = viewRegistry_m.contains(label); + if (refreshOnly && !forceHostCopy_m[label]) { + return; + } + + catalystInfo_m << level4 << "::Execute()::ExecVizChannel(" << label << ") | Type : ParticleBase).name() + << "," + << particle_dim_v + << ",...>...> [or subclass])" << endl; + const std::string channelName = "ippl_particles_" + label; + + auto particleContainer = &entry; + const size_t localNum = particleContainer->getLocalNum(); + assert((localNum == 0 || particleContainer->R.getView().data() != nullptr) + && "A non-empty R view should not be nullptr"); + + const std::string blockName = "block_allRanks"; + // const std::string blockName = "block_rank" + std::to_string(ippl::Comm->rank()); + + // channel for this particleContainer + // channel of type mesh adheres to conduits mesh blueprint + + + auto channel = node_m["catalyst/channels/"+ channelName]; + + + + channel["type"].set_string("multimesh"); + + + auto data = channel["data/block_main"]; + auto data_help = channel["data/block_help"]; + channel["assembly/main"] = "block_main"; + channel["assembly/help"] = "block_help"; + + //////////////////////////////////////////////////////// + // Note: + // Multimesh currently seems to have caused more headaches compared to what + // we have seemed to have gained from using it. So if the bugs or inconveniences stay + // in upcoming updates for ParaView catalyst we might want to use two normal meshes + // instead. + // channel["type"].set_string("mesh"); + //////////////////////////////////////////////////////// + + auto fields = data["fields"]; + data["type"].set_string("mesh"); + + const int rank = ippl::Comm->rank(); + + using IotaView_t = Kokkos::View; + using RankView_t = Kokkos::View; + IotaView_t iota_view("iota", localNum); + RankView_t rank_id_view("rank_id_view", localNum); + if (localNum > 0) { + using HostExecSpace = Kokkos::DefaultHostExecutionSpace; + Kokkos::RangePolicy host_policy(0, localNum); + Kokkos::parallel_for("fill_iota_host", host_policy, KOKKOS_LAMBDA(const int64_t i) { + iota_view(i) = i; + }); + Kokkos::parallel_for("fill_rank_ids", host_policy, KOKKOS_LAMBDA(const int64_t i) { + rank_id_view(i) = rank; + }); + } + + viewRegistry_m.set(label + "_iota", iota_view); + viewRegistry_m.set(label + "_rank_id", rank_id_view); + + // Creates a host-accessible mirror view and copies the data from the device view to the host. + using RAttrib_t = std::remove_reference_tR)>; + using hostMirror_R_t = typename RAttrib_t::host_mirror_type; + hostMirror_R_t R_hostMirror; + + if (forceHostCopy_m[label]) { + R_hostMirror = particleContainer->R.getHostMirror(); + Kokkos::deep_copy(R_hostMirror, particleContainer->R.getView()); + viewRegistry_m.set(R_hostMirror); + } else { + R_hostMirror = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), particleContainer->R.getView()); + viewRegistry_m.set(R_hostMirror); + } + + using hostMirror_ID_t = typename std::remove_reference_tID)>::host_mirror_type; + hostMirror_ID_t ID_hostMirror; + if constexpr (T::EnableIDs) { + if (forceHostCopy_m[label]) { + ID_hostMirror = particleContainer->ID.getHostMirror(); + Kokkos::deep_copy(ID_hostMirror, particleContainer->ID.getView()); + viewRegistry_m.set(label, ID_hostMirror); + } else { + ID_hostMirror = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), particleContainer->ID.getView()); + viewRegistry_m.set(label, ID_hostMirror); + } + } + + using PLayout_t = T::Layout_t; + // using vector_t = T::Layout_t::vector_type; + // using value_t = T::Layout_t::value_type; + using R_elem_t = std::remove_pointer_t; //avoids padding etc (?) Rattrib_t + static constexpr size_t R_stride_bytes = sizeof(R_elem_t); + + /* checks if playput it spatial or pure layout */ + if constexpr (has_getRegionLayout_v){ + + using RLayout_t = PLayout_t::RegionLayout_t; + using NDRegion_t = RLayout_t::NDRegion_t; + constexpr unsigned dim_ = PLayout_t::dim; + const NDRegion_t ndr = particleContainer->getLayout().getRegionLayout().getDomain(); + + /* HELPER COORDINATES TO PASS THE BOUNDING BOX in vtk format*/ + /* HELPER TOPOLOGY TO PASS THE BOUNDING BOX (??even needed??) in vtk format */ + data_help["coordsets/bound_helper_coords/type"].set_string("uniform"); + data_help["topologies/bound_helper_topo/coordset"].set_string("bound_helper_coords"); + data_help["topologies/bound_helper_topo/type"].set_string("uniform"); + /* create unfirom coordinate mesh only consisting of the corner points of the domain */ + { + data_help["coordsets/bound_helper_coords/dims/i"].set(2); + data_help["coordsets/bound_helper_coords/spacing/dx"].set( ndr[0].max() - ndr[0].min() ); + data_help["coordsets/bound_helper_coords/origin/x"].set( ndr[0].min() ); + // data_help["topologies/bound_helper_topo/origin/x"].set( ndr[0].min() ); + } + if constexpr(dim_ >= 2){ + data_help["coordsets/bound_helper_coords/dims/j"].set(2); + data_help["coordsets/bound_helper_coords/spacing/dy"].set( ndr[1].max()- ndr[1].min() ); + data_help["coordsets/bound_helper_coords/origin/y"].set( ndr[1].min() ); + // data_help["topologies/bound_helper_topo/origin/y"].set( ndr[1].min() ); + } + if constexpr(dim_ >= 3){ + data_help["coordsets/bound_helper_coords/dims/k"].set(2); + data_help["coordsets/bound_helper_coords/spacing/dz"].set( ndr[2].max()- ndr[2].min() ); + data_help["coordsets/bound_helper_coords/origin/z"].set( ndr[2].min() ); + // data_help["topologies/bound_helper_topo/origin/z"].set( ndr[2].min() ); + } + } + // else { + /* will use raw particle data instead .. */ + // } + + + /* ATTRIBUTES HARDCODED IN PARTICELBASE are identity ID and position R */ + /* EXPLICIT COORDINATES -> EACH PARTICLE'S POSITION */ + + ///////////////////////////////////////////////////////////// + // Note: + // For debuggng purposes checking distribution of particles on ranks... + // #if defined(MPI_VERSION) + // // MPI_Barrier(MPI_COMM_WORLD); + // // if(ippl::Comm->rank()==0) catalystInfo_m << level4 <<"[Rank 0] Local Particles: " << localNum << endl; + // MPI_Barrier(MPI_COMM_WORLD); + // catalystWarn_m << "Local Particles: " << localNum << endl; + // MPI_Barrier(MPI_COMM_WORLD); + // #endif + ///////////////////////////////////////////////////////////// + + data["coordsets/p_explicit_coords/type"].set_string("explicit"); + /* unstructured topology relying on per rank unique particle ID */ + data["topologies/p_unstructured_topo/coordset"].set_string("p_explicit_coords"); + data["topologies/p_unstructured_topo/type"].set_string("unstructured"); + data["topologies/p_unstructured_topo/elements/shape"].set_string("point"); + data["topologies/p_unstructured_topo/elements/connectivity"].set_external(iota_view.data(),particleContainer->getLocalNum()); + + // left hardcodeed we already have the hostViews (instead of integrating the into the loop) + + /* Process ID ATTRIBUTE */ + auto rank_field = fields["RankID"]; + rank_field["association"].set_string("vertex"); + rank_field["topology"].set_string("p_unstructured_topo"); + rank_field["volume_dependent"].set_string("false"); + rank_field["values"].set_external(rank_id_view.data(), localNum); + data["metadata/vtk_fields/RankID/attribute_type"].set_string("ProcessIds"); + + /* Global ID ATTRIBUTE (only when particle IDs are enabled in ParticleBase) */ + if constexpr (T::EnableIDs) { + auto id_field = fields["ParticleIDs"]; + id_field["association"].set_string("vertex"); + id_field["topology"].set_string("p_unstructured_topo"); + id_field["volume_dependent"].set_string("false"); + if (localNum > 0) { + id_field["values"].set_external(ID_hostMirror.data(), localNum); + } else { + using id_value_t = typename hostMirror_ID_t::value_type; + id_field["values"].set_external(static_cast(nullptr), 0); + } + data["metadata/vtk_fields/ParticleIDs/attribute_type"].set_string("GlobalIds"); + } + + /* POSITION ATTRIBUTE */ + auto R_field = fields["position"]; + R_field["association"].set_string("vertex"); + R_field["topology"].set_string("p_unstructured_topo"); + R_field["volume_dependent"].set_string("false"); + + + constexpr unsigned ParticleDim = particle_dim_v; + static_assert(ParticleDim >= 1 && ParticleDim <= 3, + "Catalyst visualization supports particles with one to three dimensions"); + + if (localNum > 0) + { + /* COORDINATE DEFINITION... */ + data["coordsets/p_explicit_coords/values/x"].set_external(&R_hostMirror.data()[0][0], particleContainer->getLocalNum(), 0, R_stride_bytes); + if constexpr (ParticleDim >= 2) { + data["coordsets/p_explicit_coords/values/y"].set_external(&R_hostMirror.data()[0][1], particleContainer->getLocalNum(), 0, R_stride_bytes); + } + if constexpr (ParticleDim >= 3) { + data["coordsets/p_explicit_coords/values/z"].set_external(&R_hostMirror.data()[0][2], particleContainer->getLocalNum(), 0, R_stride_bytes); + } + + /* POSITION ATTRIBUTE */ + R_field["values/x"].set_external(&R_hostMirror.data()[0][0], particleContainer->getLocalNum(), 0, R_stride_bytes); + if constexpr (ParticleDim >= 2) { + R_field["values/y"].set_external(&R_hostMirror.data()[0][1], particleContainer->getLocalNum(), 0, R_stride_bytes); + } + if constexpr (ParticleDim >= 3) { + R_field["values/z"].set_external(&R_hostMirror.data()[0][2], particleContainer->getLocalNum(), 0, R_stride_bytes); + } + + /* concept for no copy in situ vis would be */ + //mesh["topologies/p_unstructured_topo/elements/connectivity"].set_external(particleContainer->ID.getView().data(),particleContainer->getLocalNum()); + }else + { + // In case a rank has no particles-> data()[0] is nulllptr dereferencing !!!!! + using component_type = typename R_elem_t::value_type; + data["coordsets/p_explicit_coords/values/x"].set_external(static_cast(nullptr), 0); + R_field["values/x"].set_external(static_cast(nullptr), 0); + if constexpr (ParticleDim >= 2) { + data["coordsets/p_explicit_coords/values/y"].set_external(static_cast(nullptr), 0); + R_field["values/y"].set_external(static_cast(nullptr), 0); + } + if constexpr (ParticleDim >= 3) { + data["coordsets/p_explicit_coords/values/z"].set_external(static_cast(nullptr), 0); + R_field["values/z"].set_external(static_cast(nullptr), 0); + } + } + + // Attribute containers are grouped by memory space, so their global ordinal is not + // stable. Identify ParticleBase's built-ins by object identity instead. + entry.forAllAttributes([&](const Attributes& atts) { + for (auto* attribute : atts) { + if (!entry.isBuiltinAttribute(attribute)) { + attribute->signConduitBlueprintNode(localNum, fields, viewRegistry_m, + catalystInfo_m, catalystWarn_m, + forceHostCopy_m[label]); + } + } + }); + //////////////////////////////////////////////////////////////////////////////////////////// + // Note: + // All ways on how to iterate over particle attributes rely on base class pointers. + // In whih case dimensions and types particle attributes is not retrievable from the + // a pointer instance. Therefore conduit maniupulation are done as a membermethod + // for paricleAttrib (overriding a virtual method in the base class). + // + // entry.template forAllAttributes( + // [&](const Attributes& atts) { + // for (auto* attribute : atts) { + //////////////////////////////////////////////////////////////////////////////////////////// + +} + + +// BASE CASE: only enabled if EntryT is NOT derived from ippl::ParticleBaseBase +template +requires (!std::derived_from, ParticleBaseBase>) + void CatalystAdaptor::ExecVizChannel( [[maybe_unused]] T&& entry, const std::string label) +{ + catalystInfo_m << level4 <<" Entry type can't be processed: ID "<< label <<" "<< typeid(std::decay_t).name() << endl; +} + + +/* SHARED_PTR DISPATCHER - automatically unwraps and dispatches to appropriate overload */ +template + void CatalystAdaptor::ExecVizChannel( const std::shared_ptr& entry,const std::string label ) +{ + if (entry) { + catalystInfo_m << level4 <<" dereferencing shared pointer and reattempting execute..." << endl; + ExecVizChannel(*entry, label ); + } else { + catalystInfo_m << level4 <<" Null shared_ptr encountered" << endl; + } +} + + +} + + +// ===================================================================================== +// STEERING: +// ===================================================================================== +#include "Stream/InSitu/CatalystAdaptorSteering.hpp" +// ===================================================================================== + +namespace ippl{ +/////////////////////////////////////////////////////////////// +// Note: +// this does not really need to be a separate function call, +// May we should just inline this into the execute functions ... +/////////////////////////////////////////////////////////////// +void CatalystAdaptor::fetchResults() { + + catalyst_status err = catalyst_results(conduit_cpp::c_node(&results_m)); + if (err != catalyst_status_ok) + { + std::cerr << "Failed to execute Catalyst-results: " << err << std::endl; + } + } + + +// ===================================================================================== +// Runtime registry based Initialize / Execute (non-templated registry) +// ===================================================================================== + +void CatalystAdaptor::Initialize( + const std::shared_ptr& visReg, + const std::shared_ptr& steerReg + ) { +if ( !visEnabled_m) return; + + + catalystInfo_m << level4 <<"::Initialize() START============================================================= 0" << endl; + + int all_ready = 1; + #if defined(MPI_VERSION) + MPI_Allreduce(MPI_IN_PLACE, &all_ready, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + #endif + catalystInfo_m << level4 <<"::InitializeRuntime() Ranks ready for catalyst int: " << all_ready << " ranks" << endl; + + visRegistry_m = visReg; + steerRegistry_m = steerReg; + + + const int fcomm = MPI_Comm_c2f(MPI_COMM_WORLD); + const int64_t fcomm64 = static_cast(fcomm); + node_m["catalyst/mpi_comm"].set(fcomm64); + + + setNodeScript( node_m["catalyst/scripts/script/filename"], //where in node_m + "CATALYST_PIPELINE_PATH", // environment override + resourceDir_m / "pipeline_default.py") //default + ; + conduit_cpp::Node args = node_m["catalyst/scripts/script/args"]; + + + args.append().set_string("--channel_names"); + InitVisitor initV{*this}; + visRegistry_m->forEach(initV); + // Visitor will (also) append channel names here into the node_m (sequence of overall arguments is important!!) + + args.append().set_string("--verbosity"); + args.append().set_string(std::to_string(catalystInfo_m.getOutputLevel())); + + + args.append().set_string("--VTKextract"); + args.append().set_string(catalystVtk_m); + + args.append().set_string("--live"); + args.append().set_string(catalystLive_m); + + args.append().set_string("--steer"); + args.append().set_string(catalystSteer_m); + + + args.append().set_string("--steer_channel_names"); + + std::filesystem::path proxyPath = outputDir_m / "catalyst_proxy.xml"; + bool useExistingProxy = false; + if (const char* proxyPathEnv = std::getenv("IPPL_CATALYST_PROXY_PATH"); + proxyPathEnv && *proxyPathEnv) { + proxyPath = proxyPathEnv; + useExistingProxy = true; + } else if (const char* legacyProxyPathEnv = std::getenv("CATALYST_PROXYS_PATH"); + legacyProxyPathEnv && *legacyProxyPathEnv) { + proxyPath = legacyProxyPathEnv; + useExistingProxy = true; + catalystWarn_m << "CATALYST_PROXYS_PATH is deprecated; use " + "IPPL_CATALYST_PROXY_PATH instead." + << endl; + } + + std::string cfgYaml; + if (const char* cfg_env = std::getenv("IPPL_PROXY_CONFIG_YAML")) { + if (std::filesystem::exists(cfg_env)) { + cfgYaml = std::string(cfg_env); + } else { + catalystInfo_m << level4 <<"::Initialize() IPPL_PROXY_CONFIG_YAML set but file not found: '" << cfg_env << "', using default." << endl; + } + } + if (cfgYaml.empty()) { + auto default_cfgYaml = (resourceDir_m / "proxy_default_config.yaml").string(); + if (std::filesystem::exists(default_cfgYaml)) { + cfgYaml = std::move(default_cfgYaml); + } // else leave empty -> ProxyWriter can proceed without config + } + + + proxyWriter_m.initialize(proxyPath, cfgYaml); + if (steerEnabled_m ) { + SteerInitVisitor steerInitV{*this}; + steerRegistry_m->forEach(steerInitV); + } + const bool generateProxy = !useExistingProxy && proxyOption_m != "OFF"; + int proxyReady = 1; + + if (ippl::Comm->rank() == 0) { + if (useExistingProxy) { + proxyReady = std::filesystem::is_regular_file(proxyPath) ? 1 : 0; + } else if (generateProxy) { + proxyReady = proxyWriter_m.produceUnified( + "SteerableParameters_SCALARS", "SteerableParameters") + ? 1 + : 0; + } else { + proxyReady = std::filesystem::is_regular_file(proxyPath) ? 1 : 0; + } + } + +#if defined(MPI_VERSION) + MPI_Bcast(&proxyReady, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Barrier(MPI_COMM_WORLD); +#endif + + if ((useExistingProxy || generateProxy) && !proxyReady) { + throw IpplException( + "Stream::InSitu::CatalystAdaptor::Initialize()", + "Could not create or read Catalyst proxy XML: " + proxyPath.string()); + } + + if (proxyReady) { + node_m["catalyst/proxies/proxy_/filename"].set(proxyPath.string()); + } else if (steerEnabled_m) { + catalystWarn_m << "Steering is enabled without a Catalyst proxy XML because " + "IPPL_CATALYST_PROXY_OPTION=OFF and no existing proxy was found at " + << proxyPath.string() << endl; + } + + if (proxyOption_m == "PRODUCE_ONLY") { + throw IpplException( + "Stream::InSitu::CatalystAdaptor", + "write_proxy_only_run: proxy available at " + proxyPath.string()); + } + + catalystInfo_m << level4 <<"::Initialize() Printing Conduit `node_m` instance passed to catalyst_initialize() =>" << endl; + catalystInfo_m << level4 <second; + forceHostCopy_m[label] = true; + ExecVisitor execV{*this}; + const bool ok = visRegistry_m->forOne(label, execV); + + // Restore prior state + forceHostCopy_m[label] = tmp; + if (!ok) { + throw IpplException("Stream::InSitu::CatalystAdaptor::rememberNow", "Label not found in executable entries or has no execute callback: " + label); + } + +} + +void CatalystAdaptor::Execute( int cycle, double time, int rank /* default = ippl::Comm->rank() */) { + if ( !visEnabled_m) return; + + catalystInfo_m << level4 <<"::Execute() START =============================================================== 0" << endl; + + static IpplTimings::TimerRef TMRcatalyst_execute = IpplTimings::getTimer("catalyst_execute"); + static IpplTimings::TimerRef TMRexecVizVisitor = IpplTimings::getTimer("execVizVisitor"); + static IpplTimings::TimerRef TMRexecSteerVisitor = IpplTimings::getTimer("execSteerVisitor"); + + + auto state = node_m["catalyst/state"]; + state["cycle"].set(cycle); + state["time"].set(time); + state["domain_id"].set(rank); + + IpplTimings::startTimer(TMRexecVizVisitor); + if ( !!visEnabled_m){ + // edit forward Node: add visualisation channels + ExecVisitor execV{*this}; + visRegistry_m->forEach(execV); + } + IpplTimings::stopTimer(TMRexecVizVisitor); + + + IpplTimings::startTimer(TMRexecSteerVisitor); + if (steerEnabled_m) { + // edit forward Node: add steering channels + SteerForwardVisitor steerV{*this}; + steerRegistry_m->forEach(steerV); + } + IpplTimings::stopTimer(TMRexecSteerVisitor); + + + if(cycle == 0){ + + #if defined(MPI_VERSION) + MPI_Barrier(MPI_COMM_WORLD); + catalystInfo_m << level4 <<"::Execute() [rank = 0] Printing first Conduit Node passed from to catalyst_execute() ==>" << endl; + if(catalystInfo_m.getOutputLevel() >= 4 && ippl::Comm->rank()==0) node_m.print(); + catalystInfo_m << level4 <<"::Execute() [rank = 1] Printing first Conduit Node passed from to catalyst_execute() ==>" << endl; + MPI_Barrier(MPI_COMM_WORLD); + if(catalystInfo_m.getOutputLevel() >= 4 && ippl::Comm->rank()==1) node_m.print(); + MPI_Barrier(MPI_COMM_WORLD); + #endif + // if(level >= 5 && ippl::Comm->rank()==0) node_m.print(); + + catalystInfo_m << level4 << "::Execute() During first catalyst_execute() catalyst will\n" + << " for each passed script - in order how they were \n" + << " passed to the conduit node - run the globa scope,\n" + << " the initialize() and the execute()." << endl; + } + + + //////////////////////////////////////////////////////////////// + // Note: + // Possibly helpful for further debugging. + // + // Kokkos::fence(); + // #if defined(MPI_VERSION) + // MPI_Barrier(MPI_COMM_WORLD); + // #endif + // + // int all_ready = 1; + // #if defined(MPI_VERSION) + // MPI_Allreduce(MPI_IN_PLACE, &all_ready, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + // #endif + // catalystInfo_m << level4 <<"::Execute() All ranks ready for catalyst_execute: " << all_ready << " ranks" << endl; + //////////////////////////////////////////////////////////////// + + catalystInfo_m << level4 <<"::Execute()::catalyst_execute() ==>" << endl; + IpplTimings::startTimer(TMRcatalyst_execute); + catalyst_status err = catalyst_execute(conduit_cpp::c_node(&node_m)); + IpplTimings::stopTimer(TMRcatalyst_execute); + + //////////////////////////////////////////////////////////////////////////////// + // Note: + // catalyst execute seems to be the current bottleneck of a medium sized simulation... + //////////////////////////////////////////////////////////////////////////////// + + + if (err != catalyst_status_ok) { + std::cerr << "::Execute() Failed to execute Catalyst (runtime path): " << err << std::endl; + } + + if (steerEnabled_m) { + + static IpplTimings::TimerRef TMRfetchResult = IpplTimings::getTimer("fetchSteerParameters"); + IpplTimings::startTimer(TMRfetchResult); + + fetchResults(); + // backward Node: fetch updated steering values + SteerFetchVisitor fetchV{*this}; + steerRegistry_m->forEach(fetchV); + + IpplTimings::stopTimer(TMRfetchResult); + + + if(true){ + // if(cycle == 0){ + catalystInfo_m << level4 <<"::Execute() Printing Conduit Node received from catalyst_execute() ==>" << endl; + catalystInfo_m << level4 << results_m.to_yaml() << endl; + } + + } + + + viewRegistry_m.clear(); + ghostMaskCache_m.clear(); + node_m.reset(); + + /////////////////////////////////////////////////// + // Note: + // We deliberately don't reset results since + // 1. they will be properly overwritten by Catalyst. + // 2. If the part of the Catalyst backend crashes and the + // results are not sent back, the old results will be used, possibly + // avoiding a problems during result retrieval. + // + // results.reset(); + /////////////////////////////////////////////////// + + + catalystInfo_m << level4 <<"::Execute() DONE =============================================================== 1" << endl; + +} + + +void CatalystAdaptor::Finalize() { + if ( !visEnabled_m) return; + + conduit_cpp::Node node; + catalyst_status err = catalyst_finalize(conduit_cpp::c_node(&node_m)); + if (err != catalyst_status_ok) { + std::cerr << "Failed to finalize Catalyst: " << err << std::endl; + } +} + + + +}//ippl diff --git a/src/Stream/InSitu/CatalystAdaptorSteering.hpp b/src/Stream/InSitu/CatalystAdaptorSteering.hpp new file mode 100644 index 000000000..a53816a63 --- /dev/null +++ b/src/Stream/InSitu/CatalystAdaptorSteering.hpp @@ -0,0 +1,934 @@ +#pragma once + +#include "Stream/InSitu/CatalystAdaptor.h" +#include "Stream/InSitu/CatalystVisitors.h" // ensure AllowedSteerType_v available +#include +#include +#include +#include +#include + + + + +namespace ippl{ + +namespace detail { + + ///////////////////////////////////////////////// + // Note: + // Label sanitation: replace '/' to avoid + // unintended Conduit subtree splitting. + // TODO: Shoul be used also in Viz channels + ///////////////////////////////////////////////// + inline std::string sanitize_label(const std::string& in) { + std::string tmp; tmp.reserve(in.size()); + std::string out; out.reserve(in.size()); + for (char c : in) tmp += (c=='/' ? '_' : c); + for (char c : tmp) out += (c=='/' ? '_' : c); + + + return out; + } + + // Recursive applicators with const and non-const overloads to avoid const_cast. + // Base cases (no members left). + template + inline void apply_struct_members(Visitor&, T&, const std::string&) {} + template + inline void apply_struct_members(Visitor&, const T&, const std::string&) {} + + // Recursive steps: visit (name, pointer-to-member) then recurse. + template + inline void apply_struct_members(Visitor& vis, T& obj, const std::string& rootLabel, + NameType name, MemberPtr ptr, Rest&&... rest) { + // Build full label and dispatch to existing visitor overloads. + std::string fullLabel = rootLabel + "." + sanitize_label(std::string(name)); + vis(fullLabel, obj.*ptr); // rely on visitor operator() overload selection + apply_struct_members(vis, obj, rootLabel, std::forward(rest)...); + } + template + inline void apply_struct_members(Visitor& vis, const T& obj, const std::string& rootLabel, + NameType name, MemberPtr ptr, Rest&&... rest) { + std::string fullLabel = rootLabel + "." + sanitize_label(std::string(name)); + vis(fullLabel, obj.*ptr); + apply_struct_members(vis, obj, rootLabel, std::forward(rest)...); + } + + // Meta storage per struct type T. + template + struct StructMeta { + static inline bool registered = false; + // Dispatch lambda handles scalar/member steering. Takes a variant of visitors. + static inline std::function dispatch; + // Array-of-struct aggregation (vector) – built at registration time. + static inline std::function&, const std::string&)> dispatch_vec; + }; + +} // namespace detail + + + +// ===================================================================================== +// Registration functions +// ===================================================================================== +template +void CatalystAdaptor::RegisterStructMembers(Args&&... args) { + using DecayT = std::decay_t; + static_assert(sizeof...(Args) % 2 == 0, "RegisterStructMembers requires (name, memberPtr) pairs"); + if (detail::StructMeta::registered) return; + + // Pack names and member pointers interleaved. + auto pack = std::tuple(std::forward(args)...); + constexpr size_t N = sizeof...(Args); + constexpr size_t PairCount = N / 2; + + // Validate each (name, memberPtr) pair. + [&](std::index_sequence){ + (([] (auto name, auto memberPtr){ + using MemberType = std::decay_t().*memberPtr)>; + if constexpr (!AllowedSteerType_v) { + throw IpplException( + "CatalystAdaptor::RegisterStructMembers", + std::string("Unsupported member type for steering in struct '") + typeid(DecayT).name() + + "' member '" + std::string(name) + "'" + ); + } + })(std::get<2*I>(pack), std::get<2*I+1>(pack)), ...); + }(std::make_index_sequence{}); + + // Visitor dispatch reuses generic applicator. + detail::StructMeta::dispatch = [pack](CatalystAdaptor::SteerVisitorVariant_t var, DecayT& obj, const std::string& root) mutable { + std::visit([&](auto* active_ptr) { + if (!active_ptr) return; + auto& vis = *active_ptr; + std::apply([&](auto&&... all){ ippl::detail::apply_struct_members(vis, obj, ippl::detail::sanitize_label(root), all...); }, pack); + }, var); + }; + + // ================= Array-of-Struct (vector) support ================= + detail::StructMeta::dispatch_vec = [pack](CatalystAdaptor::SteerVisitorVariant_t var, std::vector& arr, const std::string& root) mutable { + if (arr.empty()) return; + constexpr size_t PC = PairCount; + std::visit([&](auto* active_ptr) { + if (!active_ptr) return; + auto& vis = *active_ptr; + using VisitorType = std::decay_t; + + [&](std::index_sequence){ + ( [&](){ + // using MemberPtrT = std::tuple_element_t<2*I+1, decltype(pack)>; + // MemberPtrT = // causese problems with "older" compilers... + auto mptr = std::get<2*I+1>(pack); + auto rawName = std::get<2*I>(pack); + std::string memberLabel = ippl::detail::sanitize_label(root) + '.' + ippl::detail::sanitize_label(std::string(rawName)); + using MType = std::remove_reference_t().*mptr)>; + std::vector tmp; tmp.reserve(arr.size()); + for (auto& el : arr) tmp.push_back(el.*mptr); + + vis(memberLabel, tmp); + + if constexpr (std::is_same_v) { + if (tmp.size() != arr.size()) arr.resize(tmp.size()); + // Write back member values + for (std::size_t i = 0; i < tmp.size(); ++i) { + arr[i].*mptr = tmp[i]; + } + } + }(), ... ); + }(std::make_index_sequence{}); + }, var); + }; + + detail::StructMeta::registered = true; +} + +/* DEPRECATED atm ... */ +// template +// requires (std::is_enum_v>) +// void CatalystAdaptor::RegisterEnumChoicesTyped(const std::string& label, const std::vector>& entries) { +// std::vector> conv; +// conv.reserve(entries.size()); +// for (const auto& p : entries) { +// conv.emplace_back(p.first, static_cast(p.second)); +// } +// RegisterEnumChoices(label, conv); +// } + +template +requires (std::is_enum_v>) +void CatalystAdaptor::RegisterEnumChoicesTyped(const std::vector>& entries) { + std::vector> conv; + conv.reserve(entries.size()); + for (const auto& p : entries) { + conv.emplace_back(p.first, static_cast(p.second)); + } + enumChoicesByType_m[std::type_index(typeid(E))] = std::move(conv); +} + + + + +// ===================================================================================== +// INITIALISATION: +// ===================================================================================== + +template +requires (!std::is_enum_v>) +void CatalystAdaptor::InitSteerChannel( [[maybe_unused]] const T& steerableScalarForwardpass, const std::string& label ){ + catalystInfo_m << level4 << "::Initialize()::InitSteerChannel(" << label << "): | Type: " << typeid(T).name() << endl; + // Only invoke ProxyWriter scalar include for arithmetic types; others are placeholders. + if constexpr (std::is_arithmetic_v>) { + proxyWriter_m.include(steerableScalarForwardpass, label); + } else { + catalystWarn_m << "ProxyWriter placeholder: include() for label '" << label + << "' (type=" << typeid(T).name() << ") not implemented yet (TODO)." << endl; + } + conduit_cpp::Node scriptArgs = node_m["catalyst/scripts/script/args"]; + scriptArgs.append().set_string(label); +} + +// Enum overload (explicit) to ensure proper dropdown setup even if template above is shadowed +template +requires (std::is_enum_v>) +void CatalystAdaptor::InitSteerChannel( [[maybe_unused]] const E& e, const std::string& label ){ + catalystInfo_m << level4 << "::Initialize()::InitSteerChannel(" << label << "): | Type: Enum" << endl; + auto it = enumChoices_m.find(label); + if (it != enumChoices_m.end()) { + proxyWriter_m.includeEnum(label, it->second, static_cast(e)); + } else { + // Try type-based enum choices + auto itt = enumChoicesByType_m.find(std::type_index(typeid(E))); + if (itt != enumChoicesByType_m.end()) { + proxyWriter_m.includeEnum(label, itt->second, static_cast(e)); + } else { + // Fallback to a checkbox if no choices registered + proxyWriter_m.includeBool(label, false); + } + } + conduit_cpp::Node scriptArgs = node_m["catalyst/scripts/script/args"]; + scriptArgs.append().set_string(label); +} + +// Bool-like Switch: init (checkbox in GUI) +void CatalystAdaptor::InitSteerChannel( [[maybe_unused]] const bool& sw, const std::string& label ){ + catalystInfo_m << level4 << "::Initialize()::InitSteerChannel(" << label << "): | Type: Switch" << endl; + proxyWriter_m.includeBool(label, static_cast(sw)); + + conduit_cpp::Node scriptArgs = node_m["catalyst/scripts/script/args"]; + scriptArgs.append().set_string(label); +} + +// Button-like: init (push button in GUI) +void CatalystAdaptor::InitSteerChannel( [[maybe_unused]] const ippl::Button& btn, const std::string& label ){ + catalystInfo_m << level4 << "::Initialize()::InitSteerChannel(" << label << "): | Type: Button" << endl; + proxyWriter_m.includeButton(label); + conduit_cpp::Node scriptArgs = node_m["catalyst/scripts/script/args"]; + scriptArgs.append().set_string(label); +} + +// Vector steerable overloads +template +void CatalystAdaptor::InitSteerChannel( [[maybe_unused]] const ippl::Vector& steerableVecForwardpass, const std::string& label ) +{ + catalystInfo_m << level4 << "::Initialize()::InitSteerChannel(" << label << "): | Vector<" << typeid(T).name() << "," << Dim_v << ">" << endl; + // Register this label as a vector channel in the proxy writer (limit to 3 comps in GUI) + proxyWriter_m.includeVector(label); + (void)steerableVecForwardpass; + + // Ensure the Python pipeline receives this label after `--steer_channel_names` + // so it creates the corresponding forward reader and unified sender wiring. + conduit_cpp::Node scriptArgs = node_m["catalyst/scripts/script/args"]; + scriptArgs.append().set_string(label); +} + + + +// ============================================================= +// Generic std::vector steerables (Array-of-struct members) +// Supports arithmetic, bool and ippl::Button as element types. +// ============================================================= +template +requires (std::is_arithmetic_v> || std::is_enum_v> || std::is_same_v, bool> || std::is_same_v, ippl::Button>) +void CatalystAdaptor::InitSteerChannel( [[maybe_unused]] const std::vector& arr, const std::string& label ) +{ + // Normalize: ensure 'array:' prefix for any std::vector steerable labels (user shouldn't add it) + const std::string alabel = (label.rfind("array:", 0) == 0) ? label : std::string("array:") + label; + catalystInfo_m << level4 << "::Initialize()::InitSteerChannel(" << alabel << "): | Type: std::vector size=" << arr.size() << endl; + // Derive namespace from canonical label of the form "array:." + std::string ns = alabel; + if (ns.rfind("array:", 0) == 0) ns = ns.substr(6); + auto dp = ns.find('.'); + if (dp != std::string::npos) ns = ns.substr(0, dp); + + // Register this array label with ProxyWriter so per-namespace array proxies are generated. + // The incoming label is expected to be of the form "array:." (constructed in RegisterStructMembers). + // Scalar-like element arrays use the generic include() path; ProxyWriter will parse and mark as array. + if constexpr (std::is_arithmetic_v>) { + // Use the first element as a default if available; otherwise 0. + Elem def{}; + if (!arr.empty()) def = arr.front(); + proxyWriter_m.include(def, alabel); + } else if constexpr (std::is_enum_v>) { + // Enum arrays: register as enum with choices; default from first element if present + int def = !arr.empty() ? static_cast(arr.front()) : 0; + // Prefer label-scoped choices, then type-scoped, else fallback to a checkbox + auto it = enumChoices_m.find(alabel); + if (it != enumChoices_m.end()) { + proxyWriter_m.includeEnum(alabel, it->second, def); + } else { + auto itt = enumChoicesByType_m.find(std::type_index(typeid(Elem))); + if (itt != enumChoicesByType_m.end()) { + proxyWriter_m.includeEnum(alabel, itt->second, def); + } else { + // No choices registered; use a checkbox as a minimal fallback UI + proxyWriter_m.includeBool(alabel, false); + } + } + } else if constexpr (std::is_same_v, bool>) { + bool def = !arr.empty() ? static_cast(arr.front()) : false; + proxyWriter_m.includeBool(alabel, def); + } else if constexpr (std::is_same_v, ippl::Button>) { + proxyWriter_m.includeButton(alabel); + } + + // Inform proxy writer about desired initial size for this array namespace + proxyWriter_m.setArrayInitialSize(ns, arr.size()); + + // Inform Python pipeline about each label (still needed for backward mapping) + conduit_cpp::Node scriptArgs = node_m["catalyst/scripts/script/args"]; + scriptArgs.append().set_string(alabel); + +} + +// Init: std::vector> steerables +template +void CatalystAdaptor::InitSteerChannel( [[maybe_unused]] const std::vector>& arr, const std::string& label ) +{ + // Normalize: ensure 'array:' prefix for any std::vector steerable labels + const std::string alabel = (label.rfind("array:", 0) == 0) ? label : std::string("array:") + label; + catalystInfo_m << level4 << "::Initialize()::InitSteerChannel(" << alabel << "): | Type: std::vector> size=" << arr.size() << endl; + // Derive namespace from canonical label of the form "array:." + std::string ns = alabel; + if (ns.rfind("array:", 0) == 0) ns = ns.substr(6); + auto dp = ns.find('.'); + if (dp != std::string::npos) ns = ns.substr(0, dp); + + // Register the vector array label with ProxyWriter (parsed as array namespace) + proxyWriter_m.includeVector(alabel); + + // Inform proxy writer about desired initial size for this array namespace + proxyWriter_m.setArrayInitialSize(ns, arr.size()); + + conduit_cpp::Node scriptArgs = node_m["catalyst/scripts/script/args"]; + scriptArgs.append().set_string(alabel); + +} + + + +// ===================================================================================== +// SETTING UP / FORWARDING CONDUIT NODE: +// ===================================================================================== + +// basic integer and double scalar overload +template +requires (!std::is_enum_v>) +void CatalystAdaptor::ForwardSteerChannel( const T& steerableScalarForwardpass, const std::string& steerableSuffix ) +{ + catalystInfo_m << level4 << "::Execute()::ForwardSteerChannel(" << steerableSuffix << "); | Type: " << typeid(T).name() << endl; + + auto steerableChannel = node_m["catalyst/channels/steerable_channel_0D_mesh"]; + + steerableChannel["type"].set("mesh"); + auto steerableData = steerableChannel["data"]; + steerableData["coordsets/coords/type"].set_string("explicit"); + steerableData["coordsets/coords/values/x"].set( 0 ); + + steerableData["topologies/sMesh_topo/type"].set("unstructured"); + steerableData["topologies/sMesh_topo/coordset"].set("coords"); + steerableData["topologies/sMesh_topo/elements/shape"].set("point"); + steerableData["topologies/sMesh_topo/elements/connectivity"].set( 0 ); + + + conduit_cpp::Node steerableField = steerableData["fields/steerable_field_f_" + steerableSuffix]; + steerableField["association"].set("vertex"); + steerableField["topology"].set("sMesh_topo"); + steerableField["volume_dependent"].set("false"); + + conduit_cpp::Node values = steerableField["values"]; + + // if constexpr(std::is_enum_v>){ + // values.set(static_cast(steerableScalarForwardpass)); + // } + // else + if constexpr(std::is_scalar_v){ + values.set(steerableScalarForwardpass); + } + else { + throw IpplException("Stream::InSitu::CatalystAdaptor::ForwardSteerChannel", "Unsupported steerable type for channel: " + steerableSuffix); + } +} + +// Enum overload (explicitfor clarity) +template +requires (std::is_enum_v>) +void CatalystAdaptor::ForwardSteerChannel( const E& e, const std::string& steerableSuffix ) +{ + catalystInfo_m << level4 << "::Execute()::ForwardSteerChannel(" << steerableSuffix << "); | Type: Enum" << endl; + auto steerableChannel = node_m["catalyst/channels/steerable_channel_0D_mesh"]; + steerableChannel["type"].set("mesh"); + auto steerableData = steerableChannel["data"]; + steerableData["coordsets/coords/type"].set_string("explicit"); + steerableData["coordsets/coords/values/x"].set( 0 ); + steerableData["topologies/sMesh_topo/type"].set("unstructured"); + steerableData["topologies/sMesh_topo/coordset"].set("coords"); + steerableData["topologies/sMesh_topo/elements/shape"].set("point"); + steerableData["topologies/sMesh_topo/elements/connectivity"].set( 0 ); + auto steerableField = steerableData["fields/steerable_field_f_" + steerableSuffix]; + steerableField["association"].set("vertex"); + steerableField["topology"].set("sMesh_topo"); + steerableField["volume_dependent"].set("false"); + steerableField["values"].set(static_cast(e)); +} + +// Bool-like Switch overload: forward as single scalar (0/1) +void CatalystAdaptor::ForwardSteerChannel( const bool& sw, const std::string& steerableSuffix ) +{ + catalystInfo_m << level4 << "::Execute()::ForwardSteerChannel(" << steerableSuffix << "); | Type: bool/Switch" << endl; + + auto steerableChannel = node_m["catalyst/channels/steerable_channel_0D_mesh"]; + steerableChannel["type"].set("mesh"); + auto steerableData = steerableChannel["data"]; + steerableData["coordsets/coords/type"].set_string("explicit"); + steerableData["coordsets/coords/values/x"].set( 0 ); + + steerableData["topologies/sMesh_topo/type"].set("unstructured"); + steerableData["topologies/sMesh_topo/coordset"].set("coords"); + steerableData["topologies/sMesh_topo/elements/shape"].set("point"); + steerableData["topologies/sMesh_topo/elements/connectivity"].set( 0 ); + + conduit_cpp::Node steerableField = steerableData["fields/steerable_field_f_" + steerableSuffix]; + steerableField["association"].set("vertex"); + steerableField["topology"].set("sMesh_topo"); + steerableField["volume_dependent"].set("false"); + + int boolAsInt = sw ? 1 : 0; + steerableField["values"].set(boolAsInt); +} + +// Bool-like Button overload: forward as single scalar (0/1) +void CatalystAdaptor::ForwardSteerChannel( const ippl::Button& btn, const std::string& steerableSuffix ) +{ + catalystInfo_m << level4 << "::Execute()::ForwardSteerChannel(" << steerableSuffix << "); | Type: Button" << endl; + auto steerableChannel = node_m["catalyst/channels/steerable_channel_0D_mesh"]; + steerableChannel["type"].set("mesh"); + auto steerableData = steerableChannel["data"]; + steerableData["coordsets/coords/type"].set_string("explicit"); + steerableData["coordsets/coords/values/x"].set( 0 ); + + steerableData["topologies/sMesh_topo/type"].set("unstructured"); + steerableData["topologies/sMesh_topo/coordset"].set("coords"); + steerableData["topologies/sMesh_topo/elements/shape"].set("point"); + steerableData["topologies/sMesh_topo/elements/connectivity"].set( 0 ); + + conduit_cpp::Node steerableField = steerableData["fields/steerable_field_f_" + steerableSuffix]; + steerableField["association"].set("vertex"); + steerableField["topology"].set("sMesh_topo"); + steerableField["volume_dependent"].set("false"); + steerableField["values"].set(static_cast(btn ? 1 : 0)); +} + +// ippl::Vector overload: forward single scalar conduit field with 3 components in "values" +////////////////////////////////////////////////// +// Note: +// I think x/y/z and 0/1/2 are valid syntax to send +// multicomponents forward, but in when returned in +// results the usual syntax is 0/1/2. +////////////////////////////////////////////////// +template +void CatalystAdaptor::ForwardSteerChannel( const ippl::Vector& steerableVecForwardpass, const std::string& steerableSuffix ) +{ + catalystInfo_m << level4 << "::Execute()::ForwardSteerChannel(" << steerableSuffix << "); | Vector<" << typeid(T).name() << "," << Dim_v << ">" << endl; + + auto steerableChannel = node_m["catalyst/channels/steerable_channel_0D_mesh"]; + steerableChannel["type"].set("mesh"); + auto steerableData = steerableChannel["data"]; + steerableData["coordsets/coords/type"].set_string("explicit"); + steerableData["coordsets/coords/values/x"].set(0); + + steerableData["topologies/sMesh_topo/type"].set("unstructured"); + steerableData["topologies/sMesh_topo/coordset"].set("coords"); + steerableData["topologies/sMesh_topo/elements/shape"].set("point"); + steerableData["topologies/sMesh_topo/elements/connectivity"].set(0); + + const std::string base = std::string("fields/steerable_field_f_") + steerableSuffix; + auto fnode = steerableData[base]; + fnode["association"].set_string("vertex"); + fnode["topology"].set_string("sMesh_topo"); + fnode["volume_dependent"].set_string("false"); + + if constexpr (std::is_integral_v) { + const int vx = static_cast(steerableVecForwardpass[0]); + fnode["values/x"].set(vx); + if constexpr (Dim_v >= 2) { + const int vy = static_cast(steerableVecForwardpass[1]); + fnode["values/y"].set(vy); + } + if constexpr (Dim_v >= 3) { + const int vz = static_cast(steerableVecForwardpass[2]); + fnode["values/z"].set(vz); + } + } else { + const double vx = static_cast(steerableVecForwardpass[0]); + fnode["values/x"].set(vx); + if constexpr (Dim_v >= 2) { + const double vy = static_cast(steerableVecForwardpass[1]); + fnode["values/y"].set(vy); + } + if constexpr (Dim_v >= 3) { + const double vz = static_cast(steerableVecForwardpass[2]); + fnode["values/z"].set(vz); + } + } +} + + +// std::vector overload: forward, publish as 1D mesh array under fields/steerableField_f_