stormbird 0.9.0

// Copyright (C) 2024, NTNU
// Author: Jarle Vinje Kramer <jarlekramer@gmail.com; jarle.a.kramer@ntnu.no>
// License: GPL v3.0 (see separate file LICENSE or https://www.gnu.org/licenses/gpl-3.0.html)

//! Compare the result from different solvers

use std::time::Instant;

use stormath::type_aliases::Float;
use stormath::consts::TAU;

use crate::lifting_line::prelude::*;
use crate::lifting_line::simulation_builder::{
    SimulationBuilder,
    SimulationSettings,
    QuasiSteadySettings,
    DynamicSettings,
};
use crate::line_force_model::corrections::circulation::{
    CirculationCorrectionBuilder,
    prescribed::PrescribedCirculation
};

use super::test_setup::RectangularWing;
use super::elliptic_wing_theory::EllipticWingTheory;

#[test]
fn steady_lift() {
    let aspect_ratio = 5.0;
    let cl_zero_angle = 1.2;
    let angle_of_attack = Float::from(2.0).to_radians();

    let theory = EllipticWingTheory {
        cl_2d: cl_zero_angle + TAU * angle_of_attack,
        aspect_ratio
    };

    let cl_theory = theory.cl();
    let cd_theory = theory.cd();

    let wing_builder = RectangularWing {
        aspect_ratio,
        cl_zero_angle,
        angle_of_attack,
        nr_strips: 64,
        ..Default::default()
    }.build();

    let mut prescribed_wing_builder = wing_builder.clone();

    let mut prescribed_circulation = PrescribedCirculation::default();
    prescribed_circulation.curve_fit_shape_parameters = true;

    prescribed_wing_builder.circulation_correction = CirculationCorrectionBuilder::Prescribed(
        prescribed_circulation
    );

    let steady_settings  = QuasiSteadySettings::default();
    let dynamic_settings = DynamicSettings::default();

    let nr_time_steps = 200;

    let velocity = SpatialVector::from([1.2, 0.0, 0.0]);

    let time_step = 0.25 / velocity.length();

    let mut steady_sim = SimulationBuilder::new(
        wing_builder.clone(),
        SimulationSettings::QuasiSteady(steady_settings.clone())
    ).build();

    let mut prescribed_sim = SimulationBuilder::new(
        prescribed_wing_builder.clone(),
        SimulationSettings::QuasiSteady(steady_settings)
    ).build();

    let mut dynamic_sim = SimulationBuilder::new(
        wing_builder.clone(),
        SimulationSettings::Dynamic(dynamic_settings)
    ).build();

    let force_factor = steady_sim.line_force_model.total_force_factor(velocity.length());

    let dynamic_velocity_points = dynamic_sim.get_freestream_velocity_points();
    let static_velocity_points = steady_sim.get_freestream_velocity_points();

    let dynamic_velocity_freestream: Vec<SpatialVector> = vec![velocity; dynamic_velocity_points.len()];
    let static_velocity_freestream: Vec<SpatialVector> = vec![velocity; static_velocity_points.len()];

    let start = Instant::now();
    let result_steady  = steady_sim.do_step(
        0.0,
        time_step,
        &static_velocity_freestream
    );

    println!("Time to run steady simulation: {} secs", start.elapsed().as_secs_f64());
    println!("Number of steady iterations: {}", result_steady.iterations);
    println!("Residual for steady simulation: {}", result_steady.residual);

    let cd_steady = result_steady.integrated_forces_sum()[0] / force_factor;
    let cl_steady = result_steady.integrated_forces_sum()[1] / force_factor;

    let start = Instant::now();
    let result_prescribed = prescribed_sim.do_step(
        0.0,
        time_step,
        &static_velocity_freestream
    );

    println!("Time to run prescribed simulation: {} secs", start.elapsed().as_secs_f64());
    println!("Number of prescribed iterations: {}", result_prescribed.iterations);
    println!("Residual for prescribed simulation: {}", result_prescribed.residual);

    let cd_prescribed = result_prescribed.integrated_forces_sum()[0] / force_factor;
    let cl_prescribed = result_prescribed.integrated_forces_sum()[1] / force_factor;

    let mut cd_dynamic = 0.0;
    let mut cl_dynamic = 0.0;

    let mut result_dynamic = SimulationResult::default();
    let start = Instant::now();
    for i in 0..nr_time_steps {
        let time = (i as Float) * time_step;

        result_dynamic = dynamic_sim.do_step(
            time,
            time_step,
            &dynamic_velocity_freestream
        );

        cd_dynamic = result_dynamic.integrated_forces_sum()[0] / force_factor;
        cl_dynamic = result_dynamic.integrated_forces_sum()[1] / force_factor;
    }

    println!("Time to run dynamic simulation: {} secs", start.elapsed().as_secs_f64());
    println!("Number of dynamic iterations on last time step: {}", result_dynamic.iterations);
    println!("Residual for dynamic simulation: {}", result_dynamic.residual);


    println!("Theory");
    dbg!(cd_theory, cl_theory);
    println!("Steady");
    dbg!(cd_steady, cl_steady);
    println!("Prescribed");
    dbg!(cd_prescribed, cl_prescribed);
    println!("Dynamic");
    dbg!(cd_dynamic, cl_dynamic);

    let steady_cl_error = (cl_theory - cl_steady).abs() / cl_theory.abs();
    let steady_cd_error = (cd_theory - cd_steady).abs() / cd_theory.abs();

    let prescribed_cl_error = (cl_theory - cl_prescribed).abs() / cl_theory.abs();
    let prescribed_cd_error = (cd_theory - cd_prescribed).abs() / cd_theory.abs();

    let dynamic_cl_error = (cl_theory - cl_dynamic).abs() / cl_theory.abs();
    let dynamic_cd_error = (cd_theory - cd_dynamic).abs() / cd_theory.abs();

    let allowable_cd_error = 0.06;
    let allowable_cl_error = 0.07;

    assert!(steady_cl_error < allowable_cl_error, "Steady cl error: {}", steady_cl_error);
    assert!(steady_cd_error < allowable_cd_error, "Steady cd error: {}", steady_cd_error);
    assert!(prescribed_cl_error < allowable_cl_error, "Prescribed cl error: {}", prescribed_cl_error);
    assert!(prescribed_cd_error < allowable_cd_error, "Prescribed cd error: {}", prescribed_cd_error);
    assert!(dynamic_cl_error < allowable_cl_error, "Dynamic cl error: {}", dynamic_cl_error);
    assert!(dynamic_cd_error < allowable_cd_error, "Dynamic cd error: {}", dynamic_cd_error);
}