Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/doc/api.md
fb013dc5d277reactmicp/unsaturated

Last Change
Raw File

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
c4science · Help