connor/nrlmsise00
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
23d2be94448654fee9f95d4d51ff2087d42efe06
nrlmsise00/tests/derivatives.rs
T

452 lines
18 KiB
Rust
Raw Blame History

//! That the model differentiates.
//!
//! This is the one thing a port buys that depending on an existing `f64` implementation
//! would not. Density falls by a factor of e about every fifty kilometers, so a drag
//! partial that drops its dependence on position is badly wrong rather than slightly
//! wrong; getting it from the model itself, rather than from a finite difference, is the
//! reason the model is generic over its scalar type.
//!
//! Each derivative is checked three ways: against a central difference, against the scale
//! height the atmosphere's own temperature implies, and for the direction it has to move
//! in. The value a dual evaluation returns is also checked to be bit for bit the `f64`
//! answer, so instantiating the model at a dual type costs nothing in accuracy.
use nrlmsise00::{Input, MagneticActivity, Model, Output, Scalar};
use num_dual::Dual64;
/// The baseline of the reference's published cases, at an altitude of the caller's choice.
fn input<T: Scalar>(altitude_km: T) -> Input<T> {
Input {
day_of_year: 172,
seconds_of_day: T::lift(29000.0),
altitude_km,
latitude_deg: T::lift(60.0),
longitude_deg: T::lift(-70.0),
local_solar_time_h: T::lift(16.0),
f107_average: T::lift(150.0),
f107: T::lift(150.0),
magnetic: MagneticActivity::History {
ap: [100.0; 7].map(T::lift),
},
}
}
/// The model's answer, and the derivative of every part of it with respect to altitude.
fn with_derivatives(model: &Model, altitude_km: f64) -> (Output<f64>, Output<f64>) {
let seeded = Dual64::from(altitude_km).derivative();
let dual = model.gtd7d(&input(seeded));
let value = |x: Dual64| x.re;
let slope = |x: Dual64| x.eps;
let project = |pick: &dyn Fn(Dual64) -> f64| Output {
helium: pick(dual.helium),
atomic_oxygen: pick(dual.atomic_oxygen),
molecular_nitrogen: pick(dual.molecular_nitrogen),
molecular_oxygen: pick(dual.molecular_oxygen),
argon: pick(dual.argon),
total_mass: pick(dual.total_mass),
hydrogen: pick(dual.hydrogen),
atomic_nitrogen: pick(dual.atomic_nitrogen),
anomalous_oxygen: pick(dual.anomalous_oxygen),
exospheric_temperature: pick(dual.exospheric_temperature),
temperature: pick(dual.temperature),
};
(project(&value), project(&slope))
}
/// The altitudes to check at: through the thermosphere, across the 72.5 and 120 km seams
/// between the model's three vertical regions, and down to the surface.
/// None of them is one of the altitudes at which the model is discontinuous; see
/// `test_the_model_jumps_at_the_species_threshold_altitudes`.
const ALTITUDES: [f64; 13] = [
0.0, 20.0, 50.0, 71.0, 74.0, 90.0, 110.0, 130.0, 210.0, 400.0, 600.0, 800.0, 1000.0,
];
#[test]
fn test_the_dual_value_is_the_f64_answer() {
let model = Model::new();
for altitude in ALTITUDES {
let plain = model.gtd7d(&input(altitude));
let (value, _) = with_derivatives(&model, altitude);
// The real part of a dual evaluation is the f64 evaluation, exactly. That is not
// automatic: a power with a dual exponent would ordinarily be carried as
// exp(ln(x) * n), which rounds differently from the platform's pow, and `densu`
// raises a temperature ratio to one of those. See `Scalar::powt`.
for (dual, plain) in value.densities().iter().zip(plain.densities().iter()) {
if *plain == 0.0 {
assert_eq!(*dual, 0.0);
continue;
}
let relative = (dual - plain).abs() / plain.abs();
assert!(
relative < 1.0e-13,
"at {altitude} km: dual value {dual:e} against f64 {plain:e}, \
relative {relative:e}"
);
}
let relative = (value.temperature - plain.temperature).abs() / plain.temperature;
assert!(
relative < 1.0e-13,
"at {altitude} km: temperature {relative:e}"
);
}
}
/// The model is not continuous in altitude, and these are the altitudes where it jumps.
///
/// Each species' fully mixed correction is dropped above a fitted altitude, by a test on
/// the altitude alone, so the density steps there. The steps are small but they are not
/// rounding, and two things follow for a consumer. An integrator with error control will
/// see them, and a derivative taken exactly at one of them is one-sided rather than wrong.
///
/// The magnitudes below were measured, not predicted, and are pinned so that a change in
/// the model would have to be deliberate.
#[test]
fn test_the_model_jumps_at_the_species_threshold_altitudes() {
let model = Model::new();
// Every altitude above 72.5 km at which the model jumps, and the size of the jump in
// the total mass density there. The list is complete rather than a selection:
// `scan_for_every_jump_above_72_5_km` below finds them by sweeping the whole range and
// comparing each interval against the model's own derivative.
//
// Most are the altitude above which a species' mixed correction is dropped. Three are
// not: 72.5 km is the join between the model's two halves, 123.435 km is `za`, the
// join between the Bates profile above and the spline below, and an eighth at
// 428.475 km, too small to pin here, is where atomic nitrogen's turbopause blend
// saturates and `dnet` switches to the mixed density with a step of about exp(-10).
// Annotated, because `Scalar` is blanket-implemented over `num-dual`'s trait and a
// bare float literal flowing into the generic model no longer falls back to `f64`.
let thresholds: [(f64, f64); 7] = [
(72.5, 1.0e-3),
(123.435, 1.5e-3),
(160.0, 7.1e-5),
(240.0, 2.3e-6),
(250.0, 1.9e-5),
(300.0, 1.7e-4),
(450.0, 3.5e-5),
];
for (altitude, expected) in thresholds {
let step = 1.0e-9;
let below = model.gtd7d(&input(altitude - step));
let above = model.gtd7d(&input(altitude + step));
let jump = (above.total_mass - below.total_mass).abs() / below.total_mass;
assert!(
jump > expected / 2.0 && jump < expected * 2.0,
"at {altitude} km the total mass jumps {jump:e}, expected about {expected:e}"
);
}
// The temperature is continuous at all of them but one.
for (altitude, _) in thresholds {
let step = 1.0e-9;
let below = model.gtd7d(&input(altitude - step));
let above = model.gtd7d(&input(altitude + step));
let jump = (above.temperature - below.temperature).abs() / below.temperature;
if altitude == 123.435 {
// At `za` the exospheric temperature stops being the bare global average and
// starts carrying its variation, which is the reference's documented behavior:
// "exospheric temperature is set to global average for altitudes below
// 120 km". The step is 1027.3 K to 1426.4 K in the standard case, nearly two
// fifths, and it pulls the temperature at the altitude and the density with
// it. This is the largest discontinuity in the model above 72.5 km.
assert!(
(1.0e-3..3.0e-3).contains(&jump),
"at za the temperature jumps {jump:e}"
);
let exospheric = (above.exospheric_temperature - below.exospheric_temperature).abs()
/ below.exospheric_temperature;
assert!(
(0.3..0.5).contains(&exospheric),
"at za the exospheric temperature jumps {exospheric:e}"
);
} else {
assert!(
jump < 1.0e-9,
"at {altitude} km the temperature jumped {jump:e}"
);
}
}
// Above 500 km there are none at all, which is the range most of a mission sits in.
let mut z: f64 = 500.0;
while z < 1000.0 {
let step = 1.0e-9;
let below = model.gtd7d(&input(z - step)).total_mass;
let above = model.gtd7d(&input(z + step)).total_mass;
assert!(
(above - below).abs() / below < 1.0e-9,
"the density jumps at {z} km"
);
z += 0.5;
}
}
/// The altitude derivative is the one a central difference converges to.
///
/// Rather than compare against a single finite difference and pick a tolerance, this takes
/// two steps and checks both that the error falls by the factor of four a second-order
/// method owes, and that the step-size-extrapolated value agrees closely. Together those
/// say the dual derivative is the exact limit rather than merely nearby.
#[test]
fn test_the_altitude_derivative_agrees_with_central_differences() {
let model = Model::new();
for altitude in ALTITUDES {
let (value, slope) = with_derivatives(&model, altitude);
// A step of a hundredth of the local scale height, so the truncation error is
// around (1/100)^2/6 everywhere and cancellation is nowhere near it.
let scale_height = -value.total_mass / slope.total_mass;
let h = scale_height / 100.0;
let difference = |h: f64| {
let above = model.gtd7d(&input(altitude + h));
let below = model.gtd7d(&input(altitude - h));
(above.total_mass - below.total_mass) / (2.0 * h)
};
let coarse = difference(h);
let fine = difference(h / 2.0);
// Richardson: the h^2 term cancels between the two. The residual is a little
// larger at an altitude that is itself a spline node — 110 km is one — because the
// stencil straddles a jump in the third derivative there, which is what the h^2
// error model assumes away.
let extrapolated = (fine * 4.0 - coarse) / 3.0;
let relative = (slope.total_mass - extrapolated).abs() / extrapolated.abs();
assert!(
relative < 1.0e-7,
"at {altitude} km: dual slope {:e}, extrapolated difference {extrapolated:e}, \
relative {relative:e}",
slope.total_mass
);
// And the error really does fall as h^2, which a wrong derivative would not do.
let coarse_error = (coarse - slope.total_mass).abs();
let fine_error = (fine - slope.total_mass).abs();
assert!(
fine_error * 3.0 < coarse_error && coarse_error < fine_error * 6.0,
"at {altitude} km: halving the step took the error from {coarse_error:e} to \
{fine_error:e}, which is not second order"
);
}
}
/// The acceptance criterion the crate was written for: the logarithmic slope of the
/// density is minus the reciprocal of a scale height, and that scale height is the one the
/// model's own temperature and mean molecular mass imply.
#[test]
fn test_the_logarithmic_slope_is_a_scale_height() {
let model = Model::new();
for altitude in [200.0, 300.0, 400.0, 600.0, 800.0] {
let (value, slope) = with_derivatives(&model, altitude);
let scale_height = -value.total_mass / slope.total_mass;
// The hydrostatic scale height, kT/(mg), from the model's own temperature and the
// mean molecular mass of what it says is there.
let number: f64 = value.densities()[..5].iter().sum::<f64>()
+ value.hydrogen
+ value.atomic_nitrogen
+ value.anomalous_oxygen;
let mean_mass = value.total_mass / number;
const BOLTZMANN: f64 = 1.380649e-23;
let radius = 6371.0e3 + altitude * 1000.0;
let gravity = 3.986004418e14 / (radius * radius);
let implied = BOLTZMANN * value.temperature / (mean_mass * gravity) / 1000.0;
assert!(
(10.0..200.0).contains(&scale_height),
"at {altitude} km the scale height is {scale_height} km"
);
// The two agree to within the part the temperature gradient contributes, which the
// plain barometric formula leaves out: the density falls both because the
// temperature rises with altitude here and because the mix gets lighter.
let relative = (scale_height - implied).abs() / implied;
assert!(
relative < 0.5,
"at {altitude} km: slope gives {scale_height:.1} km, kT/mg gives {implied:.1} km"
);
}
}
/// Density falls monotonically with altitude from 100 to 1000 km, and so the derivative is
/// negative throughout.
#[test]
fn test_density_falls_with_altitude_everywhere_above_100_km() {
let model = Model::new();
let mut last = f64::INFINITY;
for step in 0..=90 {
let altitude = 100.0 + f64::from(step) * 10.0;
let (value, slope) = with_derivatives(&model, altitude);
assert!(
slope.total_mass < 0.0,
"at {altitude} km the density derivative is {:e}",
slope.total_mass
);
assert!(
value.total_mass < last,
"at {altitude} km the density rose to {:e}",
value.total_mass
);
last = value.total_mass;
}
}
/// Finds every jump above 72.5 km, rather than checking the ones already known.
///
/// Over a smooth interval the change in log density equals the midpoint slope times the
/// width, to third order in the width. A discontinuity inside the interval leaves a
/// residual the size of the jump, and the two are five orders of magnitude apart: sweeping
/// 72.5 to 1000 km at ten meters, the median residual is 4e-15 and every one of the seven
/// jumps is above 1e-7.
///
/// Ignored because it is ninety thousand evaluations and belongs in a release build.
/// `cargo test --release --test derivatives -- --ignored --nocapture` runs it.
#[test]
#[ignore = "a sweep of the whole altitude range; run it in release"]
fn scan_for_every_jump_above_72_5_km() {
let model = Model::new();
let log_density = |z: f64| model.gtd7d(&input(z)).total_mass.ln();
let slope = |z: f64| {
let density = model.gtd7d(&input(Dual64::from(z).derivative())).total_mass;
density.eps / density.re
};
let width = 0.01;
let mut found = Vec::new();
let mut z = 72.5;
while z < 1000.0 {
let residual =
(log_density(z + width) - log_density(z) - width * slope(z + width / 2.0)).abs();
if residual > 1.0e-7 {
found.push(z);
}
z += width;
}
println!("{} intervals hold a jump:", found.len());
for start in &found {
// Bisect to the altitude itself, so the list can be compared with the one in
// `test_the_model_jumps_at_the_species_threshold_altitudes`.
let (mut low, mut high) = (*start, start + width);
for _ in 0..60 {
let middle = 0.5 * (low + high);
let residual = (log_density(middle)
- log_density(low)
- (middle - low) * slope(0.5 * (low + middle)))
.abs();
if residual > 1.0e-9 {
high = middle;
} else {
low = middle;
}
}
let below = model.gtd7d(&input(low)).total_mass;
let above = model.gtd7d(&input(high)).total_mass;
println!(
" {low:>10.5} km: total mass jumps {:.1e}",
(above - below).abs() / below
);
}
assert_eq!(
found.len(),
8,
"the known set is seven jumps plus the 72.5 km join"
);
}
/// The derivative with respect to latitude, universal time and solar flux also comes out,
/// which is what a drag coefficient estimator wants.
#[test]
fn test_the_other_inputs_differentiate_too() {
let model = Model::new();
let base = input::<f64>(400.0);
let check = |what: &str, seed: &dyn Fn(f64) -> Input<Dual64>, nominal: f64, h: f64| {
let slope = model.gtd7d(&seed(nominal)).total_mass.eps;
let plus = model.gtd7d(&into_f64(&seed(nominal + h))).total_mass;
let minus = model.gtd7d(&into_f64(&seed(nominal - h))).total_mass;
let difference = (plus - minus) / (2.0 * h);
let relative = (slope - difference).abs() / difference.abs();
assert!(
relative < 1.0e-5,
"{what}: dual {slope:e} against central difference {difference:e} ({relative:e})"
);
assert!(slope != 0.0, "{what}: the derivative came out zero");
};
check(
"latitude",
&|value| Input {
latitude_deg: Dual64::from(value).derivative(),
..lift(&base)
},
60.0,
0.01,
);
check(
"seconds of day",
&|value| Input {
seconds_of_day: Dual64::from(value).derivative(),
..lift(&base)
},
29000.0,
1.0,
);
check(
"f107",
&|value| Input {
f107: Dual64::from(value).derivative(),
..lift(&base)
},
150.0,
0.1,
);
}
/// The same input at a dual type, with no derivative seeded.
fn lift(from: &Input<f64>) -> Input<Dual64> {
Input {
day_of_year: from.day_of_year,
seconds_of_day: Dual64::lift(from.seconds_of_day),
altitude_km: Dual64::lift(from.altitude_km),
latitude_deg: Dual64::lift(from.latitude_deg),
longitude_deg: Dual64::lift(from.longitude_deg),
local_solar_time_h: Dual64::lift(from.local_solar_time_h),
f107_average: Dual64::lift(from.f107_average),
f107: Dual64::lift(from.f107),
magnetic: match from.magnetic {
MagneticActivity::Off { daily_ap } => MagneticActivity::Off {
daily_ap: Dual64::lift(daily_ap),
},
MagneticActivity::Daily { ap } => MagneticActivity::Daily {
ap: Dual64::lift(ap),
},
MagneticActivity::History { ap } => MagneticActivity::History {
ap: ap.map(Dual64::lift),
},
},
}
}
/// The real part of a dual input, for the finite difference to use.
fn into_f64(from: &Input<Dual64>) -> Input<f64> {
Input {
day_of_year: from.day_of_year,
seconds_of_day: from.seconds_of_day.re,
altitude_km: from.altitude_km.re,
latitude_deg: from.latitude_deg.re,
longitude_deg: from.longitude_deg.re,
local_solar_time_h: from.local_solar_time_h.re,
f107_average: from.f107_average.re,
f107: from.f107.re,
magnetic: match from.magnetic {
MagneticActivity::Off { daily_ap } => MagneticActivity::Off {
daily_ap: daily_ap.re,
},
MagneticActivity::Daily { ap } => MagneticActivity::Daily { ap: ap.re },
MagneticActivity::History { ap } => MagneticActivity::History {
ap: ap.map(|value| value.re),
},
},
}
}
Reference in new issue View Git Blame Copy Permalink