Page MenuHomec4science
Files F64176131

adimensional_system_conditions.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, May 25, 02:26

adimensional_system_conditions.cpp
View Options

#include <iostream>
#include "catch.hpp"
#include "utils/log.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
TEST_CASE("Fixed activity", "[specmicp, MiCP, program, adimensional, solver, conditions]") {
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Error;
SECTION("Carbonate speciation - fixed activity") {
specmicp::database::Database dbhandler("../data/cemdata_specmicp.js");
std::vector<std::string> to_keep = {"H2O", "H[+]", "HCO3[-]", "Na[+]"};
dbhandler.keep_only_components(to_keep);
specmicp::RawDatabasePtr thedatabase = dbhandler.get_database();
specmicp::index_t id_h = 1;
specmicp::index_t id_na = dbhandler.component_label_to_id("Na[+]");
specmicp::index_t id_hco3 = dbhandler.component_label_to_id("HCO3[-]");
specmicp::index_t id_co2 = dbhandler.aqueous_label_to_id("CO2");
specmicp::index_t id_co3 = dbhandler.aqueous_label_to_id("CO3[2-]");
specmicp::index_t id_naco3 = dbhandler.aqueous_label_to_id("NaCO3[-]");
specmicp::index_t id_nahco3 = dbhandler.aqueous_label_to_id("NaHCO3");
specmicp::index_t id_co2g = dbhandler.gas_label_to_id("CO2(g)");
specmicp::Vector total_concentrations(thedatabase->nb_component);
total_concentrations(0) = 55;
total_concentrations(id_hco3) = 0.1;
specmicp::AdimensionalSystemBC conditions(total_concentrations);
conditions.charge_keeper = id_na;
conditions.conservation_water = true;
specmicp::AdimensionalSystemSolverOptions options;
options.system_options.non_ideality = true;
options.units_set.length = specmicp::units::LengthUnit::decimeter;
options.solver_options.max_iter = 100;
options.solver_options.maxstep = 10;
options.solver_options.maxiter_maxstep = 10;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = false;
options.solver_options.penalization_factor = 0.8;
options.solver_options.factor_descent_condition = -1;
options.solver_options.fvectol = 1e-14;
specmicp::Vector x(thedatabase->nb_component+thedatabase->nb_mineral);
x << 0.8, -3.0, -2.0, -1.0;
conditions.add_fixed_activity_component(id_h, -4.0);
specmicp::AdimensionalSystemSolver solver(thedatabase, conditions, options);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x);
REQUIRE((int) perf.return_code > 0 );
specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set);
std::cout << 4.0 << " \t" << sol.pH() << " \t"
<< sol.molality_component(id_hco3) << " \t"
<< sol.molality_aqueous(id_co2) << " \t"
<< sol.molality_aqueous(id_co3) << " \t"
<< sol.molality_aqueous(id_naco3) << " \t"
<< sol.molality_aqueous(id_nahco3) << "\t"
<< sol.gas_fugacity(id_co2g)
<< std::endl;
for (double ph=4.5; ph<13; ph+=0.5)
{
conditions.fixed_activity_bc[0].log_value = -ph;
solver = specmicp::AdimensionalSystemSolver(
thedatabase, conditions, options);
perf = solver.solve(x);
REQUIRE((int) perf.return_code > 0 );
solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set);
std::cout << ph << "\t" << sol.pH() << " \t"
<< sol.molality_component(id_hco3) << " \t"
<< sol.molality_aqueous(id_co2) << " \t"
<< sol.molality_aqueous(id_co3) << "\t"
<< sol.molality_aqueous(id_naco3) << "\t"
<< sol.molality_aqueous(id_nahco3) << "\t"
<< sol.gas_fugacity(id_co2g)
<< std::endl;
}
std::cout << x << std::endl;
}
SECTION("Carbonate speciation - fixed fugacity") {
specmicp::database::Database dbhandler("../data/cemdata_specmicp.js");
std::vector<std::string> to_keep = {"H2O", "H[+]", "HCO3[-]"};
dbhandler.keep_only_components(to_keep);
specmicp::RawDatabasePtr thedatabase = dbhandler.get_database();
specmicp::index_t id_h = 1;
specmicp::index_t id_hco3 = dbhandler.component_label_to_id("HCO3[-]");
specmicp::index_t id_co2 = dbhandler.aqueous_label_to_id("CO2");
specmicp::index_t id_co3 = dbhandler.aqueous_label_to_id("CO3[2-]");;
specmicp::index_t id_co2g = dbhandler.gas_label_to_id("CO2(g)");
specmicp::Vector total_concentrations(thedatabase->nb_component);
total_concentrations(0) = 55;
total_concentrations(id_hco3) = 0.1;
specmicp::AdimensionalSystemBC conditions(total_concentrations);
conditions.charge_keeper = id_h;
conditions.conservation_water = true;
specmicp::AdimensionalSystemSolverOptions options;
options.system_options.non_ideality = true;
options.units_set.length = specmicp::units::LengthUnit::decimeter;
options.solver_options.max_iter = 100;
options.solver_options.maxstep = 10;
options.solver_options.maxiter_maxstep = 10;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = false;
options.solver_options.penalization_factor = 0.8;
options.solver_options.factor_descent_condition = -1;
options.solver_options.fvectol = 1e-14;
specmicp::Vector x(thedatabase->nb_component+thedatabase->nb_mineral);
x << 0.8, -3.0, -2.0;
conditions.add_fixed_fugacity_gas(id_co2g, id_hco3, -5);
specmicp::AdimensionalSystemSolver solver(thedatabase, conditions, options);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x);
REQUIRE((int) perf.return_code > 0 );
specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set);
std::cout << -5 << " \t" << sol.pH() << " \t"
<< sol.molality_component(id_hco3) << " \t"
<< sol.molality_aqueous(id_co2) << " \t"
<< sol.molality_aqueous(id_co3) << " \t"
<< sol.gas_fugacity(id_co2g)
<< std::endl;
for (double lfug=-4.5; lfug<=-0.5; lfug+=0.25)
{
conditions.fixed_fugacity_bc[0].log_value = lfug;
solver = specmicp::AdimensionalSystemSolver(
thedatabase, conditions, options);
perf = solver.solve(x);
REQUIRE((int) perf.return_code > 0 );
solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set);
std::cout << lfug << "\t" << sol.pH() << " \t"
<< sol.molality_component(id_hco3) << " \t"
<< sol.molality_aqueous(id_co2) << " \t"
<< sol.molality_aqueous(id_co3) << "\t"
<< sol.gas_fugacity(id_co2g)
<< std::endl;
}
// std::cout << x << std::endl;
}
}

Event Timeline

c4science · Help