Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

unsaturated

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 reduced_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_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
			//! It is called reduced becaused it does not solve explicitely for the secondary species.
			//!
			class AdimensionalSystem :
			public micpsolver::MiCPProg<AdimensionalSystem>,
			public OptionsHandler<AdimensionalSystemOptions>,
			public units::UnitBaseClass
			{
			public:
			//! \brief Create a reduced system
			//!
			//! \param ptrthermo Shared pointer to the thermodynamic data
			//! \param totconc Vector containing the total concentrations (in mol/volume) for each components
			AdimensionalSystem(
			RawDatabasePtr ptrdata,
			const Vector& total_concentration,
			const AdimensionalSystemOptions& options=AdimensionalSystemOptions()
			);

			AdimensionalSystem(
			RawDatabasePtr ptrdata,
			const Vector& totconc,
			const AdimensionalSystemSolution& previous_solution,
			const AdimensionalSystemOptions& options=AdimensionalSystemOptions()
			);

			// 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
			scalar_t residual_water(const Vector& x) const;
			//! \brief Return the residual for the conservation equation of component (i)
			//!
			//! For (i==0), water use the method : 'residual_water'
			double residual_component(const Vector& x, index_t i) const;
			//! \brief Equilibrium condition for the minerals \f$ 1 - \mathrm{SI}_l \f$
			double residual_mineral(const Vector& x, index_t m) const;
			//! \brief Charge conservation
			double residual_charge_conservation(const Vector& x) const;

			//! \brief Return the residuals
			void get_residuals(const Vector& x, Vector& residual);
			//! \brief Return the jacobian
			void get_jacobian(Vector& x,
			Matrix& jacobian);


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

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

			//! \brief Return the equation id
			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[0];}
			//! \brief Return the equation number of component 'i'
			index_t ideq_paq(index_t i) const {
			specmicp_assert(i < m_data->nb_component);
			return m_ideq[i];
			}
			//! \brief Return the equation number of mineral 'm'
			index_t ideq_min(index_t m) const {
			specmicp_assert(m < m_data->nb_mineral);
			return m_ideq[m+m_data->nb_component];
			}

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

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

			//! \brief Return the saturation of water
			scalar_t saturation_water(const Vector& x) const;
			//! \brief Return the density of water
			scalar_t density_water() const;
			//! \brief Return the saturation (mol/volume) of a mineral
			scalar_t saturation_mineral(const Vector& x, index_t mineral) const;
			//! \brief Return the volume of minerals
			scalar_t molar_volume_mineral(index_t mineral) const;
			//! \brief Return the sum of saturation of the minerals
			scalar_t sum_saturation_minerals(const Vector& x) const;
			//! \brief Return the log_10 of component concentration
			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 concentration of 'component'
			scalar_t component_molality(const Vector& x, index_t component) const
			{
			return pow10(log_component_molality(x, component));
			}

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

			//! \brief Return the concentration of secondary specis 'aqueous'
			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
			bool is_aqueous_active(index_t aqueous) const
			{
			return m_active_aqueous[aqueous];
			}
			//! \brief Return log_10(γ_i) where i is a component
			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
			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'
			scalar_t total_concentration(index_t component) const {return m_tot_conc(component);}

			//! \brief Return the fugacity of 'gas'
			scalar_t gas_fugacity(index_t gas) const {
			return m_gas_fugacity(gas);
			}
			//! \brief Return true if gas is active
			bool is_active_gas(index_t gas) const {
			return m_active_gas[gas];
			}

			//! \brief Return the saturation of the gas phase
			scalar_t saturation_gas_phase() const {
			return m_saturation_gas;
			}

			//! \brief Return the gas pressure
			scalar_t gas_total_pressure() const {
			return 1.01325e5;
			}
			//! \brief Return the temperature
			scalar_t temperature() const {
			return 273.16+25;
			}


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

			//! \brief Return the equilibrium state of the system, the 'Solution' of the speciation problem
			AdimensionalSystemSolution get_solution(const Vector& xtot, const Vector& x);


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

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

			// For the normal algorithm
			void set_secondary_concentration(const Vector& x);
			//! \brief Compute the secondary variables
			void set_secondary_variables(const Vector& x);
			//! \brief Set the gas phase saturation
			void set_saturation_gas_phase(const Vector& x);
			//! \brief Set the fugacity pressure
			void set_pressure_fugacity(const Vector& x);

			//! \brief Compute the activity model, called by the solver at the beginning of an iteration
			bool hook_start_iteration(const Vector &x, scalar_t norm_residual);

			//! \brief Return a scale factor to avoid negative mass during Newton's 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
			void reasonable_starting_guess(Vector& x);
			//! \brief A reasonable (maybe...) restarting guess
			void reasonable_restarting_guess(Vector& x);

			//! \brief Reduce the mineral system if one cannot dissolve
			void reduce_mineral_problem(Vector& true_var);
			//! \brief Reset the mineral system to the true problem
			void reset_mineral_system(Vector& true_var, Vector& output_var);

			//! \brief set the charge keeper to 'component'
			//!
			//! The equation for the charge keeper is replaced by the charge balance
			//!
			//! Set component to 'no_species' to solve charge conservation implicitely.
			//! In this case, the total concentration must conserve the charge.
			void set_charge_keeper(index_t component)
			{
			specmicp_assert(component >= 1 and component < m_data->nb_component);
			get_options().charge_keeper = component;
			}

			//! \brief Return true if 'component' is thecharge keeper
			bool is_charge_keeper(index_t component)
			{
			return ( get_options().charge_keeper != no_species
			and component == get_options().charge_keeper
			);
			}

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

			// Equation numbering
			// ##################

			void number_eq();

			// 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_tot_conc;

			index_t m_nb_tot_vars;
			index_t m_nb_free_vars;
			index_t m_nb_compl_vars;

			std::vector<index_t> m_ideq;

			Vector m_secondary_conc;

			Vector m_loggamma;
			scalar_t m_ionic_strength;

			scalar_t m_saturation_gas;
			Vector m_gas_fugacity;

			bool not_in_linesearch;

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

			};

			} // end namespace specmicp


			#endif // SPECMICP_ADIMENSIONALSYSTEM_HPP

specmicp/src/specmicp/adimensional/adimensional_system.hpp3cd1304c75d9reactmicp/unsaturated

adimensional_system.hppView Options

specmicp/src/specmicp/adimensional/adimensional_system.hpp
3cd1304c75d9reactmicp/unsaturated

adimensional_system.hpp
View Options