Diffusion SpecMiCP / ReactMiCP (chloride) Browse

variables_interface.cpp
78496b77e932chloride

variables_interface.cpp
View Options

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

			Copyright (c) 2015 F. 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:

			1. Redistributions of source code must retain the above copyright notice, this
			list of conditions and the following disclaimer.

			2. 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.

			3. Neither the name of the copyright holder 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 HOLDER 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.

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

			#include "variables_interface.hpp"

			#include "variables.hpp"
			#include "variables_box.hpp"

			#include "../../solver/staggers_base/variables_base.hpp"

			#include "../../../specmicp/adimensional/adimensional_system_solution_extractor.hpp"

			#include "../../../utils/compat.hpp"
			#include "../../../utils/log.hpp"
			#include "../../../utils/value_checker.hpp"
			#include "../../../utils/vector_checker.hpp"

			#include <stdexcept>

			// This is a boring class, with only small differences between
			// all the functions and variables. But, they make all the
			// difference at the end. The main purpose of the class is to
			// hide all the variables structure and only present a unified
			// interface to the user.

			// None of the structure defined in "variables_sub.hpp" should
			// be leaked outside this class

			namespace specmicp {
			namespace reactmicp {
			namespace systems {
			namespace unsaturated {

			using namespace specmicp::vector_checker;


			// typedef
			// -------

			using VariablesBase = solver::VariablesBase;
			using VariablesBasePtr = std::shared_ptr<VariablesBase>;

			// =================================
			//
			// Declaration of the implementation
			//
			// =================================

			// This class just store and retrieve variables group

			struct VariablesInterface::VariablesInterfaceImpl
			{
			UnsaturatedVariablesPtr m_vars;
			mesh::Mesh1DPtr m_mesh;
			database::RawDatabasePtr m_database;


			VariablesInterfaceImpl(
			mesh::Mesh1DPtr the_mesh,
			database::RawDatabasePtr the_database,
			const std::vector<index_t>& component_with_gas,
			units::UnitsSet units_set = {}
			);


			MainVariable& liquid_saturation() {
			return m_vars->get_liquid_saturation();
			}
			const MainVariable& liquid_saturation() const {
			return m_vars->get_liquid_saturation();
			}

			MainVariable& aqueous_concentration(index_t component) {
			return m_vars->get_aqueous_concentration(component);
			}
			const MainVariable& aqueous_concentration(index_t component) const {
			return m_vars->get_aqueous_concentration(component);
			}

			MainVariable& partial_pressure(index_t component) {
			if (not m_vars->component_has_gas(component)) {
			throw std::logic_error("Component has no associated gas");
			}
			return m_vars->get_pressure_main_variables(component);
			}
			const MainVariable& partial_pressure(index_t component) const {
			return m_vars->get_pressure_main_variables(component);
			}


			MainVariable& solid_concentration(index_t component) {
			return m_vars->get_solid_concentration(component);
			}
			const MainVariable& solid_concentration(index_t component) const {
			return m_vars->get_solid_concentration(component);
			}

			SecondaryTransientVariable& porosity() {
			return m_vars->get_porosity();
			}
			const SecondaryTransientVariable& porosity() const {
			return m_vars->get_porosity();
			}


			SecondaryTransientVariable& water_aqueous_concentration() {
			return m_vars->get_water_aqueous_concentration();
			}
			const SecondaryTransientVariable& water_aqueous_concentration() const {
			return m_vars->get_water_aqueous_concentration();
			}

			SecondaryVariable& liquid_diffusivity() {
			return m_vars->get_liquid_diffusivity();
			}
			const SecondaryVariable& liquid_diffusivity() const {
			return m_vars->get_liquid_diffusivity();
			}

			SecondaryVariable& liquid_permeability() {
			return m_vars->get_liquid_permeability();
			}
			const SecondaryVariable& liquid_permeability() const {
			return m_vars->get_liquid_permeability();
			}


			scalar_t& binary_gas_diffusivity(index_t component) {
			return m_vars->get_binary_gas_diffusivity(component);
			}
			scalar_t binary_gas_diffusivity(index_t component) const {
			return m_vars->get_binary_gas_diffusivity(component);
			}

			SecondaryVariable& resistance_gas_diffusivity() {
			return m_vars->get_resistance_gas_diffusivity();
			}
			const SecondaryVariable& resistance_gas_diffusivity() const {
			return m_vars->get_resistance_gas_diffusivity();
			}

			SecondaryVariable& advection_flux() {
			return m_vars->get_advection_flux();
			}
			const SecondaryVariable& advection_flux() const {
			return m_vars->get_advection_flux();
			}


			VariablesValidity check_variables();
			bool check_bounds(const Vector& to_check, scalar_t lower, scalar_t upper);

			VariablesValidity check_saturation();
			VariablesValidity check_relative();
			VariablesValidity check_porosity();
			VariablesValidity check_transport_param();
			};


			// =================================
			//
			// Implementation
			//
			// =================================

			VariablesInterface::VariablesInterface(
			mesh::Mesh1DPtr the_mesh,
			database::RawDatabasePtr the_database,
			const std::vector<index_t>& component_with_gas):
			m_impl(make_unique<VariablesInterfaceImpl>(
			the_mesh, the_database, component_with_gas))
			{

			}


			VariablesInterface::VariablesInterface(
			mesh::Mesh1DPtr the_mesh,
			database::RawDatabasePtr the_database,
			const std::vector<index_t>& component_with_gas,
			const units::UnitsSet& units_set
			):
			m_impl(make_unique<VariablesInterfaceImpl>(
			the_mesh, the_database, component_with_gas, units_set))
			{

			}

			VariablesInterface::VariablesInterfaceImpl::VariablesInterfaceImpl(
			mesh::Mesh1DPtr the_mesh,
			database::RawDatabasePtr the_database,
			const std::vector<index_t>& component_with_gas,
			units::UnitsSet units_set
			):
			m_mesh(the_mesh),
			m_database(the_database)
			{
			std::vector<bool> comp_has_gas(the_database->nb_component(), false);
			for (auto component: component_with_gas) {
			comp_has_gas[component] = true;
			}
			m_vars = std::make_shared<UnsaturatedVariables>(
			the_mesh, the_database, comp_has_gas, units_set.length);

			}

			VariablesInterface::~VariablesInterface() = default;

			std::shared_ptr<VariablesBase> VariablesInterface::get_base_variables()
			{
			return std::static_pointer_cast<VariablesBase>(m_impl->m_vars);
			}

			std::shared_ptr<UnsaturatedVariables> VariablesInterface::get_variables()
			{
			return m_impl->m_vars;
			}


			UnsaturatedVariables* VariablesInterface::get_raw_variables()
			{
			return m_impl->m_vars.get();
			}

			// Saturation

			const Vector& VariablesInterface::get_liquid_saturation() const
			{
			return m_impl->liquid_saturation().variable;
			}

			void VariablesInterface::set_liquid_saturation(index_t node, scalar_t value)
			{
			m_impl->liquid_saturation().variable(node) = value;
			m_impl->liquid_saturation().predictor(node) = value;
			}

			void VariablesInterface::set_liquid_saturation(const Vector& value)
			{
			m_impl->liquid_saturation().variable = value;
			m_impl->liquid_saturation().predictor = value;
			}

			void VariablesInterface::set_liquid_saturation(scalar_t value)
			{
			m_impl->liquid_saturation().variable.setConstant(value);
			m_impl->liquid_saturation().predictor.setConstant(value);
			}

			// Aqueous component concentration


			const Vector& VariablesInterface::get_aqueous_concentration(index_t component) const
			{
			return m_impl->aqueous_concentration(component).variable;
			}

			void VariablesInterface::set_aqueous_concentration(index_t component, index_t node, scalar_t value)
			{
			MainVariable& aq_conc = m_impl->aqueous_concentration(component);
			aq_conc.variable(node) = value;
			aq_conc.predictor(node) = value;
			}

			void VariablesInterface::set_aqueous_concentration(index_t component, const Vector& value)
			{
			MainVariable& aq_conc = m_impl->aqueous_concentration(component);
			aq_conc.variable = value;
			aq_conc.predictor = value;
			}

			void VariablesInterface::set_aqueous_concentration(index_t component, scalar_t value)
			{
			MainVariable& aq_conc = m_impl->aqueous_concentration(component);
			aq_conc.variable.setConstant(value);
			aq_conc.predictor.setConstant(value);
			}


			// Partial pressure


			const Vector& VariablesInterface::get_partial_pressure(index_t component) const
			{
			return m_impl->partial_pressure(component).variable;
			}

			void VariablesInterface::set_partial_pressure(index_t component, index_t node, scalar_t value)
			{
			MainVariable& comp_pressure = m_impl->partial_pressure(component);
			comp_pressure.variable(node) = value;
			comp_pressure.predictor(node) = value;
			}

			void VariablesInterface::set_partial_pressure(index_t component, const Vector& value)
			{
			MainVariable& comp_pressure = m_impl->partial_pressure(component);
			comp_pressure.variable = value;
			comp_pressure.predictor = value;
			}

			void VariablesInterface::set_partial_pressure(index_t component, scalar_t value)
			{
			MainVariable& comp_pressure = m_impl->partial_pressure(component);
			comp_pressure.variable.setConstant(value);
			comp_pressure.predictor.setConstant(value);
			}


			// Solid concentration


			const Vector& VariablesInterface::get_solid_concentration(index_t component) const
			{
			return m_impl->solid_concentration(component).variable;
			}

			void VariablesInterface::set_solid_concentration(index_t component, index_t node, scalar_t value)
			{
			MainVariable& solid_conc = m_impl->solid_concentration(component);
			solid_conc.variable(node) = value;
			solid_conc.predictor(node) = value;
			}

			void VariablesInterface::set_solid_concentration(index_t component, const Vector& value)
			{
			MainVariable& solid_conc = m_impl->solid_concentration(component);
			solid_conc.variable = value;
			solid_conc.predictor = value;
			}

			void VariablesInterface::set_solid_concentration(index_t component, scalar_t value)
			{
			MainVariable& solid_conc = m_impl->solid_concentration(component);
			solid_conc.variable.setConstant(value);
			solid_conc.predictor.setConstant(value);
			}


			// Porosity

			const Vector& VariablesInterface::get_porosity() const
			{
			return m_impl->porosity().variable;
			}

			void VariablesInterface::set_porosity(index_t node, scalar_t value)
			{
			m_impl->porosity().variable(node) = value;
			m_impl->porosity().predictor(node) = value;
			}

			void VariablesInterface::set_porosity(const Vector& value)
			{
			m_impl->porosity().variable = value;
			m_impl->porosity().predictor = value;
			}

			void VariablesInterface::set_porosity(scalar_t value)
			{
			m_impl->porosity().variable.setConstant(value);
			m_impl->porosity().predictor.setConstant(value);
			}

			// Water aqueous concentration


			const Vector& VariablesInterface::get_water_aqueous_concentration() const
			{
			return m_impl->water_aqueous_concentration().variable;
			}

			void VariablesInterface::set_water_aqueous_concentration(index_t node, scalar_t value)
			{
			m_impl->water_aqueous_concentration().variable(node) = value;
			m_impl->water_aqueous_concentration().predictor(node) = value;
			}

			void VariablesInterface::set_water_aqueous_concentration(const Vector& value)
			{
			m_impl->water_aqueous_concentration().variable = value;
			m_impl->water_aqueous_concentration().predictor = value;
			}

			void VariablesInterface::set_water_aqueous_concentration(scalar_t value)
			{
			m_impl->water_aqueous_concentration().variable.setConstant(value);
			m_impl->water_aqueous_concentration().predictor.setConstant(value);
			}


			// Liquid diffusion coefficient

			const Vector& VariablesInterface::get_liquid_diffusivity() const
			{
			return m_impl->liquid_diffusivity().variable;
			}

			void VariablesInterface::set_liquid_diffusivity(index_t node, scalar_t value)
			{
			m_impl->liquid_diffusivity().variable(node) = value;
			}

			void VariablesInterface::set_liquid_diffusivity(const Vector& value)
			{
			m_impl->liquid_diffusivity().variable = value;
			}

			void VariablesInterface::set_liquid_diffusivity(scalar_t value)
			{
			m_impl->liquid_diffusivity().variable.setConstant(value);
			}

			// Liquid permeability

			const Vector& VariablesInterface::get_liquid_permeability() const
			{
			return m_impl->liquid_permeability().variable;
			}

			void VariablesInterface::set_liquid_permeability(index_t node, scalar_t value)
			{
			m_impl->liquid_permeability().variable(node) = value;
			}

			void VariablesInterface::set_liquid_permeability(const Vector& value)
			{
			m_impl->liquid_permeability().variable = value;
			}

			void VariablesInterface::set_liquid_permeability(scalar_t value)
			{
			m_impl->liquid_permeability().variable.setConstant(value);
			}

			// Gas diffusivity

			scalar_t
			VariablesInterface::get_binary_gas_diffusivity(
			index_t component
			) const
			{
			return m_impl->binary_gas_diffusivity(component);
			}

			void
			VariablesInterface::set_binary_gas_diffusivity(
			index_t component,
			scalar_t value
			)
			{
			m_impl->binary_gas_diffusivity(component) = value;
			}


			const Vector&
			VariablesInterface::get_resistance_gas_diffusivity() const
			{
			return m_impl->resistance_gas_diffusivity().variable;
			}

			void
			VariablesInterface::set_resistance_gas_diffusivity(
			index_t node,
			scalar_t value
			)
			{
			m_impl->resistance_gas_diffusivity().variable(node) = value;
			}
			void
			VariablesInterface::set_resistance_gas_diffusivity(
			const Vector& value)
			{
			m_impl->resistance_gas_diffusivity().variable = value;
			}
			void
			VariablesInterface::set_resistance_gas_diffusivity(
			scalar_t value)
			{
			m_impl->resistance_gas_diffusivity().variable.setConstant(value);
			}

			const Vector& VariablesInterface::get_advection_flux() const
			{
			return m_impl->advection_flux().variable;
			}

			void VariablesInterface::set_advection_flux(index_t node, scalar_t value)
			{
			m_impl->advection_flux().variable(node) = value;
			}

			void VariablesInterface::set_advection_flux(const Vector& value)
			{
			m_impl->advection_flux().variable = value;
			}

			void VariablesInterface::set_advection_flux(scalar_t value)
			{
			m_impl->advection_flux().variable.setConstant(value);
			}

			// User Models

			const Vector& VariablesInterface::get_capillary_pressure() const
			{
			return m_impl->m_vars->get_capillary_pressure().variable;
			}

			void VariablesInterface::set_capillary_pressure_model(
			user_model_saturation_f capillary_pressure_model
			)
			{
			m_impl->m_vars->set_capillary_pressure_model(capillary_pressure_model);
			const Vector& sat = get_liquid_saturation();
			Vector& cap_pressure = m_impl->m_vars->get_capillary_pressure().variable;
			for (auto node: m_impl->m_mesh->range_nodes())
			{
			cap_pressure(node) = capillary_pressure_model(node, sat(node));
			}
			}

			const Vector& VariablesInterface::get_vapor_pressure() const
			{
			return m_impl->m_vars->get_pressure_main_variables(0).variable;
			}
			void VariablesInterface::set_vapor_pressure_model(
			user_model_saturation_f vapor_pressure_model
			)
			{
			m_impl->m_vars->set_vapor_pressure_model(vapor_pressure_model);
			}


			const Vector& VariablesInterface::get_relative_liquid_diffusivity() const
			{
			return m_impl->m_vars->get_relative_liquid_diffusivity().variable;
			}

			void VariablesInterface::set_relative_liquid_diffusivity_model(
			user_model_saturation_f relative_liquid_diffusivity_model
			)
			{
			m_impl->m_vars->set_relative_liquid_diffusivity_model(relative_liquid_diffusivity_model);
			const Vector& sat = get_liquid_saturation();
			Vector& rel_diff = m_impl->m_vars->get_relative_liquid_diffusivity().variable;
			for (auto node: m_impl->m_mesh->range_nodes())
			{
			rel_diff(node) = relative_liquid_diffusivity_model(node, sat(node));
			}
			}

			const Vector& VariablesInterface::get_relative_liquid_permeability() const
			{
			return m_impl->m_vars->get_relative_liquid_permeability().variable;
			}

			void VariablesInterface::set_relative_liquid_permeability_model(
			user_model_saturation_f relative_liquid_permeability_model
			)
			{
			m_impl->m_vars->set_relative_liquid_permeability_model(relative_liquid_permeability_model);
			const Vector& sat = get_liquid_saturation();
			Vector& rel_diff = m_impl->m_vars->get_relative_liquid_permeability().variable;
			for (auto node: m_impl->m_mesh->range_nodes())
			{
			rel_diff(node) = relative_liquid_permeability_model(node, sat(node));
			}
			}

			const Vector& VariablesInterface::get_relative_gas_diffusivity() const
			{
			return m_impl->m_vars->get_relative_gas_diffusivity().variable;
			}

			void VariablesInterface::set_relative_gas_diffusivity_model(
			user_model_saturation_f relative_gas_diffusivity_model
			)
			{
			m_impl->m_vars->set_relative_gas_diffusivity_model(relative_gas_diffusivity_model);
			const Vector& sat = get_liquid_saturation();
			Vector& rel_diff = m_impl->m_vars->get_relative_gas_diffusivity().variable;
			for (auto node: m_impl->m_mesh->range_nodes())
			{
			rel_diff(node) = relative_gas_diffusivity_model(node, sat(node));
			}
			}


			// chemistry
			void VariablesInterface::initialize_variables(index_t node,
			const AdimensionalSystemSolutionExtractor& extractor
			)
			{
			set_porosity(node, extractor.porosity());
			scalar_t saturation = extractor.saturation_water();
			set_liquid_saturation(node, extractor.saturation_water());
			scalar_t rho_l = extractor.density_water();

			set_water_aqueous_concentration(node,
			rho_l*extractor.total_aqueous_concentration(0));
			set_solid_concentration(0, node,
			extractor.total_solid_concentration(0));
			if (m_impl->m_vars->component_has_gas(0))
			{
			set_partial_pressure(0, node,
			m_impl->m_vars->get_vapor_pressure_model()(node, saturation)
			);
			}

			for (index_t component: m_impl->m_database->range_aqueous_component())
			{
			set_aqueous_concentration(component, node,
			rho_l*extractor.total_aqueous_concentration(component));
			set_solid_concentration(component, node,
			extractor.total_solid_concentration(component));
			if (m_impl->m_vars->component_has_gas(component))
			{
			set_partial_pressure(component, node,
			extractor.fugacity_gas(m_impl->m_vars->get_id_gas(component))
			*m_impl->m_vars->get_total_pressure());
			}
			}

			m_impl->m_vars->set_adim_solution(node,
			std::move(extractor.get_solution()));
			}

			VariablesValidity VariablesInterface::check_variables()
			{
			return m_impl->check_variables();
			}


			bool VariablesInterface::VariablesInterfaceImpl::check_bounds(
			const Vector& to_check,
			scalar_t lower,
			scalar_t upper
			)
			{
			bool flag = true;
			for (index_t ind=0; ind< to_check.rows(); ++ind)
			{
			const scalar_t& value = to_check(ind);
			if (value < lower or value > upper) {
			flag = false;
			break;
			}
			}
			return flag;
			}

			VariablesValidity max_error_level(VariablesValidity err1, VariablesValidity err2)
			{
			return (err1>err2)?err1:err2;
			}

			VariablesValidity VariablesInterface::VariablesInterfaceImpl::check_variables()
			{

			VariablesValidity return_code = VariablesValidity::good;

			auto tmpret = check_saturation();
			return_code = max_error_level(return_code, tmpret);

			tmpret = check_porosity();
			return_code = max_error_level(return_code, tmpret);

			tmpret = check_relative();
			return_code = max_error_level(return_code, tmpret);

			return return_code;
			}

			static bool SPECMICP_DLL_LOCAL saturation_check(const Vector& x)
			{
			return (is_bounded(x, 0.0, 1.0) && (any(x) > 0.0));
			}

			VariablesValidity
			VariablesInterface::VariablesInterfaceImpl::check_saturation()
			{
			using namespace vector_checker;
			const auto& x = m_vars->get_liquid_saturation().variable;

			if (not saturation_check(x))
			{
			ERROR << "Saturation vector is incorrect !";
			return VariablesValidity::error;
			}
			return VariablesValidity::good;
			}

			static bool SPECMICP_DLL_LOCAL porosity_check(const Vector& x)
			{
			return (is_bounded(x, 0.0, 1.0) && (any(x) > 0.0));
			}

			VariablesValidity
			VariablesInterface::VariablesInterfaceImpl::check_porosity()
			{
			using namespace vector_checker;
			const auto& x = m_vars->get_porosity().variable;

			if (not porosity_check(x)) {
			ERROR << "Porosity is incorrect";
			return VariablesValidity::error;
			}
			return VariablesValidity::good;
			}

			static bool SPECMICP_DLL_LOCAL relative_variable_check(const Vector& x)
			{
			return (is_bounded(x, 0.0, 1.0) && (any(x) > 0.0));
			}

			VariablesValidity
			VariablesInterface::VariablesInterfaceImpl::check_relative()
			{
			auto is_valid = VariablesValidity::good;

			if (not relative_variable_check(
			m_vars->get_relative_liquid_diffusivity().variable)) {
			ERROR << "Relative liquid diffusivity is incorrect";
			is_valid = VariablesValidity::error;
			}
			if (not relative_variable_check(
			m_vars->get_relative_liquid_permeability().variable)) {
			ERROR << "Relative liquid permeability is incorrect";
			is_valid = VariablesValidity::error;
			}
			if (not relative_variable_check(
			m_vars->get_relative_gas_diffusivity().variable)) {
			ERROR << "Relative gas diffusivity is incorrect";
			is_valid = VariablesValidity::error;
			}
			if (not relative_variable_check(
			m_vars->get_resistance_gas_diffusivity().variable)) {
			ERROR << "Resistance gas diffusivity is incorrect";
			is_valid = VariablesValidity::error;
			}

			return is_valid;
			}

			static bool SPECMICP_DLL_LOCAL transport_parameters_check(const Vector& x)
			{
			return (is_lower_bounded(x, 0.0) && (any(x) > 0.0));
			}

			VariablesValidity
			VariablesInterface::VariablesInterfaceImpl::check_transport_param()
			{
			auto is_valid = VariablesValidity::good;

			if (not transport_parameters_check(
			m_vars->get_liquid_diffusivity().variable)) {
			ERROR << "Liquid diffusivity is incorrect";
			is_valid = VariablesValidity::error;
			}
			if (not transport_parameters_check(
			m_vars->get_liquid_permeability().variable)) {
			ERROR << "Liquid permeability is incorrect";
			is_valid = VariablesValidity::error;
			}

			return is_valid;
			}

			} //end namespace unsaturated
			} //end namespace systems
			} //end namespace reactmicp
			} //end namespace specmicp

specmicp/src/reactmicp/systems/unsaturated/variables_interface.cpp78496b77e932chloride

variables_interface.cppView Options

specmicp/src/reactmicp/systems/unsaturated/variables_interface.cpp
78496b77e932chloride

variables_interface.cpp
View Options