diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..c1a63dc --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,46 @@ +# Checks that run on every commit. They cover the fast, purely local gates: +# formatting and linting for each language in the tree, plus the Python type +# check. Tests, native builds, and backend conformance are not here; they need +# a configured build and hardware, and run through CTest and the conformance +# scripts instead. +repos: + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v5.0.0 + hooks: + - id: check-yaml + - id: check-toml + - id: check-merge-conflict + - id: check-added-large-files + # The device photomasks under docs/tutorials/devices are legitimately + # large and are tracked deliberately. + args: [--maxkb=1024] + exclude: ^docs/tutorials/devices/ + - id: end-of-file-fixer + exclude: ^docs/tutorials/devices/ + - id: trailing-whitespace + exclude: ^docs/tutorials/devices/ + + # Pinned to the version the project depends on, so a hook run and a local + # `ruff check` enforce the same rule set. Formatting is deliberately absent: + # the linter's rules, including import order, are what this tree follows. + - repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.16.3 + hooks: + - id: ruff + args: [--fix] + + - repo: https://github.com/pre-commit/mirrors-clang-format + rev: v19.1.7 + hooks: + - id: clang-format + types_or: [c++, c, cuda, objective-c++, metal] + files: ^cpp/ + + - repo: local + hooks: + - id: pyright + name: pyright + entry: uv run pyright + language: system + pass_filenames: false + types: [python] diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d6e6e3..4b820db 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -97,11 +97,12 @@ if(CM_ENABLE_METAL) add_custom_command( OUTPUT "${CM_METAL_COUPLED_RATES_HEADER}" COMMAND "${CMAKE_COMMAND}" - "-DINPUT=${CMAKE_CURRENT_SOURCE_DIR}/cpp/metal/kernels/coupled_rates.metal" + "-DINPUT=${CMAKE_CURRENT_SOURCE_DIR}/cpp/metal/kernels/grid_transport.metal;${CMAKE_CURRENT_SOURCE_DIR}/cpp/metal/kernels/coupled_rates.metal" "-DOUTPUT=${CM_METAL_COUPLED_RATES_HEADER}" "-DSYMBOL=coupled_rates_source" -P "${CMAKE_CURRENT_SOURCE_DIR}/cmake/EmbedMetalSource.cmake" DEPENDS + cpp/metal/kernels/grid_transport.metal cpp/metal/kernels/coupled_rates.metal cmake/EmbedMetalSource.cmake VERBATIM @@ -121,11 +122,12 @@ if(CM_ENABLE_METAL) add_custom_command( OUTPUT "${CM_METAL_SIGNALS_HEADER}" COMMAND "${CMAKE_COMMAND}" - "-DINPUT=${CMAKE_CURRENT_SOURCE_DIR}/cpp/metal/kernels/signals.metal" + "-DINPUT=${CMAKE_CURRENT_SOURCE_DIR}/cpp/metal/kernels/grid_transport.metal;${CMAKE_CURRENT_SOURCE_DIR}/cpp/metal/kernels/signals.metal" "-DOUTPUT=${CM_METAL_SIGNALS_HEADER}" "-DSYMBOL=signals_source" -P "${CMAKE_CURRENT_SOURCE_DIR}/cmake/EmbedMetalSource.cmake" DEPENDS + cpp/metal/kernels/grid_transport.metal cpp/metal/kernels/signals.metal cmake/EmbedMetalSource.cmake VERBATIM diff --git a/cmake/EmbedMetalSource.cmake b/cmake/EmbedMetalSource.cmake index 7d903f9..9fb8fd7 100644 --- a/cmake/EmbedMetalSource.cmake +++ b/cmake/EmbedMetalSource.cmake @@ -1,3 +1,6 @@ +# INPUT names one Metal source, or several to concatenate in order. A Metal +# library is compiled from source at runtime with no include path, so a source +# that shares helpers with another receives them by concatenation here. if(NOT DEFINED INPUT OR NOT DEFINED OUTPUT OR NOT DEFINED SYMBOL) message(FATAL_ERROR "EmbedMetalSource.cmake requires INPUT, OUTPUT, and SYMBOL") endif() @@ -5,11 +8,16 @@ if(NOT SYMBOL MATCHES "^[A-Za-z_][A-Za-z0-9_]*$") message(FATAL_ERROR "EmbedMetalSource.cmake received an invalid C++ symbol") endif() -file(READ "${INPUT}" CM_METAL_SOURCE) get_filename_component(CM_OUTPUT_DIRECTORY "${OUTPUT}" DIRECTORY) file(MAKE_DIRECTORY "${CM_OUTPUT_DIRECTORY}") file(WRITE "${OUTPUT}" "#pragma once\n\nnamespace cm::metal {\ninline constexpr char ${SYMBOL}[] = R\"CM_METAL(") -file(APPEND "${OUTPUT}" "${CM_METAL_SOURCE}") +foreach(CM_METAL_INPUT IN LISTS INPUT) + file(READ "${CM_METAL_INPUT}" CM_METAL_SOURCE) + if(CM_METAL_SOURCE MATCHES "CM_METAL\\(" OR CM_METAL_SOURCE MATCHES "\\)CM_METAL") + message(FATAL_ERROR "Metal source ${CM_METAL_INPUT} contains the embedding delimiter") + endif() + file(APPEND "${OUTPUT}" "${CM_METAL_SOURCE}") +endforeach() file(APPEND "${OUTPUT}" [=[)CM_METAL"; } // namespace cm::metal ]=]) diff --git a/cpp/core/constraints.cpp b/cpp/core/constraints.cpp index 45c28f0..498b2cf 100644 --- a/cpp/core/constraints.cpp +++ b/cpp/core/constraints.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -44,9 +45,8 @@ void validate_sphere(const SphereConstraint& sphere) { } bool positive_finite_extents(const Vec3& half_extents) { - return std::isfinite(half_extents.x) && half_extents.x > 0.0F && - std::isfinite(half_extents.y) && half_extents.y > 0.0F && - std::isfinite(half_extents.z) && half_extents.z > 0.0F; + return std::isfinite(half_extents.x) && half_extents.x > 0.0F && std::isfinite(half_extents.y) && + half_extents.y > 0.0F && std::isfinite(half_extents.z) && half_extents.z > 0.0F; } void validate_box(const BoxConstraint& box) { @@ -94,50 +94,25 @@ void validate_constraint_state(ConstraintId next_id, std::span ids; ids.reserve(total); - ConstraintId previous_plane = invalid_constraint_id; - for (const auto& plane : planes) { - validate_plane(plane); - if (plane.id <= previous_plane || plane.id >= next_id) { - throw std::invalid_argument("checkpoint plane identifiers are not ordered and allocated"); + const auto check_ordered = [&ids, next_id](const auto& constraints, auto&& validate, + const char* kind) { + ConstraintId previous = invalid_constraint_id; + for (const auto& constraint : constraints) { + validate(constraint); + if (constraint.id <= previous || constraint.id >= next_id) { + throw std::invalid_argument(std::string("checkpoint ") + kind + + " identifiers are not ordered and allocated"); + } + if (!ids.insert(constraint.id).second) { + throw std::invalid_argument("checkpoint contains a duplicate constraint identifier"); + } + previous = constraint.id; } - if (!ids.insert(plane.id).second) { - throw std::invalid_argument("checkpoint contains a duplicate constraint identifier"); - } - previous_plane = plane.id; - } - ConstraintId previous_sphere = invalid_constraint_id; - for (const auto& sphere : spheres) { - validate_sphere(sphere); - if (sphere.id <= previous_sphere || sphere.id >= next_id) { - throw std::invalid_argument("checkpoint sphere identifiers are not ordered and allocated"); - } - if (!ids.insert(sphere.id).second) { - throw std::invalid_argument("checkpoint contains a duplicate constraint identifier"); - } - previous_sphere = sphere.id; - } - ConstraintId previous_box = invalid_constraint_id; - for (const auto& box : boxes) { - validate_box(box); - if (box.id <= previous_box || box.id >= next_id) { - throw std::invalid_argument("checkpoint box identifiers are not ordered and allocated"); - } - if (!ids.insert(box.id).second) { - throw std::invalid_argument("checkpoint contains a duplicate constraint identifier"); - } - previous_box = box.id; - } - ConstraintId previous_cylinder = invalid_constraint_id; - for (const auto& cylinder : cylinders) { - validate_cylinder(cylinder); - if (cylinder.id <= previous_cylinder || cylinder.id >= next_id) { - throw std::invalid_argument("checkpoint cylinder identifiers are not ordered and allocated"); - } - if (!ids.insert(cylinder.id).second) { - throw std::invalid_argument("checkpoint contains a duplicate constraint identifier"); - } - previous_cylinder = cylinder.id; - } + }; + check_ordered(planes, validate_plane, "plane"); + check_ordered(spheres, validate_sphere, "sphere"); + check_ordered(boxes, validate_box, "box"); + check_ordered(cylinders, validate_cylinder, "cylinder"); } std::size_t checked_offset_count(std::size_t cell_count) { diff --git a/cpp/core/contact_graph.cpp b/cpp/core/contact_graph.cpp index 22ee480..1b12b31 100644 --- a/cpp/core/contact_graph.cpp +++ b/cpp/core/contact_graph.cpp @@ -70,8 +70,8 @@ void validate_contact_parameters(const ContactParameters& parameters) { } } -std::vector find_cell_contact_candidates( - const WorldState& state, const ContactParameters& parameters) { +std::vector find_cell_contact_candidates(const WorldState& state, + const ContactParameters& parameters) { validate_contact_parameters(parameters); const auto geometry = state.geometry_state(); std::vector bounds; @@ -112,10 +112,10 @@ std::vector find_cell_contact_candidates( }); active.erase(expired.begin(), expired.end()); for (const auto* candidate : active) { - const auto overlaps_y = candidate->maximum_y >= current.minimum_y && - current.maximum_y >= candidate->minimum_y; - const auto overlaps_z = candidate->maximum_z >= current.minimum_z && - current.maximum_z >= candidate->minimum_z; + const auto overlaps_y = + candidate->maximum_y >= current.minimum_y && current.maximum_y >= candidate->minimum_y; + const auto overlaps_z = + candidate->maximum_z >= current.minimum_z && current.maximum_z >= candidate->minimum_z; if (!overlaps_y || !overlaps_z) { continue; } @@ -125,11 +125,11 @@ std::vector find_cell_contact_candidates( } active.push_back(¤t); } - std::ranges::sort(candidates, [&geometry](const ContactCandidate& left, - const ContactCandidate& right) { - return std::tuple{geometry.ids[left.first_slot], geometry.ids[left.second_slot]} < - std::tuple{geometry.ids[right.first_slot], geometry.ids[right.second_slot]}; - }); + std::ranges::sort( + candidates, [&geometry](const ContactCandidate& left, const ContactCandidate& right) { + return std::tuple{geometry.ids[left.first_slot], geometry.ids[left.second_slot]} < + std::tuple{geometry.ids[right.first_slot], geometry.ids[right.second_slot]}; + }); return candidates; } diff --git a/cpp/core/mechanics_integration.cpp b/cpp/core/mechanics_integration.cpp index 25e4a6e..a13be1d 100644 --- a/cpp/core/mechanics_integration.cpp +++ b/cpp/core/mechanics_integration.cpp @@ -20,6 +20,8 @@ bool finite(const Vec3& value) { return std::isfinite(value.x) && std::isfinite(value.y) && std::isfinite(value.z); } +} // namespace + Vec3 rotate_axis_angle(Vec3 direction, Vec3 rotation, float max_rotation) { const auto magnitude = norm(rotation); if (magnitude <= 1.0e-12F || max_rotation == 0.0F) { @@ -33,8 +35,6 @@ Vec3 rotate_axis_angle(Vec3 direction, Vec3 rotation, float max_rotation) { axis * (dot(axis, direction) * (1.0F - cosine))); } -} // namespace - void validate_mechanics_integration_parameters(const MechanicsIntegrationParameters& parameters) { if (!std::isfinite(parameters.max_rotation_radians) || parameters.max_rotation_radians < 0.0F) { throw std::invalid_argument("mechanics rotation limit must be finite and non-negative"); @@ -74,10 +74,9 @@ void integrate_mechanics_result(WorldState& state, const MechanicsSolveResult& r const auto applied_length_increment = cell.fixed ? desired_increment : std::max(0.0F, desired_increment + correction.length); const auto new_position = cell.fixed ? cell.position : cell.position + correction.translation; - const auto new_direction = - cell.fixed ? cell.direction - : rotate_axis_angle(cell.direction, correction.rotation, - parameters.max_rotation_radians); + const auto new_direction = cell.fixed ? cell.direction + : rotate_axis_angle(cell.direction, correction.rotation, + parameters.max_rotation_radians); const auto new_length = cell.length + applied_length_increment; if (!finite(new_position) || !finite(new_direction) || !std::isfinite(new_length)) { throw std::overflow_error("mechanics integration produced non-finite geometry"); diff --git a/cpp/core/signals.cpp b/cpp/core/signals.cpp index c96c014..1b5cd5c 100644 --- a/cpp/core/signals.cpp +++ b/cpp/core/signals.cpp @@ -39,7 +39,7 @@ struct AxisWeights { }; AxisWeights interpolation_axis(float position, float origin, float spacing, std::uint32_t dimension, - const char* axis) { + const char* axis, GridSampleBound bound) { if (!std::isfinite(position)) { throw std::invalid_argument("signal sample position must be finite"); } @@ -47,10 +47,14 @@ AxisWeights interpolation_axis(float position, float origin, float spacing, std: return {}; } - const auto coordinate = + auto coordinate = (static_cast(position) - static_cast(origin)) / static_cast(spacing); const auto upper_bound = static_cast(dimension - 1); - if (coordinate < 0.0 || coordinate > upper_bound) { + if (bound == GridSampleBound::clamped) { + // Clamping happens in lattice coordinates, the same space the bound is + // tested in, so a clamped position always lands inside the lattice. + coordinate = std::clamp(coordinate, 0.0, upper_bound); + } else if (coordinate < 0.0 || coordinate > upper_bound) { throw std::out_of_range(std::string("signal sample is outside the ") + axis + " grid bound"); } const auto lower = static_cast(std::floor(coordinate)); @@ -242,10 +246,9 @@ double max_reaction_loss(const SignalGridSpec& spec, std::size_t signal) { return 0.0; } const auto sites = spec.site_count(); - const auto begin = spec.reaction->loss_rates.begin() + - static_cast(signal * sites); - return static_cast( - *std::max_element(begin, begin + static_cast(sites))); + const auto begin = + spec.reaction->loss_rates.begin() + static_cast(signal * sites); + return static_cast(*std::max_element(begin, begin + static_cast(sites))); } float rms(std::span values) { @@ -261,11 +264,15 @@ float rms(std::span values) { } // namespace -SignalGridStencil signal_grid_stencil(const SignalGridSpec& spec, Vec3 position) { - spec.validate(); - const auto x = interpolation_axis(position.x, spec.origin.x, spec.spacing.x, spec.shape.x, "x"); - const auto y = interpolation_axis(position.y, spec.origin.y, spec.spacing.y, spec.shape.y, "y"); - const auto z = interpolation_axis(position.z, spec.origin.z, spec.spacing.z, spec.shape.z, "z"); +SignalGridStencil signal_grid_stencil(const SignalGridSpec& spec, Vec3 position, + GridSampleBound bound) { + spec.validate_lattice(); + const auto x = + interpolation_axis(position.x, spec.origin.x, spec.spacing.x, spec.shape.x, "x", bound); + const auto y = + interpolation_axis(position.y, spec.origin.y, spec.spacing.y, spec.shape.y, "y", bound); + const auto z = + interpolation_axis(position.z, spec.origin.z, spec.spacing.z, spec.shape.z, "z", bound); SignalGridStencil result; for (std::size_t xi = 0; xi < x.count; ++xi) { for (std::size_t yi = 0; yi < y.count; ++yi) { @@ -290,7 +297,8 @@ SignalGridStencil signal_grid_stencil(const SignalGridSpec& spec, Vec3 position) } if (dropped) { if (fluid_weight <= 0.0F) { - throw std::invalid_argument("signal sample position is inside a grid obstacle"); + result.entirely_solid = true; + return result; } for (std::size_t entry = 0; entry < result.count; ++entry) { result.weights[entry] /= fluid_weight; @@ -384,7 +392,7 @@ bool SignalGridSpec::solid_site(std::size_t site) const noexcept { return !obstacles.empty() && obstacles[site] != 0; } -void SignalGridSpec::validate() const { +void SignalGridSpec::validate_lattice() const { if (signal_count == 0) { throw std::invalid_argument("signal grid must contain at least one signal"); } @@ -397,6 +405,20 @@ void SignalGridSpec::validate() const { if (!finite(spacing) || spacing.x <= 0.0F || spacing.y <= 0.0F || spacing.z <= 0.0F) { throw std::invalid_argument("signal grid spacing must be finite and positive"); } + if (!obstacles.empty() && obstacles.size() != site_count()) { + throw std::invalid_argument("signal grid obstacle mask must cover every site"); + } + if (velocity_field.has_value()) { + const auto& field = *velocity_field; + if (field.x_faces.size() != x_face_count() || field.y_faces.size() != y_face_count() || + field.z_faces.size() != z_face_count()) { + throw std::invalid_argument("signal grid velocity field must cover every lattice face"); + } + } +} + +void SignalGridSpec::validate() const { + validate_lattice(); if (diffusion.size() != signal_count || advection.size() != signal_count) { throw std::invalid_argument("signal grid transport arrays must match signal count"); } @@ -414,9 +436,6 @@ void SignalGridSpec::validate() const { reaction->validate(level_count()); } if (!obstacles.empty()) { - if (obstacles.size() != site_count()) { - throw std::invalid_argument("signal grid obstacle mask must cover every site"); - } for (const auto value : obstacles) { if (value > 1) { throw std::invalid_argument("signal grid obstacle mask values must be 0 or 1"); @@ -426,9 +445,8 @@ void SignalGridSpec::validate() const { const auto sites = site_count(); for (std::size_t signal = 0; signal < signal_count; ++signal) { for (std::size_t site = 0; site < sites; ++site) { - if (obstacles[site] != 0 && - (reaction->source_rates[(signal * sites) + site] != 0.0F || - reaction->loss_rates[(signal * sites) + site] != 0.0F)) { + if (obstacles[site] != 0 && (reaction->source_rates[(signal * sites) + site] != 0.0F || + reaction->loss_rates[(signal * sites) + site] != 0.0F)) { throw std::invalid_argument( "signal grid affine reaction must be zero at obstacle sites"); } @@ -438,10 +456,6 @@ void SignalGridSpec::validate() const { } if (velocity_field.has_value()) { const auto& field = *velocity_field; - if (field.x_faces.size() != x_face_count() || field.y_faces.size() != y_face_count() || - field.z_faces.size() != z_face_count()) { - throw std::invalid_argument("signal grid velocity field must cover every lattice face"); - } for (const auto* faces : {&field.x_faces, &field.y_faces, &field.z_faces}) { for (const auto value : *faces) { if (!std::isfinite(value)) { @@ -572,6 +586,9 @@ std::span SignalGrid::levels() const& noexcept { return levels_; } std::vector SignalGrid::sample(Vec3 position) const { const auto stencil = signal_grid_stencil(spec_, position); + if (stencil.entirely_solid) { + throw std::invalid_argument("signal sample position is inside a grid obstacle"); + } const auto sites = spec_.site_count(); std::vector result(spec_.signal_count, 0.0F); for (std::size_t entry = 0; entry < stencil.count; ++entry) { @@ -582,6 +599,41 @@ std::vector SignalGrid::sample(Vec3 position) const { return result; } +Vec3 SignalGrid::sample_velocity(Vec3 position, GridSampleBound bound) const { + if (!spec_.velocity_field.has_value()) { + throw std::logic_error("signal grid does not declare a velocity field"); + } + const auto stencil = signal_grid_stencil(spec_, position, bound); + // The field is zero on every face of a solid site, so a stencil with no + // fluid in it samples zero: a cell that mechanics has pressed into a wall + // does not drift. + const auto& field = *spec_.velocity_field; + Vec3 result{}; + for (std::size_t entry = 0; entry < stencil.count; ++entry) { + const auto site = stencil.sites[entry]; + const auto weight = stencil.weights[entry]; + if (weight == 0.0F) { + continue; + } + const auto z = site % spec_.shape.z; + const auto y = (site / spec_.shape.z) % spec_.shape.y; + const auto x = site / (static_cast(spec_.shape.y) * spec_.shape.z); + const auto fx = static_cast(x); + const auto fy = static_cast(y); + const auto fz = static_cast(z); + result.x += weight * 0.5F * + (field.x_faces[x_face_index(spec_.shape, fx, fy, fz)] + + field.x_faces[x_face_index(spec_.shape, fx + 1, fy, fz)]); + result.y += weight * 0.5F * + (field.y_faces[y_face_index(spec_.shape, fx, fy, fz)] + + field.y_faces[y_face_index(spec_.shape, fx, fy + 1, fz)]); + result.z += weight * 0.5F * + (field.z_faces[z_face_index(spec_.shape, fx, fy, fz)] + + field.z_faces[z_face_index(spec_.shape, fx, fy, fz + 1)]); + } + return result; +} + SignalGridCheckpoint SignalGrid::checkpoint() const { validate(); return {.spec = spec_, .levels = levels_}; @@ -604,6 +656,13 @@ void SignalGrid::set_velocity_field(std::optional field spec_ = std::move(candidate.spec); } +void SignalGrid::set_reaction(std::optional reaction) { + SignalGridCheckpoint candidate{.spec = spec_, .levels = levels_}; + candidate.spec.reaction = std::move(reaction); + candidate.validate(); + spec_ = std::move(candidate.spec); +} + void SignalGrid::validate_step(float dt) const { if (!std::isfinite(dt) || dt < 0.0F) { throw std::invalid_argument("time step must be finite and non-negative"); @@ -647,10 +706,9 @@ void SignalGrid::validate_step(float dt) const { inverse_square_sum += inverse_spacing * inverse_spacing; courant_sum += std::abs(static_cast(velocity[axis])) * inverse_spacing; } - const auto factor = - static_cast(dt) * - ((2.0 * static_cast(spec_.diffusion[signal]) * inverse_square_sum) + courant_sum + - max_reaction_loss(spec_, signal)); + const auto factor = static_cast(dt) * + ((2.0 * static_cast(spec_.diffusion[signal]) * inverse_square_sum) + + courant_sum + max_reaction_loss(spec_, signal)); if (!std::isfinite(factor) || factor > 1.0) { throw std::invalid_argument("signal grid time step violates the explicit stability bound"); } @@ -792,8 +850,19 @@ SignalSolveResult signal_grid_crank_nicolson_candidate(const SignalGrid& grid, f const auto source = source_rates.empty() ? 0.0F : source_rates[index]; right_hand_side[index] = old[index] + (half_dt * old_rates[index]) + (dt * source); } - const auto threshold = - spec.solver.absolute_tolerance + (spec.solver.relative_tolerance * rms(right_hand_side)); + // The relative term scales the residual the step starts with, not the field + // it starts from. A field's own magnitude says nothing about how much of it + // this step has to change, so scaling by the field lets a small source fall + // under the threshold and be discarded; scaling by the initial residual asks + // for a fixed reduction of whatever this step actually has to resolve. + // + // A residual cannot fall below what float32 can represent for a field of + // this magnitude, and the right-hand side carries both the field and the + // operator terms of the step, so it sets that floor. An absolute tolerance + // asking for less than the floor is raised to it rather than making the + // solve unreachable. + auto threshold = std::max(spec.solver.absolute_tolerance, + std::numeric_limits::epsilon() * rms(right_hand_side)); std::vector current(old.begin(), old.end()); std::vector residual(old.size()); auto residual_rms = std::numeric_limits::infinity(); @@ -808,6 +877,9 @@ SignalSolveResult signal_grid_crank_nicolson_candidate(const SignalGrid& grid, f if (!std::isfinite(residual_rms)) { break; } + if (iterations == 0) { + threshold += spec.solver.relative_tolerance * residual_rms; + } if (residual_rms <= threshold) { return { .levels = std::move(current), diff --git a/cpp/core/simulation.cpp b/cpp/core/simulation.cpp index e96f76e..eef5a77 100644 --- a/cpp/core/simulation.cpp +++ b/cpp/core/simulation.cpp @@ -1,5 +1,6 @@ #include "cm/simulation.hpp" +#include #include #include @@ -102,6 +103,63 @@ CellId Simulation::add_cell(const CellInit& cell) { return state_.add_cell(cell) void Simulation::remove_cell(CellId id) { state_.remove_cell(id); } +void Simulation::apply_flow_drift(float dt, const MechanicsIntegrationParameters& parameters) { + if (!std::isfinite(dt) || dt < 0.0F) { + throw std::invalid_argument("time step must be finite and non-negative"); + } + validate_mechanics_integration_parameters(parameters); + if (!signal_grid_.has_value() || !signal_grid_->spec().velocity_field.has_value()) { + throw std::logic_error("flow drift requires a signal grid with a velocity field"); + } + if (dt == 0.0F || state_.empty()) { + return; + } + constexpr float degeneracy_epsilon = 1.0e-6F; + const auto geometry = state_.geometry_state(); + const auto attributes = state_.cell_attributes(); + struct DriftUpdate { + Slot slot; + Vec3 position; + Vec3 direction; + float length; + }; + std::vector updates; + updates.reserve(geometry.size()); + for (std::size_t slot = 0; slot < geometry.size(); ++slot) { + if (attributes.fixed[slot] != 0) { + continue; + } + const Vec3 center{geometry.position_x[slot], geometry.position_y[slot], + geometry.position_z[slot]}; + const Vec3 axis{geometry.direction_x[slot], geometry.direction_y[slot], + geometry.direction_z[slot]}; + const auto half_length = geometry.lengths[slot] * 0.5F; + // Rod endpoints may poke past the lattice of site centers (the mechanics + // walls, not the lattice edge, bound cells), so drift samples the nearest + // in-lattice point instead of erroring. + const auto first_velocity = + signal_grid_->sample_velocity(center - axis * half_length, GridSampleBound::clamped); + const auto second_velocity = + signal_grid_->sample_velocity(center + axis * half_length, GridSampleBound::clamped); + const auto mean_velocity = (first_velocity + second_velocity) * 0.5F; + const auto position = center + mean_velocity * dt; + auto direction = axis; + if (geometry.lengths[slot] > degeneracy_epsilon) { + const auto rotation = + cross(axis, (second_velocity - first_velocity) * (dt / geometry.lengths[slot])); + direction = rotate_axis_angle(axis, rotation, parameters.max_rotation_radians); + } + if (!std::isfinite(position.x) || !std::isfinite(position.y) || !std::isfinite(position.z) || + !std::isfinite(direction.x) || !std::isfinite(direction.y) || !std::isfinite(direction.z)) { + throw std::runtime_error("flow drift produced non-finite geometry"); + } + updates.push_back({static_cast(slot), position, direction, geometry.lengths[slot]}); + } + for (const auto& update : updates) { + state_.set_cell_geometry(update.slot, update.position, update.direction, update.length); + } +} + ConstraintId Simulation::add_plane_constraint(const PlaneConstraintInit& plane) { return constraints_.add_plane(plane); } @@ -178,6 +236,13 @@ void Simulation::set_velocity_field(std::optional field signal_grid_->set_velocity_field(std::move(field)); } +void Simulation::set_signal_reaction(std::optional reaction) { + if (!signal_grid_.has_value()) { + throw std::logic_error("simulation does not have a signal grid"); + } + signal_grid_->set_reaction(std::move(reaction)); +} + std::pair Simulation::divide(CellId parent_id, float first_fraction) { return state_.divide(parent_id, first_fraction); } diff --git a/cpp/core/world_state.cpp b/cpp/core/world_state.cpp index ea98ae1..a3a2ae1 100644 --- a/cpp/core/world_state.cpp +++ b/cpp/core/world_state.cpp @@ -475,7 +475,7 @@ void WorldState::validate() const { const std::array sizes{ position_x_.size(), position_y_.size(), position_z_.size(), direction_x_.size(), direction_y_.size(), direction_z_.size(), length_.size(), radius_.size(), - growth_rate_.size(), cell_type_.size(), fixed_.size(), + growth_rate_.size(), cell_type_.size(), fixed_.size(), }; if (!std::ranges::all_of(sizes, [expected](std::size_t size) { return size == expected; })) { throw std::logic_error("world state arrays have inconsistent lengths"); diff --git a/cpp/cpu/cpu_contacts.cpp b/cpp/cpu/cpu_contacts.cpp index 1094874..a4ce142 100644 --- a/cpp/cpu/cpu_contacts.cpp +++ b/cpp/cpu/cpu_contacts.cpp @@ -155,8 +155,7 @@ Vec3 deterministic_normal(const Capsule& first, const Capsule& second, const Poi return normalized(cross(first.axis, *least_aligned)); } -ContactGraph contacts_for_candidates(const WorldState& state, - const ContactParameters& parameters, +ContactGraph contacts_for_candidates(const WorldState& state, const ContactParameters& parameters, std::span candidates) { const auto geometry = state.geometry_state(); std::vector contacts; @@ -216,11 +215,10 @@ ContactGraph find_cell_contacts_cpu_exhaustive(const WorldState& state, candidates.reserve(pair_count); for (std::size_t first = 0; first < geometry.size(); ++first) { for (std::size_t second = first + 1; second < geometry.size(); ++second) { - candidates.push_back(geometry.ids[first] < geometry.ids[second] - ? ContactCandidate{static_cast(first), - static_cast(second)} - : ContactCandidate{static_cast(second), - static_cast(first)}); + candidates.push_back( + geometry.ids[first] < geometry.ids[second] + ? ContactCandidate{static_cast(first), static_cast(second)} + : ContactCandidate{static_cast(second), static_cast(first)}); } } return contacts_for_candidates(state, parameters, candidates); diff --git a/cpp/cpu/cpu_coupled.cpp b/cpp/cpu/cpu_coupled.cpp index 0b73526..85c1967 100644 --- a/cpp/cpu/cpu_coupled.cpp +++ b/cpp/cpu/cpu_coupled.cpp @@ -106,9 +106,12 @@ void validate_coupled_step(const WorldState& state, const SignalGrid& grid, } const auto geometry = state.geometry_state(); for (std::size_t cell = 0; cell < state.size(); ++cell) { - static_cast(signal_grid_stencil( + const auto stencil = signal_grid_stencil( grid.spec(), - {geometry.position_x[cell], geometry.position_y[cell], geometry.position_z[cell]})); + {geometry.position_x[cell], geometry.position_y[cell], geometry.position_z[cell]}); + if (stencil.entirely_solid) { + throw std::invalid_argument("signal sample position is inside a grid obstacle"); + } } } diff --git a/cpp/cuda/cuda_backend.cu b/cpp/cuda/cuda_backend.cu index ab14c98..7bbb579 100644 --- a/cpp/cuda/cuda_backend.cu +++ b/cpp/cuda/cuda_backend.cu @@ -340,8 +340,7 @@ class CudaBackend final : public ComputeBackend { copy_to_device(signal_z_faces_, spec.velocity_field->z_faces, "failed to upload CUDA signal z faces"); } - const auto has_velocity_field = - static_cast(spec.velocity_field.has_value()); + const auto has_velocity_field = static_cast(spec.velocity_field.has_value()); check_cuda(cudaMemsetAsync(signal_error_.data(), 0, sizeof(std::uint32_t), stream_), "failed to clear the CUDA signal-grid error flag"); @@ -365,9 +364,8 @@ class CudaBackend final : public ComputeBackend { signal_advection_.data(), signal_fixed_values_.data(), signal_reaction_source_.data(), signal_reaction_loss_.data(), signal_obstacles_.data(), signal_x_faces_.data(), signal_y_faces_.data(), signal_z_faces_.data(), has_velocity_field, signal_error_.data(), - boundaries, shape, - make_float4(spec.spacing.x, spec.spacing.y, spec.spacing.z, 0.0F), dt, signal_count, - level_count, crank_nicolson, stream_); + boundaries, shape, make_float4(spec.spacing.x, spec.spacing.y, spec.spacing.z, 0.0F), dt, + signal_count, level_count, crank_nicolson, stream_); check_cuda(cudaGetLastError(), "failed to launch the CUDA signal-grid kernel"); std::uint32_t error = 0; @@ -386,10 +384,9 @@ class CudaBackend final : public ComputeBackend { signal_levels_.data(), signal_output_.data(), signal_diffusion_.data(), signal_advection_.data(), signal_fixed_values_.data(), signal_reaction_source_.data(), signal_reaction_loss_.data(), signal_obstacles_.data(), signal_x_faces_.data(), - signal_y_faces_.data(), signal_z_faces_.data(), has_velocity_field, - signal_error_.data(), boundaries, shape, - make_float4(spec.spacing.x, spec.spacing.y, spec.spacing.z, 0.0F), dt, signal_count, - level_count, spec.solver); + signal_y_faces_.data(), signal_z_faces_.data(), has_velocity_field, signal_error_.data(), + boundaries, shape, make_float4(spec.spacing.x, spec.spacing.y, spec.spacing.z, 0.0F), dt, + signal_count, level_count, spec.solver); result_device = solve.first; report = solve.second; if (!report.converged) { @@ -562,8 +559,7 @@ class CudaBackend final : public ComputeBackend { copy_to_device(coupled_z_faces_, spec.velocity_field->z_faces, "failed to upload CUDA coupled z faces"); } - const auto has_velocity_field = - static_cast(spec.velocity_field.has_value()); + const auto has_velocity_field = static_cast(spec.velocity_field.has_value()); check_cuda(cudaMemsetAsync(coupled_error_.data(), 0, sizeof(std::uint32_t), stream_), "failed to clear the CUDA coupled error flag"); @@ -976,8 +972,8 @@ class CudaBackend final : public ComputeBackend { .id = cylinder.id, .kind = static_cast(ExternalConstraintKind::cylinder), .allowed_region = static_cast(cylinder.allowed_region), - .geometry = make_float4(cylinder.center.x, cylinder.center.y, cylinder.center.z, - cylinder.radius), + .geometry = + make_float4(cylinder.center.x, cylinder.center.y, cylinder.center.z, cylinder.radius), .parameters = make_float4(cylinder.half_height, 0.0F, 0.0F, cylinder.coefficient), }); } @@ -1153,8 +1149,7 @@ class CudaBackend final : public ComputeBackend { std::vector contacts; contacts.reserve(contact_count); for (std::uint32_t index = 0; index < contact_count; ++index) { - if (constraint_kinds[index] > - static_cast(ExternalConstraintKind::cylinder) || + if (constraint_kinds[index] > static_cast(ExternalConstraintKind::cylinder) || locations[index] > static_cast(RodContactLocation::interior)) { throw std::runtime_error("CUDA external-contact kernel produced an invalid tag"); } @@ -1347,25 +1342,13 @@ class CudaBackend final : public ComputeBackend { return result; } - [[nodiscard]] float signal_rhs_rms(const float* right_hand_side, std::uint32_t level_count) { - cuda::launch_signal_square_terms(right_hand_side, signal_cn_terms_.data(), level_count, - stream_); - check_cuda(cudaGetLastError(), "failed to launch the CUDA signal-norm kernel"); - return std::sqrt(reduce_signal_terms(level_count, "CUDA signal norm failed") / - static_cast(level_count)); - } - - [[nodiscard]] float signal_residual_rms(const float* current, const float* right_hand_side, - const float* diffusion, const float4* advection, - const float* fixed_values, const float* reaction_source, - const float* reaction_loss, - const std::uint8_t* obstacles, const float* x_faces, - const float* y_faces, const float* z_faces, - std::uint32_t has_velocity_field, - cuda::SignalGridBoundariesGpu boundaries, - cuda::SignalGridShapeGpu shape, float4 spacing, - float half_dt, std::uint32_t signal_count, - std::uint32_t level_count) { + [[nodiscard]] float signal_residual_rms( + const float* current, const float* right_hand_side, const float* diffusion, + const float4* advection, const float* fixed_values, const float* reaction_source, + const float* reaction_loss, const std::uint8_t* obstacles, const float* x_faces, + const float* y_faces, const float* z_faces, std::uint32_t has_velocity_field, + cuda::SignalGridBoundariesGpu boundaries, cuda::SignalGridShapeGpu shape, float4 spacing, + float half_dt, std::uint32_t signal_count, std::uint32_t level_count) { cuda::launch_signal_crank_nicolson_residual_terms( current, right_hand_side, signal_cn_terms_.data(), diffusion, advection, fixed_values, reaction_source, reaction_loss, obstacles, x_faces, y_faces, z_faces, has_velocity_field, @@ -1375,6 +1358,14 @@ class CudaBackend final : public ComputeBackend { static_cast(level_count)); } + [[nodiscard]] float signal_rhs_rms(const float* right_hand_side, std::uint32_t level_count) { + cuda::launch_signal_square_terms(right_hand_side, signal_cn_terms_.data(), level_count, + stream_); + check_cuda(cudaGetLastError(), "failed to launch the CUDA signal-norm kernel"); + return std::sqrt(reduce_signal_terms(level_count, "CUDA signal norm failed") / + static_cast(level_count)); + } + [[nodiscard]] std::pair solve_signal_crank_nicolson( const float* initial, const float* right_hand_side, const float* diffusion, const float4* advection, const float* fixed_values, const float* reaction_source, @@ -1385,14 +1376,23 @@ class CudaBackend final : public ComputeBackend { std::uint32_t level_count, const SignalSolveParameters& parameters) { ensure_signal_solve_capacity(level_count); const auto half_dt = 0.5F * dt; - const auto right_hand_side_rms = signal_rhs_rms(right_hand_side, level_count); - const auto threshold = - parameters.absolute_tolerance + (parameters.relative_tolerance * right_hand_side_rms); SignalSolveReport report; report.residual_rms = signal_residual_rms( initial, right_hand_side, diffusion, advection, fixed_values, reaction_source, - reaction_loss, obstacles, x_faces, y_faces, z_faces, has_velocity_field, boundaries, - shape, spacing, half_dt, signal_count, level_count); + reaction_loss, obstacles, x_faces, y_faces, z_faces, has_velocity_field, boundaries, shape, + spacing, half_dt, signal_count, level_count); + // The relative term scales the residual the step starts with, matching the + // CPU reference: the field's own magnitude says nothing about how much of + // it this step has to change. + // A residual cannot fall below what float32 can represent for a field of + // this magnitude, and the right-hand side carries both the field and the + // operator terms of the step, so it sets that floor. An absolute tolerance + // asking for less than the floor is raised to it rather than making the + // solve unreachable. + const auto floor = + std::numeric_limits::epsilon() * signal_rhs_rms(right_hand_side, level_count); + const auto threshold = std::max(parameters.absolute_tolerance, floor) + + (parameters.relative_tolerance * report.residual_rms); if (std::isfinite(report.residual_rms) && report.residual_rms <= threshold) { return {initial, report}; } @@ -1406,12 +1406,10 @@ class CudaBackend final : public ComputeBackend { const float* current = initial; for (std::uint32_t iteration = 1; iteration <= parameters.max_iterations; ++iteration) { float* output = current == signal_cn_a_.data() ? signal_cn_b_.data() : signal_cn_a_.data(); - cuda::launch_signal_crank_nicolson_jacobi(current, output, right_hand_side, diffusion, - advection, fixed_values, reaction_source, - reaction_loss, obstacles, x_faces, y_faces, - z_faces, has_velocity_field, error, boundaries, - shape, spacing, half_dt, signal_count, - level_count, stream_); + cuda::launch_signal_crank_nicolson_jacobi( + current, output, right_hand_side, diffusion, advection, fixed_values, reaction_source, + reaction_loss, obstacles, x_faces, y_faces, z_faces, has_velocity_field, error, + boundaries, shape, spacing, half_dt, signal_count, level_count, stream_); check_cuda(cudaGetLastError(), "failed to launch the CUDA signal Jacobi kernel"); report.residual_rms = signal_residual_rms( output, right_hand_side, diffusion, advection, fixed_values, reaction_source, diff --git a/cpp/cuda/kernels/contacts.cu b/cpp/cuda/kernels/contacts.cu index f7142fb..40a6b9a 100644 --- a/cpp/cuda/kernels/contacts.cu +++ b/cpp/cuda/kernels/contacts.cu @@ -528,13 +528,15 @@ __global__ void inclusive_scan_step(const std::uint32_t* input, std::uint32_t* o output[index] = value; } -__global__ void fill_cell_contacts( - const std::uint64_t* ids, const float4* centers, const float4* axes, const float4* geometry, - const uint2* candidates, const std::uint32_t* counts, - const std::uint32_t* inclusive_counts, std::uint64_t* first_ids, - std::uint64_t* second_ids, std::uint32_t* first_slots, std::uint32_t* second_slots, - std::uint32_t* ordinals, float4* points_on_first, float4* normals, float* separations, - float* weights, ContactParametersGpu parameters, std::uint32_t candidate_count) { +__global__ void fill_cell_contacts(const std::uint64_t* ids, const float4* centers, + const float4* axes, const float4* geometry, + const uint2* candidates, const std::uint32_t* counts, + const std::uint32_t* inclusive_counts, std::uint64_t* first_ids, + std::uint64_t* second_ids, std::uint32_t* first_slots, + std::uint32_t* second_slots, std::uint32_t* ordinals, + float4* points_on_first, float4* normals, float* separations, + float* weights, ContactParametersGpu parameters, + std::uint32_t candidate_count) { const auto pair_index = blockIdx.x * blockDim.x + threadIdx.x; if (pair_index >= candidate_count) { return; @@ -656,12 +658,11 @@ void launch_inclusive_scan_step(const std::uint32_t* input, std::uint32_t* outpu void launch_contact_fill(const std::uint64_t* ids, const float4* centers, const float4* axes, const float4* geometry, const uint2* candidates, - const std::uint32_t* counts, - const std::uint32_t* inclusive_counts, std::uint64_t* first_ids, - std::uint64_t* second_ids, std::uint32_t* first_slots, - std::uint32_t* second_slots, std::uint32_t* ordinals, - float4* points_on_first, float4* normals, float* separations, - float* weights, ContactParametersGpu parameters, + const std::uint32_t* counts, const std::uint32_t* inclusive_counts, + std::uint64_t* first_ids, std::uint64_t* second_ids, + std::uint32_t* first_slots, std::uint32_t* second_slots, + std::uint32_t* ordinals, float4* points_on_first, float4* normals, + float* separations, float* weights, ContactParametersGpu parameters, std::uint32_t candidate_count, cudaStream_t stream) { constexpr std::uint32_t threads = 256; const auto blocks = ((candidate_count - 1) / threads) + 1; diff --git a/cpp/cuda/kernels/contacts.cuh b/cpp/cuda/kernels/contacts.cuh index 5ed8205..c2d08a5 100644 --- a/cpp/cuda/kernels/contacts.cuh +++ b/cpp/cuda/kernels/contacts.cuh @@ -38,12 +38,11 @@ void launch_inclusive_scan_step(const std::uint32_t* input, std::uint32_t* outpu void launch_contact_fill(const std::uint64_t* ids, const float4* centers, const float4* axes, const float4* geometry, const uint2* candidates, - const std::uint32_t* counts, - const std::uint32_t* inclusive_counts, std::uint64_t* first_ids, - std::uint64_t* second_ids, std::uint32_t* first_slots, - std::uint32_t* second_slots, std::uint32_t* ordinals, - float4* points_on_first, float4* normals, float* separations, - float* weights, ContactParametersGpu parameters, + const std::uint32_t* counts, const std::uint32_t* inclusive_counts, + std::uint64_t* first_ids, std::uint64_t* second_ids, + std::uint32_t* first_slots, std::uint32_t* second_slots, + std::uint32_t* ordinals, float4* points_on_first, float4* normals, + float* separations, float* weights, ContactParametersGpu parameters, std::uint32_t candidate_count, cudaStream_t stream); void launch_external_contact_count(const std::uint64_t* ids, const float4* centers, diff --git a/cpp/cuda/kernels/coupled_rates.cu b/cpp/cuda/kernels/coupled_rates.cu index 04d7f73..9f720d4 100644 --- a/cpp/cuda/kernels/coupled_rates.cu +++ b/cpp/cuda/kernels/coupled_rates.cu @@ -8,16 +8,6 @@ namespace { constexpr float pi = 3.14159265358979323846F; constexpr std::uint32_t threads_per_block = 256; -__device__ std::uint32_t site_index(SignalGridShapeGpu shape, std::uint32_t x, std::uint32_t y, - std::uint32_t z) { - return x * shape.y * shape.z + y * shape.z + z; -} - -__device__ float grid_level(const float* levels, SignalGridShapeGpu shape, std::uint32_t signal, - std::uint32_t x, std::uint32_t y, std::uint32_t z) { - return levels[signal * shape.sites + site_index(shape, x, y, z)]; -} - __device__ float effective_volume(float length, float radius) { return pi * radius * radius * (length + 2.0F * radius); } @@ -98,6 +88,8 @@ __device__ float sample_signal(const float* levels, SignalGridShapeGpu shape, fl } } } + // A stencil with no fluid corner is rejected by the host's coupled-step + // validation before any kernel runs, so the fluid weight is positive here. if (dropped) { result /= fluid_weight; } @@ -107,10 +99,16 @@ __device__ float sample_signal(const float* levels, SignalGridShapeGpu shape, fl __device__ float cell_scatter_weight(float4 center, SignalGridShapeGpu shape, float4 origin, float4 spacing, const std::uint8_t* obstacles, std::uint32_t x, std::uint32_t y, std::uint32_t z) { + // A cell only scatters into the eight sites of its own stencil, and the + // weight is pure arithmetic, so testing it first keeps the obstacle mask out + // of the sites a cell cannot reach - which is nearly all of them. + const auto raw = cell_site_weight(center, shape, origin, spacing, x, y, z); + if (raw == 0.0F) { + return 0.0F; + } if (obstacles[site_index(shape, x, y, z)] != 0) { return 0.0F; } - const auto raw = cell_site_weight(center, shape, origin, spacing, x, y, z); const auto coordinate_x = axis_coordinate(center.x, origin.x, spacing.x, shape.x); const auto coordinate_y = axis_coordinate(center.y, origin.y, spacing.y, shape.y); const auto coordinate_z = axis_coordinate(center.z, origin.z, spacing.z, shape.z); @@ -142,6 +140,8 @@ __device__ float cell_scatter_weight(float4 center, SignalGridShapeGpu shape, fl } } } + // A stencil with no fluid corner is rejected by the host's coupled-step + // validation before any kernel runs, so the fluid weight is positive here. return dropped ? raw / fluid_weight : raw; } @@ -217,9 +217,8 @@ __global__ void advance_coupled_cells( const RateInstructionGpu* instructions, const std::uint32_t* species_outputs, const std::uint32_t* signal_outputs, float* workspace, const float* grid_levels, float* cell_signal_rates, const std::uint8_t* obstacles, std::uint32_t* error, - SignalGridShapeGpu shape, float4 origin, float4 spacing, float dt, - std::uint32_t species_count, std::uint32_t signal_count, std::uint32_t instruction_count, - std::uint32_t cell_count) { + SignalGridShapeGpu shape, float4 origin, float4 spacing, float dt, std::uint32_t species_count, + std::uint32_t signal_count, std::uint32_t instruction_count, std::uint32_t cell_count) { const auto cell = (blockIdx.x * blockDim.x) + threadIdx.x; if (cell >= cell_count) { return; @@ -265,29 +264,14 @@ __global__ void advance_coupled_cells( } } -__device__ float exterior_value(std::uint32_t kind, const float* fixed_values, std::uint32_t face, - std::uint32_t signal, std::uint32_t signal_count, float current, - float periodic) { - if (kind == 0) { - return current; - } - if (kind == 1) { - return periodic; - } - return fixed_values[face * signal_count + signal]; -} - -__global__ void advance_coupled_grid(const float* levels, float* output, const float* diffusion, - const float4* advection, const float* fixed_values, - const float* reaction_source, const float* reaction_loss, - const float4* centers, const float* cell_signal_rates, - const std::uint8_t* obstacles, const float* x_faces, - const float* y_faces, const float* z_faces, - std::uint32_t has_velocity_field, std::uint32_t* error, - SignalGridBoundariesGpu boundaries, SignalGridShapeGpu shape, - float4 origin, float4 spacing, float dt, - std::uint32_t signal_count, std::uint32_t cell_count, - std::uint32_t level_count, bool crank_nicolson) { +__global__ void advance_coupled_grid( + const float* levels, float* output, const float* diffusion, const float4* advection, + const float* fixed_values, const float* reaction_source, const float* reaction_loss, + const float4* centers, const float* cell_signal_rates, const std::uint8_t* obstacles, + const float* x_faces, const float* y_faces, const float* z_faces, + std::uint32_t has_velocity_field, std::uint32_t* error, SignalGridBoundariesGpu boundaries, + SignalGridShapeGpu shape, float4 origin, float4 spacing, float dt, std::uint32_t signal_count, + std::uint32_t cell_count, std::uint32_t level_count, bool crank_nicolson) { const auto index = (blockIdx.x * blockDim.x) + threadIdx.x; if (index >= level_count) { return; @@ -330,51 +314,12 @@ __global__ void advance_coupled_grid(const float* levels, float* output, const f : grid_level(levels, shape, signal, x, y, z + 1); const std::uint32_t dimensions[3]{shape.x, shape.y, shape.z}; - bool closed_lower[3]; - bool closed_upper[3]; - closed_lower[0] = - x == 0 ? (boundaries.x_lower == 0 || - (boundaries.x_lower == 1 && obstacles[site_index(shape, shape.x - 1, y, z)] != 0)) - : obstacles[site_index(shape, x - 1, y, z)] != 0; - closed_upper[0] = - x + 1 == shape.x - ? (boundaries.x_upper == 0 || - (boundaries.x_upper == 1 && obstacles[site_index(shape, 0, y, z)] != 0)) - : obstacles[site_index(shape, x + 1, y, z)] != 0; - closed_lower[1] = - y == 0 ? (boundaries.y_lower == 0 || - (boundaries.y_lower == 1 && obstacles[site_index(shape, x, shape.y - 1, z)] != 0)) - : obstacles[site_index(shape, x, y - 1, z)] != 0; - closed_upper[1] = - y + 1 == shape.y - ? (boundaries.y_upper == 0 || - (boundaries.y_upper == 1 && obstacles[site_index(shape, x, 0, z)] != 0)) - : obstacles[site_index(shape, x, y + 1, z)] != 0; - closed_lower[2] = - z == 0 ? (boundaries.z_lower == 0 || - (boundaries.z_lower == 1 && obstacles[site_index(shape, x, y, shape.z - 1)] != 0)) - : obstacles[site_index(shape, x, y, z - 1)] != 0; - closed_upper[2] = - z + 1 == shape.z - ? (boundaries.z_upper == 0 || - (boundaries.z_upper == 1 && obstacles[site_index(shape, x, y, 0)] != 0)) - : obstacles[site_index(shape, x, y, z + 1)] != 0; - float face_lower[3]; - float face_upper[3]; - if (has_velocity_field != 0) { - face_lower[0] = x_faces[x * shape.y * shape.z + y * shape.z + z]; - face_upper[0] = x_faces[(x + 1) * shape.y * shape.z + y * shape.z + z]; - face_lower[1] = y_faces[x * (shape.y + 1) * shape.z + y * shape.z + z]; - face_upper[1] = y_faces[x * (shape.y + 1) * shape.z + (y + 1) * shape.z + z]; - face_lower[2] = z_faces[x * shape.y * (shape.z + 1) + y * (shape.z + 1) + z]; - face_upper[2] = z_faces[x * shape.y * (shape.z + 1) + y * (shape.z + 1) + z + 1]; - } else { - const float velocity[3]{advection[signal].x, advection[signal].y, advection[signal].z}; - for (std::uint32_t axis = 0; axis < 3; ++axis) { - face_lower[axis] = velocity[axis]; - face_upper[axis] = velocity[axis]; - } - } + const auto faces = grid_face_state(shape, boundaries, obstacles, x_faces, y_faces, z_faces, + has_velocity_field, advection[signal], x, y, z); + const bool* closed_lower = faces.closed_lower; + const bool* closed_upper = faces.closed_upper; + const float* face_lower = faces.lower; + const float* face_upper = faces.upper; const float grid_spacing[3]{spacing.x, spacing.y, spacing.z}; float rate = 0.0F; for (std::uint32_t axis = 0; axis < 3; ++axis) { @@ -390,10 +335,10 @@ __global__ void advance_coupled_grid(const float* levels, float* output, const f const auto inverse_spacing = 1.0F / grid_spacing[axis]; rate += diffusion[signal] * (lower[axis] - 2.0F * current + upper[axis]) * inverse_spacing * inverse_spacing; - auto lower_flux = face_lower[axis] >= 0.0F ? face_lower[axis] * lower[axis] - : face_lower[axis] * current; - auto upper_flux = face_upper[axis] >= 0.0F ? face_upper[axis] * current - : face_upper[axis] * upper[axis]; + auto lower_flux = + face_lower[axis] >= 0.0F ? face_lower[axis] * lower[axis] : face_lower[axis] * current; + auto upper_flux = + face_upper[axis] >= 0.0F ? face_upper[axis] * current : face_upper[axis] * upper[axis]; if (closed_lower[axis]) { lower_flux = 0.0F; } @@ -407,7 +352,8 @@ __global__ void advance_coupled_grid(const float* levels, float* output, const f float source = 0.0F; const auto inverse_voxel_volume = 1.0F / (spacing.x * spacing.y * spacing.z); for (std::uint32_t cell = 0; cell < cell_count; ++cell) { - const auto weight = cell_scatter_weight(centers[cell], shape, origin, spacing, obstacles, x, y, z); + const auto weight = + cell_scatter_weight(centers[cell], shape, origin, spacing, obstacles, x, y, z); source += weight * cell_signal_rates[cell * signal_count + signal] * inverse_voxel_volume; } const auto transport_scale = crank_nicolson ? 0.5F * dt : dt; diff --git a/cpp/cuda/kernels/mechanics.cu b/cpp/cuda/kernels/mechanics.cu index ccfbf46..ae5536c 100644 --- a/cpp/cuda/kernels/mechanics.cu +++ b/cpp/cuda/kernels/mechanics.cu @@ -181,8 +181,7 @@ __global__ void add_mechanics_regularizer(const float4* axes, const float4* geom __global__ void initialize_mechanics_vectors(MechanicsDofsGpu* right_hand_side, MechanicsDofsGpu* solution, MechanicsDofsGpu* residual, MechanicsDofsGpu* search_direction, - const std::uint8_t* fixed, - std::uint32_t cell_count) { + const std::uint8_t* fixed, std::uint32_t cell_count) { const auto cell = blockIdx.x * blockDim.x + threadIdx.x; if (cell >= cell_count) { return; diff --git a/cpp/cuda/kernels/signals.cu b/cpp/cuda/kernels/signals.cu index 3c37eed..541cfd8 100644 --- a/cpp/cuda/kernels/signals.cu +++ b/cpp/cuda/kernels/signals.cu @@ -5,42 +5,17 @@ namespace cm::cuda { namespace { -__device__ std::uint32_t site_index(SignalGridShapeGpu shape, std::uint32_t x, std::uint32_t y, - std::uint32_t z) { - return x * shape.y * shape.z + y * shape.z + z; -} - -__device__ float grid_level(const float* levels, SignalGridShapeGpu shape, std::uint32_t signal, - std::uint32_t x, std::uint32_t y, std::uint32_t z) { - return levels[signal * shape.sites + site_index(shape, x, y, z)]; -} - -__device__ float exterior_value(std::uint32_t kind, const float* fixed_values, std::uint32_t face, - std::uint32_t signal, std::uint32_t signal_count, float current, - float periodic) { - if (kind == 0) { - return current; - } - if (kind == 1) { - return periodic; - } - return fixed_values[face * signal_count + signal]; -} - struct TransportPoint { float rate; float diagonal; }; -__device__ TransportPoint transport_point(const float* levels, const float* diffusion, - const float4* advection, const float* fixed_values, - const float* reaction_source, const float* reaction_loss, - const std::uint8_t* obstacles, const float* x_faces, - const float* y_faces, const float* z_faces, - std::uint32_t has_velocity_field, - SignalGridBoundariesGpu boundaries, - SignalGridShapeGpu shape, float4 spacing, - std::uint32_t signal_count, std::uint32_t index) { +__device__ TransportPoint transport_point( + const float* levels, const float* diffusion, const float4* advection, const float* fixed_values, + const float* reaction_source, const float* reaction_loss, const std::uint8_t* obstacles, + const float* x_faces, const float* y_faces, const float* z_faces, + std::uint32_t has_velocity_field, SignalGridBoundariesGpu boundaries, SignalGridShapeGpu shape, + float4 spacing, std::uint32_t signal_count, std::uint32_t index) { const auto signal = index / shape.sites; const auto site = index - signal * shape.sites; const auto x = site / (shape.y * shape.z); @@ -77,51 +52,12 @@ __device__ TransportPoint transport_point(const float* levels, const float* diff : grid_level(levels, shape, signal, x, y, z + 1); const std::uint32_t dimensions[3]{shape.x, shape.y, shape.z}; - bool closed_lower[3]; - bool closed_upper[3]; - closed_lower[0] = - x == 0 ? (boundaries.x_lower == 0 || - (boundaries.x_lower == 1 && obstacles[site_index(shape, shape.x - 1, y, z)] != 0)) - : obstacles[site_index(shape, x - 1, y, z)] != 0; - closed_upper[0] = - x + 1 == shape.x - ? (boundaries.x_upper == 0 || - (boundaries.x_upper == 1 && obstacles[site_index(shape, 0, y, z)] != 0)) - : obstacles[site_index(shape, x + 1, y, z)] != 0; - closed_lower[1] = - y == 0 ? (boundaries.y_lower == 0 || - (boundaries.y_lower == 1 && obstacles[site_index(shape, x, shape.y - 1, z)] != 0)) - : obstacles[site_index(shape, x, y - 1, z)] != 0; - closed_upper[1] = - y + 1 == shape.y - ? (boundaries.y_upper == 0 || - (boundaries.y_upper == 1 && obstacles[site_index(shape, x, 0, z)] != 0)) - : obstacles[site_index(shape, x, y + 1, z)] != 0; - closed_lower[2] = - z == 0 ? (boundaries.z_lower == 0 || - (boundaries.z_lower == 1 && obstacles[site_index(shape, x, y, shape.z - 1)] != 0)) - : obstacles[site_index(shape, x, y, z - 1)] != 0; - closed_upper[2] = - z + 1 == shape.z - ? (boundaries.z_upper == 0 || - (boundaries.z_upper == 1 && obstacles[site_index(shape, x, y, 0)] != 0)) - : obstacles[site_index(shape, x, y, z + 1)] != 0; - float face_lower[3]; - float face_upper[3]; - if (has_velocity_field != 0) { - face_lower[0] = x_faces[x * shape.y * shape.z + y * shape.z + z]; - face_upper[0] = x_faces[(x + 1) * shape.y * shape.z + y * shape.z + z]; - face_lower[1] = y_faces[x * (shape.y + 1) * shape.z + y * shape.z + z]; - face_upper[1] = y_faces[x * (shape.y + 1) * shape.z + (y + 1) * shape.z + z]; - face_lower[2] = z_faces[x * shape.y * (shape.z + 1) + y * (shape.z + 1) + z]; - face_upper[2] = z_faces[x * shape.y * (shape.z + 1) + y * (shape.z + 1) + z + 1]; - } else { - const float velocity[3]{advection[signal].x, advection[signal].y, advection[signal].z}; - for (std::uint32_t axis = 0; axis < 3; ++axis) { - face_lower[axis] = velocity[axis]; - face_upper[axis] = velocity[axis]; - } - } + const auto faces = grid_face_state(shape, boundaries, obstacles, x_faces, y_faces, z_faces, + has_velocity_field, advection[signal], x, y, z); + const bool* closed_lower = faces.closed_lower; + const bool* closed_upper = faces.closed_upper; + const float* face_lower = faces.lower; + const float* face_upper = faces.upper; const float grid_spacing[3]{spacing.x, spacing.y, spacing.z}; float rate = 0.0F; float diagonal = 0.0F; @@ -145,10 +81,10 @@ __device__ TransportPoint transport_point(const float* levels, const float* diff if (closed_upper[axis]) { diagonal += diffusion_scale; } - auto lower_flux = face_lower[axis] >= 0.0F ? face_lower[axis] * lower[axis] - : face_lower[axis] * current; - auto upper_flux = face_upper[axis] >= 0.0F ? face_upper[axis] * current - : face_upper[axis] * upper[axis]; + auto lower_flux = + face_lower[axis] >= 0.0F ? face_lower[axis] * lower[axis] : face_lower[axis] * current; + auto upper_flux = + face_upper[axis] >= 0.0F ? face_upper[axis] * current : face_upper[axis] * upper[axis]; if (closed_lower[axis]) { lower_flux = 0.0F; } @@ -174,10 +110,9 @@ __global__ void advance_signal_grid(const float* levels, float* output, const fl const std::uint8_t* obstacles, const float* x_faces, const float* y_faces, const float* z_faces, std::uint32_t has_velocity_field, std::uint32_t* error, - SignalGridBoundariesGpu boundaries, - SignalGridShapeGpu shape, float4 spacing, float dt, - std::uint32_t signal_count, std::uint32_t level_count, - bool crank_nicolson) { + SignalGridBoundariesGpu boundaries, SignalGridShapeGpu shape, + float4 spacing, float dt, std::uint32_t signal_count, + std::uint32_t level_count, bool crank_nicolson) { const auto index = (blockIdx.x * blockDim.x) + threadIdx.x; if (index >= level_count) { return; @@ -207,9 +142,8 @@ __global__ void signal_crank_nicolson_jacobi( const float4* advection, const float* fixed_values, const float* reaction_source, const float* reaction_loss, const std::uint8_t* obstacles, const float* x_faces, const float* y_faces, const float* z_faces, std::uint32_t has_velocity_field, - std::uint32_t* error, SignalGridBoundariesGpu boundaries, - SignalGridShapeGpu shape, float4 spacing, float half_dt, std::uint32_t signal_count, - std::uint32_t level_count) { + std::uint32_t* error, SignalGridBoundariesGpu boundaries, SignalGridShapeGpu shape, + float4 spacing, float half_dt, std::uint32_t signal_count, std::uint32_t level_count) { const auto index = (blockIdx.x * blockDim.x) + threadIdx.x; if (index >= level_count) { return; @@ -232,8 +166,8 @@ __global__ void signal_crank_nicolson_residual_terms( const float4* advection, const float* fixed_values, const float* reaction_source, const float* reaction_loss, const std::uint8_t* obstacles, const float* x_faces, const float* y_faces, const float* z_faces, std::uint32_t has_velocity_field, - SignalGridBoundariesGpu boundaries, SignalGridShapeGpu shape, - float4 spacing, float half_dt, std::uint32_t signal_count, std::uint32_t level_count) { + SignalGridBoundariesGpu boundaries, SignalGridShapeGpu shape, float4 spacing, float half_dt, + std::uint32_t signal_count, std::uint32_t level_count) { const auto index = (blockIdx.x * blockDim.x) + threadIdx.x; if (index >= level_count) { return; @@ -248,22 +182,19 @@ __global__ void signal_crank_nicolson_residual_terms( } // namespace -void launch_advance_signal_grid(const float* levels, float* output, const float* diffusion, - const float4* advection, const float* fixed_values, - const float* reaction_source, const float* reaction_loss, - const std::uint8_t* obstacles, const float* x_faces, - const float* y_faces, const float* z_faces, - std::uint32_t has_velocity_field, std::uint32_t* error, - SignalGridBoundariesGpu boundaries, SignalGridShapeGpu shape, - float4 spacing, float dt, std::uint32_t signal_count, - std::uint32_t level_count, bool crank_nicolson, - cudaStream_t stream) { +void launch_advance_signal_grid( + const float* levels, float* output, const float* diffusion, const float4* advection, + const float* fixed_values, const float* reaction_source, const float* reaction_loss, + const std::uint8_t* obstacles, const float* x_faces, const float* y_faces, const float* z_faces, + std::uint32_t has_velocity_field, std::uint32_t* error, SignalGridBoundariesGpu boundaries, + SignalGridShapeGpu shape, float4 spacing, float dt, std::uint32_t signal_count, + std::uint32_t level_count, bool crank_nicolson, cudaStream_t stream) { constexpr std::uint32_t threads_per_block = 256; const auto block_count = ((level_count - 1) / threads_per_block) + 1; advance_signal_grid<<>>( - levels, output, diffusion, advection, fixed_values, reaction_source, reaction_loss, - obstacles, x_faces, y_faces, z_faces, has_velocity_field, error, boundaries, shape, spacing, - dt, signal_count, level_count, crank_nicolson); + levels, output, diffusion, advection, fixed_values, reaction_source, reaction_loss, obstacles, + x_faces, y_faces, z_faces, has_velocity_field, error, boundaries, shape, spacing, dt, + signal_count, level_count, crank_nicolson); } void launch_signal_square_terms(const float* input, float* terms, std::uint32_t level_count, @@ -273,17 +204,14 @@ void launch_signal_square_terms(const float* input, float* terms, std::uint32_t signal_square_terms<<>>(input, terms, level_count); } -void launch_signal_crank_nicolson_jacobi(const float* current, float* output, - const float* right_hand_side, const float* diffusion, - const float4* advection, const float* fixed_values, - const float* reaction_source, const float* reaction_loss, - const std::uint8_t* obstacles, const float* x_faces, - const float* y_faces, const float* z_faces, - std::uint32_t has_velocity_field, std::uint32_t* error, - SignalGridBoundariesGpu boundaries, - SignalGridShapeGpu shape, float4 spacing, float half_dt, - std::uint32_t signal_count, std::uint32_t level_count, - cudaStream_t stream) { +void launch_signal_crank_nicolson_jacobi( + const float* current, float* output, const float* right_hand_side, const float* diffusion, + const float4* advection, const float* fixed_values, const float* reaction_source, + const float* reaction_loss, const std::uint8_t* obstacles, const float* x_faces, + const float* y_faces, const float* z_faces, std::uint32_t has_velocity_field, + std::uint32_t* error, SignalGridBoundariesGpu boundaries, SignalGridShapeGpu shape, + float4 spacing, float half_dt, std::uint32_t signal_count, std::uint32_t level_count, + cudaStream_t stream) { constexpr std::uint32_t threads_per_block = 256; const auto block_count = ((level_count - 1) / threads_per_block) + 1; signal_crank_nicolson_jacobi<<>>( diff --git a/cpp/cuda/kernels/signals.cuh b/cpp/cuda/kernels/signals.cuh index 8109a5e..088995f 100644 --- a/cpp/cuda/kernels/signals.cuh +++ b/cpp/cuda/kernels/signals.cuh @@ -22,30 +22,109 @@ struct SignalGridBoundariesGpu { std::uint32_t z_upper; }; -void launch_advance_signal_grid(const float* levels, float* output, const float* diffusion, - const float4* advection, const float* fixed_values, - const float* reaction_source, const float* reaction_loss, - const std::uint8_t* obstacles, const float* x_faces, - const float* y_faces, const float* z_faces, - std::uint32_t has_velocity_field, std::uint32_t* error, - SignalGridBoundariesGpu boundaries, SignalGridShapeGpu shape, - float4 spacing, float dt, std::uint32_t signal_count, - std::uint32_t level_count, bool crank_nicolson, - cudaStream_t stream); +// Grid geometry and transport helpers shared by the signal and coupled-rate +// kernels. Whole-program compilation gives each translation unit its own copy, +// so a device-inline definition in this header needs no device linking. + +__device__ inline std::uint32_t site_index(SignalGridShapeGpu shape, std::uint32_t x, + std::uint32_t y, std::uint32_t z) { + return x * shape.y * shape.z + y * shape.z + z; +} + +__device__ inline float grid_level(const float* levels, SignalGridShapeGpu shape, + std::uint32_t signal, std::uint32_t x, std::uint32_t y, + std::uint32_t z) { + return levels[signal * shape.sites + site_index(shape, x, y, z)]; +} + +__device__ inline float exterior_value(std::uint32_t kind, const float* fixed_values, + std::uint32_t face, std::uint32_t signal, + std::uint32_t signal_count, float current, float periodic) { + if (kind == 0) { + return current; + } + if (kind == 1) { + return periodic; + } + return fixed_values[face * signal_count + signal]; +} + +// Whether each of a site's six faces is closed to transport, and the velocity +// it carries. A face is closed by a no-flux boundary, by a periodic boundary +// that wraps onto a solid site, or by a solid neighbour. +struct GridFaceState { + bool closed_lower[3]; + bool closed_upper[3]; + float lower[3]; + float upper[3]; +}; + +__device__ inline GridFaceState grid_face_state(SignalGridShapeGpu shape, + SignalGridBoundariesGpu boundaries, + const std::uint8_t* obstacles, const float* x_faces, + const float* y_faces, const float* z_faces, + std::uint32_t has_velocity_field, float4 advection, + std::uint32_t x, std::uint32_t y, std::uint32_t z) { + GridFaceState faces{}; + faces.closed_lower[0] = + x == 0 ? (boundaries.x_lower == 0 || + (boundaries.x_lower == 1 && obstacles[site_index(shape, shape.x - 1, y, z)] != 0)) + : obstacles[site_index(shape, x - 1, y, z)] != 0; + faces.closed_upper[0] = + x + 1 == shape.x ? (boundaries.x_upper == 0 || + (boundaries.x_upper == 1 && obstacles[site_index(shape, 0, y, z)] != 0)) + : obstacles[site_index(shape, x + 1, y, z)] != 0; + faces.closed_lower[1] = + y == 0 ? (boundaries.y_lower == 0 || + (boundaries.y_lower == 1 && obstacles[site_index(shape, x, shape.y - 1, z)] != 0)) + : obstacles[site_index(shape, x, y - 1, z)] != 0; + faces.closed_upper[1] = + y + 1 == shape.y ? (boundaries.y_upper == 0 || + (boundaries.y_upper == 1 && obstacles[site_index(shape, x, 0, z)] != 0)) + : obstacles[site_index(shape, x, y + 1, z)] != 0; + faces.closed_lower[2] = + z == 0 ? (boundaries.z_lower == 0 || + (boundaries.z_lower == 1 && obstacles[site_index(shape, x, y, shape.z - 1)] != 0)) + : obstacles[site_index(shape, x, y, z - 1)] != 0; + faces.closed_upper[2] = + z + 1 == shape.z ? (boundaries.z_upper == 0 || + (boundaries.z_upper == 1 && obstacles[site_index(shape, x, y, 0)] != 0)) + : obstacles[site_index(shape, x, y, z + 1)] != 0; + if (has_velocity_field != 0) { + faces.lower[0] = x_faces[x * shape.y * shape.z + y * shape.z + z]; + faces.upper[0] = x_faces[(x + 1) * shape.y * shape.z + y * shape.z + z]; + faces.lower[1] = y_faces[x * (shape.y + 1) * shape.z + y * shape.z + z]; + faces.upper[1] = y_faces[x * (shape.y + 1) * shape.z + (y + 1) * shape.z + z]; + faces.lower[2] = z_faces[x * shape.y * (shape.z + 1) + y * (shape.z + 1) + z]; + faces.upper[2] = z_faces[x * shape.y * (shape.z + 1) + y * (shape.z + 1) + z + 1]; + } else { + const float velocity[3]{advection.x, advection.y, advection.z}; + for (std::uint32_t axis = 0; axis < 3; ++axis) { + faces.lower[axis] = velocity[axis]; + faces.upper[axis] = velocity[axis]; + } + } + return faces; +} + +void launch_advance_signal_grid( + const float* levels, float* output, const float* diffusion, const float4* advection, + const float* fixed_values, const float* reaction_source, const float* reaction_loss, + const std::uint8_t* obstacles, const float* x_faces, const float* y_faces, const float* z_faces, + std::uint32_t has_velocity_field, std::uint32_t* error, SignalGridBoundariesGpu boundaries, + SignalGridShapeGpu shape, float4 spacing, float dt, std::uint32_t signal_count, + std::uint32_t level_count, bool crank_nicolson, cudaStream_t stream); void launch_signal_square_terms(const float* input, float* terms, std::uint32_t level_count, cudaStream_t stream); -void launch_signal_crank_nicolson_jacobi(const float* current, float* output, - const float* right_hand_side, const float* diffusion, - const float4* advection, const float* fixed_values, - const float* reaction_source, const float* reaction_loss, - const std::uint8_t* obstacles, const float* x_faces, - const float* y_faces, const float* z_faces, - std::uint32_t has_velocity_field, std::uint32_t* error, - SignalGridBoundariesGpu boundaries, - SignalGridShapeGpu shape, float4 spacing, float half_dt, - std::uint32_t signal_count, std::uint32_t level_count, - cudaStream_t stream); +void launch_signal_crank_nicolson_jacobi( + const float* current, float* output, const float* right_hand_side, const float* diffusion, + const float4* advection, const float* fixed_values, const float* reaction_source, + const float* reaction_loss, const std::uint8_t* obstacles, const float* x_faces, + const float* y_faces, const float* z_faces, std::uint32_t has_velocity_field, + std::uint32_t* error, SignalGridBoundariesGpu boundaries, SignalGridShapeGpu shape, + float4 spacing, float half_dt, std::uint32_t signal_count, std::uint32_t level_count, + cudaStream_t stream); void launch_signal_crank_nicolson_residual_terms( const float* current, const float* right_hand_side, float* terms, const float* diffusion, const float4* advection, const float* fixed_values, const float* reaction_source, diff --git a/cpp/include/cm/mechanics.hpp b/cpp/include/cm/mechanics.hpp index 4abf99a..a920c40 100644 --- a/cpp/include/cm/mechanics.hpp +++ b/cpp/include/cm/mechanics.hpp @@ -58,6 +58,10 @@ struct MechanicsIntegrationParameters { void validate_mechanics_parameters(const MechanicsParameters& parameters); void validate_mechanics_integration_parameters(const MechanicsIntegrationParameters& parameters); +// The axis-angle rotation of a unit direction by a rotation vector, with the +// angle capped at ``max_rotation``. +[[nodiscard]] Vec3 rotate_axis_angle(Vec3 direction, Vec3 rotation, float max_rotation); + void integrate_mechanics_result( WorldState& state, const MechanicsSolveResult& result, const MechanicsIntegrationParameters& parameters = MechanicsIntegrationParameters{}, diff --git a/cpp/include/cm/signals.hpp b/cpp/include/cm/signals.hpp index f08d802..4c13251 100644 --- a/cpp/include/cm/signals.hpp +++ b/cpp/include/cm/signals.hpp @@ -94,6 +94,7 @@ struct SignalGridSpec { [[nodiscard]] std::size_t x_face_count() const; [[nodiscard]] std::size_t y_face_count() const; [[nodiscard]] std::size_t z_face_count() const; + void validate_lattice() const; void validate() const; }; @@ -108,9 +109,18 @@ struct SignalGridStencil { std::array sites{}; std::array weights{}; std::uint32_t count{0}; + // Every site of the stencil is solid, so it carries no fluid to interpolate + // and every weight is zero. Concentration there is undefined and sampling it + // is a model error; velocity there is zero by the field's own validation. + bool entirely_solid{false}; }; -[[nodiscard]] SignalGridStencil signal_grid_stencil(const SignalGridSpec& spec, Vec3 position); +// Whether a sample position outside the lattice of site centers is an error or +// is drawn in to the nearest in-lattice point. +enum class GridSampleBound : std::uint8_t { inside, clamped }; + +[[nodiscard]] SignalGridStencil signal_grid_stencil( + const SignalGridSpec& spec, Vec3 position, GridSampleBound bound = GridSampleBound::inside); class SignalGrid { public: @@ -121,10 +131,13 @@ class SignalGrid { [[nodiscard]] std::span levels() const& noexcept; [[nodiscard]] std::span levels() && = delete; [[nodiscard]] std::vector sample(Vec3 position) const; + [[nodiscard]] Vec3 sample_velocity(Vec3 position, + GridSampleBound bound = GridSampleBound::inside) const; [[nodiscard]] SignalGridCheckpoint checkpoint() const; void set_levels(std::span levels); void replace_levels(std::vector levels); void set_velocity_field(std::optional field); + void set_reaction(std::optional reaction); void validate_step(float dt) const; void validate() const; diff --git a/cpp/include/cm/simulation.hpp b/cpp/include/cm/simulation.hpp index 7b51898..00e8c55 100644 --- a/cpp/include/cm/simulation.hpp +++ b/cpp/include/cm/simulation.hpp @@ -31,6 +31,8 @@ class Simulation { CellId add_cell(const CellInit& cell); void remove_cell(CellId id); + void apply_flow_drift(float dt, const MechanicsIntegrationParameters& parameters = + MechanicsIntegrationParameters{}); ConstraintId add_plane_constraint(const PlaneConstraintInit& plane); ConstraintId add_sphere_constraint(const SphereConstraintInit& sphere); ConstraintId add_box_constraint(const BoxConstraintInit& box); @@ -45,6 +47,7 @@ class Simulation { void configure_signal_grid(const SignalGridSpec& spec, std::vector levels = {}); void set_signal_levels(std::span levels); void set_velocity_field(std::optional field); + void set_signal_reaction(std::optional reaction); std::pair divide(CellId parent_id, float first_fraction); std::pair divide_equal(CellId parent_id); void step(float dt); diff --git a/cpp/metal/kernels/contacts.metal b/cpp/metal/kernels/contacts.metal index 5883254..29a8c91 100644 --- a/cpp/metal/kernels/contacts.metal +++ b/cpp/metal/kernels/contacts.metal @@ -53,7 +53,7 @@ struct CenterlineMinimum { }; Capsule load_capsule(device const ulong* ids, device const float4* centers, - device const float4* axes, device const float4* geometry, uint slot) { + device const float4* axes, device const float4* geometry, uint slot) { Capsule result; result.id = ids[slot]; result.slot = slot; @@ -100,14 +100,12 @@ PointPair closest_points(const Capsule first, const Capsule second, float epsilo first_parameter = clamp(-first_projection / first_length_squared, 0.0f, 1.0f); } else { float cross_projection = dot(first_delta, second_delta); - float denominator = first_length_squared * second_length_squared - - cross_projection * cross_projection; - float parallel_tolerance = - float_epsilon * first_length_squared * second_length_squared; + float denominator = + first_length_squared * second_length_squared - cross_projection * cross_projection; + float parallel_tolerance = float_epsilon * first_length_squared * second_length_squared; if (denominator > parallel_tolerance) { first_parameter = clamp( - (cross_projection * second_projection - - first_projection * second_length_squared) / + (cross_projection * second_projection - first_projection * second_length_squared) / denominator, 0.0f, 1.0f); } @@ -428,15 +426,12 @@ ExternalEvaluation evaluate_external_constraint(const Capsule cell, return result; } -kernel void count_cell_contacts(device const ulong* ids [[buffer(0)]], - device const float4* centers [[buffer(1)]], - device const float4* axes [[buffer(2)]], - device const float4* geometry [[buffer(3)]], - device uint* counts [[buffer(4)]], - constant float4& parameters [[buffer(5)]], - constant uint& candidate_count [[buffer(6)]], - device const uint2* candidates [[buffer(7)]], - uint pair_index [[thread_position_in_grid]]) { +kernel void count_cell_contacts( + device const ulong* ids [[buffer(0)]], device const float4* centers [[buffer(1)]], + device const float4* axes [[buffer(2)]], device const float4* geometry [[buffer(3)]], + device uint* counts [[buffer(4)]], constant float4& parameters [[buffer(5)]], + constant uint& candidate_count [[buffer(6)]], device const uint2* candidates [[buffer(7)]], + uint pair_index [[thread_position_in_grid]]) { if (pair_index >= candidate_count) { return; } @@ -480,8 +475,7 @@ kernel void fill_cell_contacts( device uint* ordinals [[buffer(10)]], device float4* points_on_first [[buffer(11)]], device float4* normals [[buffer(12)]], device float* separations [[buffer(13)]], device float* weights [[buffer(14)]], constant float4& parameters [[buffer(15)]], - constant uint& candidate_count [[buffer(16)]], - device const uint2* candidates [[buffer(17)]], + constant uint& candidate_count [[buffer(16)]], device const uint2* candidates [[buffer(17)]], uint pair_index [[thread_position_in_grid]]) { if (pair_index >= candidate_count) { return; @@ -511,7 +505,8 @@ kernel void fill_cell_contacts( first_slots[output_index] = first.slot; second_slots[output_index] = second.slot; ordinals[output_index] = ordinal; - points_on_first[output_index] = float4(points.values[ordinal].first + normal * first.radius, 0.0f); + points_on_first[output_index] = + float4(points.values[ordinal].first + normal * first.radius, 0.0f); normals[output_index] = float4(normal, 0.0f); separations[output_index] = separation; weights[output_index] = weight; diff --git a/cpp/metal/kernels/coupled_rates.metal b/cpp/metal/kernels/coupled_rates.metal index a377601..83080b0 100644 --- a/cpp/metal/kernels/coupled_rates.metal +++ b/cpp/metal/kernels/coupled_rates.metal @@ -12,20 +12,7 @@ struct RateInstruction { float value; }; -struct GridShape { - uint x; - uint y; - uint z; - uint sites; -}; - -uint site_index(GridShape shape, uint x, uint y, uint z) { - return x * shape.y * shape.z + y * shape.z + z; -} - -float grid_level(device const float* levels, GridShape shape, uint signal, uint x, uint y, uint z) { - return levels[signal * shape.sites + site_index(shape, x, y, z)]; -} +; float effective_volume(float length, float radius) { return pi * radius * radius * (length + 2.0f * radius); @@ -104,6 +91,8 @@ float sample_signal(device const float* levels, GridShape shape, float4 origin, } } } + // A stencil with no fluid corner is rejected by the host's coupled-step + // validation before any kernel runs, so the fluid weight is positive here. if (dropped) { result /= fluid_weight; } @@ -112,10 +101,16 @@ float sample_signal(device const float* levels, GridShape shape, float4 origin, float cell_scatter_weight(float4 center, GridShape shape, float4 origin, float4 spacing, device const uchar* obstacles, uint x, uint y, uint z) { + // A cell only scatters into the eight sites of its own stencil, and the + // weight is pure arithmetic, so testing it first keeps the obstacle mask out + // of the sites a cell cannot reach - which is nearly all of them. + float raw = cell_site_weight(center, shape, origin, spacing, x, y, z); + if (raw == 0.0f) { + return 0.0f; + } if (obstacles[site_index(shape, x, y, z)] != 0u) { return 0.0f; } - float raw = cell_site_weight(center, shape, origin, spacing, x, y, z); float coordinate_x = axis_coordinate(center.x, origin.x, spacing.x, shape.x); float coordinate_y = axis_coordinate(center.y, origin.y, spacing.y, shape.y); float coordinate_z = axis_coordinate(center.z, origin.z, spacing.z, shape.z); @@ -147,6 +142,8 @@ float cell_scatter_weight(float4 center, GridShape shape, float4 origin, float4 } } } + // A stencil with no fluid corner is rejected by the host's coupled-step + // validation before any kernel runs, so the fluid weight is positive here. return dropped ? raw / fluid_weight : raw; } @@ -228,8 +225,8 @@ kernel void advance_coupled_cells( constant float4& origin [[buffer(14)]], constant float4& spacing [[buffer(15)]], constant float& dt [[buffer(16)]], constant uint& species_count [[buffer(17)]], constant uint& signal_count [[buffer(18)]], constant uint& instruction_count [[buffer(19)]], - constant uint& cell_count [[buffer(20)]], - device const uchar* obstacles [[buffer(21)]], uint cell [[thread_position_in_grid]]) { + constant uint& cell_count [[buffer(20)]], device const uchar* obstacles [[buffer(21)]], + uint cell [[thread_position_in_grid]]) { if (cell >= cell_count) { return; } @@ -275,17 +272,6 @@ kernel void advance_coupled_cells( } } -float exterior_value(uint kind, device const float* fixed_values, uint face, uint signal, - uint signal_count, float current, float periodic) { - if (kind == 0u) { - return current; - } - if (kind == 1u) { - return periodic; - } - return fixed_values[face * signal_count + signal]; -} - kernel void advance_coupled_grid( device const float* levels [[buffer(0)]], device float* output [[buffer(1)]], device const float* diffusion [[buffer(2)]], device const float4* advection [[buffer(3)]], @@ -297,11 +283,10 @@ kernel void advance_coupled_grid( constant uint& cell_count [[buffer(14)]], constant uint& level_count [[buffer(15)]], constant uint& crank_nicolson [[buffer(16)]], device const float* reaction_source [[buffer(17)]], - device const float* reaction_loss [[buffer(18)]], - device const uchar* obstacles [[buffer(19)]], + device const float* reaction_loss [[buffer(18)]], device const uchar* obstacles [[buffer(19)]], device const float* x_faces [[buffer(20)]], device const float* y_faces [[buffer(21)]], - device const float* z_faces [[buffer(22)]], - constant uint& has_velocity_field [[buffer(23)]], uint index [[thread_position_in_grid]]) { + device const float* z_faces [[buffer(22)]], constant uint& has_velocity_field [[buffer(23)]], + uint index [[thread_position_in_grid]]) { if (index >= level_count) { return; } @@ -343,51 +328,10 @@ kernel void advance_coupled_grid( : grid_level(levels, shape, signal, x, y, z + 1u); uint3 dimensions = uint3(shape.x, shape.y, shape.z); - bool3 closed_lower; - bool3 closed_upper; - closed_lower.x = - x == 0u ? (boundary_kinds[0] == 0u || - (boundary_kinds[0] == 1u && obstacles[site_index(shape, shape.x - 1u, y, z)] != 0u)) - : obstacles[site_index(shape, x - 1u, y, z)] != 0u; - closed_upper.x = - x + 1u == shape.x - ? (boundary_kinds[1] == 0u || - (boundary_kinds[1] == 1u && obstacles[site_index(shape, 0u, y, z)] != 0u)) - : obstacles[site_index(shape, x + 1u, y, z)] != 0u; - closed_lower.y = - y == 0u ? (boundary_kinds[2] == 0u || - (boundary_kinds[2] == 1u && obstacles[site_index(shape, x, shape.y - 1u, z)] != 0u)) - : obstacles[site_index(shape, x, y - 1u, z)] != 0u; - closed_upper.y = - y + 1u == shape.y - ? (boundary_kinds[3] == 0u || - (boundary_kinds[3] == 1u && obstacles[site_index(shape, x, 0u, z)] != 0u)) - : obstacles[site_index(shape, x, y + 1u, z)] != 0u; - closed_lower.z = - z == 0u ? (boundary_kinds[4] == 0u || - (boundary_kinds[4] == 1u && obstacles[site_index(shape, x, y, shape.z - 1u)] != 0u)) - : obstacles[site_index(shape, x, y, z - 1u)] != 0u; - closed_upper.z = - z + 1u == shape.z - ? (boundary_kinds[5] == 0u || - (boundary_kinds[5] == 1u && obstacles[site_index(shape, x, y, 0u)] != 0u)) - : obstacles[site_index(shape, x, y, z + 1u)] != 0u; - float face_lower[3]; - float face_upper[3]; - if (has_velocity_field != 0u) { - face_lower[0] = x_faces[x * shape.y * shape.z + y * shape.z + z]; - face_upper[0] = x_faces[(x + 1u) * shape.y * shape.z + y * shape.z + z]; - face_lower[1] = y_faces[x * (shape.y + 1u) * shape.z + y * shape.z + z]; - face_upper[1] = y_faces[x * (shape.y + 1u) * shape.z + (y + 1u) * shape.z + z]; - face_lower[2] = z_faces[x * shape.y * (shape.z + 1u) + y * (shape.z + 1u) + z]; - face_upper[2] = z_faces[x * shape.y * (shape.z + 1u) + y * (shape.z + 1u) + z + 1u]; - } else { - float3 velocity = advection[signal].xyz; - for (uint axis = 0; axis < 3u; ++axis) { - face_lower[axis] = velocity[axis]; - face_upper[axis] = velocity[axis]; - } - } + GridFaceState faces = grid_face_state(shape, boundary_kinds, obstacles, x_faces, y_faces, z_faces, + has_velocity_field, advection[signal], x, y, z); + bool3 closed_lower = faces.closed_lower; + bool3 closed_upper = faces.closed_upper; float3 grid_spacing = spacing.xyz; float rate = 0.0f; for (uint axis = 0; axis < 3u; ++axis) { @@ -403,10 +347,10 @@ kernel void advance_coupled_grid( float inverse_spacing = 1.0f / grid_spacing[axis]; rate += diffusion[signal] * (lower[axis] - 2.0f * current + upper[axis]) * inverse_spacing * inverse_spacing; - float lower_flux = face_lower[axis] >= 0.0f ? face_lower[axis] * lower[axis] - : face_lower[axis] * current; - float upper_flux = face_upper[axis] >= 0.0f ? face_upper[axis] * current - : face_upper[axis] * upper[axis]; + float lower_flux = + faces.lower[axis] >= 0.0f ? faces.lower[axis] * lower[axis] : faces.lower[axis] * current; + float upper_flux = + faces.upper[axis] >= 0.0f ? faces.upper[axis] * current : faces.upper[axis] * upper[axis]; if (closed_lower[axis]) { lower_flux = 0.0f; } diff --git a/cpp/metal/kernels/grid_transport.metal b/cpp/metal/kernels/grid_transport.metal new file mode 100644 index 0000000..7e3d744 --- /dev/null +++ b/cpp/metal/kernels/grid_transport.metal @@ -0,0 +1,96 @@ +// Grid geometry and transport helpers shared by the signal and coupled-rate +// kernels. A Metal library is compiled from source at runtime, with no include +// path, so this fragment is concatenated ahead of each kernel source at build +// time rather than included by it. + +#include + +using namespace metal; + +struct GridShape { + uint x; + uint y; + uint z; + uint sites; +}; + +uint site_index(GridShape shape, uint x, uint y, uint z) { + return x * shape.y * shape.z + y * shape.z + z; +} + +float grid_level(device const float* levels, GridShape shape, uint signal, uint x, uint y, uint z) { + return levels[signal * shape.sites + site_index(shape, x, y, z)]; +} + +float exterior_value(uint kind, device const float* fixed_values, uint face, uint signal, + uint signal_count, float current, float periodic) { + if (kind == 0u) { + return current; + } + if (kind == 1u) { + return periodic; + } + return fixed_values[face * signal_count + signal]; +} + +// Whether each of a site's six faces is closed to transport, and the velocity +// it carries. A face is closed by a no-flux boundary, by a periodic boundary +// that wraps onto a solid site, or by a solid neighbour. +struct GridFaceState { + bool3 closed_lower; + bool3 closed_upper; + float lower[3]; + float upper[3]; +}; + +GridFaceState grid_face_state(GridShape shape, constant uint* boundary_kinds, + device const uchar* obstacles, device const float* x_faces, + device const float* y_faces, device const float* z_faces, + uint has_velocity_field, float4 advection, uint x, uint y, uint z) { + GridFaceState faces; + faces.closed_lower.x = + x == 0u + ? (boundary_kinds[0] == 0u || + (boundary_kinds[0] == 1u && obstacles[site_index(shape, shape.x - 1u, y, z)] != 0u)) + : obstacles[site_index(shape, x - 1u, y, z)] != 0u; + faces.closed_upper.x = + x + 1u == shape.x + ? (boundary_kinds[1] == 0u || + (boundary_kinds[1] == 1u && obstacles[site_index(shape, 0u, y, z)] != 0u)) + : obstacles[site_index(shape, x + 1u, y, z)] != 0u; + faces.closed_lower.y = + y == 0u + ? (boundary_kinds[2] == 0u || + (boundary_kinds[2] == 1u && obstacles[site_index(shape, x, shape.y - 1u, z)] != 0u)) + : obstacles[site_index(shape, x, y - 1u, z)] != 0u; + faces.closed_upper.y = + y + 1u == shape.y + ? (boundary_kinds[3] == 0u || + (boundary_kinds[3] == 1u && obstacles[site_index(shape, x, 0u, z)] != 0u)) + : obstacles[site_index(shape, x, y + 1u, z)] != 0u; + faces.closed_lower.z = + z == 0u + ? (boundary_kinds[4] == 0u || + (boundary_kinds[4] == 1u && obstacles[site_index(shape, x, y, shape.z - 1u)] != 0u)) + : obstacles[site_index(shape, x, y, z - 1u)] != 0u; + faces.closed_upper.z = + z + 1u == shape.z + ? (boundary_kinds[5] == 0u || + (boundary_kinds[5] == 1u && obstacles[site_index(shape, x, y, 0u)] != 0u)) + : obstacles[site_index(shape, x, y, z + 1u)] != 0u; + if (has_velocity_field != 0u) { + faces.lower[0] = x_faces[x * shape.y * shape.z + y * shape.z + z]; + faces.upper[0] = x_faces[(x + 1u) * shape.y * shape.z + y * shape.z + z]; + faces.lower[1] = y_faces[x * (shape.y + 1u) * shape.z + y * shape.z + z]; + faces.upper[1] = y_faces[x * (shape.y + 1u) * shape.z + (y + 1u) * shape.z + z]; + faces.lower[2] = z_faces[x * shape.y * (shape.z + 1u) + y * (shape.z + 1u) + z]; + faces.upper[2] = z_faces[x * shape.y * (shape.z + 1u) + y * (shape.z + 1u) + z + 1u]; + } else { + float3 velocity = advection.xyz; + for (uint axis = 0; axis < 3u; ++axis) { + faces.lower[axis] = velocity[axis]; + faces.upper[axis] = velocity[axis]; + } + } + return faces; +} diff --git a/cpp/metal/kernels/mechanics.metal b/cpp/metal/kernels/mechanics.metal index c54dbcc..b88fbe0 100644 --- a/cpp/metal/kernels/mechanics.metal +++ b/cpp/metal/kernels/mechanics.metal @@ -36,8 +36,7 @@ MechanicsDofs contact_jacobian(float3 normal, float3 arm, float3 axis, float tot float weight) { MechanicsDofs result; result.linear_length = - float4(weight * normal, - weight * dot(axis, arm) * dot(axis, normal) / total_length); + float4(weight * normal, weight * dot(axis, arm) * dot(axis, normal) / total_length); result.rotation = float4(weight * cross(arm, normal), 0.0f); return result; } @@ -48,9 +47,8 @@ kernel void build_mechanics_rows( device const uint* second_slots [[buffer(4)]], device const float4* points [[buffer(5)]], device const float4* normals [[buffer(6)]], device const float* separations [[buffer(7)]], device const float* weights [[buffer(8)]], device MechanicsDofs* first_rows [[buffer(9)]], - device MechanicsDofs* second_rows [[buffer(10)]], - device float* right_hand_side [[buffer(11)]], constant uint& contact_count [[buffer(12)]], - uint index [[thread_position_in_grid]]) { + device MechanicsDofs* second_rows [[buffer(10)]], device float* right_hand_side [[buffer(11)]], + constant uint& contact_count [[buffer(12)]], uint index [[thread_position_in_grid]]) { if (index >= contact_count) { return; } @@ -59,9 +57,8 @@ kernel void build_mechanics_rows( float weight = weights[index]; float3 normal = normals[index].xyz; float3 point = points[index].xyz; - first_rows[index] = - contact_jacobian(normal, point - centers[first].xyz, axes[first].xyz, - geometry[first].x + 2.0f * geometry[first].y, weight); + first_rows[index] = contact_jacobian(normal, point - centers[first].xyz, axes[first].xyz, + geometry[first].x + 2.0f * geometry[first].y, weight); second_rows[index] = second == 0xffffffffu ? zero_dofs() @@ -70,14 +67,15 @@ kernel void build_mechanics_rows( right_hand_side[index] = weight * separations[index]; } -kernel void apply_mechanics_b( - device const MechanicsDofs* first_rows [[buffer(0)]], - device const MechanicsDofs* second_rows [[buffer(1)]], - device const uint* first_slots [[buffer(2)]], - device const uint* second_slots [[buffer(3)]], - device const MechanicsDofs* input [[buffer(4)]], device float* row_values [[buffer(5)]], - constant uint& contact_count [[buffer(6)]], device const uchar* fixed [[buffer(7)]], - uint index [[thread_position_in_grid]]) { +kernel void apply_mechanics_b(device const MechanicsDofs* first_rows [[buffer(0)]], + device const MechanicsDofs* second_rows [[buffer(1)]], + device const uint* first_slots [[buffer(2)]], + device const uint* second_slots [[buffer(3)]], + device const MechanicsDofs* input [[buffer(4)]], + device float* row_values [[buffer(5)]], + constant uint& contact_count [[buffer(6)]], + device const uchar* fixed [[buffer(7)]], + uint index [[thread_position_in_grid]]) { if (index >= contact_count) { return; } @@ -91,14 +89,15 @@ kernel void apply_mechanics_b( } } -kernel void apply_mechanics_transpose( - device const MechanicsDofs* first_rows [[buffer(0)]], - device const MechanicsDofs* second_rows [[buffer(1)]], - device const float* row_values [[buffer(2)]], - device const uint* incidence_offsets [[buffer(3)]], - device const uint* incidence_indices [[buffer(4)]], - device const uint* first_slots [[buffer(5)]], device MechanicsDofs* output [[buffer(6)]], - constant uint& cell_count [[buffer(7)]], uint cell [[thread_position_in_grid]]) { +kernel void apply_mechanics_transpose(device const MechanicsDofs* first_rows [[buffer(0)]], + device const MechanicsDofs* second_rows [[buffer(1)]], + device const float* row_values [[buffer(2)]], + device const uint* incidence_offsets [[buffer(3)]], + device const uint* incidence_indices [[buffer(4)]], + device const uint* first_slots [[buffer(5)]], + device MechanicsDofs* output [[buffer(6)]], + constant uint& cell_count [[buffer(7)]], + uint cell [[thread_position_in_grid]]) { if (cell >= cell_count) { return; } @@ -131,8 +130,7 @@ kernel void add_mechanics_regularizer( float radius = geometry[cell].y; float mass = mu_a * total_length; float axial_inertia = 0.5f * mass * radius * radius; - float transverse_inertia = - mass * (total_length * total_length + 3.0f * radius * radius) / 12.0f; + float transverse_inertia = mass * (total_length * total_length + 3.0f * radius * radius) / 12.0f; float3 axis = axes[cell].xyz; float3 rotation = input[cell].rotation.xyz; float3 inertia_rotation = rotation * transverse_inertia + @@ -146,12 +144,13 @@ kernel void add_mechanics_regularizer( output[cell] = result; } -kernel void initialize_mechanics_vectors( - device MechanicsDofs* right_hand_side [[buffer(0)]], - device MechanicsDofs* solution [[buffer(1)]], device MechanicsDofs* residual [[buffer(2)]], - device MechanicsDofs* search_direction [[buffer(3)]], - constant uint& cell_count [[buffer(4)]], device const uchar* fixed [[buffer(5)]], - uint cell [[thread_position_in_grid]]) { +kernel void initialize_mechanics_vectors(device MechanicsDofs* right_hand_side [[buffer(0)]], + device MechanicsDofs* solution [[buffer(1)]], + device MechanicsDofs* residual [[buffer(2)]], + device MechanicsDofs* search_direction [[buffer(3)]], + constant uint& cell_count [[buffer(4)]], + device const uchar* fixed [[buffer(5)]], + uint cell [[thread_position_in_grid]]) { if (cell >= cell_count) { return; } @@ -174,20 +173,22 @@ kernel void update_mechanics_solution_residual( residual[cell] = added(residual[cell], scaled(applied[cell], -alpha)); } -kernel void update_mechanics_search_direction( - device const MechanicsDofs* residual [[buffer(0)]], - device MechanicsDofs* search_direction [[buffer(1)]], constant float& beta [[buffer(2)]], - constant uint& cell_count [[buffer(3)]], uint cell [[thread_position_in_grid]]) { +kernel void update_mechanics_search_direction(device const MechanicsDofs* residual [[buffer(0)]], + device MechanicsDofs* search_direction [[buffer(1)]], + constant float& beta [[buffer(2)]], + constant uint& cell_count [[buffer(3)]], + uint cell [[thread_position_in_grid]]) { if (cell >= cell_count) { return; } search_direction[cell] = added(residual[cell], scaled(search_direction[cell], beta)); } -kernel void subtract_mechanics_vectors( - device const MechanicsDofs* left [[buffer(0)]], - device const MechanicsDofs* right [[buffer(1)]], device MechanicsDofs* output [[buffer(2)]], - constant uint& cell_count [[buffer(3)]], uint cell [[thread_position_in_grid]]) { +kernel void subtract_mechanics_vectors(device const MechanicsDofs* left [[buffer(0)]], + device const MechanicsDofs* right [[buffer(1)]], + device MechanicsDofs* output [[buffer(2)]], + constant uint& cell_count [[buffer(3)]], + uint cell [[thread_position_in_grid]]) { if (cell >= cell_count) { return; } diff --git a/cpp/metal/kernels/signals.metal b/cpp/metal/kernels/signals.metal index eddf3b7..8e3e8d9 100644 --- a/cpp/metal/kernels/signals.metal +++ b/cpp/metal/kernels/signals.metal @@ -2,31 +2,7 @@ using namespace metal; -struct GridShape { - uint x; - uint y; - uint z; - uint sites; -}; - -uint site_index(GridShape shape, uint x, uint y, uint z) { - return x * shape.y * shape.z + y * shape.z + z; -} - -float grid_level(device const float* levels, GridShape shape, uint signal, uint x, uint y, uint z) { - return levels[signal * shape.sites + site_index(shape, x, y, z)]; -} - -float exterior_value(uint kind, device const float* fixed_values, uint face, uint signal, - uint signal_count, float current, float periodic) { - if (kind == 0u) { - return current; - } - if (kind == 1u) { - return periodic; - } - return fixed_values[face * signal_count + signal]; -} +; struct TransportPoint { float rate; @@ -36,11 +12,11 @@ struct TransportPoint { TransportPoint transport_point(device const float* levels, device const float* diffusion, device const float4* advection, device const float* fixed_values, device const float* reaction_source, - device const float* reaction_loss, - device const uchar* obstacles, device const float* x_faces, - device const float* y_faces, device const float* z_faces, - uint has_velocity_field, constant uint* boundary_kinds, - GridShape shape, float4 spacing, uint signal_count, uint index) { + device const float* reaction_loss, device const uchar* obstacles, + device const float* x_faces, device const float* y_faces, + device const float* z_faces, uint has_velocity_field, + constant uint* boundary_kinds, GridShape shape, float4 spacing, + uint signal_count, uint index) { uint signal = index / shape.sites; uint site = index - signal * shape.sites; uint x = site / (shape.y * shape.z); @@ -77,51 +53,10 @@ TransportPoint transport_point(device const float* levels, device const float* d : grid_level(levels, shape, signal, x, y, z + 1u); uint3 dimensions = uint3(shape.x, shape.y, shape.z); - bool3 closed_lower; - bool3 closed_upper; - closed_lower.x = - x == 0u ? (boundary_kinds[0] == 0u || - (boundary_kinds[0] == 1u && obstacles[site_index(shape, shape.x - 1u, y, z)] != 0u)) - : obstacles[site_index(shape, x - 1u, y, z)] != 0u; - closed_upper.x = - x + 1u == shape.x - ? (boundary_kinds[1] == 0u || - (boundary_kinds[1] == 1u && obstacles[site_index(shape, 0u, y, z)] != 0u)) - : obstacles[site_index(shape, x + 1u, y, z)] != 0u; - closed_lower.y = - y == 0u ? (boundary_kinds[2] == 0u || - (boundary_kinds[2] == 1u && obstacles[site_index(shape, x, shape.y - 1u, z)] != 0u)) - : obstacles[site_index(shape, x, y - 1u, z)] != 0u; - closed_upper.y = - y + 1u == shape.y - ? (boundary_kinds[3] == 0u || - (boundary_kinds[3] == 1u && obstacles[site_index(shape, x, 0u, z)] != 0u)) - : obstacles[site_index(shape, x, y + 1u, z)] != 0u; - closed_lower.z = - z == 0u ? (boundary_kinds[4] == 0u || - (boundary_kinds[4] == 1u && obstacles[site_index(shape, x, y, shape.z - 1u)] != 0u)) - : obstacles[site_index(shape, x, y, z - 1u)] != 0u; - closed_upper.z = - z + 1u == shape.z - ? (boundary_kinds[5] == 0u || - (boundary_kinds[5] == 1u && obstacles[site_index(shape, x, y, 0u)] != 0u)) - : obstacles[site_index(shape, x, y, z + 1u)] != 0u; - float face_lower[3]; - float face_upper[3]; - if (has_velocity_field != 0u) { - face_lower[0] = x_faces[x * shape.y * shape.z + y * shape.z + z]; - face_upper[0] = x_faces[(x + 1u) * shape.y * shape.z + y * shape.z + z]; - face_lower[1] = y_faces[x * (shape.y + 1u) * shape.z + y * shape.z + z]; - face_upper[1] = y_faces[x * (shape.y + 1u) * shape.z + (y + 1u) * shape.z + z]; - face_lower[2] = z_faces[x * shape.y * (shape.z + 1u) + y * (shape.z + 1u) + z]; - face_upper[2] = z_faces[x * shape.y * (shape.z + 1u) + y * (shape.z + 1u) + z + 1u]; - } else { - float3 velocity = advection[signal].xyz; - for (uint axis = 0; axis < 3u; ++axis) { - face_lower[axis] = velocity[axis]; - face_upper[axis] = velocity[axis]; - } - } + GridFaceState faces = grid_face_state(shape, boundary_kinds, obstacles, x_faces, y_faces, z_faces, + has_velocity_field, advection[signal], x, y, z); + bool3 closed_lower = faces.closed_lower; + bool3 closed_upper = faces.closed_upper; float3 grid_spacing = spacing.xyz; float rate = 0.0f; float diagonal = 0.0f; @@ -145,10 +80,10 @@ TransportPoint transport_point(device const float* levels, device const float* d if (closed_upper[axis]) { diagonal += diffusion_scale; } - float lower_flux = face_lower[axis] >= 0.0f ? face_lower[axis] * lower[axis] - : face_lower[axis] * current; - float upper_flux = face_upper[axis] >= 0.0f ? face_upper[axis] * current - : face_upper[axis] * upper[axis]; + float lower_flux = + faces.lower[axis] >= 0.0f ? faces.lower[axis] * lower[axis] : faces.lower[axis] * current; + float upper_flux = + faces.upper[axis] >= 0.0f ? faces.upper[axis] * current : faces.upper[axis] * upper[axis]; if (closed_lower[axis]) { lower_flux = 0.0f; } @@ -156,11 +91,11 @@ TransportPoint transport_point(device const float* levels, device const float* d upper_flux = 0.0f; } rate -= (upper_flux - lower_flux) * inverse_spacing; - if (!closed_upper[axis] && face_upper[axis] > 0.0f) { - diagonal -= face_upper[axis] * inverse_spacing; + if (!closed_upper[axis] && faces.upper[axis] > 0.0f) { + diagonal -= faces.upper[axis] * inverse_spacing; } - if (!closed_lower[axis] && face_lower[axis] < 0.0f) { - diagonal += face_lower[axis] * inverse_spacing; + if (!closed_lower[axis] && faces.lower[axis] < 0.0f) { + diagonal += faces.lower[axis] * inverse_spacing; } } rate += reaction_source[index] - reaction_loss[index] * current; @@ -177,11 +112,10 @@ kernel void advance_signal_grid( constant uint& signal_count [[buffer(10)]], constant uint& level_count [[buffer(11)]], constant uint& crank_nicolson [[buffer(12)]], device const float* reaction_source [[buffer(13)]], - device const float* reaction_loss [[buffer(14)]], - device const uchar* obstacles [[buffer(15)]], + device const float* reaction_loss [[buffer(14)]], device const uchar* obstacles [[buffer(15)]], device const float* x_faces [[buffer(16)]], device const float* y_faces [[buffer(17)]], - device const float* z_faces [[buffer(18)]], - constant uint& has_velocity_field [[buffer(19)]], uint index [[thread_position_in_grid]]) { + device const float* z_faces [[buffer(18)]], constant uint& has_velocity_field [[buffer(19)]], + uint index [[thread_position_in_grid]]) { if (index >= level_count) { return; } @@ -207,11 +141,10 @@ kernel void crank_nicolson_jacobi( constant GridShape& shape [[buffer(8)]], constant float4& spacing [[buffer(9)]], constant float& half_dt [[buffer(10)]], constant uint& signal_count [[buffer(11)]], constant uint& level_count [[buffer(12)]], device const float* reaction_source [[buffer(13)]], - device const float* reaction_loss [[buffer(14)]], - device const uchar* obstacles [[buffer(15)]], + device const float* reaction_loss [[buffer(14)]], device const uchar* obstacles [[buffer(15)]], device const float* x_faces [[buffer(16)]], device const float* y_faces [[buffer(17)]], - device const float* z_faces [[buffer(18)]], - constant uint& has_velocity_field [[buffer(19)]], uint index [[thread_position_in_grid]]) { + device const float* z_faces [[buffer(18)]], constant uint& has_velocity_field [[buffer(19)]], + uint index [[thread_position_in_grid]]) { if (index >= level_count) { return; } @@ -236,11 +169,10 @@ kernel void crank_nicolson_residual_terms( constant float4& spacing [[buffer(8)]], constant float& half_dt [[buffer(9)]], constant uint& signal_count [[buffer(10)]], constant uint& level_count [[buffer(11)]], device const float* reaction_source [[buffer(12)]], - device const float* reaction_loss [[buffer(13)]], - device const uchar* obstacles [[buffer(14)]], + device const float* reaction_loss [[buffer(13)]], device const uchar* obstacles [[buffer(14)]], device const float* x_faces [[buffer(15)]], device const float* y_faces [[buffer(16)]], - device const float* z_faces [[buffer(17)]], - constant uint& has_velocity_field [[buffer(18)]], uint index [[thread_position_in_grid]]) { + device const float* z_faces [[buffer(17)]], constant uint& has_velocity_field [[buffer(18)]], + uint index [[thread_position_in_grid]]) { if (index >= level_count) { return; } diff --git a/cpp/metal/kernels/species.metal b/cpp/metal/kernels/species.metal index 15b84ba..6474b4e 100644 --- a/cpp/metal/kernels/species.metal +++ b/cpp/metal/kernels/species.metal @@ -20,13 +20,9 @@ float effective_surface_area(float length, float radius) { return 2.0f * pi * radius * (length + 2.0f * radius); } -float evaluate_instruction(const RateInstruction instruction, - device const float* workspace, - device const float* species, - float4 center, - float4 geometry, - float growth_rate, - int cell_type) { +float evaluate_instruction(const RateInstruction instruction, device const float* workspace, + device const float* species, float4 center, float4 geometry, + float growth_rate, int cell_type) { switch (instruction.operation) { case 0: return instruction.value; @@ -82,36 +78,29 @@ float evaluate_instruction(const RateInstruction instruction, return workspace[instruction.first] == workspace[instruction.second] ? 1.0f : 0.0f; case 26: return workspace[instruction.first] != 0.0f ? workspace[instruction.second] - : workspace[instruction.third]; + : workspace[instruction.third]; default: return NAN; } } kernel void advance_species( - device float* levels [[buffer(0)]], - device const float* previous_lengths [[buffer(1)]], - device const float4* centers [[buffer(2)]], - device const float4* geometry [[buffer(3)]], - device const float* growth_rates [[buffer(4)]], - device const int* cell_types [[buffer(5)]], + device float* levels [[buffer(0)]], device const float* previous_lengths [[buffer(1)]], + device const float4* centers [[buffer(2)]], device const float4* geometry [[buffer(3)]], + device const float* growth_rates [[buffer(4)]], device const int* cell_types [[buffer(5)]], device const RateInstruction* instructions [[buffer(6)]], - device const uint* outputs [[buffer(7)]], - device float* workspace [[buffer(8)]], - device atomic_uint* error [[buffer(9)]], - constant float& dt [[buffer(10)]], - constant uint& species_count [[buffer(11)]], - constant uint& instruction_count [[buffer(12)]], - constant uint& cell_count [[buffer(13)]], - uint cell [[thread_position_in_grid]]) { + device const uint* outputs [[buffer(7)]], device float* workspace [[buffer(8)]], + device atomic_uint* error [[buffer(9)]], constant float& dt [[buffer(10)]], + constant uint& species_count [[buffer(11)]], constant uint& instruction_count [[buffer(12)]], + constant uint& cell_count [[buffer(13)]], uint cell [[thread_position_in_grid]]) { if (cell >= cell_count) { return; } uint species_offset = cell * species_count; float radius = geometry[cell].y; - float dilution = effective_volume(previous_lengths[cell], radius) / - effective_volume(geometry[cell].x, radius); + float dilution = + effective_volume(previous_lengths[cell], radius) / effective_volume(geometry[cell].x, radius); for (uint species = 0; species < species_count; ++species) { levels[species_offset + species] *= dilution; } @@ -120,9 +109,9 @@ kernel void advance_species( device float* cell_workspace = workspace + workspace_offset; device const float* cell_species = levels + species_offset; for (uint index = 0; index < instruction_count; ++index) { - float value = evaluate_instruction(instructions[index], cell_workspace, cell_species, - centers[cell], geometry[cell], growth_rates[cell], - cell_types[cell]); + float value = + evaluate_instruction(instructions[index], cell_workspace, cell_species, centers[cell], + geometry[cell], growth_rates[cell], cell_types[cell]); cell_workspace[index] = value; if (!isfinite(value)) { atomic_fetch_or_explicit(error, 1u, memory_order_relaxed); diff --git a/cpp/metal/metal_backend.mm b/cpp/metal/metal_backend.mm index e25af2f..95464d7 100644 --- a/cpp/metal/metal_backend.mm +++ b/cpp/metal/metal_backend.mm @@ -461,8 +461,7 @@ SignalSolveReport advance_signal_grid(SignalGrid& grid, float dt) override { } ensure_signal_face_capacity(largest_face_count(spec)); fill_velocity_faces(spec, signal_x_faces_, signal_y_faces_, signal_z_faces_); - const auto has_velocity_field = - static_cast(spec.velocity_field.has_value()); + const auto has_velocity_field = static_cast(spec.velocity_field.has_value()); auto* advection = static_cast(signal_advection_.contents); for (std::size_t signal = 0; signal < signal_count; ++signal) { advection[signal] = { @@ -643,8 +642,7 @@ SignalSolveReport advance_coupled(WorldState& state, SignalGrid& grid, } ensure_coupled_face_capacity(largest_face_count(spec)); fill_velocity_faces(spec, coupled_x_faces_, coupled_y_faces_, coupled_z_faces_); - const auto has_velocity_field = - static_cast(spec.velocity_field.has_value()); + const auto has_velocity_field = static_cast(spec.velocity_field.has_value()); auto* advection = static_cast(coupled_advection_.contents); for (std::size_t signal = 0; signal < signal_count_size; ++signal) { advection[signal] = {spec.advection[signal].x, spec.advection[signal].y, @@ -1048,8 +1046,8 @@ void ensure_signal_capacity(std::size_t level_count, std::size_t signal_count) { allocate_shared_buffer(device_, byte_count, "signal-grid affine sources"); signal_reaction_loss_ = allocate_shared_buffer(device_, byte_count, "signal-grid affine losses"); - signal_obstacles_ = allocate_shared_buffer(device_, signal_level_capacity_, - "signal-grid obstacles"); + signal_obstacles_ = + allocate_shared_buffer(device_, signal_level_capacity_, "signal-grid obstacles"); } if (signal_count > signal_count_capacity_) { signal_count_capacity_ = std::bit_ceil(signal_count); @@ -1119,38 +1117,14 @@ void ensure_signal_solve_capacity(std::uint32_t level_count) { return input; } - [[nodiscard]] float signal_rhs_rms(id right_hand_side, std::uint32_t level_count) { - @autoreleasepool { - id command_buffer = [queue_ commandBuffer]; - id encoder = [command_buffer computeCommandEncoder]; - if (command_buffer == nil || encoder == nil) { - throw std::runtime_error("failed to create a Metal signal norm command"); - } - [encoder setComputePipelineState:signals_square_pipeline_]; - [encoder setBuffer:right_hand_side offset:0 atIndex:0]; - [encoder setBuffer:signal_cn_terms_ offset:0 atIndex:1]; - [encoder setBytes:&level_count length:sizeof(level_count) atIndex:2]; - dispatch_1d(encoder, signals_square_pipeline_, level_count); - [encoder memoryBarrierWithScope:MTLBarrierScopeBuffers]; - const auto reduction = encode_signal_reduction(encoder, level_count); - [encoder endEncoding]; - wait_for_command(command_buffer, "Metal signal norm failed"); - const auto sum = *static_cast(reduction.contents); - return std::sqrt(sum / static_cast(level_count)); - } - } - - id encode_signal_residual(id encoder, id current, - id right_hand_side, id diffusion, - id advection, id fixed_values, - id reaction_source, id reaction_loss, - id obstacles, id x_faces, - id y_faces, id z_faces, - std::uint32_t has_velocity_field, - const std::array& boundary_kinds, - const MetalUInt4& shape, const MetalFloat4& spacing, - float half_dt, std::uint32_t signal_count, - std::uint32_t level_count) { + id encode_signal_residual( + id encoder, id current, id right_hand_side, + id diffusion, id advection, id fixed_values, + id reaction_source, id reaction_loss, id obstacles, + id x_faces, id y_faces, id z_faces, + std::uint32_t has_velocity_field, const std::array& boundary_kinds, + const MetalUInt4& shape, const MetalFloat4& spacing, float half_dt, + std::uint32_t signal_count, std::uint32_t level_count) { [encoder setComputePipelineState:signals_cn_residual_pipeline_]; [encoder setBuffer:current offset:0 atIndex:0]; [encoder setBuffer:right_hand_side offset:0 atIndex:1]; @@ -1183,8 +1157,7 @@ void ensure_signal_solve_capacity(std::uint32_t level_count) { id fixed_values, id reaction_source, id reaction_loss, id obstacles, id x_faces, id y_faces, - id z_faces, - std::uint32_t has_velocity_field, + id z_faces, std::uint32_t has_velocity_field, const std::array& boundary_kinds, const MetalUInt4& shape, const MetalFloat4& spacing, float half_dt, std::uint32_t signal_count, @@ -1206,25 +1179,54 @@ void ensure_signal_solve_capacity(std::uint32_t level_count) { } } + [[nodiscard]] float signal_rhs_rms(id right_hand_side, std::uint32_t level_count) { + @autoreleasepool { + id command_buffer = [queue_ commandBuffer]; + id encoder = [command_buffer computeCommandEncoder]; + if (command_buffer == nil || encoder == nil) { + throw std::runtime_error("failed to create a Metal signal norm command"); + } + [encoder setComputePipelineState:signals_square_pipeline_]; + [encoder setBuffer:right_hand_side offset:0 atIndex:0]; + [encoder setBuffer:signal_cn_terms_ offset:0 atIndex:1]; + [encoder setBytes:&level_count length:sizeof(level_count) atIndex:2]; + dispatch_1d(encoder, signals_square_pipeline_, level_count); + [encoder memoryBarrierWithScope:MTLBarrierScopeBuffers]; + const auto reduction = encode_signal_reduction(encoder, level_count); + [encoder endEncoding]; + wait_for_command(command_buffer, "Metal signal norm failed"); + const auto sum = *static_cast(reduction.contents); + return std::sqrt(sum / static_cast(level_count)); + } + } + [[nodiscard]] std::pair, SignalSolveReport> solve_signal_crank_nicolson( id initial, id right_hand_side, id diffusion, id advection, id fixed_values, id reaction_source, id reaction_loss, id obstacles, id x_faces, id y_faces, id z_faces, std::uint32_t has_velocity_field, - id error, - const std::array& boundary_kinds, const MetalUInt4& shape, - const MetalFloat4& spacing, float dt, std::uint32_t signal_count, std::uint32_t level_count, - const SignalSolveParameters& parameters) { + id error, const std::array& boundary_kinds, + const MetalUInt4& shape, const MetalFloat4& spacing, float dt, std::uint32_t signal_count, + std::uint32_t level_count, const SignalSolveParameters& parameters) { ensure_signal_solve_capacity(level_count); const auto half_dt = 0.5F * dt; - const auto right_hand_side_rms = signal_rhs_rms(right_hand_side, level_count); - const auto threshold = - parameters.absolute_tolerance + (parameters.relative_tolerance * right_hand_side_rms); SignalSolveReport report; report.residual_rms = signal_residual_rms( initial, right_hand_side, diffusion, advection, fixed_values, reaction_source, reaction_loss, obstacles, x_faces, y_faces, z_faces, has_velocity_field, boundary_kinds, shape, spacing, half_dt, signal_count, level_count); + // The relative term scales the residual the step starts with, matching the + // CPU reference: the field's own magnitude says nothing about how much of + // it this step has to change. + // A residual cannot fall below what float32 can represent for a field of + // this magnitude, and the right-hand side carries both the field and the + // operator terms of the step, so it sets that floor. An absolute tolerance + // asking for less than the floor is raised to it rather than making the + // solve unreachable. + const auto floor = + std::numeric_limits::epsilon() * signal_rhs_rms(right_hand_side, level_count); + const auto threshold = std::max(parameters.absolute_tolerance, floor) + + (parameters.relative_tolerance * report.residual_rms); if (std::isfinite(report.residual_rms) && report.residual_rms <= threshold) { return {initial, report}; } @@ -1270,8 +1272,8 @@ void ensure_signal_solve_capacity(std::uint32_t level_count) { [encoder memoryBarrierWithScope:MTLBarrierScopeBuffers]; const auto reduction = encode_signal_residual( encoder, output, right_hand_side, diffusion, advection, fixed_values, reaction_source, - reaction_loss, obstacles, x_faces, y_faces, z_faces, has_velocity_field, - boundary_kinds, shape, spacing, half_dt, signal_count, level_count); + reaction_loss, obstacles, x_faces, y_faces, z_faces, has_velocity_field, boundary_kinds, + shape, spacing, half_dt, signal_count, level_count); [encoder endEncoding]; wait_for_command(command_buffer, "Metal signal Jacobi iteration failed"); const auto sum = *static_cast(reduction.contents); @@ -1365,8 +1367,8 @@ void ensure_coupled_capacity(std::size_t cell_count, std::size_t species_level_c coupled_reaction_source_ = allocate_shared_buffer(device_, byte_count, "coupled affine sources"); coupled_reaction_loss_ = allocate_shared_buffer(device_, byte_count, "coupled affine losses"); - coupled_obstacles_ = allocate_shared_buffer(device_, coupled_grid_level_capacity_, - "coupled grid obstacles"); + coupled_obstacles_ = + allocate_shared_buffer(device_, coupled_grid_level_capacity_, "coupled grid obstacles"); } if (coupled_error_ == nil) { coupled_error_ = allocate_shared_buffer(device_, sizeof(std::uint32_t), "coupled error flag"); @@ -1736,8 +1738,7 @@ void fill_external_contacts(std::uint32_t cell_count, std::uint32_t constraint_c std::vector contacts; contacts.reserve(contact_count); for (std::uint32_t index = 0; index < contact_count; ++index) { - if (constraint_kinds[index] > - static_cast(ExternalConstraintKind::cylinder) || + if (constraint_kinds[index] > static_cast(ExternalConstraintKind::cylinder) || locations[index] > static_cast(RodContactLocation::interior)) { throw std::runtime_error("Metal external-contact kernel produced an invalid tag"); } @@ -1801,8 +1802,8 @@ void ensure_mechanics_capacity(std::size_t cell_count, std::size_t contact_count mechanics_incidence_offsets_ = allocate_shared_buffer(device_, (mechanics_cell_capacity_ + 1) * sizeof(std::uint32_t), "mechanics incidence offsets"); - mechanics_fixed_ = allocate_shared_buffer(device_, mechanics_cell_capacity_, - "mechanics fixed flags"); + mechanics_fixed_ = + allocate_shared_buffer(device_, mechanics_cell_capacity_, "mechanics fixed flags"); mechanics_solution_ = allocate_shared_buffer(device_, dof_bytes, "mechanics solution"); mechanics_rhs_ = allocate_shared_buffer(device_, dof_bytes, "mechanics right-hand side"); mechanics_residual_ = allocate_shared_buffer(device_, dof_bytes, "mechanics residual"); diff --git a/cpp/python/bindings.cpp b/cpp/python/bindings.cpp index 33dc5b7..2f7eba5 100644 --- a/cpp/python/bindings.cpp +++ b/cpp/python/bindings.cpp @@ -278,7 +278,7 @@ NB_MODULE(_core, module) { .def_prop_ro("instructions", [](const cm::SpeciesRatePlan& plan) { return std::vector(plan.instructions().begin(), - plan.instructions().end()); + plan.instructions().end()); }) .def_prop_ro("outputs", [](const cm::SpeciesRatePlan& plan) { @@ -297,7 +297,7 @@ NB_MODULE(_core, module) { .def_prop_ro("instructions", [](const cm::CoupledRatePlan& plan) { return std::vector(plan.instructions().begin(), - plan.instructions().end()); + plan.instructions().end()); }) .def_prop_ro("species_outputs", [](const cm::CoupledRatePlan& plan) { @@ -334,7 +334,7 @@ NB_MODULE(_core, module) { .def_prop_ro("contacts", [](const cm::ContactGraph& graph) { return std::vector(graph.contacts().begin(), - graph.contacts().end()); + graph.contacts().end()); }) .def("__len__", &cm::ContactGraph::size) .def( @@ -455,7 +455,7 @@ NB_MODULE(_core, module) { .def_prop_ro("contacts", [](const cm::ExternalContactGraph& graph) { return std::vector(graph.contacts().begin(), - graph.contacts().end()); + graph.contacts().end()); }) .def("__len__", &cm::ExternalContactGraph::size) .def( @@ -498,8 +498,8 @@ NB_MODULE(_core, module) { .def(nb::init(), "backend"_a = cm::BackendKind::cpu, "reserved_capacity"_a = 0, "species_count"_a = 0, "device_index"_a = 0) - .def(nb::init(), - "backend"_a, "checkpoint"_a, "device_index"_a = 0) + .def(nb::init(), "backend"_a, + "checkpoint"_a, "device_index"_a = 0) .def_prop_ro("backend_info", &cm::Simulation::backend_info) .def("supports", &cm::Simulation::supports, "feature"_a) .def_prop_ro("time", &cm::Simulation::time) @@ -511,6 +511,8 @@ NB_MODULE(_core, module) { .def_prop_ro("has_coupled_rate_plan", &cm::Simulation::has_coupled_rate_plan) .def("add_cell", &cm::Simulation::add_cell, "cell"_a) .def("remove_cell", &cm::Simulation::remove_cell, "id"_a) + .def("apply_flow_drift", &cm::Simulation::apply_flow_drift, "dt"_a, + "integration"_a = cm::MechanicsIntegrationParameters{}) .def("add_plane_constraint", &cm::Simulation::add_plane_constraint, "plane"_a) .def("add_sphere_constraint", &cm::Simulation::add_sphere_constraint, "sphere"_a) .def("add_box_constraint", &cm::Simulation::add_box_constraint, "box"_a) @@ -538,6 +540,7 @@ NB_MODULE(_core, module) { }, "levels"_a) .def("set_velocity_field", &cm::Simulation::set_velocity_field, "field"_a.none()) + .def("set_signal_reaction", &cm::Simulation::set_signal_reaction, "reaction"_a.none()) .def("divide", &cm::Simulation::divide, "parent_id"_a, "first_fraction"_a) .def("divide_equal", &cm::Simulation::divide_equal, "parent_id"_a) .def("step", &cm::Simulation::step, "dt"_a) diff --git a/docs/README.md b/docs/README.md index 3fe2841..767f6d4 100644 --- a/docs/README.md +++ b/docs/README.md @@ -82,3 +82,14 @@ cmake --preset cpu-debug cmake --build --preset cpu-debug ctest --preset cpu-debug ``` + +Formatting, linting, and the Python type check run as commit hooks: + +```console +uv run --with pre-commit pre-commit install +uv run --with pre-commit pre-commit run --all-files +``` + +The hooks cover the fast local gates only. Tests, native builds, and backend +conformance need a configured build and hardware, and run through CTest and the +[conformance scripts](development/validation.md). diff --git a/docs/architecture/0008-crank-nicolson-signals.md b/docs/architecture/0008-crank-nicolson-signals.md index 370c84c..f6af5f2 100644 --- a/docs/architecture/0008-crank-nicolson-signals.md +++ b/docs/architecture/0008-crank-nicolson-signals.md @@ -29,10 +29,15 @@ The linear system is solved by matrix-free Jacobi iteration over the local trans Convergence uses the RMS residual of the declared equation: ```text -residual_rms <= absolute_tolerance + relative_tolerance * rms(right_hand_side) +residual_rms <= max(absolute_tolerance, epsilon * rms(right_hand_side)) + + relative_tolerance * initial_residual_rms ``` -The grid specification records the integration kind, maximum iteration count, and both tolerances. These fields are exact checkpoint state. A solver that reaches the iteration limit or produces a non-finite residual fails the step; it never commits an unconverged field. A successful step exposes its iteration count and final residual. Final concentrations retain the engine-wide finite, non-negative invariant. Crank-Nicolson is stable for large diffusion steps but is not positivity preserving, so an oscillatory negative result is rejected. +The relative term scales the residual the step begins with, not the field it begins from. A field's magnitude says nothing about how much of it one step has to change, so scaling by the field makes the threshold grow with the background: a cell's exchange with a well-stocked field then falls under it and the solve returns the old field, discarding the source. Scaling by the initial residual instead asks for a fixed reduction of whatever this step actually has to resolve, so the accuracy a model gets is the accuracy it asked for, independent of concentration scale. + +The absolute term is raised to the residual floor of the field, `epsilon * rms(right_hand_side)` in float32. The right-hand side carries both the field and the step's operator terms, so this tracks stiffness as well as magnitude. A tolerance below the floor is unreachable, and clamping it up converges rather than iterating to the limit on rounding noise. The floor is also the resolution limit of an implicit step: a source that moves the field by less than its own float32 resolution leaves no residual to detect, and the step converges without it. A model whose exchange is that small next to its background belongs on forward Euler, which applies sources unconditionally, or on a concentration scale that resolves it. + +The grid specification records the integration kind, maximum iteration count, and both tolerances. These fields are exact checkpoint state, so a checkpoint written before this rule resumes with the tolerances it recorded and the current interpretation of them. A solver that reaches the iteration limit or produces a non-finite residual fails the step; it never commits an unconverged field. A successful step exposes its iteration count and final residual. Final concentrations retain the engine-wide finite, non-negative invariant. Crank-Nicolson is stable for large diffusion steps but is not positivity preserving, so an oscillatory negative result is rejected. ## Backend staging diff --git a/docs/architecture/0015-affine-grid-reactions.md b/docs/architecture/0015-affine-grid-reactions.md index 0c61e10..bc5d703 100644 --- a/docs/architecture/0015-affine-grid-reactions.md +++ b/docs/architecture/0015-affine-grid-reactions.md @@ -27,6 +27,19 @@ CPU, Metal, and CUDA implement the same fixed operation. Device implementations This is a generic data representation for one focused numerical operation, not a general voxel-program extension point. It does not support arbitrary functions of position or time, cross-signal reactions, nonlinear local kinetics, mutable masks, or runtime coefficient callbacks. A demonstrated need for one of those behaviors requires another named numerical contract with an explicit integration and checkpoint design. +## Runtime replacement + +A reaction field is data, so a model may replace it while a simulation runs: +`Simulation.set_signal_reaction` validates a candidate against the full grid +specification and swaps it atomically, exactly as a velocity field is swapped. +This is what lets a first-order loss that depends on cell state - an enzyme the +cells carry, degrading a signal in proportion to how much of it is there - be +carried by transport rather than scattered from the cells. Transport takes a +loss into its implicit diagonal, where a step stays positive while the loss +times the step is under two, whereas a cell-scattered source of the same +strength is explicit and needs half that step. Model code chooses the cadence, +and the field remains exact checkpoint state. + ## Consequences - Spatial reservoirs and first-order sinks remain inspectable and portable. diff --git a/docs/architecture/0021-flow-drift.md b/docs/architecture/0021-flow-drift.md new file mode 100644 index 0000000..a3ecfde --- /dev/null +++ b/docs/architecture/0021-flow-drift.md @@ -0,0 +1,63 @@ +# ADR 0021: advective flow drift on cells + +- Status: accepted +- Date: 2026-08-16 + +## Context + +The velocity field advects grid signals but exerts nothing on cells. In a flow-fed trap, +cells outside the trap should be carried downstream; in the overdamped regime a cell's +velocity relaxes to the local fluid velocity on a timescale far below one time step, so the +cell translates with the flow and a rod in shear rotates. + +## Decision + +`Simulation.apply_flow_drift(dt, integration)` advects every non-fixed cell through the grid's velocity +field by one explicit step, as an operation between growth and contact relaxation. The fluid +velocity is sampled at both capsule centerline endpoints: each stencil site's cell-centered +velocity is the mean of its two face velocities per axis, and the trilinear weights are the +signal-sampling weights, including the obstacle renormalization near walls. Endpoints are +clamped to the lattice of site centers before sampling: mechanics walls, not the lattice +edge, bound cells, so a rod whose tip pokes past the outermost site samples the nearest +in-grid point rather than erroring. Clamping happens in the lattice coordinate the bound is +tested in rather than in position space, so a clamped endpoint lies inside the lattice for +every origin and spacing. From endpoint velocities `v1` and `v2` with cylinder length `l` +and axis `a`: + +- translation is `dt * (v1 + v2) / 2`; +- the rotation vector is `dt * (a x (v2 - v1)) / l`, the least-squares rigid rotation for the + endpoint velocity difference, taken as zero when `l` is degenerate, applied as an + axis-angle rotation capped by the caller's mechanics integration rotation limit, through + the same axis-angle helper mechanical integration uses. + +Every update is validated before any world-state mutation, matching mechanical integration. +The operation requires a signal grid with a velocity field. A cell whose sampling stencil +holds no fluid samples zero rather than erroring: the field is validated zero on every face +of a solid site, so that is the field's own value there, and a cell that contact relaxation +has pressed into a wall simply does not drift. Concentration has no such value, so sampling +it inside an obstacle remains a model error. Fixed cells do not move. + +Drift composes with contact relaxation by operator splitting: drift first, then the ordinary +relaxation resolves any overlap the drift produced against walls or neighbors. The controller +applies drift when its mechanics configuration enables `flow_drift`, before the relaxation +passes. A formulation that couples drag into the relaxation right-hand side, so wall contact +forces balance fluid forces within one solve, is a candidate refinement with its own contract; +the explicit split is the reference behavior. + +The operation is host-side over committed state and identical on every backend; no kernel or +checkpoint change is involved. `flow_drift` is part of the controller's mechanics +configuration payload. + +## Validation sequence + +1. A free cell in uniform flow translates by exactly `velocity * dt` per drift call. +2. A rod spanning a shear gradient rotates toward alignment; a fixed cell does not move. +3. Drift against a wall followed by relaxation leaves the cell outside the wall. + +## Consequences + +- Washout becomes dynamic: flow carries cells to the removal predicate rather than the model + teleporting them. +- Splitting error is first order in `dt`; models choose steps so per-step drift stays small + relative to cell size, as they already do for growth. +- Cells in zero-velocity regions, including trap interiors, are unaffected. diff --git a/docs/architecture/README.md b/docs/architecture/README.md index bd74c2a..23fec13 100644 --- a/docs/architecture/README.md +++ b/docs/architecture/README.md @@ -15,6 +15,7 @@ The [numerical contract](numerical-contract.md) is the best starting point for w - [Axis-aligned box constraints](0016-box-constraints.md) - [Axis-aligned cylinder constraints](0018-cylinder-constraints.md) - [Cell removal](0020-cell-removal.md) +- [Advective flow drift on cells](0021-flow-drift.md) - [Steady Hele-Shaw-Brinkman flow solve](0022-brinkman-flow.md) - [Staggered MAC Stokes-Brinkman solve and flow benchmarks](0023-mac-stokes.md) diff --git a/docs/compatibility/legacy-example-migrations.md b/docs/compatibility/legacy-example-migrations.md index 33bf16a..f3bdbe9 100644 --- a/docs/compatibility/legacy-example-migrations.md +++ b/docs/compatibility/legacy-example-migrations.md @@ -26,6 +26,6 @@ Every translated model requests one exact mechanics pass per biological step. Th The legacy signaling callbacks expose extracellular derivatives as concentration rates and divide cell exchange by the `4 * 4 * 4 = 64` voxel volume before returning them. CellModeller2 coupled plans expose extracellular _amount_ rates; the native scatter operation performs the voxel-volume division. The migrated signal outputs therefore return the unscaled exchange amount. Intracellular equations that explicitly used `area / gridVolume` retain their division by 64. -The migrations use conventional diffusion coefficients rather than preserving the legacy implementation's accidental extra factor of one sixth. They retain the declared no-flux boundary behavior and integration choice: forward Euler for `ex3_simpleSignal.py`, and Crank-Nicolson for the other three models. The Crank-Nicolson ports set an absolute residual tolerance of `1e-12`; the engine default would accept a zero field while these models' initially small signal sources remained below its absolute threshold. +The migrations use conventional diffusion coefficients rather than preserving the legacy implementation's accidental extra factor of one sixth. They retain the declared no-flux boundary behavior and integration choice: forward Euler for `ex3_simpleSignal.py`, and Crank-Nicolson for the other three models. The Crank-Nicolson ports set an absolute residual tolerance of `1e-12`, which the engine raises to the float32 residual floor of the field being solved; convergence for these models is carried by the relative term against the residual each step starts with, so their initially small signal sources reach the field. `EdgeDetectorChamber.py` initialized `targetVol` but tested the absent `target_volume` attribute against 3.0. Its migration restores the evident intended uniform division threshold of 3.5 to 4.0. This is an explicit repair, not a claim that the original typo's behavior was reproduced. diff --git a/examples/culture_dish.py b/examples/culture_dish.py index 68754fe..a4483a9 100644 --- a/examples/culture_dish.py +++ b/examples/culture_dish.py @@ -1,7 +1,9 @@ """Colony growth confined to a round culture dish. The dish is a single inside-region cylinder constraint: its barrel is the circular -dish wall and its caps confine the colony to a monolayer. +dish wall and its caps hold the colony in a shallow layer. The caps leave a cell +diameter of play, so the colony stays essentially planar but is free to relieve +crowding in z rather than being pinned to one plane. """ from __future__ import annotations @@ -25,7 +27,7 @@ from cellmodeller2.checkpoint import CheckpointBundle, JSONValue MODEL_ID = "examples.culture-dish" -MODEL_VERSION = 1 +MODEL_VERSION = 2 DIVISION = UniformLengthDivision(3.2, 3.8, jitter_z=False) DISH_RADIUS = 30.0 @@ -59,7 +61,9 @@ def build(context: ModelContext) -> NativeController: founder_ids = [] for index in range(FOUNDER_COUNT): placement = context.rng.uniform(0.0, 2.0 * math.pi) - distance = context.rng.uniform(0.0, DISH_RADIUS * 0.5) + # The square root spreads founders uniformly over the seeded area + # rather than crowding them toward the middle. + distance = DISH_RADIUS * 0.5 * math.sqrt(context.rng.uniform(0.0, 1.0)) angle = context.rng.uniform(0.0, 2.0 * math.pi) founder = CellInit() founder.position = Vec3( diff --git a/python/src/cellmodeller2/_core.pyi b/python/src/cellmodeller2/_core.pyi index fffd835..dee19e3 100644 --- a/python/src/cellmodeller2/_core.pyi +++ b/python/src/cellmodeller2/_core.pyi @@ -579,6 +579,9 @@ class Simulation: def has_coupled_rate_plan(self) -> bool: ... def add_cell(self, cell: CellInit) -> int: ... def remove_cell(self, id: int) -> None: ... + def apply_flow_drift( + self, dt: float, integration: MechanicsIntegrationParameters = ... + ) -> None: ... def add_plane_constraint(self, plane: PlaneConstraintInit) -> int: ... def add_sphere_constraint(self, sphere: SphereConstraintInit) -> int: ... def add_box_constraint(self, box: BoxConstraintInit) -> int: ... @@ -595,6 +598,7 @@ class Simulation: def configure_signal_grid(self, spec: SignalGridSpec, levels: list[float] = ...) -> None: ... def set_signal_levels(self, levels: list[float]) -> None: ... def set_velocity_field(self, field: SignalGridVelocityField | None) -> None: ... + def set_signal_reaction(self, reaction: SignalGridAffineReaction | None) -> None: ... def divide(self, parent_id: int, first_fraction: float) -> tuple[int, int]: ... def divide_equal(self, parent_id: int) -> tuple[int, int]: ... def step(self, dt: float) -> None: ... diff --git a/python/src/cellmodeller2/controller.py b/python/src/cellmodeller2/controller.py index 3e15f4c..54b2844 100644 --- a/python/src/cellmodeller2/controller.py +++ b/python/src/cellmodeller2/controller.py @@ -185,6 +185,7 @@ class MechanicsConfig: require_convergence: bool = True constraint_activation_margin: float = 0.0 constraint_degeneracy_epsilon: float = 1.0e-6 + flow_drift: bool = False def __post_init__(self) -> None: _integer(self.passes, "mechanics.passes", 1, _UINT32_MAX) @@ -221,6 +222,8 @@ def __post_init__(self) -> None: raise ControllerStateError("mechanics constraint/contact parameters are invalid") if not isinstance(cast(object, self.require_convergence), bool): raise ControllerStateError("mechanics.require_convergence must be Boolean") + if not isinstance(cast(object, self.flow_drift), bool): + raise ControllerStateError("mechanics.flow_drift must be Boolean") def native_parameters( self, @@ -262,6 +265,9 @@ def from_json(cls, value: JSONValue) -> MechanicsConfig: require_convergence = value["require_convergence"] if not isinstance(require_convergence, bool): raise ControllerStateError("mechanics.require_convergence must be Boolean") + flow_drift = value["flow_drift"] + if not isinstance(flow_drift, bool): + raise ControllerStateError("mechanics.flow_drift must be Boolean") return cls( passes=_integer(value["passes"], "mechanics.passes", 1, _UINT32_MAX), mu_a=_finite_number(value["mu_a"], "mechanics.mu_a"), @@ -296,6 +302,7 @@ def from_json(cls, value: JSONValue) -> MechanicsConfig: value["constraint_degeneracy_epsilon"], "mechanics.constraint_degeneracy_epsilon", ), + flow_drift=flow_drift, ) @@ -461,6 +468,13 @@ def step(self, dt: float) -> None: self.simulation.step(dt) reports: list[MechanicsSolveResult] = [] + if ( + self._mechanics is not None + and self._mechanics.flow_drift + and self.simulation.cell_count != 0 + ): + _, _, integration, _ = self._mechanics.native_parameters() + self.simulation.apply_flow_drift(dt, integration) if self._mechanics is not None and self.simulation.cell_count != 0: parameters = self._mechanics.native_parameters() for _ in range(self._mechanics.passes): diff --git a/python/src/cellmodeller2/division.py b/python/src/cellmodeller2/division.py index 8e7d055..0b013d0 100644 --- a/python/src/cellmodeller2/division.py +++ b/python/src/cellmodeller2/division.py @@ -105,8 +105,11 @@ def on_division(self, step: ControllerStep, event: DivisionEvent) -> None: parent_key = str(event.parent.id) daughter_keys = {str(event.first.id), str(event.second.id)} active = {str(cell.id) for cell in step.cells} + # Plan removals apply after divisions, so cells already forgotten for + # this step's removals are still active here; targets may be a subset + # of the pre-division identities but never contain anything else. expected = (active - daughter_keys) | {parent_key} - if set(targets) != expected: + if parent_key not in targets or not set(targets) <= expected: raise ControllerStateError("division targets do not match pre-division identities") del targets[parent_key] if self.jitter_z is not None: diff --git a/python/src/cellmodeller2/py.typed b/python/src/cellmodeller2/py.typed index 8b13789..e69de29 100644 --- a/python/src/cellmodeller2/py.typed +++ b/python/src/cellmodeller2/py.typed @@ -1 +0,0 @@ - diff --git a/python/tests/test_controller.py b/python/tests/test_controller.py index e88e50e..a1c0704 100644 --- a/python/tests/test_controller.py +++ b/python/tests/test_controller.py @@ -18,9 +18,13 @@ ControllerStep, DivisionEvent, DivisionRequest, + GridBoundaryKind, + GridShape, MechanicsConfig, ModelContext, NativeController, + SignalGridSpec, + SignalGridVelocityField, Simulation, SimulationController, StepPlan, @@ -51,6 +55,37 @@ def _simulation_payload(path: Path) -> object: return document["simulation"] +def _one_cell_in_uniform_flow() -> Simulation: + sites = 9 + shape = GridShape() + shape.x, shape.y, shape.z = sites, 1, 1 + grid = SignalGridSpec() + grid.signal_count = 1 + grid.shape = shape + grid.spacing = Vec3(1.0, 1.0, 1.0) + grid.diffusion = [0.0] + grid.advection = [Vec3()] + grid.x_lower.kind = GridBoundaryKind.FIXED + grid.x_lower.values = [0.0] + grid.x_upper.kind = GridBoundaryKind.FIXED + grid.x_upper.values = [0.0] + velocity = SignalGridVelocityField() + velocity.x_faces = [0.5] * (sites + 1) + velocity.y_faces = [0.0] * (2 * sites) + velocity.z_faces = [0.0] * (2 * sites) + grid.velocity_field = velocity + + simulation = Simulation(BackendKind.CPU) + simulation.configure_signal_grid(grid) + cell = CellInit() + cell.position = Vec3(3.0, 0.0, 0.0) + cell.direction = Vec3(1.0, 0.0, 0.0) + cell.length = 1.0 + cell.radius = 0.25 + simulation.add_cell(cell) + return simulation + + def test_random_stream_round_trip_preserves_uniform_and_gaussian_draws() -> None: stream = random.Random(1729) stream.random() @@ -94,6 +129,61 @@ def invalidate_gaussian(value: dict[str, Any]) -> None: restore_random_state(value) +def test_mechanics_config_round_trip_preserves_flow_drift() -> None: + configuration = MechanicsConfig(passes=2, flow_drift=True) + + assert MechanicsConfig.from_json(configuration.to_json()) == configuration + + invalid = configuration.to_json() + invalid["flow_drift"] = 1 + with pytest.raises( + ControllerStateError, + match=r"mechanics\.flow_drift must be Boolean", + ): + MechanicsConfig.from_json(invalid) + + +def test_native_controller_resume_preserves_flow_drift_trajectory(tmp_path: Path) -> None: + def build() -> NativeController: + return NativeController( + _one_cell_in_uniform_flow(), + model_id="flow-drift-resume-test", + model_version=1, + rng=random.Random(7), + mechanics=MechanicsConfig(flow_drift=True), + ) + + uninterrupted = build() + for _ in range(4): + uninterrupted.step(0.25) + + split = build() + for _ in range(2): + split.step(0.25) + midpoint = tmp_path / "flow-midpoint.cm2.json" + save_checkpoint(split.simulation, midpoint, controller=split.controller_state()) + resumed = NativeController.from_checkpoint( + load_checkpoint_bundle(midpoint), + model_id="flow-drift-resume-test", + model_version=1, + ) + for _ in range(2): + resumed.step(0.25) + + expected = tmp_path / "flow-expected.cm2.json" + actual = tmp_path / "flow-actual.cm2.json" + save_checkpoint( + uninterrupted.simulation, + expected, + controller=uninterrupted.controller_state(), + ) + save_checkpoint(resumed.simulation, actual, controller=resumed.controller_state()) + + assert _simulation_payload(actual) == _simulation_payload(expected) + assert resumed.controller_state() == uninterrupted.controller_state() + assert math.isclose(resumed.simulation.cell(1).position.x, 3.5, abs_tol=1.0e-6) + + @pytest.mark.parametrize("backend", list(BackendKind)) def test_native_controller_composes_regulation_division_and_mechanics( backend: BackendKind, diff --git a/python/tests/test_division.py b/python/tests/test_division.py index 6902e28..1386c5d 100644 --- a/python/tests/test_division.py +++ b/python/tests/test_division.py @@ -13,6 +13,7 @@ Simulation, StepPlan, UniformLengthDivision, + Vec3, ) from cellmodeller2.checkpoint import JSONValue @@ -49,6 +50,46 @@ def regulate(step: ControllerStep) -> StepPlan: assert all(cell.direction.z == 0.0 for cell in daughters) +def test_division_and_removal_in_the_same_plan_coexist() -> None: + # Plan removals apply after divisions, so a forgotten (about-to-be-removed) + # cell is still active while division callbacks run. + simulation = Simulation(BackendKind.CPU) + divider = CellInit() + divider.length = 4.0 + divider.position = Vec3(0.0, 5.0, 0.0) + divider_id = simulation.add_cell(divider) + leaver = CellInit() + leaver.length = 1.0 + leaver_id = simulation.add_cell(leaver) + stream = random.Random(11) + state: dict[str, JSONValue] = {} + policy = UniformLengthDivision(2.0, 2.5, jitter_z=False) + policy.initialize(state, stream, (divider_id, leaver_id)) + + def regulate(step: ControllerStep) -> StepPlan: + divisions = policy.requests(step) + policy.forget(step, (leaver_id,)) + return StepPlan(divisions=divisions, removals=(leaver_id,)) + + controller = NativeController( + simulation, + model_id="division-with-removal-test", + model_version=1, + rng=stream, + regulate=regulate, + on_division=policy.on_division, + state=state, + ) + controller.step(0.0) + + targets_state = cast(dict[str, JSONValue], controller.state["length_division"]) + targets = cast(dict[str, JSONValue], targets_state["targets"]) + daughters = simulation.cells() + assert len(daughters) == 2 + assert set(targets) == {str(cell.id) for cell in daughters} + assert all(simulation.lineage_parent(cell.id) == divider_id for cell in daughters) + + def test_uniform_length_division_rejects_missing_target_state() -> None: simulation = Simulation() simulation.add_cell(CellInit()) diff --git a/python/tests/test_flow.py b/python/tests/test_flow.py index f0649cf..4a1f6d1 100644 --- a/python/tests/test_flow.py +++ b/python/tests/test_flow.py @@ -9,6 +9,7 @@ BackendKind, GridBoundaryKind, GridShape, + SignalGridAffineReaction, SignalGridSpec, SignalGridVelocityField, SignalIntegrationKind, @@ -16,7 +17,13 @@ Vec3, backend_available, ) -from cellmodeller2.flow import FlowError, colony_mobility, gap_mobility, solve_flow_field +from cellmodeller2.flow import ( + FlowError, + colony_mobility, + colony_species_density, + gap_mobility, + solve_flow_field, +) from cellmodeller2.microfluidics import TrapChannelDevice @@ -376,3 +383,74 @@ def test_partly_blocked_inlets_and_walled_off_pockets_solve() -> None: # The pocket is cut off from the flow, so it carries none. assert sealed_field.y_faces[_y_face(pocket, 3, 3, 0)] == 0.0 assert sealed_field.y_faces[_y_face(pocket, 4, 3, 0)] == 0.0 + + +def test_colony_species_density_rasterizes_one_channel() -> None: + """A per-cell rate becomes a per-volume rate through the voxel it sits in.""" + + spec = _duct(nx=3, ny=3, nz=1) + spec.spacing = Vec3(4.0, 4.0, 4.0) + + @dataclass + class _Cell: + position: Vec3 + species: list[float] + + cells = [ + _Cell(Vec3(0.0, 0.0, 0.0), [1.0, 2.0]), + _Cell(Vec3(0.0, 0.0, 0.0), [3.0, 4.0]), + _Cell(Vec3(4.0, 4.0, 0.0), [0.0, 5.0]), + _Cell(Vec3(400.0, 0.0, 0.0), [9.0, 9.0]), + ] + density = colony_species_density(spec, cells, species=1) + voxel = spec.voxel_volume + assert math.isclose(density[_site(spec, 0, 0, 0)], 6.0 / voxel) + assert math.isclose(density[_site(spec, 1, 1, 0)], 5.0 / voxel) + assert density[_site(spec, 2, 2, 0)] == 0.0 + with pytest.raises(FlowError, match="outside the cell"): + colony_species_density(spec, cells, species=5) + + +def test_a_running_simulation_swaps_its_signal_reaction() -> None: + """A loss that depends on cell state can be handed to transport each step.""" + + spec = _duct(nx=3, ny=1, nz=1) + spec.x_lower.kind = GridBoundaryKind.NO_FLUX + spec.x_lower.values = [] + spec.x_upper.kind = GridBoundaryKind.NO_FLUX + spec.x_upper.values = [] + for name in ("y_lower", "y_upper"): + boundary = getattr(spec, name) + boundary.kind = GridBoundaryKind.NO_FLUX + boundary.values = [] + setattr(spec, name, boundary) + spec.diffusion = [0.0] + spec.integration = SignalIntegrationKind.CRANK_NICOLSON + + simulation = Simulation() + simulation.configure_signal_grid(spec, [4.0, 4.0, 4.0]) + simulation.step(0.5) + assert simulation.signal_levels == [4.0, 4.0, 4.0] + + reaction = SignalGridAffineReaction() + reaction.source_rates = [0.0, 0.0, 0.0] + reaction.loss_rates = [2.0, 0.0, 0.0] + simulation.set_signal_reaction(reaction) + simulation.step(0.5) + levels = simulation.signal_levels + assert levels[0] < 4.0 + assert levels[1] == 4.0 + + # Transport takes the reaction into its implicit diagonal, so a loss that + # would remove more than a site holds in one explicit step is still stable. + # Crank-Nicolson stays positive while the loss times the step is under two. + reaction.loss_rates = [3.0, 0.0, 0.0] + simulation.set_signal_reaction(reaction) + before_strong = simulation.signal_levels[0] + simulation.step(0.5) + assert 0.0 < simulation.signal_levels[0] < before_strong * 0.5 + + simulation.set_signal_reaction(None) + before = simulation.signal_levels + simulation.step(0.5) + assert simulation.signal_levels == before diff --git a/python/tests/test_signals.py b/python/tests/test_signals.py index 45005c8..2d4fdfd 100644 --- a/python/tests/test_signals.py +++ b/python/tests/test_signals.py @@ -13,6 +13,7 @@ GridShape, RateInstruction, RateOp, + RatePlanBuilder, SignalGridSpec, SignalIntegrationKind, Simulation, @@ -73,6 +74,11 @@ def test_cpu_signal_transport_sampling_and_stability() -> None: def test_cpu_crank_nicolson_accepts_a_step_beyond_the_euler_bound() -> None: spec = _line_spec() spec.integration = SignalIntegrationKind.CRANK_NICOLSON + # The relative tolerance asks for a reduction of the residual the step + # starts with, so it sets how closely the committed field approaches the + # exact one: this solve is checked to a millionth, so it asks for rather + # better than that. + spec.solver.relative_tolerance = 1.0e-7 simulation = Simulation() simulation.configure_signal_grid(spec, [0.0, 1.0, 0.0]) @@ -86,6 +92,80 @@ def test_cpu_crank_nicolson_accepts_a_step_beyond_the_euler_bound() -> None: assert report.residual_rms <= 2.0e-5 +def _uptake_removed(background: float, integration: SignalIntegrationKind) -> float: + """Total signal a single growing cell removes from a uniform field in one step.""" + + rates = RatePlanBuilder() + uptake = -(rates.growth_rate() * rates.cell_volume()) + plan = rates.coupled_plan(0, 1, (), (uptake,)) + + shape = GridShape() + shape.x, shape.y, shape.z = 4, 4, 4 + spec = SignalGridSpec() + spec.signal_count = 1 + spec.shape = shape + spec.spacing = Vec3(4.0, 4.0, 4.0) + spec.diffusion = [0.0] + spec.advection = [Vec3()] + spec.integration = integration + + simulation = Simulation() + simulation.configure_signal_grid(spec, [background] * 64) + simulation.set_coupled_rate_plan(plan) + cell = CellInit() + cell.position = Vec3(4.0, 4.0, 4.0) + cell.length = 2.6 + cell.radius = 0.5 + cell.growth_rate = 0.667 + simulation.add_cell(cell) + + before = sum(simulation.signal_levels) + simulation.step(0.02) + return (before - sum(simulation.signal_levels)) * spec.voxel_volume + + +@pytest.mark.parametrize("background", [1.0, 10.0, 100.0]) +def test_cell_sources_reach_the_grid_at_any_resolvable_background(background: float) -> None: + """A cell's exchange with the field must survive the convergence test. + + An implicit step is accepted on a residual, and a cell's contribution is + small next to a well-stocked field. Judging that residual against the field + would let the contribution fall under the threshold and be dropped, leaving + the model silently inert, so it is judged against the residual the step + starts with instead. Forward Euler applies its sources unconditionally and + is the reference here. + """ + + explicit = _uptake_removed(background, SignalIntegrationKind.FORWARD_EULER) + implicit = _uptake_removed(background, SignalIntegrationKind.CRANK_NICOLSON) + assert math.isclose(implicit, explicit, rel_tol=1.0e-3) + assert implicit > 0.0 + + +def test_a_source_below_the_field_noise_commits_a_converged_step() -> None: + """The limit of an implicit solve: a source it cannot see is not an error. + + Convergence is decided on a residual computed in float32, so a source that + moves the field by less than its own representable resolution leaves no + residual to detect. The step still converges and commits rather than + failing; a model whose exchange is that small next to its background needs + forward Euler or a concentration scale that resolves it. + """ + + removed = _uptake_removed(1.0e4, SignalIntegrationKind.CRANK_NICOLSON) + assert removed == 0.0 + + +def test_cell_sources_match_their_declared_amount() -> None: + cell_volume = math.pi * 0.5**2 * (2.6 + 2.0 * 0.5) + expected = 0.02 * 0.667 * cell_volume + for integration in ( + SignalIntegrationKind.FORWARD_EULER, + SignalIntegrationKind.CRANK_NICOLSON, + ): + assert math.isclose(_uptake_removed(10.0, integration), expected, rel_tol=0.02) + + def test_fixed_and_periodic_boundaries_are_explicit() -> None: fixed = _line_spec(2) lower = GridBoundary() diff --git a/python/tests/test_simulation.py b/python/tests/test_simulation.py index f22d776..327d096 100644 --- a/python/tests/test_simulation.py +++ b/python/tests/test_simulation.py @@ -12,6 +12,8 @@ ConstraintRegion, ContactParameters, ExternalConstraintKind, + GridBoundaryKind, + GridShape, MechanicsIntegrationParameters, MechanicsParameters, PlaneConstraintInit, @@ -19,6 +21,8 @@ RateOp, RodContactLocation, RodEndpoint, + SignalGridSpec, + SignalGridVelocityField, Simulation, SolverBreakdown, SolverStatus, @@ -466,3 +470,131 @@ def test_native_growth_matches_cpu(backend: BackendKind) -> None: for cpu_cell, native_cell in zip(cpu.cells(), native.cells(), strict=True): assert math.isclose(cpu_cell.length, native_cell.length, abs_tol=1.0e-6) + + +def uniform_flow_grid( + *, origin: float, spacing: float, sites: int, speed: float +) -> SignalGridSpec: + """A collapsed y/z lattice carrying a uniform x flow between fixed ends.""" + + shape = GridShape() + shape.x, shape.y, shape.z = sites, 1, 1 + grid = SignalGridSpec() + grid.signal_count = 1 + grid.shape = shape + grid.origin = Vec3(origin, 0.0, 0.0) + grid.spacing = Vec3(spacing, 1.0, 1.0) + grid.diffusion = [0.0] + grid.advection = [Vec3()] + grid.x_lower.kind = GridBoundaryKind.FIXED + grid.x_lower.values = [0.0] + grid.x_upper.kind = GridBoundaryKind.FIXED + grid.x_upper.values = [0.0] + field = SignalGridVelocityField() + field.x_faces = [speed] * (sites + 1) + field.y_faces = [0.0] * (2 * sites) + field.z_faces = [0.0] * (2 * sites) + grid.velocity_field = field + return grid + + +@pytest.mark.parametrize( + ("origin", "spacing"), [(0.0, 1.0), (0.1, 0.3), (-97.5, 1.65), (0.7, 5.0)] +) +def test_flow_drift_clamps_endpoints_on_any_lattice(origin: float, spacing: float) -> None: + sites = 33 + simulation = Simulation() + simulation.configure_signal_grid(uniform_flow_grid( + origin=origin, spacing=spacing, sites=sites, speed=2.0 + )) + cell = CellInit() + cell.position = Vec3(origin + spacing * (sites - 1), 0.0, 0.0) + cell.direction = Vec3(1.0, 0.0, 0.0) + cell.length = 2.0 * spacing + cell.radius = 0.3 + cell_id = simulation.add_cell(cell) + + simulation.apply_flow_drift(0.1) + + assert math.isclose( + simulation.cell(cell_id).position.x, + origin + spacing * (sites - 1) + 0.2, + rel_tol=1.0e-5, + abs_tol=1.0e-5, + ) + + +def test_flow_drift_honors_the_mechanics_rotation_limit() -> None: + shape = GridShape() + shape.x, shape.y, shape.z = 3, 3, 1 + grid = SignalGridSpec() + grid.signal_count = 1 + grid.shape = shape + grid.spacing = Vec3(1.0, 1.0, 1.0) + grid.diffusion = [0.0] + grid.advection = [Vec3()] + grid.x_lower.kind = GridBoundaryKind.FIXED + grid.x_lower.values = [0.0] + grid.x_upper.kind = GridBoundaryKind.FIXED + grid.x_upper.values = [0.0] + field = SignalGridVelocityField() + field.x_faces = [float(y) for _ in range(4) for y in range(3)] + field.y_faces = [0.0] * 12 + field.z_faces = [0.0] * 18 + grid.velocity_field = field + + cell = CellInit() + cell.position = Vec3(1.0, 1.0, 0.0) + cell.direction = Vec3(0.0, 1.0, 0.0) + cell.length = 2.0 + cell.radius = 0.3 + + capped = Simulation() + capped.configure_signal_grid(grid) + capped_id = capped.add_cell(cell) + capped.apply_flow_drift(1.0) + limit = MechanicsIntegrationParameters().max_rotation_radians + assert math.isclose(capped.cell(capped_id).direction.x, math.sin(limit), abs_tol=1.0e-6) + + frozen = Simulation() + frozen.configure_signal_grid(grid) + frozen_id = frozen.add_cell(cell) + integration = MechanicsIntegrationParameters() + integration.max_rotation_radians = 0.0 + frozen.apply_flow_drift(1.0, integration) + assert math.isclose(frozen.cell(frozen_id).direction.x, 0.0, abs_tol=1.0e-6) + assert math.isclose(frozen.cell(frozen_id).direction.y, 1.0, abs_tol=1.0e-6) + + +def test_a_cell_inside_a_wall_samples_no_flow_and_does_not_drift() -> None: + """Mechanics can press a crowded cell into a wall; drift must survive it. + + The velocity field is zero on every face of a solid site, so a stencil with + no fluid in it samples exactly zero and the cell stays put. Concentration + there has no such value, so sampling it is still an error. + """ + + sites = 5 + grid = uniform_flow_grid(origin=0.0, spacing=1.0, sites=sites, speed=2.0) + # A solid block in the middle of the line, with the flow stopping at it. + grid.obstacles = [0, 0, 1, 0, 0] + field = SignalGridVelocityField() + field.x_faces = [2.0, 2.0, 0.0, 0.0, 2.0, 2.0] + field.y_faces = [0.0] * (2 * sites) + field.z_faces = [0.0] * (2 * sites) + grid.velocity_field = field + + simulation = Simulation() + simulation.configure_signal_grid(grid, [1.0, 1.0, 0.0, 1.0, 1.0]) + buried = CellInit() + buried.position = Vec3(2.0, 0.0, 0.0) + buried.direction = Vec3(1.0, 0.0, 0.0) + buried.length = 0.0 + buried.radius = 0.3 + buried_id = simulation.add_cell(buried) + + simulation.apply_flow_drift(0.5) + + assert simulation.cell(buried_id).position.x == 2.0 + with pytest.raises(ValueError, match="inside a grid obstacle"): + simulation.sample_signals(Vec3(2.0, 0.0, 0.0)) diff --git a/tests/cpp/signal_grid_test.cpp b/tests/cpp/signal_grid_test.cpp index 61b3a38..0f0a1e9 100644 --- a/tests/cpp/signal_grid_test.cpp +++ b/tests/cpp/signal_grid_test.cpp @@ -368,6 +368,72 @@ int main() { assert_throws([&] { invalid.validate(); }); } + { + auto spec = line_spec(3); + spec.diffusion = {0.0F}; + spec.x_lower.kind = cm::GridBoundaryKind::fixed; + spec.x_lower.values = {0.0F}; + spec.x_upper.kind = cm::GridBoundaryKind::fixed; + spec.x_upper.values = {0.0F}; + spec.velocity_field = cm::SignalGridVelocityField{ + .x_faces = {1.0F, 1.0F, 1.0F, 1.0F}, + .y_faces = std::vector(6, 0.0F), + .z_faces = std::vector(6, 0.0F), + }; + cm::Simulation simulation; + simulation.configure_signal_grid(spec); + cm::CellInit mover; + mover.position = {1.0F, 0.0F, 0.0F}; + mover.length = 0.0F; + mover.radius = 0.4F; + const auto mover_id = simulation.add_cell(mover); + cm::CellInit anchored = mover; + anchored.position = {0.5F, 0.0F, 0.0F}; + anchored.fixed = true; + const auto anchored_id = simulation.add_cell(anchored); + + simulation.apply_flow_drift(0.25F); + assert_close(simulation.cell(mover_id).position.x, 1.25F); + assert_close(simulation.cell(anchored_id).position.x, 0.5F); + } + + { + cm::SignalGridSpec spec; + spec.signal_count = 1; + spec.shape = {.x = 3, .y = 3, .z = 1}; + spec.diffusion = {0.0F}; + spec.advection = {{0.0F, 0.0F, 0.0F}}; + spec.x_lower.kind = cm::GridBoundaryKind::fixed; + spec.x_lower.values = {0.0F}; + spec.x_upper.kind = cm::GridBoundaryKind::fixed; + spec.x_upper.values = {0.0F}; + std::vector x_faces(12, 0.0F); + for (std::uint32_t fx = 0; fx < 4; ++fx) { + for (std::uint32_t y = 0; y < 3; ++y) { + x_faces[(fx * 3) + y] = 0.5F * static_cast(y); + } + } + spec.velocity_field = cm::SignalGridVelocityField{ + .x_faces = x_faces, + .y_faces = std::vector(12, 0.0F), + .z_faces = std::vector(18, 0.0F), + }; + cm::Simulation simulation; + simulation.configure_signal_grid(spec); + cm::CellInit rod; + rod.position = {1.0F, 1.0F, 0.0F}; + rod.direction = {0.0F, 1.0F, 0.0F}; + rod.length = 1.0F; + rod.radius = 0.3F; + const auto rod_id = simulation.add_cell(rod); + + simulation.apply_flow_drift(0.1F); + const auto drifted = simulation.cell(rod_id); + assert(drifted.position.x > 1.0F); + assert(drifted.direction.x > 0.04F); + assert_close(cm::norm(drifted.direction), 1.0F); + } + { auto spec = line_spec(3); spec.diffusion = {0.0F}; @@ -377,11 +443,14 @@ int main() { spec.x_upper.values = {0.0F}; cm::SignalGrid grid(spec); + assert_throws( + [&grid] { (void)grid.sample_velocity({1.5F, 0.0F, 0.0F}); }); grid.set_velocity_field(cm::SignalGridVelocityField{ .x_faces = std::vector(4, 2.0F), .y_faces = std::vector(6, 0.0F), .z_faces = std::vector(6, 0.0F), }); + assert_close(grid.sample_velocity({1.5F, 0.0F, 0.0F}).x, 2.0F); assert_throws([&grid] { grid.set_velocity_field(cm::SignalGridVelocityField{ .x_faces = std::vector(5, 0.0F), @@ -389,18 +458,116 @@ int main() { .z_faces = std::vector(6, 0.0F), }); }); + assert_close(grid.sample_velocity({1.5F, 0.0F, 0.0F}).x, 2.0F); grid.set_velocity_field(std::nullopt); + assert_throws( + [&grid] { (void)grid.sample_velocity({1.5F, 0.0F, 0.0F}); }); cm::Simulation simulation; simulation.configure_signal_grid(spec); + cm::CellInit drifter; + drifter.position = {1.5F, 0.0F, 0.0F}; + drifter.direction = {1.0F, 0.0F, 0.0F}; + drifter.length = 1.0F; + drifter.radius = 0.3F; + const auto drifter_id = simulation.add_cell(drifter); + assert_throws([&simulation] { simulation.apply_flow_drift(0.1F); }); + assert_close(simulation.cell(drifter_id).position.x, 1.5F); simulation.set_velocity_field(cm::SignalGridVelocityField{ .x_faces = std::vector(4, 2.0F), .y_faces = std::vector(6, 0.0F), .z_faces = std::vector(6, 0.0F), }); - simulation.set_velocity_field(std::nullopt); + simulation.apply_flow_drift(0.1F); + assert_close(simulation.cell(drifter_id).position.x, 1.7F); cm::Simulation bare; assert_throws([&bare] { bare.set_velocity_field(std::nullopt); }); + + // A rod whose tip pokes past the outermost site center drifts by sampling + // the nearest in-grid point instead of erroring. + cm::CellInit poking; + poking.position = {2.5F, 0.0F, 0.0F}; + poking.direction = {1.0F, 0.0F, 0.0F}; + poking.length = 1.5F; + poking.radius = 0.3F; + const auto poking_id = simulation.add_cell(poking); + simulation.apply_flow_drift(0.1F); + assert_close(simulation.cell(poking_id).position.x, 2.7F); + } + + { + // Endpoint clamping works in lattice coordinates, so an origin and spacing + // with no exact float representation still admits a rod poking past the + // outermost site center. + auto spec = line_spec(33); + spec.origin = {0.1F, 0.0F, 0.0F}; + spec.spacing = {0.3F, 1.0F, 1.0F}; + spec.diffusion = {0.0F}; + spec.x_lower.kind = cm::GridBoundaryKind::fixed; + spec.x_lower.values = {0.0F}; + spec.x_upper.kind = cm::GridBoundaryKind::fixed; + spec.x_upper.values = {0.0F}; + spec.velocity_field = cm::SignalGridVelocityField{ + .x_faces = std::vector(34, 2.0F), + .y_faces = std::vector(66, 0.0F), + .z_faces = std::vector(66, 0.0F), + }; + cm::Simulation simulation; + simulation.configure_signal_grid(spec); + cm::CellInit poking; + poking.position = {9.7F, 0.0F, 0.0F}; + poking.direction = {1.0F, 0.0F, 0.0F}; + poking.length = 2.0F; + poking.radius = 0.3F; + const auto poking_id = simulation.add_cell(poking); + simulation.apply_flow_drift(0.1F); + assert_close(simulation.cell(poking_id).position.x, 9.9F); + } + + { + // A rod spanning a shear gradient rotates toward the flow, capped by the + // caller's mechanics rotation limit. + cm::SignalGridSpec spec; + spec.signal_count = 1; + spec.shape = {.x = 3, .y = 3, .z = 1}; + spec.diffusion = {0.0F}; + spec.advection = {{0.0F, 0.0F, 0.0F}}; + spec.x_lower.kind = cm::GridBoundaryKind::fixed; + spec.x_lower.values = {0.0F}; + spec.x_upper.kind = cm::GridBoundaryKind::fixed; + spec.x_upper.values = {0.0F}; + std::vector x_faces(12, 0.0F); + for (std::uint32_t fx = 0; fx < 4; ++fx) { + for (std::uint32_t y = 0; y < 3; ++y) { + x_faces[(fx * 3) + y] = static_cast(y); + } + } + spec.velocity_field = cm::SignalGridVelocityField{ + .x_faces = x_faces, + .y_faces = std::vector(12, 0.0F), + .z_faces = std::vector(18, 0.0F), + }; + + cm::CellInit rod; + rod.position = {1.0F, 1.0F, 0.0F}; + rod.direction = {0.0F, 1.0F, 0.0F}; + rod.length = 2.0F; + rod.radius = 0.3F; + + cm::Simulation capped; + capped.configure_signal_grid(spec); + const auto capped_id = capped.add_cell(rod); + capped.apply_flow_drift(1.0F); + const auto limit = cm::MechanicsIntegrationParameters{}.max_rotation_radians; + assert_close(capped.cell(capped_id).direction.x, std::sin(limit)); + assert_close(capped.cell(capped_id).direction.y, std::cos(limit)); + + cm::Simulation frozen; + frozen.configure_signal_grid(spec); + const auto frozen_id = frozen.add_cell(rod); + frozen.apply_flow_drift(1.0F, cm::MechanicsIntegrationParameters{.max_rotation_radians = 0.0F}); + assert_close(frozen.cell(frozen_id).direction.x, 0.0F); + assert_close(frozen.cell(frozen_id).direction.y, 1.0F); } }