carbonation.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Mon, Jul 7, 11:34

carbonation.cpp
View Options

	#include <Eigen/Core>
	#include <Eigen/QR>

	#include "specmicp/adimensional/adimensional_system_solver.hpp"
	#include "specmicp/problem_solver/formulation.hpp"
	#include "specmicp/problem_solver/dissolver.hpp"
	#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"

	#include "reactmicp/systems/saturated_react/variables.hpp"
	#include "reactmicp/systems/saturated_react/equilibrium_stagger.hpp"
	#include "reactmicp/systems/saturated_react/transport_stagger.hpp"
	#include "reactmicp/systems/saturated_react/init_variables.hpp"

	#include "reactmicp/solver/reactive_transport_solver.hpp"
	#include "reactmicp/solver/reactive_transport_solver_structs.hpp"
	#include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp"
	#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
	#include "reactmicp/solver/timestepper.hpp"


	#include "physics/laws.hpp"

	#include "dfpm/mesh.hpp"

	#include "utils/log.hpp"
	#include "utils/timer.hpp"

	#include "utils/io/reactive_transport.hpp"

	#include <iostream>
	#include <fstream>

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

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

	AdimensionalSystemSolution get_cement_initial_compo(
	RawDatabasePtr the_database,
	const AdimensionalSystemConstraints& constraints,
	const AdimensionalSystemSolverOptions& options
	);

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

	AdimensionalSystemSolution get_bc_initial_compo(
	RawDatabasePtr the_database,
	const AdimensionalSystemConstraints& constraints,
	const AdimensionalSystemSolverOptions& options
	);


	AdimensionalSystemConstraints get_bc_initial_constraints(
	RawDatabasePtr the_database,
	const units::UnitsSet& units_set);

	RawDatabasePtr get_database(const std::string& path_to_database);

	Vector initial_compo(const Vector& publi);

	AdimensionalSystemSolverOptions get_specmicp_options();

	mesh::Mesh1DPtr get_mesh(scalar_t dx);

	void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options);

	void set_reactive_solver_options(reactmicp::solver::ReactiveTransportOptions& reactive_options);

	void print_minerals_profile(
	RawDatabasePtr the_database,
	SaturatedVariablesPtr variables,
	mesh::Mesh1DPtr the_mesh,
	const std::string& filepath
	);

	class CarboUspcaling: public solver::UpscalingStaggerBase
	{
	public:

	static constexpr scalar_t de_0 = 1e-11;
	static constexpr scalar_t factor = 1e4;
	static constexpr scalar_t porosity_r = 0.02;

	CarboUspcaling(RawDatabasePtr the_database,
	units::UnitsSet the_units):
	m_data(the_database),
	m_units_set(the_units),
	m_dt(HUGE_VAL)
	{

	}

	//! \brief Initialize the stagger at the beginning of the computation
	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);
	}


	}

	//! \brief Initialize the stagger at the beginning of an iteration
	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;
	}
	}

	//! \brief Solve the equation for the 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;
	}

	scalar_t coeff_diffusion_law(scalar_t porosity)
	{
	return factorde_0((porosity - porosity_r)/(0.4 - porosity_r));
	}

	void upscaling_one_node(index_t node, SaturatedVariablesPtr true_var)
	{
	AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node),
	m_data, m_units_set);

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

	int main()
	{

	specmicp::logger::ErrFile::stream() = &std::cerr;
	specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;

	RawDatabasePtr the_database = get_database("../data/cemdata_specmicp.js");

	for (auto label: the_database->labels_basis)
	std::cout << label << std::endl;

	AdimensionalSystemSolverOptions specmicp_options = get_specmicp_options();
	units::UnitsSet units_set = specmicp_options.units_set;

	AdimensionalSystemConstraints cement_constraints = get_cement_initial_constraints(the_database, units_set);

	for (auto label: the_database->labels_minerals)
	std::cout << label << std::endl;

	AdimensionalSystemSolution cement_solution = get_cement_initial_compo(
	the_database,
	cement_constraints,
	specmicp_options);

	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;

	AdimensionalSystemConstraints bc_constraints = get_bc_initial_constraints(
	the_database,
	specmicp_options.units_set
	);

	AdimensionalSystemSolution bc_solution = get_bc_initial_compo(
	the_database,
	bc_constraints,
	specmicp_options
	);

	AdimensionalSystemSolutionExtractor extractor2(bc_solution, the_database, units_set);

	std::cout << "pH : " << extractor2.pH() << std::endl;

	mesh::Mesh1DPtr the_mesh = get_mesh(0.0075);
	index_t nb_nodes = the_mesh->nb_nodes();
	std::vector<specmicp::AdimensionalSystemConstraints> list_constraints = {cement_constraints,
	bc_constraints};
	std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition = {cement_solution, bc_solution};
	std::vector<int> index_initial_state(nb_nodes, 0);
	std::vector<int> index_constraints(nb_nodes, 0);
	index_initial_state[0] = 1;
	index_constraints[0] = 1 ;

	std::vector<index_t> list_fixed_nodes = {0};


	satdiff::SaturatedVariablesPtr variables =
	satdiff::init_variables(the_mesh, the_database, units_set,
	list_fixed_nodes,
	list_initial_composition,
	index_initial_state);

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

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

	std::cout << variables->upscaling_variables() << std::endl;

	set_reactive_solver_options(solver.get_options());

	scalar_t dt = 1.0;
	index_t cnt = 0;

	Timestepper timestepper(1e-2, 10.0, 30243600, 2);

	database::Database dbhandler(the_database);
	index_t id_oh = dbhandler.component_label_to_id("HO[-]");
	index_t id_hco3 = dbhandler.component_label_to_id("HCO3[-]");
	index_t id_ca = dbhandler.component_label_to_id("Ca[2+]");

	std::string prefix = "out_carbo_";
	std::string suffix = ".dat";

	std::ofstream out_c_oh(prefix+"ph"+suffix);
	std::ofstream out_c_hco3(prefix+"c_hco3"+suffix);
	std::ofstream out_s_ca(prefix+"s_ca"+suffix);
	std::ofstream out_s_co2(prefix+"s_co2"+suffix);
	std::ofstream out_iter(prefix+"iter"+suffix);

	io::print_reactmicp_performance_long_header(&out_iter);

	// std::vector<std::ofstream> out_minerals(the_database->nb_mineral);
	// for (index_t mineral: the_database->range_mineral())
	// {
	// const std::string path = prefix+"min_"+the_database->labels_minerals[mineral]+suffix;
	// out_minerals[mineral].open(path);
	// }

	while (timestepper.get_total() < timestepper.get_total_target())
	{
	Timer step_timer;
	step_timer.start();
	reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables);
	step_timer.stop();

	io::print_reactmicp_performance_long(&out_iter, cnt, timestepper.get_total()+dt, solver.get_perfs());

	dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
	if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
	{
	dt = 0.01;
	variables->reset_main_variables();
	retcode = solver.solve_timestep(dt, variables);
	dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
	io::print_reactmicp_performance_long(&out_iter, cnt, timestepper.get_total(), solver.get_perfs());
	}



	if (cnt % 50 == 0)
	{


	out_c_oh << timestepper.get_total();
	out_c_hco3 << timestepper.get_total();
	out_s_ca << timestepper.get_total();
	out_s_co2 << timestepper.get_total();

	//for (index_t mineral: the_database->range_mineral())
	// out_minerals[mineral] << timestepper.get_total();

	for (index_t node: the_mesh->range_nodes())
	{
	AdimensionalSystemSolutionExtractor extractor(variables->equilibrium_solution(node), the_database, units_set);
	out_c_oh << "\t" << extractor.pH();
	out_c_hco3 << "\t" << variables->aqueous_concentration(node, id_hco3, variables->displacement());
	out_s_ca << "\t" << variables->solid_concentration(node, id_ca, variables->displacement());
	out_s_co2 << "\t" << variables->solid_concentration(node, id_hco3, variables->displacement());

	//for (index_t mineral: the_database->range_mineral())
	// out_minerals[mineral] << "\t" << extractor.volume_fraction_mineral(mineral);
	}
	out_c_hco3 << std::endl;
	out_c_oh << std::endl;
	out_s_ca << std::endl;
	out_s_co2 << std::endl;
	//for (index_t mineral: the_database->range_mineral())
	// out_minerals[mineral] << std::endl;

	}

	++cnt;
	}

	print_minerals_profile(the_database, variables, the_mesh, prefix+"profiles_end"+suffix);

	io::print_reactmicp_timer(&out_iter, solver.get_timer());
	}

	AdimensionalSystemSolverOptions get_specmicp_options()
	{
	AdimensionalSystemSolverOptions options;

	options.solver_options.steptol = 1e-10;
	options.solver_options.maxiter_maxstep = 100;
	options.solver_options.maxstep = 10.0;
	options.solver_options.fvectol = 1e-10;
	options.solver_options.steptol = 1e-16;
	options.solver_options.condition_limit = -1;
	options.solver_options.factor_descent_condition = -1;
	options.solver_options.threshold_cycling_linesearch = 1e-6;
	options.solver_options.non_monotone_linesearch = true;
	options.solver_options.use_scaling = true;

	options.system_options.non_ideality_tolerance = 1e-14;

	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-6;
	transport_options.step_tolerance = 1e-16;
	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-4;
	reactive_options.good_enough_tolerance = 1.0;

	reactive_options.absolute_residuals_tolerance = 1e-6;
	reactive_options.implicit_upscaling = false;
	}

	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[-]"}
	});
	the_database.swap_components(swapping);

	return the_database.get_database();

	}

	void print_minerals_profile(
	RawDatabasePtr the_database,
	SaturatedVariablesPtr variables,
	mesh::Mesh1DPtr the_mesh,
	const std::string& filepath
	)
	{
	std::ofstream ofile(filepath);

	ofile << "Radius";
	for (index_t mineral: the_database->range_mineral())
	{
	ofile << "\t" << the_database->labels_minerals[mineral];
	}
	ofile << "\t" << "Porosity";
	ofile << "\t" << "pH";
	ofile << std::endl;

	for (index_t node: the_mesh->range_nodes())
	{
	ofile << the_mesh->get_position(node);
	AdimensionalSystemSolutionExtractor extractor(variables->equilibrium_solution(node), the_database, units::UnitsSet());

	for (index_t mineral: the_database->range_mineral())
	{
	ofile << "\t" << extractor.volume_fraction_mineral(mineral);
	}
	ofile << "\t" << variables->porosity(node);
	ofile << "\t" << extractor.pH();
	ofile << std::endl;

	}

	ofile.close();
	}

	AdimensionalSystemConstraints get_cement_initial_constraints(
	RawDatabasePtr the_database,
	const units::UnitsSet& units_set)
	{
	Vector init_min_publi(5);
	init_min_publi << 12.9, 1.62, 0.605, 0.131, 0.001;
	Vector init_min;
	init_min = initial_compo(init_min_publi);
	std::cout << "Init min \n ---- \n" << init_min << "\n ---- \n" << std::endl;

	scalar_t mult = 1e-6;
	//scalar_t mult = 6.73153e-4;
	//scalar_t mult = 0.001296315;
	init_min *= 6.71146e-4;

	Formulation formulation;
	formulation.mass_solution = 1.0;
	formulation.amount_minerals = {
	{"Portlandite", init_min(0)},
	{"CSH,jennite", init_min(1)},
	{"Monosulfoaluminate", init_min(2)},
	{"Ettringite", init_min(3)},
	{"Calcite", init_min(4)}
	};
	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",
	"Monosulfoaluminate",
	"Tricarboaluminate",
	"Monocarboaluminate",
	"Hemicarboaluminate",
	"Straetlingite",
	"Gypsum",
	"Ettringite",
	"Thaumasite"
	};
	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;

	database::Database dbhandler(the_database);

	Vector total_concentrations(the_database->nb_component);
	total_concentrations.setZero();
	total_concentrations(dbhandler.component_label_to_id("K[+]")) = kcl;
	total_concentrations(dbhandler.component_label_to_id("Cl[-]")) = nacl + kcl;
	total_concentrations(dbhandler.component_label_to_id("Na[+]")) = nacl;

	AdimensionalSystemConstraints constraints;
	constraints.add_fixed_activity_component(dbhandler.component_label_to_id("HO[-]"), -10.3);
	constraints.set_charge_keeper(dbhandler.component_label_to_id("HCO3[-]"));
	constraints.total_concentrations = total_concentrations;

	constraints.set_saturated_system();

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

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

	Vector variables;
	solver.initialise_variables(
	variables,
	0.4,
	{
	{"HO[-]", -1.2},
	{"Al(OH)4[-]", -4},
	{"HCO3[-]", -5},
	{"Ca[2+]", -1.6}
	});
	micpsolver::MiCPPerformance perf = solver.solve(variables);

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


	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(scalar_t dx)
	{
	const scalar_t total_length = 7.5/2.0; //cm
	const scalar_t height = 20;

	index_t nb_nodes = total_length/dx + 1;

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

carbonation.cppNo OneTemporaryActions

File Metadata

carbonation.cppView Options

Event Timeline

carbonation.cpp
No OneTemporary
Actions

carbonation.cpp
View Options