carbonationfe.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Tue, Nov 19, 09:53

carbonationfe.cpp
View Options

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

	// Saturated carbonation of a cement past

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

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


	// Include ReactMiCP
	#include "reactmicp.hpp"

	// Physical laws
	#include "physics/laws.hpp"
	// A simple timer
	#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 :

	// Database
	// ========

	// Parse and Prepare the database
	RawDatabasePtr get_database(const std::string& path_to_database);

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

	void compo_oxyde_to_mole(Vector& compo_oxyde);

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

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

	// Options
	// =======

	// Speciation solver options
	AdimensionalSystemSolverOptions get_specmicp_options();
	// Transport options
	void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options);
	// Reactive Solver options
	void set_reactive_solver_options(reactmicp::solver::ReactiveTransportOptions& reactive_options);

	// Upscaling
	// =========

	// Upscaling is the stagger in charge of finding the porosity and transport properties
	// This is a class that is dependant on the problem to solve and should be defined by the user
	// The following class is an example on how to define such a class

	// The Upscaling class should inherit from UpscalingStaggerBase
	class CarboUspcaling: public solver::UpscalingStaggerBase
	{
	public:
	// Initiation
	// This is not call automatically so can be adjusted by the user
	CarboUspcaling(RawDatabasePtr the_database,
	units::UnitsSet the_units):
	m_data(the_database),
	m_units_set(the_units),
	m_dt(HUGE_VAL)
	{
	}


	// Initialize the stagger at the beginning of the computation
	// The porosity and transport properties should be initialized here
	void initialize(VariablesBasePtr var)
	{
	SaturatedVariablesPtr true_var = cast_var_from_base(var);
	true_var->upscaling_variables().setZero();
	for (index_t node=0; node<true_var->nb_nodes(); ++node)
	{
	upscaling_one_node(node, true_var);
	}
	}

	// Initialize the stagger at the beginning of a timestep
	// Typically do nothing but erase the "dot" variable (velocity such as the vel_porosity)
	void initialize_timestep(scalar_t dt, VariablesBasePtr var)
	{
	SaturatedVariablesPtr true_var = cast_var_from_base(var);
	m_dt = dt;

	for (index_t node=1; node<true_var->nb_nodes(); ++node)
	{
	upscaling_one_node(node, true_var);
	true_var->vel_porosity(node) = 0.0;
	}
	}

	//! This is the main function called during a timestep
	StaggerReturnCode restart_timestep(VariablesBasePtr var)
	{
	SaturatedVariablesPtr true_var = cast_var_from_base(var);
	for (index_t node=1; node<true_var->nb_nodes(); ++node)
	{
	upscaling_one_node(node, true_var);
	}
	return StaggerReturnCode::ResidualMinimized; // This is the value that should be returned if everything is ok
	}

	// A function that return the diffusion coefficient when porosity is equal to 'porosity'
	scalar_t coeff_diffusion_law(scalar_t porosity)
	{
	static constexpr scalar_t de_0 = 2.726e-12; //1e-11;
	static constexpr scalar_t factor = 1e4;
	static constexpr scalar_t porosity_r = 0.02;

	const scalar_t ratio = ((porosity - porosity_r)/(0.25 - porosity_r));
	const scalar_t d_e = factorde_0std::pow(ratio, 3.32);
	return std::min(d_e, 1e4*1e-10);
	}

	// Compute the upscaling for 'node'
	void upscaling_one_node(index_t node, SaturatedVariablesPtr true_var)
	{
	// AdimensionalSystemSolutionExtractor is the class to use to
	// extract information from a SpecMiCP solution
	// To obtain correct information the correct units must be used
	AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node),
	m_data, m_units_set);
	// We can obtain the porosity very easily :
	scalar_t porosity = extractor.porosity();
	true_var->vel_porosity(node) += (porosity - true_var->porosity(node))/m_dt;
	true_var->porosity(node) = porosity;
	true_var->diffusion_coefficient(node) = coeff_diffusion_law(porosity);
	}

	private:
	RawDatabasePtr m_data;
	units::UnitsSet m_units_set;
	index_t m_id_cshj;
	scalar_t m_initial_sat_cshj;
	scalar_t m_dt;
	};

	// ==============
	// Output
	// ==============

	// Many information is generated throughout a computation and we can't save everything
	// The reactmicp runner can call a function to output the desired information
	// Here the output function is a member function for convenience
	class OutputPolicy
	{
	public:
	// Initialise the output
	// SimulationInformation contain basic information from the simulation
	OutputPolicy(const SimulationInformation& simul_info,
	mesh::Mesh1DPtr the_mesh,
	RawDatabasePtr the_database,
	units::UnitsSet the_units):
	m_mesh(the_mesh),
	m_database(the_database),
	m_units(the_units)

	{
	// Open the files
	// Save the pH
	out_c_ph.open(simul_info.complete_filepath("ph", suffix));
	// Save the total aqueous concentration of HCO3[-]
	out_c_hco3.open(simul_info.complete_filepath("c_hco3", suffix));
	// save the total solid concentration of s_ca
	out_s_ca.open(simul_info.complete_filepath("s_ca", suffix));
	// Save the volume fraction of calcite
	out_calcite.open(simul_info.complete_filepath("calcite", suffix));
	// Save the porosity
	out_poro.open(simul_info.complete_filepath("poro", suffix));
	// save useful indices from the database
	m_id_hco3 = m_database->get_id_component("HCO3[-]");
	m_id_ca = m_database->get_id_component("Ca[2+]");
	m_id_calcite = m_database->get_id_mineral("Calcite");
	}


	void output(scalar_t time, VariablesBasePtr variables)
	{
	SaturatedVariablesPtr true_var = cast_var_from_base(variables);

	out_c_ph << time;
	out_c_hco3 << time;
	out_s_ca << time;
	out_calcite << time;
	out_poro << time;

	for (index_t node: m_mesh->range_nodes())
	{
	// Again AdimensionalSystemSolutionExtractor is used to obtain information about the speciation system
	AdimensionalSystemSolutionExtractor extractor(true_var ->equilibrium_solution(node), m_database, m_units);
	out_c_ph << "\t" << extractor.pH();
	out_c_hco3 << "\t" << true_var->aqueous_concentration(node, m_id_hco3);
	out_s_ca << "\t" << true_var->solid_concentration(node, m_id_ca);
	out_calcite << "\t" << extractor.volume_fraction_mineral(m_id_calcite);
	out_poro << "\t" << extractor.porosity();
	}
	out_c_hco3 << std::endl;
	out_c_ph << std::endl;
	out_s_ca << std::endl;
	out_calcite << std::endl;
	out_poro << std::endl;

	}

	private:
	mesh::Mesh1DPtr m_mesh;
	RawDatabasePtr m_database;
	units::UnitsSet m_units;

	std::string suffix { "dat" };

	std::ofstream out_c_ph;
	std::ofstream out_c_hco3;
	std::ofstream out_s_ca;
	std::ofstream out_calcite;
	std::ofstream out_poro;

	index_t m_id_hco3;
	index_t m_id_ca;
	index_t m_id_calcite;
	};

	// ===================
	// 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;

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

	// Obtain the database
	RawDatabasePtr the_database = get_database("../data/cemdata.js");

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

	// Set the options
	AdimensionalSystemSolverOptions specmicp_options = get_specmicp_options();
	units::UnitsSet units_set = specmicp_options.units_set;


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

	// Obtain the mesh
	mesh::Mesh1DPtr the_mesh = get_mesh(0.0050); // uniform
	//mesh::Mesh1DPtr the_mesh = get_mesh(); // non uniform
	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);

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

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


	// 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

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

	std::shared_ptr<SaturatedTransportStagger> transport_stagger =
	std::make_shared<satdiff::SaturatedTransportStagger>(variables, list_fixed_nodes);
	set_transport_options(transport_stagger->options_solver());

	UpscalingStaggerPtr upscaling_stagger =
	std::make_shared<CarboUspcaling>(the_database, units_set);

	// Solver
	// ------

	// This section initialize the reactive transport solver

	ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger);
	set_reactive_solver_options(solver.get_options());
	transport_stagger->initialize(variables);
	equilibrium_stagger->initialize(variables);
	upscaling_stagger->initialize(variables);


	// Some basic information about the simulation
	// First argument : the name of the simulation
	// Second argument : the output function will be called every 3600s
	SimulationInformation simul_info("carbo", 3600);
	simul_info.output_prefix = "out_carbo_";

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

	// Create the runner : loop through the timestep
	ReactiveTransportRunner runner(solver, 0.1, 100.0, simul_info);


	// Create the output class
	OutputPolicy output_policy(simul_info, the_mesh, the_database, specmicp_options.units_set);
	// and bind it to the runner
	runner.set_output_policy(std::bind(&OutputPolicy::output, &output_policy, std::placeholders::_1, std::placeholders::_2));

	// Solve the problem
	runner.run_until(30243600, cast_var_to_base(variables));

	// output information at the end of the timestep

	io::print_minerals_profile(the_database, variables, the_mesh, 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"));
	}


	AdimensionalSystemSolverOptions get_specmicp_options()
	{
	AdimensionalSystemSolverOptions options;


	options.solver_options.maxiter_maxstep = 100;
	options.solver_options.maxstep = 50.0;
	options.solver_options.threshold_cycling_linesearch = 1e-6;

	options.solver_options.set_tolerance(1e-8, 1e-16);
	options.solver_options.disable_condition_check();
	options.solver_options.disable_descent_direction();
	options.solver_options.enable_non_monotone_linesearch();
	options.solver_options.enable_scaling();


	options.system_options.non_ideality_tolerance = 1e-14;
	options.system_options.cutoff_total_concentration = 1e-10;
	options.system_options.scaling_electron = 1e4;
	options.units_set.length = units::LengthUnit::centimeter;

	return options;
	}

	void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options)
	{
	transport_options.maximum_iterations = 450;
	transport_options.quasi_newton = 3;
	transport_options.residuals_tolerance = 1e-8;
	transport_options.step_tolerance = 1e-14;
	transport_options.absolute_tolerance = 1e-18;

	transport_options.alpha = 1.0;
	//transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
	transport_options.sparse_solver = SparseSolver::SparseLU;
	//transport_options.sparse_solver = SparseSolver::GMRES;
	transport_options.threshold_stationary_point = 1e-4;
	}

	void set_reactive_solver_options(reactmicp::solver::ReactiveTransportOptions& reactive_options)
	{
	reactive_options.maximum_iterations = 50; //500;
	reactive_options.residuals_tolerance = 1e-2;
	reactive_options.good_enough_tolerance = 1.0;

	reactive_options.absolute_residuals_tolerance = 1e-14;
	reactive_options.implicit_upscaling = true;
	}

	RawDatabasePtr get_database(const std::string& path_to_database) {
	database::Database the_database(path_to_database);
	std::map<std::string, std::string> swapping({
	{"H[+]","HO[-]"},
	{"Si(OH)4", "SiO(OH)3[-]"},
	{"Al[3+]","Al(OH)4[-]"},
	{"Fe[2+]", "Fe(OH)4[-]"}
	});
	the_database.swap_components(swapping);
	the_database.remove_gas_phases();
	return the_database.get_database();

	}



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

	Vector oxyde_compo(5);
	oxyde_compo << 66.98, 21.66, 2.78, 2.96, 4.41;
	compo_oxyde_to_mole(oxyde_compo);
	scalar_t mult = laws::density_water(units::celsius(25.0), units_set.length, units_set.mass);
	//oxyde_compo *= 0.01098; //= 1e-2; //poro 0.4
	oxyde_compo *= 0.00967; //= 1e-2; //poro 0.47
	Formulation formulation;
	formulation.mass_solution = 1.0;
	formulation.amount_minerals = {
	{"CaO(ox)", oxyde_compo(0)},
	{"SiO2(ox)", oxyde_compo(1)},
	{"Al2O3(ox)", oxyde_compo(2)},
	{"SO3(ox)", oxyde_compo(3)},
	{"Fe2O3(ox)", oxyde_compo(4)},
	{"Calcite", oxyde_compo(0)*1e-4}
	};
	formulation.extra_components_to_keep = {"HCO3[-]"};
	formulation.concentration_aqueous =
	{
	{"Na[+]", mult*1.0298},
	{"K[+]", mult*0.0801},
	{"Cl[-]", mult*0.0001},
	{"HO[-]", mult*1.8298}
	},
	formulation.minerals_to_keep = {
	"Portlandite",
	"CSH,jennite",
	"CSH,tobermorite",
	"SiO2,am",
	"Calcite",
	"Al(OH)3,am",
	"C3AH6",
	"C4AH13",
	"C2AH8",
	"CAH10",
	"Monosulfoaluminate",
	"Tricarboaluminate",
	"Monocarboaluminate",
	"Hemicarboaluminate",
	//"Straetlingite",
	"Gypsum",
	"Ettringite",
	//"Thaumasite",
	"C3FH6",
	"C4FH13",
	"C2FH8",
	"Fe(OH)3(mic)",
	"Fe-Ettringite",
	"Fe-Monosulfate",
	"Fe-Monocarbonate",
	"Fe-Hemicarbonate",
	//"Fe-Straetlingite",
	};
	Dissolver dissolver(the_database);
	dissolver.set_units(units_set);

	AdimensionalSystemConstraints constraints;
	constraints.set_charge_keeper(1);

	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("SiO(OH)3[-]")) = 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_oxyde_to_mole(Vector& compo_oxyde)
	{
	constexpr double M_CaO = 56.08;
	constexpr double M_SiO2 = 60.09;
	constexpr double M_Al2O3 = 101.96;
	constexpr double M_SO3 = 80.06;
	constexpr double M_FE2O3 = 159.7;

	Vector molar_mass(5);
	molar_mass << 1.0/M_CaO, 1.0/M_SiO2, 1.0/M_Al2O3, 1.0/M_SO3, 1.0/M_FE2O3;
	compo_oxyde = compo_oxyde.cwiseProduct(molar_mass);

	}


	mesh::Mesh1DPtr get_mesh(scalar_t dx)
	{
	const scalar_t total_length = 0.75/2.0 + dx; //cm
	const scalar_t height = 20.0;

	index_t nb_nodes = total_length/dx + 1;

	return mesh::uniform_axisymmetric_mesh1d(nb_nodes, total_length, dx, height);
	}

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

	const scalar_t height = 20;

	// radius <<
	// 3.751,
	// 3.750,
	// 3.749,
	// 3.748,
	// 3.747,
	// 3.745,
	// 3.740,
	// 3.735,
	// 3.730,
	// 3.720,
	// 3.710,
	// 3.700,
	// 3.675,
	// 3.650,
	// 3.625,
	// 3.600,
	// 3.575,
	// 3.550,
	// 3.525,
	// 3.500,
	// 3.475,
	// 3.450,
	// 3.425,
	// 3.400,
	// 3.375,
	// 3.350,
	// 3.325,
	// 3.300,
	// 3.250,
	// 3.200,
	// 3.150,
	// 3.100,
	// 3.050,
	// 3.000,
	// 2.975,
	// 2.950,
	// 2.925,
	// 2.900,
	// 2.875,
	// 2.850,
	// 2.825,
	// 2.800,
	// 2.775,
	// 2.750,
	// 2.700,
	// 2.650;

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

carbonationfe.cppNo OneTemporaryActions

File Metadata

carbonationfe.cppView Options

Event Timeline

carbonationfe.cpp
No OneTemporary
Actions

carbonationfe.cpp
View Options