From a59a19bf1c28d6381ef4f1e1cc28c8b825d813a2 Mon Sep 17 00:00:00 2001 From: Ivan Korotkin Date: Thu, 19 Dec 2024 12:15:44 +0000 Subject: [PATCH 1/3] Update gemfile --- docs/Gemfile.lock | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/docs/Gemfile.lock b/docs/Gemfile.lock index fa4c101..1be0ab0 100644 --- a/docs/Gemfile.lock +++ b/docs/Gemfile.lock @@ -1,22 +1,22 @@ GEM remote: https://rubygems.org/ specs: - addressable (2.8.6) - public_suffix (>= 2.0.2, < 6.0) + addressable (2.8.7) + public_suffix (>= 2.0.2, < 7.0) colorator (1.1.0) - concurrent-ruby (1.2.2) + concurrent-ruby (1.3.4) em-websocket (0.5.3) eventmachine (>= 0.12.9) http_parser.rb (~> 0) eventmachine (1.2.7) - ffi (1.16.3) + ffi (1.17.0) forwardable-extended (2.6.0) - google-protobuf (3.25.1-arm64-darwin) - google-protobuf (3.25.1-x86_64-linux) + google-protobuf (3.25.5-arm64-darwin) + google-protobuf (3.25.5-x86_64-linux) http_parser.rb (0.8.0) - i18n (1.14.1) + i18n (1.14.6) concurrent-ruby (~> 1.0) - jekyll (4.3.3) + jekyll (4.3.4) addressable (~> 2.4) colorator (~> 1.0) em-websocket (~> 0.5) @@ -45,32 +45,32 @@ GEM jekyll-include-cache jekyll-seo-tag (>= 2.0) rake (>= 12.3.1) - kramdown (2.4.0) - rexml + kramdown (2.5.1) + rexml (>= 3.3.9) kramdown-parser-gfm (1.1.0) kramdown (~> 2.0) liquid (4.0.4) - listen (3.8.0) + listen (3.9.0) rb-fsevent (~> 0.10, >= 0.10.3) rb-inotify (~> 0.9, >= 0.9.10) mercenary (0.4.0) pathutil (0.16.2) forwardable-extended (~> 2.6) - public_suffix (5.0.4) - rake (13.1.0) + public_suffix (6.0.1) + rake (13.2.1) rb-fsevent (0.11.2) - rb-inotify (0.10.1) + rb-inotify (0.11.1) ffi (~> 1.0) - rexml (3.2.6) - rouge (4.2.0) + rexml (3.4.0) + rouge (4.5.1) safe_yaml (1.0.5) sass-embedded (1.69.5) google-protobuf (~> 3.23) rake (>= 13.0.0) terminal-table (3.0.2) unicode-display_width (>= 1.1.1, < 3) - unicode-display_width (2.5.0) - webrick (1.8.1) + unicode-display_width (2.6.0) + webrick (1.9.1) PLATFORMS arm64-darwin-23 From 5e099759ccbc6a2f2601f82cb344449e68633ebd Mon Sep 17 00:00:00 2001 From: Ivan Korotkin Date: Fri, 20 Dec 2024 19:52:29 +0000 Subject: [PATCH 2/3] Update version, CHANGELOG, and examples --- docs/CHANGELOG.md | 16 +++++++++++++++- docs/examples.md | 18 ++++++++++++++++++ docs/index.md | 2 +- 3 files changed, 34 insertions(+), 2 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 8233b84..bd84ae2 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -8,11 +8,25 @@ nav_order: 11 All notable user-facing changes to the `dae-cpp` project are documented in this page. -## v2.0.1 +## v2.1.0 New {: .label .label-green } +Version 2.1 allows the user to define the shape (structure) of the Jacobian matrix. I.e., instead of providing analytic Jacobian (or not providing it at all, which is slow if the system is large enough), the user can specify the positions of non-zero elements in the Jacobian. The solver will use automatic differentiation for the specified elements only. This works nearly as fast as analytic Jacobian without requiring the user to differentiate the vector function manually. + +- Added `JacobianMatrixShape` and `VectorFunctionElements` helper classes to define the Jacobian matrix shape and the vector function +- Added `JacobianCompare` class that helps the user to compare the user-defined Jacobian (either defined explicitly or using Jacobian shape) with the one computed automatically from the system RHS +- Added `jacobian_shape` and `jacobian_compare` examples +- Updated `autodiff` to `v1.1.2` +- Updated Eigen (commit from 28th March 2024, see issue [#52](https://github.com/dae-cpp/dae-cpp/issues/52)) +- Updated `googletest` from `v.1.14.x` to `v1.15.2` +- Added integration test that uses Jacobian derived from the given shape +- Added unit tests for all new classes (`JacobianMatrixShape`, `VectorFunctionElements`, `JacobianCompare`) +- Minor solver updates + +## v2.0.1 + - Added [Flame Propagation](https://github.com/dae-cpp/dae-cpp/blob/master/examples/flame_propagation/flame_propagation.cpp) example (stiff equation) - Added `daecpp::dual_type` for automatic differentiation of the vector function (used in [Flame Propagation](https://github.com/dae-cpp/dae-cpp/blob/master/examples/flame_propagation/flame_propagation.cpp) example) - Added integration test based on the "Flame Propagation" example diff --git a/docs/examples.md b/docs/examples.md index 032bbec..81162ee 100644 --- a/docs/examples.md +++ b/docs/examples.md @@ -42,3 +42,21 @@ When you light a match, the ball of flame grows rapidly until it reaches a criti Then it remains at that size because the amount of oxygen being consumed by the combustion in the interior of the ball balances the amount of oxigen available through the surface. This example solves **stiff** equation of flame propagation for the scalar variable $$y(t)$$ which represents the radius of the ball. + +## Jacobian matrix shape + +Jacobian matrix shape source file example: [jacobian_shape.cpp](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_shape/jacobian_shape.cpp) +{: .fs-5 .fw-400 } + +This example demonstrates how to define the shape (structure) of the Jacobian matrix. +I.e., instead of providing the full analytic Jacobian (or not providing it at all, which is slow if the system is big), the user can specify the positions of non-zero elements in the Jacobian. +The solver will use automatic differentiation for the specified elements only. +This works (nearly) as fast as the analytic Jacobian without requiring the user to differentiate the vector function manually. + +## Checking user-defined Jacobian + +Jacobian matrix checking source file example: [jacobian_compare.cpp](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_compare/jacobian_compare.cpp) +{: .fs-5 .fw-400 } + +In this example, we do not solve any DAEs. Instead, we define a simple vector function, the corresponding Jacobian matrix, and then we use a built-in helper class `JacobianCompare` to compare our manually derived Jacobian with the one computed algorithmically from the vector function. +We will make a few mistakes in the analytic Jacobian on purpose to see what information `JacobianCompare` can provide. diff --git a/docs/index.md b/docs/index.md index 22d9b35..ede26d2 100644 --- a/docs/index.md +++ b/docs/index.md @@ -4,7 +4,7 @@ layout: home nav_order: 1 --- -![version](https://img.shields.io/badge/version-2.0.1-blue) +![version](https://img.shields.io/badge/version-2.1.0-blue)

