feat: implement mass balance validation for Story 7.1
- Added port_mass_flows to Component trait and implements for core components. - Added System::check_mass_balance and integrated it into the solver. - Restored connect methods for ExpansionValve, Compressor, and Pipe to fix integration tests. - Updated Python and C bindings for validation errors. - Updated sprint status and story documentation.
This commit is contained in:
@@ -98,6 +98,15 @@ pub enum SolverError {
|
||||
/// Human-readable description of the system defect.
|
||||
message: String,
|
||||
},
|
||||
|
||||
/// Post-solve validation failed (e.g., mass or energy balance violation).
|
||||
#[error("Validation failed: mass error = {mass_error:.3e} kg/s, energy error = {energy_error:.3e} W")]
|
||||
Validation {
|
||||
/// Mass balance error in kg/s
|
||||
mass_error: f64,
|
||||
/// Energy balance error in W
|
||||
energy_error: f64,
|
||||
},
|
||||
}
|
||||
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
@@ -1991,10 +2000,19 @@ impl Solver for SolverStrategy {
|
||||
},
|
||||
"SolverStrategy::solve dispatching"
|
||||
);
|
||||
match self {
|
||||
let result = match self {
|
||||
SolverStrategy::NewtonRaphson(cfg) => cfg.solve(system),
|
||||
SolverStrategy::SequentialSubstitution(cfg) => cfg.solve(system),
|
||||
};
|
||||
|
||||
if let Ok(state) = &result {
|
||||
if state.is_converged() {
|
||||
// Post-solve validation checks
|
||||
system.check_mass_balance(&state.state)?;
|
||||
}
|
||||
}
|
||||
|
||||
result
|
||||
}
|
||||
|
||||
fn with_timeout(self, timeout: Duration) -> Self {
|
||||
|
||||
@@ -1760,6 +1760,34 @@ impl System {
|
||||
let _ = row_offset; // avoid unused warning
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// Verifies that global mass balance is conserved.
|
||||
///
|
||||
/// Sums the mass flow rates at the ports of each component and ensures they
|
||||
/// sum to zero within a tight tolerance (1e-9 kg/s).
|
||||
pub fn check_mass_balance(&self, state: &StateSlice) -> Result<(), crate::SolverError> {
|
||||
let tolerance = 1e-9;
|
||||
let mut total_mass_error = 0.0;
|
||||
let mut has_violation = false;
|
||||
|
||||
for (_node_idx, component, _edge_indices) in self.traverse_for_jacobian() {
|
||||
if let Ok(mass_flows) = component.port_mass_flows(state) {
|
||||
let sum: f64 = mass_flows.iter().map(|m| m.to_kg_per_s()).sum();
|
||||
if sum.abs() > tolerance {
|
||||
has_violation = true;
|
||||
total_mass_error += sum.abs();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if has_violation {
|
||||
return Err(crate::SolverError::Validation {
|
||||
mass_error: total_mass_error,
|
||||
energy_error: 0.0,
|
||||
});
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
}
|
||||
|
||||
impl Default for System {
|
||||
@@ -3529,4 +3557,73 @@ mod tests {
|
||||
assert_eq!(indices.len(), 1);
|
||||
assert_eq!(indices[0].1, 2); // 2 * edge_count = 2
|
||||
}
|
||||
|
||||
struct BadMassFlowComponent {
|
||||
ports: Vec<ConnectedPort>,
|
||||
}
|
||||
|
||||
impl Component for BadMassFlowComponent {
|
||||
fn compute_residuals(
|
||||
&self,
|
||||
_state: &SystemState,
|
||||
_residuals: &mut entropyk_components::ResidualVector,
|
||||
) -> Result<(), ComponentError> {
|
||||
Ok(())
|
||||
}
|
||||
|
||||
fn jacobian_entries(
|
||||
&self,
|
||||
_state: &SystemState,
|
||||
_jacobian: &mut JacobianBuilder,
|
||||
) -> Result<(), ComponentError> {
|
||||
Ok(())
|
||||
}
|
||||
|
||||
fn n_equations(&self) -> usize {
|
||||
0
|
||||
}
|
||||
|
||||
fn get_ports(&self) -> &[ConnectedPort] {
|
||||
&self.ports
|
||||
}
|
||||
|
||||
fn port_mass_flows(&self, _state: &SystemState) -> Result<Vec<entropyk_core::MassFlow>, ComponentError> {
|
||||
Ok(vec![
|
||||
entropyk_core::MassFlow::from_kg_per_s(1.0),
|
||||
entropyk_core::MassFlow::from_kg_per_s(-0.5), // Intentionally unbalanced
|
||||
])
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_mass_balance_violation() {
|
||||
let mut system = System::new();
|
||||
|
||||
let inlet = Port::new(
|
||||
FluidId::new("R134a"),
|
||||
Pressure::from_bar(1.0),
|
||||
Enthalpy::from_joules_per_kg(400000.0),
|
||||
);
|
||||
let outlet = Port::new(
|
||||
FluidId::new("R134a"),
|
||||
Pressure::from_bar(1.0),
|
||||
Enthalpy::from_joules_per_kg(400000.0),
|
||||
);
|
||||
let (c1, c2) = inlet.connect(outlet).unwrap();
|
||||
|
||||
let comp = Box::new(BadMassFlowComponent {
|
||||
ports: vec![c1, c2], // Just to have ports
|
||||
});
|
||||
|
||||
let n0 = system.add_component(comp);
|
||||
system.add_edge(n0, n0).unwrap(); // Self-edge to avoid isolated node
|
||||
|
||||
system.finalize().unwrap();
|
||||
|
||||
// Ensure state is appropriately sized for finalize
|
||||
let state = vec![0.0; system.full_state_vector_len()];
|
||||
let result = system.check_mass_balance(&state);
|
||||
|
||||
assert!(result.is_err());
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user