Page MenuHomec4science
Files F91093383

hdf5_chloride.cpp
No OneTemporary
Actions

File Metadata

Created
Thu, Nov 7, 20:04

hdf5_chloride.cpp
View Options

#include "catch.hpp"
#include "specmicp/problem_solver/reactant_box.hpp"
#include "specmicp/problem_solver/smart_solver.hpp"
#include "reactmicp/io/hdf5_chloride.hpp"
#include "specmicp_common/io/hdf5_timesteps.hpp"
#include "reactmicp/reactmicp_chloride.hpp"
using namespace specmicp;
using namespace specmicp::reactmicp::systems::chloride;
static mesh::Mesh1DPtr get_mesh(index_t nb_nodes)
{
mesh::Uniform1DMeshGeometry amesh;
amesh.dx = 1.0;
amesh.section = 1.0;
amesh.nb_nodes = nb_nodes;
return mesh::uniform_mesh1d(amesh);
}
static specmicp::database::RawDatabasePtr get_database()
{
static database::RawDatabasePtr raw_data {nullptr};
if (raw_data == nullptr)
{
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
thedatabase.keep_only_components({"H2O", "H[+]", "Na[+]", "Cl[-]"});
thedatabase.swap_components({{"H[+]", "HO[-]"}});
thedatabase.remove_half_cell_reactions();
raw_data = thedatabase.get_database();
raw_data->freeze_db();
}
return raw_data;
}
static specmicp::AdimensionalSystemSolution low_conc(database::RawDatabasePtr data, units::UnitsSet tunits)
{
ReactantBox box(data, tunits);
box.set_solution(0.3, "L/L");
box.set_saturated_system();
box.set_inert_volume_fraction(0.7);
box.add_aqueous_species("NaCl", 0.001, "mol/kg");
box.add_aqueous_species("NaOH", 0.01, "mol/kg");
box.set_charge_keeper("Cl[-]");
SmartAdimSolver solver(data, box, false);
if (! solver.solve()) {
throw std::runtime_error("Cannot solve init cond");
}
return solver.get_solution();
}
static specmicp::AdimensionalSystemSolution high_conc(database::RawDatabasePtr data, units::UnitsSet tunits)
{
ReactantBox box(data, tunits);
box.set_solution(0.3, "L/L");
box.set_saturated_system();
box.set_inert_volume_fraction(0.7);
box.add_aqueous_species("NaCl", 0.5, "mol/kg");
box.add_aqueous_species("NaOH", 0.01, "mol/kg");
SmartAdimSolver solver(data, box, false);
if (! solver.solve()) {
throw std::runtime_error("Cannot solve bound cond");
}
return solver.get_solution();
}
TEST_CASE("HDF5 chloride", "[hdf5,chloride]") {
const auto nb_nodes = 15;
auto tmesh = get_mesh(nb_nodes);
auto data = get_database();
auto tunits = units::SI_units;
ChlorideVariablesInterface interface(data, tmesh, tunits);
std::vector<AdimensionalSystemSolution> sols = {low_conc(data, tunits), high_conc(data, tunits)};
std::vector<index_t> indices(nb_nodes, 0);
indices[0] = 1;
interface.set_adim_solutions(indices, sols);
index_t i_na = data->get_id_component("Na[+]");
//index_t i_cl = data->get_id_component("Cl[-]");
//index_t i_ho = data->get_id_component("HO[-]");
interface.set_diffusion_coefficient_component(i_na, 0.3);
interface.set_default_diffusion_coefficient_component(0.5);
interface.set_diffusion_coefficient_component("Cl[-]", 0.2);
interface.set_default_diffusion_resistance_factor(1.0);
ChlorideVariablesPtr vars = interface.get_variables();
SECTION ("Write") {
io::ChlorideHDF5Saver saver("test_hdf5_chloride.hdf5", vars);
saver.save_timestep(0.0);
saver.save_timestep(1.0);
}
SECTION("Read") {
io::ChlorideHDF5Reader reader("test_hdf5_chloride.hdf5", data);
auto runits = reader.get_units();
CHECK(runits.length == tunits.length);
CHECK(runits.mass == tunits.mass);
CHECK(runits.quantity == tunits.quantity);
auto& timesteps = reader.get_timesteps();
CHECK(timesteps[0] == 0.0);
CHECK(timesteps[1] == 1.0);
Vector out_0 = reader.main_variable_vs_x(timesteps(1), 0, "main_variables");
CHECK(out_0(0) == vars->total_mobile_concentration(0, 2));
CHECK(out_0(4) == vars->total_mobile_concentration(4, 2));
Vector out_1 = reader.main_variable_vs_x(timesteps(1), 1, "total_immobile_concentration");
CHECK(out_1(0) == vars->total_immobile_concentration(0, 3));
CHECK(out_1(2) == vars->total_immobile_concentration(2, 3));
auto adim_0 = reader.get_adim_solution(timesteps(0), 0);
REQUIRE(adim_0.is_valid);
CHECK(adim_0.main_variables == sols[1].main_variables);
CHECK(adim_0.ionic_strength == sols[1].ionic_strength);
auto adim_1 = reader.get_adim_solution(timesteps(0), 1);
REQUIRE(adim_1.is_valid);
CHECK(adim_1.main_variables == sols[0].main_variables);
CHECK(adim_1.ionic_strength == sols[0].ionic_strength);
}
}

Event Timeline

c4science · Help