sparse_laplace/
sparse_laplace.rs

1// cargo run --example sparse_laplace
2
3extern crate ohsl;
4
5pub use ohsl::vector::{Vector, Vec64};
6pub use ohsl::matrix::{Matrix, Mat64};
7use ohsl::sparse::Sparse;
8
9use std::time::Instant;
10
11fn main() {
12    println!( "------------------ Laplace's equation ------------------" );
13    println!( "  * Solve Laplace's equation on the unit square subject");
14    println!( "  * to the given boundary conditions. The exact solution");
15    println!( "  * is u(x,y) = y/[(1+x)^2 + y^2].");
16    let start = Instant::now();
17    let n: usize = 64;                          // Number of intervals
18    let dx: f64 = 1.0 / ( n as f64 );           // Step size
19    let size: usize = ( n + 1 ) * ( n + 1 );    // Size of the linear system
20    let mut rhs = Vec64::zeros( size );         // Right hand side vector 
21
22    let mut triplets = Vec::new();
23    let mut row: usize = 0;
24    
25    // x = 0 boundary ( u = y / ( 1 + y^2 ) )
26    let i = 0;
27    for j in 0..n+1 {
28        let y = (j as f64) * dx;
29        triplets.push( ( row, i * ( n + 1 ) + j, 1.0 ) );
30        rhs[ row ] = y / ( 1.0 + y * y );
31        row += 1;
32    }
33
34    for i in 1..n {
35        let x = (i as f64) * dx;
36        // y = 0 boundary ( u = 0 )
37        let j = 0;
38        triplets.push( ( row, i * ( n + 1 ) + j, 1.0 ) );
39        rhs[ row ] = 0.0;
40        row += 1;
41        // Interior nodes
42        for j in 1..n {
43            triplets.push( ( row, ( i - 1 ) * ( n + 1 ) + j, 1.0 ) );
44            triplets.push( ( row, ( i + 1 ) * ( n + 1 ) + j, 1.0 ) );
45            triplets.push( ( row, i * ( n + 1 ) + j - 1, 1.0 ) );
46            triplets.push( ( row, i * ( n + 1 ) + j + 1, 1.0 ) );
47            triplets.push( ( row, i * ( n + 1 ) + j, - 4.0 ) );
48            rhs[ row ] = 0.0;
49            row += 1;
50        }
51        // y = 1 boundary ( u = 1 / ( ( 1 + x )^2 + 1 ) )
52        let j = n;
53        triplets.push( ( row, i * ( n + 1 ) + j, 1.0 ) );
54        rhs[ row ] = 1. / ( ( 1. + x ) * ( 1. + x ) + 1. );
55        row += 1;
56    }
57
58    // x = 1 boundary ( u = y / ( 4 + y^2 ) )
59    let i = n;
60    for j in 0..n+1 {
61        let y = (j as f64) * dx;
62        triplets.push( ( row, i * ( n + 1 ) + j, 1.0 ) );
63        rhs[ row ] = y / ( 4. + y * y );
64        row += 1;
65    }
66
67    // Exact solution u = y / ( ( 1 + x )^2 + y^2 )
68    let mut u_exact = Vec64::zeros( size );    
69    row = 0;
70    for i in 0..n+1 {
71        let x = (i as f64) * dx;
72        for j in 0..n+1 {
73            let y = (j as f64) * dx;
74            u_exact[ row ] = y / ( ( 1. + x ) * ( 1. + x ) + y * y );
75            row += 1;
76        }
77    }
78    // Create the sparse matrix from the triplets
79    let sparse = Sparse::<f64>::from_triplets( size, size, &mut triplets );
80    let duration = start.elapsed();
81    println!("  * Time to create the matrix: {:?}", duration);
82    // Solve the sparse system
83    let mut u_guess = Vector::<f64>::random( size ); // A better guess would improve convergence
84    let max_iter = 1000;
85    let tol = 1e-8;
86    let result = sparse.solve_qmr( &rhs, &mut u_guess, max_iter, tol );
87    // Output time and error
88    let duration = start.elapsed();
89    match result {
90        Ok( iter ) => println!( "  * The sparse system converged in {} iterations.", iter ),
91        Err( error ) => println!( "  * The sparse system failed to converge, error = {}", error )
92    }
93    println!("  * Time to solve the system: {:?}", duration);
94    let u_diff = u_guess - u_exact;
95    println!("  * Solution error = {}", u_diff.norm_2() );
96    println!( "-------------------- FINISHED ---------------------" );
97}
sparse_laplace/sparse_laplace.rs

sparse_laplace/
sparse_laplace.rs
