use crate::*;
use assert_approx_eq::assert_approx_eq;
use core::f64::NAN;
use mathru::{
algebra::linear::Vector,
analysis::differential_equation::ordinary::{DormandPrince, ExplicitODE, ProportionalControl},
vector,
};
use test::Bencher;
struct HarmonicOde {}
impl DifferentialEquation<2> for HarmonicOde {
fn ode_dot_y(&self, _t: f64, y: &Coord<2>) -> (Coord<2>, bool) {
let x = y.0[0];
let v = y.0[1];
(Coord::<2>([v, -x]), true)
}
}
impl ExplicitODE<f64> for HarmonicOde {
fn func(&self, _t: &f64, y: &Vector<f64>) -> Vector<f64> {
let x = *y.get(0);
let v = *y.get(1);
vector![v; -x]
}
fn time_span(&self) -> (f64, f64) {
(0., 10.)
}
fn init_cond(&self) -> Vector<f64> {
vector![1.; 0.]
}
}
struct CoupledHarmonicState {
x1: f64,
x2: f64,
v1: f64,
v2: f64,
}
impl From<[f64; 4]> for CoupledHarmonicState {
fn from(ar: [f64; 4]) -> Self {
CoupledHarmonicState {
x1: ar[0],
x2: ar[1],
v1: ar[2],
v2: ar[3],
}
}
}
impl From<CoupledHarmonicState> for [f64; 4] {
fn from(s: CoupledHarmonicState) -> Self {
[s.x1, s.x2, s.v1, s.v2]
}
}
struct CoupledHarmonicStructOde {}
impl DifferentialEquation<4> for CoupledHarmonicStructOde {
fn ode_dot_y(&self, _t: f64, y: &Coord<4>) -> (Coord<4>, bool) {
let s = CoupledHarmonicState::from(y.0);
(
Coord::<4>([s.v1, s.v2, -s.x1 + 0.1 * s.x2, -s.x2 + 0.1 * s.x1]),
true,
)
}
}
struct CoupledHarmonicArOde {}
impl DifferentialEquation<4> for CoupledHarmonicArOde {
fn ode_dot_y(&self, _t: f64, y: &Coord<4>) -> (Coord<4>, bool) {
let x1 = y.0[0];
let x2 = y.0[1];
let v1 = y.0[2];
let v2 = y.0[3];
(Coord::<4>([v1, v2, -x1 + 0.1 * x2, -x2 + 0.1 * x1]), true)
}
}
struct YoungLaplaceOde {
phi_0: f64,
startsign: f64,
}
impl DifferentialEquation<2> for YoungLaplaceOde {
fn ode_dot_y(&self, pmphi: f64, y: &Coord<2>) -> (Coord<2>, bool) {
let phi = -pmphi * self.phi_0.signum();
let x = y.0[0];
let z = y.0[1];
#[allow(clippy::float_cmp)]
{
if (x * z - phi.sin()).signum() != self.startsign {
dbg!();
return (Coord::<2>([NAN, NAN]), false);
}
}
let dxdphi = x * phi.cos() / (x * z - phi.sin()); let dzdphi = x * phi.sin() / (x * z - phi.sin()); (
Coord::<2>([-self.phi_0.signum() * dxdphi, -self.phi_0.signum() * dzdphi]),
true,
)
}
}
struct YoungLaplace {
phi_0: f64,
x_contact: f64,
z_contact: f64,
startsign: f64,
}
impl ExplicitODE<f64> for YoungLaplace {
fn func(&self, pmphi: &f64, q: &Vector<f64>) -> Vector<f64> {
let phi = -pmphi * self.phi_0.signum();
let x = *q.get(0);
let z = *q.get(1);
#[allow(clippy::float_cmp)]
{
if (x * z - phi.sin()).signum() != self.startsign {
return vector![NAN; NAN];
}
}
let dxdphi = x * phi.cos() / (x * z - phi.sin()); let dzdphi = x * phi.sin() / (x * z - phi.sin()); vector![-self.phi_0.signum() * dxdphi; -self.phi_0.signum() * dzdphi]
}
fn time_span(&self) -> (f64, f64) {
(-self.phi_0.abs(), 0.)
}
fn init_cond(&self) -> Vector<f64> {
#[allow(clippy::float_cmp)]
{
assert!(self.startsign == 1. || self.startsign == -1.); }
vector![self.x_contact; self.z_contact]
}
}
#[test]
fn compare_with_mathru() {
let phi_0 = -0.7653981633974483;
let x_contact = 0.34641208585537364;
let z_contact = 0.4145142430976436;
let h_0 = 0.09026305092159702;
let fac = 0.9;
let fac_min = 0.2; let fac_max = 10.0; let n_max = 1e9 as u32;
let abs_tol = 1e-6;
let rel_tol = 1e-3;
let solver = ProportionalControl::new(n_max, h_0, fac, fac_min, fac_max, abs_tol, rel_tol);
let ode = YoungLaplaceOde {
phi_0,
startsign: (x_contact * z_contact - phi_0.sin()).signum(),
};
let myres = solve_ivp(
&ode,
(-phi_0.abs(), 0.),
Coord([x_contact, z_contact]),
|_, _| true,
1e-6,
1e-3,
);
let problem = YoungLaplace {
phi_0,
x_contact,
z_contact,
startsign: (x_contact * z_contact - phi_0.sin()).signum(),
};
let (_pm_phi, xz) = solver.solve(&problem, &DormandPrince::default()).unwrap();
let myres_y = match myres {
IvpResult::FinalTimeReached(y) => (y),
_ => panic!(),
};
assert_approx_eq!(myres_y.0[0], xz[xz.len() - 1].get(0), 1e-3);
assert_approx_eq!(myres_y.0[1], xz[xz.len() - 1].get(1), 1e-3);
}
#[bench]
fn b_mathru_young_laplace(b: &mut Bencher) {
b.iter(|| {
let phi_0 = -0.7653981633974483;
let x_contact = 0.34641208585537364;
let z_contact = 0.4145142430976436;
let h_0 = 0.09026305092159702;
let fac = 0.9;
let fac_min = 0.2; let fac_max = 10.0; let n_max = 1e9 as u32;
let abs_tol = 1e-6;
let rel_tol = 1e-3;
let solver = ProportionalControl::new(n_max, h_0, fac, fac_min, fac_max, abs_tol, rel_tol);
let problem = YoungLaplace {
phi_0,
x_contact,
z_contact,
startsign: (x_contact * z_contact - phi_0.sin()).signum(),
};
solver.solve(&problem, &DormandPrince::default()).unwrap()
});
}
#[bench]
fn b_young_laplace(b: &mut Bencher) {
b.iter(|| {
let phi_0 = -0.7653981633974483;
let x_contact = 0.34641208585537364;
let z_contact = 0.4145142430976436;
let ode = YoungLaplaceOde {
phi_0,
startsign: (x_contact * z_contact - phi_0.sin()).signum(),
};
solve_ivp(
&ode,
(-phi_0.abs(), 0.),
Coord([x_contact, z_contact]),
|_, _| true,
1e-6,
1e-3,
)
});
}
#[bench]
fn b_mathru_harmonic(b: &mut Bencher) {
b.iter(|| {
let h_0 = 0.09026305092159702;
let fac = 0.9;
let fac_min = 0.2; let fac_max = 10.0; let n_max = 1e9 as u32;
let abs_tol = 1e-6;
let rel_tol = 1e-3;
let solver = ProportionalControl::new(n_max, h_0, fac, fac_min, fac_max, abs_tol, rel_tol);
let problem = HarmonicOde {};
solver.solve(&problem, &DormandPrince::default()).unwrap()
});
}
#[bench]
fn b_harmonic(b: &mut Bencher) {
b.iter(|| {
let ode = HarmonicOde {};
solve_ivp(&ode, (0., 10.), Coord([1., 0.]), |_, _| true, 1e-6, 1e-3)
});
}
#[bench]
fn b_coupled_array(b: &mut Bencher) {
b.iter(|| {
let ode = CoupledHarmonicArOde {};
solve_ivp(
&ode,
(0., 100.),
Coord([1., 0., 0., 0.5]),
|_, _| true,
1e-6,
1e-3,
)
});
}
#[bench]
fn b_coupled_struct(b: &mut Bencher) {
b.iter(|| {
let ode = CoupledHarmonicStructOde {};
solve_ivp(
&ode,
(0., 100.),
Coord([1., 0., 0., 0.5]),
|_, _| true,
1e-6,
1e-3,
)
});
}