dae-cpp logo From 7dab4f19c6d1ae7c192d72f4356bac961632f468 Mon Sep 17 00:00:00 2001 From: Ivan Korotkin Date: Sat, 21 Dec 2024 22:09:22 +0000 Subject: [PATCH 3/3] Documented new classes in version 2.1 --- docs/CHANGELOG.md | 8 ++-- docs/examples.md | 10 ++--- docs/index.md | 1 + docs/jacobian-matrix.md | 87 ++++++++++++++++++++++++++++++++++++++++- docs/vector-function.md | 39 ++++++++++++++++++ 5 files changed, 134 insertions(+), 11 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index bd84ae2..971b5ac 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -15,14 +15,14 @@ New Version 2.1 allows the user to define the shape (structure) of the Jacobian matrix. I.e., instead of providing analytic Jacobian (or not providing it at all, which is slow if the system is large enough), the user can specify the positions of non-zero elements in the Jacobian. The solver will use automatic differentiation for the specified elements only. This works nearly as fast as analytic Jacobian without requiring the user to differentiate the vector function manually. -- Added `JacobianMatrixShape` and `VectorFunctionElements` helper classes to define the Jacobian matrix shape and the vector function -- Added `JacobianCompare` class that helps the user to compare the user-defined Jacobian (either defined explicitly or using Jacobian shape) with the one computed automatically from the system RHS -- Added `jacobian_shape` and `jacobian_compare` examples +- Added [`daecpp::JacobianMatrixShape`](jacobian-matrix.html#jacobian-matrix-shape) and [`daecpp::VectorFunctionElements`](vector-function.html#element-by-element-vector-function-to-define-the-jacobian-shape) helper classes to define the Jacobian matrix shape and the vector function +- Added [`daecpp::JacobianCompare`](jacobian-matrix.html#jacobian-matrix-check) class that helps the user to compare the user-defined Jacobian (either defined explicitly or using Jacobian shape) with the one computed automatically from the system RHS +- Added [Jacobian shape](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_shape/jacobian_shape.cpp) and [Jacobian compare](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_compare/jacobian_compare.cpp) examples - Updated `autodiff` to `v1.1.2` - Updated Eigen (commit from 28th March 2024, see issue [#52](https://github.com/dae-cpp/dae-cpp/issues/52)) - Updated `googletest` from `v.1.14.x` to `v1.15.2` - Added integration test that uses Jacobian derived from the given shape -- Added unit tests for all new classes (`JacobianMatrixShape`, `VectorFunctionElements`, `JacobianCompare`) +- Added unit tests for all new classes (`daecpp::JacobianMatrixShape`, `daecpp::VectorFunctionElements`, `daecpp::JacobianCompare`) - Minor solver updates ## v2.0.1 diff --git a/docs/examples.md b/docs/examples.md index 81162ee..ac60054 100644 --- a/docs/examples.md +++ b/docs/examples.md @@ -48,15 +48,15 @@ This example solves **stiff** equation of flame propagation for the scalar varia Jacobian matrix shape source file example: [jacobian_shape.cpp](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_shape/jacobian_shape.cpp) {: .fs-5 .fw-400 } -This example demonstrates how to define the shape (structure) of the Jacobian matrix. +This example demonstrates how to define the [shape](jacobian-matrix.html#jacobian-matrix-shape) (structure) of the Jacobian matrix. I.e., instead of providing the full analytic Jacobian (or not providing it at all, which is slow if the system is big), the user can specify the positions of non-zero elements in the Jacobian. The solver will use automatic differentiation for the specified elements only. This works (nearly) as fast as the analytic Jacobian without requiring the user to differentiate the vector function manually. -## Checking user-defined Jacobian +## Debugging the user-defined Jacobian -Jacobian matrix checking source file example: [jacobian_compare.cpp](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_compare/jacobian_compare.cpp) +Debugging the user-defined Jacobian example: [jacobian_compare.cpp](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_compare/jacobian_compare.cpp) {: .fs-5 .fw-400 } -In this example, we do not solve any DAEs. Instead, we define a simple vector function, the corresponding Jacobian matrix, and then we use a built-in helper class `JacobianCompare` to compare our manually derived Jacobian with the one computed algorithmically from the vector function. -We will make a few mistakes in the analytic Jacobian on purpose to see what information `JacobianCompare` can provide. +In this example, we do not solve any DAEs. Instead, we define a simple vector function, the corresponding Jacobian matrix, and then we use a built-in helper class [`daecpp::JacobianCompare`](jacobian-matrix.html#jacobian-matrix-check) to compare our manually derived Jacobian with the one computed algorithmically from the vector function. +We will make a few mistakes in the analytic Jacobian on purpose to see what information [`daecpp::JacobianCompare`](jacobian-matrix.html#jacobian-matrix-check) can provide. diff --git a/docs/index.md b/docs/index.md index ede26d2..7a41d09 100644 --- a/docs/index.md +++ b/docs/index.md @@ -30,6 +30,7 @@ Eigen's sparse solver performs two steps: factorization (decomposition) of the J - Header only, no pre-compilation required. - Uses [automatic](https://en.wikipedia.org/wiki/Automatic_differentiation) (algorithmic, exact) differentiation (using [autodiff](https://autodiff.github.io/)) to compute the Jacobian matrix, if it is not provided by the user. +- The user can provide analytically derived Jacobian or the Jacobian matrix shape (positions of non-zero elements) to significantly speed up the computation for big systems. - Fourth-order implicit BDF time integrator that preserves accuracy even when the time step rapidly changes. - A very flexible and customizable variable time stepping algorithm based on the solution stability and variability. - Mass matrix can be non-static (can depend on time) and it can be singular. diff --git a/docs/jacobian-matrix.md b/docs/jacobian-matrix.md index c73b903..852b9ba 100644 --- a/docs/jacobian-matrix.md +++ b/docs/jacobian-matrix.md @@ -136,5 +136,88 @@ automaticJacobian(J, x, t); std::cout << J.dense(N) << '\n'; ``` -{: .highlight } -*User-defined* Jacobian vs *automatic* Jacobian comparison functionality (element by element) will be added in future versions of the solver. +## Jacobian matrix shape + +If defining the Jacobian matrix manually as it is described above is not a feasible task (e.g., due to a very complex non-linear RHS), the solver allows the user to specify only the positions of non-zero elements of the Jacobian matrix (i.e., the *Jacobian matrix shape*). In this case, all the derivatives will be calculated automatically with a very small computation time penalty (compared to the manually derived analytic Jacobian). + +To define the Jacobian matrix shape, the user needs to define the vector function (RHS) of the system in an element-by-element (equation-by-equation) way using a custom class derived from [`daecpp::VectorFunctionElements`](vector-function.html#element-by-element-vector-function-to-define-the-jacobian-shape). +After that, use another class derived from `daecpp::JacobianMatrixShape` to specify the positions of each non-zero element in the Jacobian matrix. + +Class `daecpp::JacobianMatrixShape` contains the following helper methods: + +- `void add_element(const daecpp::int_type ind_i, const daecpp::int_type ind_j)` to define a single non-zero element *(i, j)* +- `void add_element(const daecpp::int_type ind_i, const std::vector &ind_j)` to define an entire row *i* of non-zero elements using a vector of indices *j* +- `void clear()` to clear the array of non-zero elements +- `void reserve(const daecpp::int_type N_elements)` to reserve memory for `N_elements` non-zero elements of the Jacobian (which is recommended) + +The constructor of `daecpp::JacobianMatrixShape` takes the vector function object derived from [`daecpp::VectorFunctionElements`](vector-function.html#element-by-element-vector-function-to-define-the-jacobian-shape). + +Here is a simple example of defining 3 non-zero elements of the Jacobian: + +```cpp +struct UserDefinedJacobianShape : daecpp::JacobianMatrixShape +{ + explicit UserDefinedJacobianShape(UserDefinedRHS rhs) : daecpp::JacobianMatrixShape(rhs) + { + reserve(3); // Reserves memory for 3 non-zero elements + add_element(0, 1); // Adds element (0, 1), e.g., row 0, column 1 + add_element(1, {0, 1}); // Adds elements (1, 0) and (1, 1) + } +}; +``` + +{: .important } +Note that the RHS should be defined using [`daecpp::VectorFunctionElements`](vector-function.html#element-by-element-vector-function-to-define-the-jacobian-shape), so the solver can call specific elements of the vector function individually to perform automatic differentiation element-by-element efficiently. + +For more details, refer to the [Jacobian shape](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_shape/jacobian_shape.cpp) example. + +## Jacobian matrix check + +User-defined Jacobians defined either analytically or using Jacobian shape can be conveniently compared to fully automatic Jacobian for debug purposes. +This can be achieved by using `daecpp::JacobianCompare` helper class. +The constructor of this class takes the user-defined Jacobian (defined analytically or from the Jacobian shape) and the vector function (RHS). +Then the object of the `daecpp::JacobianCompare` class can be called as a functor with the given state `x` and, optionally, time `t`. + +Consider the following example: + +```cpp +// Fill vectors x at which the Jacobian matrix will be tested +daecpp::state_vector x0 = {1.0, 2.0, 3.0, 4.0}; +daecpp::state_vector x1 = {-0.5, -0.1, 0.5, 0.1}; + +// Check user-defined Jacobian at `x0` and `x1` and times t = 1 and 2 +{ + auto jac_comparison = daecpp::JacobianCompare(UserJacobian(), UserRHS()); + + auto N_errors_1 = jac_comparison(x0, 1.0); + auto N_errors_2 = jac_comparison(x1, 2.0); +} +``` + +The functor `jac_comparison` returns the number of errors found in the Jacobian. +It also prints on screen a summary of all inconsistencies found in the user-defined Jacobian compared to the fully automatic one. +An example of such output is given below: + +```txt +Jacobian matrix comparison summary at time t = 1: +-- Found 4 difference(s) compared to the automatic (reference) Jacobian: +---------------------------------------------------------------------------------------- + Row | Column | Reference value | User-defined value | Absolute error +------------+------------+--------------------+--------------------+-------------------- + 2 | 0 | 0.540302 | 0 | -0.540302 + 2 | 1 | 0 | 0.540302 | 0.540302 + 0 | 3 | 0.5 | 1 | 0.5 + 3 | 3 | 1 | 0 | -1 +---------------------------------------------------------------------------------------- +``` + +Here we can see that in this example: + +- element (2, 1) should be (2, 0) +- element (0, 3) is incorrect (should be 0.5, not 1) +- element (3, 3) is not defined at all + +{: .note } +Jacobian comparisons are element-by-element and hence can be slow. These comparisons are for the debug purposes only and should be removed from the production runs. + +For more details, refer to the [Jacobian compare](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_compare/jacobian_compare.cpp) example. diff --git a/docs/vector-function.md b/docs/vector-function.md index 1357282..4e14d61 100644 --- a/docs/vector-function.md +++ b/docs/vector-function.md @@ -80,3 +80,42 @@ struct UserDefinedVectorFunction } }; ``` + +## Element-by-element vector function to define the Jacobian shape + +For relatively small systems, there is no need to provide the Jacobian matrix (neither its shape), as it will be computed fast enough automatically. +If the user decides to provide manually derived Jacobian, it can be done using [Jacobian Matrix](jacobian-matrix.html) class. +However, if the user provides the [shape of the Jacobian matrix](jacobian-matrix.html#jacobian-matrix-shape), the vector function must be defined element-by-element, so the solver can call specific elements of the vector function individually to perform automatic differentiation. + +The vector function from the [Examples](#examples) section above can be defined element-by-element by inheriting the `daecpp::VectorFunctionElements` class: + +```cpp +struct UserDefinedVectorFunction : daecpp::VectorFunctionElements +{ + daecpp::dual_type equations(const daecpp::state_type &x, const double t, const daecpp::int_type i) const + { + if (i == 0) + return x[2] + 1.0; // z + 1 + else if (i == 1) + return x[0] * x[0] + x[1]; // x*x + y + else if (i == 2) + return 2.0 * t // 2*t + else + { + // This should never happen + ERROR("Index i in UserDefinedVectorFunction is out of boundaries"); + } + } +}; +``` + +Unlike `daecpp::VectorFunction`, here we define each element of the vector function individually, depending on the index `i`. +This obviously comes with some computational time penalty compared to the case where we define the entire vector function as an array. +However, defining the RHS element-by-element like in the example above will allow us to define only the positions of non-zero elements of the Jacobian without manually differentiating the entire vector function (which in some cases can be difficult and error-prone). + +{: .note } +The class `VectorFunctionElements` is abstract, where `dual_type equations(const state_type &x, const double t, const int_type i) const` function must be implemented in the derived class. + +After defining the vector function this way, define the [Jacobian matrix shape](jacobian-matrix.html#jacobian-matrix-shape). + +For more detailed example, refer to the [Jacobian shape](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_shape/jacobian_shape.cpp) example.