Diffusion SpecMiCP / ReactMiCP (configuration) Browse

api.md
87ddbd8df2daconfiguration

api.md
View Options

			SpecMiCP API {#specmicp_api_file}
			============

			[TOC]

			The SpecMiCP API is divided into three main components:

			- [The constraints](#specmicp_constraints)
			- [The solver](#specmicp_solver)
			- [The solution](#specmicp_solution)

			# The constraints {#specmicp_constraints}

			## Introduction {#constraints_intro}

			The system considered is a fixed volume with an aqueous solution and possibly a
			gaseous phase and zero, one or several solid phase(s) at equilibrium.

			The [constraints](@ref specmicp::AdimensionalSystemConstraints) are the set of
			constraints on the chemical system. They define the amount of components in the
			system and the type of system (closed, open). In a closed system the amount of
			matter is conserved. In an open system matter may be exchanged to conserve the
			activity or the fugacity of a species.

			The (open, closed) system is define for each components (see [the database
			documentation](#database)). A system may be allowed to exchanged some
			components (for example through the gaseous phases) while the amounts of other
			components are fixed (the components that are only present in the aqueous
			solution or the solid phases).

			## The Constraint class

			The constraints are stored in the `struct` [AdimensionalSystemConstraints](@ref
			specmicp::AdimensionalSystemConstraints).

			It can be build in two way

			~~~~~~~~{.cpp}
			#include <specmicp/specmicp/adimensional_system_structs.hpp>
			specmicp::AdimensionalSystemConstraints the_constraints()
			~~~~~~~~

			or

			~~~~~~~~{.cpp}
			#include <specmicp/specmicp/adimensional_system_structs.hpp>
			specmicp::Vector total_concentrations
			// set the total_concentrations vector
			specmicp::AdimensionalSystemConstraints the_constraints(total_concentrations)
			~~~~~~~~

			The next section describes how to set the vector of total concentrations.

			## Total concentrations {#constraints_totconc}

			The total concentration of the "closed-system" components are fixed and it is
			the constraint that should be set.

			The total concentrations are contained in a [vector](@ref specmicp::Vector).
			Its size is equal to the number of components in the system

			~~~~~~~{.cpp}
			specmicp::RawDatabasePtr raw_data = // obtain the database in some way
			specmicp::Vector total_concentrations(raw_data->nb_components);
			// set the values in the total_concentrations vector
			specmicp::AdimensionalSystemConstraints the_constraints(total_concentrations);
			~~~~~~~

			or

			~~~~~~~{.cpp}
			specmicp::RawDatabasePtr raw_data = // obtain the database in some way
			specmicp::Vector total_concentrations(raw_data->nb_components);
			// set the values in the total_concentrations vector
			specmicp::AdimensionalSystemConstraints the_constraints;
			the_constraints.set_total_concentrations(total_concentrations);
			~~~~~~~

			The total concentrations can be computed by hand or obtained using the
			[dissolver](specmicp::Dissolver) as explained in the next section.

			## The dissolver {#constraints_dissolver}

			In a real experiments, the total concentrations are not the relevant variables.
			Instead, the initial composition is given. The classes specmicp::Dissolver and
			specmicp::Formulation allow to compute the total concentrations and prepare the
			database for the computation given an initial composition.

			### Formulation

			The [formulation](@ref specmicp::Formulation) allow the modeler to specify an
			initial composition.

			\warning this class do NOT take handle units. The user need to ensure that the
			input units are consistent.

			The first step is to set the mass of solution. This is the mass of the aqueous
			solution per volume units. Then the user can add aqueous species or solid
			phases.

			~~~~~~{.cpp}
			specmicp::Formulation the_formulation;
			// set the mass of the solution
			the_formulation.set_mass_solution(0.5); # kg/L
			// add aqueous species
			the_formulation.add_aqueous_species("Na[+]", 0.5); // Molal (mol/kgw)
			the_formulation.add_aqueous_species("Cl[-]", 0.5); // Molal
			the_formulation.add_aqueous_species("CO2", 0.001); // Molal
			// add a solid phase
			the_formulation.add_mineral("Portlandite", 0.5); // mol/L
			~~~~~~

			The aqueous species can be components (except water and electron) or aqueous
			species.
			The solid phases can be solid phases at equilibrium or kinetics solid phases in
			the database.

			### Dissolver

			The [dissolver](@ref specmicp::Dissolver) dissolves the [formulation](@ref
			specmicp::Formulation) previously defined.
			When computing the total concentrations, this class also removes the unecessary
			components from the database.


			\warning If a component is not included in the initial composition it
			will removed from the database. If a component must be kept in the database
			for a following computation, the user need to register it in the
			[formulation](@ref specmicp::Formulation)
			~~~~~~{.cpp}
			specmicp::Formulation the_formulation;
			the_formulation.keep_component("Ca[2+]");
			~~~~~~

			The dissolver can be used as follow :

			~~~~~~{.cpp}
			specmicp::Formulation the_formulation;
			... // set the formulation
			specmicp::Dissolver dissolver(raw_data);
			specmicp::AdimensionalSystemConstraints the_constraints(
			dissolver.dissolver(formulation));
			~~~~~~

			## Electroneutrality {#constraints_chargekeeper}

			A mass conservation equation for an aqueous species may be replaced by the
			electroneutrality, or charge balance, condition.
			It is optional if all components are conserved but mandatory if the equation of
			a charged component is either [fixed fugacity](#constraints_fixedfugacity) or
			[fixed activity](#constraints_fixedactivity).

			The following snippet declares that the quantity of "HO[-]" is adjusted to
			satisfy the mass balance :

			~~~~~{.cpp}
			specmicp::RawDatabasePtr raw_data; // obtained in some way
			specmicp::AdimensionalSystemConstraints the_constraints;
			the_constraints.set_charge_keeper(
			raw_data->get_id_component("HO[-]")
			);
			~~~~~

			In practice, it was observed that the electroneutrality condition converges
			slightly better the the mass conservation equation.

			## Fixed fugacity {#constraints_fixedfugacity}

			The fugacity of a gas can be fixed.
			In this case the conservation equation for a component will be replaced by the
			fixed-fugacity constraint.
			The corresponding component need to be chosen by the modeler but should be
			fairly obvious. For exemple if the fugacity of `CO2(g)` is fixed the
			corresponding component is `HCO3[-]` or whichever carbonate species is in the
			basis.

			To set a fixed fugacity gas, the [add_fixed_fugacity_gas()](@ref
			specmicp::AdimensionalSystemConstraints::add_fixed_fugacity_gas() method must
			be called. It is used as follow :

			~~~~~{.cpp}
			specmicp::RawDatabasePtr raw_data; // obtained in some way
			specmicp::AdimensionalSystemConstraints the_constraints;
			the_constraints.add_fixed_fugacity_gas(
			raw_data->get_id_gas("CO2(g)"),
			raw_data->get_id_component("HCO3[-]"),
			-2.0
			);
			the_constraints.set_charge_keeper(
			raw_data->get_id_component("HO[-]")
			); // HCO3[-] is charged, so we need to conserve Electroneutrality
			~~~~~

			The last argument of [add_fixed_fugacity_gas()](@ref
			specmicp::AdimensionalSystemConstraints::add_fixed_fugacity_gas() is the
			\f$\log_{10}\f$ of the fugacity. In this case, it corresponds to a fugacity of
			0.01.

			If the corresponding component is charged, then the charged balanced must be
			enforced explicitely, see [Electroneutrality](#constraints_chargekeeper).

			## Fixed activity {#constraints_fixedactivity}

			A component can also have a fixed activity. The main usage of this constraints
			is to fix the pH of the system.

			The following setup set the pH to 3.5 by adding HCl to the system :

			~~~~~{.cpp}
			specmicp::RawDatabasePtr raw_data; // obtained in some way
			specmicp::AdimensionalSystemConstraints the_constraints;
			the_constraints.add_fixed_activity_component(
			raw_data->get_id_gas("H[+]"),
			-3.5
			);
			the_constraints.set_charge_keeper(
			raw_data->get_id_component("Cl[-]")
			);
			~~~~~

			If the component is charged, then the charged balanced must be
			enforced explicitely, see [Electroneutrality](#constraints_chargekeeper).

			## Fixed molality {#constraints_fixedmolality}

			A component may also have its molality fixed.
			It can be used, for example, to simulate the effect of curing on a cement sample.

			The following code set the molality of NaCl to 0.5M :

			~~~~~~{.cpp}
			specmicp::RawDatabasePtr raw_data; // obtained in some way
			specmicp::AdimensionalSystemConstraints the_constraints;
			the_constraints.add_fixed_molality_component(
			raw_data->get_id_gas("Na[+]"), std::log10(0.5));
			the_constraints.add_fixed_molality_component(
			raw_data->get_id_gas("Cl[-]"), std::log10(0.5));
			~~~~~~

			## Water conservation {#constraints_water}

			Three constraints are available for the water (liquid water, `H2O`) :

			- The amount of water is conserved (default) :
			[enable_conservation_water()](@ref
			specmicp::AdimensionalSystemConstraints::enable_conservation_water)
			- The sample is saturated : [set_saturated_system()](@ref
			specmicp::AdimensionalSystemConstraints::set_saturated_system)
			-The amount of water is not conserved :
			[disable_conservation_water](@ref
			specmicp::AdimensionalSystemConstraints::disable_conservation_water)

			It was observed that the conservation of water helps the convergence, even for
			system where the water is in large excess.
			Therefore, the last option is not recommended.

			Example : set a saturated system


			~~~~~{.cpp}
			specmicp::AdimensionalSystemConstraints the_constraints;
			the_constraints.set_saturated_system();
			~~~~~

			## Sorption model {#constraints_sorption}

			The sorption model is disabled by default. It can be activated with the
			[enable_surface_model](@ref
			specmicp::AdimensionalSystemConstraints::enable_surface_model) method. It takes
			one argument, the total concentration of sorption sites.

			~~~~~{.cpp}
			specmicp::AdimensionalSystemConstraints the_constraints;
			the_constraints.enable_surface_model(1.2); // mol/V
			~~~~~

			## Inert volume {#constraints_inert}

			The user may add a filler phase in the system. This phase does not take part in
			any reaction but is used to compute the volume fractions.

			~~~~~{.cpp}
			specmicp::AdimensionalSystemConstraints the_constraints;
			the_constraints.set_inert_volume_fraction(0.2);
			~~~~~

			# The solver {#specmicp_solver}

			## Introduction {#solver_intro}

			The solver is the core of the algorithm. It sets up the program according the
			[constraints](#constraints) given by the user, it sets up the [MiCP
			solver](@ref specmicp::micpsolver::MiCPSolver), it initializes the variables,
			run the MiCP solver and it does the post-processing before returning the answer
			to the user.

			Its main purpose is to hide the complex and delicate steps involved and to
			propose a simpler interface to the user.

			## Initializing a solver {#solver_init}

			The solver is initialized with two pieces of information :

			- The [database](@ref database)
			- A set of [constraints](#specmicp_constraints)

			In addition the user may may provide a set of [options](#solver_options).

			~~~~~~{.cpp}
			specmicp::RawDatabasePtr the_database = ... ;
			specmicp::AdimensionalSystemConstraints the_constraints = ... ;
			// Initialize the solver
			specmicp::AdimensionalSystemSolver the_solver(the_database, the_constraints);
			~~~~~~

			If no previous solution is available, the user can provide information about how
			to initialize the variables using one of the
			[initialize_variables](@ref
			specmicp::AdimensionalSystemSolver::initialise_variables) methods.

			~~~~~~{.cpp}
			specmicp::AdimensionalSystemSolver the_solver = ....;
			specmicp::Vector variables;

			the_solver.initialise_variables(variables, 0.8, -3.0);
			~~~~~~

			If a similar system was solved previously, the previous solution can also be
			provided to the solver to initialize the solution with a better starting guess.

			~~~~~~{.cpp}
			specmicp::RawDatabasePtr the_database = ... ;
			specmicp::AdimensionalSystemConstraints the_constraints = ... ;
			specmicp::AdimensionalSystemSolution previous_solution = ...;
			// Initialize the solver
			specmicp::AdimensionalSystemSolver the_solver(
			the_database, the_constraints, previous_solution);
			specmicp::Vector variables = previous_solution.main_variables;
			~~~~~~

			## Solving the problem {#solver_solve}


			### Solving the equations

			The problem can be solved by using the [solve](@ref
			specmicp::AdimensionalSystemSolver::solve) method.
			This method takes the vector of main variables as first argument. the second
			argument is optional and is a boolean.
			If true, the solver will initialize the solution.
			If a previous solution or a custom initialisation is provided, this flag should
			be set to false (default).
			However, if no previous solution is available or the user don't know a good
			starting guess, this flag should be set to true.

			With a custom initialisation :

			~~~~~~{.cpp}
			#include <specmicp/specmicp/adimensional/adimensional_system_solver.hpp>
			specmicp::RawDatabasePtr the_database = ... ;
			specmicp::AdimensionalSystemConstraints the_constraints = ... ;
			// Initialize the solver
			specmicp::AdimensionalSystemSolver the_solver(
			the_database, the_constraints,);
			specmicp::Vector variables;
			the_solver.initialise_variables(0.5, {"Ca[2+]": -3.0, "HO[-]": -2.0});
			// solve the problem
			specmicp::micpsolver::MiCPPerformance perf = the_solver.solve(variables);
			~~~~~~

			Without a custom initialization :
			~~~~~~{.cpp}
			specmicp::RawDatabasePtr the_database = ... ;
			specmicp::AdimensionalSystemConstraints the_constraints = ... ;
			// Initialize the solver
			specmicp::AdimensionalSystemSolver the_solver(
			the_database, the_constraints,);
			specmicp::Vector variables;
			// solve the problem
			specmicp::micpsolver::MiCPPerformance perf = the_solver.solve(variables, true);
			~~~~~~

			### Checking the status of the solver

			The solve function [solve](@ref specmicp::AdimensionalSystemSolver::solve)
			returns a [struct](@ref specmicp::micpsolver::MiCPPerformance) containing
			information about the performance of the MiCP solver. It contains a lot of
			information such as the final residuals, the number of iterations but more
			importantly the return code of the solver.
			To check that the problem was succesfully solved, the following code can be used
			:

			~~~~~~{.cpp}
			// solve the problem
			specmicp::micpsolver::MiCPPerformance perf = the_solver.solve(variables, true);
			if (perf.return_code < specmicp::micpsolver::MiCPSolverReturnCode::Success)
			{
			throw std::runtime_error("Error : the system couldn't be solved !");
			}
			else
			{
			std::cout << "Hip hip hip, hourra, we have a solution !" << std::endl;
			}
			~~~~~~

			If the solver is unsuccessful at first, it will modify the initial conditions
			and try to solve the problem again.
			More information about the failure can be obtained by activating the
			[logger](@ref specmicp::logger)

			### Obtaining the solution

			The vector of main variables is only a partial view of the solution as it does
			not contain any information about the secondary species.
			The complete solution is stored in a `struct` [AdimensionalSystemSolution](@ref
			specmicp::AdimensionalSystemSolution).
			The solution can be retrieved from the solver by using the
			[get_raw_solution()](@ref specmicp::AdimensionalSystemSolver) method.
			The user should avoid querying information directly from this object but should
			used the [AdimensionalSystemSolutionExtractor](#specmicp_solution) class
			instead.


			~~~~~~{.cpp}
			specmicp::RawDatabasePtr the_database = ... ;
			specmicp::AdimensionalSystemConstraints the_constraints = ... ;
			// Initialize the solver
			specmicp::AdimensionalSystemSolver the_solver(
			the_database, the_constraints,);
			specmicp::Vector variables;
			the_solver.initialise_variables(variables, ...)
			// solve the problem
			specmicp::micpsolver::MiCPPerformance perf = the_solver.solve(variables);
			if (perf.return_code < specmicp::micpsolver::MiCPSolverReturnCode::Success)
			{
			throw std::runtime_error("Error : the system couldn't be solved !");
			}
			else
			{
			std::cout << "Hip hip hip, hourra, we have a solution !" << std::endl;
			}
			// Obtain the solution
			specmicp::AdimensionalSystemSolution solution = the_solver.get_raw_solution(variables);
			~~~~~~

			## Options {#solver_options}

			The [solver](@ref specmicp::AdimensionalSystemSolver) provides many options to customize its behavior.
			A complete description is available [here](@ref specmicp::AdimensionalSystemSolverOptions).
			This `struct` contains three main components :

			- the set of units : [units_set](@ref specmicp::AdimensionalSystemSolverOptions::units_set)
			- the options for the MiCPSolver : [solver_options](@ref specmicp::micpsolver::MiCPSolverOptions)
			- the options for the System : [system_options](@ref specmicp::AdimensionalSystemOptions)

			The following example shows the main options :

			~~~~~{.cpp}
			AdimensionalSystemSolverOptions options;

			// The first argument is the residuals non_ideality_tolerance
			// The second argument is the step tolerance
			options.solver_options.set_tolerance(1e-8, 1e-14);
			// Enable the non monotone linesearch (default)
			options.solver_options.enable_non_monotone_linesearch();
			// Enable the scaling (default)
			options.solver_options.enable_scaling();
			// This options set the maximum length of a Newton step
			// A small value will results in smaller steps
			options.solver_options.set_maximum_step_length(50.0);

			// Set the tolerance for the non-ideality model
			options.system_options.non_ideality_tolerance = 1e-12;
			// Set the threshold on the component total concentration
			// If the total concentration for a component is less than
			// this threshold the speciation of this component will not be computed
			options.system_options.cutoff_total_concentration = 1e-14;

			// Set the length unit
			options.units_set.length = units::LengthUnit::centimeter;
			~~~~~

			When a methods starts with `enable_` another method starting
			with `disable_` also exists and has the opposite effect.

			## Troubleshooting {#solver_troubleshooting}

			If the solver cannot solve a problem the first thing to do is to verify the [constraints](#specmicp_constraints).
			A non-physical problem is a lot harder to solve that a well-posed problem.
			Check database, the units, the total concentrations, ...

			If you are sure that the problem is well-posed, you may try to initialize the solver with a better [starting guess](#solver_init).

			If you are still experiencing problem there is a few options that you can adjust.

			Be sure that the non-monotone linesearch and the scaling of the jacobian are enabled :

			~~~~~{.cpp}
			specmicp::AdimensionalSystemSolver the_solver(...);
			the_solver.get_options().solver_options.enable_non_monotone_linesearch();
			the_solver.get_options().solver_options.enable_scaling();
			~~~~~

			If the logger often reports this warning : `Cycling has been detected by the linesearch - Taking the full Newton step`,
			you may want to disable the non monotone linesearch or using a lower threshold for the detection of cycling :

			~~~~~{.cpp}
			specmicp::AdimensionalSystemSolver the_solver(...);
			the_solver.get_options().solver_options.threshold_cycling_linesearch = 1e-7.
			~~~~~

			This warning is raised when a cycling is detected in the non-monotone linesearch.
			However it is sometimes a false-positive, and lowering the threshold will lead to convergence.


			If the number of iteration is too high, you may try to increase the maximum Newton step length.
			This step length is used to avoid taking steps that are too big \cite Dennis1983 and could lead to divergence.
			However, the default value is quite conservative and can be increased in many case. A value between 10 and 50 is recommended.

			~~~~~{.cpp}
			specmicp::AdimensionalSystemSolver the_solver(...);
			the_solver.get_options().solver_options.set_maximum_step_length(25, 50);
			~~~~~

			The second argument is the maximum number of iterations in which the solver can restrict the step length.


			# Extracting information from the solution {#specmicp_solution}

			A [solution](#ref specmicp::AdimensionalSystemSolution) is just a raw container.
			To extract useful pieces of information from it, the [SolutionExtractor](@ref specmicp::AdimensionalSystemSolutionExtractor) should be used.
			This `class` define a set of methods that allow to retrieve the physical variables in a friendly manner.

			The followin snippet initializes an extractor :

			~~~~~{.cpp}
			specmicp::RawDatabasePtr the_database = ... ;
			specmicp::units::UnitsSet the_units = ...;
			specmicp::AdimensionalSystemSolution the_solution = ;

			specmicp::AdimensionalSystemSolutionExtractor extractor(the_solution, the_database, the_units);
			~~~~~

			To avoid unecessary copy, the [SolutionExtractor](@ref specmicp::AdimensionalSystemSolutionExtractor) class only stores
			a reference to the solution. The database and the units used to initialized the extractor should be the ones used during the computation.

			The next examples presents some of the methods available.

			~~~~~{.cpp}
			specmicp::AdimensionalSystemSolutionExtractor extractor = ...;
			std::cout << "The volume fraction of water is : " << extractor.volume_fraction_water() << std::endl;
			std::cout << "The pH is : " << extractor.pH() << std::endl;
			std::cout << "The molality of Ca[2+] is : "
			<< extractor.molality_component(raw_data->get_id_component("Ca[2]")) <<std::endl;
			std::cout << "The volume fraction of portlandite is : "
			<< extractor.volume_fraction_mineral(raw_data->get_id_mineral("Portlandite")) << std::endl;
			~~~~~

			As a convention, a method to retrieve the `<variable>` of a `<species_type>` is called `<variable>_<species_type>`.
			As of today, the methods only accepts the IDs, not the labels. The reason is primaraly efficiency.
			The user should save the index at the beginning of the computation and use them every time they need to.

			The complete set of methods available is described in the documentation of the [AdimensionalSolutionExtractor](@ref specmicp::AdimensionalSystemSolutionExtractor).

			\defgroup specmicp_api SpecMiCP API
			\brief The classes and functions that forms the specmicp/reactmicp api

specmicp/doc/api.md87ddbd8df2daconfiguration

api.mdView Options

specmicp/doc/api.md
87ddbd8df2daconfiguration

api.md
View Options