blasius/
blasius.rs

1// cargo run --example blasius
2use std::fs;
3use ohsl::ode_bvp::ODEBVP;
4use ohsl::vector::{Vec64, Vector};
5
6fn blasius_eqn( y: &Vec64, _x: &f64 ) -> Vec64 {
7    let mut f = Vec64::new( 3, 0.0 );
8    f[0] = y[1];
9    f[1] = y[2];
10    f[2] = -y[0] * y[2];
11    /*
12        (y)' = y',
13        (y')' = y'',
14        (y'')' = - y * y'' ,
15        Blasius equation y''' + y * y'' = 0
16    */
17    f
18}
19
20fn plate_bc( y: &Vec64 ) -> Vector<Option<f64>> {
21    let mut bc = Vector::<Option<f64>>::new( 3, None );
22    bc[0] = Some( y[0] );       // y(0) = 0
23    bc[1] = Some( y[1] );       // y'(0) = 0
24    bc[2] = None;               // y''(0) free
25    bc
26}
27
28fn free_bc( y: &Vec64 ) -> Vector<Option<f64>> {
29    let mut bc = Vector::<Option<f64>>::new( 3, None );
30    bc[0] = None;               // y(0) free
31    bc[1] = Some( y[1] - 1.0 ); // y'(inf) = 1
32    bc[2] = None;               // y''(inf) free
33    bc
34}
35
36fn main() {
37    println!("----- Blasius equation -----");
38
39    const ETA_INF: f64 = 20.0;
40    const N: usize = 512;
41
42    println!( " * Solving the Blasius equation f''' + f * f'' = 0 in the domain" );
43    println!( " * \u{03B7}=[0,{}] subject to the boundary conditions", ETA_INF );
44    println!( " * f(0) = 0, f'(0) = 0, f' -> 1 as \u{03B7} -> \u{03B7}_\u{221E} = {}.", ETA_INF );
45    println!( " * The number of nodes in the discretisation is N = {}.", N );
46
47    let nodes = Vec64::powspace( 0.0, ETA_INF, N, 1.2 ); // More nodes near eta=0
48    let mut ode_bvp = ODEBVP::<f64>::new( nodes.clone(), 3 );
49    ode_bvp.iterations( 20 );
50
51    // Set the initial guess for the solution
52    for i in 0..N {
53        let eta = ode_bvp.solution.coord( i );
54        let exp_minus_eta = (-eta).exp();
55        ode_bvp.solution[i][0] = eta * ( 1.0 - exp_minus_eta ); 
56        ode_bvp.solution[i][1] = 1.0 - exp_minus_eta + eta * exp_minus_eta; 
57        ode_bvp.solution[i][2] = 2.0 * exp_minus_eta - eta * exp_minus_eta;
58    }
59
60    println!( " * Solving the Blasius equation using the finite difference method..." );
61    let result = ode_bvp.solve( &blasius_eqn, &plate_bc, &free_bc, None );  
62    assert!( result.is_ok() && !result.is_err() );
63
64    let fdd_0 = ode_bvp.solution.get_interpolated_vars( 0.0 )[2];
65    println!( " * The computed value of f''(0) is {:.8}.", fdd_0 ); 
66
67    let c= 1.65519036023f64;
68    let reference = 1. / c.powf( 1.5 );
69    println!( " * The reference value of f''(0) is {:.8}.", reference );
70    println!( " * The relative error is {:.8}%.", ( fdd_0 - reference ).abs() / reference * 100.0 );  
71
72    println!( " * Finished solving the Blasius equation." );
73    fs::create_dir_all( "DATA" ).expect( "Could not create DATA output directory" );
74    let mut string = String::from( "DATA/Blasius_equation" );
75    string.push_str( format!( "_N_{}.dat", N ).as_str() );
76    ode_bvp.solution.output( &string, 6 ); 
77
78    print!( " * For a plot of the solution run: " );
79    println!( "python3 Plotting/Blasius_plot.py {}", N );
80    println!( " * and view the resulting file DATA/Blasius_equation_N_{}.png", N );
81    
82    println!( "--- FINISHED ---" );
83}
blasius/blasius.rs

blasius/
blasius.rs
