heat_equation/
heat_equation.rs

1// cargo run --example heat_equation
2use std::fs;
3use ohsl::{Vec64, Mesh2D, Tridiagonal};
4
5fn main() {
6    println!("----- Heat equation -----");
7
8    const T_MAX: f64 = 1.0;         // Maximum simulation time
9    const N: usize = 200;           // Number of time steps
10    const J: usize = 50;            // Number of space steps
11    let dt = T_MAX / N as f64;      // Temporal step size
12    let dx = 1.0 / J as f64;        // Spatial step size
13    let alpha = 1.0;                // Thermal diffusion coefficient
14
15    println!( " * Solving the heat equation u_t = \u{03B1} u_xx in the domain" );
16    println!( " * t=[0,{}] and x=[0,1] subject to the boundary conditions", T_MAX );
17    println!( " * u(0,t) = u(1,t) = 0 and the initial condition u(x,0) = 100 * sin( pi * x ).");
18    println!( " * The time step size is dt = {} and the spacial step size is dx = {}.", dt, dx );
19    println!( " * The number of time and space steps are N = {} and J = {} respectively.", N, J);
20    println!( " * The thermal diffusion coefficient is \u{03B1} = {}.", alpha );
21
22    let mu = ( alpha * dt ) / ( dx * dx ); 
23
24    // Mesh for storing the solution
25    let tnodes = Vec64::linspace( 0.0, T_MAX, N + 1 );
26    let xnodes = Vec64::linspace( 0.0, 1.0, J + 1 );
27    let mut solution = Mesh2D::<f64>::new( xnodes, tnodes, 1 );
28
29    // Matrices for Crank-Nicolson method including BCs u(0,t) = u(1,t) = 0
30    let mut implicit = Tridiagonal::with_elements( -mu, 2. + 2. * mu, -mu, J + 1 );
31    implicit[( 0, 0 )]     = 1.0;
32    implicit[( 0, 1 )]     = 0.0;
33    implicit[( J, J - 1 )] = 0.0;
34    implicit[( J, J )]     = 1.0;
35
36    let mut explicit = Tridiagonal::with_elements( mu, 2. - 2. * mu, mu, J + 1 );
37    explicit[( 0, 0 )]     =   1.0;
38    explicit[( 0, 1 )]     =   0.0;
39    explicit[( J, J - 1 )] =   0.0;
40    explicit[( J, J )]     = - 1.0;
41
42    // Initial condition u(x,0) = 100 * sin( pi * x )
43    let mut current = Vec64::zeros( J + 1 );
44    current[ 0 ] = 0.0;
45    solution[( 0, 0 )][ 0 ] = current[ 0 ];
46    for j in 1..J {
47        current[ j ] = 100.0 * ( std::f64::consts::PI * ( j as f64 ) * dx ).sin();
48        solution[( j, 0 )][ 0 ] = current[ j ];
49    }
50    current[ J ] = 0.0;
51    solution[( J, 0 )][ 0 ] = current[ J ];
52    let mut next: Vec64;
53
54    println!( " * Solving the heat equation using the Crank-Nicolson method..." );
55    // Time loop
56    for i in 1..N + 1 {
57        current = &explicit * &current;
58        next = implicit.solve( &current );
59        current = next;
60        for j in 0..J + 1 {
61            solution[( j, i )][ 0 ] = current[ j ];
62        }
63    }
64    println!( " * Finished solving the heat equation." );
65    fs::create_dir_all( "DATA" ).expect( "Could not create DATA output directory" );
66    let mut string = String::from( "DATA/Heat_equation" );
67    string.push_str( format!( "_J_{}_N_{}.dat", J, N ).as_str() );
68    solution.output( &string, 6 );
69
70    print!( " * For an animation of the solution run: " );
71    println!( "python3 Plotting/Heat_plot.py {} {}", J, N );
72    println!( " * and play the resulting file DATA/Heat_equation.avi" );
73    
74    println!( "--- FINISHED ---" );
75}
heat_equation/heat_equation.rs

heat_equation/
heat_equation.rs
