Homec4science
Diffusion SpecMiCP / ReactMiCP (yaml/safe) Browse

specmicp/tests/reactmicp/binary/unsaturated/acc_carbo_variables.cpp
efe3e6ff717dyaml/safe

Last Change
Raw File

acc_carbo_variables.cpp
View Options

#include "reactmicp_unsaturated.hpp"
#include "physics/unsaturated_laws.hpp"
#include "utils/plugins/plugin_interface.h"
#include "utils/plugins/plugin_base.hpp"
#include "utils/compat.hpp"
#include <Eigen/QR>
using namespace specmicp;
using namespace specmicp::reactmicp::solver;
using namespace specmicp::reactmicp::systems::unsaturated;
AdimensionalSystemSolution get_initial_condition(
database::RawDatabasePtr the_database,
const units::UnitsSet& units_set
);
void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species);
scalar_t vapor_pressure(index_t node, scalar_t saturation);
class SPECMICP_DLL_PUBLIC AccCarboVariablesFactory: public InitializeVariables
{
public:
static AccCarboVariablesFactory* get() {
return new AccCarboVariablesFactory();
}
virtual std::shared_ptr<UnsaturatedVariables> initialize_variables(
std::shared_ptr<database::DataContainer> the_database,
std::shared_ptr<mesh::Mesh1D> the_mesh,
const units::UnitsSet& the_units,
io::YAMLConfigHandle& configuration
) override
{
AdimensionalSystemSolution init = get_initial_condition(
the_database,
the_units);
AdimensionalSystemSolutionExtractor extr(init, the_database, the_units);
index_t id_co2 = the_database->get_id_component("CO2");
VariablesInterface variable_interface(the_mesh, the_database, {0, id_co2}, the_units);
variable_interface.set_binary_gas_diffusivity(0, 0.240);
variable_interface.set_binary_gas_diffusivity(id_co2, 0.160);
variable_interface.set_vapor_pressure_model(&vapor_pressure);
for (index_t node=1; node<the_mesh->nb_nodes(); ++node)
{
variable_interface.initialize_variables(node, extr);
}
variable_interface.set_porosity(0, 1.0);
variable_interface.set_liquid_saturation(0, 0.0);
variable_interface.set_partial_pressure(id_co2, 0, 1000);
variable_interface.set_partial_pressure(0, 0, 1000);
std::shared_ptr<UnsaturatedVariables> vars = variable_interface.get_variables();
vars->set_aqueous_scaling(0, 1.0e3/18.3e-3*1e-6);
for (auto aqc: the_database->range_aqueous_component())
{
vars->set_aqueous_scaling(aqc, 100e-6);
}
return vars;
}
};
class SPECMICP_DLL_PUBLIC AccCarboVariablesPlugin: public specmicp::plugins::PluginBase
{
public:
AccCarboVariablesPlugin():
PluginBase()
{
set_api_version(0, 1, 0);
}
bool initialize(const specmicp::plugins::PluginManagerServices& services) override
{
auto retcode = services.register_object(
"variables_factory", "acc_carbo_variables",
&AccCarboVariablesFactory::get
);
if (not retcode) {
return false;
}
return true;
}
};
SPECMICP_PLUGIN(AccCarboVariablesPlugin)
// implementation
scalar_t capillary_pressure(index_t node, scalar_t saturation)
{
if (saturation == 0.0) return 0.0;
return laws::capillary_pressure_BaroghelBouny(saturation, 37.5479e6, 2.1684);
}
scalar_t vapor_pressure(index_t node, scalar_t saturation)
{
scalar_t pv;
if (saturation == 0.0) pv = 0.0;
else if (saturation >= 1.0) pv = 3170;
else pv = 3170*laws::invert_kelvin_equation(
capillary_pressure(node, saturation));
return pv;
}
AdimensionalSystemSolution get_initial_condition(
database::RawDatabasePtr the_database,
const units::UnitsSet& units_set
)
{
Vector oxyde_compo(4);
oxyde_compo << 67.98, 22.66, 5.19, 2.96;
Vector species_compo;
compo_from_oxyde(oxyde_compo, species_compo);
scalar_t mult = 1e-6;
//species_compo *= 0.947e-2; // poro 0.47
//species_compo *= 1.08e-2; // poro 0.4
species_compo *= 1.25e-2; // poro 0.3
//species_compo *= 1.1e-2;
Formulation formulation;
formulation.mass_solution = 5e-4;
formulation.amount_minerals = {
{"C3S", species_compo(0)},
{"C2S", species_compo(1)},
{"C3A", species_compo(2)},
{"Gypsum", species_compo(3)},
{"Calcite", species_compo(0)*1e-3}
};
//formulation.extra_components_to_keep = {"HCO3[-]"};
formulation.concentration_aqueous =
{
{"NaOH", mult*0.0298},
{"KOH", mult*0.0801},
{"Cl[-]", mult*0.0001}
};
Dissolver dissolver(the_database);
AdimensionalSystemConstraints constraints;
constraints.set_total_concentrations(
dissolver.dissolve(formulation, true)
);
//constraints.set_inert_volume_fraction(0.1);
AdimensionalSystemSolver the_solver(the_database, constraints);
the_solver.get_options().units_set = units_set;
the_solver.get_options().system_options.non_ideality_tolerance = 1e-12;
micpsolver::MiCPSolverOptions* solver_options = &the_solver.get_options().solver_options;
//solver_options.maxstep = 10.0;
solver_options->maxiter_maxstep = 200;
solver_options->max_iter = 200;
solver_options->set_tolerance(1e-9);
Vector x;
the_solver.initialise_variables(x, 0.8, -4);
micpsolver::MiCPPerformance perf = the_solver.solve(x);
if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
{
throw std::runtime_error("Error : failed to solve initial conditions");
}
auto prev_sol = the_solver.get_raw_solution(x);
prev_sol.main_variables(0) *= 0.8;
constraints.set_water_partial_pressure_model(
std::bind(&vapor_pressure,
1,
std::placeholders::_1)
);
the_solver = AdimensionalSystemSolver (the_database, constraints, prev_sol);
the_solver.get_options().units_set = units_set;
the_solver.get_options().system_options.non_ideality_tolerance = 1e-12;
solver_options = &the_solver.get_options().solver_options;
//solver_options.maxstep = 10.0;
solver_options->maxiter_maxstep = 200;
solver_options->set_tolerance(1e-9);
perf = the_solver.solve(x);
if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
{
throw std::runtime_error("Error : failed to solve initial conditions");
}
return the_solver.get_raw_solution(x);
}
void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species)
{
constexpr double M_CaO = 56.08;
constexpr double M_SiO2 = 60.09;
constexpr double M_Al2O3 = 101.96;
constexpr double M_SO3 = 80.06;
Eigen::MatrixXd compo_in_oxyde(4, 4);
Vector molar_mass(4);
molar_mass << 1.0/M_CaO, 1.0/M_SiO2, 1.0/M_Al2O3, 1.0/M_SO3;
compo_oxyde = compo_oxyde.cwiseProduct(molar_mass);
//C3S C2S C3A Gypsum
compo_in_oxyde << 3, 2, 3, 1, // CaO
1, 1, 0, 0, // SiO2
0, 0, 1, 0, // Al2O3
0, 0, 0, 1; // SO3
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> solver(compo_in_oxyde);
compo_species = solver.solve(compo_oxyde);
}
c4science · Help