adimensional_system.hpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, May 9, 01:55

adimensional_system.hpp
View Options

	/*-------------------------------------------------------

	- Module : specmicp
	- File : adimensional_system.hpp
	- Author : Fabien Georget

	Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
	All rights reserved.

	Redistribution and use in source and binary forms, with or without
	modification, are permitted provided that the following conditions are met:
	* Redistributions of source code must retain the above copyright
	notice, this list of conditions and the following disclaimer.
	* Redistributions in binary form must reproduce the above copyright
	notice, this list of conditions and the following disclaimer in the
	documentation and/or other materials provided with the distribution.
	* Neither the name of the Princeton University nor the
	names of its contributors may be used to endorse or promote products
	derived from this software without specific prior written permission.

	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
	ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
	WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
	DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
	ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
	(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
	LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
	ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
	(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
	SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

	---------------------------------------------------------*/

	#ifndef SPECMICP_ADIMENSIONALSYSTEM_HPP
	#define SPECMICP_ADIMENSIONALSYSTEM_HPP

	//! \file adimensional_system.hpp The MiCP program to solve speciation

	#include "common.hpp"
	#include "database.hpp"
	#include "micpsolver/micpprog.hpp"
	#include "../simulation_box.hpp"
	#include "physics/units.hpp"
	#include "utils/options_handler.hpp"

	#include "adimensional_system_numbering.hpp"
	#include "adimensional_system_structs.hpp"

	//! \namespace specmicp main namespace for the SpecMiCP solver
	namespace specmicp {

	class AdimensionalSystemSolution;

	//! \brief The equilibrium speciation problem
	//!
	//! Represent the equilibrium speciation problem as a Mixed Complementarity program
	//!
	//! The main variables are :
	//! - \f$S^t_w\f$ the volume fraction of water (total saturation)
	//! - \f$b_i\f$ the molality of aqueous component
	//! - \f$S^t_m\f$ the volume fraction of mineral
	//!
	//! The solid phase equilibrium is solved by using a complementarity formulation
	//! \f[ S^t_m \geq 0, \; - \mathrm{SI}_m \geq 0 , \; \mathrm{and} \; - S^t_m \mathrm{SI}_m = 0 \f]
	//!
	//! The secondary variables (molalities of secondary aqueous species, gas fugacities, ionic strength, ...)
	//! are solved using a fixed point iteration at the beginning of each iteration.
	//!
	class AdimensionalSystem :
	public micpsolver::MiCPProg<AdimensionalSystem>,
	public AdimemsionalSystemNumbering,
	public OptionsHandler<AdimensionalSystemOptions>,
	public units::UnitBaseClass
	{
	public:
	//! \brief Create a reduced system
	//!
	//! \param ptrdata Shared pointer to the thermodynamic database
	//! \param constraints Constraints to apply to the system
	//! \param options Numerical options, mainly related to the secondary variables
	AdimensionalSystem(
	RawDatabasePtr ptrdata,
	const AdimensionalSystemConstraints& constraints,
	const AdimensionalSystemOptions& options=AdimensionalSystemOptions()
	);

	//! \brief Create a reduced system, initialize the system using a previous solution
	//!
	//! \param ptrdata Shared pointer to the thermodynamic database
	//! \param constraints Constraints to apply to the system
	//! \param previous_solution The previous solution to use for initialization
	//! \param options Numerical options, mainly related to the secondary variables
	AdimensionalSystem(
	RawDatabasePtr ptrdata,
	const AdimensionalSystemConstraints& constraints,
	const AdimensionalSystemSolution& previous_solution,
	const AdimensionalSystemOptions& options
	);

	// Variables
	// ==========

	//! \brief Return the total number of variables
	index_t total_variables() const {return m_nb_tot_vars;}
	//! \brief Return the number of 'free' variables (not subject to complementarity conditions)
	index_t nb_free_variables() const {return m_nb_free_vars;}
	//! \brief Return the number of variables subjected to the complementarity conditions
	index_t nb_complementarity_variables() const {return m_nb_compl_vars;}

	// The linear system
	// ==================

	//! \brief Return the residual for the water conservation equation
	//!
	//! \param x Vector of the main variables
	scalar_t residual_water(const Vector& x) const;
	//! \brief Return the residual for the conservation equation of component (i)
	//!
	//! For water (i==0) use the method : 'residual_water'
	//! \param x Vector of the main variables
	//! \param i Index of the component (in the database)
	scalar_t residual_component(const Vector& x, index_t i) const;
	//! \brief Compute the residual for the surface sorption
	//! \param x Vector of the main variables
	scalar_t residual_surface(const Vector& x) const;
	//! \brief Equilibrium condition for the minerals
	//! \param x Vector of the main variables
	//! \param i Index of the component (in the database)
	scalar_t residual_mineral(const Vector& x, index_t m) const;
	//! \brief Return the residual of the charge conservation equation
	//!
	//! This equation may be used instead of a mass conservation equation to
	//! explicitely enforce the electroneutrality.
	//! The component corresponding to this equation is called the charge keeper
	//! \param x Vector of the main variables
	scalar_t residual_charge_conservation(const Vector& x) const;
	//! \brief Return the residual corresponding to a fixed fugacity equation
	//! \param x Vector of the main variables
	//! \param component Index of the corresponding component (in the database)
	scalar_t residual_fixed_fugacity(const Vector& x, index_t component) const;
	//! \brief Return the residual corresponding to a fixed activity equation
	//!
	//! If 'component ' is charged, then a charge keeper should be set.
	//! \param x Vector of the main variables
	//! \param component Index of the corresponding component (in the database)
	scalar_t residual_fixed_activity(const Vector& x, index_t component) const;
	//! \brief Return the residuals
	//!
	//! \param x Vector of the main variables
	//! \param[out] residual Vector containing the residuals
	void get_residuals(const Vector& x, Vector& residual);
	//! \brief Return the jacobian
	//!
	//! \param x Vector of the main variables
	//! \param[out] jacobian Dense matrix representing the jacobian
	void get_jacobian(Vector& x,
	Matrix& jacobian);


	// Getters
	// #######

	// Equation id access
	// ==================

	//! \brief Return the equation id
	//! \param i Index of the main variable (in the set of the main variables)
	index_t ideq(index_t i) const {return m_ideq[i];}
	//! \brief Return the equation number for conservation of water
	index_t ideq_w() const {return m_ideq[dof_water()];}
	//! \brief Return the equation number of component 'i'
	//! \param i Index of the component (in the database)
	index_t ideq_paq(index_t i) const {
	specmicp_assert(i < m_data->nb_component);
	return m_ideq[dof_component(i)];
	}
	//! \brief Return the equation number of the free surface concentration
	index_t ideq_surf() const {
	return m_ideq[dof_surface()];
	}
	//! \brief Return the equation number of mineral 'm'
	//! \param m Index of the mineral (in the database)
	index_t ideq_min(index_t m) const {
	specmicp_assert(m < m_data->nb_mineral);
	return m_ideq[dof_mineral(m)];
	}

	// Equation type
	// -------------
	//! \brief Return the type of equation used for the solvent
	//!
	//! Can be no equation, mass conservation, saturation of the system, ...
	//!
	//! \sa #WaterEquationType, #aqueous_component_equation_type, #AqueousComponentEquationType
	WaterEquationType water_equation_type() const {
	return static_cast<WaterEquationType>(m_component_equation_type[0]);
	}
	//! \brief Return the type of equation for aqueous species 'component'
	//!
	//! For water used the method #water_equation_type
	//! \param component Index of the aqueous component (in the database)
	//!
	//! \sa #AqueousComponentEquationType, water_equation_type, #WaterEquationType
	AqueousComponentEquationType aqueous_component_equation_type(index_t component) const {
	specmicp_assert(component > 0 and component < m_data->nb_component);
	return static_cast<AqueousComponentEquationType>(m_component_equation_type[component]);
	}
	//! \brief Return true if an equation exist for this component
	//!
	//! \param component Index of the component (in the database)
	bool is_active_component(index_t component) const {
	specmicp_assert(component < m_data->nb_component and component > no_species);
	return m_ideq[component] != no_equation;
	}

	// Variables
	// =========

	// Primary
	// -------

	//! \brief Return the saturation of water
	//!
	//! \param x Vector of the main variables
	scalar_t saturation_water(const Vector& x) const;
	//! \brief Return the density of water
	scalar_t density_water() const;
	//! \brief Return the volume fraction of a mineral
	//!
	//! \param x Vector of the main variables
	//! \param mineral Index of the mineral (in the database)
	scalar_t saturation_mineral(const Vector& x, index_t mineral) const;
	//! \brief Return the volume of minerals
	//!
	//! \param mineral Index of the mineral (in the database)
	scalar_t molar_volume_mineral(index_t mineral) const;
	//! \brief Return the sum of saturation of the minerals
	//!
	//! This corresponds to the volume fraction of the solid phases (minus the inert phase)
	//! \param x Vector of the main variables
	scalar_t sum_saturation_minerals(const Vector& x) const;
	//! \brief Return the log_10 of the component molality
	//!
	//! \param x Vector of the main variables
	//! \param component Index of the aqueous component (in the database)
	scalar_t log_component_molality(const Vector& x, index_t component) const
	{
	specmicp_assert(ideq_paq(component) != no_equation and component < m_data->nb_component);
	return x(ideq_paq(component));
	}
	//! \brief Return the molality of 'component'
	//!
	//! \param x Vector of the main variables
	//! \param component Index of the aqueous component (in the database)
	scalar_t component_molality(const Vector& x, index_t component) const
	{
	return pow10(log_component_molality(x, component));
	}

	//! \brief Return the concentration of free sorption site
	scalar_t free_sorption_site_concentration(const Vector& x) const {
	return pow10(log_free_sorption_site_concentration(x));
	}
	//! \brief Return the log_10 of the free sorption site concentration
	scalar_t log_free_sorption_site_concentration(const Vector& x) const {
	specmicp_assert(ideq_surf() != no_equation);
	return x(ideq_surf());
	}

	// Secondary
	// ---------

	//! \brief Return the concentration of secondary species 'aqueous'
	//!
	//! \param aqueous Index of the secondary aqueous species (in the database)
	scalar_t secondary_molality(index_t aqueous) const
	{
	if (not m_active_aqueous[aqueous]) return 0.0;
	return m_secondary_conc(aqueous);
	}
	//! \brief Return true if 'aqueous' is active
	//!
	//! i.e. Return true if all its component are present in the system
	//! \param aqueous Index of the secondary aqueous species (in the database)
	bool is_aqueous_active(index_t aqueous) const
	{
	return m_active_aqueous[aqueous];
	}
	//! \brief Return log_10(γ_i) where i is a component
	//!
	//! γ is the activity coefficient
	//! \param component Index of the aqueous component (in the database)
	scalar_t log_gamma_component(index_t component) const {
	return m_loggamma(component);
	}
	//! \brief Return log_10(γ_j) where j is a secondary aqueous species
	//!
	//! γ is the activity coefficient
	//! \param aqueous Index of the secondary aqueous species (in the database)
	scalar_t log_gamma_secondary(index_t aqueous) const {
	return m_loggamma(m_data->nb_component + aqueous);
	}
	//! \brief Return the ionic strength
	scalar_t ionic_strength() const {
	return m_ionic_strength;
	}
	//! \brief Return the total concentration of 'component'
	//!
	//! Only use this method if the mass is conserved for 'component'
	//! \param component Index of the aqueous component (in the database)
	scalar_t total_concentration_bc(index_t component) const {
	specmicp_assert((component > 0 and aqueous_component_equation_type(component)
	== AqueousComponentEquationType::MassConservation)
	or
	( water_equation_type() == WaterEquationType::MassConservation));
	return m_fixed_values(component);
	}
	//! \brief Return the fixed activity of 'component'
	//!
	//! Only use this method if 'component' has a fixed activity constraint
	//! \param component Index of the aqueous component (in the database)
	scalar_t fixed_activity_bc(index_t component) const {
	specmicp_assert(aqueous_component_equation_type(component)
	== AqueousComponentEquationType::FixedActivity);
	return m_fixed_values(component);
	}
	//! \brief Return the fixed fugacity value for 'component'
	scalar_t fixed_fugacity_bc(index_t component) const {
	specmicp_assert(aqueous_component_equation_type(component)
	== AqueousComponentEquationType::FixedFugacity);
	return m_fixed_values(component);
	}

	// Gas
	// ---

	//! \brief Return the fugacity of 'gas'
	//!
	//! \param gas Index of the gas (in the database)
	scalar_t gas_fugacity(index_t gas) const {
	return m_gas_fugacity(gas);
	}
	//! \brief Return true if gas is active
	//!
	//! \param gas Index of the gas (in the database)
	bool is_active_gas(index_t gas) const {
	return m_active_gas[gas];
	}
	//! \brief Return the volume fraction (total saturation) of the gas phase
	scalar_t saturation_gas_phase() const {
	return m_saturation_gas;
	}

	// Sorbed species
	// --------------

	//! \brief Return the surface total concentration
	scalar_t surface_total_concentration() const {
	return m_surface_concentration;
	}
	//! \brief Return true if 'sorbed' is an active species
	//! \param sorbed Index of the sorbed species (in the database)
	bool is_active_sorbed(index_t sorbed) const {
	return m_active_sorbed[sorbed];
	}
	//! \brief Return the molality of the sorbed species 'sorbed'
	//! \param sorbed Index of the sorbed species (in the database)
	scalar_t sorbed_species_concentration(index_t sorbed) const {
	return m_sorbed_concentrations(sorbed);
	}

	// Pressure and temperature
	// ------------------------

	//! \brief Return the gas pressure
	//!
	//! This is a fixed value (1 atm)
	scalar_t gas_total_pressure() const {
	return units::convert_pressure(1.01325e5, length_unit());
	}
	//! \brief Return the temperature
	//!
	//! This is a fixed value (25°C)
	scalar_t temperature() const {
	return 273.16+25;
	}

	// Equilibrium State
	// =================

	//! \brief Return the equilibrium state of the system, the Solution of the speciation problem
	//!
	//! \param xtot the complete set of main variables
	//! \param x the reduced set of main variables, only the variables with an equation
	AdimensionalSystemSolution get_solution(Vector& xtot, const Vector& x);


	// Algorithm
	// #########

	//! \brief Compute the ionic strength
	//!
	//! \param x Vector of the main variables
	void set_ionic_strength(const Vector& x);
	//! \brief Compute the activity coefficients
	//!
	//! \param x Vector of the main variables
	void compute_log_gamma(const Vector& x);

	//! \brief Compute the secondary aqueous species molalities
	//!
	//! \param x Vector of the main variables
	void set_secondary_concentration(const Vector& x);
	//! \brief Compute the secondary variables
	//!
	//! \param x Vector of the main variables
	void set_secondary_variables(const Vector& x);
	//! \brief Compute the gas phase volume fraction
	//!
	//! \param x Vector of the main variables
	void set_saturation_gas_phase(const Vector& x);
	//! \brief Compute the fugacity for all the gas
	//!
	//! \param x Vector of the main variables
	void set_pressure_fugacity(const Vector& x);
	//! \brief Compute the sorbed species concentrations
	//!
	//! \param x Vector of the main variables
	void set_sorbed_concentrations(const Vector& x);
	//! \brief This function is called at the beginning of each iteration by the solver
	//!
	//! It does the fixed-point iteration to compute the secondary variables
	//! \param x Vector of the main variables
	//! \param norm_residual Norm of the current residuals
	bool hook_start_iteration(const Vector &x, scalar_t norm_residual);

	//! \brief Return a scale factor to avoid negative mass during Newton's iteration
	//!
	//! \param x Vector of the main variables
	//! \param update Update of the current iteration
	double max_lambda(const Vector &x, const Vector &update);
	//! \brief Return the component that does not exist in the solution (inactive dof)
	const std::vector<index_t>& get_non_active_component() {return m_nonactive_component;}
	//! \brief A reasonable (well... maybe...) starting guess
	//!
	//! \param xtot Vector of the main variables
	//! \sa #reasonable_starting_guess
	void reasonable_starting_guess(Vector& xtot);
	//! \brief A reasonable (maybe...) restarting guess
	//!
	//! \param xtot Vector of the main variables
	//! \sa #reasonable_restarting_guess
	void reasonable_restarting_guess(Vector& xtot);

	//! \brief Return true if 'component' is the charge keeper
	//!
	//! The charge keeper is the component added or removed to satisfy the charge balance
	//! \param component Index of the aqueous component (in the database)
	bool is_charge_keeper(index_t component)
	{
	return (aqueous_component_equation_type(component) == AqueousComponentEquationType::ChargeBalance);
	}

	// Private Members and attributes
	// ##############################
	private:

	// Jacobian
	// ########

	void jacobian_water(Vector& x, Matrix& jacobian);
	void jacobian_aqueous_components(Vector& x, Matrix& jacobian);
	void jacobian_minerals(Vector& x, Matrix& jacobian);
	void jacobian_surface(Vector& x, Matrix& jacobian);

	// Equation numbering
	// ##################
	//! \brief Number the equations
	void number_eq(const AdimensionalSystemConstraints& constraints);
	//! \brief Number the equations for the aqueous components
	void number_eq_aqueous_component(
	const AdimensionalSystemConstraints& constraints,
	index_t& neq
	);

	// Setter
	// ######

	// main variables
	// --------------

	// they are set by the solver

	// Secondary variables
	// -------------------

	//! \brief Return a reference to log_10(γ_i) where i is a component
	scalar_t& log_gamma_component(index_t component) {
	return m_loggamma(component);
	}
	//! \brief Return a reference to log_10(γ_j) where j is a secondary aqueous species
	scalar_t& log_gamma_secondary(index_t aqueous) {
	return m_loggamma(m_data->nb_component + aqueous);
	}
	//! \brief Return a reference to the ionic strength
	scalar_t& ionic_strength() {
	return m_ionic_strength;
	}

	// Attributes
	// ##########

	RawDatabasePtr m_data;
	Vector m_fixed_values;

	index_t m_nb_tot_vars;
	index_t m_nb_free_vars;
	index_t m_nb_compl_vars;

	std::vector<index_t> m_ideq;
	std::vector<int> m_component_equation_type;
	std::vector<index_t> m_fixed_activity_species;

	Vector m_secondary_conc;

	Vector m_loggamma;
	scalar_t m_ionic_strength;

	scalar_t m_saturation_gas;
	Vector m_gas_fugacity;

	scalar_t m_surface_concentration;
	Vector m_sorbed_concentrations;

	bool not_in_linesearch;

	std::vector<index_t> m_nonactive_component;
	std::vector<bool> m_active_aqueous;
	std::vector<bool> m_active_gas;
	std::vector<bool> m_active_sorbed;

	scalar_t m_inert_volume_fraction;

	};

	} // end namespace specmicp


	#endif // SPECMICP_ADIMENSIONALSYSTEM_HPP

adimensional_system.hppNo OneTemporaryActions

File Metadata

adimensional_system.hppView Options

Event Timeline

adimensional_system.hpp
No OneTemporary
Actions

adimensional_system.hpp
View Options