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, Default)]
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 #[default]
35 RSSCartesianStep,
36 /// An RSS state error control: when in doubt, use this error controller, especially for high accurracy.
37 ///
38 /// Here is the warning from GMAT R2016a on this error controller:
39 /// > This is a more stringent error control method than [`rss_step`] that is often used as the default in other software such as STK.
40 /// > If you set [the] accuracy to a very small number, 1e-13 for example, and set the error control to [`rss_step`], integrator
41 /// > performance will be poor, for little if any improvement in the accuracy of the orbit integration.
42 /// > For more best practices of these integrators (which clone those in GMAT), please refer to the
43 /// > [GMAT reference](https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/doc/help/src/Resource_NumericalIntegrators.xml#L1292).
44 /// > (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3004]
45 RSSState,
46 /// An RSS step error control which effectively computes the L2 norm of the provided Vector of size 3
47 ///
48 /// Note that this error controller should be preferably be used only with slices of a state with the same units.
49 /// For example, one should probably use this for position independently of using it for the velocity.
50 /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3045]
51 RSSStep,
52 /// A largest error control which effectively computes the largest error at each component
53 ///
54 /// This is a standard error computation algorithm, but it's arguably bad if the state's components have different units.
55 /// It calculates the largest local estimate of the error from the integration (`error_est`)
56 /// given the difference in the candidate state and the previous state (`state_delta`).
57 /// This error estimator is from the physical model estimator of GMAT
58 /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/PhysicalModel.cpp#L987]
59 LargestError,
60 /// A largest state error control
61 ///
62 /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3018]
63 LargestState,
64
65 /// A largest step error control which effectively computes the L1 norm of the provided Vector of size 3
66 ///
67 /// Note that this error controller should be preferably be used only with slices of a state with the same units.
68 /// For example, one should probably use this for position independently of using it for the velocity.
69 /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3033]
70 LargestStep,
71}
72
73impl ErrorControl {
74 /// Computes the actual error of the current step.
75 ///
76 /// The `error_est` is the estimated error computed from the difference in the two stages of
77 /// of the RK propagator. The `candidate` variable is the candidate state, and `cur_state` is
78 /// the current state. This function must return the error.
79 pub fn estimate<N: DimName>(
80 self,
81 error_est: &OVector<f64, N>,
82 candidate: &OVector<f64, N>,
83 cur_state: &OVector<f64, N>,
84 ) -> f64
85 where
86 DefaultAllocator: Allocator<N>,
87 {
88 match self {
89 ErrorControl::RSSCartesianState => {
90 if N::dim() >= 6 {
91 let err_radius = RSSState::estimate::<U3>(
92 &error_est.fixed_rows::<3>(0).into_owned(),
93 &candidate.fixed_rows::<3>(0).into_owned(),
94 &cur_state.fixed_rows::<3>(0).into_owned(),
95 );
96 let err_velocity = RSSState::estimate::<U3>(
97 &error_est.fixed_rows::<3>(3).into_owned(),
98 &candidate.fixed_rows::<3>(3).into_owned(),
99 &cur_state.fixed_rows::<3>(3).into_owned(),
100 );
101 err_radius.max(err_velocity)
102 } else {
103 RSSStep::estimate(error_est, candidate, cur_state)
104 }
105 }
106 ErrorControl::RSSCartesianStep => {
107 if N::dim() >= 6 {
108 let err_radius = RSSStep::estimate::<U3>(
109 &error_est.fixed_rows::<3>(0).into_owned(),
110 &candidate.fixed_rows::<3>(0).into_owned(),
111 &cur_state.fixed_rows::<3>(0).into_owned(),
112 );
113 let err_velocity = RSSStep::estimate::<U3>(
114 &error_est.fixed_rows::<3>(3).into_owned(),
115 &candidate.fixed_rows::<3>(3).into_owned(),
116 &cur_state.fixed_rows::<3>(3).into_owned(),
117 );
118 err_radius.max(err_velocity)
119 } else {
120 RSSStep::estimate(error_est, candidate, cur_state)
121 }
122 }
123 ErrorControl::RSSState => {
124 let mag = 0.5 * (candidate + cur_state).norm();
125 let err = error_est.norm();
126 if mag > REL_ERR_THRESH {
127 err / mag
128 } else {
129 err
130 }
131 }
132 ErrorControl::RSSStep => {
133 let mag = (candidate - cur_state).norm();
134 let err = error_est.norm();
135 if mag > REL_ERR_THRESH.sqrt() {
136 err / mag
137 } else {
138 err
139 }
140 }
141 ErrorControl::LargestError => {
142 let state_delta = candidate - cur_state;
143 let mut max_err = 0.0;
144 for (i, prop_err_i) in error_est.iter().enumerate() {
145 let err = if state_delta[i] > REL_ERR_THRESH {
146 (prop_err_i / state_delta[i]).abs()
147 } else {
148 prop_err_i.abs()
149 };
150 if err > max_err {
151 max_err = err;
152 }
153 }
154 max_err
155 }
156 ErrorControl::LargestState => {
157 let sum_state = candidate + cur_state;
158 let mut mag = 0.0f64;
159 let mut err = 0.0f64;
160 for i in 0..N::dim() {
161 mag += 0.5 * sum_state[i].abs();
162 err += error_est[i].abs();
163 }
164
165 if mag > REL_ERR_THRESH {
166 err / mag
167 } else {
168 err
169 }
170 }
171 ErrorControl::LargestStep => {
172 let state_delta = candidate - cur_state;
173 let mut mag = 0.0f64;
174 let mut err = 0.0f64;
175 for i in 0..N::dim() {
176 mag += state_delta[i].abs();
177 err += error_est[i].abs();
178 }
179
180 if mag > REL_ERR_THRESH {
181 err / mag
182 } else {
183 err
184 }
185 }
186 }
187 }
188}
189
190/// An RSS step error control which effectively computes the L2 norm of the provided Vector of size 3
191///
192/// Note that this error controller should be preferably be used only with slices of a state with the same units.
193/// For example, one should probably use this for position independently of using it for the velocity.
194/// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3045]
195#[derive(Clone, Copy)]
196#[allow(clippy::upper_case_acronyms)]
197struct RSSStep;
198impl RSSStep {
199 fn estimate<N: DimName>(
200 error_est: &OVector<f64, N>,
201 candidate: &OVector<f64, N>,
202 cur_state: &OVector<f64, N>,
203 ) -> f64
204 where
205 DefaultAllocator: Allocator<N>,
206 {
207 let mag = (candidate - cur_state).norm();
208 let err = error_est.norm();
209 if mag > REL_ERR_THRESH.sqrt() {
210 err / mag
211 } else {
212 err
213 }
214 }
215}
216
217/// An RSS state error control: when in doubt, use this error controller, especially for high accurracy.
218///
219/// Here is the warning from GMAT R2016a on this error controller:
220/// > This is a more stringent error control method than [`rss_step`] that is often used as the default in other software such as STK.
221/// > If you set [the] accuracy to a very small number, 1e-13 for example, and set the error control to [`rss_step`], integrator
222/// > performance will be poor, for little if any improvement in the accuracy of the orbit integration.
223/// > For more best practices of these integrators (which clone those in GMAT), please refer to the
224/// > [GMAT reference](https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/doc/help/src/Resource_NumericalIntegrators.xml#L1292).
225/// > (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3004]
226#[derive(Clone, Copy)]
227#[allow(clippy::upper_case_acronyms)]
228struct RSSState;
229impl RSSState {
230 fn estimate<N: DimName>(
231 error_est: &OVector<f64, N>,
232 candidate: &OVector<f64, N>,
233 cur_state: &OVector<f64, N>,
234 ) -> f64
235 where
236 DefaultAllocator: Allocator<N>,
237 {
238 let mag = 0.5 * (candidate + cur_state).norm();
239 let err = error_est.norm();
240 if mag > REL_ERR_THRESH {
241 err / mag
242 } else {
243 err
244 }
245 }
246}