solver/numeric/finite_difference/elliptic/problem.rs
1//! Stationary (elliptic) finite-difference problems.
2//!
3//! Unlike [`crate::numeric::finite_difference::pde::PdeProblem`], which
4//! marches backward in time from terminal data, an elliptic problem has no
5//! time variable. Its data are the boundary conditions on the domain, and its
6//! solution is the direct solve of
7//!
8//! ```text
9//! L u = f
10//! ```
11//!
12//! where `L` is the same spatial operator assembled from
13//! [`crate::numeric::finite_difference::discretization::Transport`]. In the
14//! standard elliptic sign convention the operator is `-T u + rho u` where `T`
15//! is the transport operator and `rho` is the reaction coefficient.
16
17use crate::core::grid::Grid;
18use crate::models::control::StateDerivatives;
19use crate::numeric::finite_difference::discretization::{
20 BoundaryCondition, BoundaryConditions, DimensionKind, Transport,
21};
22
23/// A stationary finite-difference problem `-T u + reaction * u = f` on a grid.
24///
25/// The operator is written as
26///
27/// ```text
28/// -T u + reaction * u = f
29/// ```
30///
31/// where `T u` is the transport term built from [`Transport`]. The negative
32/// sign follows the standard elliptic convention: for a diffusion coefficient
33/// `D`, `T u = D u''` and the problem reads `-D u'' + rho u = f`. The transport
34/// `source` is not used; use [`EllipticProblem::rhs`] for the forcing `f` and
35/// [`EllipticProblem::reaction`] for the local diagonal coefficient.
36pub trait EllipticProblem<const N: usize> {
37 /// Discretization kind for `dim`.
38 fn dimension_kind(&self, dim: usize) -> DimensionKind;
39
40 /// Transport coefficients for the current field `u` and derivatives.
41 ///
42 /// The returned `source` is ignored by the elliptic assembler.
43 fn transport(&self, state: &[f64; N], derivs: &StateDerivatives<N>) -> Transport<N>;
44
45 /// Forcing term `f(state)` on the right-hand side.
46 fn rhs(&self, state: &[f64; N]) -> f64;
47
48 /// Local zero-order coefficient multiplying `u` on the diagonal.
49 ///
50 /// Defaults to zero.
51 fn reaction(&self, _state: &[f64; N]) -> f64 {
52 0.0
53 }
54
55 /// Boundary conditions for each dimension.
56 fn boundary_conditions(&self) -> BoundaryConditions<N>;
57}
58
59/// Convenience access to the boundary condition on a single side.
60///
61/// # Examples
62///
63/// ```
64/// use solver::numeric::finite_difference::discretization::{
65/// BoundaryCondition, BoundaryConditions, DimensionKind, Transport,
66/// };
67/// use solver::numeric::finite_difference::elliptic::{
68/// EllipticProblem, EllipticProblemExt,
69/// };
70/// use solver::models::control::StateDerivatives;
71///
72/// struct Poisson;
73/// impl EllipticProblem<1> for Poisson {
74/// fn dimension_kind(&self, _dim: usize) -> DimensionKind { DimensionKind::Diffusion }
75/// fn transport(&self, _s: &[f64; 1], _d: &StateDerivatives<1>) -> Transport<1> {
76/// Transport::new([1.0], [1.0], 0.0)
77/// }
78/// fn rhs(&self, _s: &[f64; 1]) -> f64 { 1.0 }
79/// fn boundary_conditions(&self) -> BoundaryConditions<1> {
80/// BoundaryConditions::new(
81/// [BoundaryCondition::Dirichlet(0.0)],
82/// [BoundaryCondition::Dirichlet(1.0)],
83/// )
84/// }
85/// }
86/// assert!(matches!(Poisson.lower(0), BoundaryCondition::Dirichlet(_)));
87/// ```
88pub trait EllipticProblemExt<const N: usize>: EllipticProblem<N> {
89 /// Boundary condition at the lower side of `dim`.
90 fn lower(&self, dim: usize) -> BoundaryCondition {
91 self.boundary_conditions().lower[dim]
92 }
93
94 /// Boundary condition at the upper side of `dim`.
95 fn upper(&self, dim: usize) -> BoundaryCondition {
96 self.boundary_conditions().upper[dim]
97 }
98}
99
100impl<const N: usize, T: EllipticProblem<N>> EllipticProblemExt<N> for T {}
101
102/// A stationary (infinite-horizon) stochastic optimal control problem.
103///
104/// The stationary HJB equation is
105///
106/// ```text
107/// 0 = sup_u { f(x,u) + L^u V - r V }
108/// ```
109///
110/// where `f` is the running reward, `L^u V` is the controlled generator, and
111/// `r` is the discount rate. After optimizing over `u`, the equation is
112/// assembled in the same `-T V + r V = source` form as [`EllipticProblem`],
113/// with `T V + source = sup_u { f + L^u V }`.
114///
115/// This trait extends [`crate::models::control::ControlProblem`] with the
116/// discretization data needed by the stationary finite-difference solver. The
117/// `transport` method has the same semantics as
118/// [`crate::numeric::finite_difference::pde::PdeProblem::transport`]:
119/// `plus` and `minus` define `T`, and `source` is the residual that makes
120/// `T V + source` equal the optimized driver.
121pub trait EllipticControlProblem<const N: usize>:
122 crate::models::control::ControlProblem<N>
123{
124 /// Discretization kind for `dim`.
125 fn dimension_kind(&self, dim: usize) -> DimensionKind;
126
127 /// Transport coefficients and residual source for the given control and
128 /// derivatives.
129 ///
130 /// The returned `source` is used as the right-hand side of the stationary
131 /// operator. It must satisfy `T V + source = f + L^u V` at the current
132 /// derivative bundle, where `T V` is the discrete transport operator
133 /// built from `plus` and `minus`.
134 fn transport(
135 &self,
136 state: &[f64; N],
137 control: &Self::Control,
138 derivs: &StateDerivatives<N>,
139 ) -> Transport<N>;
140
141 /// Boundary conditions for each dimension.
142 ///
143 /// Defaults to zero Neumann (natural boundary) on every side.
144 fn boundary_conditions(&self) -> BoundaryConditions<N> {
145 BoundaryConditions::default()
146 }
147}
148
149/// Builds a zero field over the grid for elliptic iteration.
150pub fn zero_field<const N: usize>(grid: &Grid<N>) -> Vec<f64> {
151 vec![0.0; grid.total_size()]
152}