You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: docs/CHANGELOG.md
+15-1
Original file line number
Diff line number
Diff line change
@@ -8,11 +8,25 @@ nav_order: 11
8
8
9
9
All notable user-facing changes to the `dae-cpp` project are documented in this page.
10
10
11
-
## v2.0.1
11
+
## v2.1.0
12
12
13
13
New
14
14
{: .label .label-green }
15
15
16
+
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.
17
+
18
+
- 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
19
+
- 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
20
+
- 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
21
+
- Updated `autodiff` to `v1.1.2`
22
+
- Updated Eigen (commit from 28th March 2024, see issue [#52](https://github.com/dae-cpp/dae-cpp/issues/52))
23
+
- Updated `googletest` from `v.1.14.x` to `v1.15.2`
24
+
- Added integration test that uses Jacobian derived from the given shape
25
+
- Added unit tests for all new classes (`daecpp::JacobianMatrixShape`, `daecpp::VectorFunctionElements`, `daecpp::JacobianCompare`)
26
+
- Minor solver updates
27
+
28
+
## v2.0.1
29
+
16
30
- Added [Flame Propagation](https://github.com/dae-cpp/dae-cpp/blob/master/examples/flame_propagation/flame_propagation.cpp) example (stiff equation)
17
31
- 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)
18
32
- Added integration test based on the "Flame Propagation" example
This example demonstrates how to define the [shape](jacobian-matrix.html#jacobian-matrix-shape) (structure) of the Jacobian matrix.
52
+
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.
53
+
The solver will use automatic differentiation for the specified elements only.
54
+
This works (nearly) as fast as the analytic Jacobian without requiring the user to differentiate the vector function manually.
55
+
56
+
## Debugging the user-defined Jacobian
57
+
58
+
Debugging the user-defined Jacobian example: [jacobian_compare.cpp](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_compare/jacobian_compare.cpp)
59
+
{: .fs-5 .fw-400 }
60
+
61
+
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.
62
+
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.
@@ -30,6 +30,7 @@ Eigen's sparse solver performs two steps: factorization (decomposition) of the J
30
30
31
31
- Header only, no pre-compilation required.
32
32
- 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.
33
+
- 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.
33
34
- Fourth-order implicit BDF time integrator that preserves accuracy even when the time step rapidly changes.
34
35
- A very flexible and customizable variable time stepping algorithm based on the solution stability and variability.
35
36
- Mass matrix can be non-static (can depend on time) and it can be singular.
Copy file name to clipboardExpand all lines: docs/jacobian-matrix.md
+85-2
Original file line number
Diff line number
Diff line change
@@ -136,5 +136,88 @@ automaticJacobian(J, x, t);
136
136
std::cout << J.dense(N) << '\n';
137
137
```
138
138
139
-
{: .highlight }
140
-
*User-defined* Jacobian vs *automatic* Jacobian comparison functionality (element by element) will be added in future versions of the solver.
139
+
## Jacobian matrix shape
140
+
141
+
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).
142
+
143
+
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).
144
+
After that, use another class derived from `daecpp::JacobianMatrixShape` to specify the positions of each non-zero element in the Jacobian matrix.
145
+
146
+
Class `daecpp::JacobianMatrixShape` contains the following helper methods:
147
+
148
+
- `void add_element(const daecpp::int_type ind_i, const daecpp::int_type ind_j)` to define a single non-zero element *(i, j)*
149
+
- `void add_element(const daecpp::int_type ind_i, const std::vector<int_type> &ind_j)` to define an entire row *i* of non-zero elements using a vector of indices *j*
150
+
- `void clear()` to clear the array of non-zero elements
151
+
- `void reserve(const daecpp::int_type N_elements)` to reserve memory for `N_elements` non-zero elements of the Jacobian (which is recommended)
152
+
153
+
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).
154
+
155
+
Here is a simple example of defining 3 non-zero elements of the Jacobian:
add_element(1, {0, 1}); // Adds elements (1, 0) and (1, 1)
165
+
}
166
+
};
167
+
```
168
+
169
+
{: .important }
170
+
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.
171
+
172
+
For more details, refer to the [Jacobian shape](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_shape/jacobian_shape.cpp) example.
173
+
174
+
## Jacobian matrix check
175
+
176
+
User-defined Jacobians defined either analytically or using Jacobian shape can be conveniently compared to fully automatic Jacobian for debug purposes.
177
+
This can be achieved by using `daecpp::JacobianCompare` helper class.
178
+
The constructor of this class takes the user-defined Jacobian (defined analytically or from the Jacobian shape) and the vector function (RHS).
179
+
Then the object of the `daecpp::JacobianCompare` class can be called as a functor with the given state `x` and, optionally, time `t`.
180
+
181
+
Consider the following example:
182
+
183
+
```cpp
184
+
// Fill vectors x at which the Jacobian matrix will be tested
185
+
daecpp::state_vector x0 = {1.0, 2.0, 3.0, 4.0};
186
+
daecpp::state_vector x1 = {-0.5, -0.1, 0.5, 0.1};
187
+
188
+
// Check user-defined Jacobian at `x0` and `x1` and times t = 1 and 2
189
+
{
190
+
auto jac_comparison = daecpp::JacobianCompare(UserJacobian(), UserRHS());
191
+
192
+
auto N_errors_1 = jac_comparison(x0, 1.0);
193
+
auto N_errors_2 = jac_comparison(x1, 2.0);
194
+
}
195
+
```
196
+
197
+
The functor `jac_comparison` returns the number of errors found in the Jacobian.
198
+
It also prints on screen a summary of all inconsistencies found in the user-defined Jacobian compared to the fully automatic one.
199
+
An example of such output is given below:
200
+
201
+
```txt
202
+
Jacobian matrix comparison summary at time t = 1:
203
+
-- Found 4 difference(s) compared to the automatic (reference) Jacobian:
- element (0, 3) is incorrect (should be 0.5, not 1)
218
+
- element (3, 3) is not defined at all
219
+
220
+
{: .note }
221
+
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.
222
+
223
+
For more details, refer to the [Jacobian compare](https://github.com/dae-cpp/dae-cpp/blob/master/examples/jacobian_compare/jacobian_compare.cpp) example.
## Element-by-element vector function to define the Jacobian shape
85
+
86
+
For relatively small systems, there is no need to provide the Jacobian matrix (neither its shape), as it will be computed fast enough automatically.
87
+
If the user decides to provide manually derived Jacobian, it can be done using [Jacobian Matrix](jacobian-matrix.html) class.
88
+
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.
89
+
90
+
The vector function from the [Examples](#examples) section above can be defined element-by-element by inheriting the `daecpp::VectorFunctionElements` class:
ERROR("Index i in UserDefinedVectorFunction is out of boundaries");
107
+
}
108
+
}
109
+
};
110
+
```
111
+
112
+
Unlike `daecpp::VectorFunction`, here we define each element of the vector function individually, depending on the index `i`.
113
+
This obviously comes with some computational time penalty compared to the case where we define the entire vector function as an array.
114
+
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).
115
+
116
+
{: .note }
117
+
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.
118
+
119
+
After defining the vector function this way, define the [Jacobian matrix shape](jacobian-matrix.html#jacobian-matrix-shape).
120
+
121
+
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.
0 commit comments