///|
fn validate_constant(
data : ThermodynamicData,
source : String,
log_k : Double,
) -> Result[Unit, SolveError] {
match thermo_record(data, source) {
None => Err(DatabaseInconsistent("missing thermodynamic record: " + source))
Some(record) =>
if !finite(log_k) || (record.log_k - log_k).abs() > 1.0e-12 {
Err(
DatabaseInconsistent("logK differs from its source record: " + source),
)
} else {
Ok(())
}
}
}
///|
/// Validate the supported primary basis, reaction charge/element/alkalinity
/// balances, finite numbers, explicit thermodynamic provenance and controls.
fn validate_problem(problem : Problem) -> Result[Unit, SolveError] {
if !finite(problem.temperature_k) ||
(problem.temperature_k - 298.15).abs() > 1.0e-8 {
return Err(
InvalidInput(
"temperature_k must be 298.15; no logK temperature model is implemented",
),
)
}
if !finite(problem.volume_l) || problem.volume_l <= 0.0 {
return Err(InvalidInput("volume_l must be finite and positive"))
}
let data = problem.thermodynamics
if problem.database_version.length() == 0 ||
problem.database_version != data.version {
return Err(
DatabaseInconsistent("database_version must match thermodynamics.version"),
)
}
if data.adaptation_notice.length() == 0 ||
data.reference_url.length() == 0 ||
data.reference_sha256.length() != 64 ||
data.license.length() == 0 ||
data.license_url.length() == 0 ||
data.source_standard_state.length() == 0 ||
data.solver_standard_state.length() == 0 ||
data.limitations.length() == 0 {
return Err(
DatabaseInconsistent(
"thermodynamic source, hash, license, standard states and limitations are required",
),
)
}
if !finite(data.reference_temperature_k) ||
!finite(data.min_temperature_k) ||
!finite(data.max_temperature_k) ||
!finite(data.pressure_atm) ||
(data.reference_temperature_k - 298.15).abs() > 1.0e-8 ||
data.min_temperature_k > problem.temperature_k ||
data.max_temperature_k < problem.temperature_k ||
data.min_temperature_k > data.max_temperature_k ||
(data.pressure_atm - 1.0).abs() > 1.0e-8 {
return Err(DatabaseInconsistent("dataset must cover 298.15 K at 1 atm"))
}
let record_ids : Array[String] = []
for record in data.records {
if record.id.length() == 0 ||
record.reaction.length() == 0 ||
record.source_reaction.length() == 0 ||
record.source_lines.length() == 0 ||
record.derivation.length() == 0 ||
!finite(record.log_k) ||
record_ids.contains(record.id) {
return Err(
DatabaseInconsistent(
"invalid or duplicate thermodynamic record: " + record.id,
),
)
}
record_ids.push(record.id)
}
match activity_model_name(problem.activity_model) {
None =>
return Err(
ActivityModelOutOfRange("activity_model must be Davies or DebyeHuckel"),
)
Some(_) => ()
}
if problem.fixed_p_h is Some(p_h) {
if !finite(p_h) || p_h < 0.0 || p_h > 14.0 {
return Err(InvalidInput("fixed_p_h must be finite and within 0..14"))
}
if problem.alkalinity_mol_l is Some(_) {
return Err(
InconsistentConstraints(
"fixed_p_h and alkalinity_mol_l cannot both be hard constraints",
),
)
}
}
if problem.alkalinity_mol_l is Some(value) {
if !finite(value) {
return Err(InvalidInput("alkalinity_mol_l must be finite"))
}
}
let count = problem.elements.length() + 1
if problem.hydrogen_component < 0 || problem.hydrogen_component >= count {
return Err(
DimensionMismatch("hydrogen_component must identify a primary component"),
)
}
let names : Array[String] = []
for element in problem.elements {
if element.name.length() == 0 || names.contains(element.name) {
return Err(DatabaseInconsistent("empty or duplicate element name"))
}
if !finite(element.total_mol_l) || element.total_mol_l <= 0.0 {
return Err(
InvalidInput(
"analytical totals must be finite and strictly positive; remove absent components",
),
)
}
names.push(element.name)
}
let primary = Array::makei(count, fn(_) { -1 })
let species_names : Array[String] = []
for i = 0; i < problem.species.length(); i = i + 1 {
let species = problem.species[i]
if species.name.length() == 0 || species_names.contains(species.name) {
return Err(DatabaseInconsistent("empty or duplicate species name"))
}
species_names.push(species.name)
if species.element_stoich.length() != problem.elements.length() {
return Err(
DimensionMismatch(
"species element_stoich length does not match elements",
),
)
}
if !finite(species.alkalinity) {
return Err(DatabaseInconsistent("species alkalinity must be finite"))
}
for value in species.element_stoich {
if !finite(value) || value < 0.0 {
return Err(
DatabaseInconsistent(
"element stoichiometry must be finite and non-negative",
),
)
}
}
match validate_constant(data, species.source, species.log_k) {
Err(error) => return Err(error)
Ok(_) => ()
}
let seen_terms : Array[Int32] = []
for term in species.component_terms {
if term.component < 0 ||
term.component >= count ||
!finite(term.coefficient) ||
term.coefficient == 0.0 ||
seen_terms.contains(term.component) {
return Err(
DatabaseInconsistent(
"invalid or duplicate species component term: " + species.name,
),
)
}
seen_terms.push(term.component)
}
if species.base_component is Some(component) {
if component < 0 || component >= count {
return Err(DimensionMismatch("invalid base_component"))
}
if primary[component] != -1 ||
species.component_terms.length() != 0 ||
species.log_k != 0.0 {
return Err(
DatabaseInconsistent(
"primary components must be unique identity definitions",
),
)
}
primary[component] = i
}
}
for index in primary {
if index == -1 {
return Err(
DimensionMismatch(
"one primary per analytical element plus H+ is required",
),
)
}
}
let hydrogen = problem.species[primary[problem.hydrogen_component]]
if hydrogen.name != "H+" ||
hydrogen.charge != 1 ||
hydrogen.alkalinity != -1.0 ||
hydrogen.element_stoich.any(fn(v) { v != 0.0 }) {
return Err(
DatabaseInconsistent(
"hydrogen_component must refer to H+ with zero analytical stoichiometry",
),
)
}
let element_basis = Array::makei(problem.elements.length(), fn(_) { false })
for component = 0; component < count; component = component + 1 {
if component != problem.hydrogen_component {
let species = problem.species[primary[component]]
let mut nonzero = 0
for e = 0; e < element_basis.length(); e = e + 1 {
if species.element_stoich[e] != 0.0 {
if species.element_stoich[e] != 1.0 || element_basis[e] {
return Err(
DatabaseInconsistent(
"analytical primary basis must be one-to-one with elements",
),
)
}
element_basis[e] = true
nonzero = nonzero + 1
}
}
if nonzero != 1 {
return Err(
DatabaseInconsistent(
"each analytical primary must contain one analytical element",
),
)
}
}
}
for species in problem.species {
if species.base_component is None {
let totals = Array::makei(problem.elements.length(), fn(_) { 0.0 })
let mut charge = 0.0
let mut alkalinity = 0.0
for term in species.component_terms {
let base = problem.species[primary[term.component]]
charge = charge + term.coefficient * base.charge.to_double()
alkalinity = alkalinity + term.coefficient * base.alkalinity
for e = 0; e < totals.length(); e = e + 1 {
totals[e] = totals[e] + term.coefficient * base.element_stoich[e]
}
}
for e = 0; e < totals.length(); e = e + 1 {
if (totals[e] - species.element_stoich[e]).abs() > 1.0e-12 {
return Err(
DatabaseInconsistent(
"reaction violates analytical element balance: " + species.name,
),
)
}
}
if (charge - species.charge.to_double()).abs() > 1.0e-12 ||
(alkalinity - species.alkalinity).abs() > 1.0e-12 {
return Err(
DatabaseInconsistent(
"reaction violates charge or alkalinity balance: " + species.name,
),
)
}
}
}
let mineral_names : Array[String] = []
for mineral in problem.minerals {
if mineral.name.length() == 0 ||
mineral_names.contains(mineral.name) ||
mineral.component_terms.length() == 0 {
return Err(DatabaseInconsistent("empty or duplicate mineral definition"))
}
mineral_names.push(mineral.name)
match validate_constant(data, mineral.source, mineral.log_ksp) {
Err(error) => return Err(error)
Ok(_) => ()
}
let mut charge = 0.0
let seen : Array[Int32] = []
for term in mineral.component_terms {
if term.component < 0 ||
term.component >= count ||
!finite(term.coefficient) ||
term.coefficient == 0.0 ||
seen.contains(term.component) {
return Err(DatabaseInconsistent("invalid mineral component term"))
}
seen.push(term.component)
charge = charge +
term.coefficient *
problem.species[primary[term.component]].charge.to_double()
}
if charge.abs() > 1.0e-12 {
return Err(
DatabaseInconsistent(
"mineral dissolution must conserve charge: " + mineral.name,
),
)
}
}
Ok(())
}