Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
50 changes: 47 additions & 3 deletions apps/tests/core/scene.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -293,9 +293,7 @@ TEST_CASE("scene_commit", "[scene]")
std::ranges::transform(pos_view.begin(),
pos_view.end(),
pos_view.begin(),
[](const auto& p) {
return p + Vector3{1, 1, 1};
});
[](const auto& p) { return p + Vector3{1, 1, 1}; });

SceneSnapshotCommit commit = scene - snapshot;
REQUIRE(commit.contact_models().attribute_collection().find("topo") == nullptr);
Expand Down Expand Up @@ -358,3 +356,49 @@ TEST_CASE("scene_commit_empty", "[scene]")
.find("myattribute")
!= nullptr);
}


TEST_CASE("dynamic_dt_config", "[scene]")
{
using namespace uipc;
using namespace uipc::core;
using namespace uipc::geometry;

Scene scene;

// Default dt should be 0.01
{
auto dt_slot = scene.config().find<Float>("dt");
REQUIRE(dt_slot != nullptr);
REQUIRE(dt_slot->view()[0] == Catch::Approx(0.01));
}

// Change dt via config
{
auto dt_slot = scene.config().find<Float>("dt");
view(*dt_slot)[0] = 0.005;
}

// Read back via config
{
auto dt_slot = scene.config().find<Float>("dt");
REQUIRE(dt_slot->view()[0] == Catch::Approx(0.005));
}

// Read back via SceneVisitor
{
backend::SceneVisitor visitor{scene};
REQUIRE(visitor.dt() == Catch::Approx(0.005));
}

// Change again and verify
{
auto dt_slot = scene.config().find<Float>("dt");
view(*dt_slot)[0] = 0.02;
}

{
backend::SceneVisitor visitor{scene};
REQUIRE(visitor.dt() == Catch::Approx(0.02));
}
}
175 changes: 175 additions & 0 deletions apps/tests/sim_case/93_dynamic_dt.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,175 @@
#include <app/app.h>
#include <uipc/uipc.h>
#include <uipc/constitution/affine_body_constitution.h>

