///|
/// Central differences at fixed database, activity model and other inputs.
/// `kind` is fixed_pH, alkalinity, element_total or total_scale. Units are
/// pH, mol/L, mol/L, fraction; element_index is used only for element_total.
pub(all) struct Perturbation {
kind : String
element_index : Int32
step : Double
} derive(Debug, ToJson, FromJson)
///|
pub(all) struct SpeciesSensitivity {
name : String
minus_mol_l : Double
base_mol_l : Double
plus_mol_l : Double
d_concentration_d_input : Double
d_log10_concentration_d_input : Double
} derive(Debug, ToJson, FromJson)
///|
pub(all) struct MineralSensitivity {
mineral : String
minus_si : Double
base_si : Double
plus_si : Double
d_si_d_input : Double
phase_boundary_ambiguous : Bool
} derive(Debug, ToJson, FromJson)
///|
pub(all) struct SensitivityReport {
difference_method : String
perturbation : Perturbation
input_value : Double
input_unit : String
d_p_h_d_input : Double
d_ionic_strength_d_input : Double
species : Array[SpeciesSensitivity]
minerals : Array[MineralSensitivity]
base_result : SpeciationResult
minus_result : SpeciationResult
plus_result : SpeciationResult
diagnostics : Array[String]
} derive(Debug, ToJson, FromJson)
///|
pub(all) struct SensitivityRequest {
problem : Problem
perturbation : Perturbation
} derive(Debug, ToJson, FromJson)
///|
fn perturb_problem(
problem : Problem,
perturbation : Perturbation,
delta : Double,
) -> Problem {
let elements = problem.elements.copy()
if perturbation.kind == "total_scale" {
for i = 0; i < elements.length(); i = i + 1 {
let old = elements[i]
elements[i] = {
name: old.name,
total_mol_l: old.total_mol_l * (1.0 + delta),
}
}
}
if perturbation.kind == "element_total" {
let old = elements[perturbation.element_index]
elements[perturbation.element_index] = {
name: old.name,
total_mol_l: old.total_mol_l + delta,
}
}
let fixed_p_h = if perturbation.kind == "fixed_pH" {
match problem.fixed_p_h {
Some(value) => Some(value + delta)
None => None
}
} else {
problem.fixed_p_h
}
let alkalinity_mol_l = if perturbation.kind == "alkalinity" {
match problem.alkalinity_mol_l {
Some(value) => Some(value + delta)
None => None
}
} else {
problem.alkalinity_mol_l
}
{ ..problem, elements, fixed_p_h, alkalinity_mol_l, }
}
///|
/// Solve the centre and both perturbed inputs. Any failed solve is returned
/// as an error; a local derivative does not prove behaviour beyond the step.
pub fn sensitivity(
problem : Problem,
perturbation : Perturbation,
options~ : SolveOptions = default_options(),
) -> Result[SensitivityReport, SolveError] {
if !finite(perturbation.step) || perturbation.step <= 0.0 {
return Err(InvalidInput("sensitivity step must be finite and positive"))
}
let value = match perturbation.kind {
"fixed_pH" =>
match problem.fixed_p_h {
Some(value) => value
None =>
return Err(
InvalidInput("fixed_pH sensitivity requires a fixed_p_h boundary"),
)
}
"alkalinity" =>
match problem.alkalinity_mol_l {
Some(value) => value
None =>
return Err(
InvalidInput(
"alkalinity sensitivity requires an alkalinity boundary",
),
)
}
"element_total" => {
if perturbation.element_index < 0 ||
perturbation.element_index >= problem.elements.length() {
return Err(InvalidInput("sensitivity element_index is out of range"))
}
let value = problem.elements[perturbation.element_index].total_mol_l
if value <= perturbation.step {
return Err(
InvalidInput(
"central element_total step must be smaller than the positive total",
),
)
}
value
}
"total_scale" => {
if perturbation.step >= 1.0 || problem.elements.length() == 0 {
return Err(
InvalidInput(
"total_scale requires analytical elements and a fractional step smaller than one",
),
)
}
1.0
}
_ =>
return Err(
InvalidInput(
"sensitivity kind must be fixed_pH, alkalinity, element_total or total_scale",
),
)
}
if value + perturbation.step == value || value - perturbation.step == value {
return Err(
InvalidInput(
"sensitivity step is below floating-point resolution at the input value",
),
)
}
let base = match solve(problem, options~) {
Ok(result) => result
Err(error) => return Err(error)
}
let minus = match
solve(perturb_problem(problem, perturbation, -perturbation.step), options~) {
Ok(result) => result
Err(error) => return Err(error)
}
let plus = match
solve(perturb_problem(problem, perturbation, perturbation.step), options~) {
Ok(result) => result
Err(error) => return Err(error)
}
let divisor = 2.0 * perturbation.step
let species = Array::makei(base.species.length(), fn(i) {
let left = minus.species[i].concentration_mol_l
let right = plus.species[i].concentration_mol_l
SpeciesSensitivity::{
name: base.species[i].name,
minus_mol_l: left,
base_mol_l: base.species[i].concentration_mol_l,
plus_mol_l: right,
d_concentration_d_input: (right - left) / divisor,
d_log10_concentration_d_input: (safe_log10(right) - safe_log10(left)) /
divisor,
}
})
let minerals = Array::makei(base.saturation.length(), fn(i) {
let left = minus.saturation[i].saturation_index
let centre = base.saturation[i].saturation_index
let right = plus.saturation[i].saturation_index
MineralSensitivity::{
mineral: base.saturation[i].mineral,
minus_si: left,
base_si: centre,
plus_si: right,
d_si_d_input: (right - left) / divisor,
phase_boundary_ambiguous: left.min(centre).min(right) <= 0.0 &&
left.max(centre).max(right) >= 0.0,
}
})
let d_p_h = (plus.p_h - minus.p_h) / divisor
let d_ionic = (plus.ionic_strength - minus.ionic_strength) / divisor
if !finite(d_p_h) ||
!finite(d_ionic) ||
species.any(fn(row) {
!finite(row.d_concentration_d_input) ||
!finite(row.d_log10_concentration_d_input)
}) ||
minerals.any(fn(row) { !finite(row.d_si_d_input) }) {
return Err(
InvalidInput("sensitivity derivative overflow; choose a resolvable step"),
)
}
Ok({
difference_method: "Central difference; three complete equilibrium solves; fixed phase set",
perturbation,
input_value: value,
input_unit: if perturbation.kind == "fixed_pH" {
"pH"
} else if perturbation.kind == "total_scale" {
"fraction"
} else {
"mol/L"
},
d_p_h_d_input: d_p_h,
d_ionic_strength_d_input: d_ionic,
species,
minerals,
base_result: base,
minus_result: minus,
plus_result: plus,
diagnostics: [
"Local numerical derivative only; check step refinement. SI crossing zero is PhaseBoundaryAmbiguous, not precipitation or a phase switch.",
],
})
}