Dev/mech micro #79
Dev/mech micro #79
Conversation
- Introduce simulation_parameters struct to centralize mechanics configuration - Define interaction_config for potential-specific settings - Add support for type-based interaction parameters - Resolves #26
- Add agent_data struct for Structure-of-Arrays (SoA) storage - Implement agent proxy class for object-oriented access - Add unit tests for agent data management and accessors - Resolves #28
- Add agent_container class for managing agent lifecycle - Implement create() and remove() operations - Ensure proper synchronization with base_agent_data - Add unit tests for container operations - Resolves #29
- Add vtk_mechanics_serializer for VTK output - Implement serialization of agent positions, velocities, and forces - Add cell_data structure for aggregated cell-level data - Resolves #37
- Define abstract solver base class - Implement solver_registry for runtime backend selection - Add unit tests for registry mechanics
- Add abstract spatial_index interface - Implement uniform_grid_spatial_index for efficient neighbor queries - Add unit tests for spatial indexing
- Add potential_interface for force models - Implement Standard (PhysiCell), Morse, and Kelvin-Voigt potentials - Add unit tests for all potential types
- Implement OpenMP-parallelized solver class - Add sub-solvers for forces, motility, neighbors, and positions - Integrate potentials and spatial index into solver loop - Add solver registration logic - Resolves #54
- Add Mechanics section to README - Describe environment, solver, and serialization usage - Provide minimal code example - Resolves #45
- Move VTK serializer base class from biofvm to common - Update biofvm serializers to use common base - Enable code reuse for mechanics serializer
- Update environment to use agent_container and solver_registry - Wire up simulation loop in run_single_timestep - Update build configuration to link new components - Update environment tests
There was a problem hiding this comment.
Pull request overview
This PR implements a comprehensive micromechanics simulation module for PhysiCore, addressing issues #26, #28, #29, #37, #45, and #54. The implementation follows PhysiCore's modular architecture with pluggable solvers, providing OpenMP-parallelized cell mechanics including neighbor detection, multiple force potentials (standard, Morse, Kelvin-Voigt), and position integration.
Key Changes
- Refactored
vtk_serializer_basefromphysicore::biofvmtophysicore::commonnamespace for reusability - Added complete micromechanics module with agent-based mechanics simulation
- Implemented solver registry pattern matching BioFVM's pluggable solver design
- Added comprehensive test suite covering all major components
Reviewed changes
Copilot reviewed 61 out of 61 changed files in this pull request and generated 4 comments.
Show a summary per file
| File | Description |
|---|---|
common/include/common/vtk_serializer_base.h |
Moved VTK serializer base class to common namespace for reuse across modules |
common/CMakeLists.txt |
Split common library to separate VTK-dependent code into optional target |
reactions-diffusion/biofvm/src/vtk_serializer.h |
Updated to use common namespace for vtk_serializer_base |
reactions-diffusion/biofvm/src/vtk_agents_serializer.h |
Updated to use common namespace for vtk_serializer_base |
reactions-diffusion/biofvm/CMakeLists.txt |
Added link dependency on vtk_serializer library |
mechanics/micromechanics/include/micromechanics/environment.h |
Main simulation environment with agents, solver, spatial index |
mechanics/micromechanics/include/micromechanics/agent_data.h |
SoA storage for agent mechanical properties |
mechanics/micromechanics/include/micromechanics/cell_data.h |
Aggregated cell-level data computed from agents |
mechanics/micromechanics/include/micromechanics/solver.h |
Abstract solver interface for pluggable mechanics backends |
mechanics/micromechanics/include/micromechanics/solver_registry.h |
Registry for runtime solver selection |
mechanics/micromechanics/include/micromechanics/simulation_parameters.h |
Type-based interaction configuration system |
mechanics/micromechanics/include/micromechanics/potential_interface.h |
Interface for pairwise force potentials |
mechanics/micromechanics/include/micromechanics/spatial_index.h |
Interface for neighbor finding data structures |
mechanics/micromechanics/include/micromechanics/uniform_grid_spatial_index.h |
Grid-based spatial indexing implementation |
mechanics/micromechanics/src/environment.cpp |
Environment implementation with solver orchestration |
mechanics/micromechanics/src/uniform_grid_spatial_index.cpp |
Uniform grid neighbor search implementation |
mechanics/micromechanics/src/vtk_mechanics_serializer.cpp |
VTK output for mechanics simulation state |
mechanics/micromechanics/src/solver_registry_sole.cpp |
Singleton solver registry implementation |
mechanics/micromechanics/kernels/openmp_solver/src/openmp_solver.cpp |
Main OpenMP solver coordinating sub-solvers |
mechanics/micromechanics/kernels/openmp_solver/src/force_solver.cpp |
Pairwise force calculation with type-based potentials |
mechanics/micromechanics/kernels/openmp_solver/src/neighbor_solver.cpp |
Neighbor list building via spatial index |
mechanics/micromechanics/kernels/openmp_solver/src/position_solver.cpp |
Adams-Bashforth position integration |
mechanics/micromechanics/kernels/openmp_solver/src/motility_solver.cpp |
Cell motility with persistence and bias |
mechanics/micromechanics/kernels/openmp_solver/src/basement_membrane_solver.cpp |
Boundary repulsion forces |
mechanics/micromechanics/kernels/openmp_solver/src/spring_solver.cpp |
Spring attachment force calculation |
mechanics/micromechanics/kernels/openmp_solver/src/potentials/standard_potential.cpp |
PhysiCell-style repulsion-adhesion potential |
mechanics/micromechanics/kernels/openmp_solver/src/potentials/morse_potential.cpp |
Morse potential for smooth interactions |
mechanics/micromechanics/kernels/openmp_solver/src/potentials/kelvin_voigt_potential.cpp |
Spring-damper viscoelastic model |
mechanics/micromechanics/kernels/openmp_solver/src/register_solver.cpp |
Solver registration with global registry |
mechanics/micromechanics/tests/test_*.cpp |
Comprehensive unit tests for all components |
README.md |
Added mechanics module documentation and usage examples |
| #pragma omp parallel | ||
| { | ||
| // Each thread gets its own random generator | ||
| std::random_device rd; | ||
| std::mt19937 gen(rd()); |
There was a problem hiding this comment.
The std::random_device is being initialized inside the parallel region on each thread iteration, which is very expensive. Random device initialization should be moved outside the parallel region or use thread-local storage that persists across function calls. Consider:
#pragma omp parallel
{
// Use thread ID to seed instead
std::mt19937 gen(omp_get_thread_num() + std::random_device{}());
// ... rest of code
}This will significantly improve performance for motility calculations.
| #pragma omp parallel | |
| { | |
| // Each thread gets its own random generator | |
| std::random_device rd; | |
| std::mt19937 gen(rd()); | |
| // Generate seeds for each thread outside the parallel region | |
| int n_threads = 1; | |
| #ifdef _OPENMP | |
| n_threads = omp_get_max_threads(); | |
| #endif | |
| std::vector<std::uint32_t> seeds(n_threads); | |
| { | |
| std::random_device rd; | |
| for (int i = 0; i < n_threads; ++i) | |
| seeds[i] = rd(); | |
| } | |
| #pragma omp parallel | |
| { | |
| int thread_id = 0; | |
| #ifdef _OPENMP | |
| thread_id = omp_get_thread_num(); | |
| #endif | |
| // Each thread gets its own random generator seeded from the precomputed seeds | |
| std::mt19937 gen(seeds[thread_id]); |
| // So we use NEGATIVE force_coeff to get correct direction: | ||
| // repulsion (force_coeff > 0) → -force_coeff * dx → i moves away from j | ||
| // adhesion (force_coeff < 0) → -force_coeff * dx → i moves toward j | ||
| mech_data.forces[i * 3] -= force_coeff * dx; | ||
| mech_data.forces[i * 3 + 1] -= force_coeff * dy; | ||
| mech_data.forces[i * 3 + 2] -= force_coeff * dz; |
There was a problem hiding this comment.
[nitpick] The force accumulation in lines 130-132 uses subtraction with a negative sign (-= force_coeff * dx). While the comment explains the physics convention, this double negation pattern is error-prone. Consider inverting the sign convention at the potential calculation level to make force application more straightforward:
mech_data.forces[i * 3] += force_coeff * dx;This would require updating potential implementations to return forces with the opposite sign convention, but would make the force application code clearer.
| // So we use NEGATIVE force_coeff to get correct direction: | |
| // repulsion (force_coeff > 0) → -force_coeff * dx → i moves away from j | |
| // adhesion (force_coeff < 0) → -force_coeff * dx → i moves toward j | |
| mech_data.forces[i * 3] -= force_coeff * dx; | |
| mech_data.forces[i * 3 + 1] -= force_coeff * dy; | |
| mech_data.forces[i * 3 + 2] -= force_coeff * dz; | |
| // Now, force_coeff is signed so that: | |
| // repulsion (force_coeff > 0) → i moves away from j (opposite to dx) | |
| // adhesion (force_coeff < 0) → i moves toward j (same as dx) | |
| // So we accumulate using += for clarity: | |
| mech_data.forces[i * 3] += force_coeff * dx; | |
| mech_data.forces[i * 3 + 1] += force_coeff * dy; | |
| mech_data.forces[i * 3 + 2] += force_coeff * dz; |
README.md
Outdated
| agent->position() = { 0.0, 0.0, 0.0 }; | ||
| agent->radius() = 10.0; |
There was a problem hiding this comment.
The README example shows agent->position() = { 0.0, 0.0, 0.0 }; but position() returns a std::span<real_t> which cannot be assigned from an initializer list. The example should be:
auto pos = agent->position();
pos[0] = 0.0; pos[1] = 0.0; pos[2] = 0.0;or
std::array<real_t, 3> p = {0.0, 0.0, 0.0};
std::copy(p.begin(), p.end(), agent->position().begin());| // Calculate exp power: a * (1 - r^2/r0^2) | ||
| real_t exp_power = scaling_factor * (1.0 - (distance * distance) / (equilibrium_distance * equilibrium_distance)); |
There was a problem hiding this comment.
In the Morse potential force calculation, the implementation uses r^2 in the exponent calculation but this doesn't match the standard Morse potential formula which uses linear r. Line 46 should be:
real_t exp_power = scaling_factor * (1.0 - distance / equilibrium_distance);not
real_t exp_power = scaling_factor * (1.0 - (distance * distance) / (equilibrium_distance * equilibrium_distance));The current implementation creates a different potential shape than the documented V(r) = D * (exp(2a(1-r/r0)) - 2exp(a(1-r/r0))).
|
| * All properties are indexed by cell_id. Properties are computed | ||
| * by aggregating data from agents belonging to each cell. | ||
| */ | ||
| struct cell_data |
There was a problem hiding this comment.
can you explain why we need unordered maps to store these values?
| std::unordered_map<index_t, std::array<real_t, 3>> velocities; | ||
|
|
||
| /// Cell speed (velocity magnitude) indexed by cell_id. | ||
| std::unordered_map<index_t, real_t> speeds; |
|
|
||
| void run_single_timestep() override; | ||
| /// Cell-level data (pressure, etc.) aggregated from agents | ||
| cell_data cells; |
There was a problem hiding this comment.
I see that environment is holding just cell_data but not agent_data
There was a problem hiding this comment.
I updated micromechanics so environment.h no longer pulls in agent/container (and thus agent_data) details via the includes. All the changes are in commit ba2e953
| std::array<real_t, 3> domain_min = { -500.0, -500.0, -500.0 }; | ||
|
|
||
| /// Domain boundaries [x_max, y_max, z_max] | ||
| std::array<real_t, 3> domain_max = { 500.0, 500.0, 500.0 }; |
There was a problem hiding this comment.
do you think that these defaults are reasonable? if not, just assign them in the constructor
| * @param dz Z-component of normalized direction vector | ||
| * @param force_out Output: force magnitude (positive = repulsion, negative = attraction) | ||
| */ | ||
| virtual void calculate_pairwise_force(const environment& env, index_t agent_i, index_t agent_j, real_t distance, |
There was a problem hiding this comment.
why not returning the final force instead of passing as a parameter?
| * | ||
| * Adds config for both (type_a, type_b) and (type_b, type_a). | ||
| */ | ||
| void add_interaction(std::uint8_t type_a, std::uint8_t type_b, const interaction_config& config) |
| // Damping force: F_damp = gamma * (v_rel . n) * dt | ||
| // Get previous velocities | ||
| index_t dims = 3; | ||
| real_t dvx = mech_data.previous_velocities[agent_j * dims + 0] - mech_data.previous_velocities[agent_i * dims + 0]; | ||
| real_t dvy = mech_data.previous_velocities[agent_j * dims + 1] - mech_data.previous_velocities[agent_i * dims + 1]; | ||
| real_t dvz = mech_data.previous_velocities[agent_j * dims + 2] - mech_data.previous_velocities[agent_i * dims + 2]; |
There was a problem hiding this comment.
so do we want to include damping or not?
asmelko
left a comment
asmelko
left a comment
There was a problem hiding this comment.
I suggest dividing this PR into 3:
- serialization
- openmp_solver
- the rest that will remain this PR
README.md
Outdated
| ``` | ||
|
|
||
| The run generates `output/` with VTK files for visualization. | ||
|
|
There was a problem hiding this comment.
I would not put the readme here. lets keep it for github pages in the next PR
| target_link_libraries( | ||
| physicore.mechanics.micromechanics | ||
| PRIVATE physicore.mechanics.micromechanics_internal_iface | ||
| PUBLIC physicore::common) |
There was a problem hiding this comment.
is this part necessary? you are already including common in internal_iface
There was a problem hiding this comment.
In commit 8edff26 moved physicore::common / physicore::common::vtk_serializer off the tests-only micromechanics_internal_iface and onto the real library target in CMakeLists.txt.
| INTERFACE physicore.mechanics.micromechanics_internal_iface | ||
| physicore.mechanics.micromechanics_iface physicore::common) |
There was a problem hiding this comment.
internal_iface is ment to be internal, so only for testing.
There was a problem hiding this comment.
Changed in db5073f openmp_solver no longer links openmp_solver_internal_iface
| @@ -0,0 +1,3 @@ | |||
| # Tests for OpenMP micromechanics solver | |||
| # TODO: Add tests when test framework is set up | |||
There was a problem hiding this comment.
either include tests in this PR or move whole openmp_solver to a separate PR
| * Usage: | ||
| * auto solver = solver_registry::instance().get("openmp_solver"); | ||
| */ | ||
| class solver_registry |
There was a problem hiding this comment.
since we moved vtk serializer base to common, lets do the same with solver_registry
| return true; | ||
| }(); | ||
|
|
||
| class mock_solver : public solver |
There was a problem hiding this comment.
so you do not need to add other tests than AvailableSolvers.
If you move solver_registry to common, just move the general tests from biofvm under common tests and remove general tests from here
|
|
||
| class environment; | ||
|
|
||
| class spatial_index |
There was a problem hiding this comment.
maybe spacial index should also go to the kernel implementation?
| } | ||
| }; | ||
|
|
||
| class uniform_grid_spatial_index : public spatial_index |
There was a problem hiding this comment.
maybe spacial index should also go to the kernel implementation?
|
|
||
| using common::vtkRealArray; | ||
|
|
||
| class vtk_mechanics_serializer : public common::vtk_serializer_base, private generic_agent_solver<agent> |
There was a problem hiding this comment.
i see that you are deciding to serialize just a subset of data from agent_data. why so? I would leave the choice of serialized members to a different PR and remove it from here for now
There was a problem hiding this comment.
Chenged in commit 6403027 I made the mechanics VTK serializer write all flat agent_data fields. I will continue working on this on a different PR after the split
| simulation_parameters params; | ||
|
|
||
| /// The active solver (obtained from solver_registry) | ||
| std::unique_ptr<solver> active_solver; |
There was a problem hiding this comment.
maybe just solver is fine. calling it active_solver suggests that there is a time when solver is inactive which is not the case
marcorusc
force-pushed
the
dev/mech-micro
branch
from
f73e75d to
d9ae3a7
Compare
|
Split follow-ups:
|
Implements micromechanics solver registration via a header-only common factory registry and drops the old solver_registry_sole.cpp implementation.
Converts cell-level storage to flat SoA vectors and adds cell_generic_storage to provide cell-level accessors, updating tests accordingly.
Drops leftover speed references; speed is derived from velocity on demand. Includes a minimal OpenMP backend adjustment to avoid writing removed fields.
Moves environment domain_min/domain_max defaults from the header into the environment constructor.
|
Closes #26, #28 , #29 , #37 , #45 , #54