Page MenuHomec4science
Files F60175922

hdf5_unsaturated.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, Apr 28, 02:10

hdf5_unsaturated.cpp
View Options

#include "catch.hpp"
#include "reactmicp/io/hdf5_unsaturated.hpp"
#include "specmicp/adimensional/adimensional_system_solution.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "specmicp_database/database.hpp"
#include "specmicp/problem_solver/reactant_box.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "reactmicp/systems/unsaturated/variables_interface.hpp"
#include "dfpm/mesh.hpp"
#include <string>
#define PATH_TEST_FILES "test_hdf5_unsaturated"
using namespace specmicp;
using namespace specmicp::reactmicp::systems::unsaturated;
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[+]", "Ca[2+]", "Si(OH)4", "HCO3[-]"}
);
std::map<std::string, std::string> swap = {{"HCO3[-]", "CO2"},};
thedatabase.swap_components(swap);
raw_data = thedatabase.get_database();
raw_data->freeze_db();
}
return raw_data;
}
static AdimensionalSystemSolution
get_init_solution(
database::RawDatabasePtr the_database,
units::UnitsSet the_units
)
{
specmicp::ReactantBox reactant_box(the_database, the_units);
reactant_box.set_solution(0.5, "L/L");
reactant_box.add_solid_phase("C3S", 0.350, "L/L");
reactant_box.add_solid_phase("C2S", 0.50, "L/L");
reactant_box.add_aqueous_species("CO2", 32.72, "mol/kg");
reactant_box.set_charge_keeper("H");
auto constraints = reactant_box.get_constraints(true);
AdimensionalSystemSolver adim_solver(the_database, constraints);
adim_solver.get_options().solver_options.set_tolerance(1e-14);
adim_solver.get_options().units_set = the_units;
Vector x;
adim_solver.initialize_variables(x, 0.8, -3);
micpsolver::MiCPPerformance perf = adim_solver.solve(x);
if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
{
throw std::runtime_error("Error solving the initial problem");
}
return adim_solver.get_raw_solution(x);
}
TEST_CASE("hdf5 output", "[hdf5]") {
std::string filename = "test_hdf5_unsaturated.hdf5";
database::RawDatabasePtr raw_data = get_database();
mesh::Uniform1DMeshGeometry mesh_geom;
mesh_geom.nb_nodes = 10;
mesh_geom.dx = 1;
mesh_geom.section = 1;
mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d(mesh_geom);
units::UnitsSet the_units;
the_units.length = units::LengthUnit::centimeter;
specmicp::AdimensionalSystemSolution init_sol
= get_init_solution(raw_data, the_units);
specmicp::AdimensionalSystemSolutionExtractor extr(init_sol, raw_data, the_units);
SECTION("Writer") {
VariablesInterface vars_init(the_mesh, raw_data,
{raw_data->get_id_component_from_element("C")}, the_units);
for (auto node: the_mesh->range_nodes())
{
vars_init.initialize_variables(node, extr);
}
auto vars = vars_init.get_variables();
io::UnsaturatedHDF5Saver saver(filename, vars, the_units);
saver.save_timestep(10.0);
}
SECTION("Reader") {
specmicp::io::UnsaturatedHDF5Reader reader(filename);
// units
units::UnitsSet read_unit = reader.get_units();
CHECK(read_unit.length == the_units.length);
CHECK(read_unit.mass == the_units.mass);
CHECK(read_unit.quantity == the_units.quantity);
// adimensional solution
specmicp::AdimensionalSystemSolution sol = reader.get_adim_solution(10.0, 1);
REQUIRE(sol.is_valid);
for (auto r: range(sol.main_variables.rows()))
{
if (std::isfinite(sol.main_variables(r)))
{
CHECK(sol.main_variables(r) == Approx(init_sol.main_variables(r)).epsilon(1e-10));
}
}
CHECK(sol.log_gamma.norm() == Approx(init_sol.log_gamma.norm()).epsilon(1e-10));
CHECK(sol.gas_fugacities.norm() == Approx(init_sol.gas_fugacities.norm()).epsilon(1e-10));
CHECK(sol.secondary_molalities.norm() == Approx(init_sol.secondary_molalities.norm()).epsilon(1e-10));
// main variables
Vector liq_sat = reader.main_variable_vs_x(10.0, "H2O", "liquid_saturation");
CHECK(liq_sat.rows() == the_mesh->nb_nodes());
CHECK(liq_sat(1) == Approx(extr.saturation_water()));
Vector ca_s = reader.main_variable_vs_x(10.0, "Ca[2+]", "solid_concentration");
CHECK(ca_s.rows() == the_mesh->nb_nodes());
CHECK(ca_s(1) == Approx(extr.total_solid_concentration(raw_data->get_id_component("Ca[2+]"))));
}
}

Event Timeline

c4science · Help