diff --git a/runtime/cudaq/base_integrator.h b/runtime/cudaq/base_integrator.h index e3196bd52d..0d63acf4d5 100644 --- a/runtime/cudaq/base_integrator.h +++ b/runtime/cudaq/base_integrator.h @@ -31,13 +31,27 @@ class BaseIntegrator { virtual void post_init() = 0; public: + /// @brief Default constructor + BaseIntegrator() = default; + + /// @brief Constructor to initialize the integrator with a state and time + /// stepper. + /// @param initial_state Initial quantum state. + /// @param t0 Initial time. + /// @param stepper Time stepper instance. + BaseIntegrator(const TState &initial_state, double t0, + std::shared_ptr> stepper) + : state(initial_state), t(t0), stepper(std::move(stepper)) {} + virtual ~BaseIntegrator() = default; + /// @brief Set the initial state and time void set_state(const TState &initial_state, double t0 = 0.0) { state = initial_state; t = t0; } + /// @brief Set the system parameters (dimensions, schedule, and operators) void set_system( const std::map &dimensions, std::shared_ptr schedule, std::shared_ptr hamiltonian, @@ -48,8 +62,10 @@ class BaseIntegrator { this->collapse_operators = collapse_operators; } - virtual void integrate(double t) = 0; + /// @brief Perform integration to the target time. + virtual void integrate(double target_time) = 0; + /// @brief Get the current time and state. std::pair get_state() const { return {t, state}; } }; } // namespace cudaq diff --git a/runtime/cudaq/dynamics/CMakeLists.txt b/runtime/cudaq/dynamics/CMakeLists.txt index 636d8b4096..63129f246f 100644 --- a/runtime/cudaq/dynamics/CMakeLists.txt +++ b/runtime/cudaq/dynamics/CMakeLists.txt @@ -11,7 +11,18 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-ctad-maybe-unsupported") set(INTERFACE_POSITION_INDEPENDENT_CODE ON) set(CUDAQ_OPS_SRC - scalar_operators.cpp elementary_operators.cpp product_operators.cpp operator_sum.cpp schedule.cpp definition.cpp helpers.cpp rydberg_hamiltonian.cpp cudm_helpers.cpp cudm_state.cpp cudm_time_stepper.cpp + scalar_operators.cpp + elementary_operators.cpp + product_operators.cpp + operator_sum.cpp + schedule.cpp + definition.cpp + helpers.cpp + rydberg_hamiltonian.cpp + cudm_helpers.cpp + cudm_state.cpp + cudm_time_stepper.cpp + runge_kutta_integrator.cpp ) set(CUQUANTUM_INSTALL_PREFIX "/usr/local/lib/python3.10/dist-packages/cuquantum") diff --git a/runtime/cudaq/dynamics/runge_kutta_integrator.cpp b/runtime/cudaq/dynamics/runge_kutta_integrator.cpp new file mode 100644 index 0000000000..62633bc172 --- /dev/null +++ b/runtime/cudaq/dynamics/runge_kutta_integrator.cpp @@ -0,0 +1,65 @@ +/******************************************************************************* + * Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#include "cudaq/runge_kutta_integrator.h" +#include + +using namespace cudaq; + +namespace cudaq { +void runge_kutta_integrator::integrate(double target_time) { + if (!stepper) { + throw std::runtime_error("Time stepper is not initialized."); + } + + double dt = integrator_options["dt"]; + if (dt <= 0) { + throw std::invalid_argument("Invalid time step size for integration."); + } + + auto handle = state.get_handle(); + auto hilbertSpaceDims = state.get_hilbert_space_dims(); + + while (t < target_time) { + double step_size = std::min(dt, target_time - 1); + + std::cout << "Runge-Kutta step at time " << t << " with step size: " << step_size << std::endl; + + // Empty vectors of same size as state.get_raw_data() + std::vector> zero_state(state.get_raw_data().size(), {0.0, 0.0}); + + cudm_state k1(handle, zero_state, hilbertSpaceDims); + cudm_state k2(handle, zero_state, hilbertSpaceDims); + cudm_state k3(handle, zero_state, hilbertSpaceDims); + cudm_state k4(handle, zero_state, hilbertSpaceDims); + + if (substeps_ == 1) { + // Euler method (1st order) + k1 = stepper->compute(state, t, step_size); + state = k1; + } else if (substeps_ == 2) { + // Midpoint method (2nd order) + k1 = stepper->compute(state, t, step_size / 2.0); + k2 = stepper->compute(k1, t + step_size / 2.0, step_size); + state = (k1 + k2) * 0.5; + } else if (substeps_ == 4) { + // Runge-Kutta method (4th order) + k1 = stepper->compute(state, t, step_size / 2.0); + k2 = stepper->compute(k1, t + step_size / 2.0, step_size / 2.0); + k3 = stepper->compute(k2, t + step_size / 2.0, step_size); + k4 = stepper->compute(k3, t + step_size, step_size); + state = (k1 + k2 * 2.0 + k3 * 2.0 + k4) * (1.0 / 6.0); + } + + // Update time + t += step_size; + } + + std::cout << "Integration complete. Final time: " << t << std::endl; +} +} diff --git a/runtime/cudaq/runge_kutta_integrator.h b/runtime/cudaq/runge_kutta_integrator.h index 9914258386..0b98bb4d86 100644 --- a/runtime/cudaq/runge_kutta_integrator.h +++ b/runtime/cudaq/runge_kutta_integrator.h @@ -9,39 +9,38 @@ #pragma once #include "base_integrator.h" -#include "runge_kutta_time_stepper.h" +#include "cudaq/cudm_state.h" +#include "cudaq/cudm_time_stepper.h" #include namespace cudaq { -template -class RungeKuttaIntegrator : public BaseIntegrator { +class runge_kutta_integrator : public BaseIntegrator { public: - using DerivativeFunction = std::function; - - explicit RungeKuttaIntegrator(DerivativeFunction f) - : stepper(std::make_shared>(f)) {} - - // Initializes the integrator - void post_init() override { - if (!this->stepper) { - throw std::runtime_error("Time stepper is not set"); + /// @brief Constructor to initialize the Runge-Kutta integrator + /// @param initial_state Initial quantum state. + /// @param t0 Initial time. + /// @param stepper Time stepper instance. + /// @param substeps Number of Runge-Kutta substeps (must be 1, 2, or 4) + runge_kutta_integrator(const cudm_state &initial_state, double t0, + std::shared_ptr stepper, + int substeps = 4) + : BaseIntegrator(initial_state, t0, stepper), substeps_(substeps) { + if (substeps_ != 1 && substeps_ != 2 && substeps_ != 4) { + throw std::invalid_argument("Runge-Kutta substeps must be 1, 2, or 4."); } + post_init(); } - // Advances the system's state from current time to `t` - void integrate(double target_t) override { - if (!this->schedule || !this->hamiltonian) { - throw std::runtime_error("System is not properly set!"); - } + /// @brief Perform Runge-Kutta integration until the target time. + /// @param target_time The final time to integrate to. + void integrate(double t) override; - while (this->t < target_t) { - stepper->compute(this->state, this->t); - // Time step size - this->t += 0.01; - } - } +protected: + /// @brief Any post-initialization setup + void post_init() override {} private: - std::shared_ptr> stepper; + // Number of substeps in RK integration (1, 2, or 4) + int substeps_; }; -} // namespace cudaq \ No newline at end of file +} // namespace cudaq diff --git a/runtime/cudaq/runge_kutta_time_stepper.h b/runtime/cudaq/runge_kutta_time_stepper.h deleted file mode 100644 index 1dcd1f69cc..0000000000 --- a/runtime/cudaq/runge_kutta_time_stepper.h +++ /dev/null @@ -1,33 +0,0 @@ -/******************************************************************************* - * Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates. * - * All rights reserved. * - * * - * This source code and the accompanying materials are made available under * - * the terms of the Apache License 2.0 which accompanies this distribution. * - ******************************************************************************/ - -#include "base_time_stepper.h" -#include - -namespace cudaq { -template -class RungeKuttaTimeStepper : public BaseTimeStepper { -public: - using DerivativeFunction = std::function; - - RungeKuttaTimeStepper(DerivativeFunction f) : derivativeFunc(f) {} - - void compute(TState &state, double t, double dt = 0.01) override { - // 4th order Runge-Kutta method - TState k1 = derivativeFunc(state, t); - TState k2 = derivativeFunc(state + (dt / 2.0) * k1, t + dt / 2.0); - TState k3 = derivativeFunc(state + (dt / 2.0) * k2, t + dt / 2.0); - TState k4 = derivativeFunc(state + dt * k3, t + dt); - - state = state + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4); - } - -private: - DerivativeFunction derivativeFunc; -}; -} // namespace cudaq \ No newline at end of file diff --git a/unittests/CMakeLists.txt b/unittests/CMakeLists.txt index 0d99fd5bb2..c422dfa8ee 100644 --- a/unittests/CMakeLists.txt +++ b/unittests/CMakeLists.txt @@ -49,7 +49,6 @@ set(CUDAQ_RUNTIME_TEST_SOURCES dynamics/elementary_ops_simple.cpp dynamics/elementary_ops_arithmetic.cpp dynamics/product_operators_arithmetic.cpp - dynamics/test_runge_kutta_time_stepper.cpp dynamics/test_runge_kutta_integrator.cpp dynamics/test_helpers.cpp dynamics/rydberg_hamiltonian.cpp diff --git a/unittests/dynamics/runge_kutta_test_helpers.h b/unittests/dynamics/runge_kutta_test_helpers.h deleted file mode 100644 index 4f93ffa242..0000000000 --- a/unittests/dynamics/runge_kutta_test_helpers.h +++ /dev/null @@ -1,24 +0,0 @@ -/****************************************************************-*- C++ -*-**** - * Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates. * - * All rights reserved. * - * * - * This source code and the accompanying materials are made available under * - * the terms of the Apache License 2.0 which accompanies this distribution. * - ******************************************************************************/ - -#pragma once - -#include - -// A simple state type -using TestState = double; - -// Simple derivative function: dx/dt = -x (exponential decay) -inline TestState simple_derivative(const TestState &state, double t) { - return -state; -} - -// A complex function: dx/dt = sin(t) -inline TestState sine_derivative(const TestState &state, double t) { - return std::sin(t); -} diff --git a/unittests/dynamics/test_runge_kutta_integrator.cpp b/unittests/dynamics/test_runge_kutta_integrator.cpp index c75a7c8d6d..407c211210 100644 --- a/unittests/dynamics/test_runge_kutta_integrator.cpp +++ b/unittests/dynamics/test_runge_kutta_integrator.cpp @@ -7,7 +7,7 @@ // ******************************************************************************/ #include "cudaq/runge_kutta_integrator.h" -#include "runge_kutta_test_helpers.h" +#include "test_mocks.h" #include #include #include @@ -17,116 +17,43 @@ using namespace cudaq; // Test fixture class class RungeKuttaIntegratorTest : public ::testing::Test { protected: - RungeKuttaIntegrator *integrator; - std::shared_ptr schedule; - std::shared_ptr hamiltonian; + cudensitymatHandle_t handle_; + cudensitymatOperator_t liouvillian_; + std::shared_ptr time_stepper_; + std::unique_ptr integrator_; + std::unique_ptr state_; void SetUp() override { - integrator = new RungeKuttaIntegrator(simple_derivative); - // Initial state and time - integrator->set_state(1.0, 0.0); + // Create library handle + HANDLE_CUDM_ERROR(cudensitymatCreate(&handle_)); - // A simple step sequence for the schedule - std::vector> steps = {0.1, 0.2, 0.3, 0.4, 0.5}; + // Create a mock Liouvillian + liouvillian_ = mock_liouvillian(handle_); - // Dummy parameters - std::vector parameters = {"param1"}; + // Initialize the time stepper + time_stepper_ = std::make_shared(handle_, liouvillian_); - // A simple parameter function - auto value_function = [](const std::string ¶m, - const std::complex &step) { return step; }; + // Create initial state + state_ = std::make_unique(handle_, mock_initial_state_data(), + mock_hilbert_space_dims()); - // A valid schedule instance - schedule = std::make_shared(steps, parameters, value_function); + double t0 = 0.0; + // Initialize the integrator (using substeps = 2, for mid-point rule) + integrator_ = + std::make_unique(*state_, t0, time_stepper_, 2); - // A simple hamiltonian as an operator_sum - hamiltonian = std::make_shared(); - *hamiltonian += 0.5 * elementary_operator::identity(0); - *hamiltonian += 0.5 * elementary_operator::number(0); - - // System with valid components - integrator->set_system({{0, 2}}, schedule, hamiltonian); - } - - void TearDown() override { delete integrator; } -}; - -// Basic integration -TEST_F(RungeKuttaIntegratorTest, BasicIntegration) { - integrator->integrate(1.0); - - // Expected result: x(t) = e^(-t) - double expected = std::exp(-1.0); - - EXPECT_NEAR(integrator->get_state().second, expected, 1e-3) - << "Basic Runge-Kutta integration failed!"; -} - -// Time evolution -TEST_F(RungeKuttaIntegratorTest, TimeEvolution) { - integrator->integrate(2.0); - - double expected = 2.0; - - EXPECT_NEAR(integrator->get_state().first, expected, 1e-5) - << "Integrator did not correctly update time!"; -} - -// Large step size -TEST_F(RungeKuttaIntegratorTest, LargeStepSize) { - integrator->integrate(5.0); - - double expected = std::exp(-5.0); - - EXPECT_NEAR(integrator->get_state().second, expected, 1e-1) - << "Runge-Kutta integration failed for large step size!!"; -} - -// // Integrating Sine function -// TEST_F(RungeKuttaIntegratorTest, SineFunction) { -// integrator = new RungeKuttaIntegrator(sine_derivative); -// integrator->set_state(1.0, 0.0); -// integrator->set_system({{0, 2}}, schedule, hamiltonian); - -// integrator->integrate(M_PI / 2); - -// double expected = std::cos(M_PI / 2); - -// EXPECT_NEAR(integrator->get_state().second, expected, 1e-3) << -// "Runge-Kutta integration for sine function failed!"; -// } - -// Small step size -TEST_F(RungeKuttaIntegratorTest, SmallStepIntegration) { - integrator->set_state(1.0, 0.0); - integrator->set_system({{0, 2}}, schedule, hamiltonian); - - double step_size = 0.001; - while (integrator->get_state().first < 1.0) { - integrator->integrate(integrator->get_state().first + step_size); + ASSERT_TRUE(state_->is_initialized()); } - double expected = std::exp(-1.0); - - EXPECT_NEAR(integrator->get_state().second, expected, 5e-4) - << "Runge-Kutta integration for small step size failed!"; -} - -// Large step size -TEST_F(RungeKuttaIntegratorTest, LargeStepIntegration) { - integrator->set_state(1.0, 0.0); - integrator->set_system({{0, 2}}, schedule, hamiltonian); - - double step_size = 0.5; - double t = 0.0; - double target_t = 1.0; - while (t < target_t) { - integrator->integrate(std::min(t + step_size, target_t)); - t += step_size; + void TearDown() override { + // Clean up resources + HANDLE_CUDM_ERROR(cudensitymatDestroyOperator(liouvillian_)); + HANDLE_CUDM_ERROR(cudensitymatDestroy(handle_)); } +}; - double expected = std::exp(-1.0); - - EXPECT_NEAR(integrator->get_state().second, expected, 1e-2) - << "Runge-Kutta integration for large step size failed!"; +// Test Initialization +TEST_F(RungeKuttaIntegratorTest, Initialization) { + ASSERT_NE(integrator_, nullptr); + // ASSERT_TRUE(state_->is_initialized()); } diff --git a/unittests/dynamics/test_runge_kutta_time_stepper.cpp b/unittests/dynamics/test_runge_kutta_time_stepper.cpp deleted file mode 100644 index 4c4c7b7588..0000000000 --- a/unittests/dynamics/test_runge_kutta_time_stepper.cpp +++ /dev/null @@ -1,152 +0,0 @@ -/******************************************************************************* - * Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates. * - * All rights reserved. * - * * - * This source code and the accompanying materials are made available under * - * the terms of the Apache License 2.0 which accompanies this distribution. * - ******************************************************************************/ - -#include "cudaq/runge_kutta_time_stepper.h" -#include "runge_kutta_test_helpers.h" -#include -#include -#include - -// Test fixture class -class RungeKuttaTimeStepperTest : public ::testing::Test { -protected: - std::shared_ptr> stepper; - - void SetUp() override { - stepper = std::make_shared>( - simple_derivative); - } -}; - -// Single step integration -TEST_F(RungeKuttaTimeStepperTest, SingleStep) { - // Initial values - double state = 1.0; - double t = 0.0; - double dt = 0.1; - - stepper->compute(state, t, dt); - - // Expected result using analytical solution: x(t) = e^(-t) - double expected = std::exp(-dt); - - EXPECT_NEAR(state, expected, 1e-3) - << "Single step Runge-Kutta integration failed!"; -} - -// Multiple step integration -TEST_F(RungeKuttaTimeStepperTest, MultipleSteps) { - // Initial values - double state = 1.0; - double t = 0.0; - double dt = 0.1; - int steps = 10; - - for (int i = 0; i < steps; i++) { - stepper->compute(state, t, dt); - } - - // Expected result: x(t) = e^(-t) - double expected = std::exp(-1.0); - - EXPECT_NEAR(state, expected, 1e-2) - << "Multiple step Runge-Kutta integration failed!"; -} - -// Convergence to Analytical Solution -TEST_F(RungeKuttaTimeStepperTest, Convergence) { - // Initial values - double state = 1.0; - double t = 0.0; - double dt = 0.01; - int steps = 100; - - for (int i = 0; i < steps; i++) { - stepper->compute(state, t, dt); - } - - double expected = std::exp(-1.0); - - EXPECT_NEAR(state, expected, 1e-3) - << "Runge-Kutta integration does not converge!"; -} - -// // Integrating Sine function -// TEST_F(RungeKuttaTimeStepperTest, SineFunction) { -// auto sine_stepper = -// std::make_shared>(sine_derivative); - -// // Initial values -// double state = 0.0; -// double t = 0.0; -// double dt = 0.1; -// int steps = 10; - -// for (int i = 0; i < steps; i++) { -// sine_stepper->compute(state, t, dt); -// } - -// // Expected integral of sin(t) over [0, 1] is 1 - cos(1) -// double expected = 1 - std::cos(1); - -// EXPECT_NEAR(state, expected, 1e-2) << "Runge-Kutta integration for sine -// function failed!"; -// } - -// Handling small steps sizes -TEST_F(RungeKuttaTimeStepperTest, SmallStepSize) { - // Initial values - double state = 1.0; - double t = 0.0; - double dt = 1e-5; - int steps = 100000; - - for (int i = 0; i < steps; i++) { - stepper->compute(state, t, dt); - } - - double expected = std::exp(-1.0); - - EXPECT_NEAR(state, expected, 1e-3) - << "Runge-Kutta fails with small step sizes!"; -} - -// Handling large steps sizes -TEST_F(RungeKuttaTimeStepperTest, LargeStepSize) { - // Initial values - double state = 1.0; - double t = 0.0; - double dt = 1.0; - - stepper->compute(state, t, dt); - - double expected = std::exp(-1.0); - - EXPECT_NEAR(state, expected, 1e-1) - << "Runge-Kutta is unstable with large step sizes!"; -} - -// Constant derivative (dx/dt = 0) -TEST_F(RungeKuttaTimeStepperTest, ConstantFunction) { - auto constant_stepper = - std::make_shared>( - [](const TestState &state, double t) { return 0.0; }); - - // Initial values - double state = 5.0; - double t = 0.0; - double dt = 0.1; - int steps = 10; - - for (int i = 0; i < steps; i++) { - constant_stepper->compute(state, t, dt); - } - - EXPECT_NEAR(state, 5.0, 1e-6) - << "Runge-Kutta should not change a constant function!"; -}