Diffusion SpecMiCP / ReactMiCP (dev) Browse

carbonationfe.cpp
b174f1ddf352dev

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

specmicp/examples/reactmicp/saturated_react/carbonationfe.cppb174f1ddf352dev

carbonationfe.cppView Options

specmicp/examples/reactmicp/saturated_react/carbonationfe.cpp
b174f1ddf352dev

carbonationfe.cpp
View Options