Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

unsaturated

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 <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 {

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


			// =================================
			//
			// 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()));
			}

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

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

variables_interface.cppView Options

specmicp/src/reactmicp/systems/unsaturated/variables_interface.cpp
d536fc3e2489reactmicp/unsaturated

variables_interface.cpp
View Options