1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
//! Russell - Rust Scientific Library
//!
//! `russell_pde`: Solvers for ordinary differential equations and differential algebraic equations
//!
//! **Important:** This crate depends on external libraries (non-Rust). Thus, please check the [Installation Instructions on the GitHub Repository](https://github.com/cpmech/russell).
//!
//! To account for the EBCs, two approaches are possible:
//!
//! 1. Use the system partitioning strategy (SPS)
//! 2. Use the Lagrange multipliers method (LMM)
//!
//! ## Approach 1: System partitioning strategy (SPS)
//!
//! Consider the following partitioning of the vectors `a` and `f` and the matrix `K`:
//!
//! ```text
//! ┌ ┐ ┌ ┐ ┌ ┐
//! │ K̄ Ǩ │ │ ̄a │ │ f̄ │
//! │ │ │ │ = │ │
//! │ Ḵ ̰K │ │ ǎ │ │ f̌ │
//! └ ┘ └ ┘ └ ┘
//! K a f
//! ```
//!
//! where `ā` (a-bar) is a reduced vector containing only the unknown values (i.e., non-EBC nodes), and `ǎ` (a-check)
//! is a reduced vector containing only the prescribed values (i.e., EBC nodes). `f̄` and `f̌` are the associated reduced
//! right-hand side vectors. The `K̄` (K-bar) matrix is the reduced discrete Laplacian operator and `Ǩ` (K-check) is a
//! *correction* matrix. The `Ḵ` (K-underline) and `K̰` (K-under-tilde) matrices are often not needed.
//!
//! Thus, the linear system to be solved is:
//!
//! ```text
//! K̄ ā = f̄ - Ǩ ǎ
//! ```
//!
//! ## Approach 2: Lagrange multipliers method (LMM)
//!
//! The LMM consists of augmenting the original linear system with additional equations:
//!
//! ```text
//! ┌ ┐ ┌ ┐ ┌ ┐
//! │ K Cᵀ │ │ a │ │ f │
//! │ │ │ │ = │ │
//! │ C 0 │ │ ℓ │ │ ǎ │
//! └ ┘ └ ┘ └ ┘
//! M A F
//! ```
//!
//! where `ℓ` is the vector of Lagrange multipliers, `C` is the constraints matrix, and `ǎ` is the vector of
//! prescribed values at EBC nodes. The constraints matrix `C` has a row for each EBC (prescribed) node and a column
//! for every node. Each row in `C` has a single `1` at the column corresponding to the EBC node, and `0`s elsewhere.
//!
//! # Examples
//!
//! Solve the Poisson equation in 1D with homogeneous Dirichlet boundary conditions:
//!
//! ```text
//! -d²ϕ/dx² = 1 on x ∈ [0, 1]
//!
//! ϕ(0) = 0
//! ϕ(1) = 0
//! ```
//!
//! The analytical solution is `ϕ(x) = (x - x²) / 2`.
//!
//! ```
//! use russell_lab::approx_eq;
//! use russell_pde::{EssentialBcs1d, Fdm1d, Grid1d, NaturalBcs1d, Side, StrError};
//!
//! fn main() -> Result<(), StrError> {
//! // grid
//! let xmin = 0.0;
//! let xmax = 1.0;
//! let nx = 4;
//! let mut grid = Grid1d::new_uniform(xmin, xmax, nx)?;
//!
//! // Essential BCs
//! let mut ebcs = EssentialBcs1d::new();
//! ebcs.set(Side::Xmin, |_| 0.0);
//! ebcs.set(Side::Xmax, |_| 0.0);
//!
//! // Natural BCs (none)
//! let nbcs = NaturalBcs1d::new();
//!
//! // FDM solver
//! let kx = 1.0;
//! let fdm = Fdm1d::new(grid, ebcs, nbcs, kx)?;
//!
//! // Solve system
//! let alpha = 0.0; // Poisson
//! let source = |_| 1.0;
//! let phi = fdm.solve_sps(alpha, source)?;
//!
//! // Check
//! fdm.for_each_coord(|m, x| {
//! let analytical = x * (1.0 - x) / 2.0;
//! approx_eq(phi[m], analytical, 1e-14);
//! });
//! Ok(())
//! }
//! ```
/// Defines the error output as a static string
pub type StrError = &'static str;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
pub use *;
// run code from README file
;