nyx_space/propagators/error_ctrl.rs
1/*
2 Nyx, blazing fast astrodynamics
3 Copyright (C) 2018-onwards Christopher Rabotin <christopher.rabotin@gmail.com>
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU Affero General Public License as published
7 by the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU Affero General Public License for more details.
14
15 You should have received a copy of the GNU Affero General Public License
16 along with this program. If not, see <https://www.gnu.org/licenses/>.
17*/
18
19use serde::{Deserialize, Serialize};
20
21use crate::linalg::allocator::Allocator;
22use crate::linalg::{DefaultAllocator, DimName, OVector, U3};
23
24// This determines when to take into consideration the magnitude of the state_delta and
25// prevents dividing by too small of a number.
26const REL_ERR_THRESH: f64 = 0.1;
27
28/// The Error Control manages how a propagator computes the error in the current step.
29#[derive(Copy, Clone, Debug, PartialEq, Eq, Serialize, Deserialize)]
30pub enum ErrorControl {
31 /// An RSS state error control which effectively for the provided vector composed of two vectors of the same unit, both of size 3 (e.g. position + velocity).
32 RSSCartesianState,
33 /// An RSS step error control which effectively for the provided vector composed of two vectors of the same unit, both of size 3 (e.g. position + velocity).
34 RSSCartesianStep,
35 /// An RSS state error control: when in doubt, use this error controller, especially for high accurracy.
36 ///
37 /// Here is the warning from GMAT R2016a on this error controller:
38 /// > This is a more stringent error control method than [`rss_step`] that is often used as the default in other software such as STK.
39 /// > If you set [the] accuracy to a very small number, 1e-13 for example, and set the error control to [`rss_step`], integrator
40 /// > performance will be poor, for little if any improvement in the accuracy of the orbit integration.
41 /// > For more best practices of these integrators (which clone those in GMAT), please refer to the
42 /// > [GMAT reference](https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/doc/help/src/Resource_NumericalIntegrators.xml#L1292).
43 /// > (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3004]
44 RSSState,
45 /// An RSS step error control which effectively computes the L2 norm of the provided Vector of size 3
46 ///
47 /// Note that this error controller should be preferably be used only with slices of a state with the same units.
48 /// For example, one should probably use this for position independently of using it for the velocity.
49 /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3045]
50 RSSStep,
51 /// A largest error control which effectively computes the largest error at each component
52 ///
53 /// This is a standard error computation algorithm, but it's arguably bad if the state's components have different units.
54 /// It calculates the largest local estimate of the error from the integration (`error_est`)
55 /// given the difference in the candidate state and the previous state (`state_delta`).
56 /// This error estimator is from the physical model estimator of GMAT
57 /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/PhysicalModel.cpp#L987]
58 LargestError,
59 /// A largest state error control
60 ///
61 /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3018]
62 LargestState,
63
64 /// A largest step error control which effectively computes the L1 norm of the provided Vector of size 3
65 ///
66 /// Note that this error controller should be preferably be used only with slices of a state with the same units.
67 /// For example, one should probably use this for position independently of using it for the velocity.
68 /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3033]
69 LargestStep,
70}
71
72impl ErrorControl {
73 /// Computes the actual error of the current step.
74 ///
75 /// The `error_est` is the estimated error computed from the difference in the two stages of
76 /// of the RK propagator. The `candidate` variable is the candidate state, and `cur_state` is
77 /// the current state. This function must return the error.
78 pub fn estimate<N: DimName>(
79 self,
80 error_est: &OVector<f64, N>,
81 candidate: &OVector<f64, N>,
82 cur_state: &OVector<f64, N>,
83 ) -> f64
84 where
85 DefaultAllocator: Allocator<N>,
86 {
87 match self {
88 ErrorControl::RSSCartesianState => {
89 if N::dim() >= 6 {
90 let err_radius = RSSState::estimate::<U3>(
91 &error_est.fixed_rows::<3>(0).into_owned(),
92 &candidate.fixed_rows::<3>(0).into_owned(),
93 &cur_state.fixed_rows::<3>(0).into_owned(),
94 );
95 let err_velocity = RSSState::estimate::<U3>(
96 &error_est.fixed_rows::<3>(3).into_owned(),
97 &candidate.fixed_rows::<3>(3).into_owned(),
98 &cur_state.fixed_rows::<3>(3).into_owned(),
99 );
100 err_radius.max(err_velocity)
101 } else {
102 RSSStep::estimate(error_est, candidate, cur_state)
103 }
104 }
105 ErrorControl::RSSCartesianStep => {
106 if N::dim() >= 6 {
107 let err_radius = RSSStep::estimate::<U3>(
108 &error_est.fixed_rows::<3>(0).into_owned(),
109 &candidate.fixed_rows::<3>(0).into_owned(),
110 &cur_state.fixed_rows::<3>(0).into_owned(),
111 );
112 let err_velocity = RSSStep::estimate::<U3>(
113 &error_est.fixed_rows::<3>(3).into_owned(),
114 &candidate.fixed_rows::<3>(3).into_owned(),
115 &cur_state.fixed_rows::<3>(3).into_owned(),
116 );
117 err_radius.max(err_velocity)
118 } else {
119 RSSStep::estimate(error_est, candidate, cur_state)
120 }
121 }
122 ErrorControl::RSSState => {
123 let mag = 0.5 * (candidate + cur_state).norm();
124 let err = error_est.norm();
125 if mag > REL_ERR_THRESH {
126 err / mag
127 } else {
128 err
129 }
130 }
131 ErrorControl::RSSStep => {
132 let mag = (candidate - cur_state).norm();
133 let err = error_est.norm();
134 if mag > REL_ERR_THRESH.sqrt() {
135 err / mag
136 } else {
137 err
138 }
139 }
140 ErrorControl::LargestError => {
141 let state_delta = candidate - cur_state;
142 let mut max_err = 0.0;
143 for (i, prop_err_i) in error_est.iter().enumerate() {
144 let err = if state_delta[i] > REL_ERR_THRESH {
145 (prop_err_i / state_delta[i]).abs()
146 } else {
147 prop_err_i.abs()
148 };
149 if err > max_err {
150 max_err = err;
151 }
152 }
153 max_err
154 }
155 ErrorControl::LargestState => {
156 let sum_state = candidate + cur_state;
157 let mut mag = 0.0f64;
158 let mut err = 0.0f64;
159 for i in 0..N::dim() {
160 mag += 0.5 * sum_state[i].abs();
161 err += error_est[i].abs();
162 }
163
164 if mag > REL_ERR_THRESH {
165 err / mag
166 } else {
167 err
168 }
169 }
170 ErrorControl::LargestStep => {
171 let state_delta = candidate - cur_state;
172 let mut mag = 0.0f64;
173 let mut err = 0.0f64;
174 for i in 0..N::dim() {
175 mag += state_delta[i].abs();
176 err += error_est[i].abs();
177 }
178
179 if mag > REL_ERR_THRESH {
180 err / mag
181 } else {
182 err
183 }
184 }
185 }
186 }
187}
188
189impl Default for ErrorControl {
190 fn default() -> Self {
191 Self::RSSCartesianStep
192 }
193}
194
195/// An RSS step error control which effectively computes the L2 norm of the provided Vector of size 3
196///
197/// Note that this error controller should be preferably be used only with slices of a state with the same units.
198/// For example, one should probably use this for position independently of using it for the velocity.
199/// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3045]
200#[derive(Clone, Copy)]
201#[allow(clippy::upper_case_acronyms)]
202struct RSSStep;
203impl RSSStep {
204 fn estimate<N: DimName>(
205 error_est: &OVector<f64, N>,
206 candidate: &OVector<f64, N>,
207 cur_state: &OVector<f64, N>,
208 ) -> f64
209 where
210 DefaultAllocator: Allocator<N>,
211 {
212 let mag = (candidate - cur_state).norm();
213 let err = error_est.norm();
214 if mag > REL_ERR_THRESH.sqrt() {
215 err / mag
216 } else {
217 err
218 }
219 }
220}
221
222/// An RSS state error control: when in doubt, use this error controller, especially for high accurracy.
223///
224/// Here is the warning from GMAT R2016a on this error controller:
225/// > This is a more stringent error control method than [`rss_step`] that is often used as the default in other software such as STK.
226/// > If you set [the] accuracy to a very small number, 1e-13 for example, and set the error control to [`rss_step`], integrator
227/// > performance will be poor, for little if any improvement in the accuracy of the orbit integration.
228/// > For more best practices of these integrators (which clone those in GMAT), please refer to the
229/// > [GMAT reference](https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/doc/help/src/Resource_NumericalIntegrators.xml#L1292).
230/// > (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3004]
231#[derive(Clone, Copy)]
232#[allow(clippy::upper_case_acronyms)]
233struct RSSState;
234impl RSSState {
235 fn estimate<N: DimName>(
236 error_est: &OVector<f64, N>,
237 candidate: &OVector<f64, N>,
238 cur_state: &OVector<f64, N>,
239 ) -> f64
240 where
241 DefaultAllocator: Allocator<N>,
242 {
243 let mag = 0.5 * (candidate + cur_state).norm();
244 let err = error_est.norm();
245 if mag > REL_ERR_THRESH {
246 err / mag
247 } else {
248 err
249 }
250 }
251}