carbonation.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, Jul 3, 08:24

carbonation.cpp
View Options

	// ======================================================

	// Leaching of a cement paste in CO2-saturated water

	// =======================================================

	// This file can be used to learn how to use ReactMiCP


	// Include ReactMiCP
	#include "reactmicp.hpp"

	// Other includes that would be needed later :
	#include <Eigen/QR>

	#include "physics/laws.hpp"
	#include "utils/timer.hpp"
	#include <functional>
	#include <iostream>
	#include <fstream>


	// we use the following namespaces
	using namespace specmicp;
	using namespace reactmicp;
	using namespace reactmicp::systems;

	using namespace specmicp::reactmicp::systems::satdiff;
	using namespace reactmicp::solver;

	// We declare some functions that we will be using :

	// Boundary and initial conditions
	// ===============================

	// Return the chemistry constraints for the inital condition
	AdimensionalSystemConstraints get_cement_initial_constraints(
	RawDatabasePtr the_database,
	const units::UnitsSet& units_set);
	// return the chemistry constraints for the boundary condition
	AdimensionalSystemConstraints get_bc_initial_constraints(
	RawDatabasePtr the_database,
	const units::UnitsSet& units_set);


	Vector initial_compo(const Vector& publi);

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

	// The mesh
	// ========

	// There is two version
	// a uniform mesh (cf config)
	// and a non uniform mesh
	mesh::Mesh1DPtr get_mesh();

	// ===================
	// MAIN
	// ===================

	int main()
	{
	// initialize eigen
	// This is needed if OpenMP is used
	Eigen::initParallel();

	// Setup the log
	// The output will be on cerr
	specmicp::logger::ErrFile::stream() = &std::cerr;
	// And only Warning and error messages will be printed
	specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;


	// Configuration

	auto config = specmicp::io::parse_config_file("carbonation.yaml");

	// The database
	// -------------

	// Obtain the database
	RawDatabasePtr the_database = specmicp::io::configure_database(config);

	// SpecMiCP options
	// ----------------

	// Set the options
	AdimensionalSystemSolverOptions specmicp_options;
	specmicp::io::configure_specmicp_options(specmicp_options, config);
	units::UnitsSet units_set = specmicp_options.units_set;


	// The mesh
	// --------

	// Obtain the mesh
	mesh::Mesh1DPtr the_mesh = specmicp::io::configure_mesh(config);
	index_t nb_nodes = the_mesh->nb_nodes(); // just boilerplate to make things easier

	// Print the mesh
	io::print_mesh("out_carbo_mesh.dat", the_mesh, units::LengthUnit::centimeter);


	// Initial and boundary conditions
	// --------------------------------

	// Obtain the initial constraints
	// Note : the database is modified at this point
	// Only the relevant components are kept at this point
	AdimensionalSystemConstraints cement_constraints = get_cement_initial_constraints(the_database, units_set);

	// Freeze the database
	the_database->freeze_db();


	// Obtain the initial condition
	AdimensionalSystemSolution cement_solution = solve_equilibrium(
	the_database,
	cement_constraints,
	specmicp_options);

	if(not the_database->is_valid())
	throw std::runtime_error("The database has been changed during the initial condition computation");

	database::Database db_manager(the_database);
	db_manager.save("out_carbo_db.js");


	// Just to check that everything is ok !
	AdimensionalSystemSolutionExtractor extractor(cement_solution, the_database, units_set);
	std::cout << "Porosity : " << extractor.porosity() << std::endl;
	std::cout << "Water volume fraction " << extractor.volume_fraction_water() << std::endl;



	//return EXIT_SUCCESS;

	// Obtain the boundary conditions
	AdimensionalSystemConstraints bc_constraints = get_bc_initial_constraints(
	the_database,
	specmicp_options.units_set
	);
	AdimensionalSystemSolution bc_solution = solve_equilibrium(
	the_database,
	bc_constraints,
	specmicp_options
	);

	if(not the_database->is_valid())
	throw std::runtime_error("The database has been changed during the boundary conditions computation");


	// The constraints
	// --------------

	// Boundary condition
	// 'list_fixed_nodes' lists the node with fixed composition
	std::vector<index_t> list_fixed_nodes = {0};


	// list_constraints lists the different speciation constraints in the system
	// Note : the total concentrations are computed from the initial solution,
	// and thus different total concentrations do not constitute a reason to create different constraints
	std::vector<specmicp::AdimensionalSystemConstraints> list_constraints = {cement_constraints,
	bc_constraints};

	// index_constraints links a node to the set of constraints that will be used
	std::vector<int> index_constraints(nb_nodes, 0);
	index_constraints[0] = 1; // Boundary conditions

	// This is the list of initial composition
	std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition = {cement_solution, bc_solution};
	// This list links a node to its initial conditions
	std::vector<int> index_initial_state(nb_nodes, 0);
	index_initial_state[0] = 1; // boundary condition

	// Variables
	// -----------

	// Create the main set of variables
	// They will be initialized using the list of initial composition
	satdiff::SaturatedVariablesPtr variables =
	satdiff::init_variables(the_mesh, the_database, units_set,
	list_fixed_nodes,
	list_initial_composition,
	index_initial_state);

	// Staggers
	// ---------

	// This section initialize the staggers

	// The equilibrium stagger
	// - - - - - - - - - - - -

	std::shared_ptr<EquilibriumStagger> equilibrium_stagger =
	std::make_shared<satdiff::EquilibriumStagger>(list_constraints, index_constraints, specmicp_options);

	// The transport stagger
	// - - - - - - - - - - -

	std::vector<index_t> list_fixed_component {0, 1, the_database->get_id_component("Si(OH)4")};
	std::map<index_t, index_t> slave_nodes {};


	std::shared_ptr<SaturatedTransportStagger> transport_stagger =
	std::make_shared<satdiff::SaturatedTransportStagger>(
	variables, list_fixed_nodes, slave_nodes, list_fixed_component
	);

	specmicp::io::configure_transport_options(transport_stagger->options_solver(), config);

	// The upscaling stagger
	// - - - - - - - - - - -

	CappedPowerLawParameters upscaling_param;
	upscaling_param.d_eff_0 = 2.726e-12*1e4;
	upscaling_param.cap = 1e10*1e6;
	upscaling_param.exponent = 3.32;
	upscaling_param.porosity_0 = 0.25;
	upscaling_param.porosity_res = 0.02;

	CappedPowerLaw upscaling_law(upscaling_param);

	UpscalingStaggerPtr upscaling_stagger =
	std::make_shared<DiffusiveUpscalingStagger>(upscaling_law.get_law(),the_database, units_set);

	// Solver
	// ------

	// This section initialize the reactive transport solver

	ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger);
	specmicp::io::configure_reactmicp_options(solver.get_options(), config);
	transport_stagger->initialize(variables);
	equilibrium_stagger->initialize(variables);
	upscaling_stagger->initialize(variables);


	// Some basic information about the simulation
	SimulationInformation simul_info = specmicp::io::configure_simulation_information(config);

	// Runner
	// -------

	// Create the runner : loop through the timestep

	specmicp::io::check_mandatory_yaml_node(config, "run", "__main__");
	scalar_t min_dt = specmicp::io::get_yaml_mandatory<scalar_t>(config["run"], "minimum_dt", "run");
	scalar_t max_dt = specmicp::io::get_yaml_mandatory<scalar_t>(config["run"], "maximum_dt", "run");
	scalar_t total = specmicp::io::get_yaml_mandatory<scalar_t>(config["run"], "total_execution_time", "run");
	ReactiveTransportRunner runner(solver, min_dt, max_dt, simul_info);

	// Output
	// ------

	io::OutputNodalVariables output_policy(the_database, the_mesh, specmicp_options.units_set);
	io::configure_reactmicp_output(output_policy, simul_info, config);

	// and bind it to the runner
	runner.set_output_policy(output_policy.get_output_for_reactmicp());

	// Solve the problem
	// -----------------

	runner.run_until(total, cast_var_to_base(variables));


	if(not the_database->is_valid())
	throw std::runtime_error("The database has been changed during the computation");

	AdimensionalSystemSolutionSaver saver(the_database);
	saver.save_solutions(simul_info.complete_filepath("solutions", "dat"), variables->equilibrium_solutions());

	// output information at the end of the timestep

	io::print_minerals_profile(the_database, variables, the_mesh, specmicp_options.units_set,
	simul_info.complete_filepath("profiles_end", "dat"));
	io::print_components_total_aqueous_concentration(the_database, variables,
	the_mesh, specmicp_options.units_set,
	simul_info.complete_filepath("c_profiles_end", "dat"));
	io::print_components_total_solid_concentration(the_database, variables,
	the_mesh, specmicp_options.units_set,
	simul_info.complete_filepath("s_profiles_end", "dat"));

	return EXIT_SUCCESS;
	}

	AdimensionalSystemConstraints get_cement_initial_constraints(
	RawDatabasePtr the_database,
	const units::UnitsSet& units_set)
	{

	Vector oxyde_compo(4);
	oxyde_compo << 66.98, 21.66, 7.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 = 1.0;
	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*1.0298},
	{"KOH", mult*0.0801},
	{"Cl[-]", mult*0.0001}
	};
	Dissolver dissolver(the_database);
	dissolver.set_units(units_set);

	AdimensionalSystemConstraints constraints;
	constraints.set_charge_keeper(the_database->get_id_component("HO[-]"));

	constraints.total_concentrations = dissolver.dissolve(formulation);

	std::cout << constraints.total_concentrations << std::endl;

	constraints.set_saturated_system();

	return constraints;
	}

	AdimensionalSystemConstraints get_bc_initial_constraints(
	RawDatabasePtr the_database,
	const units::UnitsSet& units_set)
	{
	const scalar_t rho_w = laws::density_water(units::celsius(25.0), units_set.length, units_set.mass);
	const scalar_t nacl = 0.3*rho_w;
	const scalar_t kcl = 0.28*rho_w;

	Vector total_concentrations;
	total_concentrations.setZero(the_database->nb_component());
	total_concentrations(the_database->get_id_component("K[+]")) = kcl;
	total_concentrations(the_database->get_id_component("Cl[-]")) = nacl + kcl;
	total_concentrations(the_database->get_id_component("Na[+]")) = nacl;
	//total_concentrations(the_database->get_id_component("Si(OH)4")) = 0.0001*rho_w;

	AdimensionalSystemConstraints constraints;
	constraints.add_fixed_activity_component(the_database->get_id_component("HO[-]"), -10.3);
	constraints.set_charge_keeper(the_database->get_id_component("HCO3[-]"));
	constraints.total_concentrations = total_concentrations;

	constraints.set_saturated_system();
	constraints.disable_surface_model();
	constraints.fixed_fugacity_cs ={};

	return constraints;
	}

	AdimensionalSystemSolution get_bc_initial_compo(
	RawDatabasePtr the_database,
	const AdimensionalSystemConstraints& constraints,
	const AdimensionalSystemSolverOptions& options
	)
	{
	Vector variables;
	AdimensionalSystemSolver solver(the_database,
	constraints,
	options);

	micpsolver::MiCPPerformance perf = solver.solve(variables, true);

	if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet)
	throw std::runtime_error("Failed to solve initial composition");

	std::cout << variables << std::endl;

	return solver.get_raw_solution(variables);
	}

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

	}


	Vector initial_compo(const Vector& publi)
	{
	Matrix publi_stochio(5,5);
	publi_stochio << 1, 0, 0, 0, 0,
	9, 6, 0, 0, 0,
	3, 0, 1, 1, 0,
	3, 0, 1, 3, 0,
	0, 0, 0, 0, 1;

	Matrix cemdata_stochio(5, 5);
	cemdata_stochio << 1, 0, 0, 0, 0,
	1.0/0.6, 1.0, 0, 0, 0,
	3, 0, 1, 1, 0,
	3, 0, 1, 3, 0,
	0, 0, 0, 0, 1;

	Vector oxyde;
	Eigen::ColPivHouseholderQR<Matrix> solver(publi_stochio);
	oxyde = solver.solve(publi);

	Vector for_cemdata = cemdata_stochio*oxyde;
	return for_cemdata;
	}

	mesh::Mesh1DPtr get_mesh()
	{
	Vector radius(71);

	const scalar_t height = 20;

	radius <<
	3.751,
	3.750,
	3.745,
	3.740,
	3.735,
	3.730,
	3.720,
	3.710,
	3.700,
	3.690,
	3.680,
	3.670,
	3.660,
	3.645,
	3.630,
	3.615,
	3.600,
	3.580,
	3.560,
	3.540,
	3.520,
	3.500,
	3.475,
	3.450,
	3.425,
	3.400,
	3.375,
	3.350,
	3.325,
	3.300,
	3.275,
	3.250,
	3.225,
	3.200,
	3.175,
	3.150,
	3.125,
	3.100,
	3.050,
	3.025,
	3.000,
	2.950,
	2.900,
	2.850,
	2.800,
	2.750,
	2.700,
	2.650,
	2.600,
	2.550,
	2.500,
	2.450,
	2.400,
	2.350,
	2.300,
	2.250,
	2.200,
	2.150,
	2.100,
	2.050,
	2.000,
	1.950,
	1.900,
	1.850,
	1.800,
	1.750,
	1.700,
	1.650,
	1.600,
	1.550,
	1.500
	;
	radius = radius/10;
	return mesh::axisymmetric_mesh1d(radius, height);
	}

carbonation.cppNo OneTemporaryActions

File Metadata

carbonation.cppView Options

Event Timeline

carbonation.cpp
No OneTemporary
Actions

carbonation.cpp
View Options