OpenMP solvers

.devcontainer/physicore-dev/devcontainer.json

@@ -8,10 +8,11 @@
 	],
 	"remoteEnv": {
 		"VCPKG_ROOT": "${containerWorkspaceFolder}/vcpkg",
-		"PATH": "${containerEnv:PATH}:${containerWorkspaceFolder}/vcpkg:/usr/lib/llvm-20/bin"
+		"PATH": "${containerEnv:PATH}:${containerWorkspaceFolder}/vcpkg:/usr/lib/llvm-20/bin",
+		"APT_KEY_DONT_WARN_ON_DANGEROUS_USAGE": "1"
 	},
 	"initializeCommand": "mkdir -p ${localEnv:HOME}/.cache/vcpkg",
-	"onCreateCommand": ".devcontainer/configure.sh",
+	"postCreateCommand": ".devcontainer/configure.sh",
 	"features": {
 		"ghcr.io/devcontainers/features/python:1": {
 			"toolsToInstall": "cmakelang"

.devcontainer/physicore-dev/dockerfile

@@ -1,9 +1,14 @@
 FROM gcc:15.2.0
 
-# Install LLVM 20
-RUN apt-get update && \
+# Configure APT to allow unsigned repositories due to LLVM key expiration
+# This is needed because the LLVM 20 signing key uses SHA1 which is no longer secure as of 2026-02-01
+RUN echo "APT::Get::AllowUnauthenticated \"true\";" > /etc/apt/apt.conf.d/99-allow-unsigned && \
+    echo "Acquire::AllowInsecureRepositories \"true\";" >> /etc/apt/apt.conf.d/99-allow-unsigned && \
+    echo "Acquire::AllowDowngradeToInsecureRepositories \"true\";" >> /etc/apt/apt.conf.d/99-allow-unsigned
+
+# Remove the problematic LLVM repository entirely to avoid signature issues during feature installations
+# Users can still install LLVM manually if needed
+RUN apt-get update || true && \
     apt-get install -y --no-install-recommends \
         lsb-release wget && \
-    wget https://apt.llvm.org/llvm.sh && \
-    chmod +x ./llvm.sh && ./llvm.sh 20 all && rm llvm.sh && \
     apt-get clean && rm -rf /var/lib/apt/lists/*

.github/workflows/cmake-multi-platform.yml

@@ -39,6 +39,19 @@ jobs:
         sudo apt-get update
         sudo apt-get install -y libomp-dev
 
+    - name: Install libomp (macOS)
+      if: matrix.os == 'macos-latest'
+      run: |
+        brew install libomp
+        echo "OpenMP_ROOT=$(brew --prefix libomp)" >> "$GITHUB_ENV"
+
+    - name: Install Ninja (Windows)
+      if: runner.os == 'Windows'
+      shell: pwsh
+      run: |
+        choco install ninja -y --no-progress
+        ninja --version
+
     - name: Install CUDA Toolkit
       uses: ./.github/actions/cuda-toolkit
 

.github/workflows/release.yml

@@ -50,6 +50,13 @@ jobs:
           submodules: recursive
           fetch-depth: 0
 
+      - name: Install Ninja (Windows)
+        if: runner.os == 'Windows'
+        shell: pwsh
+        run: |
+          choco install ninja -y --no-progress
+          ninja --version
+
       - name: Install CUDA Toolkit
         uses: ./.github/actions/cuda-toolkit
 

CMakePresets.json

@@ -91,6 +91,7 @@
             "cacheVariables": {
                 "CMAKE_C_COMPILER": "cc",
                 "CMAKE_CXX_COMPILER": "c++",
+                "OpenMP_ROOT": "$env{OpenMP_ROOT}",
                 "CMAKE_C_FLAGS_DEBUG_INIT": "-Wall -Wextra -Wpedantic -fsanitize=address -fsanitize=leak -fsanitize=pointer-compare -fsanitize=pointer-subtract -fsanitize=undefined",
                 "CMAKE_CXX_FLAGS_DEBUG_INIT": "-Wall -Wextra -Wpedantic -fsanitize=address -fsanitize=leak -fsanitize=pointer-compare -fsanitize=pointer-subtract -fsanitize=undefined",
                 "CMAKE_C_FLAGS_RELWITHDEBINFO_INIT": "-Wall -Wextra -Wpedantic -fsanitize=thread",

Testing/Temporary/CTestCostData.txt

@@ -0,0 +1 @@
+---

Testing/Temporary/LastTest.log

@@ -0,0 +1,3 @@
+Start testing: Feb 24 13:42 UTC
+----------------------------------------------------------
+End testing: Feb 24 13:42 UTC

common/CMakeLists.txt

@@ -1,9 +1,13 @@
-add_library(common INTERFACE)
+add_library(common)
 
 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)

common/include/common/cartesian_mesh.h

@@ -0,0 +1,62 @@
+#pragma once
+
+#include <array>
+#include <span>
+
+#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<sindex_t, 3> bounding_box_mins,
+							std::array<sindex_t, 3> bounding_box_maxs, std::array<index_t, 3> 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<index_t, 3> voxel_position(std::span<const real_t> 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<real_t, 3> voxel_center(std::array<index_t, 3> 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<sindex_t, 3> bounding_box_mins; ///< Minimum domain coordinates
+	std::array<sindex_t, 3> bounding_box_maxs; ///< Maximum domain coordinates
+	std::array<index_t, 3> voxel_shape;		   ///< Size of each voxel
+	std::array<index_t, 3> grid_shape;		   ///< Number of voxels per dimension
+};
+
+} // namespace physicore

common/include/common/random.h

@@ -0,0 +1,19 @@
+#pragma once
+
+#include <common/types.h>
+
+namespace physicore {
+
+class random
+{
+public:
+	static random& instance();
+
+	real_t uniform(const real_t min = 0, const real_t max = 1);
+
+	real_t normal(const real_t mean = 0, const real_t std = 1);
+
+	void set_seed(unsigned int seed);
+};
+
+} // namespace physicore

common/include/common/registry.h

@@ -0,0 +1,56 @@
+#pragma once
+
+#include <cassert>
+#include <concepts>
+#include <functional>
+#include <memory>
+#include <string>
+#include <unordered_map>
+
+namespace physicore::common {
+template <class Base>
+class factory_registry
+{
+public:
+	using ptr_t = std::unique_ptr<Base>;
+	using factory_func_t = std::function<ptr_t()>;
+	using map_t = std::unordered_map<std::string, factory_func_t>;
+
+	bool register_factory(std::string name, factory_func_t&& f);
+	ptr_t get(const std::string& name);
+
+private:
+	map_t factories_;
+};
+
+template <class Base>
+bool factory_registry<Base>::register_factory(std::string name, factory_func_t&& f)
+{
+	auto [it, emplaced] = factories_.try_emplace(std::move(name), std::move(f));
+	return emplaced;
+}
+
+template <class Base>
+typename factory_registry<Base>::ptr_t factory_registry<Base>::get(const std::string& name)
+{
+	auto it = factories_.find(name);
+	if (it == factories_.end())
+	{
+		assert(false);
+		return nullptr;
+	}
+
+	return it->second();
+}
+
+template <class Derived, class Registry, class Base>
+	requires std::derived_from<Derived, Base>
+struct registry_adder
+{
+	explicit registry_adder(std::string name)
+	{
+		Registry::instance().register_factory(std::move(name), []() { return std::make_unique<Derived>(); });
+	}
+};
+
+} // namespace physicore::common

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

common/src/cartesian_mesh.cpp

@@ -0,0 +1,78 @@
+#include "common/cartesian_mesh.h"
+
+#include <cassert>
+
+namespace physicore {
+
+cartesian_mesh::cartesian_mesh(index_t dims, std::array<sindex_t, 3> bounding_box_mins,
+							   std::array<sindex_t, 3> bounding_box_maxs, std::array<index_t, 3> 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<index_t, 3> cartesian_mesh::voxel_position(std::span<const real_t> 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<real_t, 3> cartesian_mesh::voxel_center(std::array<index_t, 3> 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

common/src/random.cpp

@@ -0,0 +1,52 @@
+#ifdef _OPENMP
+	#include <omp.h>
+#endif
+
+#include <random>
+#include <vector>
+
+#include <common/random.h>
+
+
+thread_local std::mt19937 generator;
+
+physicore::random& physicore::random::instance()
+{
+	static random instance;
+	return instance;
+}
+
+physicore::real_t physicore::random::uniform(const physicore::real_t min, const physicore::real_t max)
+{
+	std::uniform_real_distribution<physicore::real_t> distribution(min, max);
+	return distribution(generator);
+}
+
+physicore::real_t physicore::random::normal(const physicore::real_t mean, const physicore::real_t std)
+{
+	std::normal_distribution<physicore::real_t> distribution(mean, std);
+	return distribution(generator);
+}
+
+void physicore::random::set_seed(unsigned int seed)
+{
+#ifdef _OPENMP
+	std::vector<unsigned int> 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<unsigned int> seeds(omp_get_num_threads());
+	seq.generate(seeds.begin(), seeds.end());
+
+	#pragma omp parallel
+	{
+		int id = omp_get_thread_num();
+		generator.seed(seeds[id]);
+	}
+#else
+	generator.seed(seed);
+#endif
+}