use std::vec;
use crate::Thermodynamics::ChemEquilibrium::ClassicalThermodynamics::ThermodynamicCalculations;
use crate::Thermodynamics::User_PhaseOrSolution::SubstanceSystemFactory;
use crate::Thermodynamics::User_PhaseOrSolution::{CustomSubstance, SubstancesContainer};
use RustedSciThe::symbolic::symbolic_engine::Expr;
use approx::assert_relative_eq;
use std::collections::HashMap;
pub fn SubsData_examples(thermotask: usize) {
match thermotask {
3 => {
use crate::Thermodynamics::ChemEquilibrium::ClassicalThermodynamics::ThermodynamicCalculations;
use crate::Thermodynamics::User_PhaseOrSolution::{
SubstanceSystemFactory, SubstancesContainer,
};
use crate::Thermodynamics::User_PhaseOrSolution::ThermodynamicsCalculatorTrait;
use nalgebra::{DMatrix, DVector};
let subs = vec!["H2".to_string(), "O2".to_string()];
let T = 400.0;
let P = 101325.0;
let mut n = HashMap::new();
let map_of_moles_num =
HashMap::from([("H2".to_string(), 0.5), ("O2".to_string(), 0.5)]);
n.insert(None, (Some(1.0), (Some(map_of_moles_num))));
let mut n_sym: HashMap<Option<String>, (Option<Expr>, Option<Vec<Expr>>)> =
HashMap::new();
n_sym.insert(
None,
(
Some(Expr::Var("Np".to_string())),
Some(vec![
Expr::Var("N0".to_string()),
Expr::Var("N1".to_string()),
]),
),
);
let search_in_NIST = false;
let explicit_search_insructions = None;
let library_priorities = vec!["NASA_gas".to_string()];
let permitted_libraries = vec!["NASA_gas".to_string()];
let container = SubstancesContainer::SinglePhase(subs.clone());
let mut customsubs = SubstanceSystemFactory::create_system(
container,
None,
library_priorities,
permitted_libraries,
explicit_search_insructions,
search_in_NIST,
)
.unwrap();
let mut nv = HashMap::new();
nv.insert(None, (Some(1.0), (Some(vec![0.5, 0.5]))));
let dG = customsubs.calcutate_Gibbs_free_energy(T, P, nv.clone());
let dG_sym = customsubs.calculate_Gibbs_sym(T);
assert!(dG.is_ok());
assert!(dG_sym.is_ok());
let td = customsubs.create_thermodynamics(Some(T), P, Some(n.clone()), None);
assert!(td.is_ok());
let mut td = td.unwrap();
td.set_P_to_sym();
let bind = td.subdata.get_dG_sym();
let g_sym = bind.get(&None).unwrap();
let g_sym = g_sym.get(&"H2".to_string()).unwrap().clone();
let g_sym_fun = g_sym.lambdify_wrapped();
let g_sym_value = g_sym_fun(vec![400.0, 0.5, 0.5]);
let g_sym_value_from_sym = g_sym_fun(vec![400.0, 0.5, 0.5]);
assert_relative_eq!(g_sym_value, g_sym_value_from_sym, epsilon = 1e-6);
td.calculate_elem_composition_and_molar_mass(None);
td.map_of_concentration =
HashMap::from([("H2".to_string(), 0.5), ("O2".to_string(), 0.5)]);
let elem_composition_matrix = td.elem_composition_matrix.clone().unwrap();
let vec_of_elements = &td.unique_elements;
assert_eq!(vec_of_elements.len(), 2);
assert_eq!(
vec_of_elements.clone(),
vec!["H".to_string(), "O".to_string()]
);
let matrix_data = vec![
2.0, 0.0, 0.0, 2.0, ];
let matrix_expected = DMatrix::from_row_slice(2, 2, &matrix_data);
assert_eq!(elem_composition_matrix, matrix_expected);
td.initial_composition().unwrap();
let initial_composition = &td.initial_vector_of_elements;
println!("initial_composition: {:?}", initial_composition);
assert_eq!(initial_composition, &vec![1.0, 1.0]);
td.composition_equations().unwrap();
td.composition_equation_sym().unwrap();
let composition_equations: &Vec<Box<dyn Fn(DVector<f64>) -> f64>> =
&td.solver.elements_conditions;
let composition_equation_sym = &td.solver.elements_conditions_sym;
let n = DVector::from_vec(vec![0.5, 0.5]);
for (i, ref comp_eq) in composition_equations.iter().enumerate() {
let eq_for_element = comp_eq(n.clone());
let eq_for_element_sym = composition_equation_sym[i].clone();
let eq_for_element_sym_value =
eq_for_element_sym.lambdify_wrapped()(vec![0.5, 0.5]);
assert_relative_eq!(eq_for_element, eq_for_element_sym_value, epsilon = 1e-6);
}
println!("composition_equations: {:?}", composition_equation_sym);
td.create_indexed_variables();
let Lamda_sym = td.solver.Lambda.clone();
let n = td.solver.n.clone();
println!("Lamda_sym: {:?}", Lamda_sym);
println!("n: {:?}", n);
assert_eq!(Lamda_sym[0], Expr::Var("Lambda0".to_string()));
assert_eq!(Lamda_sym[1], Expr::Var("Lambda1".to_string()));
assert_eq!(n[0], Expr::Var("N0".to_string()));
assert_eq!(n[1], Expr::Var("N1".to_string()));
td.set_T_to_sym();
td.create_nonlinear_system_sym().unwrap();
let eq_mu = td.solver.eq_mu.clone();
assert_eq!(eq_mu.len(), 2);
println!("eq_mu: {:?}", eq_mu);
let vars_extracted = td.symbolic_variables_extract().unwrap();
println!("vars_extracted: {:?}", vars_extracted);
td.pretty_print_substances_verbose().unwrap();
td.pretty_print_Lagrange_equations().unwrap();
td.create_nonlinear_system_fun().unwrap();
let mut sym_vars = td.solver.Lambda.clone();
sym_vars.extend(td.solver.n.clone());
let sym_vars: Vec<String> = sym_vars.iter().map(|x| x.to_string()).collect();
let mut sym_vars: Vec<&str> = sym_vars.iter().map(|x| x.as_str()).collect();
sym_vars.extend(vec!["Np"]);
println!("list of symbolic variables: {:?}", sym_vars);
let eq_s = &td.solver.eq_mu.clone();
let eq_mu_fun = &td.solver.eq_mu_fun;
let eq_mu_fun_value = eq_mu_fun(T, Some(vec![0.5, 0.5]), None, vec![0.5, 0.5]);
for (i, eq) in eq_s.clone().iter().enumerate() {
let subs = &td.vec_of_subs[i];
println!("eq_{}: {}, {:?}", subs, eq, sym_vars);
let eq_sym_fun = eq.lambdify_wrapped();
let eq_sym_fun_value = eq_sym_fun(vec![0.5, 0.5, 1.0]);
assert_relative_eq!(eq_sym_fun_value, eq_mu_fun_value[i], epsilon = 1e-6);
}
for (i, eq) in eq_s.iter().enumerate() {
let subs = &td.vec_of_subs[i];
println!("eq_{}: {}, {:?}", subs, eq, &sym_vars);
let eq_sym_fun = eq
.clone()
.lambdify_borrowed_thread_safe(sym_vars.as_slice());
let eq_sym_fun_value = eq_sym_fun(&[0.5, 0.5, 0.5, 0.5, 1.0]);
assert_relative_eq!(eq_sym_fun_value, eq_mu_fun_value[i], epsilon = 1e-6);
}
td.create_sum_of_mole_numbers_sym().unwrap();
td.form_full_system_sym().unwrap();
let full_system = &td.solver.full_system_sym;
println!("full_system: {:?}", full_system);
println!("eq_s: {:?}", eq_s);
td.pretty_print_full_system();
}
4 => {
let subs = vec![
"H2".to_string(),
"O2".to_string(),
"H2O".to_string(),
"O".to_string(),
"H".to_string(),
];
let T = 400.0;
let P = 101325.0;
let search_in_NIST = false;
let explicit_search_insructions = None;
let library_priorities = vec!["NASA_gas".to_string()];
let permitted_libraries = vec!["NUIG".to_string()];
let container = SubstancesContainer::SinglePhase(subs.clone());
let mut customsubs = SubstanceSystemFactory::create_system(
container,
None,
library_priorities,
permitted_libraries,
explicit_search_insructions,
search_in_NIST,
)
.unwrap();
let mut n = HashMap::new();
let map_of_concentration = HashMap::from([
("H2".to_string(), 0.5),
("O2".to_string(), 0.5),
("H2O".to_string(), 0.1),
]);
n.insert(None, (Some(1.0), Some(map_of_concentration)));
let td = customsubs.create_thermodynamics(Some(T), P, Some(n), None);
assert!(td.is_ok());
let mut td = td.unwrap();
td.set_P_to_sym();
td.initial_composition().unwrap();
td.create_indexed_variables();
println!("Lambda: {:#?} \n", td.solver.Lambda);
println!("n: {:#?} \n", td.solver.n);
println!("n_sym: {:#?} \n", td.solver.Np);
td.calculate_elem_composition_and_molar_mass(None);
println!(
"Initial vector of elements: {:#?} \n",
td.initial_vector_of_elements
);
println!(
"composition: {}, ncols = {}\n",
&td.clone().elem_composition_matrix.unwrap(),
&td.clone().elem_composition_matrix.unwrap().ncols()
);
td.composition_equations().unwrap();
td.composition_equation_sym().unwrap();
td.set_T_to_sym();
td.create_nonlinear_system_sym().unwrap();
td.create_nonlinear_system_fun().unwrap();
td.create_sum_of_mole_numbers_sym().unwrap();
td.form_full_system_sym().unwrap();
td.pretty_print_full_system();
}
5 => {
let gas_subs = vec![
"CO".to_string(),
"O".to_string(),
"CO2".to_string(),
"O2".to_string(),
];
let map_of_subs = HashMap::from([
("gas".to_string(), gas_subs.clone()),
("solid".to_string(), vec!["C".to_string()]),
]);
let T = 400.0;
let P = 101325.0;
let search_in_NIST = false;
let explicit_search_insructions = None;
let library_priorities = vec!["NASA_gas".to_string()];
let permitted_libraries = vec!["NUIG".to_string()];
let container = SubstancesContainer::MultiPhase(map_of_subs.clone());
let mut customsubs = SubstanceSystemFactory::create_system(
container,
None,
library_priorities,
permitted_libraries,
explicit_search_insructions,
search_in_NIST,
)
.unwrap();
println!("{:#?} \n", customsubs);
match customsubs {
CustomSubstance::PhaseOrSolution(ref phase) => {
let subs_data = &phase.subs_data;
assert!(subs_data.get(&Some("gas".to_string())).is_some());
assert!(subs_data.get(&Some("solid".to_string())).is_some());
}
_ => panic!(),
}
let mut n = HashMap::new();
let map_of_gas = HashMap::from([("O2".to_string(), 0.5)]);
let map_of_solid = HashMap::from([("C".to_string(), 0.5)]);
n.insert(Some("gas".to_string()), (Some(1.0), Some(map_of_gas)));
n.insert(Some("solid".to_string()), (Some(1.0), Some(map_of_solid)));
let td = customsubs.create_thermodynamics(Some(T), P, Some(n), None);
assert!(td.is_ok());
let mut td = td.unwrap();
td.set_P_to_sym();
println!("\n \n Lambda: {:#?} \n", td.solver.Lambda);
println!("n: {:#?} \n", td.solver.n);
println!("n_sym: {:#?} \n", td.solver.Np);
td.calculate_elem_composition_and_molar_mass(None);
println!(
"Initial vector of elements: {:#?} \n",
td.initial_vector_of_elements
);
println!(
"composition: {}, ncols = {}\n",
&td.clone().elem_composition_matrix.unwrap(),
&td.clone().elem_composition_matrix.unwrap().ncols()
);
td.composition_equations().unwrap();
td.composition_equation_sym().unwrap();
td.set_T_to_sym();
td.create_nonlinear_system_sym().unwrap();
td.create_nonlinear_system_fun().unwrap();
td.create_sum_of_mole_numbers_sym().unwrap();
td.form_full_system_sym().unwrap();
td.pretty_print_full_system();
use crate::Thermodynamics::ChemEquilibrium::ClassicalThermodynamicsSolver::{
NRParams, SolverType,
};
let params = NRParams {
initial_guess: None,
tolerance: 1e-4,
max_iterations: 100,
..Default::default()
};
td.solver.set_solver_type(SolverType::NewtonRaphson(params));
td.solver.solve();
}
6 => {
let subs = vec!["CO2".to_string(), "CO".to_string(), "O2".to_string()];
let T = 273.15;
let P = 101325.0;
let search_in_NIST = false;
let explicit_search_insructions = None;
let library_priorities = vec!["NASA_gas".to_string()];
let permitted_libraries = vec!["NUIG".to_string()];
let container = SubstancesContainer::SinglePhase(subs.clone());
let mut customsubs = SubstanceSystemFactory::create_system(
container,
None,
library_priorities,
permitted_libraries,
explicit_search_insructions,
search_in_NIST,
)
.unwrap();
let mut n = HashMap::new();
let map_of_concentration = HashMap::from([
("CO".to_string(), 0.4999),
("CO2".to_string(), 0.5),
("O2".to_string(), 0.0001),
]);
n.insert(None, (Some(1.0), Some(map_of_concentration)));
let td = customsubs.create_thermodynamics(Some(T), P, Some(n), None);
assert!(td.is_ok());
let mut td = td.unwrap();
td.set_P_to_sym();
td.initial_composition().unwrap();
td.create_indexed_variables();
println!("Lambda: {:#?} \n", td.solver.Lambda);
println!("n: {:#?} \n", td.solver.n);
println!("n_sym: {:#?} \n", td.solver.Np);
td.calculate_elem_composition_and_molar_mass(None);
td.create_full_system_of_equations_with_const_T().unwrap();
use crate::Thermodynamics::ChemEquilibrium::ClassicalThermodynamicsSolver::{
NRParams, SolverType,
};
let params = NRParams {
initial_guess: None,
tolerance: 5.0 * 1e-3,
max_iterations: 200,
damping_factor: Some(0.005),
..Default::default()
};
td.solver.set_solver_type(SolverType::NewtonRaphson(params));
td.solver.solve();
let result = td.solver.map_of_solutions;
println!("result: {:#?}", result);
}
_ => {
panic!("Invalid test case");
}
}
}