Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65095036
acc_carbo_variables.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, May 31, 18:48
Size
10 KB
Mime Type
text/x-c++
Expires
Sun, Jun 2, 18:48 (2 d)
Engine
blob
Format
Raw Data
Handle
18003255
Attached To
rSPECMICP SpecMiCP / ReactMiCP
acc_carbo_variables.cpp
View Options
#include "reactmicp/reactmicp_unsaturated.hpp"
#include "specmicp_common/physics/unsaturated_laws.hpp"
#include "specmicp_common/plugins/plugin_interface.h"
#include "specmicp_common/plugins/plugin_base.hpp"
#include "specmicp_common/compat.hpp"
#include "specmicp_common/io/safe_config.hpp"
#include "specmicp/problem_solver/reactant_box.hpp"
#include "specmicp/problem_solver/smart_solver.hpp"
#include <iostream>
#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
);
AdimensionalSystemSolution carbonate_initial_condition(
AdimensionalSystemSolution& init_cond,
database::RawDatabasePtr the_database,
const units::UnitsSet& units_set
);
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);
AdimensionalSystemSolution init_carb = carbonate_initial_condition(
init,
the_database,
the_units);
AdimensionalSystemSolutionExtractor extr_carb(init_carb, the_database, the_units);
std::cout<< "Hop : " << init.main_variables(0) << std::endl;
std::cout<< "Hop : " << extr.saturation_water() << std::endl;
std::cout<< "Hop : " << extr.porosity() << std::endl;
//throw std::runtime_error("that's enough");
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);
variable_interface.initialize_variables(0, extr_carb);
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);
const auto p_co2 = configuration.get_required_attribute<scalar_t>("bc_p_co2");
variable_interface.set_partial_pressure(id_co2, 0, p_co2);
const auto p_h2o = configuration.get_required_attribute<scalar_t>("bc_p_h2o");
variable_interface.set_partial_pressure(0, 0, p_h2o);
std::shared_ptr<UnsaturatedVariables> vars = variable_interface.get_variables();
const auto scal_sat = configuration.get_required_attribute<scalar_t>(
"scaling_sat");
//vars->set_aqueous_scaling(0, scal_sat);
vars->set_aqueous_scaling(0, extr_carb.total_aqueous_concentration(0));
const auto scal_aqueous = configuration.get_required_attribute<scalar_t>(
"scaling_aqueous");
for (auto aqc: the_database->range_aqueous_component())
{
//vars->set_aqueous_scaling(aqc, scal_aqueous);
vars->set_aqueous_scaling(aqc, std::max(1e-3, std::abs(extr_carb.total_aqueous_concentration(aqc))));
}
const auto scal_gas = configuration.get_required_attribute<scalar_t>(
"scaling_gas");
for (auto akk: the_database->range_component())
{
vars->set_gaseous_scaling(akk, scal_gas);
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(
"reactmicp/unsaturated/variables", "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));
std::cout << "saturation " << saturation << " - pv " << pv << std::endl;
return pv;
}
AdimensionalSystemSolution get_initial_condition(
database::RawDatabasePtr the_database,
const units::UnitsSet& units_set
)
{
ReactantBox reactant_box(the_database, units_set);
reactant_box.set_solution(0.57, "kg/L");
scalar_t factor = 1.0/100*1.2;
reactant_box.add_solid_phase("CaO(ox)", factor*67.87, "kg/L");
reactant_box.add_solid_phase("SiO2(ox)", factor*22.66, "kg/L");
reactant_box.add_solid_phase("Al2O3(ox)", factor* 5.19, "kg/L");
reactant_box.add_solid_phase("SO3(ox)", factor* 2.96, "kg/L");
reactant_box.add_solid_phase("Calcite", factor*1.0, "kg/L");
//reactant_box.add_component("CO2");
Vector tot_conc = reactant_box.get_total_concentration(true);// also simplify the db
std::cout << "\n ------ \n " << tot_conc << "\n ------ \n" << std::endl;
AdimensionalSystemConstraints constraints;
constraints.set_total_concentrations(tot_conc);
//constraints.set_saturated_system();
constraints.set_inert_volume_fraction(0.0);
constraints.set_water_partial_pressure_model(
std::bind(&vapor_pressure,
1,
std::placeholders::_1)
);
//constraints.charge_keeper = the_database->get_id_component_from_element("H");
//constraints.set_inert_volume_fraction(0.1);
AdimensionalSystemSolverOptions opts;
opts.units_set = units_set;
opts.system_options.non_ideality_tolerance = 1e-14;
opts.system_options.restart_water_volume_fraction = 0.8;
opts.system_options.cutoff_total_concentration = 1e-14;
opts.solver_options.set_tolerance(1e-10, 1e-12);
//micpsolver::MiCPSolverOptions& solver_options = opts.solver_options;
//solver_options.maxstep = 10.0;
//solver_options.maxiter_maxstep = 200;
//solver_options.max_iter = 200;
//solver_options.set_tolerance(1e-6, 1e-12);
SmartAdimSolver smart_solver(the_database, constraints, opts);
smart_solver.solve();
if (not smart_solver.is_solved())
{
throw std::runtime_error("Error : failed to solve initial conditions");
}
AdimensionalSystemSolution sol = smart_solver.get_solution();
// std::cout << "normal solver " << std::endl;
// AdimensionalSystemSolver solver(the_database, constraints, opts);
// Vector x;
// solver.set_options(opts);
// solver.initialise_variables(x, 0.5, -6);
// micpsolver::MiCPPerformance perf = solver.solve(x, false);
// std::cout << perf.current_residual << " - "
// << perf.current_update << " - "
// << perf.nb_iterations
// << std::endl;
// if (perf.return_code <= micpsolver::MiCPSolverReturnCode::Success) {
// throw std::runtime_error("Error : failed to solve initial conditions");
// }
// AdimensionalSystemSolution sol = solver.get_raw_solution(x);
std::cout << sol.main_variables << std::endl;
AdimensionalSystemSolutionExtractor extr(sol, the_database, units_set);
for (index_t idc: the_database->range_component())
{
std::cout << the_database->get_label_component(idc) << " : "
<< extr.total_concentration(idc) << " - "
<< tot_conc(idc) << " - "
<< (extr.total_concentration(idc) - tot_conc(idc))/tot_conc(idc) << std::endl;
}
//throw std::runtime_error("That's enough");
return sol;
}
AdimensionalSystemSolution carbonate_initial_condition(
AdimensionalSystemSolution& init_cond,
database::RawDatabasePtr the_database,
const units::UnitsSet& units_set
)
{
AdimensionalSystemSolutionExtractor extr(init_cond, the_database, units_set);
scalar_t init_porosity = extr.porosity();
AdimensionalSystemConstraints constraints;
Vector tot_conc = extr.total_concentrations();
index_t id_ca = the_database->get_id_component_from_element("Ca");
index_t id_co2 = the_database->get_id_component_from_element("C");
tot_conc(id_co2) = tot_conc(id_ca);
constraints.set_total_concentrations(tot_conc);
constraints.set_saturated_system();
AdimensionalSystemSolverOptions opts;
opts.units_set = units_set;
opts.system_options.non_ideality_tolerance = 1e-14;
opts.system_options.restart_water_volume_fraction = 0.8;
opts.system_options.cutoff_total_concentration = 1e-14;
opts.solver_options.set_tolerance(1e-10, 1e-12);
AdimensionalSystemSolver solver(the_database, constraints, init_cond, opts);
Vector x = init_cond.main_variables;
auto perf = solver.solve(x, false);
if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
{
throw std::runtime_error("Failed to solve carbonated initial condition");
}
AdimensionalSystemSolution solution = solver.get_raw_solution(x);
// set the same saturation for both cases
AdimensionalSystemSolutionExtractor ext2r(solution, the_database, units_set);
scalar_t porosity = ext2r.porosity();
solution.main_variables(0) = porosity/init_porosity * init_cond.main_variables(0);
return solution;
}
Event Timeline
Log In to Comment