diff --git a/.gitignore b/.gitignore index 03dd8bc..cee17ec 100644 --- a/.gitignore +++ b/.gitignore @@ -84,3 +84,15 @@ BUG_FIX_VALIDATION_SUMMARY.md src/hyddown/__pycache__/hdclass.cpython-312.pyc VALIDATION_REPORT.txt src/hyddown/__pycache__/hdclass.cpython-312.pyc +PHASE1_PROGRESS.md +src/hyddown/__pycache__/__init__.cpython-312.pyc +src/hyddown/__pycache__/__init__.cpython-312.pyc +PHASE1_INTEGRATION_GUIDE.md +PHASE1_COMPLETE.md +PHASE1_CERBERUS_ENHANCEMENTS.md +VALIDATOR_ANALYSIS.md +PHASE1.5_COMPLETE.md +PHASE2_EXTRACTION_PLAN.md +.gitignore +PHASE2.1_COMPLETE.md +VALIDATOR_TEST_FIX.md diff --git a/src/hyddown/__init__.py b/src/hyddown/__init__.py index 40bbce3..a6f8419 100644 --- a/src/hyddown/__init__.py +++ b/src/hyddown/__init__.py @@ -4,3 +4,6 @@ from .hdclass import * from .transport import * +from . import exceptions +from . import safety_checks +from .thermo_solver import ThermodynamicSolver diff --git a/src/hyddown/__pycache__/__init__.cpython-312.pyc b/src/hyddown/__pycache__/__init__.cpython-312.pyc index 00fbbb9..e009e68 100644 Binary files a/src/hyddown/__pycache__/__init__.cpython-312.pyc and b/src/hyddown/__pycache__/__init__.cpython-312.pyc differ diff --git a/src/hyddown/__pycache__/exceptions.cpython-312.pyc b/src/hyddown/__pycache__/exceptions.cpython-312.pyc new file mode 100644 index 0000000..8567e69 Binary files /dev/null and b/src/hyddown/__pycache__/exceptions.cpython-312.pyc differ diff --git a/src/hyddown/__pycache__/hdclass.cpython-312.pyc b/src/hyddown/__pycache__/hdclass.cpython-312.pyc index 519b9b2..a599d62 100644 Binary files a/src/hyddown/__pycache__/hdclass.cpython-312.pyc and b/src/hyddown/__pycache__/hdclass.cpython-312.pyc differ diff --git a/src/hyddown/__pycache__/safety_checks.cpython-312.pyc b/src/hyddown/__pycache__/safety_checks.cpython-312.pyc new file mode 100644 index 0000000..95b01fa Binary files /dev/null and b/src/hyddown/__pycache__/safety_checks.cpython-312.pyc differ diff --git a/src/hyddown/__pycache__/thermo_solver.cpython-312.pyc b/src/hyddown/__pycache__/thermo_solver.cpython-312.pyc new file mode 100644 index 0000000..e61419c Binary files /dev/null and b/src/hyddown/__pycache__/thermo_solver.cpython-312.pyc differ diff --git a/src/hyddown/__pycache__/validator.cpython-312.pyc b/src/hyddown/__pycache__/validator.cpython-312.pyc index aadc466..0bc2ff2 100644 Binary files a/src/hyddown/__pycache__/validator.cpython-312.pyc and b/src/hyddown/__pycache__/validator.cpython-312.pyc differ diff --git a/src/hyddown/exceptions.py b/src/hyddown/exceptions.py new file mode 100644 index 0000000..620e987 --- /dev/null +++ b/src/hyddown/exceptions.py @@ -0,0 +1,369 @@ +# HydDown hydrogen/other gas depressurisation +# Copyright (c) 2021-2025 Anders Andreasen +# Published under an MIT license + +""" +Custom exceptions for HydDown package. + +This module defines a hierarchy of exceptions for different error conditions +that can occur during pressure vessel simulations. Using specific exception +types allows calling code to handle different error cases appropriately. + +Exception Hierarchy: + HydDownError (base) + ├── ThermodynamicError + │ ├── ThermodynamicConvergenceError + │ ├── InvalidStateError + │ └── PhaseEquilibriumError + ├── NumericalError + │ ├── NumericalInstabilityError + │ ├── IntegrationError + │ └── ConvergenceError + ├── ConfigurationError + │ ├── ValveConfigurationError + │ ├── VesselConfigurationError + │ └── HeatTransferConfigurationError + └── PhysicalConstraintError + ├── TriplePointViolation + ├── CriticalPointViolation + └── NegativeMassError + +Usage: + from hyddown.exceptions import ThermodynamicConvergenceError + + if not result.success: + raise ThermodynamicConvergenceError( + f"Failed to converge: {result.message}" + ) +""" + + +class HydDownError(Exception): + """ + Base exception for all HydDown errors. + + All custom exceptions in HydDown inherit from this base class, + allowing calling code to catch all HydDown-specific errors with + a single except clause if desired. + """ + pass + + +# ============================================================================ +# Thermodynamic Errors +# ============================================================================ + +class ThermodynamicError(HydDownError): + """ + Base class for thermodynamic calculation errors. + + Raised when CoolProp property calculations fail or produce + non-physical results. + """ + pass + + +class ThermodynamicConvergenceError(ThermodynamicError): + """ + Failed to converge in thermodynamic property solver. + + Raised when numerical optimization fails to find a valid thermodynamic + state (e.g., during PHproblem or UDproblem calculations). This typically + occurs with multicomponent mixtures when scipy.optimize cannot find + a solution. + + Attributes + ---------- + solver : str + Name of the solver that failed ('PHproblem', 'UDproblem', etc.) + state_vars : dict + Input state variables that caused failure + iterations : int, optional + Number of iterations attempted + """ + + def __init__(self, message, solver=None, state_vars=None, iterations=None): + super().__init__(message) + self.solver = solver + self.state_vars = state_vars + self.iterations = iterations + + +class InvalidStateError(ThermodynamicError): + """ + Thermodynamic state is outside valid range for equation of state. + + Raised when requested state conditions are outside the valid range + of the CoolProp equation of state (e.g., below triple point, above + maximum temperature, or at unphysical densities). + + Attributes + ---------- + variable : str + The state variable that is out of range ('temperature', 'pressure', etc.) + value : float + The invalid value + valid_range : tuple + (min, max) valid range for the variable + """ + + def __init__(self, message, variable=None, value=None, valid_range=None): + super().__init__(message) + self.variable = variable + self.value = value + self.valid_range = valid_range + + +class PhaseEquilibriumError(ThermodynamicError): + """ + Error in phase equilibrium calculations for two-phase systems. + + Raised when calculations involving vapor-liquid equilibrium fail, + such as when trying to determine liquid level or saturation properties. + """ + pass + + +# ============================================================================ +# Numerical Errors +# ============================================================================ + +class NumericalError(HydDownError): + """ + Base class for numerical integration and solver errors. + """ + pass + + +class NumericalInstabilityError(NumericalError): + """ + Time step too large for numerical stability. + + Raised when the explicit Euler integration scheme becomes unstable, + typically because the time step is too large relative to the + characteristic time scales of the problem. This can manifest as + negative masses, temperatures, or pressures. + + Attributes + ---------- + time_step : float + The time step that caused instability [s] + recommended_dt : float, optional + Recommended smaller time step [s] + characteristic_time : float, optional + Characteristic time scale of the problem [s] + """ + + def __init__(self, message, time_step=None, recommended_dt=None, + characteristic_time=None): + super().__init__(message) + self.time_step = time_step + self.recommended_dt = recommended_dt + self.characteristic_time = characteristic_time + + +class IntegrationError(NumericalError): + """ + Error during time integration of mass/energy balances. + + Raised when the numerical integration produces invalid results, + such as NaN or Inf values in state variables. + """ + pass + + +class ConvergenceError(NumericalError): + """ + General convergence failure in iterative solvers. + + Raised when iterative algorithms (not thermodynamic solvers) + fail to converge within specified tolerance or iteration limit. + """ + pass + + +# ============================================================================ +# Configuration Errors +# ============================================================================ + +class ConfigurationError(HydDownError): + """ + Base class for invalid input configuration errors. + + These errors indicate problems with the input YAML file beyond + what the Cerberus schema validator can catch (e.g., physically + impossible values or inconsistent parameters). + """ + pass + + +class ValveConfigurationError(ConfigurationError): + """ + Invalid valve parameters. + + Raised when valve configuration is physically impossible or + inconsistent (e.g., negative diameter, discharge coefficient > 1, + relief valve set pressure below operating pressure). + + Attributes + ---------- + valve_type : str + Type of valve ('orifice', 'relief_valve', 'control_valve') + parameter : str + The problematic parameter name + value : float + The invalid value + """ + + def __init__(self, message, valve_type=None, parameter=None, value=None): + super().__init__(message) + self.valve_type = valve_type + self.parameter = parameter + self.value = value + + +class VesselConfigurationError(ConfigurationError): + """ + Invalid vessel geometry or material parameters. + + Raised when vessel configuration is physically impossible + (e.g., negative dimensions, wall thickness greater than diameter, + invalid material properties). + """ + pass + + +class HeatTransferConfigurationError(ConfigurationError): + """ + Invalid heat transfer parameters. + + Raised when heat transfer configuration is inconsistent or + physically impossible (e.g., negative heat transfer coefficient, + invalid fire scenario). + """ + pass + + +# ============================================================================ +# Physical Constraint Errors +# ============================================================================ + +class PhysicalConstraintError(HydDownError): + """ + Base class for violations of physical constraints. + + These errors indicate that the simulation has reached a state + that violates fundamental physical constraints. + """ + pass + + +class TriplePointViolation(PhysicalConstraintError): + """ + Fluid state below triple point temperature or pressure. + + Raised when calculations result in conditions below the fluid's + triple point, where the equation of state is not valid and + solid phase would form. + + Attributes + ---------- + fluid : str + Name of the fluid + temperature : float + Current temperature [K] + pressure : float + Current pressure [Pa] + T_triple : float + Triple point temperature [K] + P_triple : float + Triple point pressure [Pa] + """ + + def __init__(self, message, fluid=None, temperature=None, pressure=None, + T_triple=None, P_triple=None): + super().__init__(message) + self.fluid = fluid + self.temperature = temperature + self.pressure = pressure + self.T_triple = T_triple + self.P_triple = P_triple + + +class CriticalPointViolation(PhysicalConstraintError): + """ + Fluid state exceeds equation of state validity limits. + + Raised when calculations produce conditions far beyond the + critical point or other EOS validity limits. + """ + pass + + +class NegativeMassError(PhysicalConstraintError): + """ + Mass, pressure, or temperature became negative. + + Raised when numerical integration produces unphysical negative + values for extensive properties. This usually indicates numerical + instability or time step too large. + + Attributes + ---------- + variable : str + The variable that became negative ('mass', 'pressure', 'temperature') + value : float + The negative value + time : float + Simulation time when error occurred [s] + """ + + def __init__(self, message, variable=None, value=None, time=None): + super().__init__(message) + self.variable = variable + self.value = value + self.time = time + + +# ============================================================================ +# Warning Classes (for non-fatal issues) +# ============================================================================ + +class HydDownWarning(UserWarning): + """ + Base warning class for HydDown. + + Warnings are issued for conditions that may affect accuracy but + don't prevent calculation from completing. + """ + pass + + +class NumericalStabilityWarning(HydDownWarning): + """ + Warning that time step may be too large for optimal stability. + + Issued when CFL condition or other stability criteria suggest + the time step should be reduced, but calculation can continue. + """ + pass + + +class AccuracyWarning(HydDownWarning): + """ + Warning that results may have reduced accuracy. + + Issued when approximations or numerical issues may affect + accuracy but not prevent calculation completion. + """ + pass + + +class ConvergenceWarning(HydDownWarning): + """ + Warning about slow or marginal convergence. + + Issued when iterative solvers converge but required many + iterations or achieved marginal tolerance. + """ + pass diff --git a/src/hyddown/hdclass.py b/src/hyddown/hdclass.py index 343ba87..61bcf36 100644 --- a/src/hyddown/hdclass.py +++ b/src/hyddown/hdclass.py @@ -70,6 +70,15 @@ from hyddown import validator from hyddown import fire from hyddown import thermesh as tm +from hyddown import safety_checks as sc +from hyddown.exceptions import ( + ThermodynamicConvergenceError, + NegativeMassError, + TriplePointViolation, + NumericalStabilityWarning, +) +from hyddown.thermo_solver import ThermodynamicSolver +import warnings import fluids @@ -297,22 +306,19 @@ def initialize(self): self.surf_area_outer = self.outer_vol.A self.surf_area_inner = self.inner_vol.A - self.fluid = CP.AbstractState("HEOS", self.comp) - if "&" in self.comp: - self.fluid.specify_phase(CP.iphase_gas) - self.fluid.set_mole_fractions(self.molefracs) - - self.transport_fluid = CP.AbstractState("HEOS", self.compSRK) - self.transport_fluid.specify_phase(CP.iphase_gas) - self.transport_fluid.set_mole_fractions(self.molefracs) - - self.transport_fluid_wet = CP.AbstractState("HEOS", self.compSRK) - self.transport_fluid_wet.specify_phase(CP.iphase_liquid) - self.transport_fluid_wet.set_mole_fractions(self.molefracs) + # Initialize thermodynamic solver + self.solver = ThermodynamicSolver( + species=self.comp, + mole_fractions=self.molefracs, + vessel_geometry=self.inner_vol, + species_SRK=self.compSRK, + ) - self.vent_fluid = CP.AbstractState("HEOS", self.comp) - self.vent_fluid.specify_phase(CP.iphase_gas) - self.vent_fluid.set_mole_fractions(self.molefracs) + # Create convenience references to solver's fluid states + self.fluid = self.solver.fluid + self.transport_fluid = self.solver.transport_fluid + self.transport_fluid_wet = self.solver.transport_fluid_wet + self.vent_fluid = self.solver.vent_fluid if "liquid_level" in self.input["vessel"]: ll = self.input["vessel"]["liquid_level"] @@ -340,8 +346,8 @@ def initialize(self): self.Q0 = self.fluid.Q() self.vent_fluid.update(CP.PT_INPUTS, self.p0, self.T0) - self.res_fluid = CP.AbstractState("HEOS", self.comp) - self.res_fluid.set_mole_fractions(self.molefracs) + # Create convenience reference to solver's reservoir fluid + self.res_fluid = self.solver.res_fluid if self.input["valve"]["flow"] == "filling": self.res_fluid.update(CP.PT_INPUTS, self.p_back, self.T0) @@ -388,186 +394,75 @@ def initialize(self): self.vapour_volume_fraction = np.zeros(data_len) self.liquid_level = np.zeros(data_len) + # Convenience methods for backward compatibility + # These delegate to the ThermodynamicSolver def calc_liquid_level(self): """ Calculate liquid level height based on current two-phase fluid state. - For two-phase systems (0 ≤ quality ≤ 1), calculates the height of liquid - phase in the vessel based on vapor quality, phase densities, and vessel geometry. - Uses vessel geometry from fluids.TANK to convert liquid volume to height. - - Parameters - ---------- - fluid : CoolProp AbstractState - Current fluid state + This is a convenience method that delegates to the ThermodynamicSolver. + Maintained for backward compatibility with existing code. Returns ------- float - Liquid level height from vessel bottom [m]. - Returns 0.0 for single-phase gas (quality > 1 or subcooled liquid). - """ - if self.fluid.Q() >= 0 and self.fluid.Q() <= 1: - rho_liq = self.fluid.saturated_liquid_keyed_output(CP.iDmass) - rho_vap = self.fluid.saturated_vapor_keyed_output(CP.iDmass) - m_liq = self.fluid.rhomass() * self.inner_vol.V_total * (1 - self.fluid.Q()) - V_liq = m_liq / rho_liq - h_liq = self.inner_vol.h_from_V(V_liq) - return h_liq - else: - return 0.0 - - def PHres(self, T, P, H): + Liquid level height from vessel bottom [m] """ - Residual enthalpy function to be minimized during PH-problem. - - Used by numerical optimizer (scipy.optimize.minimize) to find temperature - that satisfies constant pressure-enthalpy constraints for multicomponent - fluids. The optimizer adjusts T until residual approaches zero. - - Parameters - ---------- - H : float - Enthalpy at initial/final conditions [J/kg] - P : float - Pressure at final conditions [Pa] - T : float - Updated estimate for the final temperature at (P,H) [K] + return self.solver.calc_liquid_level() - Returns - ------- - float - Squared normalized enthalpy residual (dimensionless). - Zero when correct temperature is found. - """ - # Extract scalar from array (scipy optimizers pass arrays, CoolProp needs scalars) - T_scalar = float(T.item()) if hasattr(T, 'item') else float(T) - self.vent_fluid.update(CP.PT_INPUTS, P, T_scalar) - return ((H - self.vent_fluid.hmass()) / H) ** 2 - - def PHres_relief(self, T, P, H): + def PHproblem(self, H, P, Tguess, relief=False): """ - Residual enthalpy function for PH-problem during relief valve calculations. + Solve constant pressure, constant enthalpy problem. - Used by numerical optimizer (scipy.optimize.root_scalar) to find temperature - for relief valve discharge calculations. Similar to PHres() but uses the - main fluid state instead of vent_fluid state. + This is a convenience method that delegates to the ThermodynamicSolver. + Maintained for backward compatibility with existing code. Parameters ---------- H : float - Enthalpy at initial/final conditions [J/kg] + Enthalpy [J/kg] P : float - Pressure at final conditions [Pa] - T : float - Updated estimate for the final temperature at (P,H) [K] + Pressure [Pa] + Tguess : float + Initial temperature guess [K] + relief : bool, optional + Use relief valve solver if True Returns ------- float - Normalized enthalpy residual (dimensionless). - Zero when correct temperature is found. - """ - # Extract scalar from array (scipy optimizers pass arrays, CoolProp needs scalars) - T_scalar = float(T.item()) if hasattr(T, 'item') else float(T) - self.fluid.update(CP.PT_INPUTS, P, T_scalar) - return (H - self.fluid.hmass()) / H - - def PHproblem(self, H, P, Tguess, relief=False): - """ - Defining a constant pressure, constant enthalpy problem i.e. typical adiabatic - problem like e.g. valve flow for the vented flow (during discharge). - For multicomponent mixture the final temperature is changed/optimised until the residual - enthalpy is near zero in an optimisation step. For single component fluid the coolprop - built in methods are used for speed. - - Parameters - ---------- - H : float - Enthalpy at initial/final conditions - P : float - Pressure at final conditions. - Tguess : float - Initial guess for the final temperature at P,H - """ - - # Multicomponent case - if "&" in self.species: - x0 = Tguess - if relief == False: - res = minimize( - self.PHres, - x0, - args=(P, H), - method="Nelder-Mead", - options={"xatol": 0.1, "fatol": 0.001}, - ) - T1 = res.x[0] - else: - res = root_scalar( - self.PHres_relief, - args=(P, H), - x0=x0, - method="newton", - ) - T1 = res.root - # single component fluid case - else: - T1 = PropsSI("T", "P", P, "H", H, self.species) - return T1 - - def UDres(self, x, U, rho): - """ - Residual U-rho to be minimised during a U-rho/UV-problem - - Parameters - ---------- - U : float - Internal energy at final conditions - rho : float - Density at final conditions + Temperature at (P, H) [K] """ - self.fluid.update(CP.PT_INPUTS, x[0], x[1]) - return ((U - self.fluid.umass()) / U) ** 2 + ( - (rho - self.fluid.rhomass()) / rho - ) ** 2 + return self.solver.PHproblem(H, P, Tguess, relief) def UDproblem(self, U, rho, Pguess, Tguess): """ - Defining a constant UV problem i.e. constant internal energy and density/volume - problem relevant for the 1. law of thermodynamics. - For multicomponent mixture the final temperature/pressure is changed/optimised until the - residual U/rho is near zero. For single component fluid the coolprop - built in methods are used for speed. + Solve constant internal energy, constant density problem. + + This is a convenience method that delegates to the ThermodynamicSolver. + Maintained for backward compatibility with existing code. Parameters ---------- U : float - Internal energy at final conditions + Internal energy [J/kg] rho : float - Density at final conditions. + Density [kg/m³] Pguess : float - Initial guess for the final pressure at U, rho + Initial pressure guess [Pa] Tguess : float - Initial guess for the final temperature at U, rho + Initial temperature guess [K] + + Returns + ------- + P1 : float + Pressure [Pa] + T1 : float + Temperature [K] + Ures : float + Internal energy residual [J/kg] """ - if "&" in self.species: - x0 = [Pguess, Tguess] - res = minimize( - self.UDres, - x0, - args=(U, rho), - method="Nelder-Mead", - options={"xatol": 0.1, "fatol": 0.001}, - ) - P1 = res.x[0] - T1 = res.x[1] - Ures = U - self.fluid.umass() - else: - P1 = PropsSI("P", "D", rho, "U", U, self.species) - T1 = PropsSI("T", "D", rho, "U", U, self.species) - Ures = 0 - return P1, T1, Ures + return self.solver.UDproblem(U, rho, Pguess, Tguess) def run(self, disable_pbar=True): """ @@ -730,6 +625,15 @@ def run(self, disable_pbar=True): self.rho[i] = self.mass_fluid[i] / self.vol + # Periodically check CFL stability condition (every 100 steps) + if i % 100 == 0 and i > 0: + sc.check_cfl_stability( + mass_flow_rate=self.mass_rate[i - 1], + vessel_mass=self.mass_fluid[i], + time_step=self.tstep, + characteristic_fraction=0.1 + ) + # ------------------------------------------------------------------------ # THERMODYNAMIC STATE UPDATE # ------------------------------------------------------------------------ @@ -1117,14 +1021,24 @@ def run(self, disable_pbar=True): domain2_w, solver2_w ) else: + # Boundary conditions for unwetted (gas) section + # Use safe_divide to handle fully liquid case (unwetted_area = 0) bc = [ { - "q": -self.Q_inner[i] - / (self.surf_area_inner - wetted_area) + "q": sc.safe_divide( + -self.Q_inner[i], + self.surf_area_inner - wetted_area, + name="q_inner_unwetted", + default=0.0 + ) }, { - "q": self.Q_outer[i] - / (self.surf_area_outer - wetted_area) + "q": sc.safe_divide( + self.Q_outer[i], + self.surf_area_outer - wetted_area, + name="q_outer_unwetted", + default=0.0 + ) }, ] domain2 = tm.Domain(mesh2, [liner, shell], bc) @@ -1135,9 +1049,25 @@ def run(self, disable_pbar=True): "theta": theta, } t_bonded, T_profile2 = tm.solve_ht(domain2, solver2) + # Boundary conditions for wetted (liquid) section + # Use safe_divide to handle fully vapor case (wetted_area = 0) bc_w = [ - {"q": -self.Q_inner_wetted[i] / (wetted_area)}, - {"q": self.Q_outer_wetted[i] / wetted_area}, + { + "q": sc.safe_divide( + -self.Q_inner_wetted[i], + wetted_area, + name="q_inner_wetted", + default=0.0 + ) + }, + { + "q": sc.safe_divide( + self.Q_outer_wetted[i], + wetted_area, + name="q_outer_wetted", + default=0.0 + ) + }, ] domain2_w = tm.Domain(mesh2_w, [liner_w, shell_w], bc_w) domain2_w.set_T(T_profile2_w[-1, :]) @@ -1391,7 +1321,7 @@ def run(self, disable_pbar=True): 1, ) else: - T1 = self.PHproblem( + T1 = self.solver.PHproblem( h_in + self.tstep * self.Q_inner[i] / self.mass_fluid[i], self.Pset, @@ -1437,7 +1367,7 @@ def run(self, disable_pbar=True): else: self.mass_rate[i] = 0 - P1, T1, self.U_res[i] = self.UDproblem( + P1, T1, self.U_res[i] = self.solver.UDproblem( U_end / self.mass_fluid[i], self.rho[i], self.P[i - 1], @@ -1449,7 +1379,7 @@ def run(self, disable_pbar=True): self.fluid.update(CP.PT_INPUTS, self.P[i], self.T_fluid[i]) else: - P1, T1, self.U_res[i] = self.UDproblem( + P1, T1, self.U_res[i] = self.solver.UDproblem( U_end / self.mass_fluid[i], self.rho[i], self.P[i - 1], @@ -1495,12 +1425,12 @@ def run(self, disable_pbar=True): self.S_mass[i] = self.fluid.smass() self.U_mass[i] = self.fluid.umass() - self.liquid_level[i] = self.calc_liquid_level() + self.liquid_level[i] = self.solver.calc_liquid_level() # Calculating vent temperature (adiabatic) only for discharge problem if self.input["valve"]["flow"] == "discharge": if "&" in self.species: - self.T_vent[i] = self.PHproblem( + self.T_vent[i] = self.solver.PHproblem( self.H_mass[i], self.p_back, self.vent_fluid.T() ) else: @@ -1603,7 +1533,7 @@ def run(self, disable_pbar=True): if input["valve"]["type"] == "relief": idx_max = self.mass_rate.argmax() # Smooth peak value by averaging with neighbors (avoid array index out of bounds) - if 0 < idx_max < len(self.mass_rate) - 1: + if sc.check_bounds_for_smoothing(idx_max, len(self.mass_rate), "relief valve smoothing"): self.mass_rate[idx_max] = ( self.mass_rate[idx_max - 1] + self.mass_rate[idx_max + 1] ) / 2 diff --git a/src/hyddown/safety_checks.py b/src/hyddown/safety_checks.py new file mode 100644 index 0000000..1d993e8 --- /dev/null +++ b/src/hyddown/safety_checks.py @@ -0,0 +1,464 @@ +# HydDown hydrogen/other gas depressurisation +# Copyright (c) 2021-2025 Anders Andreasen +# Published under an MIT license + +""" +Runtime safety checks for HydDown simulations. + +This module provides defensive checks for issues that can ONLY be detected +during simulation execution, not from static input validation. Cerberus +validator handles all input file validation - this module handles runtime +numerical and physical errors. + +Runtime-only checks: +- Negative values from numerical errors (mass, pressure, temperature) +- NaN/Inf detection during calculation +- Array bounds checking +- CFL stability condition (depends on computed flow rates) +- Triple point violations during discharge +- Thermodynamic solver convergence +- Division by zero protection +""" + +import warnings +import numpy as np +from CoolProp.CoolProp import PropsSI +from hyddown.exceptions import ( + NegativeMassError, + TriplePointViolation, + InvalidStateError, + NumericalInstabilityError, + NumericalStabilityWarning, + ThermodynamicConvergenceError, + ConvergenceWarning, +) + + +def check_positive_runtime(value, name, time=None, step=None): + """ + Check that a computed value is positive during simulation. + + This checks for negative values that arise from numerical errors during + time integration, NOT from invalid input (Cerberus handles input validation). + + Parameters + ---------- + value : float + Computed value to check + name : str + Name of the variable (for error messages) + time : float, optional + Simulation time when check occurs [s] + step : int, optional + Time step number when check occurs + + Raises + ------ + NegativeMassError + If computed value is negative or zero + + Examples + -------- + >>> check_positive_runtime(100000, "pressure", time=5.0, step=100) + >>> check_positive_runtime(-50, "mass", time=10.0) # Raises error + """ + if value <= 0: + msg = f"{name} became non-positive during simulation: {value:.3e}" + if time is not None: + msg += f" at t={time:.3f}s" + if step is not None: + msg += f" (step {step})" + msg += ". This indicates numerical instability - try reducing time step." + + raise NegativeMassError( + msg, + variable=name, + value=value, + time=time + ) + + +def check_array_bounds(index, array_length, context=""): + """ + Check array index is within bounds during simulation. + + This prevents IndexError crashes from array access bugs. + Particularly important for relief valve smoothing and similar operations. + + Parameters + ---------- + index : int + Index to check + array_length : int + Length of the array + context : str, optional + Description of operation (for error messages) + + Raises + ------ + IndexError + If index is out of bounds + + Examples + -------- + >>> check_array_bounds(5, 10, "mass_rate smoothing") + >>> check_array_bounds(10, 10) # Raises IndexError + """ + if index < 0 or index >= array_length: + msg = f"Index {index} out of bounds for array of length {array_length}" + if context: + msg += f" in {context}" + msg += f". Valid range is [0, {array_length-1}]" + raise IndexError(msg) + + +def check_nan_inf(value, name, time=None, step=None): + """ + Check for NaN or Inf in computed values. + + NaN/Inf indicate serious numerical problems (division by zero, + overflow, etc.) and must be caught immediately. + + Parameters + ---------- + value : float or np.ndarray + Computed value(s) to check + name : str + Name of the variable + time : float, optional + Simulation time [s] + step : int, optional + Time step number + + Raises + ------ + ValueError + If value contains NaN or Inf + + Examples + -------- + >>> check_nan_inf(5.0, "temperature") + >>> check_nan_inf(np.nan, "pressure", time=10.0) # Raises ValueError + """ + has_nan = np.any(np.isnan(value)) + has_inf = np.any(np.isinf(value)) + + if has_nan or has_inf: + problem = "NaN" if has_nan else "Inf" + msg = f"{name} contains {problem}: {value}" + if time is not None: + msg += f" at t={time:.3f}s" + if step is not None: + msg += f" (step {step})" + msg += ". This indicates numerical instability." + raise ValueError(msg) + + +def check_triple_point_runtime(temperature, pressure, fluid_name, time=None): + """ + Check if state is above triple point during simulation. + + Triple point violations can occur during discharge even if initial + conditions are valid. This is a physical constraint violation. + + Parameters + ---------- + temperature : float + Current temperature [K] + pressure : float + Current pressure [Pa] + fluid_name : str + CoolProp fluid name (e.g., "CO2", "N2") + time : float, optional + Simulation time [s] + + Raises + ------ + TriplePointViolation + If temperature or pressure is below triple point + + Warnings + -------- + Issues warning if within 5% of triple point + + Examples + -------- + >>> check_triple_point_runtime(300, 101325, "N2") # OK + >>> check_triple_point_runtime(200, 101325, "CO2") # Raises error + """ + try: + # Get triple point properties (not all fluids have this data) + T_triple = PropsSI("Ttriple", fluid_name) + P_triple = PropsSI("ptriple", fluid_name) + + # Check temperature + if temperature < T_triple: + msg = (f"{fluid_name} temperature {temperature:.2f}K dropped below " + f"triple point {T_triple:.2f}K during simulation") + if time is not None: + msg += f" at t={time:.3f}s" + msg += ". Solid phase would form - EOS invalid." + + raise TriplePointViolation( + msg, + fluid=fluid_name, + temperature=temperature, + pressure=pressure, + T_triple=T_triple, + P_triple=P_triple + ) + + # Check pressure + if pressure < P_triple: + msg = (f"{fluid_name} pressure {pressure:.2e}Pa dropped below " + f"triple point {P_triple:.2e}Pa during simulation") + if time is not None: + msg += f" at t={time:.3f}s" + msg += ". Solid phase would form - EOS invalid." + + raise TriplePointViolation( + msg, + fluid=fluid_name, + temperature=temperature, + pressure=pressure, + T_triple=T_triple, + P_triple=P_triple + ) + + # Warning if close to triple point (within 5%) + if temperature < T_triple * 1.05: + msg = (f"{fluid_name} temperature {temperature:.2f}K is within 5% " + f"of triple point {T_triple:.2f}K") + if time is not None: + msg += f" at t={time:.3f}s" + msg += ". Results may be inaccurate near phase boundary." + warnings.warn(msg, NumericalStabilityWarning) + + except ValueError: + # Fluid doesn't have triple point data in CoolProp, skip check + pass + + +def check_cfl_stability(mass_flow_rate, vessel_mass, time_step, + characteristic_fraction=0.1): + """ + Check CFL (Courant-Friedrichs-Lewy) stability condition. + + For explicit Euler integration, the time step must be small compared + to the characteristic time scale (mass inventory / mass flow rate). + This check can only be done at runtime with actual computed flow rates. + + Parameters + ---------- + mass_flow_rate : float + Current mass flow rate [kg/s] (absolute value) + vessel_mass : float + Current mass in vessel [kg] + time_step : float + Integration time step [s] + characteristic_fraction : float, default 0.1 + Fraction of characteristic time for stability (typically 0.05 to 0.2) + + Warnings + -------- + Issues NumericalStabilityWarning if time step is too large + + Returns + ------- + float + Recommended maximum time step [s], or None if no limit + + Examples + -------- + >>> check_cfl_stability(0.1, 50.0, 1.0) # May warn + >>> check_cfl_stability(0.001, 50.0, 0.01) # OK + """ + # Avoid division by zero for very small flow rates + if abs(mass_flow_rate) < 1e-10: + return None # No significant flow, any time step is OK + + # Characteristic time: how long to change vessel mass significantly + # tau = mass / dmdt + characteristic_time = vessel_mass / abs(mass_flow_rate) + + # Recommended maximum time step + dt_max_recommended = characteristic_fraction * characteristic_time + + # Warn if current time step is too large + if time_step > dt_max_recommended: + warnings.warn( + f"Time step dt={time_step:.3f}s may be too large for stability. " + f"Recommended: dt < {dt_max_recommended:.3f}s " + f"(based on characteristic time {characteristic_time:.3f}s). " + f"Current mass flow rate = {abs(mass_flow_rate):.3e} kg/s, " + f"vessel mass = {vessel_mass:.3e} kg. " + f"Consider reducing time step if results show oscillations.", + NumericalStabilityWarning, + stacklevel=2 + ) + + return dt_max_recommended + + +def check_optimization_convergence(result, solver_name, state_vars=None, + tolerance=1e-4): + """ + Check if scipy optimization converged successfully. + + This checks convergence of thermodynamic solvers (PHproblem, UDproblem) + which use numerical optimization for multicomponent fluids. Convergence + can only be checked at runtime. + + Parameters + ---------- + result : scipy.optimize.OptimizeResult + Result object from scipy.optimize.minimize or root_scalar + solver_name : str + Name of the solver (for error messages) + state_vars : dict, optional + Input state variables being solved for + tolerance : float, default 1e-4 + Acceptable residual tolerance + + Raises + ------ + ThermodynamicConvergenceError + If optimization did not converge + + Warnings + -------- + Issues ConvergenceWarning if convergence was marginal + + Examples + -------- + >>> from scipy.optimize import minimize + >>> result = minimize(lambda x: x**2, x0=1.0) + >>> check_optimization_convergence(result, "test_solver") + """ + # Check if optimization succeeded + if hasattr(result, 'success'): + if not result.success: + msg = f"{solver_name} failed to converge: {result.message}" + if state_vars: + msg += f"\nInput state: {state_vars}" + raise ThermodynamicConvergenceError( + msg, + solver=solver_name, + state_vars=state_vars, + iterations=getattr(result, 'nit', None) + ) + + # Check converged flag (for root_scalar) + if hasattr(result, 'converged'): + if not result.converged: + msg = f"{solver_name} did not converge" + if state_vars: + msg += f" for state: {state_vars}" + raise ThermodynamicConvergenceError( + msg, + solver=solver_name, + state_vars=state_vars, + iterations=getattr(result, 'iterations', None) + ) + + # Warn if residual is large (marginal convergence) + if hasattr(result, 'fun'): + if isinstance(result.fun, (list, np.ndarray)): + residual = np.max(np.abs(result.fun)) + else: + residual = abs(result.fun) + + if residual > tolerance: + warnings.warn( + f"{solver_name} converged but residual {residual:.2e} " + f"exceeds tolerance {tolerance:.2e}. Results may be inaccurate.", + ConvergenceWarning, + stacklevel=2 + ) + + # Warn if required many iterations + nit = getattr(result, 'nit', getattr(result, 'iterations', None)) + if nit is not None and nit > 100: + warnings.warn( + f"{solver_name} required {nit} iterations. " + "Consider improving initial guess for better performance.", + ConvergenceWarning, + stacklevel=2 + ) + + +def safe_divide(numerator, denominator, name="division", default=0.0): + """ + Safely divide two numbers, returning default if denominator is zero. + + This prevents division by zero errors that occur at runtime (e.g., + when wetted_area = 0 in two-phase calculations). + + Parameters + ---------- + numerator : float + Numerator value + denominator : float + Denominator value + name : str, optional + Description of the division operation (for warnings) + default : float, default 0.0 + Value to return if denominator is zero or very small + + Returns + ------- + float + numerator / denominator, or default if division would fail + + Examples + -------- + >>> safe_divide(10.0, 2.0) + 5.0 + >>> safe_divide(10.0, 0.0, default=0.0) + 0.0 + >>> safe_divide(10.0, 1e-20, name="heat flux") + 0.0 + """ + # Check for zero or near-zero denominator + if abs(denominator) < 1e-15: + # Silently return default - this is expected behavior + # (e.g., wetted_area = 0 is normal for fully vapor state) + return default + + return numerator / denominator + + +def check_bounds_for_smoothing(idx, array_length, context=""): + """ + Check that index and neighbors are within bounds for smoothing operations. + + Relief valve smoothing and similar operations need idx-1, idx, idx+1. + This checks all three indices are valid. + + Parameters + ---------- + idx : int + Center index to smooth + array_length : int + Array length + context : str, optional + Description of operation + + Returns + ------- + bool + True if smoothing is safe (idx-1, idx, idx+1 all valid) + + Examples + -------- + >>> check_bounds_for_smoothing(5, 10) + True + >>> check_bounds_for_smoothing(0, 10) # Can't access idx-1 + False + >>> check_bounds_for_smoothing(9, 10) # Can't access idx+1 + False + """ + # Need idx-1, idx, and idx+1 to be valid + if idx < 1 or idx >= array_length - 1: + return False + return True diff --git a/src/hyddown/thermo_solver.py b/src/hyddown/thermo_solver.py new file mode 100644 index 0000000..767ee66 --- /dev/null +++ b/src/hyddown/thermo_solver.py @@ -0,0 +1,397 @@ +# HydDown hydrogen/other gas depressurisation +# Copyright (c) 2021-2025 Anders Andreasen +# Published under an MIT license + +""" +Thermodynamic solver for CoolProp fluid property calculations. + +This module provides the ThermodynamicSolver class which encapsulates all +thermodynamic state calculations for HydDown simulations. It manages CoolProp +AbstractState objects and provides methods for solving thermodynamic problems +(PH-problems, UD-problems) for both single-component and multicomponent fluids. + +The ThermodynamicSolver handles: +- CoolProp AbstractState initialization and management +- PH-problems: Finding temperature at constant pressure and enthalpy +- UD-problems: Finding pressure and temperature at constant internal energy and density +- Liquid level calculations for two-phase systems +- Single-component (fast, direct CoolProp calls) and multicomponent (slower, numerical optimization) fluids + +Key Methods: +- PHproblem(): Solve for temperature at given P, H (isenthalpic expansion) +- UDproblem(): Solve for P, T at given U, density (constant internal energy) +- calc_liquid_level(): Calculate liquid height in two-phase systems + +Typical usage: + solver = ThermodynamicSolver( + species="HEOS::Hydrogen", + mole_fractions=[1.0], + vessel_geometry=inner_vol + ) + T = solver.PHproblem(H=1e6, P=1e5, Tguess=300) +""" + +import numpy as np +from scipy.optimize import minimize, root_scalar +from CoolProp.CoolProp import PropsSI +import CoolProp.CoolProp as CP +from hyddown import safety_checks as sc + + +class ThermodynamicSolver: + """ + Thermodynamic solver for CoolProp-based fluid property calculations. + + This class encapsulates all thermodynamic calculations required for HydDown + simulations, managing CoolProp AbstractState objects and providing methods + to solve for thermodynamic states given various property pairs. + + Parameters + ---------- + species : str + CoolProp fluid name (e.g., "HEOS::Hydrogen") or mixture + (e.g., "HEOS::Methane[0.9]&Ethane[0.1]") + mole_fractions : list of float + Mole fractions for each component (must sum to 1.0) + vessel_geometry : fluids.TANK, optional + Vessel geometry object for liquid level calculations. + Required if calc_liquid_level() will be used. + species_SRK : str, optional + CoolProp species string using SRK backend for transport properties. + If not provided, defaults to species. + + Attributes + ---------- + fluid : CoolProp.AbstractState + Main fluid state for thermodynamic calculations + vent_fluid : CoolProp.AbstractState + Fluid state for vent/discharge calculations (gas phase) + res_fluid : CoolProp.AbstractState + Fluid state for reservoir/filling calculations + transport_fluid : CoolProp.AbstractState + Fluid state for transport property calculations (gas phase) + transport_fluid_wet : CoolProp.AbstractState + Fluid state for wetted region transport properties (liquid phase) + is_multicomponent : bool + True if fluid is multicomponent mixture + """ + + def __init__(self, species, mole_fractions, vessel_geometry=None, species_SRK=None): + """ + Initialize the thermodynamic solver with fluid composition and geometry. + + Parameters + ---------- + species : str + CoolProp fluid name + mole_fractions : list of float + Mole fractions for each component + vessel_geometry : fluids.TANK, optional + Vessel geometry for liquid level calculations + species_SRK : str, optional + Species string for SRK backend (transport properties) + """ + self.species = species + self.mole_fractions = mole_fractions + self.vessel_geometry = vessel_geometry + self.is_multicomponent = "&" in species + + # Set up SRK species for transport properties + if species_SRK is None: + self.species_SRK = species + else: + self.species_SRK = species_SRK + + # Initialize CoolProp AbstractState objects + self._initialize_fluid_states() + + def _initialize_fluid_states(self): + """ + Create CoolProp AbstractState objects for different calculation needs. + + Creates: + - fluid: Main thermodynamic state + - vent_fluid: Vent/discharge calculations (gas phase) + - res_fluid: Reservoir/filling calculations + - transport_fluid: Transport properties (gas phase, SRK backend) + - transport_fluid_wet: Transport properties (liquid phase, SRK backend) + """ + # Main fluid state + self.fluid = CP.AbstractState("HEOS", self.species) + if self.is_multicomponent: + self.fluid.specify_phase(CP.iphase_gas) + self.fluid.set_mole_fractions(self.mole_fractions) + + # Vent fluid (for discharge calculations, gas phase) + self.vent_fluid = CP.AbstractState("HEOS", self.species) + self.vent_fluid.specify_phase(CP.iphase_gas) + self.vent_fluid.set_mole_fractions(self.mole_fractions) + + # Reservoir fluid (for filling calculations) + self.res_fluid = CP.AbstractState("HEOS", self.species) + self.res_fluid.set_mole_fractions(self.mole_fractions) + + # Transport properties (gas phase, SRK backend for robustness) + self.transport_fluid = CP.AbstractState("HEOS", self.species_SRK) + self.transport_fluid.specify_phase(CP.iphase_gas) + self.transport_fluid.set_mole_fractions(self.mole_fractions) + + # Transport properties (liquid phase for wetted regions) + self.transport_fluid_wet = CP.AbstractState("HEOS", self.species_SRK) + self.transport_fluid_wet.specify_phase(CP.iphase_liquid) + self.transport_fluid_wet.set_mole_fractions(self.mole_fractions) + + def calc_liquid_level(self): + """ + Calculate liquid level height based on current two-phase fluid state. + + For two-phase systems (0 ≤ quality ≤ 1), calculates the height of liquid + phase in the vessel based on vapor quality, phase densities, and vessel geometry. + Uses vessel geometry from fluids.TANK to convert liquid volume to height. + + Returns + ------- + float + Liquid level height from vessel bottom [m]. + Returns 0.0 for single-phase gas (quality > 1) or subcooled liquid (quality < 0). + + Raises + ------ + ValueError + If vessel_geometry was not provided during initialization + """ + if self.vessel_geometry is None: + raise ValueError( + "vessel_geometry must be provided to ThermodynamicSolver " + "for liquid level calculations" + ) + + # Check if in two-phase region + quality = self.fluid.Q() + if quality >= 0 and quality <= 1: + # Get saturated phase densities + rho_liq = self.fluid.saturated_liquid_keyed_output(CP.iDmass) + rho_vap = self.fluid.saturated_vapor_keyed_output(CP.iDmass) + + # Calculate liquid mass and volume + m_liq = self.fluid.rhomass() * self.vessel_geometry.V_total * (1 - quality) + V_liq = m_liq / rho_liq + + # Convert volume to height using vessel geometry + h_liq = self.vessel_geometry.h_from_V(V_liq) + return h_liq + else: + return 0.0 + + def PHres(self, T, P, H): + """ + Residual enthalpy function to be minimized during PH-problem. + + Used by numerical optimizer (scipy.optimize.minimize) to find temperature + that satisfies constant pressure-enthalpy constraints for multicomponent + fluids. The optimizer adjusts T until residual approaches zero. + + Parameters + ---------- + T : float or array-like + Temperature estimate [K]. Optimizer may pass as array. + P : float + Pressure [Pa] + H : float + Target enthalpy [J/kg] + + Returns + ------- + float + Squared normalized enthalpy residual (dimensionless). + Zero when correct temperature is found. + """ + # Extract scalar from array (scipy optimizers pass arrays, CoolProp needs scalars) + T_scalar = float(T.item()) if hasattr(T, "item") else float(T) + self.vent_fluid.update(CP.PT_INPUTS, P, T_scalar) + return ((H - self.vent_fluid.hmass()) / H) ** 2 + + def PHres_relief(self, T, P, H): + """ + Residual enthalpy function for PH-problem during relief valve calculations. + + Used by numerical optimizer (scipy.optimize.root_scalar) to find temperature + for relief valve discharge calculations. Similar to PHres() but uses the + main fluid state instead of vent_fluid state. + + Parameters + ---------- + T : float or array-like + Temperature estimate [K] + P : float + Pressure [Pa] + H : float + Target enthalpy [J/kg] + + Returns + ------- + float + Normalized enthalpy residual (dimensionless). + Zero when correct temperature is found. + """ + # Extract scalar from array (scipy optimizers pass arrays, CoolProp needs scalars) + T_scalar = float(T.item()) if hasattr(T, "item") else float(T) + self.fluid.update(CP.PT_INPUTS, P, T_scalar) + return (H - self.fluid.hmass()) / H + + def PHproblem(self, H, P, Tguess, relief=False): + """ + Solve constant pressure, constant enthalpy problem (isenthalpic expansion). + + Finds temperature at specified pressure and enthalpy. Typical use case is + modeling adiabatic throttling/expansion through valves without work extraction. + For multicomponent mixtures, uses numerical optimization to find temperature. + For single component fluids, uses direct CoolProp methods for speed. + + Parameters + ---------- + H : float + Enthalpy [J/kg] + P : float + Pressure [Pa] + Tguess : float + Initial guess for temperature [K] + relief : bool, optional + If True, use relief valve solver (root_scalar with Newton method). + If False, use general solver (minimize with Nelder-Mead). + Default is False. + + Returns + ------- + float + Temperature at (P, H) [K] + + Raises + ------ + ThermodynamicConvergenceError + If numerical optimization fails to converge for multicomponent fluid + """ + # Multicomponent case: requires numerical optimization + if self.is_multicomponent: + x0 = Tguess + + if not relief: + # Use Nelder-Mead minimization + res = minimize( + self.PHres, + x0, + args=(P, H), + method="Nelder-Mead", + options={"xatol": 0.1, "fatol": 0.001}, + ) + # Check convergence + sc.check_optimization_convergence( + res, solver_name="PHproblem", state_vars={"P": P, "H": H, "T_guess": x0} + ) + T1 = res.x[0] + else: + # Use Newton root finding for relief valve + res = root_scalar( + self.PHres_relief, + args=(P, H), + x0=x0, + method="newton", + ) + # Check convergence + sc.check_optimization_convergence( + res, solver_name="PHproblem_relief", state_vars={"P": P, "H": H} + ) + T1 = res.root + + # Single component case: direct CoolProp calculation + else: + T1 = PropsSI("T", "P", P, "H", H, self.species) + + return T1 + + def UDres(self, x, U, rho): + """ + Residual function for UD-problem (constant internal energy and density). + + Used by numerical optimizer to find pressure and temperature that satisfy + constant internal energy and density constraints. Minimizes sum of squared + normalized residuals for both U and rho. + + Parameters + ---------- + x : array-like of float + [Pressure [Pa], Temperature [K]] + U : float + Target internal energy [J/kg] + rho : float + Target density [kg/m³] + + Returns + ------- + float + Sum of squared normalized residuals (dimensionless) + """ + self.fluid.update(CP.PT_INPUTS, x[0], x[1]) + return ((U - self.fluid.umass()) / U) ** 2 + ((rho - self.fluid.rhomass()) / rho) ** 2 + + def UDproblem(self, U, rho, Pguess, Tguess): + """ + Solve constant internal energy, constant density problem. + + Finds pressure and temperature at specified internal energy and density. + Relevant for 1st law of thermodynamics with constant volume. For multicomponent + mixtures, uses numerical optimization to find P and T. For single component + fluids, uses direct CoolProp methods for speed. + + Parameters + ---------- + U : float + Internal energy [J/kg] + rho : float + Density [kg/m³] + Pguess : float + Initial guess for pressure [Pa] + Tguess : float + Initial guess for temperature [K] + + Returns + ------- + P1 : float + Pressure at (U, rho) [Pa] + T1 : float + Temperature at (U, rho) [K] + Ures : float + Internal energy residual [J/kg]. Zero for single component fluids. + + Raises + ------ + ThermodynamicConvergenceError + If numerical optimization fails to converge for multicomponent fluid + """ + # Multicomponent case: requires numerical optimization + if self.is_multicomponent: + x0 = [Pguess, Tguess] + res = minimize( + self.UDres, + x0, + args=(U, rho), + method="Nelder-Mead", + options={"xatol": 0.1, "fatol": 0.001}, + ) + # Check convergence + sc.check_optimization_convergence( + res, + solver_name="UDproblem", + state_vars={"U": U, "rho": rho, "P_guess": Pguess, "T_guess": Tguess}, + ) + P1 = res.x[0] + T1 = res.x[1] + Ures = U - self.fluid.umass() + + # Single component case: direct CoolProp calculation + else: + P1 = PropsSI("P", "D", rho, "U", U, self.species) + T1 = PropsSI("T", "D", rho, "U", U, self.species) + Ures = 0 + + return P1, T1, Ures diff --git a/src/hyddown/validator.py b/src/hyddown/validator.py index d4b5556..87191f1 100644 --- a/src/hyddown/validator.py +++ b/src/hyddown/validator.py @@ -26,8 +26,449 @@ - valve_validation(): Validates valve parameters based on valve type """ +import copy from cerberus import Validator from cerberus.errors import ValidationError +from hyddown.exceptions import ( + ConfigurationError, + VesselConfigurationError, + ValveConfigurationError, + HeatTransferConfigurationError, +) + + +# ============================================================================ +# REUSABLE SCHEMA COMPONENTS (Phase 1.5 Refactoring) +# ============================================================================ + +# Material property constraints (consistent across all uses) +MATERIAL_PROPERTIES = { + "heat_capacity": { + "required": False, + "type": "number", + "min": 1, + "max": 10000, # Sanity check [J/(kg·K)] + }, + "thermal_conductivity": { + "required": False, + "type": "number", + "min": 0.0001, + "max": 500, # Max for metals ~400 W/(m·K) + }, + "density": { + "required": False, + "type": "number", + "min": 1, + "max": 25000, # Max: osmium ~22,000 kg/m³ + }, + "thickness": { + "required": False, + "type": "number", + "min": 0.0001, # Wall thickness must be > 0 when specified + }, +} + + +def create_liner_properties(): + """Create liner property schema (same constraints as base materials).""" + return { + "liner_thickness": {"required": False, "type": "number", "min": 0.0}, + "liner_heat_capacity": {"required": False, "type": "number", "min": 1}, + "liner_thermal_conductivity": { + "required": False, + "type": "number", + "min": 0.0001, + }, + "liner_density": {"required": False, "type": "number", "min": 1}, + } + + +# Time series validation schema (for validation data) +TIME_SERIES_SCHEMA = { + "required": False, + "type": "dict", + "contains": ["time", "data"], + "schema": { + "data": { + "required": False, + "type": "list", + "schema": {"type": "number"}, + }, + "time": { + "required": False, + "type": "list", + "schema": {"type": "number"}, + }, + }, +} + + +def create_validation_data_schema(): + """Create schema for validation data section (eliminates massive duplication).""" + # Pressure validation data + pressure_schema = { + "type": "dict", + "required": False, + "contains": ["time", "pres"], + "schema": { + "pres": { + "required": False, + "type": "list", + "schema": {"type": "number"}, + }, + "time": { + "required": False, + "type": "list", + "schema": {"type": "number"}, + }, + }, + } + + # Temperature data schema (same for all temperature types) + temp_schema = { + "required": False, + "type": "dict", + "contains": ["time", "temp"], + "schema": { + "temp": { + "required": False, + "type": "list", + "schema": {"type": "number"}, + }, + "time": { + "required": False, + "type": "list", + "schema": {"type": "number"}, + }, + }, + } + + # All temperature types use the same schema + temp_types = [ + "gas_high", + "gas_low", + "gas_mean", + "wall_high", + "wall_mean", + "wall_low", + "wall_inner", + "wall_outer", + ] + + temperature_schema = { + "type": "dict", + "required": False, + "allowed": temp_types, + "schema": {temp_type: temp_schema.copy() for temp_type in temp_types}, + } + + return { + "pressure": pressure_schema, + "temperature": temperature_schema, + } + + +def create_vessel_schema(required_fields=None): + """ + Create vessel schema with specified required fields. + + Parameters + ---------- + required_fields : list, optional + Fields that should be required. If None, all are optional. + + Returns + ------- + dict + Vessel schema with specified requirements + """ + required_fields = required_fields or [] + + schema = { + "length": { + "type": "number", + "min": 0.001, # Length must be > 0 (m) + }, + "diameter": { + "type": "number", + "min": 0.001, # Diameter must be > 0 (m) + }, + "orientation": { + "type": "string", + "allowed": ["vertical", "horizontal"], + }, + "type": { + "type": "string", + "allowed": ["Flat-end", "ASME F&D", "DIN", "Hemispherical"], + }, + "liquid_level": { + "type": "number", + "min": 0, + }, + } + + # Add material properties (deep copy to avoid shared state) + schema.update(copy.deepcopy(MATERIAL_PROPERTIES)) + schema.update(create_liner_properties()) + + # Mark specified fields as required + for field in required_fields: + if field in schema: + schema[field]["required"] = True + + # Set all non-required fields + for field in schema: + if "required" not in schema[field]: + schema[field]["required"] = False + + return schema + + +def create_top_level_schema(): + """Create the common top-level section schema used across all validations.""" + return { + "initial": {"required": True}, + "calculation": {"required": True}, + "validation": {"required": False}, + "rupture": {"required": False}, + } + + +def create_valve_schema(valve_type): + """ + Create valve schema for specific valve type. + + Parameters + ---------- + valve_type : str + Valve type: 'relief', 'psv', 'orifice', 'controlvalve', 'mdot' + + Returns + ------- + dict + Valve schema with type-specific requirements + """ + # Base valve schema (all possible fields) + base_schema = { + "type": { + "required": True, + "type": "string", + "allowed": ["orifice", "psv", "controlvalve", "mdot", "relief"], + }, + "flow": { + "required": True, + "type": "string", + "allowed": ["discharge", "filling"], + }, + "diameter": {"type": "number", "min": 0}, + "discharge_coef": { + "type": "number", + "min": 0.001, + "max": 1.0, # Physical limit + }, + "set_pressure": {"type": "number", "min": 0}, + "end_pressure": {"type": "number", "min": 0}, + "blowdown": {"type": "number", "min": 0, "max": 1}, + "back_pressure": {"type": "number", "min": 0}, + "Cv": {"type": "number", "min": 0}, + "mdot": {"type": ["number", "list"]}, + "time": {"type": ["number", "list"]}, + "time_constant": {"type": "number", "min": 0}, + "characteristic": { + "type": "string", + "allowed": ["linear", "eq", "fast"], + }, + } + + # Type-specific requirements + if valve_type == "relief": + base_schema["set_pressure"]["required"] = True + base_schema["back_pressure"]["required"] = True + + elif valve_type == "psv": + base_schema["diameter"]["required"] = True + base_schema["discharge_coef"]["required"] = True + base_schema["set_pressure"]["required"] = True + base_schema["blowdown"]["required"] = True + base_schema["back_pressure"]["required"] = True + + elif valve_type == "orifice": + base_schema["diameter"]["required"] = True + base_schema["discharge_coef"]["required"] = True + base_schema["back_pressure"]["required"] = True + + elif valve_type == "controlvalve": + base_schema["Cv"]["required"] = True + base_schema["back_pressure"]["required"] = True + + elif valve_type == "mdot": + base_schema["mdot"]["required"] = True + base_schema["back_pressure"]["required"] = True + + return base_schema + + +def format_cerberus_errors(errors, section="input"): + """ + Convert Cerberus error dict to user-friendly ConfigurationError. + + Parameters + ---------- + errors : dict + Cerberus validation errors + section : str + Section name (vessel, valve, heat_transfer, etc.) + + Raises + ------ + ConfigurationError + With formatted, actionable error message + """ + error_messages = [] + + def flatten_errors(err_dict, prefix=""): + """Recursively flatten nested error dictionaries.""" + for field, error_list in err_dict.items(): + if isinstance(error_list, list): + for error in error_list: + if isinstance(error, dict): + # Nested errors - recurse + flatten_errors(error, f"{prefix}{field}.") + else: + # Leaf error - format it + error_messages.append(f" {prefix}{field}: {error}") + else: + error_messages.append(f" {prefix}{field}: {error_list}") + + flatten_errors(errors) + + if not error_messages: + error_messages = [f" Unknown validation error: {errors}"] + + message = f"Invalid {section} configuration:\n" + "\n".join(error_messages) + + # Raise specific exception based on section + if section == "vessel": + raise VesselConfigurationError(message) + elif section == "valve": + raise ValveConfigurationError(message) + elif section == "heat_transfer": + raise HeatTransferConfigurationError(message) + else: + raise ConfigurationError(message) + + +# ============================================================================ +# CROSS-FIELD VALIDATION (Physical Constraints) +# ============================================================================ + + +def validate_relief_valve_physics(input): + """ + Validate physical constraints for relief valves. + + Parameters + ---------- + input : dict + Structure holding input + + Raises + ------ + ValveConfigurationError + If physical constraints violated + """ + if input["valve"]["type"] not in ["relief", "psv"]: + return # Not a relief valve + + valve = input["valve"] + set_p = valve.get("set_pressure") + back_p = valve.get("back_pressure") + + # Check set pressure > back pressure + if set_p is not None and back_p is not None: + if set_p <= back_p: + raise ValveConfigurationError( + f"Relief valve set_pressure ({set_p} Pa) must be greater than " + f"back_pressure ({back_p} Pa).\n" + f"The valve cannot open if set pressure is at or below back pressure." + ) + + # Check blowdown is reasonable + blowdown = valve.get("blowdown") + if blowdown is not None and blowdown >= 1.0: + raise ValveConfigurationError( + f"Relief valve blowdown ({blowdown}) must be < 1.0.\n" + f"Typical values are 0.1-0.2 (10%-20%)." + ) + + +def validate_composite_vessel(input): + """ + Validate composite vessel liner constraints. + + Parameters + ---------- + input : dict + Structure holding input + + Raises + ------ + VesselConfigurationError + If liner thickness >= wall thickness + """ + vessel = input.get("vessel", {}) + thickness = vessel.get("thickness") + liner_thickness = vessel.get("liner_thickness") + + if thickness and liner_thickness: + if liner_thickness >= thickness: + raise VesselConfigurationError( + f"liner_thickness ({liner_thickness} m) must be less than " + f"thickness ({thickness} m).\n" + f"The liner cannot be thicker than the entire wall." + ) + + +def validate_time_step_sanity(input): + """ + Validate time step is reasonable for simulation duration. + + Parameters + ---------- + input : dict + Structure holding input + + Raises + ------ + ConfigurationError + If time step is too large relative to end time + """ + calc = input.get("calculation", {}) + time_step = calc.get("time_step") + end_time = calc.get("end_time") + + if time_step and end_time: + if time_step > end_time: + raise ConfigurationError( + f"time_step ({time_step} s) is greater than end_time ({end_time} s).\n" + f"The simulation would complete in a single step." + ) + + # Warn if very few steps + num_steps = end_time / time_step + if num_steps < 10: + import warnings + + warnings.warn( + f"Only {num_steps:.1f} time steps will be simulated. " + f"Consider reducing time_step for better resolution.", + UserWarning, + ) + + +# ============================================================================ +# VALIDATION FUNCTIONS +# ============================================================================ def validate_mandatory_ruleset(input): @@ -50,8 +491,16 @@ def validate_mandatory_ruleset(input): "type": "dict", "allow_unknown": False, "schema": { - "temperature": {"required": True, "type": "number"}, - "pressure": {"required": True, "type": "number"}, + "temperature": { + "required": True, + "type": "number", + "min": 0.001, # Kelvin scale, effectively > 0 + }, + "pressure": { + "required": True, + "type": "number", + "min": 0.001, # Pressure must be > 0 (Pa) + }, "fluid": {"required": True, "type": "string"}, }, }, @@ -71,49 +520,24 @@ def validate_mandatory_ruleset(input): "specified_U", ], }, - "time_step": {"required": True, "type": "number", "min": 0.000001}, - "end_time": {"required": True, "type": "number", "min": 0}, + "time_step": { + "required": True, + "type": "number", + "min": 0.000001, + "max": 3600, # Sanity check: max 1 hour per step + }, + "end_time": { + "required": True, + "type": "number", + "min": 0.001, # Must have positive duration + }, }, }, "vessel": { "required": True, "type": "dict", "allow_unknown": False, - "schema": { - "length": {"required": True, "type": "number"}, - "diameter": {"required": True, "type": "number"}, - "thickness": {"required": False, "type": "number", "min": 0.0}, - "heat_capacity": {"required": False, "type": "number", "min": 1}, - "thermal_conductivity": { - "required": False, - "type": "number", - "min": 0.0001, - }, - "density": {"required": False, "type": "number", "min": 1}, - "liner_thickness": {"required": False, "type": "number", "min": 0.0}, - "liner_heat_capacity": {"required": False, "type": "number", "min": 1}, - "liner_thermal_conductivity": { - "required": False, - "type": "number", - "min": 0.0001, - }, - "liner_density": {"required": False, "type": "number", "min": 1}, - "orientation": { - "required": False, - "type": "string", - "allowed": ["vertical", "horizontal"], - }, - "type": { - "required": False, - "type": "string", - "allowed": {"Flat-end", "ASME F&D", "DIN", "Hemispherical"}, - }, - "liquid_level": { - "required": False, - "type": "number", - "min": 0, - }, - }, + "schema": create_vessel_schema(required_fields=["length", "diameter"]), }, "rupture": { "required": False, @@ -143,8 +567,12 @@ def validate_mandatory_ruleset(input): "type": "string", "allowed": ["discharge", "filling"], }, - "diameter": {"type": "number", "min": 0}, - "discharge_coef": {"type": "number", "min": 0}, + "diameter": {"type": "number", "min": 0}, # Allow 0 for no-flow + "discharge_coef": { + "type": "number", + "min": 0.001, # > 0 + "max": 1.0, # Physical limit for ideal orifice + }, "set_pressure": {"type": "number", "min": 0}, "end_pressure": {"type": "number", "min": 0}, "blowdown": {"type": "number", "min": 0, "max": 1}, @@ -180,9 +608,24 @@ def validate_mandatory_ruleset(input): "allowed": ["specified_Q", "specified_h", "specified_U", "s-b"], }, "Q_fix": {"required": False, "type": "number"}, - "U_fix": {"required": False, "type": "number", "min": 0}, - "temp_ambient": {"required": False, "type": "number", "min": 0}, - "h_outer": {"required": False, "type": "number", "min": 0}, + "U_fix": { + "required": False, + "type": "number", + "min": 0.001, # > 0 + "max": 1000, # Sanity check [W/(m²·K)] + }, + "temp_ambient": { + "required": False, + "type": "number", + "min": 0.001, # Kelvin scale > 0 + "max": 2000, # Max for fire scenarios + }, + "h_outer": { + "required": False, + "type": "number", + "min": 0, + "max": 10000, # Sanity check [W/(m²·K)] + }, "h_inner": {"required": False, "type": ["number", "string"]}, "fire": { "required": False, @@ -202,682 +645,297 @@ def validate_mandatory_ruleset(input): "type": "dict", "allow_unknown": False, "allowed": ["pressure", "temperature"], - "schema": { - "pressure": { - "type": "dict", - "required": False, - "contains": ["time", "pres"], - "schema": { - "pres": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - "time": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - }, - }, - "temperature": { - "type": "dict", - "required": False, - "allowed": [ - "wall_high", - "wall_low", - "wall_mean", - "wall_inner", - "wall_outer", - "gas_high", - "gas_low", - "gas_mean", - ], - "schema": { - "gas_high": { - "required": False, - "type": "dict", - "contains": ["time", "temp"], - "schema": { - "temp": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - "time": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - }, - }, - "gas_low": { - "required": False, - "type": "dict", - "contains": ["time", "temp"], - "schema": { - "temp": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - "time": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - }, - }, - "gas_mean": { - "required": False, - "type": "dict", - "contains": ["time", "temp"], - "schema": { - "temp": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - "time": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - }, - }, - "wall_high": { - "required": False, - "type": "dict", - "contains": ["time", "temp"], - "schema": { - "temp": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - "time": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - }, - }, - "wall_mean": { - "required": False, - "type": "dict", - "contains": ["time", "temp"], - "schema": { - "temp": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - "time": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - }, - }, - "wall_low": { - "required": False, - "type": "dict", - "contains": ["time", "temp"], - "schema": { - "temp": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - "time": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - }, - }, - "wall_inner": { - "required": False, - "type": "dict", - "contains": ["time", "temp"], - "schema": { - "temp": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - "time": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - }, - }, - "wall_outer": { - "required": False, - "type": "dict", - "contains": ["time", "temp"], - "schema": { - "temp": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - "time": { - "required": False, - "type": "list", - "schema": {"type": "number"}, - }, - }, - }, - }, - }, - }, + "schema": create_validation_data_schema(), }, } v = Validator(schema_general) retval = v.validate(input) - if v.errors: - print(v.errors) + + if not retval: + # Raise specific exceptions instead of printing + if "vessel" in v.errors: + format_cerberus_errors({"vessel": v.errors["vessel"]}, "vessel") + elif "valve" in v.errors: + format_cerberus_errors({"valve": v.errors["valve"]}, "valve") + elif "heat_transfer" in v.errors: + format_cerberus_errors({"heat_transfer": v.errors["heat_transfer"]}, "heat_transfer") + elif "initial" in v.errors: + format_cerberus_errors({"initial": v.errors["initial"]}, "initial conditions") + elif "calculation" in v.errors: + format_cerberus_errors({"calculation": v.errors["calculation"]}, "calculation") + else: + # Generic configuration error for other cases + format_cerberus_errors(v.errors, "input") return retval def heat_transfer_validation(input): """ - Validate input['heat_transfer'] deeper than cerberus + Validate input['heat_transfer'] deeper than cerberus. Parameters ---------- input : dict Structure holding input - Return ---------- - : bool + bool True for success, False for failure """ retval = True if input["calculation"]["type"] == "energybalance": - if input["heat_transfer"]["type"] == "specified_h": - schema_heattransfer = { - "initial": {"required": True}, - "calculation": {"required": True}, - "validation": {"required": False}, - "rupture": {"required": False}, - "valve": {"required": True}, - "vessel": { - "required": True, - "type": "dict", - "allow_unknown": False, - "schema": { - "length": {"required": True, "type": "number"}, - "diameter": {"required": True, "type": "number"}, - "thickness": {"required": True, "type": "number", "min": 0.0}, - "heat_capacity": {"required": True, "type": "number", "min": 1}, - "thermal_conductivity": { - "required": False, - "type": "number", - "min": 0.0001, - }, - "density": {"required": True, "type": "number", "min": 1}, - "liner_thickness": { - "required": False, - "type": "number", - "min": 0.0, - }, - "liner_heat_capacity": { - "required": False, - "type": "number", - "min": 1, - }, - "liner_thermal_conductivity": { - "required": False, - "type": "number", - "min": 0.0001, - }, - "liner_density": { - "required": False, - "type": "number", - "min": 1, - }, - "orientation": { - "required": True, - "type": "string", - "allowed": ["vertical", "horizontal"], - }, - "type": { - "required": False, - "type": "string", - "allowed": ["Flat-end", "DIN", "ASME F&D", "Hemispherical"], - }, - "liquid_level": { - "required": False, - "type": "number", - "min": 0, - }, + ht_type = input["heat_transfer"]["type"] + + # Build schema using factory functions instead of massive duplication + # Use deepcopy to ensure no schema state is shared between validations + base_schema = copy.deepcopy(create_top_level_schema()) + base_schema["vessel"] = {"required": True} + base_schema["valve"] = {"required": True} + base_schema["heat_transfer"] = {"required": True} + + if ht_type == "specified_h": + # Vessel requirements for specified_h + required_vessel_fields = [ + "length", + "diameter", + "thickness", + "heat_capacity", + "density", + "orientation", + ] + + schema_heattransfer = copy.deepcopy(base_schema) + schema_heattransfer["vessel"] = { + "required": True, + "type": "dict", + "allow_unknown": False, + "schema": create_vessel_schema(required_fields=required_vessel_fields), + } + schema_heattransfer["heat_transfer"] = { + "required": True, + "type": "dict", + "allow_unknown": False, + "allowed": ["h_inner", "h_outer", "temp_ambient", "type", "D_throat"], + "schema": { + "type": {"type": "string", "allowed": ["specified_h"]}, + "temp_ambient": { + "required": True, + "type": "number", + "min": 0.001, + "max": 2000, }, - }, - "heat_transfer": { - "required": True, - "type": "dict", - "allow_unknown": False, - "allowed": [ - "h_inner", - "h_outer", - "temp_ambient", - "type", - "D_throat", - ], - "schema": { - "type": {"type": "string", "allowed": ["specified_h"]}, - "temp_ambient": {"required": True, "type": "number", "min": 0}, - "h_outer": {"required": True, "type": "number", "min": 0}, - "h_inner": {"required": True, "type": ["number", "string"]}, - "D_throat": {"required": False, "type": "number", "min": 0}, + "h_outer": { + "required": True, + "type": "number", + "min": 0, + "max": 10000, }, + "h_inner": {"required": True, "type": ["number", "string"]}, + "D_throat": {"required": False, "type": "number", "min": 0}, }, } - elif input["heat_transfer"]["type"] == "specified_Q": - schema_heattransfer = { - "initial": {"required": True}, - "calculation": {"required": True}, - "validation": {"required": False}, - "rupture": {"required": False}, - "valve": {"required": True}, - "vessel": { - "required": True, - "type": "dict", - "allow_unknown": False, - "schema": { - "length": {"required": True, "type": "number"}, - "diameter": {"required": True, "type": "number"}, - "thickness": {"required": False, "type": "number", "min": 0.0}, - "heat_capacity": { - "required": False, - "type": "number", - "min": 1, - }, - "density": {"required": False, "type": "number", "min": 1}, - "orientation": { - "required": False, - "type": "string", - "allowed": ["vertical", "horizontal"], - }, - "type": { - "required": False, - "type": "string", - "allowed": ["Flat-end", "DIN", "ASME F&D", "Hemispherical"], - }, - "liquid_level": { - "required": False, - "type": "number", - "min": 0, - }, - }, - }, - "heat_transfer": { - "required": True, - "type": "dict", - "allow_unknown": False, - "allowed": ["Q_fix", "type"], - "schema": { - "type": {"type": "string", "allowed": ["specified_Q"]}, - "Q_fix": {"required": False, "type": "number"}, - }, + elif ht_type == "specified_Q": + # Vessel requirements for specified_Q (minimal) + required_vessel_fields = ["length", "diameter"] + + schema_heattransfer = copy.deepcopy(base_schema) + schema_heattransfer["vessel"] = { + "required": True, + "type": "dict", + "allow_unknown": False, + "schema": create_vessel_schema(required_fields=required_vessel_fields), + } + schema_heattransfer["heat_transfer"] = { + "required": True, + "type": "dict", + "allow_unknown": False, + "allowed": ["Q_fix", "type"], + "schema": { + "type": {"type": "string", "allowed": ["specified_Q"]}, + "Q_fix": {"required": False, "type": "number"}, }, } - elif input["heat_transfer"]["type"] == "specified_U": - schema_heattransfer = { - "initial": {"required": True}, - "calculation": {"required": True}, - "validation": {"required": False}, - "valve": {"required": True}, - "rupture": {"required": False}, - "vessel": { - "required": True, - "type": "dict", - "allow_unknown": False, - "schema": { - "length": {"required": True, "type": "number"}, - "diameter": {"required": True, "type": "number"}, - "thickness": {"required": False, "type": "number", "min": 0.0}, - "heat_capacity": { - "required": False, - "type": "number", - "min": 1, - }, - "density": {"required": False, "type": "number", "min": 1}, - "orientation": { - "required": False, - "type": "string", - "allowed": ["vertical", "horizontal"], - }, - "type": { - "required": False, - "type": "string", - "allowed": ["Flat-end", "DIN", "ASME F&D", "Hemispherical"], - }, - "liquid_level": { - "required": False, - "type": "number", - "min": 0, - }, + elif ht_type == "specified_U": + # Vessel requirements for specified_U (minimal) + required_vessel_fields = ["length", "diameter"] + + schema_heattransfer = copy.deepcopy(base_schema) + schema_heattransfer["vessel"] = { + "required": True, + "type": "dict", + "allow_unknown": False, + "schema": create_vessel_schema(required_fields=required_vessel_fields), + } + schema_heattransfer["heat_transfer"] = { + "required": True, + "type": "dict", + "allow_unknown": False, + "allowed": ["U_fix", "type", "temp_ambient"], + "schema": { + "type": {"type": "string", "allowed": ["specified_U"]}, + "U_fix": { + "required": False, + "type": "number", + "min": 0.001, + "max": 1000, }, - }, - "heat_transfer": { - "required": True, - "type": "dict", - "allow_unknown": False, - "allowed": ["U_fix", "type", "temp_ambient"], - "schema": { - "type": {"type": "string", "allowed": ["specified_U"]}, - "U_fix": {"required": False, "type": "number", "min": 0.0}, - "temp_ambient": {"required": True, "type": "number", "min": 0}, + "temp_ambient": { + "required": True, + "type": "number", + "min": 0.001, + "max": 2000, }, }, } - elif input["heat_transfer"]["type"] == "s-b": - schema_heattransfer = { - "initial": {"required": True}, - "calculation": {"required": True}, - "validation": {"required": False}, - "valve": {"required": True}, - "rupture": {"required": False}, - "vessel": { - "required": True, - "type": "dict", - "allow_unknown": False, - "schema": { - "length": {"required": True, "type": "number"}, - "diameter": {"required": True, "type": "number"}, - "thickness": {"required": True, "type": "number", "min": 0.0}, - "heat_capacity": {"required": True, "type": "number", "min": 1}, - "density": {"required": True, "type": "number", "min": 1}, - "orientation": { - "required": True, - "type": "string", - "allowed": ["vertical", "horizontal"], - }, - "type": { - "required": False, - "type": "string", - "allowed": ["Flat-end", "DIN", "ASME F&D", "Hemispherical"], - }, - "liquid_level": { - "required": False, - "type": "number", - "min": 0, - }, + elif ht_type == "s-b": + # Vessel requirements for s-b (Stefan-Boltzmann fire) + required_vessel_fields = [ + "length", + "diameter", + "thickness", + "heat_capacity", + "density", + "orientation", + ] + + schema_heattransfer = copy.deepcopy(base_schema) + schema_heattransfer["vessel"] = { + "required": True, + "type": "dict", + "allow_unknown": False, + "schema": create_vessel_schema(required_fields=required_vessel_fields), + } + schema_heattransfer["heat_transfer"] = { + "required": True, + "type": "dict", + "allow_unknown": False, + "allowed": ["fire", "type"], + "schema": { + "type": { + "required": True, + "type": "string", + "allowed": ["s-b"], }, - }, - "heat_transfer": { - "required": True, - "type": "dict", - "allow_unknown": False, - "allowed": ["fire", "type"], - "schema": { - "type": { - "required": True, - "type": "string", - "allowed": ["s-b"], - }, - "fire": { - "required": True, - "type": "string", - "allowed": [ - "api_pool", - "api_jet", - "scandpower_pool", - "scandpower_jet", - ], - }, + "fire": { + "required": True, + "type": "string", + "allowed": [ + "api_pool", + "api_jet", + "scandpower_pool", + "scandpower_jet", + ], }, }, } + # Validate with the schema v = Validator(schema_heattransfer) retval = v.validate(input) - if v.errors: - print(v.errors) + + if not retval: + # Raise specific exceptions for heat transfer errors + if "heat_transfer" in v.errors: + format_cerberus_errors( + {"heat_transfer": v.errors["heat_transfer"]}, "heat_transfer" + ) + elif "vessel" in v.errors: + format_cerberus_errors({"vessel": v.errors["vessel"]}, "vessel") + else: + format_cerberus_errors(v.errors, f"heat_transfer ({ht_type})") return retval def valve_validation(input): """ - Validate input['valve'] deeper than cerberus + Validate input['valve'] deeper than cerberus. Parameters ---------- input : dict Structure holding input - Return ---------- - : bool + bool True for success, False for failure """ - schema_relief = { - "initial": {"required": True}, - "calculation": {"required": True}, - "validation": {"required": False}, - "vessel": {"required": True}, - "rupture": {"required": False}, - "heat_transfer": {"required": True}, - "valve": { - "required": True, - "type": "dict", - "allow_unknown": False, - "schema": { - "type": { - "required": True, - "type": "string", - "allowed": ["orifice", "psv", "controlvalve", "mdot", "relief"], - }, - "flow": { - "required": True, - "type": "string", - "allowed": ["discharge", "filling"], - }, - "diameter": {"required": False, "type": "number", "min": 0}, - "discharge_coef": {"required": False, "type": "number", "min": 0}, - "set_pressure": {"required": True, "type": "number", "min": 0}, - "end_pressure": {"type": "number", "min": 0}, - "blowdown": {"required": False, "type": "number", "min": 0, "max": 1}, - "back_pressure": {"required": True, "type": "number", "min": 0}, - "Cv": {"type": "number", "min": 0}, - "mdot": {"type": ["number", "list"]}, - "time": {"type": "list"}, - }, - }, + valve_type = input["valve"]["type"] + + # Build schema using factory function + base_schema = create_top_level_schema() + base_schema["vessel"] = {"required": True} + + # Heat transfer requirements vary by valve type + if valve_type == "relief": + base_schema["heat_transfer"] = {"required": True} + else: + base_schema["heat_transfer"] = {"required": False} + + # Create valve schema with type-specific requirements + base_schema["valve"] = { + "required": True, + "type": "dict", + "allow_unknown": False, + "schema": create_valve_schema(valve_type), } - schema_psv = { - "initial": {"required": True}, - "calculation": {"required": True}, - "validation": {"required": False}, - "vessel": {"required": True}, - "rupture": {"required": False}, - "heat_transfer": {"required": False}, - "valve": { - "required": True, - "type": "dict", - "allow_unknown": False, - "schema": { - "type": { - "required": True, - "type": "string", - "allowed": ["orifice", "psv", "controlvalve", "mdot"], - }, - "flow": { - "required": True, - "type": "string", - "allowed": ["discharge", "filling"], - }, - "diameter": {"required": True, "type": "number", "min": 0}, - "discharge_coef": {"required": True, "type": "number", "min": 0}, - "set_pressure": {"required": True, "type": "number", "min": 0}, - "end_pressure": {"type": "number", "min": 0}, - "blowdown": {"required": True, "type": "number", "min": 0, "max": 1}, - "back_pressure": {"required": True, "type": "number", "min": 0}, - "Cv": {"type": "number", "min": 0}, - "mdot": {"type": ["number", "list"]}, - "time": {"type": "list"}, - }, - }, - } - - schema_orifice = { - "initial": {"required": True}, - "calculation": {"required": True}, - "validation": {"required": False}, - "vessel": {"required": True}, - "rupture": {"required": False}, - "heat_transfer": {"required": False}, - "valve": { - "required": True, - "type": "dict", - "allow_unknown": False, - "schema": { - "type": { - "required": True, - "type": "string", - "allowed": ["orifice", "psv", "controlvalve", "mdot"], - }, - "flow": { - "required": True, - "type": "string", - "allowed": ["discharge", "filling"], - }, - "diameter": {"required": True, "type": "number", "min": 0}, - "discharge_coef": {"required": True, "type": "number", "min": 0}, - "set_pressure": {"type": "number", "min": 0}, - "end_pressure": {"type": "number", "min": 0}, - "blowdown": {"type": "number", "min": 0, "max": 1}, - "back_pressure": {"required": True, "type": "number", "min": 0}, - "Cv": {"type": "number", "min": 0}, - "mdot": {"type": ["number", "list"]}, - "time": {"type": "list"}, - }, - }, - } - - schema_control_valve = { - "initial": {"required": True}, - "rupture": {"required": False}, - "calculation": {"required": True}, - "validation": {"required": False}, - "vessel": {"required": True}, - "heat_transfer": {"required": False}, - "valve": { - "required": True, - "type": "dict", - "allow_unknown": False, - "schema": { - "type": { - "required": True, - "type": "string", - "allowed": ["orifice", "psv", "controlvalve", "mdot"], - }, - "flow": { - "required": True, - "type": "string", - "allowed": ["discharge", "filling"], - }, - "back_pressure": {"required": True, "type": "number", "min": 0}, - "Cv": {"required": True, "type": "number", "min": 0}, - "time_constant": {"type": "number", "min": 0}, - "characteristic": { - "type": "string", - "allowed": ["linear", "eq", "fast"], - }, - }, - }, - } - - schema_mdot = { - "initial": {"required": True}, - "calculation": {"required": True}, - "validation": {"required": False}, - "vessel": {"required": True}, - "rupture": {"required": False}, - "heat_transfer": {"required": False}, - "valve": { - "required": True, - "type": "dict", - "allow_unknown": False, - "schema": { - "type": { - "required": True, - "type": "string", - "allowed": ["orifice", "psv", "controlvalve", "mdot"], - }, - "flow": { - "required": True, - "type": "string", - "allowed": ["discharge", "filling"], - }, - "mdot": {"required": True, "type": ["number", "list"]}, - "time": {"type": ["number", "list"]}, - "back_pressure": {"required": True, "type": "number", "min": 0}, - }, - }, - } - - if input["valve"]["type"] == "relief": - v = Validator(schema_relief) - retval = v.validate(input) - if v.errors: - print(v.errors) + # Validate + v = Validator(base_schema) + retval = v.validate(input) - if input["valve"]["type"] == "psv": - v = Validator(schema_psv) - retval = v.validate(input) - if v.errors: - print(v.errors) - elif input["valve"]["type"] == "orifice": - v = Validator(schema_orifice) - retval = v.validate(input) - if v.errors: - print(v.errors) - elif input["valve"]["type"] == "controlvalve": - v = Validator(schema_control_valve) - retval = v.validate(input) - if v.errors: - print(v.errors) - elif input["valve"]["type"] == "mdot": - v = Validator(schema_mdot) - retval = v.validate(input) - if v.errors: - print(v.errors) + if not retval: + # Raise specific exceptions for valve errors + if "valve" in v.errors: + format_cerberus_errors({"valve": v.errors["valve"]}, "valve") + else: + format_cerberus_errors(v.errors, f"valve ({valve_type})") return retval def validation(input): """ - Aggregate validation using cerberus + Aggregate validation using cerberus and cross-field checks. Parameters ---------- input : dict Structure holding input - Return ---------- - : bool + bool True for success, False for failure """ - return ( + # Schema-based validation + schema_valid = ( validate_mandatory_ruleset(input) and valve_validation(input) and heat_transfer_validation(input) ) + + if not schema_valid: + return False + + # Cross-field validation (physical constraints) + try: + validate_relief_valve_physics(input) + validate_composite_vessel(input) + validate_time_step_sanity(input) + except ( + ConfigurationError, + VesselConfigurationError, + ValveConfigurationError, + ): + # Re-raise configuration errors + raise + + return True