diff --git a/.github/workflows/sonarqube.yml b/.github/workflows/sonarqube.yml index 250dc705..0510cc2d 100644 --- a/.github/workflows/sonarqube.yml +++ b/.github/workflows/sonarqube.yml @@ -38,14 +38,32 @@ jobs: LLVM_PROFILE_FILE="biofvm.profraw" build/llvm-cov-debug/reactions-diffusion/biofvm/tests/reactions-diffusion.biofvm.tests LLVM_PROFILE_FILE="micromechanics.profraw" build/llvm-cov-debug/mechanics/micromechanics/tests/mechanics.micromechanics.tests LLVM_PROFILE_FILE="physicell.profraw" build/llvm-cov-debug/mechanics/physicell/tests/mechanics.physicell.tests + LLVM_PROFILE_FILE="physicell.openmp_solver.profraw" build/llvm-cov-debug/mechanics/physicell/kernels/openmp_solver/tests/mechanics.physicell.kernels.openmp_solver.tests LLVM_PROFILE_FILE="biofvm.openmp_solver.profraw" build/llvm-cov-debug/reactions-diffusion/biofvm/kernels/openmp_solver/tests/reactions-diffusion.biofvm.kernels.openmp_solver.tests LLVM_PROFILE_FILE="biofvm.thrust_solver.profraw" build/llvm-cov-debug/reactions-diffusion/biofvm/kernels/thrust_solver/tests/reactions-diffusion.biofvm.kernels.tbb_thrust_solver.tests - name: Merge counters from different runs and convert the result into an llvm-cov-compatible form run: | - llvm-profdata-18 merge -o merged.profdata common.profraw biofvm.profraw micromechanics.profraw physicell.profraw biofvm.openmp_solver.profraw biofvm.thrust_solver.profraw + llvm-profdata-18 merge -o merged.profdata \ + common.profraw \ + biofvm.profraw \ + micromechanics.profraw \ + physicell.profraw \ + physicell.openmp_solver.profraw \ + biofvm.openmp_solver.profraw \ + biofvm.thrust_solver.profraw - name: Collect the coverage run: | - llvm-cov-18 show --show-branches=count --instr-profile merged.profdata build/llvm-cov-debug/common/tests/common.tests -object build/llvm-cov-debug/reactions-diffusion/biofvm/tests/reactions-diffusion.biofvm.tests -object build/llvm-cov-debug/mechanics/micromechanics/tests/mechanics.micromechanics.tests -object build/llvm-cov-debug/reactions-diffusion/biofvm/kernels/openmp_solver/tests/reactions-diffusion.biofvm.kernels.openmp_solver.tests -object build/llvm-cov-debug/reactions-diffusion/biofvm/kernels/thrust_solver/tests/reactions-diffusion.biofvm.kernels.tbb_thrust_solver.tests -object build/llvm-cov-debug/mechanics/physicell/tests/mechanics.physicell.tests > coverage.txt + llvm-cov-18 show \ + --show-branches=count \ + --instr-profile merged.profdata \ + build/llvm-cov-debug/common/tests/common.tests \ + -object build/llvm-cov-debug/reactions-diffusion/biofvm/tests/reactions-diffusion.biofvm.tests \ + -object build/llvm-cov-debug/mechanics/micromechanics/tests/mechanics.micromechanics.tests \ + -object build/llvm-cov-debug/mechanics/physicell/tests/mechanics.physicell.tests \ + -object build/llvm-cov-debug/mechanics/physicell/kernels/openmp_solver/tests/mechanics.physicell.kernels.openmp_solver.tests \ + -object build/llvm-cov-debug/reactions-diffusion/biofvm/kernels/openmp_solver/tests/reactions-diffusion.biofvm.kernels.openmp_solver.tests \ + -object build/llvm-cov-debug/reactions-diffusion/biofvm/kernels/thrust_solver/tests/reactions-diffusion.biofvm.kernels.tbb_thrust_solver.tests \ + > coverage.txt - name: SonarQube Scan uses: SonarSource/sonarqube-scan-action@v6 env: @@ -56,7 +74,7 @@ jobs: --define sonar.organization=bsc-life --define sonar.cfamily.compile-commands=${{ env.BUILD_WRAPPER_OUT_DIR }}/compile_commands.json --define sonar.cfamily.llvm-cov.reportPath=coverage.txt - --define sonar.sources=phenotype/physicore/src,common/include,reactions-diffusion/biofvm/src,reactions-diffusion/biofvm/include,mechanics/micromechanics/src,mechanics/micromechanics/include,mechanics/physicell/src,mechanics/physicell/include,reactions-diffusion/biofvm/kernels/openmp_solver/src,reactions-diffusion/biofvm/kernels/openmp_solver/examples,reactions-diffusion/biofvm/kernels/openmp_solver/include,reactions-diffusion/biofvm/kernels/thrust_solver/src,reactions-diffusion/biofvm/kernels/thrust_solver/examples,reactions-diffusion/biofvm/kernels/thrust_solver/include - --define sonar.tests=common/tests,reactions-diffusion/biofvm/tests,mechanics/micromechanics/tests,mechanics/physicell/tests,reactions-diffusion/biofvm/kernels/openmp_solver/tests,reactions-diffusion/biofvm/kernels/thrust_solver/tests + --define sonar.sources=phenotype/physicore/src,common/include,reactions-diffusion/biofvm/src,reactions-diffusion/biofvm/include,mechanics/micromechanics/src,mechanics/micromechanics/include,mechanics/physicell/src,mechanics/physicell/include,mechanics/physicell/kernels/openmp_solver/src,mechanics/physicell/kernels/openmp_solver/include,reactions-diffusion/biofvm/kernels/openmp_solver/src,reactions-diffusion/biofvm/kernels/openmp_solver/examples,reactions-diffusion/biofvm/kernels/openmp_solver/include,reactions-diffusion/biofvm/kernels/thrust_solver/src,reactions-diffusion/biofvm/kernels/thrust_solver/examples,reactions-diffusion/biofvm/kernels/thrust_solver/include + --define sonar.tests=common/tests,reactions-diffusion/biofvm/tests,mechanics/micromechanics/tests,mechanics/physicell/tests,mechanics/physicell/kernels/openmp_solver/tests,reactions-diffusion/biofvm/kernels/openmp_solver/tests,reactions-diffusion/biofvm/kernels/thrust_solver/tests --define sonar.coverage.exclusions=reactions-diffusion/biofvm/kernels/openmp_solver/examples/*,reactions-diffusion/biofvm/kernels/thrust_solver/examples/* --define sonar.cpd.exclusions=reactions-diffusion/biofvm/kernels/openmp_solver/examples/*,reactions-diffusion/biofvm/kernels/thrust_solver/examples/* diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index ecf680fb..4bfb1e01 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -1,9 +1,17 @@ -add_library(common INTERFACE) +add_library(common OBJECT) + +if(${BUILD_SHARED_LIBS}) + set_target_properties(common PROPERTIES POSITION_INDEPENDENT_CODE ON) +endif() file(GLOB PUBLIC_INC "include/common/*.h") +file(GLOB SOURCES "src/*.cpp") + target_sources( - common PUBLIC FILE_SET HEADERS BASE_DIRS - "${CMAKE_CURRENT_SOURCE_DIR}/include" FILES ${PUBLIC_INC}) + common + PUBLIC FILE_SET HEADERS BASE_DIRS "${CMAKE_CURRENT_SOURCE_DIR}/include" FILES + ${PUBLIC_INC} + PRIVATE ${SOURCES}) if(PHYSICORE_BUILD_TESTS) add_subdirectory(tests) diff --git a/common/include/common/cartesian_mesh.h b/common/include/common/cartesian_mesh.h new file mode 100644 index 00000000..2a18367b --- /dev/null +++ b/common/include/common/cartesian_mesh.h @@ -0,0 +1,62 @@ +#pragma once + +#include +#include + +#include "types.h" + +namespace physicore { + +/** + * @brief Uniform Cartesian mesh for spatial domain discretization. + * + * Provides axis-aligned structured grid with uniform spacing. + * Used by both reaction-diffusion and mechanics modules for domain partitioning. + * + * Supports 1D, 2D, and 3D domains with configurable voxel sizes and bounding boxes. + * Enables efficient spatial queries: position-to-voxel mapping and voxel linearization. + * + * @see reactions-diffusion/biofvm for diffusion solver usage + * @see mechanics/physicell for mechanics engine usage + */ +class cartesian_mesh +{ +public: + /// @brief Constructor + /// @param dims Number of spatial dimensions (1, 2, or 3) + /// @param bounding_box_mins Minimum coordinates of domain bounds + /// @param bounding_box_maxs Maximum coordinates of domain bounds + /// @param voxel_shape Size of each voxel in each dimension + explicit cartesian_mesh(index_t dims, std::array bounding_box_mins, + std::array bounding_box_maxs, std::array voxel_shape); + + /// @brief Get total number of voxels in the mesh + [[nodiscard]] std::size_t voxel_count() const; + + /// @brief Get volume of a single voxel + [[nodiscard]] index_t voxel_volume() const; + + /// @brief Find voxel indices containing the given position + /// @param position Spatial coordinates (must match dims) + /// @return Voxel indices in x, y, z order + [[nodiscard]] std::array voxel_position(std::span position) const; + + /// @brief Get the center position of a voxel + /// @param position Voxel indices + /// @return Center coordinates in x, y, z order + [[nodiscard]] std::array voxel_center(std::array position) const; + + /// @brief Convert 3D voxel indices to linear index + /// @param x, y, z Voxel indices + /// @return Linear index for use in flat arrays + [[nodiscard]] std::size_t linearize(index_t x, index_t y, index_t z) const; + + // Member data (public for direct access, following biofvm pattern) + index_t dims; ///< Number of spatial dimensions + std::array bounding_box_mins; ///< Minimum domain coordinates + std::array bounding_box_maxs; ///< Maximum domain coordinates + std::array voxel_shape; ///< Size of each voxel + std::array grid_shape; ///< Number of voxels per dimension +}; + +} // namespace physicore diff --git a/common/include/common/macros.h b/common/include/common/macros.h new file mode 100644 index 00000000..31be2b22 --- /dev/null +++ b/common/include/common/macros.h @@ -0,0 +1,13 @@ +#pragma once + +namespace physicore { + +#if defined(_MSC_VER) + #define PHYSICORE_RESTRICT __restrict +#elif defined(__GNUC__) || defined(__clang__) + #define PHYSICORE_RESTRICT __restrict__ +#else + #define PHYSICORE_RESTRICT +#endif + +} // namespace physicore diff --git a/common/include/common/registry.h b/common/include/common/registry.h new file mode 100644 index 00000000..14a61d16 --- /dev/null +++ b/common/include/common/registry.h @@ -0,0 +1,54 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace physicore { +template +class factory_registry +{ +public: + using ptr_t = std::unique_ptr; + using factory_func_t = std::function; + using map_t = std::unordered_map; + + bool register_factory(std::string name, factory_func_t&& f); + ptr_t get(const std::string& name); + +private: + map_t factories_; +}; + +template +bool factory_registry::register_factory(std::string name, factory_func_t&& f) +{ + auto [it, emplaced] = factories_.try_emplace(std::move(name), std::move(f)); + return emplaced; +} + +template +typename factory_registry::ptr_t factory_registry::get(const std::string& name) +{ + auto it = factories_.find(name); + if (it == factories_.end()) + { + return nullptr; + } + + return it->second(); +} + +template +struct generic_registry_adder +{ + explicit generic_registry_adder(std::string name) + { + Registry::instance().register_factory(std::move(name), []() { return std::make_unique(); }); + } +}; + +} // namespace physicore diff --git a/common/include/common/types.h b/common/include/common/types.h index 1ef4e590..492ff856 100644 --- a/common/include/common/types.h +++ b/common/include/common/types.h @@ -8,4 +8,13 @@ using real_t = double; using index_t = std::uint64_t; using sindex_t = std::int64_t; + +#if defined(_MSC_VER) + #define PHYSICORE_RESTRICT __restrict +#elif defined(__GNUC__) || defined(__clang__) + #define PHYSICORE_RESTRICT __restrict__ +#else + #define PHYSICORE_RESTRICT +#endif + } // namespace physicore diff --git a/common/src/cartesian_mesh.cpp b/common/src/cartesian_mesh.cpp new file mode 100644 index 00000000..5721e2d1 --- /dev/null +++ b/common/src/cartesian_mesh.cpp @@ -0,0 +1,78 @@ +#include "common/cartesian_mesh.h" + +#include + +namespace physicore { + +cartesian_mesh::cartesian_mesh(index_t dims, std::array bounding_box_mins, + std::array bounding_box_maxs, std::array voxel_shape) + : dims(dims), + bounding_box_mins(bounding_box_mins), + bounding_box_maxs(bounding_box_maxs), + voxel_shape(voxel_shape), + grid_shape({ 1, 1, 1 }) +{ + if (dims >= 1) + { + grid_shape[0] = (bounding_box_maxs[0] - bounding_box_mins[0] + voxel_shape[0] - 1) / voxel_shape[0]; + } + if (dims >= 2) + { + grid_shape[1] = (bounding_box_maxs[1] - bounding_box_mins[1] + voxel_shape[1] - 1) / voxel_shape[1]; + } + if (dims >= 3) + { + grid_shape[2] = (bounding_box_maxs[2] - bounding_box_mins[2] + voxel_shape[2] - 1) / voxel_shape[2]; + } +} + +std::size_t cartesian_mesh::voxel_count() const +{ + return (std::size_t)grid_shape[0] * (std::size_t)grid_shape[1] * (std::size_t)grid_shape[2]; +} + +index_t cartesian_mesh::voxel_volume() const { return voxel_shape[0] * voxel_shape[1] * voxel_shape[2]; } + +std::array cartesian_mesh::voxel_position(std::span position) const +{ + assert(position.size() == (size_t)dims); + assert(position[0] <= bounding_box_maxs[0] && position[0] >= bounding_box_mins[0]); + if (dims >= 2) + assert(position[1] <= bounding_box_maxs[1] && position[1] >= bounding_box_mins[1]); + if (dims >= 3) + assert(position[2] <= bounding_box_maxs[2] && position[2] >= bounding_box_mins[2]); + + switch (position.size()) + { + case 1: + return { (index_t)((position[0] - (real_t)bounding_box_mins[0]) / (real_t)voxel_shape[0]), 0, 0 }; + case 2: + return { (index_t)((position[0] - (real_t)bounding_box_mins[0]) / (real_t)voxel_shape[0]), + (index_t)((position[1] - (real_t)bounding_box_mins[1]) / (real_t)voxel_shape[1]), 0 }; + case 3: + return { (index_t)((position[0] - (real_t)bounding_box_mins[0]) / (real_t)voxel_shape[0]), + (index_t)((position[1] - (real_t)bounding_box_mins[1]) / (real_t)voxel_shape[1]), + (index_t)((position[2] - (real_t)bounding_box_mins[2]) / (real_t)voxel_shape[2]) }; + default: + assert(false); // Should never reach here + return { 0, 0, 0 }; + } +} + +std::array cartesian_mesh::voxel_center(std::array position) const +{ + assert(position[0] < grid_shape[0]); + assert(position[1] < grid_shape[1]); + assert(position[2] < grid_shape[2]); + + return { (real_t)(position[0] * voxel_shape[0] + bounding_box_mins[0]) + ((real_t)voxel_shape[0] / (real_t)2.0), + (real_t)(position[1] * voxel_shape[1] + bounding_box_mins[1]) + ((real_t)voxel_shape[1] / (real_t)2.0), + (real_t)(position[2] * voxel_shape[2] + bounding_box_mins[2]) + ((real_t)voxel_shape[2] / (real_t)2.0) }; +} + +std::size_t cartesian_mesh::linearize(index_t x, index_t y, index_t z) const +{ + return x + y * grid_shape[0] + z * grid_shape[0] * grid_shape[1]; +} + +} // namespace physicore diff --git a/mechanics/physicell/CMakeLists.txt b/mechanics/physicell/CMakeLists.txt index 2a02e7c1..1eae4f60 100644 --- a/mechanics/physicell/CMakeLists.txt +++ b/mechanics/physicell/CMakeLists.txt @@ -2,12 +2,25 @@ add_library(mechanics.physicell) add_library(physicore::mechanics::physicell ALIAS mechanics.physicell) add_library(mechanics.physicell_internal_iface INTERFACE) +if(${BUILD_SHARED_LIBS}) + set_target_properties( + mechanics.physicell + PROPERTIES CXX_VISIBILITY_PRESET hidden + VISIBILITY_INLINES_HIDDEN ON + POSITION_INDEPENDENT_CODE ON + INSTALL_RPATH_USE_LINK_PATH TRUE) +endif() + file(GLOB PUBLIC_INC "include/physicell/*.h") target_sources( mechanics.physicell - PRIVATE src/environment.cpp src/config_reader.cpp src/agent_data.cpp - src/mechanical_agent.cpp src/vtk_serializer_base.cpp + PRIVATE src/environment.cpp + src/config_reader.cpp + src/agent_data.cpp + src/mechanical_agent.cpp + src/vtk_serializer_base.cpp src/vtk_agents_serializer.cpp + src/solver_registry.cpp PUBLIC FILE_SET HEADERS BASE_DIRS "${CMAKE_CURRENT_SOURCE_DIR}/include" FILES ${PUBLIC_INC}) @@ -28,6 +41,8 @@ target_link_libraries( PUBLIC common PRIVATE mechanics.physicell_internal_iface pugixml::pugixml VTK::IOXML) +add_subdirectory(kernels/openmp_solver) + if(PHYSICORE_BUILD_TESTS) add_subdirectory(tests) endif() diff --git a/mechanics/physicell/include/physicell/agent_data.h b/mechanics/physicell/include/physicell/agent_data.h index b9105246..94f2b7c0 100644 --- a/mechanics/physicell/include/physicell/agent_data.h +++ b/mechanics/physicell/include/physicell/agent_data.h @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -43,6 +44,7 @@ struct mechanics_properties */ struct motility_properties { + using direction_update_func = std::function; // TODO: it needs access to substrate concentrations // Migration parameters std::vector is_motile; std::vector persistence_time; @@ -51,6 +53,7 @@ struct motility_properties std::vector migration_bias; std::vector motility_vector; // dims per agent std::vector restrict_to_2d; + std::vector direction_update_funcs; // Chemotaxis parameters std::vector chemotaxis_index; diff --git a/mechanics/physicell/include/physicell/environment.h b/mechanics/physicell/include/physicell/environment.h index 3908bb97..500e9e52 100644 --- a/mechanics/physicell/include/physicell/environment.h +++ b/mechanics/physicell/include/physicell/environment.h @@ -1,26 +1,43 @@ #pragma once #include +#include +#include #include #include +#include "mechanical_agent_container.h" #include "serializer.h" +#include "solver.h" namespace physicore::mechanics::physicell { class environment : public timestep_executor { - real_t timestep; - public: - explicit environment(real_t timestep) : timestep(timestep) {} + environment(real_t timestep, index_t dims, index_t agent_types_count, index_t substrates_count); void run_single_timestep() override; void serialize_state(real_t current_time) override; + real_t mechanics_timestep; + bool automated_spring_adhesion = false; + bool virtual_wall_at_domain_edges = true; + serializer_ptr serializer; + solver_ptr solver; + + std::unique_ptr agents; + + // Mesh data for spatial queries + void set_mesh(cartesian_mesh mesh) { this->mesh_ = std::move(mesh); } + const cartesian_mesh& get_mesh() const; + bool has_mesh() const { return mesh_.has_value(); } + +private: + std::optional mesh_; }; } // namespace physicore::mechanics::physicell diff --git a/mechanics/physicell/include/physicell/physicell_export.h b/mechanics/physicell/include/physicell/physicell_export.h new file mode 100644 index 00000000..b4aee060 --- /dev/null +++ b/mechanics/physicell/include/physicell/physicell_export.h @@ -0,0 +1,7 @@ +#pragma once + +// Placeholder export macro for shared-library builds. +// If/when CMake generates an export header, it can replace this file. +#ifndef PHYSICELL_EXPORT + #define PHYSICELL_EXPORT +#endif diff --git a/mechanics/physicell/include/physicell/solver.h b/mechanics/physicell/include/physicell/solver.h new file mode 100644 index 00000000..91b1b361 --- /dev/null +++ b/mechanics/physicell/include/physicell/solver.h @@ -0,0 +1,24 @@ +#pragma once + +#include + +#include +#include + +namespace physicore::mechanics::physicell { + +class environment; + +class PHYSICELL_EXPORT solver +{ +public: + virtual void initialize(environment& e) = 0; + + virtual void solve(environment& e, index_t iterations) = 0; + + virtual ~solver() = default; +}; + +using solver_ptr = std::unique_ptr; + +} // namespace physicore::mechanics::physicell diff --git a/mechanics/physicell/include/physicell/solver_registry.h b/mechanics/physicell/include/physicell/solver_registry.h new file mode 100644 index 00000000..635f1398 --- /dev/null +++ b/mechanics/physicell/include/physicell/solver_registry.h @@ -0,0 +1,23 @@ +#pragma once + +#include + +#include "solver.h" + +namespace physicore::mechanics::physicell { + +class solver_registry : public physicore::factory_registry +{ +public: + using base_t = physicore::factory_registry; + + static solver_registry& instance(); +}; + +template +concept derived_from_solver = std::derived_from; + +template +using registry_adder = physicore::generic_registry_adder; + +} // namespace physicore::mechanics::physicell diff --git a/mechanics/physicell/kernels/.clang-tidy b/mechanics/physicell/kernels/.clang-tidy new file mode 100644 index 00000000..c67794d3 --- /dev/null +++ b/mechanics/physicell/kernels/.clang-tidy @@ -0,0 +1,6 @@ +InheritParentConfig: true +Checks: > + performance-*, + -performance-avoid-endl, + + -readability-function-cognitive-complexity, diff --git a/mechanics/physicell/kernels/openmp_solver/CMakeLists.txt b/mechanics/physicell/kernels/openmp_solver/CMakeLists.txt new file mode 100644 index 00000000..a422b80f --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/CMakeLists.txt @@ -0,0 +1,50 @@ +add_library(mechanics.physicell.kernels.openmp_solver_internal_iface INTERFACE) +add_library(mechanics.physicell.kernels.openmp_solver OBJECT) + +if(${BUILD_SHARED_LIBS}) + set_target_properties(mechanics.physicell.kernels.openmp_solver + PROPERTIES POSITION_INDEPENDENT_CODE ON) +endif() + +find_package(OpenMP 4) + +target_sources( + mechanics.physicell.kernels.openmp_solver + PRIVATE src/openmp_solver.cpp src/register_solver.cpp src/position_solver.cpp + src/common_solver.cpp src/random.cpp + PUBLIC FILE_SET + HEADERS + BASE_DIRS + "${CMAKE_CURRENT_SOURCE_DIR}/include" + FILES + include/physicell/openmp_solver/register_solver.h + include/physicell/openmp_solver/position_solver.h) + +target_link_libraries(mechanics.physicell.kernels.openmp_solver_internal_iface + INTERFACE mechanics.physicell_internal_iface common) + +target_include_directories( + mechanics.physicell.kernels.openmp_solver_internal_iface + INTERFACE $ + $ + $) + +if(OpenMP_CXX_FOUND) + target_link_libraries(mechanics.physicell.kernels.openmp_solver_internal_iface + INTERFACE OpenMP::OpenMP_CXX) +endif() + +target_link_libraries( + mechanics.physicell.kernels.openmp_solver + PRIVATE mechanics.physicell.kernels.openmp_solver_internal_iface) + +target_link_libraries(mechanics.physicell + PRIVATE mechanics.physicell.kernels.openmp_solver) + +if(PHYSICORE_BUILD_EXAMPLES AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/examples") + add_subdirectory(examples) +endif() + +if(PHYSICORE_BUILD_TESTS) + add_subdirectory(tests) +endif() diff --git a/mechanics/physicell/kernels/openmp_solver/include/physicell/openmp_solver/position_solver.h b/mechanics/physicell/kernels/openmp_solver/include/physicell/openmp_solver/position_solver.h new file mode 100644 index 00000000..e857adbc --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/include/physicell/openmp_solver/position_solver.h @@ -0,0 +1,29 @@ +#pragma once + +#include + +#include +#include +#include +#include + + +namespace physicore::mechanics::physicell::kernels::openmp_solver { + +class position_solver : private generic_agent_solver +{ +public: + void update_cell_forces(environment& e); + + void update_cell_neighbors(environment& e, const cartesian_mesh& mesh); + + void update_motility(environment& e); + + void update_basement_membrane_interactions(environment& e, const cartesian_mesh& mesh); + + void update_spring_attachments(environment& e); + + void update_positions(environment& e); +}; + +} // namespace physicore::mechanics::physicell::kernels::openmp_solver diff --git a/mechanics/physicell/kernels/openmp_solver/include/physicell/openmp_solver/register_solver.h b/mechanics/physicell/kernels/openmp_solver/include/physicell/openmp_solver/register_solver.h new file mode 100644 index 00000000..cef1d260 --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/include/physicell/openmp_solver/register_solver.h @@ -0,0 +1,7 @@ +#pragma once + +namespace physicore::mechanics::physicell::kernels::openmp_solver { + +void attach_to_registry(); + +} // namespace physicore::mechanics::physicell::kernels::openmp_solver diff --git a/mechanics/physicell/kernels/openmp_solver/src/common_solver.cpp b/mechanics/physicell/kernels/openmp_solver/src/common_solver.cpp new file mode 100644 index 00000000..c774c452 --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/src/common_solver.cpp @@ -0,0 +1,17 @@ +#include "common_solver.h" + +#include + +namespace physicore::mechanics::physicell::kernels::openmp_solver { + +std::size_t common_solver::get_mesh_index(const voxel_pos_t& position, const cartesian_mesh& mesh) +{ + return mesh.linearize(position[0], position[1], position[2]); +} + +common_solver::voxel_pos_t common_solver::get_mesh_position(const real_t* position, const cartesian_mesh& mesh) +{ + return mesh.voxel_position(std::span(position, static_cast(mesh.dims))); +} + +} // namespace physicore::mechanics::physicell::kernels::openmp_solver diff --git a/mechanics/physicell/kernels/openmp_solver/src/common_solver.h b/mechanics/physicell/kernels/openmp_solver/src/common_solver.h new file mode 100644 index 00000000..4b7ba95b --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/src/common_solver.h @@ -0,0 +1,72 @@ +#pragma once + +#include +#include +#include +#include + +#include + +namespace physicore::mechanics::physicell::kernels::openmp_solver { + +class common_solver +{ +public: + using voxel_pos_t = std::array; + + static std::size_t get_mesh_index(const voxel_pos_t& position, const cartesian_mesh& mesh); + + static voxel_pos_t get_mesh_position(const real_t* position, const cartesian_mesh& mesh); + + template + static void for_each_in_mech_neighborhood_symmetric(const cartesian_mesh& mesh, + const std::vector>& cells_in_voxels, + const voxel_pos_t& position, index_t i, func_t&& f) + { + for_each_in_mech_neighborhood(mesh, cells_in_voxels, position, i, [&](index_t cell_idx) { + if (cell_idx > i) + f(cell_idx); + }); + } + + template + static void for_each_in_mech_neighborhood(const cartesian_mesh& mesh, + const std::vector>& cells_in_voxels, + const voxel_pos_t& position, index_t i, func_t&& f) + { + assert(cells_in_voxels.size() == mesh.voxel_count()); + + const sindex_t grid_x = static_cast(mesh.grid_shape[0]); + const sindex_t grid_y = static_cast(mesh.grid_shape[1]); + const sindex_t grid_z = static_cast(mesh.grid_shape[2]); + + for (sindex_t z = static_cast(position[2]) - 1; z <= static_cast(position[2]) + 1; z++) + { + if (z < 0 || z >= grid_z) + continue; + + for (sindex_t y = static_cast(position[1]) - 1; y <= static_cast(position[1]) + 1; y++) + { + if (y < 0 || y >= grid_y) + continue; + + for (sindex_t x = static_cast(position[0]) - 1; x <= static_cast(position[0]) + 1; + x++) + { + if (x < 0 || x >= grid_x) + continue; + + const auto voxel_idx = + mesh.linearize(static_cast(x), static_cast(y), static_cast(z)); + for (const auto cell_idx : cells_in_voxels[voxel_idx]) + { + if (i != cell_idx) + f(cell_idx); + } + } + } + } + } +}; + +} // namespace physicore::mechanics::physicell::kernels::openmp_solver diff --git a/mechanics/physicell/kernels/openmp_solver/src/openmp_solver.cpp b/mechanics/physicell/kernels/openmp_solver/src/openmp_solver.cpp new file mode 100644 index 00000000..6f539779 --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/src/openmp_solver.cpp @@ -0,0 +1,38 @@ +#include "openmp_solver.h" + +#include + +namespace physicore::mechanics::physicell::kernels::openmp_solver { + +void openmp_solver::initialize(environment& e) +{ + (void)e; + initialized = true; +} + +void openmp_solver::solve(environment& e, index_t iterations) +{ + if (!initialized) + { + initialize(e); + } + + position_solver mechanics_position_solver; + + for (index_t i = 0; i < iterations; ++i) + { + mechanics_position_solver.update_cell_neighbors(e, e.get_mesh()); + + mechanics_position_solver.update_cell_forces(e); + + mechanics_position_solver.update_motility(e); + + mechanics_position_solver.update_basement_membrane_interactions(e, e.get_mesh()); + + mechanics_position_solver.update_spring_attachments(e); + + mechanics_position_solver.update_positions(e); + } +} + +} // namespace physicore::mechanics::physicell::kernels::openmp_solver diff --git a/mechanics/physicell/kernels/openmp_solver/src/openmp_solver.h b/mechanics/physicell/kernels/openmp_solver/src/openmp_solver.h new file mode 100644 index 00000000..f830375f --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/src/openmp_solver.h @@ -0,0 +1,16 @@ +#pragma once + +#include + +namespace physicore::mechanics::physicell::kernels::openmp_solver { + +class openmp_solver : public solver +{ + bool initialized = false; + +public: + void initialize(environment& e) override; + void solve(environment& e, index_t iterations) override; +}; + +} // namespace physicore::mechanics::physicell::kernels::openmp_solver diff --git a/mechanics/physicell/kernels/openmp_solver/src/position_solver.cpp b/mechanics/physicell/kernels/openmp_solver/src/position_solver.cpp new file mode 100644 index 00000000..32b81647 --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/src/position_solver.cpp @@ -0,0 +1,571 @@ +#include "../include/physicell/openmp_solver/position_solver.h" + +#include +#include +#include +#include + +#include +#include +#include + +#include "common_solver.h" +#include "random.h" +#include "solver_helper.h" + +using namespace physicore::mechanics::physicell; + +namespace physicore::mechanics::physicell::kernels::openmp_solver { +constexpr real_t simple_pressure_coefficient = 36.64504274775163; // 1 / (12 * (1 - sqrt(pi/(2*sqrt(3))))^2) + +namespace { + +void clear_simple_pressure(real_t* PHYSICORE_RESTRICT simple_pressure, index_t count) +{ +#pragma omp for + for (index_t i = 0; i < count; i++) + { + simple_pressure[i] = 0; + } +} + +template +void solve_pair(index_t lhs, index_t rhs, index_t cell_defs_count, real_t* PHYSICORE_RESTRICT velocity, + real_t* PHYSICORE_RESTRICT simple_pressure, const real_t* PHYSICORE_RESTRICT position, + const real_t* PHYSICORE_RESTRICT radius, const real_t* PHYSICORE_RESTRICT cell_cell_repulsion_strength, + const real_t* PHYSICORE_RESTRICT cell_cell_adhesion_strength, + const real_t* PHYSICORE_RESTRICT relative_maximum_adhesion_distance, + const real_t* PHYSICORE_RESTRICT cell_adhesion_affinity, + const index_t* PHYSICORE_RESTRICT cell_definition_index) +{ + std::array position_difference {}; + + const real_t distance = + std::max(position_helper::difference_and_distance(position + lhs * dims, position + rhs * dims, + position_difference.data()), + 0.00001); + + // compute repulsion + real_t repulsion {}; + { + const real_t repulsive_distance = radius[lhs] + radius[rhs]; + + repulsion = 1 - distance / repulsive_distance; + + repulsion = repulsion < 0 ? 0 : repulsion; + + repulsion *= repulsion; + + // update simple pressure + simple_pressure[lhs] += repulsion * simple_pressure_coefficient; + simple_pressure[rhs] += repulsion * simple_pressure_coefficient; + + repulsion *= std::sqrt(cell_cell_repulsion_strength[lhs] * cell_cell_repulsion_strength[rhs]); + } + + // compute adhesion + real_t adhesion {}; + { + const real_t adhesion_distance = relative_maximum_adhesion_distance[lhs] * radius[lhs] + + relative_maximum_adhesion_distance[rhs] * radius[rhs]; + + adhesion = 1 - distance / adhesion_distance; + + adhesion = adhesion < 0 ? 0 : adhesion; + + adhesion *= adhesion; + + const index_t lhs_cell_def_index = cell_definition_index[lhs]; + const index_t rhs_cell_def_index = cell_definition_index[rhs]; + + adhesion *= std::sqrt(cell_cell_adhesion_strength[lhs] * cell_cell_adhesion_strength[rhs] + * cell_adhesion_affinity[lhs * cell_defs_count + rhs_cell_def_index] + * cell_adhesion_affinity[rhs * cell_defs_count + lhs_cell_def_index]); + } + + real_t const force = (repulsion - adhesion) / distance; + + position_helper::update_velocity(velocity + lhs * dims, position_difference.data(), force); +} + +template +void update_cell_forces_single(index_t i, index_t cell_def_count, real_t* PHYSICORE_RESTRICT velocity, + real_t* PHYSICORE_RESTRICT simple_pressure, const real_t* PHYSICORE_RESTRICT position, + const real_t* PHYSICORE_RESTRICT radius, + const real_t* PHYSICORE_RESTRICT cell_cell_repulsion_strength, + const real_t* PHYSICORE_RESTRICT cell_cell_adhesion_strength, + const real_t* PHYSICORE_RESTRICT relative_maximum_adhesion_distance, + const index_t* PHYSICORE_RESTRICT cell_definition_index, + const real_t* PHYSICORE_RESTRICT cell_adhesion_affinities, + std::vector* PHYSICORE_RESTRICT neighbors) +{ + for (const index_t j : neighbors[i]) + { + solve_pair(i, j, cell_def_count, velocity, simple_pressure, position, radius, + cell_cell_repulsion_strength, cell_cell_adhesion_strength, relative_maximum_adhesion_distance, + cell_adhesion_affinities, cell_definition_index); + } +} + +template +void update_cell_forces_internal(index_t agents_count, index_t cell_def_count, real_t* PHYSICORE_RESTRICT velocity, + real_t* PHYSICORE_RESTRICT simple_pressure, const real_t* PHYSICORE_RESTRICT position, + const real_t* PHYSICORE_RESTRICT radius, + const real_t* PHYSICORE_RESTRICT cell_cell_repulsion_strength, + const real_t* PHYSICORE_RESTRICT cell_cell_adhesion_strength, + const real_t* PHYSICORE_RESTRICT relative_maximum_adhesion_distance, + const index_t* PHYSICORE_RESTRICT cell_definition_index, + const real_t* PHYSICORE_RESTRICT cell_adhesion_affinities, + const std::uint8_t* PHYSICORE_RESTRICT is_movable, + std::vector* PHYSICORE_RESTRICT neighbors) +{ +#pragma omp for + for (index_t i = 0; i < agents_count; i++) + { + if (is_movable[i] == 0) + continue; + + update_cell_forces_single(i, cell_def_count, velocity, simple_pressure, position, radius, + cell_cell_repulsion_strength, cell_cell_adhesion_strength, + relative_maximum_adhesion_distance, cell_definition_index, + cell_adhesion_affinities, neighbors); + } +} + +template +void update_cell_neighbors_single(environment& e, index_t i, const real_t* PHYSICORE_RESTRICT position, + const real_t* PHYSICORE_RESTRICT radius, + const real_t* PHYSICORE_RESTRICT relative_maximum_adhesion_distance, + std::vector* PHYSICORE_RESTRICT neighbors, const cartesian_mesh& mesh, + const std::vector>& cells_in_voxels) +{ + (void)e; + common_solver::for_each_in_mech_neighborhood( + mesh, cells_in_voxels, common_solver::get_mesh_position(position + dims * i, mesh), i, [=](index_t j) { + const real_t adhesion_distance = + relative_maximum_adhesion_distance[i] * radius[i] + relative_maximum_adhesion_distance[j] * radius[j]; + + const real_t distance = position_helper::distance(position + i * dims, position + j * dims); + + if (distance <= adhesion_distance) + { + neighbors[i].push_back(j); + } + }); +} + +template +void update_cell_neighbors_internal(environment& e, index_t agents_count, const real_t* PHYSICORE_RESTRICT position, + const real_t* PHYSICORE_RESTRICT radius, + const real_t* PHYSICORE_RESTRICT relative_maximum_adhesion_distance, + const std::uint8_t* PHYSICORE_RESTRICT is_movable, + std::vector* PHYSICORE_RESTRICT neighbors, const cartesian_mesh& mesh, + const std::vector>& cells_in_voxels) +{ +#pragma omp for + for (index_t i = 0; i < agents_count; i++) + { + if (is_movable[i] == 0) + continue; + + update_cell_neighbors_single(e, i, position, radius, relative_maximum_adhesion_distance, neighbors, mesh, + cells_in_voxels); + } +} + +template +void update_motility_single( + index_t i, real_t time_step, real_t* PHYSICORE_RESTRICT motility_vector, real_t* PHYSICORE_RESTRICT velocity, + const real_t* PHYSICORE_RESTRICT persistence_time, const real_t* PHYSICORE_RESTRICT migration_bias, + real_t* PHYSICORE_RESTRICT migration_bias_direction, const std::uint8_t* PHYSICORE_RESTRICT restrict_to_2d, + const std::uint8_t* PHYSICORE_RESTRICT is_motile, const real_t* PHYSICORE_RESTRICT migration_speed, + const motility_properties::direction_update_func* PHYSICORE_RESTRICT update_migration_bias_direction_f, + const index_t* PHYSICORE_RESTRICT cell_definition_index) +{ + if (is_motile[i] == 0) + return; + + if (physicore::random::uniform() < time_step / persistence_time[i]) + { + std::array random_walk {}; + + position_helper::random_walk(restrict_to_2d, random_walk.data()); + + if (update_migration_bias_direction_f != nullptr && update_migration_bias_direction_f[i]) + { + update_migration_bias_direction_f[i](cell_definition_index[i]); + } + + position_helper::update_motility_vector(motility_vector + i * dims, random_walk.data(), + migration_bias_direction + i * dims, migration_bias[i]); + + position_helper::normalize_and_scale(motility_vector + i * dims, migration_speed[i]); + } + + position_helper::add(velocity + i * dims, motility_vector + i * dims); +} + +template +void update_motility_internal( + index_t agents_count, real_t time_step, real_t* PHYSICORE_RESTRICT motility_vector, + real_t* PHYSICORE_RESTRICT velocity, const real_t* PHYSICORE_RESTRICT persistence_time, + const real_t* PHYSICORE_RESTRICT migration_bias, real_t* PHYSICORE_RESTRICT migration_bias_direction, + const std::uint8_t* PHYSICORE_RESTRICT restrict_to_2d, const std::uint8_t* PHYSICORE_RESTRICT is_motile, + const real_t* PHYSICORE_RESTRICT migration_speed, + const motility_properties::direction_update_func* PHYSICORE_RESTRICT update_migration_bias_direction_f, + const index_t* PHYSICORE_RESTRICT cell_definition_index) +{ +#pragma omp for + for (index_t i = 0; i < agents_count; i++) + { + update_motility_single(i, time_step, motility_vector, velocity, persistence_time, migration_bias, + migration_bias_direction, restrict_to_2d, is_motile, migration_speed, + update_migration_bias_direction_f, cell_definition_index); + } +} + +template +void update_basement_membrane_interactions_single(index_t i, real_t* PHYSICORE_RESTRICT velocity, + const real_t* PHYSICORE_RESTRICT position, + const real_t* PHYSICORE_RESTRICT radius, + const real_t* PHYSICORE_RESTRICT cell_BM_repulsion_strength, + const cartesian_mesh& mesh) +{ + position_helper::update_membrane_velocities(velocity + i * dims, position + i * dims, mesh, radius[i], + cell_BM_repulsion_strength[i]); +} + +template +void update_basement_membrane_interactions_internal(index_t agents_count, real_t* PHYSICORE_RESTRICT velocity, + const real_t* PHYSICORE_RESTRICT position, + const real_t* PHYSICORE_RESTRICT radius, + const real_t* PHYSICORE_RESTRICT cell_BM_repulsion_strength, + const std::uint8_t* PHYSICORE_RESTRICT is_movable, + const cartesian_mesh& mesh) +{ +#pragma omp for + for (index_t i = 0; i < agents_count; i++) + { + if (is_movable[i] == 0) + continue; + + update_basement_membrane_interactions_single(i, velocity, position, radius, cell_BM_repulsion_strength, + mesh); + } +} + +template +void spring_contract_function(index_t agents_count, index_t cell_defs_count, real_t* PHYSICORE_RESTRICT velocity, + const index_t* PHYSICORE_RESTRICT cell_definition_index, + const real_t* PHYSICORE_RESTRICT attachment_elastic_constant, + const real_t* PHYSICORE_RESTRICT cell_adhesion_affinity, + const real_t* PHYSICORE_RESTRICT position, + const std::uint8_t* PHYSICORE_RESTRICT is_movable, + std::vector* PHYSICORE_RESTRICT springs) +{ +#pragma omp for + for (index_t this_cell_index = 0; this_cell_index < agents_count; this_cell_index++) + { + if (is_movable[this_cell_index] == 0) + continue; + + for (std::size_t j = 0; j < springs[this_cell_index].size(); j++) + { + const index_t other_cell_index = springs[this_cell_index][j]; + + const index_t this_cell_def_index = cell_definition_index[this_cell_index]; + const index_t other_cell_def_index = cell_definition_index[other_cell_index]; + + const real_t adhesion = + sqrt(attachment_elastic_constant[this_cell_index] * attachment_elastic_constant[other_cell_index] + * cell_adhesion_affinity[this_cell_index * cell_defs_count + other_cell_def_index] + * cell_adhesion_affinity[other_cell_index * cell_defs_count + this_cell_def_index]); + + std::array difference {}; + + position_helper::subtract(difference.data(), position + other_cell_index * dims, + position + this_cell_index * dims); + + position_helper::update_velocity(velocity + this_cell_index * dims, difference.data(), adhesion); + } + } +} + +constexpr index_t erased_spring = -1; + +void update_spring_attachments_internal(index_t agents_count, real_t time_step, index_t cell_defs_count, + const real_t* PHYSICORE_RESTRICT detachment_rate, + const real_t* PHYSICORE_RESTRICT attachment_rate, + const real_t* PHYSICORE_RESTRICT cell_adhesion_affinities, + const index_t* PHYSICORE_RESTRICT maximum_number_of_attachments, + const index_t* PHYSICORE_RESTRICT cell_definition_index, + const std::vector* PHYSICORE_RESTRICT neighbors, + std::vector* PHYSICORE_RESTRICT springs) +{ +// mark springs for detachment +#pragma omp for + for (index_t this_cell_index = 0; this_cell_index < agents_count; this_cell_index++) + { + for (index_t j = 0; j < (index_t)springs[this_cell_index].size(); j++) + { + if (physicore::random::uniform() <= detachment_rate[this_cell_index] * time_step) + { +#pragma omp critical + { + const index_t other_cell_index = springs[this_cell_index][j]; + + if (other_cell_index != erased_spring) + { + springs[this_cell_index][j] = erased_spring; + + *std::ranges::find(springs[other_cell_index], this_cell_index) = erased_spring; + } + } + } + } + } + +// remove marked springs +#pragma omp for + for (index_t this_cell_index = 0; this_cell_index < agents_count; this_cell_index++) + { + auto removed = std::ranges::remove(springs[this_cell_index], erased_spring); + + springs[this_cell_index].erase(removed.begin(), removed.end()); + } + + // attach cells to springs + +#pragma omp for + for (index_t this_cell_index = 0; this_cell_index < agents_count; this_cell_index++) + { + for (std::size_t j = 0; j < neighbors[this_cell_index].size(); j++) + { + const index_t other_cell_index = neighbors[this_cell_index][j]; + + if (other_cell_index < this_cell_index) + continue; + + const real_t affinity_l = + cell_adhesion_affinities[this_cell_index * cell_defs_count + cell_definition_index[other_cell_index]]; + + const real_t attachment_prob_l = attachment_rate[this_cell_index] * time_step * affinity_l; + + const real_t affinity_r = + cell_adhesion_affinities[other_cell_index * cell_defs_count + cell_definition_index[this_cell_index]]; + + const real_t attachment_prob_r = attachment_rate[other_cell_index] * time_step * affinity_r; + + if (physicore::random::uniform() <= attachment_prob_l || physicore::random::uniform() <= attachment_prob_r) + { +#pragma omp critical + { + if ((index_t)springs[this_cell_index].size() < maximum_number_of_attachments[this_cell_index] + && (index_t)springs[other_cell_index].size() < maximum_number_of_attachments[other_cell_index]) + { + springs[this_cell_index].push_back(other_cell_index); + springs[other_cell_index].push_back(this_cell_index); + } + } + } + } + } +} + +template +void update_positions_internal(index_t agents_count, real_t time_step, real_t* PHYSICORE_RESTRICT position, + real_t* PHYSICORE_RESTRICT velocity, real_t* PHYSICORE_RESTRICT previous_velocity, + const std::uint8_t* PHYSICORE_RESTRICT is_movable) +{ +#pragma omp for + for (index_t i = 0; i < agents_count; i++) + { + if (!is_movable[i]) + continue; + + const real_t factor = time_step * 1.5; + const real_t previous_factor = time_step * -0.5; + + for (index_t d = 0; d < dims; d++) + { + position[i * dims + d] += + velocity[i * dims + d] * factor + previous_velocity[i * dims + d] * previous_factor; + + previous_velocity[i * dims + d] = velocity[i * dims + d]; + velocity[i * dims + d] = 0; + } + } +} + +} // namespace + +void position_solver::update_cell_forces(environment& e) +{ + auto& data = retrieve_agent_data(*e.agents); + const index_t dims = data.base_data.dims; + + clear_simple_pressure(data.state_data.simple_pressure.data(), data.agents_count); + + if (dims == 1) + update_cell_forces_internal<1>(data.agents_count, data.agent_types_count, data.velocity.data(), + data.state_data.simple_pressure.data(), data.base_data.positions.data(), + data.radius.data(), data.mechanics_data.cell_cell_repulsion_strength.data(), + data.mechanics_data.cell_cell_adhesion_strength.data(), + data.mechanics_data.relative_maximum_adhesion_distance.data(), + data.state_data.agent_type_index.data(), + data.mechanics_data.cell_adhesion_affinities.data(), + data.state_data.is_movable.data(), data.state_data.neighbors.data()); + else if (dims == 2) + update_cell_forces_internal<2>(data.agents_count, data.agent_types_count, data.velocity.data(), + data.state_data.simple_pressure.data(), data.base_data.positions.data(), + data.radius.data(), data.mechanics_data.cell_cell_repulsion_strength.data(), + data.mechanics_data.cell_cell_adhesion_strength.data(), + data.mechanics_data.relative_maximum_adhesion_distance.data(), + data.state_data.agent_type_index.data(), + data.mechanics_data.cell_adhesion_affinities.data(), + data.state_data.is_movable.data(), data.state_data.neighbors.data()); + else if (dims == 3) + update_cell_forces_internal<3>(data.agents_count, data.agent_types_count, data.velocity.data(), + data.state_data.simple_pressure.data(), data.base_data.positions.data(), + data.radius.data(), data.mechanics_data.cell_cell_repulsion_strength.data(), + data.mechanics_data.cell_cell_adhesion_strength.data(), + data.mechanics_data.relative_maximum_adhesion_distance.data(), + data.state_data.agent_type_index.data(), + data.mechanics_data.cell_adhesion_affinities.data(), + data.state_data.is_movable.data(), data.state_data.neighbors.data()); +} + +void position_solver::update_cell_neighbors(environment& e, const cartesian_mesh& mesh) +{ + auto& data = retrieve_agent_data(*e.agents); + const index_t dims = data.base_data.dims; + + std::vector> cells_in_voxels(mesh.voxel_count()); + for (index_t i = 0; i < data.agents_count; i++) + { + const auto voxel_pos = common_solver::get_mesh_position(data.base_data.positions.data() + dims * i, mesh); + const auto voxel_idx = common_solver::get_mesh_index(voxel_pos, mesh); + cells_in_voxels[voxel_idx].push_back(i); + } + + // clear neighbors +#pragma omp for + for (index_t i = 0; i < data.agents_count; i++) + data.state_data.neighbors[i].clear(); + + if (dims == 1) + update_cell_neighbors_internal<1>(e, data.agents_count, data.base_data.positions.data(), data.radius.data(), + data.mechanics_data.relative_maximum_adhesion_distance.data(), + data.state_data.is_movable.data(), data.state_data.neighbors.data(), mesh, + cells_in_voxels); + else if (dims == 2) + update_cell_neighbors_internal<2>(e, data.agents_count, data.base_data.positions.data(), data.radius.data(), + data.mechanics_data.relative_maximum_adhesion_distance.data(), + data.state_data.is_movable.data(), data.state_data.neighbors.data(), mesh, + cells_in_voxels); + else if (dims == 3) + update_cell_neighbors_internal<3>(e, data.agents_count, data.base_data.positions.data(), data.radius.data(), + data.mechanics_data.relative_maximum_adhesion_distance.data(), + data.state_data.is_movable.data(), data.state_data.neighbors.data(), mesh, + cells_in_voxels); +} + +void position_solver::update_motility(environment& e) +{ + auto& data = retrieve_agent_data(*e.agents); + + if (data.base_data.dims == 1) + update_motility_internal<1>( + data.agents_count, e.mechanics_timestep, data.motility_data.motility_vector.data(), data.velocity.data(), + data.motility_data.persistence_time.data(), data.motility_data.migration_bias.data(), + data.motility_data.migration_bias_direction.data(), data.motility_data.restrict_to_2d.data(), + data.motility_data.is_motile.data(), data.motility_data.migration_speed.data(), + data.motility_data.direction_update_funcs.data(), data.state_data.agent_type_index.data()); + else if (data.base_data.dims == 2) + update_motility_internal<2>( + data.agents_count, e.mechanics_timestep, data.motility_data.motility_vector.data(), data.velocity.data(), + data.motility_data.persistence_time.data(), data.motility_data.migration_bias.data(), + data.motility_data.migration_bias_direction.data(), data.motility_data.restrict_to_2d.data(), + data.motility_data.is_motile.data(), data.motility_data.migration_speed.data(), + data.motility_data.direction_update_funcs.data(), data.state_data.agent_type_index.data()); + else if (data.base_data.dims == 3) + update_motility_internal<3>( + data.agents_count, e.mechanics_timestep, data.motility_data.motility_vector.data(), data.velocity.data(), + data.motility_data.persistence_time.data(), data.motility_data.migration_bias.data(), + data.motility_data.migration_bias_direction.data(), data.motility_data.restrict_to_2d.data(), + data.motility_data.is_motile.data(), data.motility_data.migration_speed.data(), + data.motility_data.direction_update_funcs.data(), data.state_data.agent_type_index.data()); +} + +void position_solver::update_basement_membrane_interactions(environment& e, const cartesian_mesh& mesh) +{ + if (!e.virtual_wall_at_domain_edges) // note: where do we include this + return; + + auto& data = retrieve_agent_data(*e.agents); + + if (data.base_data.dims == 1) + update_basement_membrane_interactions_internal<1>( + data.agents_count, data.velocity.data(), data.base_data.positions.data(), data.radius.data(), + data.mechanics_data.cell_BM_repulsion_strength.data(), data.state_data.is_movable.data(), mesh); + else if (data.base_data.dims == 2) + update_basement_membrane_interactions_internal<2>( + data.agents_count, data.velocity.data(), data.base_data.positions.data(), data.radius.data(), + data.mechanics_data.cell_BM_repulsion_strength.data(), data.state_data.is_movable.data(), mesh); + else if (data.base_data.dims == 3) + update_basement_membrane_interactions_internal<3>( + data.agents_count, data.velocity.data(), data.base_data.positions.data(), data.radius.data(), + data.mechanics_data.cell_BM_repulsion_strength.data(), data.state_data.is_movable.data(), mesh); +} + +void position_solver::update_spring_attachments(environment& e) +{ + if (!e.automated_spring_adhesion) + return; + + auto& data = retrieve_agent_data(*e.agents); + const index_t dims = data.base_data.dims; + + update_spring_attachments_internal( + data.agents_count, e.mechanics_timestep, data.agent_types_count, data.mechanics_data.detachment_rate.data(), + data.mechanics_data.attachment_rate.data(), data.mechanics_data.cell_adhesion_affinities.data(), + data.mechanics_data.maximum_number_of_attachments.data(), data.state_data.agent_type_index.data(), + data.state_data.neighbors.data(), data.state_data.springs.data()); + + if (dims == 1) + spring_contract_function<1>( + data.agents_count, data.agent_types_count, data.velocity.data(), data.state_data.agent_type_index.data(), + data.mechanics_data.attachment_elastic_constant.data(), data.mechanics_data.cell_adhesion_affinities.data(), + data.base_data.positions.data(), data.state_data.is_movable.data(), data.state_data.springs.data()); + else if (dims == 2) + spring_contract_function<2>( + data.agents_count, data.agent_types_count, data.velocity.data(), data.state_data.agent_type_index.data(), + data.mechanics_data.attachment_elastic_constant.data(), data.mechanics_data.cell_adhesion_affinities.data(), + data.base_data.positions.data(), data.state_data.is_movable.data(), data.state_data.springs.data()); + else if (dims == 3) + spring_contract_function<3>( + data.agents_count, data.agent_types_count, data.velocity.data(), data.state_data.agent_type_index.data(), + data.mechanics_data.attachment_elastic_constant.data(), data.mechanics_data.cell_adhesion_affinities.data(), + data.base_data.positions.data(), data.state_data.is_movable.data(), data.state_data.springs.data()); +} + +void position_solver::update_positions(environment& e) +{ + auto& data = retrieve_agent_data(*e.agents); + const index_t dims = data.base_data.dims; + + if (dims == 1) + update_positions_internal<1>(data.agents_count, e.mechanics_timestep, data.base_data.positions.data(), + data.velocity.data(), data.previous_velocity.data(), + data.state_data.is_movable.data()); + else if (dims == 2) + update_positions_internal<2>(data.agents_count, e.mechanics_timestep, data.base_data.positions.data(), + data.velocity.data(), data.previous_velocity.data(), + data.state_data.is_movable.data()); + else if (dims == 3) + update_positions_internal<3>(data.agents_count, e.mechanics_timestep, data.base_data.positions.data(), + data.velocity.data(), data.previous_velocity.data(), + data.state_data.is_movable.data()); +} + +} // namespace physicore::mechanics::physicell::kernels::openmp_solver diff --git a/mechanics/physicell/kernels/openmp_solver/src/random.cpp b/mechanics/physicell/kernels/openmp_solver/src/random.cpp new file mode 100644 index 00000000..90e079bb --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/src/random.cpp @@ -0,0 +1,50 @@ +#ifdef _OPENMP + #include +#endif + +#include + +#include "random.h" + +namespace { +std::mt19937& thread_local_generator() +{ + static thread_local std::mt19937 generator; + return generator; +} +} // namespace + +physicore::real_t physicore::random::uniform(const physicore::real_t min, const physicore::real_t max) +{ + std::uniform_real_distribution distribution(min, max); + return distribution(thread_local_generator()); +} + +physicore::real_t physicore::random::normal(const physicore::real_t mean, const physicore::real_t std) +{ + std::normal_distribution distribution(mean, std); + return distribution(thread_local_generator()); +} + +void physicore::random::set_seed(unsigned int seed) +{ +#ifdef _OPENMP + std::vector initial_sequence(omp_get_num_threads()); + + for (int i = 0; i < omp_get_num_threads(); i++) + initial_sequence[i] = seed + i; + + std::seed_seq seq(initial_sequence.begin(), initial_sequence.end()); + + std::vector seeds(omp_get_num_threads()); + seq.generate(seeds.begin(), seeds.end()); + + #pragma omp parallel + { + int id = omp_get_thread_num(); + thread_local_generator().seed(seeds[id]); + } +#else + thread_local_generator().seed(seed); +#endif +} diff --git a/mechanics/physicell/kernels/openmp_solver/src/random.h b/mechanics/physicell/kernels/openmp_solver/src/random.h new file mode 100644 index 00000000..ca896584 --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/src/random.h @@ -0,0 +1,17 @@ +#pragma once + +#include + +namespace physicore { + +class random +{ +public: + static real_t uniform(real_t min = 0, real_t max = 1); + + static real_t normal(real_t mean = 0, real_t std = 1); + + static void set_seed(unsigned int seed); +}; + +} // namespace physicore diff --git a/mechanics/physicell/kernels/openmp_solver/src/register_solver.cpp b/mechanics/physicell/kernels/openmp_solver/src/register_solver.cpp new file mode 100644 index 00000000..6459171d --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/src/register_solver.cpp @@ -0,0 +1,10 @@ +#include "physicell/openmp_solver/register_solver.h" + +#include + +#include "openmp_solver.h" + +void physicore::mechanics::physicell::kernels::openmp_solver::attach_to_registry() +{ + static const physicore::mechanics::physicell::registry_adder openmp_solver_adder("openmp_solver"); +} diff --git a/mechanics/physicell/kernels/openmp_solver/src/solver_helper.h b/mechanics/physicell/kernels/openmp_solver/src/solver_helper.h new file mode 100644 index 00000000..303195f4 --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/src/solver_helper.h @@ -0,0 +1,294 @@ +#pragma once + +#include +#include + +#include "random.h" + +constexpr physicore::real_t zero_threshold = 1e-16; + +namespace physicore::mechanics::physicell::kernels::openmp_solver { + + +inline void update_membrane_velocity(real_t position, real_t bounding_box, real_t sign, real_t radius, + real_t repulsion_strength, real_t& velocity) +{ + real_t distance = bounding_box - position; + distance = (distance < 0) ? -distance : distance; + distance = (distance < 0.00001) ? 0.00001 : distance; + + // BITHACK VERSION: TO BE TESTED + // distance = (distance ^ (distance >> 63)) - (distance >> 63); + // distance += (0.00001 - distance) & -(distance < 0.00001); + + real_t repulsion = 1 - distance / radius; + repulsion = repulsion < 0 ? 0 : repulsion; + + repulsion *= repulsion * repulsion_strength * sign; + + velocity += repulsion * distance; +} + +template +struct position_helper +{}; + +template <> +struct position_helper<1> +{ + static inline real_t distance(const real_t* PHYSICORE_RESTRICT lhs, const real_t* PHYSICORE_RESTRICT rhs) + { + return std::abs(lhs[0] - rhs[0]); + } + + static inline real_t difference_and_distance(const real_t* PHYSICORE_RESTRICT lhs, + const real_t* PHYSICORE_RESTRICT rhs, + real_t* PHYSICORE_RESTRICT difference) + { + difference[0] = lhs[0] - rhs[0]; + + return std::abs(difference[0]); + } + + static inline void update_velocities(real_t* PHYSICORE_RESTRICT lhs, real_t* PHYSICORE_RESTRICT rhs, + const real_t* PHYSICORE_RESTRICT difference, const real_t force) + { + lhs[0] += force * difference[0]; + rhs[0] -= force * difference[0]; + } + + static inline void update_velocity(real_t* PHYSICORE_RESTRICT velocity, const real_t* PHYSICORE_RESTRICT difference, + const real_t force) + { + velocity[0] += force * difference[0]; + } + + static void random_walk(bool, real_t* PHYSICORE_RESTRICT walk) + { + real_t rand = random::uniform(); + walk[0] = rand < 0.5 ? -1 : 1; + } + + static inline void update_motility_vector(real_t* PHYSICORE_RESTRICT motility_vector, + const real_t* PHYSICORE_RESTRICT walk, + const real_t* PHYSICORE_RESTRICT migration_bias_direction, + const real_t migration_bias) + { + motility_vector[0] = (1 - migration_bias) * walk[0] + migration_bias * migration_bias_direction[0]; + } + + static inline void normalize_and_scale(real_t* PHYSICORE_RESTRICT vector, real_t scale) + { + real_t length = std::abs(vector[0]); + + vector[0] = length > zero_threshold ? vector[0] * scale / length : 0; + } + + static inline void update_membrane_velocities(real_t* PHYSICORE_RESTRICT velocity, + const real_t* PHYSICORE_RESTRICT position, const cartesian_mesh& mesh, + const real_t radius, const real_t repulsion_strength) + { + update_membrane_velocity(position[0], mesh.bounding_box_mins[0], 1, radius, repulsion_strength, velocity[0]); + update_membrane_velocity(position[0], mesh.bounding_box_maxs[0], -1, radius, repulsion_strength, velocity[0]); + } + + static inline void add(real_t* PHYSICORE_RESTRICT lhs, const real_t* PHYSICORE_RESTRICT rhs) { lhs[0] += rhs[0]; } + + static inline void subtract(real_t* PHYSICORE_RESTRICT dst, const real_t* PHYSICORE_RESTRICT lhs, + const real_t* PHYSICORE_RESTRICT rhs) + { + dst[0] = lhs[0] - rhs[0]; + } +}; + +template <> +struct position_helper<2> +{ + static inline real_t distance(const real_t* PHYSICORE_RESTRICT lhs, const real_t* PHYSICORE_RESTRICT rhs) + { + return std::sqrt((lhs[0] - rhs[0]) * (lhs[0] - rhs[0]) + (lhs[1] - rhs[1]) * (lhs[1] - rhs[1])); + } + + static inline real_t difference_and_distance(const real_t* PHYSICORE_RESTRICT lhs, + const real_t* PHYSICORE_RESTRICT rhs, + real_t* PHYSICORE_RESTRICT difference) + { + difference[0] = lhs[0] - rhs[0]; + difference[1] = lhs[1] - rhs[1]; + + return std::sqrt(difference[0] * difference[0] + difference[1] * difference[1]); + } + + static inline void update_velocities(real_t* PHYSICORE_RESTRICT lhs, real_t* PHYSICORE_RESTRICT rhs, + const real_t* PHYSICORE_RESTRICT difference, const real_t force) + { + lhs[0] += force * difference[0]; + lhs[1] += force * difference[1]; + + rhs[0] -= force * difference[0]; + rhs[1] -= force * difference[1]; + } + + static inline void update_velocity(real_t* PHYSICORE_RESTRICT velocity, const real_t* PHYSICORE_RESTRICT difference, + const real_t force) + { + velocity[0] += force * difference[0]; + velocity[1] += force * difference[1]; + } + + static void random_walk(bool, real_t* PHYSICORE_RESTRICT walk) + { + real_t theta = random::uniform(0, 2 * std::numbers::pi_v); + walk[0] = std::cos(theta); + walk[1] = std::sin(theta); + } + + static inline void update_motility_vector(real_t* PHYSICORE_RESTRICT motility_vector, + const real_t* PHYSICORE_RESTRICT walk, + const real_t* PHYSICORE_RESTRICT migration_bias_direction, + const real_t migration_bias) + { + motility_vector[0] = (1 - migration_bias) * walk[0] + migration_bias * migration_bias_direction[0]; + motility_vector[1] = (1 - migration_bias) * walk[1] + migration_bias * migration_bias_direction[1]; + } + + static inline void normalize_and_scale(real_t* PHYSICORE_RESTRICT vector, real_t scale) + { + real_t length = std::sqrt(vector[0] * vector[0] + vector[1] * vector[1]); + + vector[0] = length > zero_threshold ? vector[0] * scale / length : 0; + vector[1] = length > zero_threshold ? vector[1] * scale / length : 0; + } + + static inline void update_membrane_velocities(real_t* PHYSICORE_RESTRICT velocity, + const real_t* PHYSICORE_RESTRICT position, const cartesian_mesh& mesh, + const real_t radius, const real_t repulsion_strength) + { + update_membrane_velocity(position[0], mesh.bounding_box_mins[0], 1, radius, repulsion_strength, velocity[0]); + update_membrane_velocity(position[0], mesh.bounding_box_maxs[0], -1, radius, repulsion_strength, velocity[0]); + update_membrane_velocity(position[1], mesh.bounding_box_mins[1], 1, radius, repulsion_strength, velocity[1]); + update_membrane_velocity(position[1], mesh.bounding_box_maxs[1], -1, radius, repulsion_strength, velocity[1]); + } + + static inline void add(real_t* PHYSICORE_RESTRICT lhs, const real_t* PHYSICORE_RESTRICT rhs) + { + lhs[0] += rhs[0]; + lhs[1] += rhs[1]; + } + + static inline void subtract(real_t* PHYSICORE_RESTRICT dst, const real_t* PHYSICORE_RESTRICT lhs, + const real_t* PHYSICORE_RESTRICT rhs) + { + dst[0] = lhs[0] - rhs[0]; + dst[1] = lhs[1] - rhs[1]; + } +}; + +template <> +struct position_helper<3> +{ + static inline real_t distance(const real_t* PHYSICORE_RESTRICT lhs, const real_t* PHYSICORE_RESTRICT rhs) + { + return std::sqrt((lhs[0] - rhs[0]) * (lhs[0] - rhs[0]) + (lhs[1] - rhs[1]) * (lhs[1] - rhs[1]) + + (lhs[2] - rhs[2]) * (lhs[2] - rhs[2])); + } + + static inline real_t difference_and_distance(const real_t* PHYSICORE_RESTRICT lhs, + const real_t* PHYSICORE_RESTRICT rhs, + real_t* PHYSICORE_RESTRICT difference) + { + difference[0] = lhs[0] - rhs[0]; + difference[1] = lhs[1] - rhs[1]; + difference[2] = lhs[2] - rhs[2]; + + return std::sqrt(difference[0] * difference[0] + difference[1] * difference[1] + difference[2] * difference[2]); + } + + static inline void update_velocities(real_t* PHYSICORE_RESTRICT lhs, real_t* PHYSICORE_RESTRICT rhs, + const real_t* PHYSICORE_RESTRICT difference, const real_t force) + { + lhs[0] += force * difference[0]; + lhs[1] += force * difference[1]; + lhs[2] += force * difference[2]; + + rhs[0] -= force * difference[0]; + rhs[1] -= force * difference[1]; + rhs[2] -= force * difference[2]; + } + + static inline void update_velocity(real_t* PHYSICORE_RESTRICT velocity, const real_t* PHYSICORE_RESTRICT difference, + const real_t force) + { + velocity[0] += force * difference[0]; + velocity[1] += force * difference[1]; + velocity[2] += force * difference[2]; + } + + static void random_walk(bool restrict_to_2d, real_t* PHYSICORE_RESTRICT walk) + { + if (restrict_to_2d) + { + position_helper<2>::random_walk(true, walk); + walk[2] = 0; + } + else + { + const real_t theta = random::uniform(0, 2 * std::numbers::pi_v); + const real_t z = random::uniform(-1, 1); + const real_t r = std::sqrt(1 - z * z); + + walk[0] = std::cos(theta) * r; + walk[1] = std::sin(theta) * r; + walk[2] = z; + } + } + + static inline void update_motility_vector(real_t* PHYSICORE_RESTRICT motility_vector, + const real_t* PHYSICORE_RESTRICT walk, + const real_t* PHYSICORE_RESTRICT migration_bias_direction, + const real_t migration_bias) + { + motility_vector[0] = (1 - migration_bias) * walk[0] + migration_bias * migration_bias_direction[0]; + motility_vector[1] = (1 - migration_bias) * walk[1] + migration_bias * migration_bias_direction[1]; + motility_vector[2] = (1 - migration_bias) * walk[2] + migration_bias * migration_bias_direction[2]; + } + + static inline void normalize_and_scale(real_t* PHYSICORE_RESTRICT vector, real_t scale) + { + real_t length = std::sqrt(vector[0] * vector[0] + vector[1] * vector[1] + vector[2] * vector[2]); + + vector[0] = length > zero_threshold ? vector[0] * scale / length : 0; + vector[1] = length > zero_threshold ? vector[1] * scale / length : 0; + vector[2] = length > zero_threshold ? vector[2] * scale / length : 0; + } + + static inline void update_membrane_velocities(real_t* PHYSICORE_RESTRICT velocity, + const real_t* PHYSICORE_RESTRICT position, const cartesian_mesh& mesh, + const real_t radius, const real_t repulsion_strength) + { + update_membrane_velocity(position[0], mesh.bounding_box_mins[0], 1, radius, repulsion_strength, velocity[0]); + update_membrane_velocity(position[0], mesh.bounding_box_maxs[0], -1, radius, repulsion_strength, velocity[0]); + update_membrane_velocity(position[1], mesh.bounding_box_mins[1], 1, radius, repulsion_strength, velocity[1]); + update_membrane_velocity(position[1], mesh.bounding_box_maxs[1], -1, radius, repulsion_strength, velocity[1]); + update_membrane_velocity(position[2], mesh.bounding_box_mins[2], 1, radius, repulsion_strength, velocity[2]); + update_membrane_velocity(position[2], mesh.bounding_box_maxs[2], -1, radius, repulsion_strength, velocity[2]); + } + + static inline void add(real_t* PHYSICORE_RESTRICT lhs, const real_t* PHYSICORE_RESTRICT rhs) + { + lhs[0] += rhs[0]; + lhs[1] += rhs[1]; + lhs[2] += rhs[2]; + } + + static inline void subtract(real_t* PHYSICORE_RESTRICT dst, const real_t* PHYSICORE_RESTRICT lhs, + const real_t* PHYSICORE_RESTRICT rhs) + { + dst[0] = lhs[0] - rhs[0]; + dst[1] = lhs[1] - rhs[1]; + dst[2] = lhs[2] - rhs[2]; + } +}; + + +}; // namespace physicore::mechanics::physicell::kernels::openmp_solver diff --git a/mechanics/physicell/kernels/openmp_solver/tests/CMakeLists.txt b/mechanics/physicell/kernels/openmp_solver/tests/CMakeLists.txt new file mode 100644 index 00000000..64a42c3c --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/tests/CMakeLists.txt @@ -0,0 +1,12 @@ +find_package(GTest REQUIRED) + +file(GLOB TEST_SOURCES *.cpp) +add_executable(mechanics.physicell.kernels.openmp_solver.tests ${TEST_SOURCES}) + +target_link_libraries( + mechanics.physicell.kernels.openmp_solver.tests + PRIVATE mechanics.physicell.kernels.openmp_solver_internal_iface + mechanics.physicell GTest::gtest GTest::gtest_main) + +include(GoogleTest) +gtest_discover_tests(mechanics.physicell.kernels.openmp_solver.tests) diff --git a/mechanics/physicell/kernels/openmp_solver/tests/test_update_cell_forces.cpp b/mechanics/physicell/kernels/openmp_solver/tests/test_update_cell_forces.cpp new file mode 100644 index 00000000..96ab9da7 --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/tests/test_update_cell_forces.cpp @@ -0,0 +1,635 @@ +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "physicell/environment.h" + +using namespace physicore::mechanics::physicell; +using physicore::cartesian_mesh; +using physicore::index_t; +using physicore::real_t; +using physicore::sindex_t; + +namespace { + +class agent_retriever : public physicore::generic_agent_solver +{}; + +mechanical_agent_data& retrieve_environment_agent_data(environment& env) +{ + if (env.agents == nullptr) + { + throw std::runtime_error("environment has no agents"); + } + + return agent_retriever().retrieve_agent_data(*env.agents); +} + +constexpr real_t kTolerance = static_cast(1e-6); + +kernels::openmp_solver::position_solver& position_solver_instance() +{ + static kernels::openmp_solver::position_solver solver; + return solver; +} + +// Helpers for test orchestration +void add_agent(environment& env, std::initializer_list pos, real_t radius = 1, index_t type = 0, + std::uint8_t movable = 1) +{ + ASSERT_TRUE(env.agents != nullptr); + env.agents->create(); + + auto& data = retrieve_environment_agent_data(env); + const index_t dims = data.base_data.dims; + const index_t idx = data.agents_count - 1; + + std::vector p(pos); + p.resize(static_cast(dims), 0); + for (index_t d = 0; d < dims; ++d) + data.base_data.positions[idx * dims + d] = p[d]; + + data.radius[idx] = radius; + data.state_data.agent_type_index[idx] = type; + data.state_data.is_movable[idx] = movable; + + data.mechanics_data.cell_cell_repulsion_strength[idx] = 1; + data.mechanics_data.cell_cell_adhesion_strength[idx] = 1; + data.mechanics_data.relative_maximum_adhesion_distance[idx] = 1; + + for (index_t t = 0; t < std::max(data.agent_types_count, 1); ++t) + data.mechanics_data.cell_adhesion_affinities[idx * data.agent_types_count + t] = 1; +} + +void clear_kinematics_and_pressure(environment& env) +{ + auto& data = retrieve_environment_agent_data(env); + std::ranges::fill(data.velocity, static_cast(0)); + std::ranges::fill(data.previous_velocity, static_cast(0)); + std::ranges::fill(data.state_data.simple_pressure, static_cast(0)); +} + +void connect_pair(environment& env) +{ + auto& data = retrieve_environment_agent_data(env); + data.state_data.neighbors[0] = { 1 }; + data.state_data.neighbors[1] = { 0 }; +} + +void run_update_cell_forces(environment& env) +{ + clear_kinematics_and_pressure(env); + position_solver_instance().update_cell_forces(env); +} + +real_t velocity_component(environment& env, index_t agent, index_t dim) +{ + auto& data = retrieve_environment_agent_data(env); + return data.velocity[agent * data.base_data.dims + dim]; +} + +void set_uniform_affinity(environment& env, real_t value) +{ + auto& data = retrieve_environment_agent_data(env); + for (index_t agent = 0; agent < data.agents_count; ++agent) + for (index_t type = 0; type < data.agent_types_count; ++type) + data.mechanics_data.cell_adhesion_affinities[agent * data.agent_types_count + type] = value; +} + +void set_uniform_position(environment& env, index_t agent, real_t value) +{ + auto& data = retrieve_environment_agent_data(env); + const index_t dims = data.base_data.dims; + for (index_t dim = 0; dim < dims; ++dim) + data.base_data.positions[agent * dims + dim] = value; +} + +cartesian_mesh make_mesh(index_t dims) +{ + if (dims == 1) + return cartesian_mesh(dims, std::array { 0, 0, 0 }, std::array { 20, 0, 0 }, + std::array { 10, 1, 1 }); + if (dims == 2) + return cartesian_mesh(dims, std::array { 0, 0, 0 }, std::array { 20, 20, 0 }, + std::array { 10, 10, 1 }); + return cartesian_mesh(dims, std::array { 0, 0, 0 }, std::array { 20, 20, 20 }, + std::array { 10, 10, 10 }); +} + +std::vector compute_expected_velocities(environment& env) +{ + auto& data = retrieve_environment_agent_data(env); + const index_t dims = data.base_data.dims; + std::vector expected_velocities(static_cast(data.agents_count * dims), 0); + + for (index_t i = 0; i < data.agents_count; ++i) + { + if (data.state_data.is_movable[i] == 0) + continue; + + for (const index_t j : data.state_data.neighbors[i]) + { + std::vector diff(static_cast(dims), 0); + real_t distance_squared = 0; + + for (index_t dim = 0; dim < dims; ++dim) + { + diff[dim] = data.base_data.positions[i * dims + dim] - data.base_data.positions[j * dims + dim]; + distance_squared += diff[dim] * diff[dim]; + } + + const real_t distance = std::max(std::sqrt(distance_squared), static_cast(0.00001)); + + real_t repulsion {}; + { + const real_t repulsive_distance = data.radius[i] + data.radius[j]; + + repulsion = 1 - distance / repulsive_distance; + repulsion = repulsion < 0 ? 0 : repulsion; + repulsion *= repulsion; + repulsion *= std::sqrt(data.mechanics_data.cell_cell_repulsion_strength[i] + * data.mechanics_data.cell_cell_repulsion_strength[j]); + } + + real_t adhesion {}; + { + const real_t adhesion_distance = + data.mechanics_data.relative_maximum_adhesion_distance[i] * data.radius[i] + + data.mechanics_data.relative_maximum_adhesion_distance[j] * data.radius[j]; + + adhesion = 1 - distance / adhesion_distance; + adhesion = adhesion < 0 ? 0 : adhesion; + adhesion *= adhesion; + + const index_t lhs_type = data.state_data.agent_type_index[i]; + const index_t rhs_type = data.state_data.agent_type_index[j]; + + adhesion *= + std::sqrt(data.mechanics_data.cell_cell_adhesion_strength[i] + * data.mechanics_data.cell_cell_adhesion_strength[j] + * data.mechanics_data.cell_adhesion_affinities[i * data.agent_types_count + rhs_type] + * data.mechanics_data.cell_adhesion_affinities[j * data.agent_types_count + lhs_type]); + } + + const real_t force = (repulsion - adhesion) / distance; + + for (index_t dim = 0; dim < dims; ++dim) + expected_velocities[i * dims + dim] += force * diff[dim]; + } + } + + return expected_velocities; +} + +} // namespace + +TEST(UpdateCellForcesTest, NoAgentsDoesNotCrash) +{ + environment env(0.1, 2, 1, 1); + clear_kinematics_and_pressure(env); + position_solver_instance().update_cell_forces(env); + + auto& data = retrieve_environment_agent_data(env); + EXPECT_EQ(data.agents_count, 0); + EXPECT_TRUE(data.velocity.empty()); + EXPECT_TRUE(data.state_data.simple_pressure.empty()); +} + +TEST(UpdateCellForcesTest, SingleAgentNoNeighborsLeavesZeroForces) +{ + environment env(0.1, 2, 1, 1); + add_agent(env, { 0, 0 }); + clear_kinematics_and_pressure(env); + + position_solver_instance().update_cell_forces(env); + + auto& data = retrieve_environment_agent_data(env); + EXPECT_EQ(data.agents_count, 1); + EXPECT_FLOAT_EQ(data.velocity[0], 0); + EXPECT_FLOAT_EQ(data.velocity[1], 0); + EXPECT_FLOAT_EQ(data.state_data.simple_pressure[0], 0); +} + +TEST(UpdateCellForcesTest, TwoAgentsRepelSymmetrically) +{ + environment env(0.1, 2, 1, 1); + add_agent(env, { 0, 0 }); + add_agent(env, { 0.5, 0 }); + + auto& data = retrieve_environment_agent_data(env); + data.radius[0] = 1; + data.radius[1] = 1; + data.mechanics_data.cell_cell_adhesion_strength[0] = 0; + data.mechanics_data.cell_cell_adhesion_strength[1] = 0; + data.state_data.neighbors[0] = { 1 }; + data.state_data.neighbors[1] = { 0 }; + + clear_kinematics_and_pressure(env); + position_solver_instance().update_cell_forces(env); + + EXPECT_LT(data.velocity[0], 0); + EXPECT_GT(data.velocity[2], 0); + EXPECT_NEAR(data.velocity[0] + data.velocity[2], 0, kTolerance); + EXPECT_NEAR(data.velocity[1], 0, kTolerance); + EXPECT_NEAR(data.velocity[3], 0, kTolerance); + EXPECT_GT(data.state_data.simple_pressure[0], 0); + EXPECT_GT(data.state_data.simple_pressure[1], 0); +} + +TEST(UpdateCellForcesTest, OverlappingAgentsProduceFiniteVelocities) +{ + environment env(0.1, 2, 1, 1); + add_agent(env, { 0, 0 }); + add_agent(env, { 0, 0 }); + + auto& data = retrieve_environment_agent_data(env); + data.radius[0] = 1; + data.radius[1] = 1; + data.mechanics_data.cell_cell_adhesion_strength[0] = 0; + data.mechanics_data.cell_cell_adhesion_strength[1] = 0; + data.state_data.neighbors[0] = { 1 }; + + clear_kinematics_and_pressure(env); + position_solver_instance().update_cell_forces(env); + + for (real_t const v : data.velocity) + EXPECT_TRUE(std::isfinite(v)); + for (real_t const p : data.state_data.simple_pressure) + EXPECT_TRUE(std::isfinite(p)); +} + +TEST(UpdateCellForcesTest, AdhesionDependsOnAffinities) +{ + environment env(0.1, 2, 2, 1); + add_agent(env, { 0, 0 }, 1, 0); + add_agent(env, { 1, 0 }, 1, 1); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.cell_cell_repulsion_strength[0] = 0; + data.mechanics_data.cell_cell_repulsion_strength[1] = 0; + data.mechanics_data.relative_maximum_adhesion_distance[0] = 2; + data.mechanics_data.relative_maximum_adhesion_distance[1] = 2; + data.state_data.neighbors[0] = { 1 }; + data.state_data.neighbors[1] = { 0 }; + + data.mechanics_data.cell_adhesion_affinities[0 * 2 + 1] = 1; + data.mechanics_data.cell_adhesion_affinities[1 * 2 + 0] = 1; + clear_kinematics_and_pressure(env); + position_solver_instance().update_cell_forces(env); + EXPECT_GT(data.velocity[0], 0); + EXPECT_LT(data.velocity[2], 0); + + data.mechanics_data.cell_adhesion_affinities[0 * 2 + 1] = 0; + clear_kinematics_and_pressure(env); + position_solver_instance().update_cell_forces(env); + EXPECT_NEAR(data.velocity[0], 0, kTolerance); + EXPECT_NEAR(data.velocity[2], 0, kTolerance); +} + +class SolvePairComplexTest : public ::testing::TestWithParam +{}; + +TEST(SolvePair, RepulsiveForce1D_Overlapping) +{ + environment env(0.1, 1, 1, 1); + add_agent(env, { 0 }); + add_agent(env, { 0.5 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.cell_cell_adhesion_strength[0] = 0; + data.mechanics_data.cell_cell_adhesion_strength[1] = 0; + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_LT(velocity_component(env, 0, 0), 0); + EXPECT_GT(velocity_component(env, 1, 0), 0); + EXPECT_NEAR(velocity_component(env, 0, 0) + velocity_component(env, 1, 0), 0, kTolerance); + EXPECT_GT(data.state_data.simple_pressure[0], 0); + EXPECT_GT(data.state_data.simple_pressure[1], 0); +} + +TEST(SolvePair, NoForce1D_FarApart) +{ + environment env(0.1, 1, 1, 1); + add_agent(env, { 0 }); + add_agent(env, { 5 }); + connect_pair(env); + + run_update_cell_forces(env); + + auto& data = retrieve_environment_agent_data(env); + EXPECT_NEAR(velocity_component(env, 0, 0), 0, kTolerance); + EXPECT_NEAR(velocity_component(env, 1, 0), 0, kTolerance); + EXPECT_NEAR(data.state_data.simple_pressure[0], 0, kTolerance); + EXPECT_NEAR(data.state_data.simple_pressure[1], 0, kTolerance); +} + +TEST(SolvePair, AdhesiveForce1D_InAdhesionRange) +{ + environment env(0.1, 1, 1, 1); + add_agent(env, { 0 }); + add_agent(env, { 3 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.relative_maximum_adhesion_distance[0] = 2; + data.mechanics_data.relative_maximum_adhesion_distance[1] = 2; + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_GT(velocity_component(env, 0, 0), 0); + EXPECT_LT(velocity_component(env, 1, 0), 0); + EXPECT_NEAR(velocity_component(env, 0, 0) + velocity_component(env, 1, 0), 0, kTolerance); + EXPECT_NEAR(data.state_data.simple_pressure[0], 0, kTolerance); + EXPECT_NEAR(data.state_data.simple_pressure[1], 0, kTolerance); +} + +TEST(SolvePair, NewtonsThirdLaw1D_ForceSymmetry) +{ + environment env(0.1, 1, 1, 1); + add_agent(env, { 0 }); + add_agent(env, { 0.5 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.cell_cell_adhesion_strength[0] = 0; + data.mechanics_data.cell_cell_adhesion_strength[1] = 0; + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_NEAR(velocity_component(env, 0, 0), -velocity_component(env, 1, 0), kTolerance); +} + +TEST(SolvePair, SimplePressure1D_Accumulates) +{ + environment env(0.1, 1, 1, 1); + add_agent(env, { 0 }); + add_agent(env, { 0.5 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.cell_cell_adhesion_strength[0] = 0; + data.mechanics_data.cell_cell_adhesion_strength[1] = 0; + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_GT(data.state_data.simple_pressure[0], 0); + EXPECT_GT(data.state_data.simple_pressure[1], 0); + EXPECT_NEAR(data.state_data.simple_pressure[0], data.state_data.simple_pressure[1], kTolerance); +} + +TEST(SolvePair, ZeroRepulsion1D_NoRepulsiveForce) +{ + environment env(0.1, 1, 1, 1); + add_agent(env, { 0 }); + add_agent(env, { 0.5 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.cell_cell_repulsion_strength[0] = 0; + data.mechanics_data.cell_cell_repulsion_strength[1] = 0; + data.mechanics_data.cell_cell_adhesion_strength[0] = 0; + data.mechanics_data.cell_cell_adhesion_strength[1] = 0; + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_NEAR(velocity_component(env, 0, 0), 0, kTolerance); + EXPECT_NEAR(velocity_component(env, 1, 0), 0, kTolerance); + EXPECT_GT(data.state_data.simple_pressure[0], 0); + EXPECT_GT(data.state_data.simple_pressure[1], 0); +} + +TEST(SolvePair, RepulsiveForce2D_Overlapping) +{ + environment env(0.1, 2, 1, 1); + add_agent(env, { 0, 0 }); + add_agent(env, { 0.5, 0.5 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.cell_cell_adhesion_strength[0] = 0; + data.mechanics_data.cell_cell_adhesion_strength[1] = 0; + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_LT(velocity_component(env, 0, 0), 0); + EXPECT_LT(velocity_component(env, 0, 1), 0); + EXPECT_GT(velocity_component(env, 1, 0), 0); + EXPECT_GT(velocity_component(env, 1, 1), 0); + EXPECT_NEAR(velocity_component(env, 0, 0) + velocity_component(env, 1, 0), 0, kTolerance); + EXPECT_NEAR(velocity_component(env, 0, 1) + velocity_component(env, 1, 1), 0, kTolerance); +} + +TEST(SolvePair, AdhesiveForce2D_InAdhesionRange) +{ + environment env(0.1, 2, 1, 1); + add_agent(env, { 0, 0 }); + add_agent(env, { 3, 4 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.relative_maximum_adhesion_distance[0] = 3; + data.mechanics_data.relative_maximum_adhesion_distance[1] = 3; + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_GT(velocity_component(env, 0, 0), 0); + EXPECT_GT(velocity_component(env, 0, 1), 0); + EXPECT_LT(velocity_component(env, 1, 0), 0); + EXPECT_LT(velocity_component(env, 1, 1), 0); +} + +TEST(SolvePair, NewtonsThirdLaw2D_ForceSymmetry) +{ + environment env(0.1, 2, 1, 1); + add_agent(env, { 0, 0 }); + add_agent(env, { 0.5, 0.5 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.cell_cell_adhesion_strength[0] = 0; + data.mechanics_data.cell_cell_adhesion_strength[1] = 0; + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_NEAR(velocity_component(env, 0, 0), -velocity_component(env, 1, 0), kTolerance); + EXPECT_NEAR(velocity_component(env, 0, 1), -velocity_component(env, 1, 1), kTolerance); +} + +TEST(SolvePair, RepulsiveForce3D_Overlapping) +{ + environment env(0.1, 3, 1, 1); + add_agent(env, { 0, 0, 0 }); + add_agent(env, { 0.5, 0.5, 0.5 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.cell_cell_adhesion_strength[0] = 0; + data.mechanics_data.cell_cell_adhesion_strength[1] = 0; + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_LT(velocity_component(env, 0, 0), 0); + EXPECT_LT(velocity_component(env, 0, 1), 0); + EXPECT_LT(velocity_component(env, 0, 2), 0); + EXPECT_GT(velocity_component(env, 1, 0), 0); + EXPECT_GT(velocity_component(env, 1, 1), 0); + EXPECT_GT(velocity_component(env, 1, 2), 0); +} + +TEST(SolvePair, AdhesiveForce3D_DirectionalAccuracy) +{ + environment env(0.1, 3, 1, 1); + add_agent(env, { 0, 0, 0 }); + add_agent(env, { 1, 2, 2 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.relative_maximum_adhesion_distance[0] = 2; + data.mechanics_data.relative_maximum_adhesion_distance[1] = 2; + connect_pair(env); + + run_update_cell_forces(env); + + const real_t vx = velocity_component(env, 0, 0); + const real_t vy = velocity_component(env, 0, 1); + const real_t vz = velocity_component(env, 0, 2); + EXPECT_GT(vx, 0); + EXPECT_GT(vy, 0); + EXPECT_GT(vz, 0); + EXPECT_NEAR(vy / vx, static_cast(2), kTolerance); + EXPECT_NEAR(vz / vx, static_cast(2), kTolerance); +} + +TEST(SolvePair, NewtonsThirdLaw3D_ForceSymmetry) +{ + environment env(0.1, 3, 1, 1); + add_agent(env, { 0, 0, 0 }); + add_agent(env, { 0.5, 0.5, 0.5 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.cell_cell_adhesion_strength[0] = 0; + data.mechanics_data.cell_cell_adhesion_strength[1] = 0; + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_NEAR(velocity_component(env, 0, 0), -velocity_component(env, 1, 0), kTolerance); + EXPECT_NEAR(velocity_component(env, 0, 1), -velocity_component(env, 1, 1), kTolerance); + EXPECT_NEAR(velocity_component(env, 0, 2), -velocity_component(env, 1, 2), kTolerance); +} + +TEST(SolvePair, ZeroDistance1D_Minimum) +{ + environment env(0.1, 1, 1, 1); + add_agent(env, { 0 }); + add_agent(env, { 0 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.cell_cell_adhesion_strength[0] = 0; + data.mechanics_data.cell_cell_adhesion_strength[1] = 0; + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_TRUE(std::isfinite(velocity_component(env, 0, 0))); + EXPECT_TRUE(std::isfinite(velocity_component(env, 1, 0))); + EXPECT_TRUE(std::isfinite(data.state_data.simple_pressure[0])); + EXPECT_TRUE(std::isfinite(data.state_data.simple_pressure[1])); +} + + +TEST(SolvePair, ZeroAffinity1D_NoAdhesion) +{ + environment env(0.1, 1, 1, 1); + add_agent(env, { 0 }); + add_agent(env, { 3 }); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.relative_maximum_adhesion_distance[0] = 2; + data.mechanics_data.relative_maximum_adhesion_distance[1] = 2; + set_uniform_affinity(env, 0); + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_NEAR(velocity_component(env, 0, 0), 0, kTolerance); + EXPECT_NEAR(velocity_component(env, 1, 0), 0, kTolerance); + EXPECT_NEAR(data.state_data.simple_pressure[0], 0, kTolerance); + EXPECT_NEAR(data.state_data.simple_pressure[1], 0, kTolerance); +} + +TEST(SolvePair, DifferentCellTypes1D_AffinityLookup) +{ + environment env(0.1, 1, 2, 1); + add_agent(env, { 0 }, 1, 0); + add_agent(env, { 3 }, 1, 1); + + auto& data = retrieve_environment_agent_data(env); + data.mechanics_data.relative_maximum_adhesion_distance[0] = 2; + data.mechanics_data.relative_maximum_adhesion_distance[1] = 2; + set_uniform_affinity(env, 0); + data.mechanics_data.cell_adhesion_affinities[0 * data.agent_types_count + 1] = 1; + data.mechanics_data.cell_adhesion_affinities[1 * data.agent_types_count + 0] = 1; + connect_pair(env); + + run_update_cell_forces(env); + + EXPECT_GT(velocity_component(env, 0, 0), 0); + EXPECT_LT(velocity_component(env, 1, 0), 0); + EXPECT_NEAR(velocity_component(env, 0, 0) + velocity_component(env, 1, 0), 0, kTolerance); +} + +TEST_P(SolvePairComplexTest, Complex) +{ + const index_t dims = GetParam(); + const real_t radius = dims == 1 ? static_cast(2) : static_cast(4); + + environment env(0.1, dims, 2, 1); + auto mesh = make_mesh(dims); + + add_agent(env, {}, radius, 0); + add_agent(env, {}, radius, 1); + add_agent(env, {}, radius, 0); + + auto& data = retrieve_environment_agent_data(env); + for (index_t agent = 0; agent < data.agents_count; ++agent) + data.mechanics_data.relative_maximum_adhesion_distance[agent] = 1; + + set_uniform_position(env, 0, 1); + set_uniform_position(env, 1, 4); + set_uniform_position(env, 2, 7); + + clear_kinematics_and_pressure(env); + +#pragma omp parallel + { + position_solver_instance().update_cell_neighbors(env, mesh); + position_solver_instance().update_cell_forces(env); + } + + EXPECT_EQ(data.state_data.neighbors[0].size(), static_cast(1)); + EXPECT_EQ(data.state_data.neighbors[1].size(), static_cast(2)); + EXPECT_EQ(data.state_data.neighbors[2].size(), static_cast(1)); + + const auto expected_velocities = compute_expected_velocities(env); + + for (index_t dim = 0; dim < dims; ++dim) + EXPECT_NEAR(data.velocity[dim], expected_velocities[dim], kTolerance); + + for (index_t dim = 0; dim < dims; ++dim) + EXPECT_NEAR(data.velocity[dims + dim], expected_velocities[dims + dim], kTolerance); + + for (index_t dim = 0; dim < dims; ++dim) + EXPECT_NEAR(data.velocity[2 * dims + dim], expected_velocities[2 * dims + dim], kTolerance); +} + +INSTANTIATE_TEST_SUITE_P(AllDims, SolvePairComplexTest, ::testing::Values(1, 2, 3)); diff --git a/mechanics/physicell/kernels/openmp_solver/tests/test_update_cell_neighbors.cpp b/mechanics/physicell/kernels/openmp_solver/tests/test_update_cell_neighbors.cpp new file mode 100644 index 00000000..513d9018 --- /dev/null +++ b/mechanics/physicell/kernels/openmp_solver/tests/test_update_cell_neighbors.cpp @@ -0,0 +1,148 @@ +#include +#include +#include + +#include +#include +#include +#include + +#include "physicell/environment.h" + +using namespace physicore::mechanics::physicell; +using physicore::cartesian_mesh; +using physicore::index_t; +using physicore::real_t; +using physicore::sindex_t; + +namespace { + +class agent_retriever : public physicore::generic_agent_solver +{}; + +mechanical_agent_data& retrieve_environment_agent_data(environment& env) +{ + if (env.agents == nullptr) + { + throw std::runtime_error("environment has no agents"); + } + + return agent_retriever().retrieve_agent_data(*env.agents); +} + +kernels::openmp_solver::position_solver& position_solver_instance() +{ + static kernels::openmp_solver::position_solver solver; + return solver; +} + +void add_agent(environment& env, std::initializer_list pos, real_t radius = 1, std::uint8_t movable = 1, + real_t rel_max_adhesion_dist = 1) +{ + ASSERT_TRUE(env.agents != nullptr); + env.agents->create(); + + auto& data = retrieve_environment_agent_data(env); + const index_t dims = data.base_data.dims; + const index_t idx = data.agents_count - 1; + + std::vector p(pos); + p.resize(static_cast(dims), 0); + for (index_t d = 0; d < dims; ++d) + data.base_data.positions[idx * dims + d] = p[d]; + + data.radius[idx] = radius; + data.state_data.is_movable[idx] = movable; + data.mechanics_data.relative_maximum_adhesion_distance[idx] = rel_max_adhesion_dist; + + data.mechanics_data.cell_cell_repulsion_strength[idx] = 1; + data.mechanics_data.cell_cell_adhesion_strength[idx] = 1; + + for (index_t t = 0; t < std::max(data.agent_types_count, 1); ++t) + data.mechanics_data.cell_adhesion_affinities[idx * data.agent_types_count + t] = 1; +} + +std::vector sorted_neighbors(environment& env, index_t i) +{ + auto neighbors = retrieve_environment_agent_data(env).state_data.neighbors[i]; + std::ranges::sort(neighbors); + return neighbors; +} + +cartesian_mesh make_mesh(index_t dims) +{ + return cartesian_mesh(dims, std::array { 0, 0, 0 }, std::array { 10, 10, 0 }, + std::array { 2, 2, 1 }); +} + +} // namespace + +TEST(UpdateCellNeighborsTest, NoAgentsDoesNotCrash) +{ + environment env(0.1, 2, 1, 1); + auto mesh = make_mesh(2); + position_solver_instance().update_cell_neighbors(env, mesh); + EXPECT_EQ(retrieve_environment_agent_data(env).agents_count, 0); +} + +TEST(UpdateCellNeighborsTest, SingleAgentHasNoNeighbors) +{ + environment env(0.1, 2, 1, 1); + add_agent(env, { 0, 0 }); + + auto mesh = make_mesh(2); + position_solver_instance().update_cell_neighbors(env, mesh); + + ASSERT_EQ(retrieve_environment_agent_data(env).agents_count, 1); + EXPECT_TRUE(retrieve_environment_agent_data(env).state_data.neighbors[0].empty()); +} + +TEST(UpdateCellNeighborsTest, DistanceEqualThresholdCountsAsNeighbor) +{ + environment env(0.1, 2, 1, 1); + add_agent(env, { 0, 0 }, 1, 1, 1); + add_agent(env, { 2, 0 }, 1, 1, 1); // adhesion_distance = 1*1 + 1*1 = 2 + + auto mesh = make_mesh(2); + position_solver_instance().update_cell_neighbors(env, mesh); + + EXPECT_EQ(sorted_neighbors(env, 0), (std::vector { 1 })); + EXPECT_EQ(sorted_neighbors(env, 1), (std::vector { 0 })); +} + +TEST(UpdateCellNeighborsTest, DistanceAboveThresholdIsNotNeighbor) +{ + environment env(0.1, 2, 1, 1); + add_agent(env, { 0, 0 }, 1, 1, 1); + add_agent(env, { 2.0001, 0 }, 1, 1, 1); + + auto mesh = make_mesh(2); + position_solver_instance().update_cell_neighbors(env, mesh); + + EXPECT_TRUE(retrieve_environment_agent_data(env).state_data.neighbors[0].empty()); + EXPECT_TRUE(retrieve_environment_agent_data(env).state_data.neighbors[1].empty()); +} + +TEST(UpdateCellNeighborsTest, ClearsPreviousNeighborsAndRespectsMovableFlag) +{ + environment env(0.1, 2, 1, 1); + add_agent(env, { 0, 0 }, 1, 1, 1); // movable + add_agent(env, { 1, 0 }, 1, 0, 1); // immovable but within threshold + add_agent(env, { 10, 0 }, 1, 1, 1); + + auto& data = retrieve_environment_agent_data(env); + data.state_data.neighbors[0] = { 2, 12345 }; // garbage to prove clear + data.state_data.neighbors[1] = { 0, 2 }; // should be cleared, then skipped (immovable) + + auto mesh = make_mesh(2); + position_solver_instance().update_cell_neighbors(env, mesh); + + EXPECT_EQ(sorted_neighbors(env, 0), (std::vector { 1 })); + EXPECT_TRUE(retrieve_environment_agent_data(env).state_data.neighbors[1].empty()); + EXPECT_TRUE(retrieve_environment_agent_data(env).state_data.neighbors[2].empty()); + + // Integration: neighbor list should drive a non-zero force. + std::ranges::fill(data.velocity, static_cast(0)); + position_solver_instance().update_cell_forces(env); + EXPECT_NEAR(data.velocity[0], -data.velocity[2], 1e-6); +} diff --git a/mechanics/physicell/src/agent_data.cpp b/mechanics/physicell/src/agent_data.cpp index 8b1796d3..b7070a24 100644 --- a/mechanics/physicell/src/agent_data.cpp +++ b/mechanics/physicell/src/agent_data.cpp @@ -79,6 +79,7 @@ void mechanical_agent_data::remove_at(index_t position) base_storage_t::move_vector(&motility_data.chemotactic_sensitivities[position * substrates_count], &motility_data.chemotactic_sensitivities[last * substrates_count], substrates_count); + motility_data.direction_update_funcs[position] = std::move(motility_data.direction_update_funcs[last]); // State properties state_data.neighbors[position] = std::move(state_data.neighbors[last]); @@ -123,6 +124,7 @@ void mechanical_agent_data::resize_storage() motility_data.migration_bias.resize(agents_count); motility_data.motility_vector.resize(agents_count * base_data.dims); motility_data.restrict_to_2d.resize(agents_count); + motility_data.direction_update_funcs.resize(agents_count); motility_data.chemotaxis_index.resize(agents_count); motility_data.chemotaxis_direction.resize(agents_count); motility_data.chemotactic_sensitivities.resize(agents_count * substrates_count); diff --git a/mechanics/physicell/src/environment.cpp b/mechanics/physicell/src/environment.cpp index d1054902..c8f11502 100644 --- a/mechanics/physicell/src/environment.cpp +++ b/mechanics/physicell/src/environment.cpp @@ -1,14 +1,22 @@ #include "environment.h" +#include +#include + +#include + using namespace physicore::mechanics::physicell; -void environment::run_single_timestep() +environment::environment(real_t timestep, index_t dims, index_t agent_types_count, index_t substrates_count) + : mechanics_timestep(timestep) { - // Implementation of the single timestep logic - // This is a placeholder; actual implementation will depend on the specific requirements - (void)timestep; + auto base_data = std::make_unique(dims); + auto data = std::make_unique(*base_data, agent_types_count, substrates_count); + agents = std::make_unique(std::move(base_data), std::move(data)); } +void environment::run_single_timestep() { solver->solve(*this, 1); } + void environment::serialize_state(real_t current_time) { if (serializer) @@ -16,3 +24,12 @@ void environment::serialize_state(real_t current_time) serializer->serialize(current_time); } } + +const physicore::cartesian_mesh& environment::get_mesh() const +{ + if (!mesh_) + { + throw std::runtime_error("environment has no mesh"); + } + return *mesh_; +} diff --git a/mechanics/physicell/src/register_solver.cpp b/mechanics/physicell/src/register_solver.cpp new file mode 100644 index 00000000..6459171d --- /dev/null +++ b/mechanics/physicell/src/register_solver.cpp @@ -0,0 +1,10 @@ +#include "physicell/openmp_solver/register_solver.h" + +#include + +#include "openmp_solver.h" + +void physicore::mechanics::physicell::kernels::openmp_solver::attach_to_registry() +{ + static const physicore::mechanics::physicell::registry_adder openmp_solver_adder("openmp_solver"); +} diff --git a/mechanics/physicell/src/solver_registry.cpp b/mechanics/physicell/src/solver_registry.cpp new file mode 100644 index 00000000..90497f2b --- /dev/null +++ b/mechanics/physicell/src/solver_registry.cpp @@ -0,0 +1,20 @@ +#include "physicell/solver_registry.h" + +#include "physicell/openmp_solver/register_solver.h" + +namespace physicore::mechanics::physicell { + +solver_registry& solver_registry::instance() +{ + static solver_registry r; + return r; +} + +struct attachment_point +{ + attachment_point() { kernels::openmp_solver::attach_to_registry(); } +}; + +static const attachment_point ap; + +} // namespace physicore::mechanics::physicell diff --git a/mechanics/physicell/tests/test_environment.cpp b/mechanics/physicell/tests/test_environment.cpp index 1c69fe62..7c096cf8 100644 --- a/mechanics/physicell/tests/test_environment.cpp +++ b/mechanics/physicell/tests/test_environment.cpp @@ -1,12 +1,97 @@ +#include +#include + #include -#include "environment.h" +#include "physicell/environment.h" +#include "physicell/serializer.h" +#include "physicell/solver.h" + +namespace physicore::mechanics::physicell::tests { +using namespace physicore; using namespace physicore::mechanics::physicell; -TEST(EnvironmentTest, RunSingleTimestep) +// Mock solver for testing +class MockSolver : public solver +{ +public: + void initialize(environment& /* e */) override {} + void solve(environment& /* e */, index_t /* iterations */) override { solve_called = true; } + bool solve_called = false; +}; + +// Mock serializer for testing +class MockSerializer : public serializer +{ +public: + void serialize(real_t /* current_time */) override { serialize_called = true; } + bool serialize_called = false; +}; + +class EnvironmentTest : public ::testing::Test +{}; + +// Test constructor and basic initialization +TEST_F(EnvironmentTest, ConstructorInitializesEnvironment) { - environment env(0.1); + environment env(0.1, 3, 2, 3); + EXPECT_DOUBLE_EQ(env.mechanics_timestep, 0.1); + EXPECT_FALSE(env.automated_spring_adhesion); + EXPECT_TRUE(env.virtual_wall_at_domain_edges); + ASSERT_NE(env.agents, nullptr); + EXPECT_EQ(env.agents->size(), 0); +} + +// Test run_single_timestep calls solver and serialize_state calls serializer +TEST_F(EnvironmentTest, RunSingleTimestepAndSerializeState) +{ + environment env(0.1, 3, 2, 3); + + auto mock_solver = std::make_unique(); + auto mock_serializer = std::make_unique(); + + auto* solver_ptr = mock_solver.get(); + auto* serializer_ptr = mock_serializer.get(); + + env.solver = std::move(mock_solver); + env.serializer = std::move(mock_serializer); + env.run_single_timestep(); - // No assertion, just check it runs without error + EXPECT_TRUE(solver_ptr->solve_called); + + env.serialize_state(1.5); + EXPECT_TRUE(serializer_ptr->serialize_called); } + +// Test set_mesh and has_mesh methods +TEST_F(EnvironmentTest, SetMeshAndHasMesh) +{ + environment env(0.1, 3, 2, 3); + + // Initially has_mesh returns false + EXPECT_FALSE(env.has_mesh()); + + // Set a mesh and verify has_mesh returns true + auto mesh = physicore::cartesian_mesh(3, { -100, -100, -100 }, { 100, 100, 100 }, { 20, 20, 20 }); + env.set_mesh(mesh); + EXPECT_TRUE(env.has_mesh()); +} + +// Test mesh operations and exception handling +TEST_F(EnvironmentTest, MeshOperationsAndErrorHandling) +{ + environment env(0.1, 3, 2, 3); + + // Test exception when getting mesh without setting it + EXPECT_THROW(env.get_mesh(), std::runtime_error); + + // Set a mesh and verify retrieval + auto mesh = physicore::cartesian_mesh(3, { -100, -100, -100 }, { 100, 100, 100 }, { 20, 20, 20 }); + env.set_mesh(mesh); + + const auto& retrieved = env.get_mesh(); + EXPECT_EQ(retrieved.dims, 3); +} + +} // namespace physicore::mechanics::physicell::tests diff --git a/mechanics/physicell/tests/test_registry.cpp b/mechanics/physicell/tests/test_registry.cpp new file mode 100644 index 00000000..6abc10c8 --- /dev/null +++ b/mechanics/physicell/tests/test_registry.cpp @@ -0,0 +1,33 @@ +#include + +#include "physicell/solver.h" +#include "physicell/solver_registry.h" + +namespace physicore::mechanics::physicell::tests { + +using namespace physicore; +using namespace physicore::mechanics::physicell; + +class RegistryTest : public ::testing::Test +{}; + +// Test registry covers the if body in registry.h::37 (key not found case) +TEST_F(RegistryTest, GetNonExistentSolverReturnsNullptr) +{ + // Try to get a solver that doesn't exist in the registry + // This exercises the if (it == factories_.end()) branch at registry.h::37 + auto result = solver_registry::instance().get("nonexistent_solver_xyz"); + EXPECT_EQ(result, nullptr); +} + +// Test registry covers the else case (key found) +TEST_F(RegistryTest, GetExistentSolverReturnsValidInstance) +{ + // Get the registered openmp solver + // This exercises the else branch where the factory is found and called + auto result = solver_registry::instance().get("openmp_solver"); + ASSERT_NE(result, nullptr); + EXPECT_TRUE(dynamic_cast(result.get()) != nullptr); +} + +} // namespace physicore::mechanics::physicell::tests diff --git a/mechanics/physicell/tests/test_solver_registry.cpp b/mechanics/physicell/tests/test_solver_registry.cpp new file mode 100644 index 00000000..68e8d661 --- /dev/null +++ b/mechanics/physicell/tests/test_solver_registry.cpp @@ -0,0 +1,16 @@ +#include +#include + +#include +#include + +#include "physicell/solver_registry.h" + +using namespace physicore; +using namespace physicore::mechanics::physicell; + +TEST(SolverRegistryTest, OpenMPSolverIsRegistered) +{ + auto solver_instance = solver_registry::instance().get("openmp_solver"); + ASSERT_NE(solver_instance, nullptr); +} diff --git a/mechanics/physicell/tests/test_vtk_agents_serializer.cpp b/mechanics/physicell/tests/test_vtk_agents_serializer.cpp index bd41b503..3dcbe34b 100644 --- a/mechanics/physicell/tests/test_vtk_agents_serializer.cpp +++ b/mechanics/physicell/tests/test_vtk_agents_serializer.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include #include diff --git a/reactions-diffusion/biofvm/CMakeLists.txt b/reactions-diffusion/biofvm/CMakeLists.txt index ae6dbbb5..e2b60050 100644 --- a/reactions-diffusion/biofvm/CMakeLists.txt +++ b/reactions-diffusion/biofvm/CMakeLists.txt @@ -22,7 +22,6 @@ target_sources( reactions-diffusion.biofvm PRIVATE src/microenvironment.cpp src/microenvironment_builder.cpp - src/mesh.cpp src/vtk_serializer.cpp src/vtk_serializer_base.cpp src/config_reader.cpp @@ -47,7 +46,8 @@ target_sources( INTERFACE FILE_SET HEADERS BASE_DIRS "${CMAKE_CURRENT_SOURCE_DIR}/include" "${CMAKE_CURRENT_BINARY_DIR}/include") -target_compile_definitions(reactions-diffusion.biofvm_iface INTERFACE BIOFVM_STATIC_DEFINE) +target_compile_definitions(reactions-diffusion.biofvm_iface + INTERFACE BIOFVM_STATIC_DEFINE) target_link_libraries(reactions-diffusion.biofvm_internal_iface INTERFACE common) diff --git a/reactions-diffusion/biofvm/include/biofvm/mesh.h b/reactions-diffusion/biofvm/include/biofvm/mesh.h deleted file mode 100644 index 43e16e36..00000000 --- a/reactions-diffusion/biofvm/include/biofvm/mesh.h +++ /dev/null @@ -1,34 +0,0 @@ -#pragma once - -#include -#include - -#include -#include - -namespace physicore::biofvm { - -struct BIOFVM_EXPORT cartesian_mesh -{ - index_t dims; // 1 or 2 or 3 - - std::array bounding_box_mins; // [x_min, y_min, z_min] - std::array bounding_box_maxs; // [x_max, y_max, z_max] - - std::array voxel_shape; // [dx, dy, dz] - std::array grid_shape; // [x_size, y_size, z_size] - - cartesian_mesh(index_t dims, std::array bounding_box_mins, std::array bounding_box_maxs, - std::array voxel_shape); - - std::size_t voxel_count() const; - index_t voxel_volume() const; - - std::size_t linearize(index_t x, index_t y, index_t z) const; - - std::array voxel_position(std::span position) const; - - std::array voxel_center(std::array position) const; -}; - -} // namespace physicore::biofvm diff --git a/reactions-diffusion/biofvm/include/biofvm/microenvironment.h b/reactions-diffusion/biofvm/include/biofvm/microenvironment.h index d9002d00..886131bf 100644 --- a/reactions-diffusion/biofvm/include/biofvm/microenvironment.h +++ b/reactions-diffusion/biofvm/include/biofvm/microenvironment.h @@ -6,12 +6,12 @@ #include #include +#include #include #include #include "agent_container.h" #include "bulk_functor.h" -#include "mesh.h" #include "serializer.h" #include "solver.h" diff --git a/reactions-diffusion/biofvm/include/biofvm/microenvironment_builder.h b/reactions-diffusion/biofvm/include/biofvm/microenvironment_builder.h index 9476e5b1..9472a597 100644 --- a/reactions-diffusion/biofvm/include/biofvm/microenvironment_builder.h +++ b/reactions-diffusion/biofvm/include/biofvm/microenvironment_builder.h @@ -3,10 +3,10 @@ #include #include +#include #include #include "bulk_functor.h" -#include "mesh.h" #include "microenvironment.h" namespace physicore::biofvm { diff --git a/reactions-diffusion/biofvm/include/biofvm/solver_registry.h b/reactions-diffusion/biofvm/include/biofvm/solver_registry.h index 69bb5433..f70974a4 100644 --- a/reactions-diffusion/biofvm/include/biofvm/solver_registry.h +++ b/reactions-diffusion/biofvm/include/biofvm/solver_registry.h @@ -1,27 +1,14 @@ #pragma once -#include -#include -#include -#include -#include +#include #include "solver.h" namespace physicore::biofvm { -class solver_registry +class solver_registry : public physicore::factory_registry { public: - using solver_factory_func_t = std::function; - using registry_map_t = std::unordered_map; - - registry_map_t factory_registry; - - bool register_factory(std::string solver_name, solver_factory_func_t&& f); - - solver_ptr get(const std::string& solver_name); - static solver_registry& instance(); }; @@ -29,13 +16,5 @@ template concept derived_from_solver = std::derived_from; template -struct registry_adder -{ - explicit registry_adder(std::string solver_name) - { - solver_registry::instance().register_factory(std::move(solver_name), - []() { return std::make_unique(); }); - } -}; - +using registry_adder = physicore::generic_registry_adder; } // namespace physicore::biofvm diff --git a/reactions-diffusion/biofvm/src/solver_registry.cpp b/reactions-diffusion/biofvm/src/solver_registry.cpp index 99a6be20..b8f27a58 100644 --- a/reactions-diffusion/biofvm/src/solver_registry.cpp +++ b/reactions-diffusion/biofvm/src/solver_registry.cpp @@ -1,7 +1,5 @@ #include "solver_registry.h" -#include - #include "openmp_solver/register_solver.h" #if defined(PHYSICORE_HAS_TBB_THRUST) or defined(PHYSICORE_HAS_CUDA_THRUST) @@ -10,31 +8,12 @@ using namespace physicore::biofvm; -bool solver_registry::register_factory(std::string solver_name, solver_factory_func_t&& f) -{ - auto [it, emplaced] = factory_registry.try_emplace(std::move(solver_name), std::move(f)); - - return emplaced; -} - -std::unique_ptr solver_registry::get(const std::string& solver_name) -{ - if (!factory_registry.contains(solver_name)) - { - assert(false); - return nullptr; - } - - return factory_registry[solver_name](); -} - solver_registry& solver_registry::instance() { static solver_registry r; return r; } - struct attachment_point { attachment_point()