diff --git a/tests/reactmicp/binary/unsaturated/acc_carbo.yml b/tests/reactmicp/binary/unsaturated/acc_carbo.yml index f4acfe4..9f0f265 100644 --- a/tests/reactmicp/binary/unsaturated/acc_carbo.yml +++ b/tests/reactmicp/binary/unsaturated/acc_carbo.yml @@ -1,98 +1,103 @@ --- plugins: to_load: - acc_carbo_ups.so - acc_carbo_variables.so database: path: @TEST_CEMDATA_PATH@ swap_components: - out: HCO3[-] in: CO2 remove_half_cells: true units: length: centimeter mesh: type: uniform_mesh uniform_mesh: dx: 0.05 section: 5.0 nb_nodes: 25 variables: type: plugin plugin_name: acc_carbo_variables + bc_p_co2: 1000 + bc_p_h2o: 1500 + scaling_sat: 1.0e3/18.3e-3*1e-6 + scaling_aqueous: 100e-6 + boundary_conditions: gas_nodes: 0 staggers: upscaling: type: plugin plugin_name: acc_carbo_ups chemistry: type: equilibrium equilibrium_options: default: residual_tolerance: 1e-8 step_tolerance: 1e-10 maximum_step_length: 10 maximum_iterations: 200 maximum_step_maximum_iterations: 200 non_ideality_tolerance: 1e-14 cutoff_total_concentration: 1e-14 transport: # general options merge_saturation: true merge_aqueous: true # dfpm options default_options: residual_tolerance: 1e-8 step_tolerance: 1e-12 absolute_tolerance: 1e-14 quasi_newton: 3 linesearch: backtracking sparse_solver: LU saturation_options: threshold_stationary: 1e-12 maximum_step_length: 0.1 maximum_step_maximum_iterations: 500 maximum_iterations: 500 #aqueous_options: #gas_options: reactmicp: residual_tolerance: 1e-2 good_enough_tolerance: 0.99 maximum_iterations: 100 step_tolerance: 1e-14 implicit_upscaling: true simulation: name: acc_carbo output_prefix: acc_carbo_ print_iter_info: true output_step: 3600 run: run_until: 30*24*3600 timestepper: minimum_dt: 0.01 maximum_dt: 1.0 restart_dt: 0.1 output: type: hdf5 file: /tmp/acc_carbo.hdf5 diff --git a/tests/reactmicp/binary/unsaturated/acc_carbo_variables.cpp b/tests/reactmicp/binary/unsaturated/acc_carbo_variables.cpp index 8a51a49..db97e5e 100644 --- a/tests/reactmicp/binary/unsaturated/acc_carbo_variables.cpp +++ b/tests/reactmicp/binary/unsaturated/acc_carbo_variables.cpp @@ -1,231 +1,239 @@ #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 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 initialize_variables( std::shared_ptr the_database, std::shared_ptr 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; nodenb_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); + const auto p_co2 = configuration.get_required_attribute("bc_p_co2"); + variable_interface.set_partial_pressure(id_co2, 0, p_co2); + const auto p_h2o = configuration.get_required_attribute("bc_p_h2o"); + variable_interface.set_partial_pressure(0, 0, p_h2o); std::shared_ptr vars = variable_interface.get_variables(); - vars->set_aqueous_scaling(0, 1.0e3/18.3e-3*1e-6); + const auto scal_sat = configuration.get_required_attribute( + "scaling_sat"); + vars->set_aqueous_scaling(0, scal_sat); + const auto scal_aqueous = configuration.get_required_attribute( + "scaling_aqueous"); for (auto aqc: the_database->range_aqueous_component()) { - vars->set_aqueous_scaling(aqc, 100e-6); + vars->set_aqueous_scaling(aqc, scal_aqueous); } 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 solver(compo_in_oxyde); compo_species = solver.solve(compo_oxyde); }