TEST_CASE("93_dynamic_dt", "[abd][dynamic_dt]")
{
using namespace uipc;
using namespace uipc::core;
using namespace uipc::geometry;
using namespace uipc::constitution;

auto output_path = AssetDir::output_path(UIPC_RELATIVE_SOURCE_FILE);

Engine engine{"cuda", output_path};
World world{engine};

auto config = test::Scene::default_config();
config["gravity"] = Vector3{0, -9.8, 0};
config["contact"]["enable"] = false;
config["dt"] = 0.01;
test::Scene::dump_config(config, output_path);

Scene scene{config};
{
AffineBodyConstitution abd;

auto object = scene.objects().create("falling_tet");

vector<Vector4i> Ts = {Vector4i{0, 1, 2, 3}};
vector<Vector3> Vs = {Vector3{0, 1, 0},
Vector3{0, 0, 1},
Vector3{-std::sqrt(3) / 2, 0, -0.5},
Vector3{std::sqrt(3) / 2, 0, -0.5}};

std::transform(
Vs.begin(), Vs.end(), Vs.begin(), [](auto& v) { return v * 0.3; });

auto mesh = tetmesh(Vs, Ts);
label_surface(mesh);
label_triangle_orient(mesh);

mesh.instances().resize(1);
abd.apply_to(mesh, 100.0_MPa);

auto trans_view = view(mesh.transforms());
auto is_fixed = mesh.instances().find<IndexT>(builtin::is_fixed);
auto is_fixed_view = view(*is_fixed);

Transform t = Transform::Identity();
t.translation() = Vector3::UnitY() * 2;
trans_view[0] = t.matrix();
is_fixed_view[0] = 0;

object->geometries().create(mesh);
}

world.init(scene);
REQUIRE(world.is_valid());

auto geo_slots = scene.geometries().find(0);
auto geo_slot = geo_slots.geometry;

auto read_center_y = [&]() -> Float
{
auto geo = geo_slot->geometry().as<SimplicialComplex>();
auto trans_view = geo->transforms().view();
Matrix4x4 T = trans_view[0];
return T(1, 3);
};

Float prev_y = read_center_y();

// Run frame 1 with dt=0.01
world.advance();
REQUIRE(world.is_valid());
world.retrieve();
Float y_after_1 = read_center_y();
Float dy_frame1 = prev_y - y_after_1;
REQUIRE(dy_frame1 > 0.0);

// Now change dt to 0.005 and reset to the same starting position
// by running a fresh scene. Instead, we compare the first-frame displacement
// at two different dt values. Run additional frames to let it settle.
//
// Simpler approach: compare dy/dt ratio.
// At frame 1 with dt=0.01: dy1 = 0.5 * g * dt^2 = 0.5 * 9.8 * 0.0001 = 0.00049
// With dt=0.005: dy = 0.5 * g * dt^2 = 0.5 * 9.8 * 0.000025 = 0.0001225
// ratio should be ~4x

// Reset: create a second simulation with dt=0.005 from the start
Engine engine2{"cuda", output_path};
World world2{engine2};

auto config2 = test::Scene::default_config();
config2["gravity"] = Vector3{0, -9.8, 0};
config2["contact"]["enable"] = false;
config2["dt"] = 0.005;

Scene scene2{config2};
{
AffineBodyConstitution abd2;

auto object2 = scene2.objects().create("falling_tet");

vector<Vector4i> Ts2 = {Vector4i{0, 1, 2, 3}};
vector<Vector3> Vs2 = {Vector3{0, 1, 0},
Vector3{0, 0, 1},
Vector3{-std::sqrt(3) / 2, 0, -0.5},
Vector3{std::sqrt(3) / 2, 0, -0.5}};

std::transform(
Vs2.begin(), Vs2.end(), Vs2.begin(), [](auto& v) { return v * 0.3; });

auto mesh2 = tetmesh(Vs2, Ts2);
label_surface(mesh2);
label_triangle_orient(mesh2);

mesh2.instances().resize(1);
abd2.apply_to(mesh2, 100.0_MPa);

auto trans_view2 = view(mesh2.transforms());
auto is_fixed2 = mesh2.instances().find<IndexT>(builtin::is_fixed);
auto is_fixed_view2 = view(*is_fixed2);

Transform t2 = Transform::Identity();
t2.translation() = Vector3::UnitY() * 2;
trans_view2[0] = t2.matrix();
is_fixed_view2[0] = 0;

object2->geometries().create(mesh2);
}

world2.init(scene2);
REQUIRE(world2.is_valid());

auto geo_slots2 = scene2.geometries().find(0);
auto geo_slot2 = geo_slots2.geometry;

auto read_center_y2 = [&]() -> Float
{
auto geo2 = geo_slot2->geometry().as<SimplicialComplex>();
auto trans_view = geo2->transforms().view();
Matrix4x4 T = trans_view[0];
return T(1, 3);
};

Float prev_y2 = read_center_y2();

world2.advance();
REQUIRE(world2.is_valid());
world2.retrieve();
Float y_after_1_small = read_center_y2();
Float dy_frame1_small = prev_y2 - y_after_1_small;
REQUIRE(dy_frame1_small > 0.0);

// With dt halved, first-frame displacement (0.5*g*dt^2) should be ~4x smaller
REQUIRE(dy_frame1_small < dy_frame1);

// Now test dynamic dt change: change dt in scene2 from 0.005 to 0.01
{
auto dt_slot = scene2.config().find<Float>("dt");
view(*dt_slot)[0] = 0.01;
}

// Run one more frame with the changed dt
Float prev_y2_before = read_center_y2();
world2.advance();
REQUIRE(world2.is_valid());
world2.retrieve();
Float y_after_2_changed = read_center_y2();
Float dy_frame2_changed = prev_y2_before - y_after_2_changed;

// After changing dt from 0.005 to 0.01, the displacement should increase
REQUIRE(dy_frame2_changed > dy_frame1_small);
}
1 change: 1 addition & 0 deletions include/uipc/backend/visitors/scene_visitor.h
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ class UIPC_CORE_API SceneVisitor

span<IndexT> pending_destroy_ids() const noexcept;
const geometry::AttributeCollection& config() const noexcept;
Float dt() const noexcept;

const core::ConstitutionTabular& constitution_tabular() const noexcept;
core::ConstitutionTabular& constitution_tabular() noexcept;
Expand Down
14 changes: 8 additions & 6 deletions src/backends/cuda/active_set_system/global_active_set_manager.cu
Original file line number Diff line number Diff line change
Expand Up @@ -557,7 +557,7 @@ void GlobalActiveSetManager::Impl::update_lambda()
{
auto PT = PTs(idx);
auto mu = min(min(mu_vertices(PT(0)), mu_vertices(PT(1))),
min(mu_vertices(PT(2)), mu_vertices(PT(3))));
min(mu_vertices(PT(2)), mu_vertices(PT(3))));
auto d_grad = PT_d_grad(idx);
auto d = d0(idx), &lambda = PT_lambda(idx), d_shift = 0.0;
auto& cnt = PT_cnt(idx);
Expand Down Expand Up @@ -598,7 +598,7 @@ void GlobalActiveSetManager::Impl::update_lambda()
{
auto EE = EEs(idx);
auto mu = min(min(mu_vertices(EE(0)), mu_vertices(EE(1))),
min(mu_vertices(EE(2)), mu_vertices(EE(3))));
min(mu_vertices(EE(2)), mu_vertices(EE(3))));
auto d_grad = EE_d_grad(idx);
auto d = d0(idx), &lambda = EE_lambda(idx), d_shift = 0.0;
auto& cnt = EE_cnt(idx);
Expand Down Expand Up @@ -851,16 +851,18 @@ muda::BufferView<Float> GlobalActiveSetManager::StiffnessEstimateInfo::mu_vertic

