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*/
18use serde::{Deserialize, Serialize};
19use serde_dhall::StaticType;
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, StaticType)]
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 { err / mag } else { err }
127 }
128 ErrorControl::RSSStep => {
129 let mag = (candidate - cur_state).norm();
130 let err = error_est.norm();
131 if mag > REL_ERR_THRESH.sqrt() {
132 err / mag
133 } else {
134 err
135 }
136 }
137 ErrorControl::LargestError => {
138 let state_delta = candidate - cur_state;
139 let mut max_err = 0.0;
140 for (i, prop_err_i) in error_est.iter().enumerate() {
141 let err = if state_delta[i] > REL_ERR_THRESH {
142 (prop_err_i / state_delta[i]).abs()
143 } else {
144 prop_err_i.abs()
145 };
146 if err > max_err {
147 max_err = err;
148 }
149 }
150 max_err
151 }
152 ErrorControl::LargestState => {
153 let sum_state = candidate + cur_state;
154 let mut mag = 0.0f64;
155 let mut err = 0.0f64;
156 for i in 0..N::dim() {
157 mag += 0.5 * sum_state[i].abs();
158 err += error_est[i].abs();
159 }
160
161 if mag > REL_ERR_THRESH { err / mag } else { err }
162 }
163 ErrorControl::LargestStep => {
164 let state_delta = candidate - cur_state;
165 let mut mag = 0.0f64;
166 let mut err = 0.0f64;
167 for i in 0..N::dim() {
168 mag += state_delta[i].abs();
169 err += error_est[i].abs();
170 }
171
172 if mag > REL_ERR_THRESH { err / mag } else { err }
173 }
174 }
175 }
176}
177
178/// An RSS step error control which effectively computes the L2 norm of the provided Vector of size 3
179///
180/// Note that this error controller should be preferably be used only with slices of a state with the same units.
181/// For example, one should probably use this for position independently of using it for the velocity.
182/// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3045]
183#[derive(Clone, Copy)]
184#[allow(clippy::upper_case_acronyms)]
185struct RSSStep;
186impl RSSStep {
187 fn estimate<N: DimName>(
188 error_est: &OVector<f64, N>,
189 candidate: &OVector<f64, N>,
190 cur_state: &OVector<f64, N>,
191 ) -> f64
192 where
193 DefaultAllocator: Allocator<N>,
194 {
195 let mag = (candidate - cur_state).norm();
196 let err = error_est.norm();
197 if mag > REL_ERR_THRESH.sqrt() {
198 err / mag
199 } else {
200 err
201 }
202 }
203}
204
205/// An RSS state error control: when in doubt, use this error controller, especially for high accurracy.
206///
207/// Here is the warning from GMAT R2016a on this error controller:
208/// > This is a more stringent error control method than [`rss_step`] that is often used as the default in other software such as STK.
209/// > If you set [the] accuracy to a very small number, 1e-13 for example, and set the error control to [`rss_step`], integrator
210/// > performance will be poor, for little if any improvement in the accuracy of the orbit integration.
211/// > For more best practices of these integrators (which clone those in GMAT), please refer to the
212/// > [GMAT reference](https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/doc/help/src/Resource_NumericalIntegrators.xml#L1292).
213/// > (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3004]
214#[derive(Clone, Copy)]
215#[allow(clippy::upper_case_acronyms)]
216struct RSSState;
217impl RSSState {
218 fn estimate<N: DimName>(
219 error_est: &OVector<f64, N>,
220 candidate: &OVector<f64, N>,
221 cur_state: &OVector<f64, N>,
222 ) -> f64
223 where
224 DefaultAllocator: Allocator<N>,
225 {
226 let mag = 0.5 * (candidate + cur_state).norm();
227 let err = error_est.norm();
228 if mag > REL_ERR_THRESH { err / mag } else { err }
229 }
230}