Diffusion SpecMiCP / ReactMiCP (dev) Browse

acc_carbo.cpp
dee7f564bf24dev

acc_carbo.cpp
View Options

			/* =============================================================================

			Copyright (c) 2014 - 2016
			F. Georget <fabieng@princeton.edu> Princeton University
			All rights reserved.


			Redistribution and use in source and binary forms, with or without modification,
			are permitted provided that the following conditions are met:

			1. Redistributions of source code must retain the above copyright notice,
			this list of conditions and the following disclaimer.

			2. Redistributions in binary form must reproduce the above copyright notice,
			this list of conditions and the following disclaimer in the documentation
			and/or other materials provided with the distribution.

			3. Neither the name of the copyright holder nor the names of its
			contributors may be used to endorse or promote products derived from this
			software without specific prior written permission.

			* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
			AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
			IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
			DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
			ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
			(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
			LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
			ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
			(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
			SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *

			============================================================================= */

			//! \file acc_carbo.cpp
			//! \brief Accelerated carbonation example
			//!
			//! This is an example to show the features of the unsaturated
			//! reactive transport solver
			//! It is a simple study of the accelerated carbonation of cement paste
			//! (with drying)
			//!
			//!

			#include "reactmicp/reactmicp_unsaturated.hpp"

			#include "reactmicp/systems/unsaturated/variables.hpp"

			#include "reactmicp/io/hdf5_unsaturated.hpp"
			#include "dfpm/io/hdf5_mesh.hpp"
			#include "specmicp_database/io/hdf5_database.hpp"

			#include "dfpm/solver/parabolic_structs.hpp"

			#include "specmicp_common/physics/unsaturated_laws.hpp"

			#include "specmicp_common/cli/parser.hpp"


			#include <Eigen/QR>
			#include <iostream>

			using namespace specmicp;
			using namespace specmicp::reactmicp;
			using namespace specmicp::reactmicp::systems::unsaturated;

			using VariablesBase = solver::VariablesBase;
			using StaggerReturnCode = solver::StaggerReturnCode;


			class UpscalingStagger;

			database::RawDatabasePtr get_database();

			AdimensionalSystemSolution get_initial_condition(
			database::RawDatabasePtr the_database,
			const units::UnitsSet& units_set
			);


			void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species);

			// Upscaling stagger -> user models

			class UpscalingStagger: public solver::UpscalingStaggerBase
			{
			public:
			UpscalingStagger() {}

			~UpscalingStagger() {}

			//! \brief Initialize the stagger at the beginning of the computation
			//!
			//! \param var a shared_ptr to the variables
			virtual void initialize(VariablesBase * const var) override;

			//! \brief Initialize the stagger at the beginning of an iteration
			//!
			//! This is where the predictor can be saved, the first trivial iteration done, ...
			//!
			//! \param dt the duration of the timestep
			//! \param var a shared_ptr to the variables
			virtual void initialize_timestep(
			scalar_t dt,
			VariablesBase * const var) override;

			//! \brief Solve the equation for the timestep
			//!
			//! \param var a shared_ptr to the variables
			virtual StaggerReturnCode restart_timestep(
			VariablesBase * const var) override;

			// User models
			scalar_t capillary_pressure(index_t node, scalar_t saturation);

			scalar_t vapor_pressure(index_t node, scalar_t saturation);

			scalar_t relative_liquid_diffusivity(index_t node, scalar_t saturation);

			scalar_t relative_liquid_permeability(index_t node, scalar_t saturation);

			scalar_t relative_gas_diffusivity(index_t node, scalar_t saturation);

			scalar_t intrinsic_permeability(scalar_t porosity);

			scalar_t intrinsic_liquid_diffusivity(scalar_t porosity);

			scalar_t resistance_gas_diffusivity(scalar_t porosity);

			private:
			scalar_t m_a {37.5479e6};
			scalar_t m_b {2.1684};
			scalar_t m_alpha {7.27};
			scalar_t m_beta {0.440};
			scalar_t m_pv_sat {3170};
			scalar_t m_exp_r_l_d {3.00};
			scalar_t m_exp_i_g_d {2.74};
			scalar_t m_exp_r_g_d {4.20};
			scalar_t m_exp_r_l_k {0.4191};

			scalar_t m_k_0 {1e-21*1e4};
			scalar_t m_d_0 {2.3e-13*1e4};
			scalar_t m_phi_0 {0.2};
			};

			// implementations

			int main(int argc, char* argv[])
			{
			cli::CommandLineParser cli_parser;
			cli_parser.add_option('n', "min_dt", 0.01, "minimum timestep (s)");
			cli_parser.add_option('m', "max_dt", 10.0, "maximum timestep (s)");
			cli_parser.add_option('d', "dx", 0.05, "space step (cm)");
			cli_parser.add_option('r', "run_until", 10.024.03600.0, "run simulation until <value> s");
			cli_parser.add_option('p', "p_h20", 1000.0, "vapor pressure at the boundary");
			cli_parser.add_option('c', "p_co2", 500.0, "vapor pressure at the boundary");

			cli_parser.set_help_message("Atmospheric carbonation");

			cli_parser.parse(argc, argv);

			init_logger(&std::cout, logger::Warning);

			index_t nb_nodes = 25;
			scalar_t dx = cli_parser.get_option<cli::ValueType::floating>("dx"); // cm
			scalar_t cross_section = 5; // cm^2

			units::UnitsSet units_set;
			units_set.length = units::LengthUnit::centimeter;

			specmicp::RawDatabasePtr the_database = get_database();

			AdimensionalSystemSolution init = get_initial_condition(
			the_database,
			units_set);

			AdimensionalSystemSolutionExtractor extr(init, the_database, units_set);

			std::cout << extr.porosity() << " - " << extr.saturation_water() << std::endl;

			//return EXIT_SUCCESS;

			mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d({dx, nb_nodes, cross_section});

			index_t id_co2 = the_database->get_id_component("CO2");

			VariablesInterface variable_interface(the_mesh, the_database, {0, id_co2}, units_set);

			variable_interface.set_binary_gas_diffusivity(0, 0.240);
			variable_interface.set_binary_gas_diffusivity(id_co2, 0.160);

			std::shared_ptr<UpscalingStagger> upscaling_stagger =
			std::make_shared<UpscalingStagger>();

			auto* call_ptr = upscaling_stagger.get();
			variable_interface.set_vapor_pressure_model(
			[call_ptr](index_t node, scalar_t sat) -> scalar_t {
			return call_ptr->vapor_pressure(node, sat);
			});

			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,
			cli_parser.get_option<cli::ValueType::floating>("p_co2"));
			variable_interface.set_partial_pressure(0, 0,
			cli_parser.get_option<cli::ValueType::floating>("p_h20"));

			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);
			}
			//vars->set_gaseous_scaling(0, 1.0/vars->get_rt());
			//vars->set_gaseous_scaling(id_co2, 1.0/vars->get_rt());

			std::cout << "Youpla : " << vars->get_rt() << std::endl;

			auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), the_database->nb_component());
			bcs->add_gas_node(0);
			bcs->fork_constraint(0, "gas_node");

			auto chem_options = EquilibriumOptionsVector::make(the_mesh->nb_nodes());
			auto& def_chem_opts = chem_options->get("default");
			def_chem_opts.units_set = units_set;

			micpsolver::MiCPSolverOptions& copts = def_chem_opts.solver_options;
			copts.set_maximum_iterations(200);
			copts.set_maximum_step_length(10, 200);
			copts.set_tolerance(1e-8, 1e-12);
			def_chem_opts.system_options.non_ideality_tolerance = 1e-14;


			std::shared_ptr<EquilibriumStagger> chemistry_stagger =
			EquilibriumStagger::make(vars, bcs, chem_options);

			upscaling_stagger->initialize(vars.get());
			vars->set_relative_variables(0);
			chemistry_stagger->initialize(vars.get());

			std::shared_ptr<UnsaturatedTransportStagger> transport_stagger =
			UnsaturatedTransportStagger::make(vars, bcs);

			transport_stagger->initialize(vars.get());

			bcs->get_constraint("gas_node").set_water_partial_pressure_model(
			std::bind(&UpscalingStagger::vapor_pressure,
			upscaling_stagger.get(),
			1,
			std::placeholders::_1)
			);

			auto check = variable_interface.check_variables();
			if (check >= VariablesValidity::error) {
			std::cout << variable_interface.get_liquid_saturation() << std::endl;
			throw std::runtime_error("Invalid variables");
			}

			solver::ReactiveTransportSolver react_solver(transport_stagger, chemistry_stagger, upscaling_stagger);
			solver::ReactiveTransportOptions& opts = react_solver.get_options();
			opts.residuals_tolerance = 1e-2;
			opts.good_enough_tolerance = 0.99;
			opts.maximum_iterations = 100;
			transport_stagger->get_saturation_options()->maximum_iterations = 500;
			transport_stagger->get_saturation_options()->step_tolerance = 1e-10;
			transport_stagger->get_saturation_options()->residuals_tolerance = 1e-8;

			for (auto ind: the_database->range_aqueous_component())
			{
			auto* opts = transport_stagger->get_aqueous_options(ind);
			opts->step_tolerance = 1e-10;
			opts->residuals_tolerance = 1e-8;
			}

			//transport_stagger->get_gas_options(3)->maximum_iterations = 5000;
			opts.step_tolerance = 1e-10;

			io::UnsaturatedHDF5Saver saver("/tmp/out_acc_carbo.hdf5", vars, units_set);
			saver.save_timestep(0.0);

			solver::SimulationInformation simul_info("acc_carbo", 500);

			auto out_pol = [&saver](
			scalar_t timestep,
			solver::VariablesBasePtr _) {
			saver.save_timestep(timestep);
			};

			database::Database db_handler(the_database);
			db_handler.save("acc_carbo_data.yml");

			//solver::ReactiveTransportReturnCode retcode = react_solver.solve_timestep(0.01, vars);
			//retcode = react_solver.solve_timestep(1.0, vars);

			//std::cout << (int) retcode << std::endl;
			// retcode = react_solver.solve_timestep(5.0, vars);
			// std::cout << (int) retcode << std::endl;
			// retcode = react_solver.solve_timestep(10.0, vars);
			// std::cout << (int) retcode << std::endl;
			// retcode = react_solver.solve_timestep(10.0, vars);
			// std::cout << (int) retcode << std::endl;

			solver::ReactiveTransportRunner runner(
			react_solver,
			cli_parser.get_option<cli::ValueType::floating>("min_dt"),
			cli_parser.get_option<cli::ValueType::floating>("max_dt"),
			simul_info);
			runner.set_output_policy(out_pol);


			runner.run_until(cli_parser.get_option<cli::ValueType::floating>("run_until"),
			vars);

			std::cout << vars->get_liquid_saturation().variable << std::endl;
			std::cout << vars->get_liquid_saturation()(1) - vars->get_liquid_saturation()(2) << std::endl;

			std::cout << vars->get_pressure_main_variables(id_co2).variable << std::endl;
			}

			void UpscalingStagger::initialize(VariablesBase * const var)
			{
			UnsaturatedVariables * const true_var = static_cast<UnsaturatedVariables * const>(var);

			// set relative models

			true_var->set_capillary_pressure_model(
			[this](index_t x, scalar_t sat){
			return this->capillary_pressure(x, sat);
			});
			true_var->set_relative_liquid_permeability_model(
			[this](index_t x, scalar_t sat){
			return this->relative_liquid_permeability(x, sat);
			});
			true_var->set_relative_liquid_diffusivity_model(
			[this](index_t x, scalar_t sat){
			return this->relative_liquid_diffusivity(x, sat);
			});
			true_var->set_relative_gas_diffusivity_model(
			[this](index_t x, scalar_t sat){
			return this->relative_gas_diffusivity(x, sat);
			});
			true_var->set_vapor_pressure_model(
			[this](index_t x, scalar_t sat){
			return this->vapor_pressure(x, sat);
			});

			// initialization

			SecondaryTransientVariable& porosity = true_var->get_porosity();
			SecondaryVariable& permeability = true_var->get_liquid_permeability();
			SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity();
			SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity();

			true_var->set_relative_variables();
			gas_diffusivity(0) = resistance_gas_diffusivity(1.0);
			true_var->get_relative_gas_diffusivity()(0) = relative_gas_diffusivity(0, 1.0);

			for (index_t node=1; node<true_var->get_mesh()->nb_nodes(); ++node)
			{
			scalar_t phi = porosity(node);

			permeability(node) = intrinsic_permeability(phi);
			liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi);
			gas_diffusivity(node) = resistance_gas_diffusivity(phi);

			//true_var->set_relative_variables(node);
			}
			}

			scalar_t UpscalingStagger::capillary_pressure(index_t node, scalar_t saturation)
			{
			if (saturation == 0.0) return 0.0;

			return laws::capillary_pressure_BaroghelBouny(saturation, m_a, m_b);
			}

			scalar_t UpscalingStagger::vapor_pressure(index_t node, scalar_t saturation)
			{
			scalar_t pv;
			if (saturation == 0.0) pv = 0.0;
			else if (saturation >= 1.0) pv = m_pv_sat;
			else pv = m_pv_sat*laws::invert_kelvin_equation(
			capillary_pressure(node, saturation));
			return pv;
			}

			scalar_t UpscalingStagger::relative_liquid_diffusivity(
			index_t node,
			scalar_t saturation
			)
			{
			if (saturation == 0.0) return 0.0;
			return std::pow(saturation, m_exp_r_l_d);
			}

			scalar_t UpscalingStagger::relative_liquid_permeability(
			index_t node,
			scalar_t saturation
			)
			{
			if (saturation == 0.0) return 0.0;

			return laws::relative_liquid_permeability_Mualem(saturation, 1.0/m_b);
			}

			scalar_t UpscalingStagger::relative_gas_diffusivity(index_t node, scalar_t saturation)
			{
			if (saturation == 0.0) return 1.0;

			return std::pow(1.0 - saturation, m_exp_r_g_d);
			}

			scalar_t UpscalingStagger::intrinsic_permeability(scalar_t porosity)
			{
			scalar_t tmp_1 = porosity / m_phi_0;
			scalar_t tmp_2 = (1.0 - m_phi_0) / (1.0 - porosity);
			tmp_1 = std::pow(tmp_1, 3);
			tmp_2 = std::pow(tmp_2, 2);
			return m_k_0 * tmp_1 * tmp_2;
			}

			scalar_t UpscalingStagger::intrinsic_liquid_diffusivity(scalar_t porosity)
			{
			return m_d_0 * std::exp(-9.95*porosity);
			}

			scalar_t UpscalingStagger::resistance_gas_diffusivity(scalar_t porosity)
			{
			return std::pow(porosity, m_exp_i_g_d);
			}

			void UpscalingStagger::initialize_timestep(
			scalar_t dt,
			VariablesBase * const var
			)
			{
			}

			StaggerReturnCode UpscalingStagger::restart_timestep(VariablesBase * const var)
			{
			UnsaturatedVariables* true_var = static_cast<UnsaturatedVariables*>(var);

			SecondaryTransientVariable& porosity = true_var->get_porosity();
			SecondaryVariable& permeability = true_var->get_liquid_permeability();
			SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity();
			SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity();


			for (index_t node=1; node<true_var->get_mesh()->nb_nodes(); ++node)
			{
			scalar_t phi = porosity(node);

			permeability(node) = intrinsic_permeability(phi);
			liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi);
			gas_diffusivity(node) = resistance_gas_diffusivity(phi);

			true_var->set_relative_variables(node);
			}

			return StaggerReturnCode::ResidualMinimized;
			}

			specmicp::RawDatabasePtr get_database()
			{
			specmicp::database::Database thedatabase(CEMDATA_PATH);
			std::map<std::string, std::string> swap = {{"HCO3[-]", "CO2"},};
			thedatabase.swap_components(swap);
			thedatabase.remove_half_cell_reactions();
			database::RawDatabasePtr raw_data = thedatabase.get_database();
			return raw_data;
			}

			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);

			UpscalingStagger upsstag_init; // just for vapor pressure model

			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.initialize_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(&UpscalingStagger::vapor_pressure,
			&upsstag_init,
			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);

			}

specmicp/examples/reactmicp/unsaturated/acc_carbo.cppdee7f564bf24dev

acc_carbo.cppView Options

specmicp/examples/reactmicp/unsaturated/acc_carbo.cpp
dee7f564bf24dev

acc_carbo.cpp
View Options