Float GlobalActiveSetManager::StiffnessEstimateInfo::dt() const noexcept
{
return m_impl->dt;
return m_impl->dt_attr->view()[0];
}

void GlobalActiveSetManager::Impl::init(WorldVisitor& world)
{
auto config = world.scene().config();
dt = config.find<Float>("dt")->view()[0];
auto config = world.scene().config();
dt_attr = config.find<Float>("dt");
UIPC_ASSERT(dt_attr, "Scene config must have a 'dt' attribute.");
decay_factor = config.find<Float>("contact/al-ipc/decay_factor")->view()[0];
toi_threshold = config.find<Float>("contact/al-ipc/toi_threshold")->view()[0];
alpha_lower_bound = config.find<Float>("contact/al-ipc/alpha_lower_bound")->view()[0];
alpha_lower_bound =
config.find<Float>("contact/al-ipc/alpha_lower_bound")->view()[0];
energy_enabled = true;
}

Expand Down
11 changes: 6 additions & 5 deletions src/backends/cuda/active_set_system/global_active_set_manager.h
Original file line number Diff line number Diff line change
Expand Up @@ -106,11 +106,12 @@ class GlobalActiveSetManager final : public SimSystem

muda::DeviceBuffer<Vector3> non_penetrate_positions;

Float decay_factor, dt;
Float toi_threshold;
Float alpha_lower_bound;
bool energy_enabled;
bool should_discard_friction_candidates = false;
Float decay_factor;
S<const geometry::AttributeSlot<Float>> dt_attr;
Float toi_threshold;
Float alpha_lower_bound;
bool energy_enabled;
bool should_discard_friction_candidates = false;

Float m_reserve_ratio = 1.5;

Expand Down
6 changes: 4 additions & 2 deletions src/backends/cuda/affine_body/abd_linear_subsystem.cu
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,8 @@ void ABDLinearSubsystem::do_build(DiagLinearSubsystem::BuildInfo& info)
{
m_impl.affine_body_dynamics = require<AffineBodyDynamics>();
m_impl.affine_body_vertex_reporter = require<AffineBodyVertexReporter>();
auto attr = world().scene().config().find<Float>("dt");
m_impl.dt = attr->view()[0];
m_impl.dt_attr = world().scene().config().find<Float>("dt");
UIPC_ASSERT(m_impl.dt_attr, "Scene config must have a 'dt' attribute.");

m_impl.dytopo_effect_receiver = find<ABDDyTopoEffectReceiver>();
}
Expand Down Expand Up @@ -201,6 +201,8 @@ void ABDLinearSubsystem::Impl::_assemble_kinetic_shape(IndexT& hess_offset,
{
using namespace muda;

Float dt = dt_attr->view()[0];

// Collect Kinetic
ABDLinearSubsystem::ComputeGradientHessianInfo this_info{
info.gradient_only(), body_id_to_kinetic_gradient, body_id_to_kinetic_hessian, dt};
Expand Down
12 changes: 6 additions & 6 deletions src/backends/cuda/affine_body/abd_linear_subsystem.h
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,9 @@ class ABDLinearSubsystem final : public DiagLinearSubsystem
private:
friend class ABDLinearSubsystem;
friend class ABDLinearSubsystemReporter;
SizeT m_gradient_count = 0;
SizeT m_hessian_count = 0;
bool m_gradient_only = false;
SizeT m_gradient_count = 0;
SizeT m_hessian_count = 0;
bool m_gradient_only = false;
mutable bool m_gradient_only_checked = false;
};

Expand Down Expand Up @@ -99,8 +99,8 @@ class ABDLinearSubsystem final : public DiagLinearSubsystem
void _assemble_reporters(IndexT& offset, GlobalLinearSystem::DiagInfo& info);
void _assemble_dytopo_effect(IndexT& offset, GlobalLinearSystem::DiagInfo& info);

void accuracy_check(GlobalLinearSystem::AccuracyInfo& info);
void retrieve_solution(GlobalLinearSystem::SolutionInfo& info);
void accuracy_check(GlobalLinearSystem::AccuracyInfo& info);
void retrieve_solution(GlobalLinearSystem::SolutionInfo& info);
Float diag_norm();
Float mass_norm();

Expand Down Expand Up @@ -129,7 +129,7 @@ class ABDLinearSubsystem final : public DiagLinearSubsystem
muda::DeviceBuffer<Float> block_norm;
muda::DeviceVar<Float> reduced_norm;

Float dt = 0.0f; // time step, used in assemble
S<const geometry::AttributeSlot<Float>> dt_attr;
};

private:
Expand Down
Loading
Loading